# | Line 54 | Line 54 | |
---|---|---|
54 | #include "math/Vector3.hpp" | |
55 | #include "primitives/Molecule.hpp" | |
56 | #include "primitives/StuntDouble.hpp" | |
57 | – | #include "UseTheForce/fCutoffPolicy.h" |
58 | – | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" |
57 | #include "UseTheForce/doForces_interface.h" | |
58 | #include "UseTheForce/DarkSide/neighborLists_interface.h" | |
61 | – | #include "UseTheForce/DarkSide/switcheroo_interface.h" |
59 | #include "utils/MemoryUtils.hpp" | |
60 | #include "utils/simError.h" | |
61 | #include "selection/SelectionManager.hpp" | |
62 | #include "io/ForceFieldOptions.hpp" | |
63 | #include "UseTheForce/ForceField.hpp" | |
64 | + | #include "nonbonded/SwitchingFunction.hpp" |
65 | ||
66 | ||
67 | #ifdef IS_MPI | |
# | Line 71 | Line 69 | |
69 | #include "UseTheForce/DarkSide/simParallel_interface.h" | |
70 | #endif | |
71 | ||
72 | + | using namespace std; |
73 | namespace OpenMD { | |
75 | – | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { |
76 | – | std::map<int, std::set<int> >::iterator i = container.find(index); |
77 | – | std::set<int> result; |
78 | – | if (i != container.end()) { |
79 | – | result = i->second; |
80 | – | } |
81 | – | |
82 | – | return result; |
83 | – | } |
74 | ||
75 | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | |
76 | forceField_(ff), simParams_(simParams), | |
# | Line 90 | Line 80 | namespace OpenMD { | |
80 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), | |
81 | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | |
82 | nConstraints_(0), sman_(NULL), fortranInitialized_(false), | |
83 | < | calcBoxDipole_(false), useAtomicVirial_(true) { |
84 | < | |
85 | < | |
86 | < | MoleculeStamp* molStamp; |
87 | < | int nMolWithSameStamp; |
88 | < | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
89 | < | int nGroups = 0; //total cutoff groups defined in meta-data file |
90 | < | CutoffGroupStamp* cgStamp; |
91 | < | RigidBodyStamp* rbStamp; |
92 | < | int nRigidAtoms = 0; |
93 | < | |
94 | < | std::vector<Component*> components = simParams->getComponents(); |
83 | > | calcBoxDipole_(false), useAtomicVirial_(true) { |
84 | > | |
85 | > | MoleculeStamp* molStamp; |
86 | > | int nMolWithSameStamp; |
87 | > | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
88 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
89 | > | CutoffGroupStamp* cgStamp; |
90 | > | RigidBodyStamp* rbStamp; |
91 | > | int nRigidAtoms = 0; |
92 | > | |
93 | > | vector<Component*> components = simParams->getComponents(); |
94 | > | |
95 | > | for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
96 | > | molStamp = (*i)->getMoleculeStamp(); |
97 | > | nMolWithSameStamp = (*i)->getNMol(); |
98 | ||
99 | < | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
100 | < | molStamp = (*i)->getMoleculeStamp(); |
101 | < | nMolWithSameStamp = (*i)->getNMol(); |
102 | < | |
103 | < | addMoleculeStamp(molStamp, nMolWithSameStamp); |
104 | < | |
105 | < | //calculate atoms in molecules |
106 | < | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
107 | < | |
108 | < | //calculate atoms in cutoff groups |
109 | < | int nAtomsInGroups = 0; |
110 | < | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); |
118 | < | |
119 | < | for (int j=0; j < nCutoffGroupsInStamp; j++) { |
120 | < | cgStamp = molStamp->getCutoffGroupStamp(j); |
121 | < | nAtomsInGroups += cgStamp->getNMembers(); |
122 | < | } |
123 | < | |
124 | < | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; |
125 | < | |
126 | < | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; |
127 | < | |
128 | < | //calculate atoms in rigid bodies |
129 | < | int nAtomsInRigidBodies = 0; |
130 | < | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); |
131 | < | |
132 | < | for (int j=0; j < nRigidBodiesInStamp; j++) { |
133 | < | rbStamp = molStamp->getRigidBodyStamp(j); |
134 | < | nAtomsInRigidBodies += rbStamp->getNMembers(); |
135 | < | } |
136 | < | |
137 | < | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; |
138 | < | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; |
139 | < | |
99 | > | addMoleculeStamp(molStamp, nMolWithSameStamp); |
100 | > | |
101 | > | //calculate atoms in molecules |
102 | > | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
103 | > | |
104 | > | //calculate atoms in cutoff groups |
105 | > | int nAtomsInGroups = 0; |
106 | > | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); |
107 | > | |
108 | > | for (int j=0; j < nCutoffGroupsInStamp; j++) { |
109 | > | cgStamp = molStamp->getCutoffGroupStamp(j); |
110 | > | nAtomsInGroups += cgStamp->getNMembers(); |
111 | } | |
112 | < | |
113 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff |
114 | < | //group therefore the total number of cutoff groups in the system is |
115 | < | //equal to the total number of atoms minus number of atoms belong to |
116 | < | //cutoff group defined in meta-data file plus the number of cutoff |
117 | < | //groups defined in meta-data file |
118 | < | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
119 | < | |
120 | < | //every free atom (atom does not belong to rigid bodies) is an |
121 | < | //integrable object therefore the total number of integrable objects |
122 | < | //in the system is equal to the total number of atoms minus number of |
123 | < | //atoms belong to rigid body defined in meta-data file plus the number |
124 | < | //of rigid bodies defined in meta-data file |
125 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
126 | < | + nGlobalRigidBodies_; |
127 | < | |
128 | < | nGlobalMols_ = molStampIds_.size(); |
129 | < | molToProcMap_.resize(nGlobalMols_); |
130 | < | } |
131 | < | |
112 | > | |
113 | > | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; |
114 | > | |
115 | > | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; |
116 | > | |
117 | > | //calculate atoms in rigid bodies |
118 | > | int nAtomsInRigidBodies = 0; |
119 | > | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); |
120 | > | |
121 | > | for (int j=0; j < nRigidBodiesInStamp; j++) { |
122 | > | rbStamp = molStamp->getRigidBodyStamp(j); |
123 | > | nAtomsInRigidBodies += rbStamp->getNMembers(); |
124 | > | } |
125 | > | |
126 | > | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; |
127 | > | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; |
128 | > | |
129 | > | } |
130 | > | |
131 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff |
132 | > | //group therefore the total number of cutoff groups in the system is |
133 | > | //equal to the total number of atoms minus number of atoms belong to |
134 | > | //cutoff group defined in meta-data file plus the number of cutoff |
135 | > | //groups defined in meta-data file |
136 | > | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
137 | > | |
138 | > | //every free atom (atom does not belong to rigid bodies) is an |
139 | > | //integrable object therefore the total number of integrable objects |
140 | > | //in the system is equal to the total number of atoms minus number of |
141 | > | //atoms belong to rigid body defined in meta-data file plus the number |
142 | > | //of rigid bodies defined in meta-data file |
143 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
144 | > | + nGlobalRigidBodies_; |
145 | > | |
146 | > | nGlobalMols_ = molStampIds_.size(); |
147 | > | molToProcMap_.resize(nGlobalMols_); |
148 | > | } |
149 | > | |
150 | SimInfo::~SimInfo() { | |
151 | < | std::map<int, Molecule*>::iterator i; |
151 | > | map<int, Molecule*>::iterator i; |
152 | for (i = molecules_.begin(); i != molecules_.end(); ++i) { | |
153 | delete i->second; | |
154 | } | |
# | Line 170 | Line 159 | namespace OpenMD { | |
159 | delete forceField_; | |
160 | } | |
161 | ||
173 | – | int SimInfo::getNGlobalConstraints() { |
174 | – | int nGlobalConstraints; |
175 | – | #ifdef IS_MPI |
176 | – | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
177 | – | MPI_COMM_WORLD); |
178 | – | #else |
179 | – | nGlobalConstraints = nConstraints_; |
180 | – | #endif |
181 | – | return nGlobalConstraints; |
182 | – | } |
162 | ||
163 | bool SimInfo::addMolecule(Molecule* mol) { | |
164 | MoleculeIterator i; | |
165 | < | |
165 | > | |
166 | i = molecules_.find(mol->getGlobalIndex()); | |
167 | if (i == molecules_.end() ) { | |
168 | < | |
169 | < | molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); |
170 | < | |
168 | > | |
169 | > | molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); |
170 | > | |
171 | nAtoms_ += mol->getNAtoms(); | |
172 | nBonds_ += mol->getNBonds(); | |
173 | nBends_ += mol->getNBends(); | |
# | Line 198 | Line 177 | namespace OpenMD { | |
177 | nIntegrableObjects_ += mol->getNIntegrableObjects(); | |
178 | nCutoffGroups_ += mol->getNCutoffGroups(); | |
179 | nConstraints_ += mol->getNConstraintPairs(); | |
180 | < | |
180 | > | |
181 | addInteractionPairs(mol); | |
182 | < | |
182 | > | |
183 | return true; | |
184 | } else { | |
185 | return false; | |
186 | } | |
187 | } | |
188 | < | |
188 | > | |
189 | bool SimInfo::removeMolecule(Molecule* mol) { | |
190 | MoleculeIterator i; | |
191 | i = molecules_.find(mol->getGlobalIndex()); | |
# | Line 234 | Line 213 | namespace OpenMD { | |
213 | } else { | |
214 | return false; | |
215 | } | |
237 | – | |
238 | – | |
216 | } | |
217 | ||
218 | ||
# | Line 253 | Line 230 | namespace OpenMD { | |
230 | void SimInfo::calcNdf() { | |
231 | int ndf_local; | |
232 | MoleculeIterator i; | |
233 | < | std::vector<StuntDouble*>::iterator j; |
233 | > | vector<StuntDouble*>::iterator j; |
234 | Molecule* mol; | |
235 | StuntDouble* integrableObject; | |
236 | ||
# | Line 304 | Line 281 | namespace OpenMD { | |
281 | int ndfRaw_local; | |
282 | ||
283 | MoleculeIterator i; | |
284 | < | std::vector<StuntDouble*>::iterator j; |
284 | > | vector<StuntDouble*>::iterator j; |
285 | Molecule* mol; | |
286 | StuntDouble* integrableObject; | |
287 | ||
# | Line 353 | Line 330 | namespace OpenMD { | |
330 | ||
331 | void SimInfo::addInteractionPairs(Molecule* mol) { | |
332 | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | |
333 | < | std::vector<Bond*>::iterator bondIter; |
334 | < | std::vector<Bend*>::iterator bendIter; |
335 | < | std::vector<Torsion*>::iterator torsionIter; |
336 | < | std::vector<Inversion*>::iterator inversionIter; |
333 | > | vector<Bond*>::iterator bondIter; |
334 | > | vector<Bend*>::iterator bendIter; |
335 | > | vector<Torsion*>::iterator torsionIter; |
336 | > | vector<Inversion*>::iterator inversionIter; |
337 | Bond* bond; | |
338 | Bend* bend; | |
339 | Torsion* torsion; | |
# | Line 374 | Line 351 | namespace OpenMD { | |
351 | // always be excluded. These are done at the bottom of this | |
352 | // function. | |
353 | ||
354 | < | std::map<int, std::set<int> > atomGroups; |
354 | > | map<int, set<int> > atomGroups; |
355 | Molecule::RigidBodyIterator rbIter; | |
356 | RigidBody* rb; | |
357 | Molecule::IntegrableObjectIterator ii; | |
# | Line 386 | Line 363 | namespace OpenMD { | |
363 | ||
364 | if (integrableObject->isRigidBody()) { | |
365 | rb = static_cast<RigidBody*>(integrableObject); | |
366 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
367 | < | std::set<int> rigidAtoms; |
366 | > | vector<Atom*> atoms = rb->getAtoms(); |
367 | > | set<int> rigidAtoms; |
368 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
369 | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | |
370 | } | |
371 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
372 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
372 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
373 | } | |
374 | } else { | |
375 | < | std::set<int> oneAtomSet; |
375 | > | set<int> oneAtomSet; |
376 | oneAtomSet.insert(integrableObject->getGlobalIndex()); | |
377 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
377 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
378 | } | |
379 | } | |
380 | ||
# | Line 500 | Line 477 | namespace OpenMD { | |
477 | ||
478 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | |
479 | rb = mol->nextRigidBody(rbIter)) { | |
480 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
480 | > | vector<Atom*> atoms = rb->getAtoms(); |
481 | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | |
482 | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | |
483 | a = atoms[i]->getGlobalIndex(); | |
# | Line 514 | Line 491 | namespace OpenMD { | |
491 | ||
492 | void SimInfo::removeInteractionPairs(Molecule* mol) { | |
493 | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | |
494 | < | std::vector<Bond*>::iterator bondIter; |
495 | < | std::vector<Bend*>::iterator bendIter; |
496 | < | std::vector<Torsion*>::iterator torsionIter; |
497 | < | std::vector<Inversion*>::iterator inversionIter; |
494 | > | vector<Bond*>::iterator bondIter; |
495 | > | vector<Bend*>::iterator bendIter; |
496 | > | vector<Torsion*>::iterator torsionIter; |
497 | > | vector<Inversion*>::iterator inversionIter; |
498 | Bond* bond; | |
499 | Bend* bend; | |
500 | Torsion* torsion; | |
# | Line 527 | Line 504 | namespace OpenMD { | |
504 | int c; | |
505 | int d; | |
506 | ||
507 | < | std::map<int, std::set<int> > atomGroups; |
507 | > | map<int, set<int> > atomGroups; |
508 | Molecule::RigidBodyIterator rbIter; | |
509 | RigidBody* rb; | |
510 | Molecule::IntegrableObjectIterator ii; | |
# | Line 539 | Line 516 | namespace OpenMD { | |
516 | ||
517 | if (integrableObject->isRigidBody()) { | |
518 | rb = static_cast<RigidBody*>(integrableObject); | |
519 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
520 | < | std::set<int> rigidAtoms; |
519 | > | vector<Atom*> atoms = rb->getAtoms(); |
520 | > | set<int> rigidAtoms; |
521 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
522 | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | |
523 | } | |
524 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
525 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
525 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
526 | } | |
527 | } else { | |
528 | < | std::set<int> oneAtomSet; |
528 | > | set<int> oneAtomSet; |
529 | oneAtomSet.insert(integrableObject->getGlobalIndex()); | |
530 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
530 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
531 | } | |
532 | } | |
533 | ||
# | Line 653 | Line 630 | namespace OpenMD { | |
630 | ||
631 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | |
632 | rb = mol->nextRigidBody(rbIter)) { | |
633 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
633 | > | vector<Atom*> atoms = rb->getAtoms(); |
634 | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | |
635 | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | |
636 | a = atoms[i]->getGlobalIndex(); | |
# | Line 676 | Line 653 | namespace OpenMD { | |
653 | molStampIds_.insert(molStampIds_.end(), nmol, curStampId); | |
654 | } | |
655 | ||
679 | – | void SimInfo::update() { |
656 | ||
657 | < | setupSimType(); |
657 | > | /** |
658 | > | * update |
659 | > | * |
660 | > | * Performs the global checks and variable settings after the objects have been |
661 | > | * created. |
662 | > | * |
663 | > | */ |
664 | > | void SimInfo::update() { |
665 | > | |
666 | > | setupSimVariables(); |
667 | > | setupCutoffs(); |
668 | > | setupSwitching(); |
669 | > | setupElectrostatics(); |
670 | > | setupNeighborlists(); |
671 | ||
672 | #ifdef IS_MPI | |
673 | setupFortranParallel(); | |
674 | #endif | |
686 | – | |
675 | setupFortranSim(); | |
676 | + | fortranInitialized_ = true; |
677 | ||
689 | – | //setup fortran force field |
690 | – | /** @deprecate */ |
691 | – | int isError = 0; |
692 | – | |
693 | – | setupCutoff(); |
694 | – | |
695 | – | setupElectrostaticSummationMethod( isError ); |
696 | – | setupSwitchingFunction(); |
697 | – | setupAccumulateBoxDipole(); |
698 | – | |
699 | – | if(isError){ |
700 | – | sprintf( painCave.errMsg, |
701 | – | "ForceField error: There was an error initializing the forceField in fortran.\n" ); |
702 | – | painCave.isFatal = 1; |
703 | – | simError(); |
704 | – | } |
705 | – | |
678 | calcNdf(); | |
679 | calcNdfRaw(); | |
680 | calcNdfTrans(); | |
709 | – | |
710 | – | fortranInitialized_ = true; |
681 | } | |
682 | < | |
683 | < | std::set<AtomType*> SimInfo::getUniqueAtomTypes() { |
682 | > | |
683 | > | set<AtomType*> SimInfo::getSimulatedAtomTypes() { |
684 | SimInfo::MoleculeIterator mi; | |
685 | Molecule* mol; | |
686 | Molecule::AtomIterator ai; | |
687 | Atom* atom; | |
688 | < | std::set<AtomType*> atomTypes; |
689 | < | |
690 | < | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
721 | < | |
688 | > | set<AtomType*> atomTypes; |
689 | > | |
690 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
691 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
692 | atomTypes.insert(atom->getAtomType()); | |
693 | < | } |
694 | < | |
726 | < | } |
727 | < | |
693 | > | } |
694 | > | } |
695 | return atomTypes; | |
696 | } | |
697 | ||
698 | < | void SimInfo::setupSimType() { |
699 | < | std::set<AtomType*>::iterator i; |
700 | < | std::set<AtomType*> atomTypes; |
701 | < | atomTypes = getUniqueAtomTypes(); |
698 | > | /** |
699 | > | * setupCutoffs |
700 | > | * |
701 | > | * Sets the values of cutoffRadius and cutoffMethod |
702 | > | * |
703 | > | * cutoffRadius : realType |
704 | > | * If the cutoffRadius was explicitly set, use that value. |
705 | > | * If the cutoffRadius was not explicitly set: |
706 | > | * Are there electrostatic atoms? Use 12.0 Angstroms. |
707 | > | * No electrostatic atoms? Poll the atom types present in the |
708 | > | * simulation for suggested cutoff values (e.g. 2.5 * sigma). |
709 | > | * Use the maximum suggested value that was found. |
710 | > | * |
711 | > | * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) |
712 | > | * If cutoffMethod was explicitly set, use that choice. |
713 | > | * If cutoffMethod was not explicitly set, use SHIFTED_FORCE |
714 | > | */ |
715 | > | void SimInfo::setupCutoffs() { |
716 | ||
717 | < | int useLennardJones = 0; |
718 | < | int useElectrostatic = 0; |
719 | < | int useEAM = 0; |
720 | < | int useSC = 0; |
721 | < | int useCharge = 0; |
722 | < | int useDirectional = 0; |
723 | < | int useDipole = 0; |
724 | < | int useGayBerne = 0; |
725 | < | int useSticky = 0; |
726 | < | int useStickyPower = 0; |
727 | < | int useShape = 0; |
728 | < | int useFLARB = 0; //it is not in AtomType yet |
729 | < | int useDirectionalAtom = 0; |
730 | < | int useElectrostatics = 0; |
731 | < | //usePBC and useRF are from simParams |
732 | < | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
733 | < | int useRF; |
734 | < | int useSF; |
735 | < | int useSP; |
736 | < | int useBoxDipole; |
717 | > | if (simParams_->haveCutoffRadius()) { |
718 | > | cutoffRadius_ = simParams_->getCutoffRadius(); |
719 | > | } else { |
720 | > | if (usesElectrostaticAtoms_) { |
721 | > | sprintf(painCave.errMsg, |
722 | > | "SimInfo: No value was set for the cutoffRadius.\n" |
723 | > | "\tOpenMD will use a default value of 12.0 angstroms" |
724 | > | "\tfor the cutoffRadius.\n"); |
725 | > | painCave.isFatal = 0; |
726 | > | painCave.severity = OPENMD_INFO; |
727 | > | simError(); |
728 | > | cutoffRadius_ = 12.0; |
729 | > | } else { |
730 | > | RealType thisCut; |
731 | > | set<AtomType*>::iterator i; |
732 | > | set<AtomType*> atomTypes; |
733 | > | atomTypes = getSimulatedAtomTypes(); |
734 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
735 | > | thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); |
736 | > | cutoffRadius_ = max(thisCut, cutoffRadius_); |
737 | > | } |
738 | > | sprintf(painCave.errMsg, |
739 | > | "SimInfo: No value was set for the cutoffRadius.\n" |
740 | > | "\tOpenMD will use %lf angstroms.\n", |
741 | > | cutoffRadius_); |
742 | > | painCave.isFatal = 0; |
743 | > | painCave.severity = OPENMD_INFO; |
744 | > | simError(); |
745 | > | } |
746 | > | } |
747 | ||
748 | < | std::string myMethod; |
748 | > | InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); |
749 | ||
750 | < | // set the useRF logical |
751 | < | useRF = 0; |
752 | < | useSF = 0; |
753 | < | useSP = 0; |
754 | < | useBoxDipole = 0; |
750 | > | map<string, CutoffMethod> stringToCutoffMethod; |
751 | > | stringToCutoffMethod["HARD"] = HARD; |
752 | > | stringToCutoffMethod["SWITCHING_FUNCTION"] = SWITCHING_FUNCTION; |
753 | > | stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; |
754 | > | stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; |
755 | > | |
756 | > | if (simParams_->haveCutoffMethod()) { |
757 | > | string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); |
758 | > | map<string, CutoffMethod>::iterator i; |
759 | > | i = stringToCutoffMethod.find(cutMeth); |
760 | > | if (i == stringToCutoffMethod.end()) { |
761 | > | sprintf(painCave.errMsg, |
762 | > | "SimInfo: Could not find chosen cutoffMethod %s\n" |
763 | > | "\tShould be one of: " |
764 | > | "HARD, SWITCHING_FUNCTION, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", |
765 | > | cutMeth.c_str()); |
766 | > | painCave.isFatal = 1; |
767 | > | painCave.severity = OPENMD_ERROR; |
768 | > | simError(); |
769 | > | } else { |
770 | > | cutoffMethod_ = i->second; |
771 | > | } |
772 | > | } else { |
773 | > | sprintf(painCave.errMsg, |
774 | > | "SimInfo: No value was set for the cutoffMethod.\n" |
775 | > | "\tOpenMD will use SHIFTED_FORCE.\n"); |
776 | > | painCave.isFatal = 0; |
777 | > | painCave.severity = OPENMD_INFO; |
778 | > | simError(); |
779 | > | cutoffMethod_ = SHIFTED_FORCE; |
780 | > | } |
781 | ||
782 | + | InteractionManager::Instance()->setCutoffMethod(cutoffMethod_); |
783 | + | } |
784 | + | |
785 | + | /** |
786 | + | * setupSwitching |
787 | + | * |
788 | + | * Sets the values of switchingRadius and |
789 | + | * If the switchingRadius was explicitly set, use that value (but check it) |
790 | + | * If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ |
791 | + | */ |
792 | + | void SimInfo::setupSwitching() { |
793 | + | |
794 | + | if (simParams_->haveSwitchingRadius()) { |
795 | + | switchingRadius_ = simParams_->getSwitchingRadius(); |
796 | + | if (switchingRadius_ > cutoffRadius_) { |
797 | + | sprintf(painCave.errMsg, |
798 | + | "SimInfo: switchingRadius (%f) is larger than cutoffRadius(%f)\n", |
799 | + | switchingRadius_, cutoffRadius_); |
800 | + | painCave.isFatal = 1; |
801 | + | painCave.severity = OPENMD_ERROR; |
802 | + | simError(); |
803 | + | } |
804 | + | } else { |
805 | + | switchingRadius_ = 0.85 * cutoffRadius_; |
806 | + | sprintf(painCave.errMsg, |
807 | + | "SimInfo: No value was set for the switchingRadius.\n" |
808 | + | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
809 | + | "\tswitchingRadius = %f. for this simulation\n", switchingRadius_); |
810 | + | painCave.isFatal = 0; |
811 | + | painCave.severity = OPENMD_WARNING; |
812 | + | simError(); |
813 | + | } |
814 | + | |
815 | + | InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); |
816 | ||
817 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
818 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
819 | < | toUpper(myMethod); |
820 | < | if (myMethod == "REACTION_FIELD"){ |
821 | < | useRF = 1; |
822 | < | } else if (myMethod == "SHIFTED_FORCE"){ |
823 | < | useSF = 1; |
824 | < | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
825 | < | useSP = 1; |
817 | > | SwitchingFunctionType ft; |
818 | > | |
819 | > | if (simParams_->haveSwitchingFunctionType()) { |
820 | > | string funcType = simParams_->getSwitchingFunctionType(); |
821 | > | toUpper(funcType); |
822 | > | if (funcType == "CUBIC") { |
823 | > | ft = cubic; |
824 | > | } else { |
825 | > | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
826 | > | ft = fifth_order_poly; |
827 | > | } else { |
828 | > | // throw error |
829 | > | sprintf( painCave.errMsg, |
830 | > | "SimInfo : Unknown switchingFunctionType. (Input file specified %s .)\n" |
831 | > | "\tswitchingFunctionType must be one of: " |
832 | > | "\"cubic\" or \"fifth_order_polynomial\".", |
833 | > | funcType.c_str() ); |
834 | > | painCave.isFatal = 1; |
835 | > | painCave.severity = OPENMD_ERROR; |
836 | > | simError(); |
837 | > | } |
838 | } | |
839 | } | |
840 | < | |
841 | < | if (simParams_->haveAccumulateBoxDipole()) |
842 | < | if (simParams_->getAccumulateBoxDipole()) |
780 | < | useBoxDipole = 1; |
840 | > | |
841 | > | InteractionManager::Instance()->setSwitchingFunctionType(ft); |
842 | > | } |
843 | ||
844 | + | /** |
845 | + | * setupSkinThickness |
846 | + | * |
847 | + | * If the skinThickness was explicitly set, use that value (but check it) |
848 | + | * If the skinThickness was not explicitly set: use 1.0 angstroms |
849 | + | */ |
850 | + | void SimInfo::setupSkinThickness() { |
851 | + | if (simParams_->haveSkinThickness()) { |
852 | + | skinThickness_ = simParams_->getSkinThickness(); |
853 | + | } else { |
854 | + | skinThickness_ = 1.0; |
855 | + | sprintf(painCave.errMsg, |
856 | + | "SimInfo Warning: No value was set for the skinThickness.\n" |
857 | + | "\tOpenMD will use a default value of %f Angstroms\n" |
858 | + | "\tfor this simulation\n", skinThickness_); |
859 | + | painCave.isFatal = 0; |
860 | + | simError(); |
861 | + | } |
862 | + | } |
863 | + | |
864 | + | void SimInfo::setupSimType() { |
865 | + | set<AtomType*>::iterator i; |
866 | + | set<AtomType*> atomTypes; |
867 | + | atomTypes = getSimulatedAtomTypes(); |
868 | + | |
869 | useAtomicVirial_ = simParams_->getUseAtomicVirial(); | |
870 | ||
871 | + | int usesElectrostatic = 0; |
872 | + | int usesMetallic = 0; |
873 | + | int usesDirectional = 0; |
874 | //loop over all of the atom types | |
875 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
876 | < | useLennardJones |= (*i)->isLennardJones(); |
877 | < | useElectrostatic |= (*i)->isElectrostatic(); |
878 | < | useEAM |= (*i)->isEAM(); |
789 | < | useSC |= (*i)->isSC(); |
790 | < | useCharge |= (*i)->isCharge(); |
791 | < | useDirectional |= (*i)->isDirectional(); |
792 | < | useDipole |= (*i)->isDipole(); |
793 | < | useGayBerne |= (*i)->isGayBerne(); |
794 | < | useSticky |= (*i)->isSticky(); |
795 | < | useStickyPower |= (*i)->isStickyPower(); |
796 | < | useShape |= (*i)->isShape(); |
876 | > | usesElectrostatic |= (*i)->isElectrostatic(); |
877 | > | usesMetallic |= (*i)->isMetal(); |
878 | > | usesDirectional |= (*i)->isDirectional(); |
879 | } | |
880 | ||
799 | – | if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { |
800 | – | useDirectionalAtom = 1; |
801 | – | } |
802 | – | |
803 | – | if (useCharge || useDipole) { |
804 | – | useElectrostatics = 1; |
805 | – | } |
806 | – | |
881 | #ifdef IS_MPI | |
882 | int temp; | |
883 | + | temp = usesDirectional; |
884 | + | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
885 | ||
886 | < | temp = usePBC; |
887 | < | MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
886 | > | temp = usesMetallic; |
887 | > | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
888 | ||
889 | < | temp = useDirectionalAtom; |
890 | < | MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
815 | < | |
816 | < | temp = useLennardJones; |
817 | < | MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
818 | < | |
819 | < | temp = useElectrostatics; |
820 | < | MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
821 | < | |
822 | < | temp = useCharge; |
823 | < | MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
824 | < | |
825 | < | temp = useDipole; |
826 | < | MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
827 | < | |
828 | < | temp = useSticky; |
829 | < | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
830 | < | |
831 | < | temp = useStickyPower; |
832 | < | MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
833 | < | |
834 | < | temp = useGayBerne; |
835 | < | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
836 | < | |
837 | < | temp = useEAM; |
838 | < | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
839 | < | |
840 | < | temp = useSC; |
841 | < | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
842 | < | |
843 | < | temp = useShape; |
844 | < | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
845 | < | |
846 | < | temp = useFLARB; |
847 | < | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
848 | < | |
849 | < | temp = useRF; |
850 | < | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
851 | < | |
852 | < | temp = useSF; |
853 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
854 | < | |
855 | < | temp = useSP; |
856 | < | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
857 | < | |
858 | < | temp = useBoxDipole; |
859 | < | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
860 | < | |
861 | < | temp = useAtomicVirial_; |
862 | < | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
863 | < | |
889 | > | temp = usesElectrostatic; |
890 | > | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
891 | #endif | |
892 | < | |
893 | < | fInfo_.SIM_uses_PBC = usePBC; |
894 | < | fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; |
895 | < | fInfo_.SIM_uses_LennardJones = useLennardJones; |
896 | < | fInfo_.SIM_uses_Electrostatics = useElectrostatics; |
897 | < | fInfo_.SIM_uses_Charges = useCharge; |
871 | < | fInfo_.SIM_uses_Dipoles = useDipole; |
872 | < | fInfo_.SIM_uses_Sticky = useSticky; |
873 | < | fInfo_.SIM_uses_StickyPower = useStickyPower; |
874 | < | fInfo_.SIM_uses_GayBerne = useGayBerne; |
875 | < | fInfo_.SIM_uses_EAM = useEAM; |
876 | < | fInfo_.SIM_uses_SC = useSC; |
877 | < | fInfo_.SIM_uses_Shapes = useShape; |
878 | < | fInfo_.SIM_uses_FLARB = useFLARB; |
879 | < | fInfo_.SIM_uses_RF = useRF; |
880 | < | fInfo_.SIM_uses_SF = useSF; |
881 | < | fInfo_.SIM_uses_SP = useSP; |
882 | < | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
883 | < | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
892 | > | fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; |
893 | > | fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; |
894 | > | fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; |
895 | > | fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; |
896 | > | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; |
897 | > | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; |
898 | } | |
899 | ||
900 | void SimInfo::setupFortranSim() { | |
901 | int isError; | |
902 | int nExclude, nOneTwo, nOneThree, nOneFour; | |
903 | < | std::vector<int> fortranGlobalGroupMembership; |
903 | > | vector<int> fortranGlobalGroupMembership; |
904 | ||
905 | + | notifyFortranSkinThickness(&skinThickness_); |
906 | + | |
907 | + | int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
908 | + | int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
909 | + | notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); |
910 | + | |
911 | isError = 0; | |
912 | ||
913 | //globalGroupMembership_ is filled by SimCreator | |
# | Line 896 | Line 916 | namespace OpenMD { | |
916 | } | |
917 | ||
918 | //calculate mass ratio of cutoff group | |
919 | < | std::vector<RealType> mfact; |
919 | > | vector<RealType> mfact; |
920 | SimInfo::MoleculeIterator mi; | |
921 | Molecule* mol; | |
922 | Molecule::CutoffGroupIterator ci; | |
# | Line 923 | Line 943 | namespace OpenMD { | |
943 | } | |
944 | ||
945 | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) | |
946 | < | std::vector<int> identArray; |
946 | > | vector<int> identArray; |
947 | ||
948 | //to avoid memory reallocation, reserve enough space identArray | |
949 | identArray.reserve(getNAtoms()); | |
# | Line 936 | Line 956 | namespace OpenMD { | |
956 | ||
957 | //fill molMembershipArray | |
958 | //molMembershipArray is filled by SimCreator | |
959 | < | std::vector<int> molMembershipArray(nGlobalAtoms_); |
959 | > | vector<int> molMembershipArray(nGlobalAtoms_); |
960 | for (int i = 0; i < nGlobalAtoms_; i++) { | |
961 | molMembershipArray[i] = globalMolMembership_[i] + 1; | |
962 | } | |
# | Line 989 | Line 1009 | namespace OpenMD { | |
1009 | void SimInfo::setupFortranParallel() { | |
1010 | #ifdef IS_MPI | |
1011 | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex | |
1012 | < | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
1013 | < | std::vector<int> localToGlobalCutoffGroupIndex; |
1012 | > | vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
1013 | > | vector<int> localToGlobalCutoffGroupIndex; |
1014 | SimInfo::MoleculeIterator mi; | |
1015 | Molecule::AtomIterator ai; | |
1016 | Molecule::CutoffGroupIterator ci; | |
# | Line 1042 | Line 1062 | namespace OpenMD { | |
1062 | #endif | |
1063 | } | |
1064 | ||
1045 | – | void SimInfo::setupCutoff() { |
1046 | – | |
1047 | – | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
1065 | ||
1049 | – | // Check the cutoff policy |
1050 | – | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default |
1051 | – | |
1052 | – | // Set LJ shifting bools to false |
1053 | – | ljsp_ = 0; |
1054 | – | ljsf_ = 0; |
1055 | – | |
1056 | – | std::string myPolicy; |
1057 | – | if (forceFieldOptions_.haveCutoffPolicy()){ |
1058 | – | myPolicy = forceFieldOptions_.getCutoffPolicy(); |
1059 | – | }else if (simParams_->haveCutoffPolicy()) { |
1060 | – | myPolicy = simParams_->getCutoffPolicy(); |
1061 | – | } |
1062 | – | |
1063 | – | if (!myPolicy.empty()){ |
1064 | – | toUpper(myPolicy); |
1065 | – | if (myPolicy == "MIX") { |
1066 | – | cp = MIX_CUTOFF_POLICY; |
1067 | – | } else { |
1068 | – | if (myPolicy == "MAX") { |
1069 | – | cp = MAX_CUTOFF_POLICY; |
1070 | – | } else { |
1071 | – | if (myPolicy == "TRADITIONAL") { |
1072 | – | cp = TRADITIONAL_CUTOFF_POLICY; |
1073 | – | } else { |
1074 | – | // throw error |
1075 | – | sprintf( painCave.errMsg, |
1076 | – | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
1077 | – | painCave.isFatal = 1; |
1078 | – | simError(); |
1079 | – | } |
1080 | – | } |
1081 | – | } |
1082 | – | } |
1083 | – | notifyFortranCutoffPolicy(&cp); |
1084 | – | |
1085 | – | // Check the Skin Thickness for neighborlists |
1086 | – | RealType skin; |
1087 | – | if (simParams_->haveSkinThickness()) { |
1088 | – | skin = simParams_->getSkinThickness(); |
1089 | – | notifyFortranSkinThickness(&skin); |
1090 | – | } |
1091 | – | |
1092 | – | // Check if the cutoff was set explicitly: |
1093 | – | if (simParams_->haveCutoffRadius()) { |
1094 | – | rcut_ = simParams_->getCutoffRadius(); |
1095 | – | if (simParams_->haveSwitchingRadius()) { |
1096 | – | rsw_ = simParams_->getSwitchingRadius(); |
1097 | – | } else { |
1098 | – | if (fInfo_.SIM_uses_Charges | |
1099 | – | fInfo_.SIM_uses_Dipoles | |
1100 | – | fInfo_.SIM_uses_RF) { |
1101 | – | |
1102 | – | rsw_ = 0.85 * rcut_; |
1103 | – | sprintf(painCave.errMsg, |
1104 | – | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1105 | – | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
1106 | – | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1107 | – | painCave.isFatal = 0; |
1108 | – | simError(); |
1109 | – | } else { |
1110 | – | rsw_ = rcut_; |
1111 | – | sprintf(painCave.errMsg, |
1112 | – | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1113 | – | "\tOpenMD will use the same value as the cutoffRadius.\n" |
1114 | – | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1115 | – | painCave.isFatal = 0; |
1116 | – | simError(); |
1117 | – | } |
1118 | – | } |
1119 | – | |
1120 | – | if (simParams_->haveElectrostaticSummationMethod()) { |
1121 | – | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1122 | – | toUpper(myMethod); |
1123 | – | |
1124 | – | if (myMethod == "SHIFTED_POTENTIAL") { |
1125 | – | ljsp_ = 1; |
1126 | – | } else if (myMethod == "SHIFTED_FORCE") { |
1127 | – | ljsf_ = 1; |
1128 | – | } |
1129 | – | } |
1130 | – | |
1131 | – | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1132 | – | |
1133 | – | } else { |
1134 | – | |
1135 | – | // For electrostatic atoms, we'll assume a large safe value: |
1136 | – | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
1137 | – | sprintf(painCave.errMsg, |
1138 | – | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
1139 | – | "\tOpenMD will use a default value of 15.0 angstroms" |
1140 | – | "\tfor the cutoffRadius.\n"); |
1141 | – | painCave.isFatal = 0; |
1142 | – | simError(); |
1143 | – | rcut_ = 15.0; |
1144 | – | |
1145 | – | if (simParams_->haveElectrostaticSummationMethod()) { |
1146 | – | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1147 | – | toUpper(myMethod); |
1148 | – | |
1149 | – | // For the time being, we're tethering the LJ shifted behavior to the |
1150 | – | // electrostaticSummationMethod keyword options |
1151 | – | if (myMethod == "SHIFTED_POTENTIAL") { |
1152 | – | ljsp_ = 1; |
1153 | – | } else if (myMethod == "SHIFTED_FORCE") { |
1154 | – | ljsf_ = 1; |
1155 | – | } |
1156 | – | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1157 | – | if (simParams_->haveSwitchingRadius()){ |
1158 | – | sprintf(painCave.errMsg, |
1159 | – | "SimInfo Warning: A value was set for the switchingRadius\n" |
1160 | – | "\teven though the electrostaticSummationMethod was\n" |
1161 | – | "\tset to %s\n", myMethod.c_str()); |
1162 | – | painCave.isFatal = 1; |
1163 | – | simError(); |
1164 | – | } |
1165 | – | } |
1166 | – | } |
1167 | – | |
1168 | – | if (simParams_->haveSwitchingRadius()){ |
1169 | – | rsw_ = simParams_->getSwitchingRadius(); |
1170 | – | } else { |
1171 | – | sprintf(painCave.errMsg, |
1172 | – | "SimCreator Warning: No value was set for switchingRadius.\n" |
1173 | – | "\tOpenMD will use a default value of\n" |
1174 | – | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
1175 | – | painCave.isFatal = 0; |
1176 | – | simError(); |
1177 | – | rsw_ = 0.85 * rcut_; |
1178 | – | } |
1179 | – | |
1180 | – | Electrostatic::setElectrostaticCutoffRadius(rcut_, rsw_); |
1181 | – | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1182 | – | |
1183 | – | } else { |
1184 | – | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so |
1185 | – | // We'll punt and let fortran figure out the cutoffs later. |
1186 | – | |
1187 | – | notifyFortranYouAreOnYourOwn(); |
1188 | – | |
1189 | – | } |
1190 | – | } |
1191 | – | } |
1192 | – | |
1193 | – | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
1194 | – | |
1195 | – | int errorOut; |
1196 | – | ElectrostaticSummationMethod esm = NONE; |
1197 | – | ElectrostaticScreeningMethod sm = UNDAMPED; |
1198 | – | RealType alphaVal; |
1199 | – | RealType dielectric; |
1200 | – | |
1201 | – | errorOut = isError; |
1202 | – | |
1203 | – | if (simParams_->haveElectrostaticSummationMethod()) { |
1204 | – | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1205 | – | toUpper(myMethod); |
1206 | – | if (myMethod == "NONE") { |
1207 | – | esm = NONE; |
1208 | – | } else { |
1209 | – | if (myMethod == "SWITCHING_FUNCTION") { |
1210 | – | esm = SWITCHING_FUNCTION; |
1211 | – | } else { |
1212 | – | if (myMethod == "SHIFTED_POTENTIAL") { |
1213 | – | esm = SHIFTED_POTENTIAL; |
1214 | – | } else { |
1215 | – | if (myMethod == "SHIFTED_FORCE") { |
1216 | – | esm = SHIFTED_FORCE; |
1217 | – | } else { |
1218 | – | if (myMethod == "REACTION_FIELD") { |
1219 | – | esm = REACTION_FIELD; |
1220 | – | dielectric = simParams_->getDielectric(); |
1221 | – | if (!simParams_->haveDielectric()) { |
1222 | – | // throw warning |
1223 | – | sprintf( painCave.errMsg, |
1224 | – | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1225 | – | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1226 | – | painCave.isFatal = 0; |
1227 | – | simError(); |
1228 | – | } |
1229 | – | } else { |
1230 | – | // throw error |
1231 | – | sprintf( painCave.errMsg, |
1232 | – | "SimInfo error: Unknown electrostaticSummationMethod.\n" |
1233 | – | "\t(Input file specified %s .)\n" |
1234 | – | "\telectrostaticSummationMethod must be one of: \"none\",\n" |
1235 | – | "\t\"shifted_potential\", \"shifted_force\", or \n" |
1236 | – | "\t\"reaction_field\".\n", myMethod.c_str() ); |
1237 | – | painCave.isFatal = 1; |
1238 | – | simError(); |
1239 | – | } |
1240 | – | } |
1241 | – | } |
1242 | – | } |
1243 | – | } |
1244 | – | } |
1245 | – | |
1246 | – | if (simParams_->haveElectrostaticScreeningMethod()) { |
1247 | – | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
1248 | – | toUpper(myScreen); |
1249 | – | if (myScreen == "UNDAMPED") { |
1250 | – | sm = UNDAMPED; |
1251 | – | } else { |
1252 | – | if (myScreen == "DAMPED") { |
1253 | – | sm = DAMPED; |
1254 | – | if (!simParams_->haveDampingAlpha()) { |
1255 | – | // first set a cutoff dependent alpha value |
1256 | – | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1257 | – | alphaVal = 0.5125 - rcut_* 0.025; |
1258 | – | // for values rcut > 20.5, alpha is zero |
1259 | – | if (alphaVal < 0) alphaVal = 0; |
1260 | – | |
1261 | – | // throw warning |
1262 | – | sprintf( painCave.errMsg, |
1263 | – | "SimInfo warning: dampingAlpha was not specified in the input file.\n" |
1264 | – | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1265 | – | painCave.isFatal = 0; |
1266 | – | simError(); |
1267 | – | } else { |
1268 | – | alphaVal = simParams_->getDampingAlpha(); |
1269 | – | } |
1270 | – | |
1271 | – | } else { |
1272 | – | // throw error |
1273 | – | sprintf( painCave.errMsg, |
1274 | – | "SimInfo error: Unknown electrostaticScreeningMethod.\n" |
1275 | – | "\t(Input file specified %s .)\n" |
1276 | – | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" |
1277 | – | "or \"damped\".\n", myScreen.c_str() ); |
1278 | – | painCave.isFatal = 1; |
1279 | – | simError(); |
1280 | – | } |
1281 | – | } |
1282 | – | } |
1283 | – | |
1284 | – | |
1285 | – | Electrostatic::setElectrostaticSummationMethod( esm ); |
1286 | – | Electrostatic::setElectrostaticScreeningMethod( sm ); |
1287 | – | Electrostatic::setDampingAlpha( alphaVal ); |
1288 | – | Electrostatic::setReactionFieldDielectric( dielectric ); |
1289 | – | initFortranFF( &errorOut ); |
1290 | – | } |
1291 | – | |
1066 | void SimInfo::setupSwitchingFunction() { | |
1293 | – | int ft = CUBIC; |
1067 | ||
1295 | – | if (simParams_->haveSwitchingFunctionType()) { |
1296 | – | std::string funcType = simParams_->getSwitchingFunctionType(); |
1297 | – | toUpper(funcType); |
1298 | – | if (funcType == "CUBIC") { |
1299 | – | ft = CUBIC; |
1300 | – | } else { |
1301 | – | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1302 | – | ft = FIFTH_ORDER_POLY; |
1303 | – | } else { |
1304 | – | // throw error |
1305 | – | sprintf( painCave.errMsg, |
1306 | – | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1307 | – | painCave.isFatal = 1; |
1308 | – | simError(); |
1309 | – | } |
1310 | – | } |
1311 | – | } |
1312 | – | |
1313 | – | // send switching function notification to switcheroo |
1314 | – | setFunctionType(&ft); |
1315 | – | |
1068 | } | |
1069 | ||
1070 | void SimInfo::setupAccumulateBoxDipole() { | |
# | Line 1320 | Line 1072 | namespace OpenMD { | |
1072 | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | |
1073 | if ( simParams_->haveAccumulateBoxDipole() ) | |
1074 | if ( simParams_->getAccumulateBoxDipole() ) { | |
1323 | – | setAccumulateBoxDipole(); |
1075 | calcBoxDipole_ = true; | |
1076 | } | |
1077 | ||
# | Line 1330 | Line 1081 | namespace OpenMD { | |
1081 | properties_.addProperty(genData); | |
1082 | } | |
1083 | ||
1084 | < | void SimInfo::removeProperty(const std::string& propName) { |
1084 | > | void SimInfo::removeProperty(const string& propName) { |
1085 | properties_.removeProperty(propName); | |
1086 | } | |
1087 | ||
# | Line 1338 | Line 1089 | namespace OpenMD { | |
1089 | properties_.clearProperties(); | |
1090 | } | |
1091 | ||
1092 | < | std::vector<std::string> SimInfo::getPropertyNames() { |
1092 | > | vector<string> SimInfo::getPropertyNames() { |
1093 | return properties_.getPropertyNames(); | |
1094 | } | |
1095 | ||
1096 | < | std::vector<GenericData*> SimInfo::getProperties() { |
1096 | > | vector<GenericData*> SimInfo::getProperties() { |
1097 | return properties_.getProperties(); | |
1098 | } | |
1099 | ||
1100 | < | GenericData* SimInfo::getPropertyByName(const std::string& propName) { |
1100 | > | GenericData* SimInfo::getPropertyByName(const string& propName) { |
1101 | return properties_.getPropertyByName(propName); | |
1102 | } | |
1103 | ||
# | Line 1429 | Line 1180 | namespace OpenMD { | |
1180 | ||
1181 | } | |
1182 | ||
1183 | < | std::ostream& operator <<(std::ostream& o, SimInfo& info) { |
1183 | > | ostream& operator <<(ostream& o, SimInfo& info) { |
1184 | ||
1185 | return o; | |
1186 | } | |
# | Line 1472 | Line 1223 | namespace OpenMD { | |
1223 | ||
1224 | ||
1225 | [ Ixx -Ixy -Ixz ] | |
1226 | < | J =| -Iyx Iyy -Iyz | |
1226 | > | J =| -Iyx Iyy -Iyz | |
1227 | [ -Izx -Iyz Izz ] | |
1228 | */ | |
1229 | ||
# | Line 1579 | Line 1330 | namespace OpenMD { | |
1330 | return IOIndexToIntegrableObject.at(index); | |
1331 | } | |
1332 | ||
1333 | < | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1333 | > | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { |
1334 | IOIndexToIntegrableObject= v; | |
1335 | } | |
1336 | ||
# | Line 1621 | Line 1372 | namespace OpenMD { | |
1372 | return; | |
1373 | } | |
1374 | /* | |
1375 | < | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1375 | > | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { |
1376 | assert( v.size() == nAtoms_ + nRigidBodies_); | |
1377 | sdByGlobalIndex_ = v; | |
1378 | } | |
# | Line 1631 | Line 1382 | namespace OpenMD { | |
1382 | return sdByGlobalIndex_.at(index); | |
1383 | } | |
1384 | */ | |
1385 | + | int SimInfo::getNGlobalConstraints() { |
1386 | + | int nGlobalConstraints; |
1387 | + | #ifdef IS_MPI |
1388 | + | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
1389 | + | MPI_COMM_WORLD); |
1390 | + | #else |
1391 | + | nGlobalConstraints = nConstraints_; |
1392 | + | #endif |
1393 | + | return nGlobalConstraints; |
1394 | + | } |
1395 | + | |
1396 | }//end namespace OpenMD | |
1397 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |