# | 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" |
59 | – | #include "UseTheForce/doForces_interface.h" |
60 | – | #include "UseTheForce/DarkSide/neighborLists_interface.h" |
61 | – | #include "UseTheForce/DarkSide/switcheroo_interface.h" |
57 | #include "utils/MemoryUtils.hpp" | |
58 | #include "utils/simError.h" | |
59 | #include "selection/SelectionManager.hpp" | |
60 | #include "io/ForceFieldOptions.hpp" | |
61 | #include "UseTheForce/ForceField.hpp" | |
62 | < | |
68 | < | |
62 | > | #include "nonbonded/SwitchingFunction.hpp" |
63 | #ifdef IS_MPI | |
64 | < | #include "UseTheForce/mpiComponentPlan.h" |
65 | < | #include "UseTheForce/DarkSide/simParallel_interface.h" |
72 | < | #endif |
64 | > | #include <mpi.h> |
65 | > | #endif |
66 | ||
67 | + | using namespace std; |
68 | 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 | – | } |
69 | ||
70 | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | |
71 | forceField_(ff), simParams_(simParams), | |
# | Line 89 | Line 74 | namespace OpenMD { | |
74 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
75 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), | |
76 | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | |
77 | < | nConstraints_(0), sman_(NULL), fortranInitialized_(false), |
78 | < | calcBoxDipole_(false), useAtomicVirial_(true) { |
79 | < | |
80 | < | |
81 | < | MoleculeStamp* molStamp; |
82 | < | int nMolWithSameStamp; |
83 | < | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
84 | < | int nGroups = 0; //total cutoff groups defined in meta-data file |
85 | < | CutoffGroupStamp* cgStamp; |
86 | < | RigidBodyStamp* rbStamp; |
87 | < | int nRigidAtoms = 0; |
88 | < | |
89 | < | std::vector<Component*> components = simParams->getComponents(); |
77 | > | nConstraints_(0), sman_(NULL), topologyDone_(false), |
78 | > | calcBoxDipole_(false), useAtomicVirial_(true) { |
79 | > | |
80 | > | MoleculeStamp* molStamp; |
81 | > | int nMolWithSameStamp; |
82 | > | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
83 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
84 | > | CutoffGroupStamp* cgStamp; |
85 | > | RigidBodyStamp* rbStamp; |
86 | > | int nRigidAtoms = 0; |
87 | > | |
88 | > | vector<Component*> components = simParams->getComponents(); |
89 | > | |
90 | > | for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
91 | > | molStamp = (*i)->getMoleculeStamp(); |
92 | > | nMolWithSameStamp = (*i)->getNMol(); |
93 | ||
94 | < | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
95 | < | molStamp = (*i)->getMoleculeStamp(); |
96 | < | nMolWithSameStamp = (*i)->getNMol(); |
97 | < | |
98 | < | addMoleculeStamp(molStamp, nMolWithSameStamp); |
99 | < | |
100 | < | //calculate atoms in molecules |
101 | < | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
102 | < | |
103 | < | //calculate atoms in cutoff groups |
104 | < | int nAtomsInGroups = 0; |
105 | < | 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 | < | |
94 | > | addMoleculeStamp(molStamp, nMolWithSameStamp); |
95 | > | |
96 | > | //calculate atoms in molecules |
97 | > | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
98 | > | |
99 | > | //calculate atoms in cutoff groups |
100 | > | int nAtomsInGroups = 0; |
101 | > | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); |
102 | > | |
103 | > | for (int j=0; j < nCutoffGroupsInStamp; j++) { |
104 | > | cgStamp = molStamp->getCutoffGroupStamp(j); |
105 | > | nAtomsInGroups += cgStamp->getNMembers(); |
106 | } | |
107 | < | |
108 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff |
109 | < | //group therefore the total number of cutoff groups in the system is |
110 | < | //equal to the total number of atoms minus number of atoms belong to |
111 | < | //cutoff group defined in meta-data file plus the number of cutoff |
112 | < | //groups defined in meta-data file |
113 | < | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
114 | < | |
115 | < | //every free atom (atom does not belong to rigid bodies) is an |
116 | < | //integrable object therefore the total number of integrable objects |
117 | < | //in the system is equal to the total number of atoms minus number of |
118 | < | //atoms belong to rigid body defined in meta-data file plus the number |
119 | < | //of rigid bodies defined in meta-data file |
120 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
121 | < | + nGlobalRigidBodies_; |
122 | < | |
123 | < | nGlobalMols_ = molStampIds_.size(); |
158 | < | molToProcMap_.resize(nGlobalMols_); |
107 | > | |
108 | > | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; |
109 | > | |
110 | > | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; |
111 | > | |
112 | > | //calculate atoms in rigid bodies |
113 | > | int nAtomsInRigidBodies = 0; |
114 | > | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); |
115 | > | |
116 | > | for (int j=0; j < nRigidBodiesInStamp; j++) { |
117 | > | rbStamp = molStamp->getRigidBodyStamp(j); |
118 | > | nAtomsInRigidBodies += rbStamp->getNMembers(); |
119 | > | } |
120 | > | |
121 | > | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; |
122 | > | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; |
123 | > | |
124 | } | |
125 | + | |
126 | + | //every free atom (atom does not belong to cutoff groups) is a cutoff |
127 | + | //group therefore the total number of cutoff groups in the system is |
128 | + | //equal to the total number of atoms minus number of atoms belong to |
129 | + | //cutoff group defined in meta-data file plus the number of cutoff |
130 | + | //groups defined in meta-data file |
131 | ||
132 | + | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
133 | + | |
134 | + | //every free atom (atom does not belong to rigid bodies) is an |
135 | + | //integrable object therefore the total number of integrable objects |
136 | + | //in the system is equal to the total number of atoms minus number of |
137 | + | //atoms belong to rigid body defined in meta-data file plus the number |
138 | + | //of rigid bodies defined in meta-data file |
139 | + | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
140 | + | + nGlobalRigidBodies_; |
141 | + | |
142 | + | nGlobalMols_ = molStampIds_.size(); |
143 | + | molToProcMap_.resize(nGlobalMols_); |
144 | + | } |
145 | + | |
146 | SimInfo::~SimInfo() { | |
147 | < | std::map<int, Molecule*>::iterator i; |
147 | > | map<int, Molecule*>::iterator i; |
148 | for (i = molecules_.begin(); i != molecules_.end(); ++i) { | |
149 | delete i->second; | |
150 | } | |
# | Line 170 | Line 155 | namespace OpenMD { | |
155 | delete forceField_; | |
156 | } | |
157 | ||
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 | – | } |
158 | ||
159 | bool SimInfo::addMolecule(Molecule* mol) { | |
160 | MoleculeIterator i; | |
161 | < | |
161 | > | |
162 | i = molecules_.find(mol->getGlobalIndex()); | |
163 | if (i == molecules_.end() ) { | |
164 | < | |
165 | < | molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); |
166 | < | |
164 | > | |
165 | > | molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); |
166 | > | |
167 | nAtoms_ += mol->getNAtoms(); | |
168 | nBonds_ += mol->getNBonds(); | |
169 | nBends_ += mol->getNBends(); | |
# | Line 198 | Line 173 | namespace OpenMD { | |
173 | nIntegrableObjects_ += mol->getNIntegrableObjects(); | |
174 | nCutoffGroups_ += mol->getNCutoffGroups(); | |
175 | nConstraints_ += mol->getNConstraintPairs(); | |
176 | < | |
176 | > | |
177 | addInteractionPairs(mol); | |
178 | < | |
178 | > | |
179 | return true; | |
180 | } else { | |
181 | return false; | |
182 | } | |
183 | } | |
184 | < | |
184 | > | |
185 | bool SimInfo::removeMolecule(Molecule* mol) { | |
186 | MoleculeIterator i; | |
187 | i = molecules_.find(mol->getGlobalIndex()); | |
# | Line 234 | Line 209 | namespace OpenMD { | |
209 | } else { | |
210 | return false; | |
211 | } | |
237 | – | |
238 | – | |
212 | } | |
213 | ||
214 | ||
# | Line 253 | Line 226 | namespace OpenMD { | |
226 | void SimInfo::calcNdf() { | |
227 | int ndf_local; | |
228 | MoleculeIterator i; | |
229 | < | std::vector<StuntDouble*>::iterator j; |
229 | > | vector<StuntDouble*>::iterator j; |
230 | Molecule* mol; | |
231 | StuntDouble* integrableObject; | |
232 | ||
# | Line 299 | Line 272 | namespace OpenMD { | |
272 | #endif | |
273 | return fdf_; | |
274 | } | |
275 | + | |
276 | + | unsigned int SimInfo::getNLocalCutoffGroups(){ |
277 | + | int nLocalCutoffAtoms = 0; |
278 | + | Molecule* mol; |
279 | + | MoleculeIterator mi; |
280 | + | CutoffGroup* cg; |
281 | + | Molecule::CutoffGroupIterator ci; |
282 | ||
283 | + | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
284 | + | |
285 | + | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
286 | + | cg = mol->nextCutoffGroup(ci)) { |
287 | + | nLocalCutoffAtoms += cg->getNumAtom(); |
288 | + | |
289 | + | } |
290 | + | } |
291 | + | |
292 | + | return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; |
293 | + | } |
294 | + | |
295 | void SimInfo::calcNdfRaw() { | |
296 | int ndfRaw_local; | |
297 | ||
298 | MoleculeIterator i; | |
299 | < | std::vector<StuntDouble*>::iterator j; |
299 | > | vector<StuntDouble*>::iterator j; |
300 | Molecule* mol; | |
301 | StuntDouble* integrableObject; | |
302 | ||
# | Line 353 | Line 345 | namespace OpenMD { | |
345 | ||
346 | void SimInfo::addInteractionPairs(Molecule* mol) { | |
347 | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | |
348 | < | std::vector<Bond*>::iterator bondIter; |
349 | < | std::vector<Bend*>::iterator bendIter; |
350 | < | std::vector<Torsion*>::iterator torsionIter; |
351 | < | std::vector<Inversion*>::iterator inversionIter; |
348 | > | vector<Bond*>::iterator bondIter; |
349 | > | vector<Bend*>::iterator bendIter; |
350 | > | vector<Torsion*>::iterator torsionIter; |
351 | > | vector<Inversion*>::iterator inversionIter; |
352 | Bond* bond; | |
353 | Bend* bend; | |
354 | Torsion* torsion; | |
# | Line 374 | Line 366 | namespace OpenMD { | |
366 | // always be excluded. These are done at the bottom of this | |
367 | // function. | |
368 | ||
369 | < | std::map<int, std::set<int> > atomGroups; |
369 | > | map<int, set<int> > atomGroups; |
370 | Molecule::RigidBodyIterator rbIter; | |
371 | RigidBody* rb; | |
372 | Molecule::IntegrableObjectIterator ii; | |
# | Line 386 | Line 378 | namespace OpenMD { | |
378 | ||
379 | if (integrableObject->isRigidBody()) { | |
380 | rb = static_cast<RigidBody*>(integrableObject); | |
381 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
382 | < | std::set<int> rigidAtoms; |
381 | > | vector<Atom*> atoms = rb->getAtoms(); |
382 | > | set<int> rigidAtoms; |
383 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
384 | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | |
385 | } | |
386 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
387 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
387 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
388 | } | |
389 | } else { | |
390 | < | std::set<int> oneAtomSet; |
390 | > | set<int> oneAtomSet; |
391 | oneAtomSet.insert(integrableObject->getGlobalIndex()); | |
392 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
392 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
393 | } | |
394 | } | |
395 | ||
# | Line 500 | Line 492 | namespace OpenMD { | |
492 | ||
493 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | |
494 | rb = mol->nextRigidBody(rbIter)) { | |
495 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
495 | > | vector<Atom*> atoms = rb->getAtoms(); |
496 | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | |
497 | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | |
498 | a = atoms[i]->getGlobalIndex(); | |
# | Line 514 | Line 506 | namespace OpenMD { | |
506 | ||
507 | void SimInfo::removeInteractionPairs(Molecule* mol) { | |
508 | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | |
509 | < | std::vector<Bond*>::iterator bondIter; |
510 | < | std::vector<Bend*>::iterator bendIter; |
511 | < | std::vector<Torsion*>::iterator torsionIter; |
512 | < | std::vector<Inversion*>::iterator inversionIter; |
509 | > | vector<Bond*>::iterator bondIter; |
510 | > | vector<Bend*>::iterator bendIter; |
511 | > | vector<Torsion*>::iterator torsionIter; |
512 | > | vector<Inversion*>::iterator inversionIter; |
513 | Bond* bond; | |
514 | Bend* bend; | |
515 | Torsion* torsion; | |
# | Line 527 | Line 519 | namespace OpenMD { | |
519 | int c; | |
520 | int d; | |
521 | ||
522 | < | std::map<int, std::set<int> > atomGroups; |
522 | > | map<int, set<int> > atomGroups; |
523 | Molecule::RigidBodyIterator rbIter; | |
524 | RigidBody* rb; | |
525 | Molecule::IntegrableObjectIterator ii; | |
# | Line 539 | Line 531 | namespace OpenMD { | |
531 | ||
532 | if (integrableObject->isRigidBody()) { | |
533 | rb = static_cast<RigidBody*>(integrableObject); | |
534 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
535 | < | std::set<int> rigidAtoms; |
534 | > | vector<Atom*> atoms = rb->getAtoms(); |
535 | > | set<int> rigidAtoms; |
536 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
537 | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | |
538 | } | |
539 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
540 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
540 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
541 | } | |
542 | } else { | |
543 | < | std::set<int> oneAtomSet; |
543 | > | set<int> oneAtomSet; |
544 | oneAtomSet.insert(integrableObject->getGlobalIndex()); | |
545 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
545 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
546 | } | |
547 | } | |
548 | ||
# | Line 653 | Line 645 | namespace OpenMD { | |
645 | ||
646 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | |
647 | rb = mol->nextRigidBody(rbIter)) { | |
648 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
648 | > | vector<Atom*> atoms = rb->getAtoms(); |
649 | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | |
650 | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | |
651 | a = atoms[i]->getGlobalIndex(); | |
# | Line 676 | Line 668 | namespace OpenMD { | |
668 | molStampIds_.insert(molStampIds_.end(), nmol, curStampId); | |
669 | } | |
670 | ||
679 | – | void SimInfo::update() { |
671 | ||
672 | < | setupSimType(); |
673 | < | |
674 | < | #ifdef IS_MPI |
675 | < | setupFortranParallel(); |
676 | < | #endif |
677 | < | |
678 | < | setupFortranSim(); |
679 | < | |
680 | < | //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 | < | |
672 | > | /** |
673 | > | * update |
674 | > | * |
675 | > | * Performs the global checks and variable settings after the |
676 | > | * objects have been created. |
677 | > | * |
678 | > | */ |
679 | > | void SimInfo::update() { |
680 | > | setupSimVariables(); |
681 | calcNdf(); | |
682 | calcNdfRaw(); | |
683 | calcNdfTrans(); | |
709 | – | |
710 | – | fortranInitialized_ = true; |
684 | } | |
685 | < | |
686 | < | std::set<AtomType*> SimInfo::getUniqueAtomTypes() { |
685 | > | |
686 | > | /** |
687 | > | * getSimulatedAtomTypes |
688 | > | * |
689 | > | * Returns an STL set of AtomType* that are actually present in this |
690 | > | * simulation. Must query all processors to assemble this information. |
691 | > | * |
692 | > | */ |
693 | > | set<AtomType*> SimInfo::getSimulatedAtomTypes() { |
694 | SimInfo::MoleculeIterator mi; | |
695 | Molecule* mol; | |
696 | Molecule::AtomIterator ai; | |
697 | Atom* atom; | |
698 | < | std::set<AtomType*> atomTypes; |
699 | < | |
698 | > | set<AtomType*> atomTypes; |
699 | > | |
700 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
701 | < | |
702 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
701 | > | for(atom = mol->beginAtom(ai); atom != NULL; |
702 | > | atom = mol->nextAtom(ai)) { |
703 | atomTypes.insert(atom->getAtomType()); | |
704 | < | } |
705 | < | |
706 | < | } |
704 | > | } |
705 | > | } |
706 | > | |
707 | > | #ifdef IS_MPI |
708 | ||
709 | < | return atomTypes; |
710 | < | } |
730 | < | |
731 | < | void SimInfo::setupSimType() { |
732 | < | std::set<AtomType*>::iterator i; |
733 | < | std::set<AtomType*> atomTypes; |
734 | < | atomTypes = getUniqueAtomTypes(); |
709 | > | // loop over the found atom types on this processor, and add their |
710 | > | // numerical idents to a vector: |
711 | ||
712 | < | int useLennardJones = 0; |
713 | < | int useElectrostatic = 0; |
714 | < | int useEAM = 0; |
715 | < | int useSC = 0; |
740 | < | int useCharge = 0; |
741 | < | int useDirectional = 0; |
742 | < | int useDipole = 0; |
743 | < | int useGayBerne = 0; |
744 | < | int useSticky = 0; |
745 | < | int useStickyPower = 0; |
746 | < | int useShape = 0; |
747 | < | int useFLARB = 0; //it is not in AtomType yet |
748 | < | int useDirectionalAtom = 0; |
749 | < | int useElectrostatics = 0; |
750 | < | //usePBC and useRF are from simParams |
751 | < | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
752 | < | int useRF; |
753 | < | int useSF; |
754 | < | int useSP; |
755 | < | int useBoxDipole; |
712 | > | vector<int> foundTypes; |
713 | > | set<AtomType*>::iterator i; |
714 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) |
715 | > | foundTypes.push_back( (*i)->getIdent() ); |
716 | ||
717 | < | std::string myMethod; |
717 | > | // count_local holds the number of found types on this processor |
718 | > | int count_local = foundTypes.size(); |
719 | ||
720 | < | // set the useRF logical |
760 | < | useRF = 0; |
761 | < | useSF = 0; |
762 | < | useSP = 0; |
763 | < | useBoxDipole = 0; |
720 | > | int nproc = MPI::COMM_WORLD.Get_size(); |
721 | ||
722 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
723 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
724 | < | toUpper(myMethod); |
725 | < | if (myMethod == "REACTION_FIELD"){ |
726 | < | useRF = 1; |
727 | < | } else if (myMethod == "SHIFTED_FORCE"){ |
728 | < | useSF = 1; |
729 | < | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
730 | < | useSP = 1; |
731 | < | } |
722 | > | // we need arrays to hold the counts and displacement vectors for |
723 | > | // all processors |
724 | > | vector<int> counts(nproc, 0); |
725 | > | vector<int> disps(nproc, 0); |
726 | > | |
727 | > | // fill the counts array |
728 | > | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], |
729 | > | 1, MPI::INT); |
730 | > | |
731 | > | // use the processor counts to compute the displacement array |
732 | > | disps[0] = 0; |
733 | > | int totalCount = counts[0]; |
734 | > | for (int iproc = 1; iproc < nproc; iproc++) { |
735 | > | disps[iproc] = disps[iproc-1] + counts[iproc-1]; |
736 | > | totalCount += counts[iproc]; |
737 | } | |
738 | + | |
739 | + | // we need a (possibly redundant) set of all found types: |
740 | + | vector<int> ftGlobal(totalCount); |
741 | ||
742 | < | if (simParams_->haveAccumulateBoxDipole()) |
743 | < | if (simParams_->getAccumulateBoxDipole()) |
744 | < | useBoxDipole = 1; |
742 | > | // now spray out the foundTypes to all the other processors: |
743 | > | MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, |
744 | > | &ftGlobal[0], &counts[0], &disps[0], |
745 | > | MPI::INT); |
746 | ||
747 | < | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
747 | > | vector<int>::iterator j; |
748 | ||
749 | < | //loop over all of the atom types |
750 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
751 | < | useLennardJones |= (*i)->isLennardJones(); |
786 | < | useElectrostatic |= (*i)->isElectrostatic(); |
787 | < | useEAM |= (*i)->isEAM(); |
788 | < | useSC |= (*i)->isSC(); |
789 | < | useCharge |= (*i)->isCharge(); |
790 | < | useDirectional |= (*i)->isDirectional(); |
791 | < | useDipole |= (*i)->isDipole(); |
792 | < | useGayBerne |= (*i)->isGayBerne(); |
793 | < | useSticky |= (*i)->isSticky(); |
794 | < | useStickyPower |= (*i)->isStickyPower(); |
795 | < | useShape |= (*i)->isShape(); |
796 | < | } |
749 | > | // foundIdents is a stl set, so inserting an already found ident |
750 | > | // will have no effect. |
751 | > | set<int> foundIdents; |
752 | ||
753 | < | if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { |
754 | < | useDirectionalAtom = 1; |
755 | < | } |
753 | > | for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) |
754 | > | foundIdents.insert((*j)); |
755 | > | |
756 | > | // now iterate over the foundIdents and get the actual atom types |
757 | > | // that correspond to these: |
758 | > | set<int>::iterator it; |
759 | > | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
760 | > | atomTypes.insert( forceField_->getAtomType((*it)) ); |
761 | > | |
762 | > | #endif |
763 | ||
764 | < | if (useCharge || useDipole) { |
765 | < | useElectrostatics = 1; |
804 | < | } |
764 | > | return atomTypes; |
765 | > | } |
766 | ||
767 | + | void SimInfo::setupSimVariables() { |
768 | + | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
769 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
770 | + | calcBoxDipole_ = false; |
771 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
772 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
773 | + | calcBoxDipole_ = true; |
774 | + | } |
775 | + | |
776 | + | set<AtomType*>::iterator i; |
777 | + | set<AtomType*> atomTypes; |
778 | + | atomTypes = getSimulatedAtomTypes(); |
779 | + | int usesElectrostatic = 0; |
780 | + | int usesMetallic = 0; |
781 | + | int usesDirectional = 0; |
782 | + | //loop over all of the atom types |
783 | + | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
784 | + | usesElectrostatic |= (*i)->isElectrostatic(); |
785 | + | usesMetallic |= (*i)->isMetal(); |
786 | + | usesDirectional |= (*i)->isDirectional(); |
787 | + | } |
788 | + | |
789 | #ifdef IS_MPI | |
790 | int temp; | |
791 | + | temp = usesDirectional; |
792 | + | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
793 | + | |
794 | + | temp = usesMetallic; |
795 | + | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
796 | + | |
797 | + | temp = usesElectrostatic; |
798 | + | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
799 | + | #else |
800 | ||
801 | < | temp = usePBC; |
802 | < | MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
801 | > | usesDirectionalAtoms_ = usesDirectional; |
802 | > | usesMetallicAtoms_ = usesMetallic; |
803 | > | usesElectrostaticAtoms_ = usesElectrostatic; |
804 | ||
805 | < | temp = useDirectionalAtom; |
806 | < | MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
805 | > | #endif |
806 | > | |
807 | > | requiresPrepair_ = usesMetallicAtoms_ ? true : false; |
808 | > | requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; |
809 | > | requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; |
810 | > | } |
811 | ||
815 | – | temp = useLennardJones; |
816 | – | MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
812 | ||
813 | < | temp = useElectrostatics; |
814 | < | MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
813 | > | vector<int> SimInfo::getGlobalAtomIndices() { |
814 | > | SimInfo::MoleculeIterator mi; |
815 | > | Molecule* mol; |
816 | > | Molecule::AtomIterator ai; |
817 | > | Atom* atom; |
818 | ||
819 | < | temp = useCharge; |
822 | < | MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
823 | < | |
824 | < | temp = useDipole; |
825 | < | MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
826 | < | |
827 | < | temp = useSticky; |
828 | < | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
829 | < | |
830 | < | temp = useStickyPower; |
831 | < | MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
819 | > | vector<int> GlobalAtomIndices(getNAtoms(), 0); |
820 | ||
821 | < | temp = useGayBerne; |
822 | < | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
821 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
822 | > | |
823 | > | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
824 | > | GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); |
825 | > | } |
826 | > | } |
827 | > | return GlobalAtomIndices; |
828 | > | } |
829 | ||
836 | – | temp = useEAM; |
837 | – | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
830 | ||
831 | < | temp = useSC; |
832 | < | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
833 | < | |
834 | < | temp = useShape; |
835 | < | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
831 | > | vector<int> SimInfo::getGlobalGroupIndices() { |
832 | > | SimInfo::MoleculeIterator mi; |
833 | > | Molecule* mol; |
834 | > | Molecule::CutoffGroupIterator ci; |
835 | > | CutoffGroup* cg; |
836 | ||
837 | < | temp = useFLARB; |
838 | < | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
839 | < | |
840 | < | temp = useRF; |
841 | < | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
842 | < | |
843 | < | temp = useSF; |
844 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
845 | < | |
846 | < | temp = useSP; |
847 | < | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
848 | < | |
857 | < | temp = useBoxDipole; |
858 | < | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
859 | < | |
860 | < | temp = useAtomicVirial_; |
861 | < | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
862 | < | |
863 | < | #endif |
864 | < | fInfo_.SIM_uses_PBC = usePBC; |
865 | < | fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; |
866 | < | fInfo_.SIM_uses_LennardJones = useLennardJones; |
867 | < | fInfo_.SIM_uses_Electrostatics = useElectrostatics; |
868 | < | fInfo_.SIM_uses_Charges = useCharge; |
869 | < | fInfo_.SIM_uses_Dipoles = useDipole; |
870 | < | fInfo_.SIM_uses_Sticky = useSticky; |
871 | < | fInfo_.SIM_uses_StickyPower = useStickyPower; |
872 | < | fInfo_.SIM_uses_GayBerne = useGayBerne; |
873 | < | fInfo_.SIM_uses_EAM = useEAM; |
874 | < | fInfo_.SIM_uses_SC = useSC; |
875 | < | fInfo_.SIM_uses_Shapes = useShape; |
876 | < | fInfo_.SIM_uses_FLARB = useFLARB; |
877 | < | fInfo_.SIM_uses_RF = useRF; |
878 | < | fInfo_.SIM_uses_SF = useSF; |
879 | < | fInfo_.SIM_uses_SP = useSP; |
880 | < | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
881 | < | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
837 | > | vector<int> GlobalGroupIndices; |
838 | > | |
839 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
840 | > | |
841 | > | //local index of cutoff group is trivial, it only depends on the |
842 | > | //order of travesing |
843 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
844 | > | cg = mol->nextCutoffGroup(ci)) { |
845 | > | GlobalGroupIndices.push_back(cg->getGlobalIndex()); |
846 | > | } |
847 | > | } |
848 | > | return GlobalGroupIndices; |
849 | } | |
850 | ||
851 | < | void SimInfo::setupFortranSim() { |
852 | < | int isError; |
851 | > | |
852 | > | void SimInfo::prepareTopology() { |
853 | int nExclude, nOneTwo, nOneThree, nOneFour; | |
887 | – | std::vector<int> fortranGlobalGroupMembership; |
888 | – | |
889 | – | isError = 0; |
854 | ||
891 | – | //globalGroupMembership_ is filled by SimCreator |
892 | – | for (int i = 0; i < nGlobalAtoms_; i++) { |
893 | – | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
894 | – | } |
895 | – | |
855 | //calculate mass ratio of cutoff group | |
897 | – | std::vector<RealType> mfact; |
856 | SimInfo::MoleculeIterator mi; | |
857 | Molecule* mol; | |
858 | Molecule::CutoffGroupIterator ci; | |
# | Line 903 | Line 861 | namespace OpenMD { | |
861 | Atom* atom; | |
862 | RealType totalMass; | |
863 | ||
864 | < | //to avoid memory reallocation, reserve enough space for mfact |
865 | < | mfact.reserve(getNCutoffGroups()); |
864 | > | /** |
865 | > | * The mass factor is the relative mass of an atom to the total |
866 | > | * mass of the cutoff group it belongs to. By default, all atoms |
867 | > | * are their own cutoff groups, and therefore have mass factors of |
868 | > | * 1. We need some special handling for massless atoms, which |
869 | > | * will be treated as carrying the entire mass of the cutoff |
870 | > | * group. |
871 | > | */ |
872 | > | massFactors_.clear(); |
873 | > | massFactors_.resize(getNAtoms(), 1.0); |
874 | ||
875 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
876 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
876 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
877 | > | cg = mol->nextCutoffGroup(ci)) { |
878 | ||
879 | totalMass = cg->getMass(); | |
880 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
881 | // Check for massless groups - set mfact to 1 if true | |
882 | < | if (totalMass != 0) |
883 | < | mfact.push_back(atom->getMass()/totalMass); |
882 | > | if (totalMass != 0) |
883 | > | massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; |
884 | else | |
885 | < | mfact.push_back( 1.0 ); |
885 | > | massFactors_[atom->getLocalIndex()] = 1.0; |
886 | } | |
887 | } | |
888 | } | |
889 | ||
890 | < | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) |
924 | < | std::vector<int> identArray; |
890 | > | // Build the identArray_ |
891 | ||
892 | < | //to avoid memory reallocation, reserve enough space identArray |
893 | < | identArray.reserve(getNAtoms()); |
928 | < | |
892 | > | identArray_.clear(); |
893 | > | identArray_.reserve(getNAtoms()); |
894 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
895 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
896 | < | identArray.push_back(atom->getIdent()); |
896 | > | identArray_.push_back(atom->getIdent()); |
897 | } | |
898 | } | |
934 | – | |
935 | – | //fill molMembershipArray |
936 | – | //molMembershipArray is filled by SimCreator |
937 | – | std::vector<int> molMembershipArray(nGlobalAtoms_); |
938 | – | for (int i = 0; i < nGlobalAtoms_; i++) { |
939 | – | molMembershipArray[i] = globalMolMembership_[i] + 1; |
940 | – | } |
899 | ||
900 | < | //setup fortran simulation |
900 | > | //scan topology |
901 | ||
902 | nExclude = excludedInteractions_.getSize(); | |
903 | nOneTwo = oneTwoInteractions_.getSize(); | |
# | Line 951 | Line 909 | namespace OpenMD { | |
909 | int* oneThreeList = oneThreeInteractions_.getPairList(); | |
910 | int* oneFourList = oneFourInteractions_.getPairList(); | |
911 | ||
912 | < | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], |
955 | < | &nExclude, excludeList, |
956 | < | &nOneTwo, oneTwoList, |
957 | < | &nOneThree, oneThreeList, |
958 | < | &nOneFour, oneFourList, |
959 | < | &molMembershipArray[0], &mfact[0], &nCutoffGroups_, |
960 | < | &fortranGlobalGroupMembership[0], &isError); |
961 | < | |
962 | < | if( isError ){ |
963 | < | |
964 | < | sprintf( painCave.errMsg, |
965 | < | "There was an error setting the simulation information in fortran.\n" ); |
966 | < | painCave.isFatal = 1; |
967 | < | painCave.severity = OPENMD_ERROR; |
968 | < | simError(); |
969 | < | } |
970 | < | |
971 | < | |
972 | < | sprintf( checkPointMsg, |
973 | < | "succesfully sent the simulation information to fortran.\n"); |
974 | < | |
975 | < | errorCheckPoint(); |
976 | < | |
977 | < | // Setup number of neighbors in neighbor list if present |
978 | < | if (simParams_->haveNeighborListNeighbors()) { |
979 | < | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
980 | < | setNeighbors(&nlistNeighbors); |
981 | < | } |
982 | < | |
983 | < | |
984 | < | } |
985 | < | |
986 | < | |
987 | < | void SimInfo::setupFortranParallel() { |
988 | < | #ifdef IS_MPI |
989 | < | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex |
990 | < | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
991 | < | std::vector<int> localToGlobalCutoffGroupIndex; |
992 | < | SimInfo::MoleculeIterator mi; |
993 | < | Molecule::AtomIterator ai; |
994 | < | Molecule::CutoffGroupIterator ci; |
995 | < | Molecule* mol; |
996 | < | Atom* atom; |
997 | < | CutoffGroup* cg; |
998 | < | mpiSimData parallelData; |
999 | < | int isError; |
1000 | < | |
1001 | < | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
1002 | < | |
1003 | < | //local index(index in DataStorge) of atom is important |
1004 | < | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
1005 | < | localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; |
1006 | < | } |
1007 | < | |
1008 | < | //local index of cutoff group is trivial, it only depends on the order of travesing |
1009 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
1010 | < | localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); |
1011 | < | } |
1012 | < | |
1013 | < | } |
1014 | < | |
1015 | < | //fill up mpiSimData struct |
1016 | < | parallelData.nMolGlobal = getNGlobalMolecules(); |
1017 | < | parallelData.nMolLocal = getNMolecules(); |
1018 | < | parallelData.nAtomsGlobal = getNGlobalAtoms(); |
1019 | < | parallelData.nAtomsLocal = getNAtoms(); |
1020 | < | parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); |
1021 | < | parallelData.nGroupsLocal = getNCutoffGroups(); |
1022 | < | parallelData.myNode = worldRank; |
1023 | < | MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); |
1024 | < | |
1025 | < | //pass mpiSimData struct and index arrays to fortran |
1026 | < | setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), |
1027 | < | &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), |
1028 | < | &localToGlobalCutoffGroupIndex[0], &isError); |
1029 | < | |
1030 | < | if (isError) { |
1031 | < | sprintf(painCave.errMsg, |
1032 | < | "mpiRefresh errror: fortran didn't like something we gave it.\n"); |
1033 | < | painCave.isFatal = 1; |
1034 | < | simError(); |
1035 | < | } |
1036 | < | |
1037 | < | sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
1038 | < | errorCheckPoint(); |
1039 | < | |
1040 | < | #endif |
1041 | < | } |
1042 | < | |
1043 | < | void SimInfo::setupCutoff() { |
1044 | < | |
1045 | < | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
1046 | < | |
1047 | < | // Check the cutoff policy |
1048 | < | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default |
1049 | < | |
1050 | < | // Set LJ shifting bools to false |
1051 | < | ljsp_ = 0; |
1052 | < | ljsf_ = 0; |
1053 | < | |
1054 | < | std::string myPolicy; |
1055 | < | if (forceFieldOptions_.haveCutoffPolicy()){ |
1056 | < | myPolicy = forceFieldOptions_.getCutoffPolicy(); |
1057 | < | }else if (simParams_->haveCutoffPolicy()) { |
1058 | < | myPolicy = simParams_->getCutoffPolicy(); |
1059 | < | } |
1060 | < | |
1061 | < | if (!myPolicy.empty()){ |
1062 | < | toUpper(myPolicy); |
1063 | < | if (myPolicy == "MIX") { |
1064 | < | cp = MIX_CUTOFF_POLICY; |
1065 | < | } else { |
1066 | < | if (myPolicy == "MAX") { |
1067 | < | cp = MAX_CUTOFF_POLICY; |
1068 | < | } else { |
1069 | < | if (myPolicy == "TRADITIONAL") { |
1070 | < | cp = TRADITIONAL_CUTOFF_POLICY; |
1071 | < | } else { |
1072 | < | // throw error |
1073 | < | sprintf( painCave.errMsg, |
1074 | < | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
1075 | < | painCave.isFatal = 1; |
1076 | < | simError(); |
1077 | < | } |
1078 | < | } |
1079 | < | } |
1080 | < | } |
1081 | < | notifyFortranCutoffPolicy(&cp); |
1082 | < | |
1083 | < | // Check the Skin Thickness for neighborlists |
1084 | < | RealType skin; |
1085 | < | if (simParams_->haveSkinThickness()) { |
1086 | < | skin = simParams_->getSkinThickness(); |
1087 | < | notifyFortranSkinThickness(&skin); |
1088 | < | } |
1089 | < | |
1090 | < | // Check if the cutoff was set explicitly: |
1091 | < | if (simParams_->haveCutoffRadius()) { |
1092 | < | rcut_ = simParams_->getCutoffRadius(); |
1093 | < | if (simParams_->haveSwitchingRadius()) { |
1094 | < | rsw_ = simParams_->getSwitchingRadius(); |
1095 | < | } else { |
1096 | < | if (fInfo_.SIM_uses_Charges | |
1097 | < | fInfo_.SIM_uses_Dipoles | |
1098 | < | fInfo_.SIM_uses_RF) { |
1099 | < | |
1100 | < | rsw_ = 0.85 * rcut_; |
1101 | < | sprintf(painCave.errMsg, |
1102 | < | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1103 | < | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
1104 | < | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1105 | < | painCave.isFatal = 0; |
1106 | < | simError(); |
1107 | < | } else { |
1108 | < | rsw_ = rcut_; |
1109 | < | sprintf(painCave.errMsg, |
1110 | < | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1111 | < | "\tOpenMD will use the same value as the cutoffRadius.\n" |
1112 | < | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1113 | < | painCave.isFatal = 0; |
1114 | < | simError(); |
1115 | < | } |
1116 | < | } |
1117 | < | |
1118 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
1119 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1120 | < | toUpper(myMethod); |
1121 | < | |
1122 | < | if (myMethod == "SHIFTED_POTENTIAL") { |
1123 | < | ljsp_ = 1; |
1124 | < | } else if (myMethod == "SHIFTED_FORCE") { |
1125 | < | ljsf_ = 1; |
1126 | < | } |
1127 | < | } |
1128 | < | |
1129 | < | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1130 | < | |
1131 | < | } else { |
1132 | < | |
1133 | < | // For electrostatic atoms, we'll assume a large safe value: |
1134 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
1135 | < | sprintf(painCave.errMsg, |
1136 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
1137 | < | "\tOpenMD will use a default value of 15.0 angstroms" |
1138 | < | "\tfor the cutoffRadius.\n"); |
1139 | < | painCave.isFatal = 0; |
1140 | < | simError(); |
1141 | < | rcut_ = 15.0; |
1142 | < | |
1143 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
1144 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1145 | < | toUpper(myMethod); |
1146 | < | |
1147 | < | // For the time being, we're tethering the LJ shifted behavior to the |
1148 | < | // electrostaticSummationMethod keyword options |
1149 | < | if (myMethod == "SHIFTED_POTENTIAL") { |
1150 | < | ljsp_ = 1; |
1151 | < | } else if (myMethod == "SHIFTED_FORCE") { |
1152 | < | ljsf_ = 1; |
1153 | < | } |
1154 | < | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1155 | < | if (simParams_->haveSwitchingRadius()){ |
1156 | < | sprintf(painCave.errMsg, |
1157 | < | "SimInfo Warning: A value was set for the switchingRadius\n" |
1158 | < | "\teven though the electrostaticSummationMethod was\n" |
1159 | < | "\tset to %s\n", myMethod.c_str()); |
1160 | < | painCave.isFatal = 1; |
1161 | < | simError(); |
1162 | < | } |
1163 | < | } |
1164 | < | } |
1165 | < | |
1166 | < | if (simParams_->haveSwitchingRadius()){ |
1167 | < | rsw_ = simParams_->getSwitchingRadius(); |
1168 | < | } else { |
1169 | < | sprintf(painCave.errMsg, |
1170 | < | "SimCreator Warning: No value was set for switchingRadius.\n" |
1171 | < | "\tOpenMD will use a default value of\n" |
1172 | < | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
1173 | < | painCave.isFatal = 0; |
1174 | < | simError(); |
1175 | < | rsw_ = 0.85 * rcut_; |
1176 | < | } |
1177 | < | |
1178 | < | Electrostatic::setElectrostaticCutoffRadius(rcut_, rsw_); |
1179 | < | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1180 | < | |
1181 | < | } else { |
1182 | < | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so |
1183 | < | // We'll punt and let fortran figure out the cutoffs later. |
1184 | < | |
1185 | < | notifyFortranYouAreOnYourOwn(); |
1186 | < | |
1187 | < | } |
1188 | < | } |
1189 | < | } |
1190 | < | |
1191 | < | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
1192 | < | |
1193 | < | int errorOut; |
1194 | < | ElectrostaticSummationMethod esm = NONE; |
1195 | < | ElectrostaticScreeningMethod sm = UNDAMPED; |
1196 | < | RealType alphaVal; |
1197 | < | RealType dielectric; |
1198 | < | |
1199 | < | errorOut = isError; |
1200 | < | |
1201 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
1202 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1203 | < | toUpper(myMethod); |
1204 | < | if (myMethod == "NONE") { |
1205 | < | esm = NONE; |
1206 | < | } else { |
1207 | < | if (myMethod == "SWITCHING_FUNCTION") { |
1208 | < | esm = SWITCHING_FUNCTION; |
1209 | < | } else { |
1210 | < | if (myMethod == "SHIFTED_POTENTIAL") { |
1211 | < | esm = SHIFTED_POTENTIAL; |
1212 | < | } else { |
1213 | < | if (myMethod == "SHIFTED_FORCE") { |
1214 | < | esm = SHIFTED_FORCE; |
1215 | < | } else { |
1216 | < | if (myMethod == "REACTION_FIELD") { |
1217 | < | esm = REACTION_FIELD; |
1218 | < | dielectric = simParams_->getDielectric(); |
1219 | < | if (!simParams_->haveDielectric()) { |
1220 | < | // throw warning |
1221 | < | sprintf( painCave.errMsg, |
1222 | < | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1223 | < | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1224 | < | painCave.isFatal = 0; |
1225 | < | simError(); |
1226 | < | } |
1227 | < | } else { |
1228 | < | // throw error |
1229 | < | sprintf( painCave.errMsg, |
1230 | < | "SimInfo error: Unknown electrostaticSummationMethod.\n" |
1231 | < | "\t(Input file specified %s .)\n" |
1232 | < | "\telectrostaticSummationMethod must be one of: \"none\",\n" |
1233 | < | "\t\"shifted_potential\", \"shifted_force\", or \n" |
1234 | < | "\t\"reaction_field\".\n", myMethod.c_str() ); |
1235 | < | painCave.isFatal = 1; |
1236 | < | simError(); |
1237 | < | } |
1238 | < | } |
1239 | < | } |
1240 | < | } |
1241 | < | } |
1242 | < | } |
1243 | < | |
1244 | < | if (simParams_->haveElectrostaticScreeningMethod()) { |
1245 | < | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
1246 | < | toUpper(myScreen); |
1247 | < | if (myScreen == "UNDAMPED") { |
1248 | < | sm = UNDAMPED; |
1249 | < | } else { |
1250 | < | if (myScreen == "DAMPED") { |
1251 | < | sm = DAMPED; |
1252 | < | if (!simParams_->haveDampingAlpha()) { |
1253 | < | // first set a cutoff dependent alpha value |
1254 | < | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1255 | < | alphaVal = 0.5125 - rcut_* 0.025; |
1256 | < | // for values rcut > 20.5, alpha is zero |
1257 | < | if (alphaVal < 0) alphaVal = 0; |
1258 | < | |
1259 | < | // throw warning |
1260 | < | sprintf( painCave.errMsg, |
1261 | < | "SimInfo warning: dampingAlpha was not specified in the input file.\n" |
1262 | < | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1263 | < | painCave.isFatal = 0; |
1264 | < | simError(); |
1265 | < | } else { |
1266 | < | alphaVal = simParams_->getDampingAlpha(); |
1267 | < | } |
1268 | < | |
1269 | < | } else { |
1270 | < | // throw error |
1271 | < | sprintf( painCave.errMsg, |
1272 | < | "SimInfo error: Unknown electrostaticScreeningMethod.\n" |
1273 | < | "\t(Input file specified %s .)\n" |
1274 | < | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" |
1275 | < | "or \"damped\".\n", myScreen.c_str() ); |
1276 | < | painCave.isFatal = 1; |
1277 | < | simError(); |
1278 | < | } |
1279 | < | } |
1280 | < | } |
1281 | < | |
1282 | < | |
1283 | < | Electrostatic::setElectrostaticSummationMethod( esm ); |
1284 | < | Electrostatic::setElectrostaticScreeningMethod( sm ); |
1285 | < | Electrostatic::setDampingAlpha( alphaVal ); |
1286 | < | Electrostatic::setReactionFieldDielectric( dielectric ); |
1287 | < | initFortranFF( &errorOut ); |
912 | > | topologyDone_ = true; |
913 | } | |
914 | ||
1290 | – | void SimInfo::setupSwitchingFunction() { |
1291 | – | int ft = CUBIC; |
1292 | – | |
1293 | – | if (simParams_->haveSwitchingFunctionType()) { |
1294 | – | std::string funcType = simParams_->getSwitchingFunctionType(); |
1295 | – | toUpper(funcType); |
1296 | – | if (funcType == "CUBIC") { |
1297 | – | ft = CUBIC; |
1298 | – | } else { |
1299 | – | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1300 | – | ft = FIFTH_ORDER_POLY; |
1301 | – | } else { |
1302 | – | // throw error |
1303 | – | sprintf( painCave.errMsg, |
1304 | – | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
1305 | – | painCave.isFatal = 1; |
1306 | – | simError(); |
1307 | – | } |
1308 | – | } |
1309 | – | } |
1310 | – | |
1311 | – | // send switching function notification to switcheroo |
1312 | – | setFunctionType(&ft); |
1313 | – | |
1314 | – | } |
1315 | – | |
1316 | – | void SimInfo::setupAccumulateBoxDipole() { |
1317 | – | |
1318 | – | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1319 | – | if ( simParams_->haveAccumulateBoxDipole() ) |
1320 | – | if ( simParams_->getAccumulateBoxDipole() ) { |
1321 | – | setAccumulateBoxDipole(); |
1322 | – | calcBoxDipole_ = true; |
1323 | – | } |
1324 | – | |
1325 | – | } |
1326 | – | |
915 | void SimInfo::addProperty(GenericData* genData) { | |
916 | properties_.addProperty(genData); | |
917 | } | |
918 | ||
919 | < | void SimInfo::removeProperty(const std::string& propName) { |
919 | > | void SimInfo::removeProperty(const string& propName) { |
920 | properties_.removeProperty(propName); | |
921 | } | |
922 | ||
# | Line 1336 | Line 924 | namespace OpenMD { | |
924 | properties_.clearProperties(); | |
925 | } | |
926 | ||
927 | < | std::vector<std::string> SimInfo::getPropertyNames() { |
927 | > | vector<string> SimInfo::getPropertyNames() { |
928 | return properties_.getPropertyNames(); | |
929 | } | |
930 | ||
931 | < | std::vector<GenericData*> SimInfo::getProperties() { |
931 | > | vector<GenericData*> SimInfo::getProperties() { |
932 | return properties_.getProperties(); | |
933 | } | |
934 | ||
935 | < | GenericData* SimInfo::getPropertyByName(const std::string& propName) { |
935 | > | GenericData* SimInfo::getPropertyByName(const string& propName) { |
936 | return properties_.getPropertyByName(propName); | |
937 | } | |
938 | ||
# | Line 1358 | Line 946 | namespace OpenMD { | |
946 | Molecule* mol; | |
947 | RigidBody* rb; | |
948 | Atom* atom; | |
949 | + | CutoffGroup* cg; |
950 | SimInfo::MoleculeIterator mi; | |
951 | Molecule::RigidBodyIterator rbIter; | |
952 | < | Molecule::AtomIterator atomIter;; |
952 | > | Molecule::AtomIterator atomIter; |
953 | > | Molecule::CutoffGroupIterator cgIter; |
954 | ||
955 | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
956 | ||
# | Line 1371 | Line 961 | namespace OpenMD { | |
961 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
962 | rb->setSnapshotManager(sman_); | |
963 | } | |
964 | + | |
965 | + | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { |
966 | + | cg->setSnapshotManager(sman_); |
967 | + | } |
968 | } | |
969 | ||
970 | } | |
# | Line 1427 | Line 1021 | namespace OpenMD { | |
1021 | ||
1022 | } | |
1023 | ||
1024 | < | std::ostream& operator <<(std::ostream& o, SimInfo& info) { |
1024 | > | ostream& operator <<(ostream& o, SimInfo& info) { |
1025 | ||
1026 | return o; | |
1027 | } | |
# | Line 1577 | Line 1171 | namespace OpenMD { | |
1171 | return IOIndexToIntegrableObject.at(index); | |
1172 | } | |
1173 | ||
1174 | < | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1174 | > | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { |
1175 | IOIndexToIntegrableObject= v; | |
1176 | } | |
1177 | ||
# | Line 1619 | Line 1213 | namespace OpenMD { | |
1213 | return; | |
1214 | } | |
1215 | /* | |
1216 | < | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1216 | > | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { |
1217 | assert( v.size() == nAtoms_ + nRigidBodies_); | |
1218 | sdByGlobalIndex_ = v; | |
1219 | } | |
# | Line 1629 | Line 1223 | namespace OpenMD { | |
1223 | return sdByGlobalIndex_.at(index); | |
1224 | } | |
1225 | */ | |
1226 | + | int SimInfo::getNGlobalConstraints() { |
1227 | + | int nGlobalConstraints; |
1228 | + | #ifdef IS_MPI |
1229 | + | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
1230 | + | MPI_COMM_WORLD); |
1231 | + | #else |
1232 | + | nGlobalConstraints = nConstraints_; |
1233 | + | #endif |
1234 | + | return nGlobalConstraints; |
1235 | + | } |
1236 | + | |
1237 | }//end namespace OpenMD | |
1238 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |