# | Line 48 | Line 48 | |
---|---|---|
48 | ||
49 | #include <algorithm> | |
50 | #include <set> | |
51 | + | #include <map> |
52 | ||
53 | #include "brains/SimInfo.hpp" | |
54 | #include "math/Vector3.hpp" | |
55 | #include "primitives/Molecule.hpp" | |
56 | + | #include "primitives/StuntDouble.hpp" |
57 | #include "UseTheForce/fCutoffPolicy.h" | |
58 | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" | |
59 | + | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" |
60 | + | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" |
61 | #include "UseTheForce/doForces_interface.h" | |
62 | + | #include "UseTheForce/DarkSide/neighborLists_interface.h" |
63 | #include "UseTheForce/DarkSide/electrostatic_interface.h" | |
64 | < | #include "UseTheForce/notifyCutoffs_interface.h" |
64 | > | #include "UseTheForce/DarkSide/switcheroo_interface.h" |
65 | #include "utils/MemoryUtils.hpp" | |
66 | #include "utils/simError.h" | |
67 | #include "selection/SelectionManager.hpp" | |
68 | + | #include "io/ForceFieldOptions.hpp" |
69 | + | #include "UseTheForce/ForceField.hpp" |
70 | ||
71 | + | |
72 | #ifdef IS_MPI | |
73 | #include "UseTheForce/mpiComponentPlan.h" | |
74 | #include "UseTheForce/DarkSide/simParallel_interface.h" | |
75 | #endif | |
76 | ||
77 | namespace oopse { | |
78 | + | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { |
79 | + | std::map<int, std::set<int> >::iterator i = container.find(index); |
80 | + | std::set<int> result; |
81 | + | if (i != container.end()) { |
82 | + | result = i->second; |
83 | + | } |
84 | ||
85 | < | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, |
86 | < | ForceField* ff, Globals* simParams) : |
87 | < | stamps_(stamps), forceField_(ff), simParams_(simParams), |
88 | < | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
85 | > | return result; |
86 | > | } |
87 | > | |
88 | > | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : |
89 | > | forceField_(ff), simParams_(simParams), |
90 | > | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
91 | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | |
92 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
93 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), | |
94 | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), | |
95 | < | sman_(NULL), fortranInitialized_(false) { |
95 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { |
96 | ||
81 | – | |
82 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; |
97 | MoleculeStamp* molStamp; | |
98 | int nMolWithSameStamp; | |
99 | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | |
# | Line 87 | Line 101 | namespace oopse { | |
101 | CutoffGroupStamp* cgStamp; | |
102 | RigidBodyStamp* rbStamp; | |
103 | int nRigidAtoms = 0; | |
104 | < | |
105 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { |
106 | < | molStamp = i->first; |
107 | < | nMolWithSameStamp = i->second; |
104 | > | std::vector<Component*> components = simParams->getComponents(); |
105 | > | |
106 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
107 | > | molStamp = (*i)->getMoleculeStamp(); |
108 | > | nMolWithSameStamp = (*i)->getNMol(); |
109 | ||
110 | addMoleculeStamp(molStamp, nMolWithSameStamp); | |
111 | ||
112 | //calculate atoms in molecules | |
113 | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | |
114 | ||
100 | – | |
115 | //calculate atoms in cutoff groups | |
116 | int nAtomsInGroups = 0; | |
117 | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | |
118 | ||
119 | for (int j=0; j < nCutoffGroupsInStamp; j++) { | |
120 | < | cgStamp = molStamp->getCutoffGroup(j); |
120 | > | cgStamp = molStamp->getCutoffGroupStamp(j); |
121 | nAtomsInGroups += cgStamp->getNMembers(); | |
122 | } | |
123 | ||
# | Line 116 | Line 130 | namespace oopse { | |
130 | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | |
131 | ||
132 | for (int j=0; j < nRigidBodiesInStamp; j++) { | |
133 | < | rbStamp = molStamp->getRigidBody(j); |
133 | > | rbStamp = molStamp->getRigidBodyStamp(j); |
134 | nAtomsInRigidBodies += rbStamp->getNMembers(); | |
135 | } | |
136 | ||
# | Line 155 | Line 169 | namespace oopse { | |
169 | } | |
170 | molecules_.clear(); | |
171 | ||
158 | – | delete stamps_; |
172 | delete sman_; | |
173 | delete simParams_; | |
174 | delete forceField_; | |
# | Line 262 | Line 275 | namespace oopse { | |
275 | } | |
276 | } | |
277 | ||
278 | < | }//end for (integrableObject) |
279 | < | }// end for (mol) |
278 | > | } |
279 | > | } |
280 | ||
281 | // n_constraints is local, so subtract them on each processor | |
282 | ndf_local -= nConstraints_; | |
# | Line 280 | Line 293 | namespace oopse { | |
293 | ||
294 | } | |
295 | ||
296 | + | int SimInfo::getFdf() { |
297 | + | #ifdef IS_MPI |
298 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
299 | + | #else |
300 | + | fdf_ = fdf_local; |
301 | + | #endif |
302 | + | return fdf_; |
303 | + | } |
304 | + | |
305 | void SimInfo::calcNdfRaw() { | |
306 | int ndfRaw_local; | |
307 | ||
# | Line 342 | Line 364 | namespace oopse { | |
364 | int b; | |
365 | int c; | |
366 | int d; | |
367 | + | |
368 | + | std::map<int, std::set<int> > atomGroups; |
369 | + | |
370 | + | Molecule::RigidBodyIterator rbIter; |
371 | + | RigidBody* rb; |
372 | + | Molecule::IntegrableObjectIterator ii; |
373 | + | StuntDouble* integrableObject; |
374 | ||
375 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
376 | + | integrableObject = mol->nextIntegrableObject(ii)) { |
377 | + | |
378 | + | if (integrableObject->isRigidBody()) { |
379 | + | rb = static_cast<RigidBody*>(integrableObject); |
380 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
381 | + | std::set<int> rigidAtoms; |
382 | + | for (int i = 0; i < atoms.size(); ++i) { |
383 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
384 | + | } |
385 | + | for (int i = 0; i < atoms.size(); ++i) { |
386 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
387 | + | } |
388 | + | } else { |
389 | + | std::set<int> oneAtomSet; |
390 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
391 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
392 | + | } |
393 | + | } |
394 | + | |
395 | + | |
396 | + | |
397 | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
398 | a = bond->getAtomA()->getGlobalIndex(); | |
399 | b = bond->getAtomB()->getGlobalIndex(); | |
# | Line 353 | Line 404 | namespace oopse { | |
404 | a = bend->getAtomA()->getGlobalIndex(); | |
405 | b = bend->getAtomB()->getGlobalIndex(); | |
406 | c = bend->getAtomC()->getGlobalIndex(); | |
407 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
408 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
409 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
410 | ||
411 | < | exclude_.addPair(a, b); |
412 | < | exclude_.addPair(a, c); |
413 | < | exclude_.addPair(b, c); |
411 | > | exclude_.addPairs(rigidSetA, rigidSetB); |
412 | > | exclude_.addPairs(rigidSetA, rigidSetC); |
413 | > | exclude_.addPairs(rigidSetB, rigidSetC); |
414 | > | |
415 | > | //exclude_.addPair(a, b); |
416 | > | //exclude_.addPair(a, c); |
417 | > | //exclude_.addPair(b, c); |
418 | } | |
419 | ||
420 | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
# | Line 364 | Line 422 | namespace oopse { | |
422 | b = torsion->getAtomB()->getGlobalIndex(); | |
423 | c = torsion->getAtomC()->getGlobalIndex(); | |
424 | d = torsion->getAtomD()->getGlobalIndex(); | |
425 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
426 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
427 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
428 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
429 | ||
430 | + | exclude_.addPairs(rigidSetA, rigidSetB); |
431 | + | exclude_.addPairs(rigidSetA, rigidSetC); |
432 | + | exclude_.addPairs(rigidSetA, rigidSetD); |
433 | + | exclude_.addPairs(rigidSetB, rigidSetC); |
434 | + | exclude_.addPairs(rigidSetB, rigidSetD); |
435 | + | exclude_.addPairs(rigidSetC, rigidSetD); |
436 | + | |
437 | + | /* |
438 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
439 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
440 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
441 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
442 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
443 | + | exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
444 | + | |
445 | + | |
446 | exclude_.addPair(a, b); | |
447 | exclude_.addPair(a, c); | |
448 | exclude_.addPair(a, d); | |
449 | exclude_.addPair(b, c); | |
450 | exclude_.addPair(b, d); | |
451 | exclude_.addPair(c, d); | |
452 | + | */ |
453 | } | |
454 | ||
376 | – | Molecule::RigidBodyIterator rbIter; |
377 | – | RigidBody* rb; |
455 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
456 | std::vector<Atom*> atoms = rb->getAtoms(); | |
457 | for (int i = 0; i < atoms.size() -1 ; ++i) { | |
# | Line 399 | Line 476 | namespace oopse { | |
476 | int b; | |
477 | int c; | |
478 | int d; | |
479 | + | |
480 | + | std::map<int, std::set<int> > atomGroups; |
481 | + | |
482 | + | Molecule::RigidBodyIterator rbIter; |
483 | + | RigidBody* rb; |
484 | + | Molecule::IntegrableObjectIterator ii; |
485 | + | StuntDouble* integrableObject; |
486 | ||
487 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
488 | + | integrableObject = mol->nextIntegrableObject(ii)) { |
489 | + | |
490 | + | if (integrableObject->isRigidBody()) { |
491 | + | rb = static_cast<RigidBody*>(integrableObject); |
492 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
493 | + | std::set<int> rigidAtoms; |
494 | + | for (int i = 0; i < atoms.size(); ++i) { |
495 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
496 | + | } |
497 | + | for (int i = 0; i < atoms.size(); ++i) { |
498 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
499 | + | } |
500 | + | } else { |
501 | + | std::set<int> oneAtomSet; |
502 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
503 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
504 | + | } |
505 | + | } |
506 | + | |
507 | + | |
508 | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
509 | a = bond->getAtomA()->getGlobalIndex(); | |
510 | b = bond->getAtomB()->getGlobalIndex(); | |
# | Line 411 | Line 516 | namespace oopse { | |
516 | b = bend->getAtomB()->getGlobalIndex(); | |
517 | c = bend->getAtomC()->getGlobalIndex(); | |
518 | ||
519 | < | exclude_.removePair(a, b); |
520 | < | exclude_.removePair(a, c); |
521 | < | exclude_.removePair(b, c); |
519 | > | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
520 | > | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
521 | > | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
522 | > | |
523 | > | exclude_.removePairs(rigidSetA, rigidSetB); |
524 | > | exclude_.removePairs(rigidSetA, rigidSetC); |
525 | > | exclude_.removePairs(rigidSetB, rigidSetC); |
526 | > | |
527 | > | //exclude_.removePair(a, b); |
528 | > | //exclude_.removePair(a, c); |
529 | > | //exclude_.removePair(b, c); |
530 | } | |
531 | ||
532 | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
# | Line 421 | Line 534 | namespace oopse { | |
534 | b = torsion->getAtomB()->getGlobalIndex(); | |
535 | c = torsion->getAtomC()->getGlobalIndex(); | |
536 | d = torsion->getAtomD()->getGlobalIndex(); | |
537 | + | |
538 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
539 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
540 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
541 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
542 | + | |
543 | + | exclude_.removePairs(rigidSetA, rigidSetB); |
544 | + | exclude_.removePairs(rigidSetA, rigidSetC); |
545 | + | exclude_.removePairs(rigidSetA, rigidSetD); |
546 | + | exclude_.removePairs(rigidSetB, rigidSetC); |
547 | + | exclude_.removePairs(rigidSetB, rigidSetD); |
548 | + | exclude_.removePairs(rigidSetC, rigidSetD); |
549 | + | |
550 | + | /* |
551 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
552 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
553 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
554 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
555 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
556 | + | exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
557 | ||
558 | + | |
559 | exclude_.removePair(a, b); | |
560 | exclude_.removePair(a, c); | |
561 | exclude_.removePair(a, d); | |
562 | exclude_.removePair(b, c); | |
563 | exclude_.removePair(b, d); | |
564 | exclude_.removePair(c, d); | |
565 | + | */ |
566 | } | |
567 | ||
433 | – | Molecule::RigidBodyIterator rbIter; |
434 | – | RigidBody* rb; |
568 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
569 | std::vector<Atom*> atoms = rb->getAtoms(); | |
570 | for (int i = 0; i < atoms.size() -1 ; ++i) { | |
# | Line 470 | Line 603 | namespace oopse { | |
603 | /** @deprecate */ | |
604 | int isError = 0; | |
605 | ||
606 | + | setupCutoff(); |
607 | + | |
608 | setupElectrostaticSummationMethod( isError ); | |
609 | + | setupSwitchingFunction(); |
610 | + | setupAccumulateBoxDipole(); |
611 | ||
612 | if(isError){ | |
613 | sprintf( painCave.errMsg, | |
# | Line 478 | Line 615 | namespace oopse { | |
615 | painCave.isFatal = 1; | |
616 | simError(); | |
617 | } | |
481 | – | |
482 | – | |
483 | – | setupCutoff(); |
618 | ||
619 | calcNdf(); | |
620 | calcNdfRaw(); | |
# | Line 515 | Line 649 | namespace oopse { | |
649 | int useLennardJones = 0; | |
650 | int useElectrostatic = 0; | |
651 | int useEAM = 0; | |
652 | + | int useSC = 0; |
653 | int useCharge = 0; | |
654 | int useDirectional = 0; | |
655 | int useDipole = 0; | |
# | Line 528 | Line 663 | namespace oopse { | |
663 | //usePBC and useRF are from simParams | |
664 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
665 | int useRF; | |
666 | < | int useDW; |
666 | > | int useSF; |
667 | > | int useSP; |
668 | > | int useBoxDipole; |
669 | std::string myMethod; | |
670 | ||
671 | // set the useRF logical | |
672 | useRF = 0; | |
673 | < | useDW = 0; |
673 | > | useSF = 0; |
674 | > | useSP = 0; |
675 | ||
676 | ||
677 | if (simParams_->haveElectrostaticSummationMethod()) { | |
678 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
679 | toUpper(myMethod); | |
680 | < | if (myMethod == "REACTION_FIELD") { |
681 | < | useRF=1; |
682 | < | } else { |
683 | < | if (myMethod == "DAMPED_WOLF") { |
684 | < | useDW = 1; |
685 | < | } |
680 | > | if (myMethod == "REACTION_FIELD"){ |
681 | > | useRF = 1; |
682 | > | } else if (myMethod == "SHIFTED_FORCE"){ |
683 | > | useSF = 1; |
684 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
685 | > | useSP = 1; |
686 | } | |
687 | } | |
688 | + | |
689 | + | if (simParams_->haveAccumulateBoxDipole()) |
690 | + | if (simParams_->getAccumulateBoxDipole()) |
691 | + | useBoxDipole = 1; |
692 | ||
693 | //loop over all of the atom types | |
694 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
695 | useLennardJones |= (*i)->isLennardJones(); | |
696 | useElectrostatic |= (*i)->isElectrostatic(); | |
697 | useEAM |= (*i)->isEAM(); | |
698 | + | useSC |= (*i)->isSC(); |
699 | useCharge |= (*i)->isCharge(); | |
700 | useDirectional |= (*i)->isDirectional(); | |
701 | useDipole |= (*i)->isDipole(); | |
# | Line 603 | Line 746 | namespace oopse { | |
746 | temp = useEAM; | |
747 | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
748 | ||
749 | + | temp = useSC; |
750 | + | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
751 | + | |
752 | temp = useShape; | |
753 | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
754 | ||
# | Line 612 | Line 758 | namespace oopse { | |
758 | temp = useRF; | |
759 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
760 | ||
761 | < | temp = useDW; |
762 | < | MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
761 | > | temp = useSF; |
762 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
763 | ||
764 | + | temp = useSP; |
765 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
766 | + | |
767 | + | temp = useBoxDipole; |
768 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
769 | + | |
770 | #endif | |
771 | ||
772 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 627 | Line 779 | namespace oopse { | |
779 | fInfo_.SIM_uses_StickyPower = useStickyPower; | |
780 | fInfo_.SIM_uses_GayBerne = useGayBerne; | |
781 | fInfo_.SIM_uses_EAM = useEAM; | |
782 | + | fInfo_.SIM_uses_SC = useSC; |
783 | fInfo_.SIM_uses_Shapes = useShape; | |
784 | fInfo_.SIM_uses_FLARB = useFLARB; | |
785 | fInfo_.SIM_uses_RF = useRF; | |
786 | < | fInfo_.SIM_uses_DampedWolf = useDW; |
787 | < | |
788 | < | if( myMethod == "REACTION_FIELD") { |
636 | < | |
637 | < | if (simParams_->haveDielectric()) { |
638 | < | fInfo_.dielect = simParams_->getDielectric(); |
639 | < | } else { |
640 | < | sprintf(painCave.errMsg, |
641 | < | "SimSetup Error: No Dielectric constant was set.\n" |
642 | < | "\tYou are trying to use Reaction Field without" |
643 | < | "\tsetting a dielectric constant!\n"); |
644 | < | painCave.isFatal = 1; |
645 | < | simError(); |
646 | < | } |
647 | < | } |
648 | < | |
786 | > | fInfo_.SIM_uses_SF = useSF; |
787 | > | fInfo_.SIM_uses_SP = useSP; |
788 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
789 | } | |
790 | ||
791 | void SimInfo::setupFortranSim() { | |
# | Line 662 | Line 802 | namespace oopse { | |
802 | } | |
803 | ||
804 | //calculate mass ratio of cutoff group | |
805 | < | std::vector<double> mfact; |
805 | > | std::vector<RealType> mfact; |
806 | SimInfo::MoleculeIterator mi; | |
807 | Molecule* mol; | |
808 | Molecule::CutoffGroupIterator ci; | |
809 | CutoffGroup* cg; | |
810 | Molecule::AtomIterator ai; | |
811 | Atom* atom; | |
812 | < | double totalMass; |
812 | > | RealType totalMass; |
813 | ||
814 | //to avoid memory reallocation, reserve enough space for mfact | |
815 | mfact.reserve(getNCutoffGroups()); | |
# | Line 730 | Line 870 | namespace oopse { | |
870 | "succesfully sent the simulation information to fortran.\n"); | |
871 | MPIcheckPoint(); | |
872 | #endif // is_mpi | |
873 | + | |
874 | + | // Setup number of neighbors in neighbor list if present |
875 | + | if (simParams_->haveNeighborListNeighbors()) { |
876 | + | setNeighbors(simParams_->getNeighborListNeighbors()); |
877 | + | } |
878 | + | |
879 | + | |
880 | } | |
881 | ||
882 | ||
# | Line 792 | Line 939 | namespace oopse { | |
939 | ||
940 | #endif | |
941 | ||
942 | < | double SimInfo::calcMaxCutoffRadius() { |
942 | > | void SimInfo::setupCutoff() { |
943 | > | |
944 | > | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
945 | ||
946 | + | // Check the cutoff policy |
947 | + | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default |
948 | ||
949 | < | std::set<AtomType*> atomTypes; |
950 | < | std::set<AtomType*>::iterator i; |
951 | < | std::vector<double> cutoffRadius; |
952 | < | |
953 | < | //get the unique atom types |
803 | < | atomTypes = getUniqueAtomTypes(); |
804 | < | |
805 | < | //query the max cutoff radius among these atom types |
806 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
807 | < | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); |
949 | > | std::string myPolicy; |
950 | > | if (forceFieldOptions_.haveCutoffPolicy()){ |
951 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); |
952 | > | }else if (simParams_->haveCutoffPolicy()) { |
953 | > | myPolicy = simParams_->getCutoffPolicy(); |
954 | } | |
955 | ||
956 | < | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); |
811 | < | #ifdef IS_MPI |
812 | < | //pick the max cutoff radius among the processors |
813 | < | #endif |
814 | < | |
815 | < | return maxCutoffRadius; |
816 | < | } |
817 | < | |
818 | < | void SimInfo::getCutoff(double& rcut, double& rsw) { |
819 | < | |
820 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
821 | < | |
822 | < | if (!simParams_->haveCutoffRadius()){ |
823 | < | sprintf(painCave.errMsg, |
824 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
825 | < | "\tOOPSE will use a default value of 15.0 angstroms" |
826 | < | "\tfor the cutoffRadius.\n"); |
827 | < | painCave.isFatal = 0; |
828 | < | simError(); |
829 | < | rcut = 15.0; |
830 | < | } else{ |
831 | < | rcut = simParams_->getCutoffRadius(); |
832 | < | } |
833 | < | |
834 | < | if (!simParams_->haveSwitchingRadius()){ |
835 | < | sprintf(painCave.errMsg, |
836 | < | "SimCreator Warning: No value was set for switchingRadius.\n" |
837 | < | "\tOOPSE will use a default value of\n" |
838 | < | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
839 | < | painCave.isFatal = 0; |
840 | < | simError(); |
841 | < | rsw = 0.85 * rcut; |
842 | < | } else{ |
843 | < | rsw = simParams_->getSwitchingRadius(); |
844 | < | } |
845 | < | |
846 | < | } else { |
847 | < | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in |
848 | < | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used |
849 | < | |
850 | < | if (simParams_->haveCutoffRadius()) { |
851 | < | rcut = simParams_->getCutoffRadius(); |
852 | < | } else { |
853 | < | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system |
854 | < | rcut = calcMaxCutoffRadius(); |
855 | < | } |
856 | < | |
857 | < | if (simParams_->haveSwitchingRadius()) { |
858 | < | rsw = simParams_->getSwitchingRadius(); |
859 | < | } else { |
860 | < | rsw = rcut; |
861 | < | } |
862 | < | |
863 | < | } |
864 | < | } |
865 | < | |
866 | < | void SimInfo::setupCutoff() { |
867 | < | getCutoff(rcut_, rsw_); |
868 | < | double rnblist = rcut_ + 1; // skin of neighbor list |
869 | < | |
870 | < | //Pass these cutoff radius etc. to fortran. This function should be called once and only once |
871 | < | |
872 | < | int cp = TRADITIONAL_CUTOFF_POLICY; |
873 | < | if (simParams_->haveCutoffPolicy()) { |
874 | < | std::string myPolicy = simParams_->getCutoffPolicy(); |
956 | > | if (!myPolicy.empty()){ |
957 | toUpper(myPolicy); | |
958 | if (myPolicy == "MIX") { | |
959 | cp = MIX_CUTOFF_POLICY; | |
# | Line 890 | Line 972 | namespace oopse { | |
972 | } | |
973 | } | |
974 | } | |
975 | < | } |
976 | < | |
975 | > | } |
976 | > | notifyFortranCutoffPolicy(&cp); |
977 | ||
978 | + | // Check the Skin Thickness for neighborlists |
979 | + | RealType skin; |
980 | if (simParams_->haveSkinThickness()) { | |
981 | < | double skinThickness = simParams_->getSkinThickness(); |
982 | < | } |
981 | > | skin = simParams_->getSkinThickness(); |
982 | > | notifyFortranSkinThickness(&skin); |
983 | > | } |
984 | > | |
985 | > | // Check if the cutoff was set explicitly: |
986 | > | if (simParams_->haveCutoffRadius()) { |
987 | > | rcut_ = simParams_->getCutoffRadius(); |
988 | > | if (simParams_->haveSwitchingRadius()) { |
989 | > | rsw_ = simParams_->getSwitchingRadius(); |
990 | > | } else { |
991 | > | if (fInfo_.SIM_uses_Charges | |
992 | > | fInfo_.SIM_uses_Dipoles | |
993 | > | fInfo_.SIM_uses_RF) { |
994 | > | |
995 | > | rsw_ = 0.85 * rcut_; |
996 | > | sprintf(painCave.errMsg, |
997 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
998 | > | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" |
999 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1000 | > | painCave.isFatal = 0; |
1001 | > | simError(); |
1002 | > | } else { |
1003 | > | rsw_ = rcut_; |
1004 | > | sprintf(painCave.errMsg, |
1005 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1006 | > | "\tOOPSE will use the same value as the cutoffRadius.\n" |
1007 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1008 | > | painCave.isFatal = 0; |
1009 | > | simError(); |
1010 | > | } |
1011 | > | } |
1012 | > | |
1013 | > | notifyFortranCutoffs(&rcut_, &rsw_); |
1014 | > | |
1015 | > | } else { |
1016 | > | |
1017 | > | // For electrostatic atoms, we'll assume a large safe value: |
1018 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
1019 | > | sprintf(painCave.errMsg, |
1020 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
1021 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
1022 | > | "\tfor the cutoffRadius.\n"); |
1023 | > | painCave.isFatal = 0; |
1024 | > | simError(); |
1025 | > | rcut_ = 15.0; |
1026 | > | |
1027 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
1028 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1029 | > | toUpper(myMethod); |
1030 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1031 | > | if (simParams_->haveSwitchingRadius()){ |
1032 | > | sprintf(painCave.errMsg, |
1033 | > | "SimInfo Warning: A value was set for the switchingRadius\n" |
1034 | > | "\teven though the electrostaticSummationMethod was\n" |
1035 | > | "\tset to %s\n", myMethod.c_str()); |
1036 | > | painCave.isFatal = 1; |
1037 | > | simError(); |
1038 | > | } |
1039 | > | } |
1040 | > | } |
1041 | > | |
1042 | > | if (simParams_->haveSwitchingRadius()){ |
1043 | > | rsw_ = simParams_->getSwitchingRadius(); |
1044 | > | } else { |
1045 | > | sprintf(painCave.errMsg, |
1046 | > | "SimCreator Warning: No value was set for switchingRadius.\n" |
1047 | > | "\tOOPSE will use a default value of\n" |
1048 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
1049 | > | painCave.isFatal = 0; |
1050 | > | simError(); |
1051 | > | rsw_ = 0.85 * rcut_; |
1052 | > | } |
1053 | > | notifyFortranCutoffs(&rcut_, &rsw_); |
1054 | > | } else { |
1055 | > | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so |
1056 | > | // We'll punt and let fortran figure out the cutoffs later. |
1057 | > | |
1058 | > | notifyFortranYouAreOnYourOwn(); |
1059 | ||
1060 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
1061 | < | // also send cutoff notification to electrostatics |
902 | < | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
1060 | > | } |
1061 | > | } |
1062 | } | |
1063 | ||
1064 | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | |
1065 | ||
1066 | int errorOut; | |
1067 | int esm = NONE; | |
1068 | < | double alphaVal; |
1069 | < | double dielectric; |
1070 | < | |
1068 | > | int sm = UNDAMPED; |
1069 | > | RealType alphaVal; |
1070 | > | RealType dielectric; |
1071 | > | |
1072 | errorOut = isError; | |
913 | – | alphaVal = simParams_->getDampingAlpha(); |
914 | – | dielectric = simParams_->getDielectric(); |
1073 | ||
1074 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1075 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
# | Line 919 | Line 1077 | namespace oopse { | |
1077 | if (myMethod == "NONE") { | |
1078 | esm = NONE; | |
1079 | } else { | |
1080 | < | if (myMethod == "UNDAMPED_WOLF") { |
1081 | < | esm = UNDAMPED_WOLF; |
1080 | > | if (myMethod == "SWITCHING_FUNCTION") { |
1081 | > | esm = SWITCHING_FUNCTION; |
1082 | } else { | |
1083 | < | if (myMethod == "DAMPED_WOLF") { |
1084 | < | esm = DAMPED_WOLF; |
1085 | < | if (!simParams_->haveDampingAlpha()) { |
1086 | < | //throw error |
1087 | < | sprintf( painCave.errMsg, |
930 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Damped Wolf Method.", alphaVal); |
931 | < | painCave.isFatal = 0; |
932 | < | simError(); |
933 | < | } |
934 | < | } else { |
935 | < | if (myMethod == "REACTION_FIELD") { |
936 | < | esm = REACTION_FIELD; |
1083 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
1084 | > | esm = SHIFTED_POTENTIAL; |
1085 | > | } else { |
1086 | > | if (myMethod == "SHIFTED_FORCE") { |
1087 | > | esm = SHIFTED_FORCE; |
1088 | } else { | |
1089 | < | // throw error |
1090 | < | sprintf( painCave.errMsg, |
1091 | < | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"undamped_wolf\", \"damped_wolf\", or \"reaction_field\".", myMethod.c_str() ); |
1092 | < | painCave.isFatal = 1; |
1093 | < | simError(); |
1094 | < | } |
1095 | < | } |
1089 | > | if (myMethod == "REACTION_FIELD") { |
1090 | > | esm = REACTION_FIELD; |
1091 | > | dielectric = simParams_->getDielectric(); |
1092 | > | if (!simParams_->haveDielectric()) { |
1093 | > | // throw warning |
1094 | > | sprintf( painCave.errMsg, |
1095 | > | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1096 | > | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1097 | > | painCave.isFatal = 0; |
1098 | > | simError(); |
1099 | > | } |
1100 | > | } else { |
1101 | > | // throw error |
1102 | > | sprintf( painCave.errMsg, |
1103 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" |
1104 | > | "\t(Input file specified %s .)\n" |
1105 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" |
1106 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" |
1107 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); |
1108 | > | painCave.isFatal = 1; |
1109 | > | simError(); |
1110 | > | } |
1111 | > | } |
1112 | > | } |
1113 | > | } |
1114 | > | } |
1115 | > | } |
1116 | > | |
1117 | > | if (simParams_->haveElectrostaticScreeningMethod()) { |
1118 | > | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
1119 | > | toUpper(myScreen); |
1120 | > | if (myScreen == "UNDAMPED") { |
1121 | > | sm = UNDAMPED; |
1122 | > | } else { |
1123 | > | if (myScreen == "DAMPED") { |
1124 | > | sm = DAMPED; |
1125 | > | if (!simParams_->haveDampingAlpha()) { |
1126 | > | // first set a cutoff dependent alpha value |
1127 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1128 | > | alphaVal = 0.5125 - rcut_* 0.025; |
1129 | > | // for values rcut > 20.5, alpha is zero |
1130 | > | if (alphaVal < 0) alphaVal = 0; |
1131 | > | |
1132 | > | // throw warning |
1133 | > | sprintf( painCave.errMsg, |
1134 | > | "SimInfo warning: dampingAlpha was not specified in the input file.\n" |
1135 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1136 | > | painCave.isFatal = 0; |
1137 | > | simError(); |
1138 | > | } else { |
1139 | > | alphaVal = simParams_->getDampingAlpha(); |
1140 | > | } |
1141 | > | |
1142 | > | } else { |
1143 | > | // throw error |
1144 | > | sprintf( painCave.errMsg, |
1145 | > | "SimInfo error: Unknown electrostaticScreeningMethod.\n" |
1146 | > | "\t(Input file specified %s .)\n" |
1147 | > | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" |
1148 | > | "or \"damped\".\n", myScreen.c_str() ); |
1149 | > | painCave.isFatal = 1; |
1150 | > | simError(); |
1151 | } | |
1152 | } | |
1153 | } | |
1154 | + | |
1155 | // let's pass some summation method variables to fortran | |
1156 | setElectrostaticSummationMethod( &esm ); | |
1157 | < | setDampedWolfAlpha( &alphaVal ); |
1157 | > | setFortranElectrostaticMethod( &esm ); |
1158 | > | setScreeningMethod( &sm ); |
1159 | > | setDampingAlpha( &alphaVal ); |
1160 | setReactionFieldDielectric( &dielectric ); | |
1161 | < | initFortranFF( &esm, &errorOut ); |
1161 | > | initFortranFF( &errorOut ); |
1162 | } | |
1163 | ||
1164 | + | void SimInfo::setupSwitchingFunction() { |
1165 | + | int ft = CUBIC; |
1166 | + | |
1167 | + | if (simParams_->haveSwitchingFunctionType()) { |
1168 | + | std::string funcType = simParams_->getSwitchingFunctionType(); |
1169 | + | toUpper(funcType); |
1170 | + | if (funcType == "CUBIC") { |
1171 | + | ft = CUBIC; |
1172 | + | } else { |
1173 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1174 | + | ft = FIFTH_ORDER_POLY; |
1175 | + | } else { |
1176 | + | // throw error |
1177 | + | sprintf( painCave.errMsg, |
1178 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1179 | + | painCave.isFatal = 1; |
1180 | + | simError(); |
1181 | + | } |
1182 | + | } |
1183 | + | } |
1184 | + | |
1185 | + | // send switching function notification to switcheroo |
1186 | + | setFunctionType(&ft); |
1187 | + | |
1188 | + | } |
1189 | + | |
1190 | + | void SimInfo::setupAccumulateBoxDipole() { |
1191 | + | |
1192 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1193 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1194 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1195 | + | setAccumulateBoxDipole(); |
1196 | + | calcBoxDipole_ = true; |
1197 | + | } |
1198 | + | |
1199 | + | } |
1200 | + | |
1201 | void SimInfo::addProperty(GenericData* genData) { | |
1202 | properties_.addProperty(genData); | |
1203 | } | |
# | Line 1008 | Line 1254 | namespace oopse { | |
1254 | Molecule* mol; | |
1255 | ||
1256 | Vector3d comVel(0.0); | |
1257 | < | double totalMass = 0.0; |
1257 | > | RealType totalMass = 0.0; |
1258 | ||
1259 | ||
1260 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1261 | < | double mass = mol->getMass(); |
1261 | > | RealType mass = mol->getMass(); |
1262 | totalMass += mass; | |
1263 | comVel += mass * mol->getComVel(); | |
1264 | } | |
1265 | ||
1266 | #ifdef IS_MPI | |
1267 | < | double tmpMass = totalMass; |
1267 | > | RealType tmpMass = totalMass; |
1268 | Vector3d tmpComVel(comVel); | |
1269 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1270 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1269 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1270 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1271 | #endif | |
1272 | ||
1273 | comVel /= totalMass; | |
# | Line 1034 | Line 1280 | namespace oopse { | |
1280 | Molecule* mol; | |
1281 | ||
1282 | Vector3d com(0.0); | |
1283 | < | double totalMass = 0.0; |
1283 | > | RealType totalMass = 0.0; |
1284 | ||
1285 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1286 | < | double mass = mol->getMass(); |
1286 | > | RealType mass = mol->getMass(); |
1287 | totalMass += mass; | |
1288 | com += mass * mol->getCom(); | |
1289 | } | |
1290 | ||
1291 | #ifdef IS_MPI | |
1292 | < | double tmpMass = totalMass; |
1292 | > | RealType tmpMass = totalMass; |
1293 | Vector3d tmpCom(com); | |
1294 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1295 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1294 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1295 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1296 | #endif | |
1297 | ||
1298 | com /= totalMass; | |
# | Line 1070 | Line 1316 | namespace oopse { | |
1316 | Molecule* mol; | |
1317 | ||
1318 | ||
1319 | < | double totalMass = 0.0; |
1319 | > | RealType totalMass = 0.0; |
1320 | ||
1321 | ||
1322 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1323 | < | double mass = mol->getMass(); |
1323 | > | RealType mass = mol->getMass(); |
1324 | totalMass += mass; | |
1325 | com += mass * mol->getCom(); | |
1326 | comVel += mass * mol->getComVel(); | |
1327 | } | |
1328 | ||
1329 | #ifdef IS_MPI | |
1330 | < | double tmpMass = totalMass; |
1330 | > | RealType tmpMass = totalMass; |
1331 | Vector3d tmpCom(com); | |
1332 | Vector3d tmpComVel(comVel); | |
1333 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1334 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1335 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1333 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1334 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1335 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1336 | #endif | |
1337 | ||
1338 | com /= totalMass; | |
# | Line 1105 | Line 1351 | namespace oopse { | |
1351 | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | |
1352 | ||
1353 | ||
1354 | < | double xx = 0.0; |
1355 | < | double yy = 0.0; |
1356 | < | double zz = 0.0; |
1357 | < | double xy = 0.0; |
1358 | < | double xz = 0.0; |
1359 | < | double yz = 0.0; |
1354 | > | RealType xx = 0.0; |
1355 | > | RealType yy = 0.0; |
1356 | > | RealType zz = 0.0; |
1357 | > | RealType xy = 0.0; |
1358 | > | RealType xz = 0.0; |
1359 | > | RealType yz = 0.0; |
1360 | Vector3d com(0.0); | |
1361 | Vector3d comVel(0.0); | |
1362 | ||
# | Line 1122 | Line 1368 | namespace oopse { | |
1368 | Vector3d thisq(0.0); | |
1369 | Vector3d thisv(0.0); | |
1370 | ||
1371 | < | double thisMass = 0.0; |
1371 | > | RealType thisMass = 0.0; |
1372 | ||
1373 | ||
1374 | ||
# | Line 1160 | Line 1406 | namespace oopse { | |
1406 | #ifdef IS_MPI | |
1407 | Mat3x3d tmpI(inertiaTensor); | |
1408 | Vector3d tmpAngMom; | |
1409 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1410 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1409 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1410 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1411 | #endif | |
1412 | ||
1413 | return; | |
# | Line 1182 | Line 1428 | namespace oopse { | |
1428 | Vector3d thisr(0.0); | |
1429 | Vector3d thisp(0.0); | |
1430 | ||
1431 | < | double thisMass; |
1431 | > | RealType thisMass; |
1432 | ||
1433 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1434 | thisMass = mol->getMass(); | |
# | Line 1195 | Line 1441 | namespace oopse { | |
1441 | ||
1442 | #ifdef IS_MPI | |
1443 | Vector3d tmpAngMom; | |
1444 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1444 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1445 | #endif | |
1446 | ||
1447 | return angularMomentum; | |
1448 | } | |
1449 | ||
1450 | < | |
1450 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { |
1451 | > | return IOIndexToIntegrableObject.at(index); |
1452 | > | } |
1453 | > | |
1454 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1455 | > | IOIndexToIntegrableObject= v; |
1456 | > | } |
1457 | > | |
1458 | > | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1459 | > | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1460 | > | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1461 | > | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1462 | > | */ |
1463 | > | void SimInfo::getGyrationalVolume(RealType &volume){ |
1464 | > | Mat3x3d intTensor; |
1465 | > | RealType det; |
1466 | > | Vector3d dummyAngMom; |
1467 | > | RealType sysconstants; |
1468 | > | RealType geomCnst; |
1469 | > | |
1470 | > | geomCnst = 3.0/2.0; |
1471 | > | /* Get the inertial tensor and angular momentum for free*/ |
1472 | > | getInertiaTensor(intTensor,dummyAngMom); |
1473 | > | |
1474 | > | det = intTensor.determinant(); |
1475 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1476 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1477 | > | return; |
1478 | > | } |
1479 | > | |
1480 | > | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1481 | > | Mat3x3d intTensor; |
1482 | > | Vector3d dummyAngMom; |
1483 | > | RealType sysconstants; |
1484 | > | RealType geomCnst; |
1485 | > | |
1486 | > | geomCnst = 3.0/2.0; |
1487 | > | /* Get the inertial tensor and angular momentum for free*/ |
1488 | > | getInertiaTensor(intTensor,dummyAngMom); |
1489 | > | |
1490 | > | detI = intTensor.determinant(); |
1491 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1492 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1493 | > | return; |
1494 | > | } |
1495 | > | /* |
1496 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1497 | > | assert( v.size() == nAtoms_ + nRigidBodies_); |
1498 | > | sdByGlobalIndex_ = v; |
1499 | > | } |
1500 | > | |
1501 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1502 | > | //assert(index < nAtoms_ + nRigidBodies_); |
1503 | > | return sdByGlobalIndex_.at(index); |
1504 | > | } |
1505 | > | */ |
1506 | }//end namespace oopse | |
1507 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |