# | 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 | > | useAtomicVirial_(true) { |
97 | ||
82 | – | |
83 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; |
98 | MoleculeStamp* molStamp; | |
99 | int nMolWithSameStamp; | |
100 | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | |
# | Line 88 | Line 102 | namespace oopse { | |
102 | CutoffGroupStamp* cgStamp; | |
103 | RigidBodyStamp* rbStamp; | |
104 | int nRigidAtoms = 0; | |
105 | < | |
106 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { |
107 | < | molStamp = i->first; |
108 | < | nMolWithSameStamp = i->second; |
105 | > | std::vector<Component*> components = simParams->getComponents(); |
106 | > | |
107 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
108 | > | molStamp = (*i)->getMoleculeStamp(); |
109 | > | nMolWithSameStamp = (*i)->getNMol(); |
110 | ||
111 | addMoleculeStamp(molStamp, nMolWithSameStamp); | |
112 | ||
113 | //calculate atoms in molecules | |
114 | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | |
115 | ||
101 | – | |
116 | //calculate atoms in cutoff groups | |
117 | int nAtomsInGroups = 0; | |
118 | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | |
119 | ||
120 | for (int j=0; j < nCutoffGroupsInStamp; j++) { | |
121 | < | cgStamp = molStamp->getCutoffGroup(j); |
121 | > | cgStamp = molStamp->getCutoffGroupStamp(j); |
122 | nAtomsInGroups += cgStamp->getNMembers(); | |
123 | } | |
124 | ||
# | Line 117 | Line 131 | namespace oopse { | |
131 | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | |
132 | ||
133 | for (int j=0; j < nRigidBodiesInStamp; j++) { | |
134 | < | rbStamp = molStamp->getRigidBody(j); |
134 | > | rbStamp = molStamp->getRigidBodyStamp(j); |
135 | nAtomsInRigidBodies += rbStamp->getNMembers(); | |
136 | } | |
137 | ||
# | Line 156 | Line 170 | namespace oopse { | |
170 | } | |
171 | molecules_.clear(); | |
172 | ||
159 | – | delete stamps_; |
173 | delete sman_; | |
174 | delete simParams_; | |
175 | delete forceField_; | |
# | Line 263 | Line 276 | namespace oopse { | |
276 | } | |
277 | } | |
278 | ||
279 | < | }//end for (integrableObject) |
280 | < | }// end for (mol) |
279 | > | } |
280 | > | } |
281 | ||
282 | // n_constraints is local, so subtract them on each processor | |
283 | ndf_local -= nConstraints_; | |
# | Line 281 | Line 294 | namespace oopse { | |
294 | ||
295 | } | |
296 | ||
297 | + | int SimInfo::getFdf() { |
298 | + | #ifdef IS_MPI |
299 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
300 | + | #else |
301 | + | fdf_ = fdf_local; |
302 | + | #endif |
303 | + | return fdf_; |
304 | + | } |
305 | + | |
306 | void SimInfo::calcNdfRaw() { | |
307 | int ndfRaw_local; | |
308 | ||
# | Line 343 | Line 365 | namespace oopse { | |
365 | int b; | |
366 | int c; | |
367 | int d; | |
368 | + | |
369 | + | std::map<int, std::set<int> > atomGroups; |
370 | + | |
371 | + | Molecule::RigidBodyIterator rbIter; |
372 | + | RigidBody* rb; |
373 | + | Molecule::IntegrableObjectIterator ii; |
374 | + | StuntDouble* integrableObject; |
375 | ||
376 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
377 | + | integrableObject = mol->nextIntegrableObject(ii)) { |
378 | + | |
379 | + | if (integrableObject->isRigidBody()) { |
380 | + | rb = static_cast<RigidBody*>(integrableObject); |
381 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
382 | + | std::set<int> rigidAtoms; |
383 | + | for (int i = 0; i < atoms.size(); ++i) { |
384 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
385 | + | } |
386 | + | for (int i = 0; i < atoms.size(); ++i) { |
387 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
388 | + | } |
389 | + | } else { |
390 | + | std::set<int> oneAtomSet; |
391 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
392 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
393 | + | } |
394 | + | } |
395 | + | |
396 | + | |
397 | + | |
398 | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
399 | a = bond->getAtomA()->getGlobalIndex(); | |
400 | b = bond->getAtomB()->getGlobalIndex(); | |
# | Line 354 | Line 405 | namespace oopse { | |
405 | a = bend->getAtomA()->getGlobalIndex(); | |
406 | b = bend->getAtomB()->getGlobalIndex(); | |
407 | c = bend->getAtomC()->getGlobalIndex(); | |
408 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
409 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
410 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
411 | ||
412 | < | exclude_.addPair(a, b); |
413 | < | exclude_.addPair(a, c); |
414 | < | exclude_.addPair(b, c); |
412 | > | exclude_.addPairs(rigidSetA, rigidSetB); |
413 | > | exclude_.addPairs(rigidSetA, rigidSetC); |
414 | > | exclude_.addPairs(rigidSetB, rigidSetC); |
415 | > | |
416 | > | //exclude_.addPair(a, b); |
417 | > | //exclude_.addPair(a, c); |
418 | > | //exclude_.addPair(b, c); |
419 | } | |
420 | ||
421 | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
# | Line 365 | Line 423 | namespace oopse { | |
423 | b = torsion->getAtomB()->getGlobalIndex(); | |
424 | c = torsion->getAtomC()->getGlobalIndex(); | |
425 | d = torsion->getAtomD()->getGlobalIndex(); | |
426 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
427 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
428 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
429 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
430 | ||
431 | + | exclude_.addPairs(rigidSetA, rigidSetB); |
432 | + | exclude_.addPairs(rigidSetA, rigidSetC); |
433 | + | exclude_.addPairs(rigidSetA, rigidSetD); |
434 | + | exclude_.addPairs(rigidSetB, rigidSetC); |
435 | + | exclude_.addPairs(rigidSetB, rigidSetD); |
436 | + | exclude_.addPairs(rigidSetC, rigidSetD); |
437 | + | |
438 | + | /* |
439 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
440 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
441 | + | exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
442 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
443 | + | exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
444 | + | exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
445 | + | |
446 | + | |
447 | exclude_.addPair(a, b); | |
448 | exclude_.addPair(a, c); | |
449 | exclude_.addPair(a, d); | |
450 | exclude_.addPair(b, c); | |
451 | exclude_.addPair(b, d); | |
452 | exclude_.addPair(c, d); | |
453 | + | */ |
454 | } | |
455 | ||
377 | – | Molecule::RigidBodyIterator rbIter; |
378 | – | RigidBody* rb; |
456 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
457 | std::vector<Atom*> atoms = rb->getAtoms(); | |
458 | for (int i = 0; i < atoms.size() -1 ; ++i) { | |
# | Line 400 | Line 477 | namespace oopse { | |
477 | int b; | |
478 | int c; | |
479 | int d; | |
480 | + | |
481 | + | std::map<int, std::set<int> > atomGroups; |
482 | + | |
483 | + | Molecule::RigidBodyIterator rbIter; |
484 | + | RigidBody* rb; |
485 | + | Molecule::IntegrableObjectIterator ii; |
486 | + | StuntDouble* integrableObject; |
487 | ||
488 | + | for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; |
489 | + | integrableObject = mol->nextIntegrableObject(ii)) { |
490 | + | |
491 | + | if (integrableObject->isRigidBody()) { |
492 | + | rb = static_cast<RigidBody*>(integrableObject); |
493 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
494 | + | std::set<int> rigidAtoms; |
495 | + | for (int i = 0; i < atoms.size(); ++i) { |
496 | + | rigidAtoms.insert(atoms[i]->getGlobalIndex()); |
497 | + | } |
498 | + | for (int i = 0; i < atoms.size(); ++i) { |
499 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
500 | + | } |
501 | + | } else { |
502 | + | std::set<int> oneAtomSet; |
503 | + | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
504 | + | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
505 | + | } |
506 | + | } |
507 | + | |
508 | + | |
509 | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { | |
510 | a = bond->getAtomA()->getGlobalIndex(); | |
511 | b = bond->getAtomB()->getGlobalIndex(); | |
# | Line 412 | Line 517 | namespace oopse { | |
517 | b = bend->getAtomB()->getGlobalIndex(); | |
518 | c = bend->getAtomC()->getGlobalIndex(); | |
519 | ||
520 | < | exclude_.removePair(a, b); |
521 | < | exclude_.removePair(a, c); |
522 | < | exclude_.removePair(b, c); |
520 | > | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
521 | > | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
522 | > | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
523 | > | |
524 | > | exclude_.removePairs(rigidSetA, rigidSetB); |
525 | > | exclude_.removePairs(rigidSetA, rigidSetC); |
526 | > | exclude_.removePairs(rigidSetB, rigidSetC); |
527 | > | |
528 | > | //exclude_.removePair(a, b); |
529 | > | //exclude_.removePair(a, c); |
530 | > | //exclude_.removePair(b, c); |
531 | } | |
532 | ||
533 | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { | |
# | Line 422 | Line 535 | namespace oopse { | |
535 | b = torsion->getAtomB()->getGlobalIndex(); | |
536 | c = torsion->getAtomC()->getGlobalIndex(); | |
537 | d = torsion->getAtomD()->getGlobalIndex(); | |
538 | + | |
539 | + | std::set<int> rigidSetA = getRigidSet(a, atomGroups); |
540 | + | std::set<int> rigidSetB = getRigidSet(b, atomGroups); |
541 | + | std::set<int> rigidSetC = getRigidSet(c, atomGroups); |
542 | + | std::set<int> rigidSetD = getRigidSet(d, atomGroups); |
543 | + | |
544 | + | exclude_.removePairs(rigidSetA, rigidSetB); |
545 | + | exclude_.removePairs(rigidSetA, rigidSetC); |
546 | + | exclude_.removePairs(rigidSetA, rigidSetD); |
547 | + | exclude_.removePairs(rigidSetB, rigidSetC); |
548 | + | exclude_.removePairs(rigidSetB, rigidSetD); |
549 | + | exclude_.removePairs(rigidSetC, rigidSetD); |
550 | ||
551 | + | /* |
552 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); |
553 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); |
554 | + | exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); |
555 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); |
556 | + | exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); |
557 | + | exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); |
558 | + | |
559 | + | |
560 | exclude_.removePair(a, b); | |
561 | exclude_.removePair(a, c); | |
562 | exclude_.removePair(a, d); | |
563 | exclude_.removePair(b, c); | |
564 | exclude_.removePair(b, d); | |
565 | exclude_.removePair(c, d); | |
566 | + | */ |
567 | } | |
568 | ||
434 | – | Molecule::RigidBodyIterator rbIter; |
435 | – | RigidBody* rb; |
569 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
570 | std::vector<Atom*> atoms = rb->getAtoms(); | |
571 | for (int i = 0; i < atoms.size() -1 ; ++i) { | |
# | Line 471 | Line 604 | namespace oopse { | |
604 | /** @deprecate */ | |
605 | int isError = 0; | |
606 | ||
607 | + | setupCutoff(); |
608 | + | |
609 | setupElectrostaticSummationMethod( isError ); | |
610 | + | setupSwitchingFunction(); |
611 | + | setupAccumulateBoxDipole(); |
612 | ||
613 | if(isError){ | |
614 | sprintf( painCave.errMsg, | |
# | Line 479 | Line 616 | namespace oopse { | |
616 | painCave.isFatal = 1; | |
617 | simError(); | |
618 | } | |
482 | – | |
483 | – | |
484 | – | setupCutoff(); |
619 | ||
620 | calcNdf(); | |
621 | calcNdfRaw(); | |
# | Line 516 | Line 650 | namespace oopse { | |
650 | int useLennardJones = 0; | |
651 | int useElectrostatic = 0; | |
652 | int useEAM = 0; | |
653 | + | int useSC = 0; |
654 | int useCharge = 0; | |
655 | int useDirectional = 0; | |
656 | int useDipole = 0; | |
# | Line 530 | Line 665 | namespace oopse { | |
665 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
666 | int useRF; | |
667 | int useSF; | |
668 | + | int useSP; |
669 | + | int useBoxDipole; |
670 | + | |
671 | std::string myMethod; | |
672 | ||
673 | // set the useRF logical | |
674 | useRF = 0; | |
675 | useSF = 0; | |
676 | + | useSP = 0; |
677 | ||
678 | ||
679 | if (simParams_->haveElectrostaticSummationMethod()) { | |
680 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
681 | toUpper(myMethod); | |
682 | < | if (myMethod == "REACTION_FIELD") { |
683 | < | useRF=1; |
684 | < | } else { |
685 | < | if (myMethod == "SHIFTED_FORCE") { |
686 | < | useSF = 1; |
687 | < | } |
682 | > | if (myMethod == "REACTION_FIELD"){ |
683 | > | useRF = 1; |
684 | > | } else if (myMethod == "SHIFTED_FORCE"){ |
685 | > | useSF = 1; |
686 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
687 | > | useSP = 1; |
688 | } | |
689 | } | |
690 | + | |
691 | + | if (simParams_->haveAccumulateBoxDipole()) |
692 | + | if (simParams_->getAccumulateBoxDipole()) |
693 | + | useBoxDipole = 1; |
694 | ||
695 | + | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
696 | + | |
697 | //loop over all of the atom types | |
698 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
699 | useLennardJones |= (*i)->isLennardJones(); | |
700 | useElectrostatic |= (*i)->isElectrostatic(); | |
701 | useEAM |= (*i)->isEAM(); | |
702 | + | useSC |= (*i)->isSC(); |
703 | useCharge |= (*i)->isCharge(); | |
704 | useDirectional |= (*i)->isDirectional(); | |
705 | useDipole |= (*i)->isDipole(); | |
# | Line 604 | Line 750 | namespace oopse { | |
750 | temp = useEAM; | |
751 | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
752 | ||
753 | + | temp = useSC; |
754 | + | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
755 | + | |
756 | temp = useShape; | |
757 | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
758 | ||
# | Line 614 | Line 763 | namespace oopse { | |
763 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
764 | ||
765 | temp = useSF; | |
766 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
766 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
767 | ||
768 | + | temp = useSP; |
769 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
770 | + | |
771 | + | temp = useBoxDipole; |
772 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
773 | + | |
774 | + | temp = useAtomicVirial_; |
775 | + | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
776 | + | |
777 | #endif | |
778 | ||
779 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 628 | Line 786 | namespace oopse { | |
786 | fInfo_.SIM_uses_StickyPower = useStickyPower; | |
787 | fInfo_.SIM_uses_GayBerne = useGayBerne; | |
788 | fInfo_.SIM_uses_EAM = useEAM; | |
789 | + | fInfo_.SIM_uses_SC = useSC; |
790 | fInfo_.SIM_uses_Shapes = useShape; | |
791 | fInfo_.SIM_uses_FLARB = useFLARB; | |
792 | fInfo_.SIM_uses_RF = useRF; | |
793 | fInfo_.SIM_uses_SF = useSF; | |
794 | < | |
795 | < | if( myMethod == "REACTION_FIELD") { |
796 | < | |
638 | < | if (simParams_->haveDielectric()) { |
639 | < | fInfo_.dielect = simParams_->getDielectric(); |
640 | < | } else { |
641 | < | sprintf(painCave.errMsg, |
642 | < | "SimSetup Error: No Dielectric constant was set.\n" |
643 | < | "\tYou are trying to use Reaction Field without" |
644 | < | "\tsetting a dielectric constant!\n"); |
645 | < | painCave.isFatal = 1; |
646 | < | simError(); |
647 | < | } |
648 | < | } |
649 | < | |
794 | > | fInfo_.SIM_uses_SP = useSP; |
795 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
796 | > | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
797 | } | |
798 | ||
799 | void SimInfo::setupFortranSim() { | |
# | Line 663 | Line 810 | namespace oopse { | |
810 | } | |
811 | ||
812 | //calculate mass ratio of cutoff group | |
813 | < | std::vector<double> mfact; |
813 | > | std::vector<RealType> mfact; |
814 | SimInfo::MoleculeIterator mi; | |
815 | Molecule* mol; | |
816 | Molecule::CutoffGroupIterator ci; | |
817 | CutoffGroup* cg; | |
818 | Molecule::AtomIterator ai; | |
819 | Atom* atom; | |
820 | < | double totalMass; |
820 | > | RealType totalMass; |
821 | ||
822 | //to avoid memory reallocation, reserve enough space for mfact | |
823 | mfact.reserve(getNCutoffGroups()); | |
# | Line 731 | Line 878 | namespace oopse { | |
878 | "succesfully sent the simulation information to fortran.\n"); | |
879 | MPIcheckPoint(); | |
880 | #endif // is_mpi | |
881 | + | |
882 | + | // Setup number of neighbors in neighbor list if present |
883 | + | if (simParams_->haveNeighborListNeighbors()) { |
884 | + | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
885 | + | setNeighbors(&nlistNeighbors); |
886 | + | } |
887 | + | |
888 | + | |
889 | } | |
890 | ||
891 | ||
# | Line 793 | Line 948 | namespace oopse { | |
948 | ||
949 | #endif | |
950 | ||
951 | < | double SimInfo::calcMaxCutoffRadius() { |
797 | < | |
798 | < | |
799 | < | std::set<AtomType*> atomTypes; |
800 | < | std::set<AtomType*>::iterator i; |
801 | < | std::vector<double> cutoffRadius; |
802 | < | |
803 | < | //get the unique atom types |
804 | < | atomTypes = getUniqueAtomTypes(); |
805 | < | |
806 | < | //query the max cutoff radius among these atom types |
807 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
808 | < | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); |
809 | < | } |
810 | < | |
811 | < | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); |
812 | < | #ifdef IS_MPI |
813 | < | //pick the max cutoff radius among the processors |
814 | < | #endif |
815 | < | |
816 | < | return maxCutoffRadius; |
817 | < | } |
818 | < | |
819 | < | void SimInfo::getCutoff(double& rcut, double& rsw) { |
951 | > | void SimInfo::setupCutoff() { |
952 | ||
953 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
822 | < | |
823 | < | if (!simParams_->haveCutoffRadius()){ |
824 | < | sprintf(painCave.errMsg, |
825 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
826 | < | "\tOOPSE will use a default value of 15.0 angstroms" |
827 | < | "\tfor the cutoffRadius.\n"); |
828 | < | painCave.isFatal = 0; |
829 | < | simError(); |
830 | < | rcut = 15.0; |
831 | < | } else{ |
832 | < | rcut = simParams_->getCutoffRadius(); |
833 | < | } |
953 | > | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
954 | ||
955 | < | if (!simParams_->haveSwitchingRadius()){ |
956 | < | sprintf(painCave.errMsg, |
837 | < | "SimCreator Warning: No value was set for switchingRadius.\n" |
838 | < | "\tOOPSE will use a default value of\n" |
839 | < | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
840 | < | painCave.isFatal = 0; |
841 | < | simError(); |
842 | < | rsw = 0.85 * rcut; |
843 | < | } else{ |
844 | < | rsw = simParams_->getSwitchingRadius(); |
845 | < | } |
955 | > | // Check the cutoff policy |
956 | > | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default |
957 | ||
958 | < | } else { |
959 | < | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in |
960 | < | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used |
850 | < | |
851 | < | if (simParams_->haveCutoffRadius()) { |
852 | < | rcut = simParams_->getCutoffRadius(); |
853 | < | } else { |
854 | < | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system |
855 | < | rcut = calcMaxCutoffRadius(); |
856 | < | } |
958 | > | // Set LJ shifting bools to false |
959 | > | ljsp_ = false; |
960 | > | ljsf_ = false; |
961 | ||
962 | < | if (simParams_->haveSwitchingRadius()) { |
963 | < | rsw = simParams_->getSwitchingRadius(); |
964 | < | } else { |
965 | < | rsw = rcut; |
966 | < | } |
863 | < | |
962 | > | std::string myPolicy; |
963 | > | if (forceFieldOptions_.haveCutoffPolicy()){ |
964 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); |
965 | > | }else if (simParams_->haveCutoffPolicy()) { |
966 | > | myPolicy = simParams_->getCutoffPolicy(); |
967 | } | |
865 | – | } |
968 | ||
969 | < | void SimInfo::setupCutoff() { |
868 | < | getCutoff(rcut_, rsw_); |
869 | < | double rnblist = rcut_ + 1; // skin of neighbor list |
870 | < | |
871 | < | //Pass these cutoff radius etc. to fortran. This function should be called once and only once |
872 | < | |
873 | < | int cp = TRADITIONAL_CUTOFF_POLICY; |
874 | < | if (simParams_->haveCutoffPolicy()) { |
875 | < | std::string myPolicy = simParams_->getCutoffPolicy(); |
969 | > | if (!myPolicy.empty()){ |
970 | toUpper(myPolicy); | |
971 | if (myPolicy == "MIX") { | |
972 | cp = MIX_CUTOFF_POLICY; | |
# | Line 891 | Line 985 | namespace oopse { | |
985 | } | |
986 | } | |
987 | } | |
988 | < | } |
988 | > | } |
989 | > | notifyFortranCutoffPolicy(&cp); |
990 | ||
991 | < | |
991 | > | // Check the Skin Thickness for neighborlists |
992 | > | RealType skin; |
993 | if (simParams_->haveSkinThickness()) { | |
994 | < | double skinThickness = simParams_->getSkinThickness(); |
995 | < | } |
994 | > | skin = simParams_->getSkinThickness(); |
995 | > | notifyFortranSkinThickness(&skin); |
996 | > | } |
997 | > | |
998 | > | // Check if the cutoff was set explicitly: |
999 | > | if (simParams_->haveCutoffRadius()) { |
1000 | > | rcut_ = simParams_->getCutoffRadius(); |
1001 | > | if (simParams_->haveSwitchingRadius()) { |
1002 | > | rsw_ = simParams_->getSwitchingRadius(); |
1003 | > | } else { |
1004 | > | if (fInfo_.SIM_uses_Charges | |
1005 | > | fInfo_.SIM_uses_Dipoles | |
1006 | > | fInfo_.SIM_uses_RF) { |
1007 | > | |
1008 | > | rsw_ = 0.85 * rcut_; |
1009 | > | sprintf(painCave.errMsg, |
1010 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1011 | > | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" |
1012 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1013 | > | painCave.isFatal = 0; |
1014 | > | simError(); |
1015 | > | } else { |
1016 | > | rsw_ = rcut_; |
1017 | > | sprintf(painCave.errMsg, |
1018 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1019 | > | "\tOOPSE will use the same value as the cutoffRadius.\n" |
1020 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1021 | > | painCave.isFatal = 0; |
1022 | > | simError(); |
1023 | > | } |
1024 | > | } |
1025 | ||
1026 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
1027 | < | // also send cutoff notification to electrostatics |
1028 | < | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
1026 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
1027 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1028 | > | toUpper(myMethod); |
1029 | > | |
1030 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
1031 | > | ljsp_ = true; |
1032 | > | } else if (myMethod == "SHIFTED_FORCE") { |
1033 | > | ljsf_ = true; |
1034 | > | } |
1035 | > | } |
1036 | > | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1037 | > | |
1038 | > | } else { |
1039 | > | |
1040 | > | // For electrostatic atoms, we'll assume a large safe value: |
1041 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
1042 | > | sprintf(painCave.errMsg, |
1043 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
1044 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
1045 | > | "\tfor the cutoffRadius.\n"); |
1046 | > | painCave.isFatal = 0; |
1047 | > | simError(); |
1048 | > | rcut_ = 15.0; |
1049 | > | |
1050 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
1051 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1052 | > | toUpper(myMethod); |
1053 | > | |
1054 | > | // For the time being, we're tethering the LJ shifted behavior to the |
1055 | > | // electrostaticSummationMethod keyword options |
1056 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
1057 | > | ljsp_ = true; |
1058 | > | } else if (myMethod == "SHIFTED_FORCE") { |
1059 | > | ljsf_ = true; |
1060 | > | } |
1061 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1062 | > | if (simParams_->haveSwitchingRadius()){ |
1063 | > | sprintf(painCave.errMsg, |
1064 | > | "SimInfo Warning: A value was set for the switchingRadius\n" |
1065 | > | "\teven though the electrostaticSummationMethod was\n" |
1066 | > | "\tset to %s\n", myMethod.c_str()); |
1067 | > | painCave.isFatal = 1; |
1068 | > | simError(); |
1069 | > | } |
1070 | > | } |
1071 | > | } |
1072 | > | |
1073 | > | if (simParams_->haveSwitchingRadius()){ |
1074 | > | rsw_ = simParams_->getSwitchingRadius(); |
1075 | > | } else { |
1076 | > | sprintf(painCave.errMsg, |
1077 | > | "SimCreator Warning: No value was set for switchingRadius.\n" |
1078 | > | "\tOOPSE will use a default value of\n" |
1079 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
1080 | > | painCave.isFatal = 0; |
1081 | > | simError(); |
1082 | > | rsw_ = 0.85 * rcut_; |
1083 | > | } |
1084 | > | |
1085 | > | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1086 | > | |
1087 | > | } else { |
1088 | > | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so |
1089 | > | // We'll punt and let fortran figure out the cutoffs later. |
1090 | > | |
1091 | > | notifyFortranYouAreOnYourOwn(); |
1092 | > | |
1093 | > | } |
1094 | > | } |
1095 | } | |
1096 | ||
1097 | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | |
# | Line 908 | Line 1099 | namespace oopse { | |
1099 | int errorOut; | |
1100 | int esm = NONE; | |
1101 | int sm = UNDAMPED; | |
1102 | < | double alphaVal; |
1103 | < | double dielectric; |
1104 | < | |
1102 | > | RealType alphaVal; |
1103 | > | RealType dielectric; |
1104 | > | |
1105 | errorOut = isError; | |
915 | – | alphaVal = simParams_->getDampingAlpha(); |
916 | – | dielectric = simParams_->getDielectric(); |
1106 | ||
1107 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1108 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
# | Line 930 | Line 1119 | namespace oopse { | |
1119 | if (myMethod == "SHIFTED_FORCE") { | |
1120 | esm = SHIFTED_FORCE; | |
1121 | } else { | |
1122 | < | if (myMethod == "REACTION_FIELD") { |
1122 | > | if (myMethod == "REACTION_FIELD") { |
1123 | esm = REACTION_FIELD; | |
1124 | + | dielectric = simParams_->getDielectric(); |
1125 | + | if (!simParams_->haveDielectric()) { |
1126 | + | // throw warning |
1127 | + | sprintf( painCave.errMsg, |
1128 | + | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1129 | + | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1130 | + | painCave.isFatal = 0; |
1131 | + | simError(); |
1132 | + | } |
1133 | } else { | |
1134 | // throw error | |
1135 | sprintf( painCave.errMsg, | |
1136 | < | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() ); |
1136 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" |
1137 | > | "\t(Input file specified %s .)\n" |
1138 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" |
1139 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" |
1140 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); |
1141 | painCave.isFatal = 1; | |
1142 | simError(); | |
1143 | } | |
# | Line 954 | Line 1156 | namespace oopse { | |
1156 | if (myScreen == "DAMPED") { | |
1157 | sm = DAMPED; | |
1158 | if (!simParams_->haveDampingAlpha()) { | |
1159 | < | //throw error |
1159 | > | // first set a cutoff dependent alpha value |
1160 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1161 | > | alphaVal = 0.5125 - rcut_* 0.025; |
1162 | > | // for values rcut > 20.5, alpha is zero |
1163 | > | if (alphaVal < 0) alphaVal = 0; |
1164 | > | |
1165 | > | // throw warning |
1166 | sprintf( painCave.errMsg, | |
1167 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); |
1167 | > | "SimInfo warning: dampingAlpha was not specified in the input file.\n" |
1168 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1169 | painCave.isFatal = 0; | |
1170 | simError(); | |
1171 | + | } else { |
1172 | + | alphaVal = simParams_->getDampingAlpha(); |
1173 | } | |
1174 | + | |
1175 | } else { | |
1176 | // throw error | |
1177 | sprintf( painCave.errMsg, | |
1178 | < | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
1178 | > | "SimInfo error: Unknown electrostaticScreeningMethod.\n" |
1179 | > | "\t(Input file specified %s .)\n" |
1180 | > | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" |
1181 | > | "or \"damped\".\n", myScreen.c_str() ); |
1182 | painCave.isFatal = 1; | |
1183 | simError(); | |
1184 | } | |
# | Line 972 | Line 1187 | namespace oopse { | |
1187 | ||
1188 | // let's pass some summation method variables to fortran | |
1189 | setElectrostaticSummationMethod( &esm ); | |
1190 | + | setFortranElectrostaticMethod( &esm ); |
1191 | setScreeningMethod( &sm ); | |
1192 | setDampingAlpha( &alphaVal ); | |
1193 | setReactionFieldDielectric( &dielectric ); | |
1194 | < | initFortranFF( &esm, &errorOut ); |
1194 | > | initFortranFF( &errorOut ); |
1195 | } | |
1196 | ||
1197 | + | void SimInfo::setupSwitchingFunction() { |
1198 | + | int ft = CUBIC; |
1199 | + | |
1200 | + | if (simParams_->haveSwitchingFunctionType()) { |
1201 | + | std::string funcType = simParams_->getSwitchingFunctionType(); |
1202 | + | toUpper(funcType); |
1203 | + | if (funcType == "CUBIC") { |
1204 | + | ft = CUBIC; |
1205 | + | } else { |
1206 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1207 | + | ft = FIFTH_ORDER_POLY; |
1208 | + | } else { |
1209 | + | // throw error |
1210 | + | sprintf( painCave.errMsg, |
1211 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1212 | + | painCave.isFatal = 1; |
1213 | + | simError(); |
1214 | + | } |
1215 | + | } |
1216 | + | } |
1217 | + | |
1218 | + | // send switching function notification to switcheroo |
1219 | + | setFunctionType(&ft); |
1220 | + | |
1221 | + | } |
1222 | + | |
1223 | + | void SimInfo::setupAccumulateBoxDipole() { |
1224 | + | |
1225 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1226 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1227 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1228 | + | setAccumulateBoxDipole(); |
1229 | + | calcBoxDipole_ = true; |
1230 | + | } |
1231 | + | |
1232 | + | } |
1233 | + | |
1234 | void SimInfo::addProperty(GenericData* genData) { | |
1235 | properties_.addProperty(genData); | |
1236 | } | |
# | Line 1034 | Line 1287 | namespace oopse { | |
1287 | Molecule* mol; | |
1288 | ||
1289 | Vector3d comVel(0.0); | |
1290 | < | double totalMass = 0.0; |
1290 | > | RealType totalMass = 0.0; |
1291 | ||
1292 | ||
1293 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1294 | < | double mass = mol->getMass(); |
1294 | > | RealType mass = mol->getMass(); |
1295 | totalMass += mass; | |
1296 | comVel += mass * mol->getComVel(); | |
1297 | } | |
1298 | ||
1299 | #ifdef IS_MPI | |
1300 | < | double tmpMass = totalMass; |
1300 | > | RealType tmpMass = totalMass; |
1301 | Vector3d tmpComVel(comVel); | |
1302 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1303 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1302 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1303 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1304 | #endif | |
1305 | ||
1306 | comVel /= totalMass; | |
# | Line 1060 | Line 1313 | namespace oopse { | |
1313 | Molecule* mol; | |
1314 | ||
1315 | Vector3d com(0.0); | |
1316 | < | double totalMass = 0.0; |
1316 | > | RealType totalMass = 0.0; |
1317 | ||
1318 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1319 | < | double mass = mol->getMass(); |
1319 | > | RealType mass = mol->getMass(); |
1320 | totalMass += mass; | |
1321 | com += mass * mol->getCom(); | |
1322 | } | |
1323 | ||
1324 | #ifdef IS_MPI | |
1325 | < | double tmpMass = totalMass; |
1325 | > | RealType tmpMass = totalMass; |
1326 | Vector3d tmpCom(com); | |
1327 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1328 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1327 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1328 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1329 | #endif | |
1330 | ||
1331 | com /= totalMass; | |
# | Line 1096 | Line 1349 | namespace oopse { | |
1349 | Molecule* mol; | |
1350 | ||
1351 | ||
1352 | < | double totalMass = 0.0; |
1352 | > | RealType totalMass = 0.0; |
1353 | ||
1354 | ||
1355 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1356 | < | double mass = mol->getMass(); |
1356 | > | RealType mass = mol->getMass(); |
1357 | totalMass += mass; | |
1358 | com += mass * mol->getCom(); | |
1359 | comVel += mass * mol->getComVel(); | |
1360 | } | |
1361 | ||
1362 | #ifdef IS_MPI | |
1363 | < | double tmpMass = totalMass; |
1363 | > | RealType tmpMass = totalMass; |
1364 | Vector3d tmpCom(com); | |
1365 | Vector3d tmpComVel(comVel); | |
1366 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1367 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1368 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1366 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1367 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1368 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1369 | #endif | |
1370 | ||
1371 | com /= totalMass; | |
# | Line 1131 | Line 1384 | namespace oopse { | |
1384 | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | |
1385 | ||
1386 | ||
1387 | < | double xx = 0.0; |
1388 | < | double yy = 0.0; |
1389 | < | double zz = 0.0; |
1390 | < | double xy = 0.0; |
1391 | < | double xz = 0.0; |
1392 | < | double yz = 0.0; |
1387 | > | RealType xx = 0.0; |
1388 | > | RealType yy = 0.0; |
1389 | > | RealType zz = 0.0; |
1390 | > | RealType xy = 0.0; |
1391 | > | RealType xz = 0.0; |
1392 | > | RealType yz = 0.0; |
1393 | Vector3d com(0.0); | |
1394 | Vector3d comVel(0.0); | |
1395 | ||
# | Line 1148 | Line 1401 | namespace oopse { | |
1401 | Vector3d thisq(0.0); | |
1402 | Vector3d thisv(0.0); | |
1403 | ||
1404 | < | double thisMass = 0.0; |
1404 | > | RealType thisMass = 0.0; |
1405 | ||
1406 | ||
1407 | ||
# | Line 1186 | Line 1439 | namespace oopse { | |
1439 | #ifdef IS_MPI | |
1440 | Mat3x3d tmpI(inertiaTensor); | |
1441 | Vector3d tmpAngMom; | |
1442 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1443 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1442 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1443 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1444 | #endif | |
1445 | ||
1446 | return; | |
# | Line 1208 | Line 1461 | namespace oopse { | |
1461 | Vector3d thisr(0.0); | |
1462 | Vector3d thisp(0.0); | |
1463 | ||
1464 | < | double thisMass; |
1464 | > | RealType thisMass; |
1465 | ||
1466 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1467 | thisMass = mol->getMass(); | |
# | Line 1221 | Line 1474 | namespace oopse { | |
1474 | ||
1475 | #ifdef IS_MPI | |
1476 | Vector3d tmpAngMom; | |
1477 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1477 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1478 | #endif | |
1479 | ||
1480 | return angularMomentum; | |
1481 | } | |
1482 | ||
1483 | < | |
1483 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { |
1484 | > | return IOIndexToIntegrableObject.at(index); |
1485 | > | } |
1486 | > | |
1487 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1488 | > | IOIndexToIntegrableObject= v; |
1489 | > | } |
1490 | > | |
1491 | > | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1492 | > | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1493 | > | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1494 | > | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1495 | > | */ |
1496 | > | void SimInfo::getGyrationalVolume(RealType &volume){ |
1497 | > | Mat3x3d intTensor; |
1498 | > | RealType det; |
1499 | > | Vector3d dummyAngMom; |
1500 | > | RealType sysconstants; |
1501 | > | RealType geomCnst; |
1502 | > | |
1503 | > | geomCnst = 3.0/2.0; |
1504 | > | /* Get the inertial tensor and angular momentum for free*/ |
1505 | > | getInertiaTensor(intTensor,dummyAngMom); |
1506 | > | |
1507 | > | det = intTensor.determinant(); |
1508 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1509 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1510 | > | return; |
1511 | > | } |
1512 | > | |
1513 | > | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1514 | > | Mat3x3d intTensor; |
1515 | > | Vector3d dummyAngMom; |
1516 | > | RealType sysconstants; |
1517 | > | RealType geomCnst; |
1518 | > | |
1519 | > | geomCnst = 3.0/2.0; |
1520 | > | /* Get the inertial tensor and angular momentum for free*/ |
1521 | > | getInertiaTensor(intTensor,dummyAngMom); |
1522 | > | |
1523 | > | detI = intTensor.determinant(); |
1524 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1525 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1526 | > | return; |
1527 | > | } |
1528 | > | /* |
1529 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1530 | > | assert( v.size() == nAtoms_ + nRigidBodies_); |
1531 | > | sdByGlobalIndex_ = v; |
1532 | > | } |
1533 | > | |
1534 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1535 | > | //assert(index < nAtoms_ + nRigidBodies_); |
1536 | > | return sdByGlobalIndex_.at(index); |
1537 | > | } |
1538 | > | */ |
1539 | }//end namespace oopse | |
1540 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |