# | Line 55 | Line 55 | |
---|---|---|
55 | #include "primitives/Molecule.hpp" | |
56 | #include "primitives/StuntDouble.hpp" | |
57 | #include "UseTheForce/fCutoffPolicy.h" | |
58 | – | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
59 | – | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" |
58 | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" | |
59 | #include "UseTheForce/doForces_interface.h" | |
60 | #include "UseTheForce/DarkSide/neighborLists_interface.h" | |
63 | – | #include "UseTheForce/DarkSide/electrostatic_interface.h" |
61 | #include "UseTheForce/DarkSide/switcheroo_interface.h" | |
62 | #include "utils/MemoryUtils.hpp" | |
63 | #include "utils/simError.h" | |
64 | #include "selection/SelectionManager.hpp" | |
65 | #include "io/ForceFieldOptions.hpp" | |
66 | #include "UseTheForce/ForceField.hpp" | |
67 | + | #include "nonbonded/InteractionManager.hpp" |
68 | ||
69 | ||
70 | #ifdef IS_MPI | |
# | Line 74 | Line 72 | |
72 | #include "UseTheForce/DarkSide/simParallel_interface.h" | |
73 | #endif | |
74 | ||
75 | + | using namespace std; |
76 | namespace OpenMD { | |
78 | – | std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) { |
79 | – | std::map<int, std::set<int> >::iterator i = container.find(index); |
80 | – | std::set<int> result; |
81 | – | if (i != container.end()) { |
82 | – | result = i->second; |
83 | – | } |
84 | – | |
85 | – | return result; |
86 | – | } |
77 | ||
78 | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | |
79 | forceField_(ff), simParams_(simParams), | |
# | Line 93 | Line 83 | namespace OpenMD { | |
83 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), | |
84 | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | |
85 | nConstraints_(0), sman_(NULL), fortranInitialized_(false), | |
86 | < | calcBoxDipole_(false), useAtomicVirial_(true) { |
87 | < | |
88 | < | |
89 | < | MoleculeStamp* molStamp; |
90 | < | int nMolWithSameStamp; |
91 | < | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
92 | < | int nGroups = 0; //total cutoff groups defined in meta-data file |
93 | < | CutoffGroupStamp* cgStamp; |
94 | < | RigidBodyStamp* rbStamp; |
95 | < | int nRigidAtoms = 0; |
96 | < | |
97 | < | std::vector<Component*> components = simParams->getComponents(); |
86 | > | calcBoxDipole_(false), useAtomicVirial_(true) { |
87 | > | |
88 | > | MoleculeStamp* molStamp; |
89 | > | int nMolWithSameStamp; |
90 | > | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
91 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
92 | > | CutoffGroupStamp* cgStamp; |
93 | > | RigidBodyStamp* rbStamp; |
94 | > | int nRigidAtoms = 0; |
95 | > | |
96 | > | vector<Component*> components = simParams->getComponents(); |
97 | > | |
98 | > | for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
99 | > | molStamp = (*i)->getMoleculeStamp(); |
100 | > | nMolWithSameStamp = (*i)->getNMol(); |
101 | ||
102 | < | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
103 | < | molStamp = (*i)->getMoleculeStamp(); |
104 | < | nMolWithSameStamp = (*i)->getNMol(); |
105 | < | |
106 | < | addMoleculeStamp(molStamp, nMolWithSameStamp); |
107 | < | |
108 | < | //calculate atoms in molecules |
109 | < | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
110 | < | |
111 | < | //calculate atoms in cutoff groups |
112 | < | int nAtomsInGroups = 0; |
113 | < | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); |
121 | < | |
122 | < | for (int j=0; j < nCutoffGroupsInStamp; j++) { |
123 | < | cgStamp = molStamp->getCutoffGroupStamp(j); |
124 | < | nAtomsInGroups += cgStamp->getNMembers(); |
125 | < | } |
126 | < | |
127 | < | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; |
128 | < | |
129 | < | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; |
130 | < | |
131 | < | //calculate atoms in rigid bodies |
132 | < | int nAtomsInRigidBodies = 0; |
133 | < | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); |
134 | < | |
135 | < | for (int j=0; j < nRigidBodiesInStamp; j++) { |
136 | < | rbStamp = molStamp->getRigidBodyStamp(j); |
137 | < | nAtomsInRigidBodies += rbStamp->getNMembers(); |
138 | < | } |
139 | < | |
140 | < | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; |
141 | < | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; |
142 | < | |
102 | > | addMoleculeStamp(molStamp, nMolWithSameStamp); |
103 | > | |
104 | > | //calculate atoms in molecules |
105 | > | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
106 | > | |
107 | > | //calculate atoms in cutoff groups |
108 | > | int nAtomsInGroups = 0; |
109 | > | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); |
110 | > | |
111 | > | for (int j=0; j < nCutoffGroupsInStamp; j++) { |
112 | > | cgStamp = molStamp->getCutoffGroupStamp(j); |
113 | > | nAtomsInGroups += cgStamp->getNMembers(); |
114 | } | |
115 | < | |
116 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff |
117 | < | //group therefore the total number of cutoff groups in the system is |
118 | < | //equal to the total number of atoms minus number of atoms belong to |
119 | < | //cutoff group defined in meta-data file plus the number of cutoff |
120 | < | //groups defined in meta-data file |
121 | < | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
122 | < | |
123 | < | //every free atom (atom does not belong to rigid bodies) is an |
124 | < | //integrable object therefore the total number of integrable objects |
125 | < | //in the system is equal to the total number of atoms minus number of |
126 | < | //atoms belong to rigid body defined in meta-data file plus the number |
127 | < | //of rigid bodies defined in meta-data file |
128 | < | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
129 | < | + nGlobalRigidBodies_; |
130 | < | |
131 | < | nGlobalMols_ = molStampIds_.size(); |
132 | < | molToProcMap_.resize(nGlobalMols_); |
133 | < | } |
134 | < | |
115 | > | |
116 | > | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; |
117 | > | |
118 | > | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; |
119 | > | |
120 | > | //calculate atoms in rigid bodies |
121 | > | int nAtomsInRigidBodies = 0; |
122 | > | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); |
123 | > | |
124 | > | for (int j=0; j < nRigidBodiesInStamp; j++) { |
125 | > | rbStamp = molStamp->getRigidBodyStamp(j); |
126 | > | nAtomsInRigidBodies += rbStamp->getNMembers(); |
127 | > | } |
128 | > | |
129 | > | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; |
130 | > | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; |
131 | > | |
132 | > | } |
133 | > | |
134 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff |
135 | > | //group therefore the total number of cutoff groups in the system is |
136 | > | //equal to the total number of atoms minus number of atoms belong to |
137 | > | //cutoff group defined in meta-data file plus the number of cutoff |
138 | > | //groups defined in meta-data file |
139 | > | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
140 | > | |
141 | > | //every free atom (atom does not belong to rigid bodies) is an |
142 | > | //integrable object therefore the total number of integrable objects |
143 | > | //in the system is equal to the total number of atoms minus number of |
144 | > | //atoms belong to rigid body defined in meta-data file plus the number |
145 | > | //of rigid bodies defined in meta-data file |
146 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
147 | > | + nGlobalRigidBodies_; |
148 | > | |
149 | > | nGlobalMols_ = molStampIds_.size(); |
150 | > | molToProcMap_.resize(nGlobalMols_); |
151 | > | } |
152 | > | |
153 | SimInfo::~SimInfo() { | |
154 | < | std::map<int, Molecule*>::iterator i; |
154 | > | map<int, Molecule*>::iterator i; |
155 | for (i = molecules_.begin(); i != molecules_.end(); ++i) { | |
156 | delete i->second; | |
157 | } | |
# | Line 173 | Line 162 | namespace OpenMD { | |
162 | delete forceField_; | |
163 | } | |
164 | ||
176 | – | int SimInfo::getNGlobalConstraints() { |
177 | – | int nGlobalConstraints; |
178 | – | #ifdef IS_MPI |
179 | – | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
180 | – | MPI_COMM_WORLD); |
181 | – | #else |
182 | – | nGlobalConstraints = nConstraints_; |
183 | – | #endif |
184 | – | return nGlobalConstraints; |
185 | – | } |
165 | ||
166 | bool SimInfo::addMolecule(Molecule* mol) { | |
167 | MoleculeIterator i; | |
168 | < | |
168 | > | |
169 | i = molecules_.find(mol->getGlobalIndex()); | |
170 | if (i == molecules_.end() ) { | |
171 | < | |
172 | < | molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); |
173 | < | |
171 | > | |
172 | > | molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); |
173 | > | |
174 | nAtoms_ += mol->getNAtoms(); | |
175 | nBonds_ += mol->getNBonds(); | |
176 | nBends_ += mol->getNBends(); | |
# | Line 201 | Line 180 | namespace OpenMD { | |
180 | nIntegrableObjects_ += mol->getNIntegrableObjects(); | |
181 | nCutoffGroups_ += mol->getNCutoffGroups(); | |
182 | nConstraints_ += mol->getNConstraintPairs(); | |
183 | < | |
183 | > | |
184 | addInteractionPairs(mol); | |
185 | < | |
185 | > | |
186 | return true; | |
187 | } else { | |
188 | return false; | |
189 | } | |
190 | } | |
191 | < | |
191 | > | |
192 | bool SimInfo::removeMolecule(Molecule* mol) { | |
193 | MoleculeIterator i; | |
194 | i = molecules_.find(mol->getGlobalIndex()); | |
# | Line 237 | Line 216 | namespace OpenMD { | |
216 | } else { | |
217 | return false; | |
218 | } | |
240 | – | |
241 | – | |
219 | } | |
220 | ||
221 | ||
# | Line 256 | Line 233 | namespace OpenMD { | |
233 | void SimInfo::calcNdf() { | |
234 | int ndf_local; | |
235 | MoleculeIterator i; | |
236 | < | std::vector<StuntDouble*>::iterator j; |
236 | > | vector<StuntDouble*>::iterator j; |
237 | Molecule* mol; | |
238 | StuntDouble* integrableObject; | |
239 | ||
# | Line 307 | Line 284 | namespace OpenMD { | |
284 | int ndfRaw_local; | |
285 | ||
286 | MoleculeIterator i; | |
287 | < | std::vector<StuntDouble*>::iterator j; |
287 | > | vector<StuntDouble*>::iterator j; |
288 | Molecule* mol; | |
289 | StuntDouble* integrableObject; | |
290 | ||
# | Line 356 | Line 333 | namespace OpenMD { | |
333 | ||
334 | void SimInfo::addInteractionPairs(Molecule* mol) { | |
335 | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | |
336 | < | std::vector<Bond*>::iterator bondIter; |
337 | < | std::vector<Bend*>::iterator bendIter; |
338 | < | std::vector<Torsion*>::iterator torsionIter; |
339 | < | std::vector<Inversion*>::iterator inversionIter; |
336 | > | vector<Bond*>::iterator bondIter; |
337 | > | vector<Bend*>::iterator bendIter; |
338 | > | vector<Torsion*>::iterator torsionIter; |
339 | > | vector<Inversion*>::iterator inversionIter; |
340 | Bond* bond; | |
341 | Bend* bend; | |
342 | Torsion* torsion; | |
# | Line 377 | Line 354 | namespace OpenMD { | |
354 | // always be excluded. These are done at the bottom of this | |
355 | // function. | |
356 | ||
357 | < | std::map<int, std::set<int> > atomGroups; |
357 | > | map<int, set<int> > atomGroups; |
358 | Molecule::RigidBodyIterator rbIter; | |
359 | RigidBody* rb; | |
360 | Molecule::IntegrableObjectIterator ii; | |
# | Line 389 | Line 366 | namespace OpenMD { | |
366 | ||
367 | if (integrableObject->isRigidBody()) { | |
368 | rb = static_cast<RigidBody*>(integrableObject); | |
369 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
370 | < | std::set<int> rigidAtoms; |
369 | > | vector<Atom*> atoms = rb->getAtoms(); |
370 | > | set<int> rigidAtoms; |
371 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
372 | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | |
373 | } | |
374 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
375 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
375 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
376 | } | |
377 | } else { | |
378 | < | std::set<int> oneAtomSet; |
378 | > | set<int> oneAtomSet; |
379 | oneAtomSet.insert(integrableObject->getGlobalIndex()); | |
380 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
380 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
381 | } | |
382 | } | |
383 | ||
# | Line 503 | Line 480 | namespace OpenMD { | |
480 | ||
481 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | |
482 | rb = mol->nextRigidBody(rbIter)) { | |
483 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
483 | > | vector<Atom*> atoms = rb->getAtoms(); |
484 | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | |
485 | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | |
486 | a = atoms[i]->getGlobalIndex(); | |
# | Line 517 | Line 494 | namespace OpenMD { | |
494 | ||
495 | void SimInfo::removeInteractionPairs(Molecule* mol) { | |
496 | ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); | |
497 | < | std::vector<Bond*>::iterator bondIter; |
498 | < | std::vector<Bend*>::iterator bendIter; |
499 | < | std::vector<Torsion*>::iterator torsionIter; |
500 | < | std::vector<Inversion*>::iterator inversionIter; |
497 | > | vector<Bond*>::iterator bondIter; |
498 | > | vector<Bend*>::iterator bendIter; |
499 | > | vector<Torsion*>::iterator torsionIter; |
500 | > | vector<Inversion*>::iterator inversionIter; |
501 | Bond* bond; | |
502 | Bend* bend; | |
503 | Torsion* torsion; | |
# | Line 530 | Line 507 | namespace OpenMD { | |
507 | int c; | |
508 | int d; | |
509 | ||
510 | < | std::map<int, std::set<int> > atomGroups; |
510 | > | map<int, set<int> > atomGroups; |
511 | Molecule::RigidBodyIterator rbIter; | |
512 | RigidBody* rb; | |
513 | Molecule::IntegrableObjectIterator ii; | |
# | Line 542 | Line 519 | namespace OpenMD { | |
519 | ||
520 | if (integrableObject->isRigidBody()) { | |
521 | rb = static_cast<RigidBody*>(integrableObject); | |
522 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
523 | < | std::set<int> rigidAtoms; |
522 | > | vector<Atom*> atoms = rb->getAtoms(); |
523 | > | set<int> rigidAtoms; |
524 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
525 | rigidAtoms.insert(atoms[i]->getGlobalIndex()); | |
526 | } | |
527 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
528 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
528 | > | atomGroups.insert(map<int, set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); |
529 | } | |
530 | } else { | |
531 | < | std::set<int> oneAtomSet; |
531 | > | set<int> oneAtomSet; |
532 | oneAtomSet.insert(integrableObject->getGlobalIndex()); | |
533 | < | atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
533 | > | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
534 | } | |
535 | } | |
536 | ||
# | Line 656 | Line 633 | namespace OpenMD { | |
633 | ||
634 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | |
635 | rb = mol->nextRigidBody(rbIter)) { | |
636 | < | std::vector<Atom*> atoms = rb->getAtoms(); |
636 | > | vector<Atom*> atoms = rb->getAtoms(); |
637 | for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) { | |
638 | for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) { | |
639 | a = atoms[i]->getGlobalIndex(); | |
# | Line 682 | Line 659 | namespace OpenMD { | |
659 | void SimInfo::update() { | |
660 | ||
661 | setupSimType(); | |
662 | + | setupCutoffRadius(); |
663 | + | setupSwitchingRadius(); |
664 | + | setupCutoffMethod(); |
665 | + | setupSkinThickness(); |
666 | + | setupSwitchingFunction(); |
667 | + | setupAccumulateBoxDipole(); |
668 | ||
669 | #ifdef IS_MPI | |
670 | setupFortranParallel(); | |
671 | #endif | |
689 | – | |
672 | setupFortranSim(); | |
673 | + | fortranInitialized_ = true; |
674 | ||
692 | – | //setup fortran force field |
693 | – | /** @deprecate */ |
694 | – | int isError = 0; |
695 | – | |
696 | – | setupCutoff(); |
697 | – | |
698 | – | setupElectrostaticSummationMethod( isError ); |
699 | – | setupSwitchingFunction(); |
700 | – | setupAccumulateBoxDipole(); |
701 | – | |
702 | – | if(isError){ |
703 | – | sprintf( painCave.errMsg, |
704 | – | "ForceField error: There was an error initializing the forceField in fortran.\n" ); |
705 | – | painCave.isFatal = 1; |
706 | – | simError(); |
707 | – | } |
708 | – | |
675 | calcNdf(); | |
676 | calcNdfRaw(); | |
677 | calcNdfTrans(); | |
712 | – | |
713 | – | fortranInitialized_ = true; |
678 | } | |
679 | < | |
680 | < | std::set<AtomType*> SimInfo::getUniqueAtomTypes() { |
679 | > | |
680 | > | set<AtomType*> SimInfo::getSimulatedAtomTypes() { |
681 | SimInfo::MoleculeIterator mi; | |
682 | Molecule* mol; | |
683 | Molecule::AtomIterator ai; | |
684 | Atom* atom; | |
685 | < | std::set<AtomType*> atomTypes; |
686 | < | |
685 | > | set<AtomType*> atomTypes; |
686 | > | |
687 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
688 | < | |
688 | > | |
689 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
690 | atomTypes.insert(atom->getAtomType()); | |
691 | } | |
692 | < | |
692 | > | |
693 | } | |
694 | < | |
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 | > | * setupCutoffRadius |
700 | > | * |
701 | > | * If the cutoffRadius was explicitly set, use that value. |
702 | > | * If the cutoffRadius was not explicitly set: |
703 | > | * Are there electrostatic atoms? Use 12.0 Angstroms. |
704 | > | * No electrostatic atoms? Poll the atom types present in the |
705 | > | * simulation for suggested cutoff values (e.g. 2.5 * sigma). |
706 | > | * Use the maximum suggested value that was found. |
707 | > | */ |
708 | > | void SimInfo::setupCutoffRadius() { |
709 | ||
710 | < | int useLennardJones = 0; |
711 | < | int useElectrostatic = 0; |
712 | < | int useEAM = 0; |
713 | < | int useSC = 0; |
714 | < | int useCharge = 0; |
715 | < | int useDirectional = 0; |
716 | < | int useDipole = 0; |
717 | < | int useGayBerne = 0; |
718 | < | int useSticky = 0; |
719 | < | int useStickyPower = 0; |
720 | < | int useShape = 0; |
721 | < | int useFLARB = 0; //it is not in AtomType yet |
722 | < | int useDirectionalAtom = 0; |
723 | < | int useElectrostatics = 0; |
724 | < | //usePBC and useRF are from simParams |
725 | < | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
726 | < | int useRF; |
727 | < | int useSF; |
728 | < | int useSP; |
729 | < | int useBoxDipole; |
730 | < | |
731 | < | std::string myMethod; |
732 | < | |
733 | < | // set the useRF logical |
734 | < | useRF = 0; |
735 | < | useSF = 0; |
736 | < | useSP = 0; |
737 | < | useBoxDipole = 0; |
710 | > | if (simParams_->haveCutoffRadius()) { |
711 | > | cutoffRadius_ = simParams_->getCutoffRadius(); |
712 | > | } else { |
713 | > | if (usesElectrostaticAtoms_) { |
714 | > | sprintf(painCave.errMsg, |
715 | > | "SimInfo Warning: No value was set for the cutoffRadius.\n" |
716 | > | "\tOpenMD will use a default value of 12.0 angstroms" |
717 | > | "\tfor the cutoffRadius.\n"); |
718 | > | painCave.isFatal = 0; |
719 | > | simError(); |
720 | > | cutoffRadius_ = 12.0; |
721 | > | } else { |
722 | > | RealType thisCut; |
723 | > | set<AtomType*>::iterator i; |
724 | > | set<AtomType*> atomTypes; |
725 | > | atomTypes = getSimulatedAtomTypes(); |
726 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
727 | > | thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); |
728 | > | cutoffRadius_ = max(thisCut, cutoffRadius_); |
729 | > | } |
730 | > | sprintf(painCave.errMsg, |
731 | > | "SimInfo Warning: No value was set for the cutoffRadius.\n" |
732 | > | "\tOpenMD will use %lf angstroms.\n", |
733 | > | cutoffRadius_); |
734 | > | painCave.isFatal = 0; |
735 | > | simError(); |
736 | > | } |
737 | > | } |
738 | ||
739 | + | InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); |
740 | + | } |
741 | + | |
742 | + | /** |
743 | + | * setupSwitchingRadius |
744 | + | * |
745 | + | * If the switchingRadius was explicitly set, use that value (but check it) |
746 | + | * If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ |
747 | + | */ |
748 | + | void SimInfo::setupSwitchingRadius() { |
749 | + | |
750 | + | if (simParams_->haveSwitchingRadius()) { |
751 | + | switchingRadius_ = simParams_->getSwitchingRadius(); |
752 | + | if (switchingRadius_ > cutoffRadius_) { |
753 | + | sprintf(painCave.errMsg, |
754 | + | "SimInfo Error: switchingRadius (%f) is larger than cutoffRadius(%f)\n", |
755 | + | switchingRadius_, cutoffRadius_); |
756 | + | painCave.isFatal = 1; |
757 | + | simError(); |
758 | ||
769 | – | if (simParams_->haveElectrostaticSummationMethod()) { |
770 | – | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
771 | – | toUpper(myMethod); |
772 | – | if (myMethod == "REACTION_FIELD"){ |
773 | – | useRF = 1; |
774 | – | } else if (myMethod == "SHIFTED_FORCE"){ |
775 | – | useSF = 1; |
776 | – | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
777 | – | useSP = 1; |
759 | } | |
760 | < | } |
761 | < | |
762 | < | if (simParams_->haveAccumulateBoxDipole()) |
763 | < | if (simParams_->getAccumulateBoxDipole()) |
764 | < | useBoxDipole = 1; |
760 | > | } else { |
761 | > | switchingRadius_ = 0.85 * cutoffRadius_; |
762 | > | sprintf(painCave.errMsg, |
763 | > | "SimInfo Warning: No value was set for the switchingRadius.\n" |
764 | > | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
765 | > | "\tswitchingRadius = %f. for this simulation\n", switchingRadius_); |
766 | > | painCave.isFatal = 0; |
767 | > | simError(); |
768 | > | } |
769 | > | InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); |
770 | > | } |
771 | ||
772 | + | /** |
773 | + | * setupSkinThickness |
774 | + | * |
775 | + | * If the skinThickness was explicitly set, use that value (but check it) |
776 | + | * If the skinThickness was not explicitly set: use 1.0 angstroms |
777 | + | */ |
778 | + | void SimInfo::setupSkinThickness() { |
779 | + | if (simParams_->haveSkinThickness()) { |
780 | + | skinThickness_ = simParams_->getSkinThickness(); |
781 | + | } else { |
782 | + | skinThickness_ = 1.0; |
783 | + | sprintf(painCave.errMsg, |
784 | + | "SimInfo Warning: No value was set for the skinThickness.\n" |
785 | + | "\tOpenMD will use a default value of %f Angstroms\n" |
786 | + | "\tfor this simulation\n", skinThickness_); |
787 | + | painCave.isFatal = 0; |
788 | + | simError(); |
789 | + | } |
790 | + | } |
791 | + | |
792 | + | void SimInfo::setupSimType() { |
793 | + | set<AtomType*>::iterator i; |
794 | + | set<AtomType*> atomTypes; |
795 | + | atomTypes = getSimulatedAtomTypes(); |
796 | + | |
797 | useAtomicVirial_ = simParams_->getUseAtomicVirial(); | |
798 | ||
799 | + | int usesElectrostatic = 0; |
800 | + | int usesMetallic = 0; |
801 | + | int usesDirectional = 0; |
802 | //loop over all of the atom types | |
803 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
804 | < | useLennardJones |= (*i)->isLennardJones(); |
805 | < | useElectrostatic |= (*i)->isElectrostatic(); |
806 | < | useEAM |= (*i)->isEAM(); |
792 | < | useSC |= (*i)->isSC(); |
793 | < | useCharge |= (*i)->isCharge(); |
794 | < | useDirectional |= (*i)->isDirectional(); |
795 | < | useDipole |= (*i)->isDipole(); |
796 | < | useGayBerne |= (*i)->isGayBerne(); |
797 | < | useSticky |= (*i)->isSticky(); |
798 | < | useStickyPower |= (*i)->isStickyPower(); |
799 | < | useShape |= (*i)->isShape(); |
804 | > | usesElectrostatic |= (*i)->isElectrostatic(); |
805 | > | usesMetallic |= (*i)->isMetal(); |
806 | > | usesDirectional |= (*i)->isDirectional(); |
807 | } | |
808 | ||
802 | – | if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { |
803 | – | useDirectionalAtom = 1; |
804 | – | } |
805 | – | |
806 | – | if (useCharge || useDipole) { |
807 | – | useElectrostatics = 1; |
808 | – | } |
809 | – | |
809 | #ifdef IS_MPI | |
810 | int temp; | |
811 | + | temp = usesDirectional; |
812 | + | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
813 | ||
814 | < | temp = usePBC; |
815 | < | MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
814 | > | temp = usesMetallic; |
815 | > | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
816 | ||
817 | < | temp = useDirectionalAtom; |
818 | < | MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
818 | < | |
819 | < | temp = useLennardJones; |
820 | < | MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
821 | < | |
822 | < | temp = useElectrostatics; |
823 | < | MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
824 | < | |
825 | < | temp = useCharge; |
826 | < | MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
827 | < | |
828 | < | temp = useDipole; |
829 | < | MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
830 | < | |
831 | < | temp = useSticky; |
832 | < | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
833 | < | |
834 | < | temp = useStickyPower; |
835 | < | MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
836 | < | |
837 | < | temp = useGayBerne; |
838 | < | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
839 | < | |
840 | < | temp = useEAM; |
841 | < | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
842 | < | |
843 | < | temp = useSC; |
844 | < | MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
845 | < | |
846 | < | temp = useShape; |
847 | < | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
848 | < | |
849 | < | temp = useFLARB; |
850 | < | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
851 | < | |
852 | < | temp = useRF; |
853 | < | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
854 | < | |
855 | < | temp = useSF; |
856 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
857 | < | |
858 | < | temp = useSP; |
859 | < | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
860 | < | |
861 | < | temp = useBoxDipole; |
862 | < | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
863 | < | |
864 | < | temp = useAtomicVirial_; |
865 | < | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
866 | < | |
817 | > | temp = usesElectrostatic; |
818 | > | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
819 | #endif | |
820 | < | |
821 | < | fInfo_.SIM_uses_PBC = usePBC; |
822 | < | fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; |
823 | < | fInfo_.SIM_uses_LennardJones = useLennardJones; |
824 | < | fInfo_.SIM_uses_Electrostatics = useElectrostatics; |
825 | < | fInfo_.SIM_uses_Charges = useCharge; |
874 | < | fInfo_.SIM_uses_Dipoles = useDipole; |
875 | < | fInfo_.SIM_uses_Sticky = useSticky; |
876 | < | fInfo_.SIM_uses_StickyPower = useStickyPower; |
877 | < | fInfo_.SIM_uses_GayBerne = useGayBerne; |
878 | < | fInfo_.SIM_uses_EAM = useEAM; |
879 | < | fInfo_.SIM_uses_SC = useSC; |
880 | < | fInfo_.SIM_uses_Shapes = useShape; |
881 | < | fInfo_.SIM_uses_FLARB = useFLARB; |
882 | < | fInfo_.SIM_uses_RF = useRF; |
883 | < | fInfo_.SIM_uses_SF = useSF; |
884 | < | fInfo_.SIM_uses_SP = useSP; |
885 | < | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
886 | < | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
820 | > | fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; |
821 | > | fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; |
822 | > | fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; |
823 | > | fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; |
824 | > | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; |
825 | > | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; |
826 | } | |
827 | ||
828 | void SimInfo::setupFortranSim() { | |
829 | int isError; | |
830 | int nExclude, nOneTwo, nOneThree, nOneFour; | |
831 | < | std::vector<int> fortranGlobalGroupMembership; |
831 | > | vector<int> fortranGlobalGroupMembership; |
832 | ||
833 | + | notifyFortranSkinThickness(&skinThickness_); |
834 | + | |
835 | + | int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
836 | + | int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
837 | + | notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); |
838 | + | |
839 | isError = 0; | |
840 | ||
841 | //globalGroupMembership_ is filled by SimCreator | |
# | Line 899 | Line 844 | namespace OpenMD { | |
844 | } | |
845 | ||
846 | //calculate mass ratio of cutoff group | |
847 | < | std::vector<RealType> mfact; |
847 | > | vector<RealType> mfact; |
848 | SimInfo::MoleculeIterator mi; | |
849 | Molecule* mol; | |
850 | Molecule::CutoffGroupIterator ci; | |
# | Line 926 | Line 871 | namespace OpenMD { | |
871 | } | |
872 | ||
873 | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) | |
874 | < | std::vector<int> identArray; |
874 | > | vector<int> identArray; |
875 | ||
876 | //to avoid memory reallocation, reserve enough space identArray | |
877 | identArray.reserve(getNAtoms()); | |
# | Line 939 | Line 884 | namespace OpenMD { | |
884 | ||
885 | //fill molMembershipArray | |
886 | //molMembershipArray is filled by SimCreator | |
887 | < | std::vector<int> molMembershipArray(nGlobalAtoms_); |
887 | > | vector<int> molMembershipArray(nGlobalAtoms_); |
888 | for (int i = 0; i < nGlobalAtoms_; i++) { | |
889 | molMembershipArray[i] = globalMolMembership_[i] + 1; | |
890 | } | |
# | Line 992 | Line 937 | namespace OpenMD { | |
937 | void SimInfo::setupFortranParallel() { | |
938 | #ifdef IS_MPI | |
939 | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex | |
940 | < | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
941 | < | std::vector<int> localToGlobalCutoffGroupIndex; |
940 | > | vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
941 | > | vector<int> localToGlobalCutoffGroupIndex; |
942 | SimInfo::MoleculeIterator mi; | |
943 | Molecule::AtomIterator ai; | |
944 | Molecule::CutoffGroupIterator ci; | |
# | Line 1045 | Line 990 | namespace OpenMD { | |
990 | #endif | |
991 | } | |
992 | ||
1048 | – | void SimInfo::setupCutoff() { |
1049 | – | |
1050 | – | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
993 | ||
994 | < | // Check the cutoff policy |
995 | < | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default |
1054 | < | |
1055 | < | // Set LJ shifting bools to false |
1056 | < | ljsp_ = 0; |
1057 | < | ljsf_ = 0; |
1058 | < | |
1059 | < | std::string myPolicy; |
1060 | < | if (forceFieldOptions_.haveCutoffPolicy()){ |
1061 | < | myPolicy = forceFieldOptions_.getCutoffPolicy(); |
1062 | < | }else if (simParams_->haveCutoffPolicy()) { |
1063 | < | myPolicy = simParams_->getCutoffPolicy(); |
1064 | < | } |
1065 | < | |
1066 | < | if (!myPolicy.empty()){ |
1067 | < | toUpper(myPolicy); |
1068 | < | if (myPolicy == "MIX") { |
1069 | < | cp = MIX_CUTOFF_POLICY; |
1070 | < | } else { |
1071 | < | if (myPolicy == "MAX") { |
1072 | < | cp = MAX_CUTOFF_POLICY; |
1073 | < | } else { |
1074 | < | if (myPolicy == "TRADITIONAL") { |
1075 | < | cp = TRADITIONAL_CUTOFF_POLICY; |
1076 | < | } else { |
1077 | < | // throw error |
1078 | < | sprintf( painCave.errMsg, |
1079 | < | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
1080 | < | painCave.isFatal = 1; |
1081 | < | simError(); |
1082 | < | } |
1083 | < | } |
1084 | < | } |
1085 | < | } |
1086 | < | notifyFortranCutoffPolicy(&cp); |
1087 | < | |
1088 | < | // Check the Skin Thickness for neighborlists |
1089 | < | RealType skin; |
1090 | < | if (simParams_->haveSkinThickness()) { |
1091 | < | skin = simParams_->getSkinThickness(); |
1092 | < | notifyFortranSkinThickness(&skin); |
1093 | < | } |
1094 | < | |
1095 | < | // Check if the cutoff was set explicitly: |
1096 | < | if (simParams_->haveCutoffRadius()) { |
1097 | < | rcut_ = simParams_->getCutoffRadius(); |
1098 | < | if (simParams_->haveSwitchingRadius()) { |
1099 | < | rsw_ = simParams_->getSwitchingRadius(); |
1100 | < | } else { |
1101 | < | if (fInfo_.SIM_uses_Charges | |
1102 | < | fInfo_.SIM_uses_Dipoles | |
1103 | < | fInfo_.SIM_uses_RF) { |
1104 | < | |
1105 | < | rsw_ = 0.85 * rcut_; |
1106 | < | sprintf(painCave.errMsg, |
1107 | < | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1108 | < | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
1109 | < | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1110 | < | painCave.isFatal = 0; |
1111 | < | simError(); |
1112 | < | } else { |
1113 | < | rsw_ = rcut_; |
1114 | < | sprintf(painCave.errMsg, |
1115 | < | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1116 | < | "\tOpenMD will use the same value as the cutoffRadius.\n" |
1117 | < | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1118 | < | painCave.isFatal = 0; |
1119 | < | simError(); |
1120 | < | } |
1121 | < | } |
1122 | < | |
1123 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
1124 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1125 | < | toUpper(myMethod); |
1126 | < | |
1127 | < | if (myMethod == "SHIFTED_POTENTIAL") { |
1128 | < | ljsp_ = 1; |
1129 | < | } else if (myMethod == "SHIFTED_FORCE") { |
1130 | < | ljsf_ = 1; |
1131 | < | } |
1132 | < | } |
1133 | < | |
1134 | < | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1135 | < | |
1136 | < | } else { |
1137 | < | |
1138 | < | // For electrostatic atoms, we'll assume a large safe value: |
1139 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
1140 | < | sprintf(painCave.errMsg, |
1141 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
1142 | < | "\tOpenMD will use a default value of 15.0 angstroms" |
1143 | < | "\tfor the cutoffRadius.\n"); |
1144 | < | painCave.isFatal = 0; |
1145 | < | simError(); |
1146 | < | rcut_ = 15.0; |
1147 | < | |
1148 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
1149 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1150 | < | toUpper(myMethod); |
1151 | < | |
1152 | < | // For the time being, we're tethering the LJ shifted behavior to the |
1153 | < | // electrostaticSummationMethod keyword options |
1154 | < | if (myMethod == "SHIFTED_POTENTIAL") { |
1155 | < | ljsp_ = 1; |
1156 | < | } else if (myMethod == "SHIFTED_FORCE") { |
1157 | < | ljsf_ = 1; |
1158 | < | } |
1159 | < | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1160 | < | if (simParams_->haveSwitchingRadius()){ |
1161 | < | sprintf(painCave.errMsg, |
1162 | < | "SimInfo Warning: A value was set for the switchingRadius\n" |
1163 | < | "\teven though the electrostaticSummationMethod was\n" |
1164 | < | "\tset to %s\n", myMethod.c_str()); |
1165 | < | painCave.isFatal = 1; |
1166 | < | simError(); |
1167 | < | } |
1168 | < | } |
1169 | < | } |
1170 | < | |
1171 | < | if (simParams_->haveSwitchingRadius()){ |
1172 | < | rsw_ = simParams_->getSwitchingRadius(); |
1173 | < | } else { |
1174 | < | sprintf(painCave.errMsg, |
1175 | < | "SimCreator Warning: No value was set for switchingRadius.\n" |
1176 | < | "\tOpenMD will use a default value of\n" |
1177 | < | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
1178 | < | painCave.isFatal = 0; |
1179 | < | simError(); |
1180 | < | rsw_ = 0.85 * rcut_; |
1181 | < | } |
1182 | < | |
1183 | < | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1184 | < | |
1185 | < | } else { |
1186 | < | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so |
1187 | < | // We'll punt and let fortran figure out the cutoffs later. |
1188 | < | |
1189 | < | notifyFortranYouAreOnYourOwn(); |
1190 | < | |
1191 | < | } |
1192 | < | } |
1193 | < | } |
1194 | < | |
1195 | < | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
1196 | < | |
1197 | < | int errorOut; |
1198 | < | int esm = NONE; |
1199 | < | int sm = UNDAMPED; |
1200 | < | RealType alphaVal; |
1201 | < | RealType dielectric; |
994 | > | void SimInfo::setupSwitchingFunction() { |
995 | > | int ft = CUBIC; |
996 | ||
997 | < | errorOut = isError; |
998 | < | |
1205 | < | if (simParams_->haveElectrostaticSummationMethod()) { |
1206 | < | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1207 | < | toUpper(myMethod); |
1208 | < | if (myMethod == "NONE") { |
1209 | < | esm = NONE; |
1210 | < | } else { |
1211 | < | if (myMethod == "SWITCHING_FUNCTION") { |
1212 | < | esm = SWITCHING_FUNCTION; |
1213 | < | } else { |
1214 | < | if (myMethod == "SHIFTED_POTENTIAL") { |
1215 | < | esm = SHIFTED_POTENTIAL; |
1216 | < | } else { |
1217 | < | if (myMethod == "SHIFTED_FORCE") { |
1218 | < | esm = SHIFTED_FORCE; |
1219 | < | } else { |
1220 | < | if (myMethod == "REACTION_FIELD") { |
1221 | < | esm = REACTION_FIELD; |
1222 | < | dielectric = simParams_->getDielectric(); |
1223 | < | if (!simParams_->haveDielectric()) { |
1224 | < | // throw warning |
1225 | < | sprintf( painCave.errMsg, |
1226 | < | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1227 | < | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1228 | < | painCave.isFatal = 0; |
1229 | < | simError(); |
1230 | < | } |
1231 | < | } else { |
1232 | < | // throw error |
1233 | < | sprintf( painCave.errMsg, |
1234 | < | "SimInfo error: Unknown electrostaticSummationMethod.\n" |
1235 | < | "\t(Input file specified %s .)\n" |
1236 | < | "\telectrostaticSummationMethod must be one of: \"none\",\n" |
1237 | < | "\t\"shifted_potential\", \"shifted_force\", or \n" |
1238 | < | "\t\"reaction_field\".\n", myMethod.c_str() ); |
1239 | < | painCave.isFatal = 1; |
1240 | < | simError(); |
1241 | < | } |
1242 | < | } |
1243 | < | } |
1244 | < | } |
1245 | < | } |
1246 | < | } |
1247 | < | |
1248 | < | if (simParams_->haveElectrostaticScreeningMethod()) { |
1249 | < | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
1250 | < | toUpper(myScreen); |
1251 | < | if (myScreen == "UNDAMPED") { |
1252 | < | sm = UNDAMPED; |
1253 | < | } else { |
1254 | < | if (myScreen == "DAMPED") { |
1255 | < | sm = DAMPED; |
1256 | < | if (!simParams_->haveDampingAlpha()) { |
1257 | < | // first set a cutoff dependent alpha value |
1258 | < | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1259 | < | alphaVal = 0.5125 - rcut_* 0.025; |
1260 | < | // for values rcut > 20.5, alpha is zero |
1261 | < | if (alphaVal < 0) alphaVal = 0; |
1262 | < | |
1263 | < | // throw warning |
1264 | < | sprintf( painCave.errMsg, |
1265 | < | "SimInfo warning: dampingAlpha was not specified in the input file.\n" |
1266 | < | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1267 | < | painCave.isFatal = 0; |
1268 | < | simError(); |
1269 | < | } else { |
1270 | < | alphaVal = simParams_->getDampingAlpha(); |
1271 | < | } |
1272 | < | |
1273 | < | } else { |
1274 | < | // throw error |
1275 | < | sprintf( painCave.errMsg, |
1276 | < | "SimInfo error: Unknown electrostaticScreeningMethod.\n" |
1277 | < | "\t(Input file specified %s .)\n" |
1278 | < | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" |
1279 | < | "or \"damped\".\n", myScreen.c_str() ); |
1280 | < | painCave.isFatal = 1; |
1281 | < | simError(); |
1282 | < | } |
1283 | < | } |
1284 | < | } |
1285 | < | |
1286 | < | // let's pass some summation method variables to fortran |
1287 | < | setElectrostaticSummationMethod( &esm ); |
1288 | < | setFortranElectrostaticMethod( &esm ); |
1289 | < | setScreeningMethod( &sm ); |
1290 | < | setDampingAlpha( &alphaVal ); |
1291 | < | setReactionFieldDielectric( &dielectric ); |
1292 | < | initFortranFF( &errorOut ); |
1293 | < | } |
1294 | < | |
1295 | < | void SimInfo::setupSwitchingFunction() { |
1296 | < | int ft = CUBIC; |
1297 | < | |
1298 | < | if (simParams_->haveSwitchingFunctionType()) { |
1299 | < | std::string funcType = simParams_->getSwitchingFunctionType(); |
997 | > | if (simParams_->haveSwitchingFunctionType()) { |
998 | > | string funcType = simParams_->getSwitchingFunctionType(); |
999 | toUpper(funcType); | |
1000 | if (funcType == "CUBIC") { | |
1001 | ft = CUBIC; | |
# | Line 1323 | Line 1022 | namespace OpenMD { | |
1022 | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | |
1023 | if ( simParams_->haveAccumulateBoxDipole() ) | |
1024 | if ( simParams_->getAccumulateBoxDipole() ) { | |
1326 | – | setAccumulateBoxDipole(); |
1025 | calcBoxDipole_ = true; | |
1026 | } | |
1027 | ||
# | Line 1333 | Line 1031 | namespace OpenMD { | |
1031 | properties_.addProperty(genData); | |
1032 | } | |
1033 | ||
1034 | < | void SimInfo::removeProperty(const std::string& propName) { |
1034 | > | void SimInfo::removeProperty(const string& propName) { |
1035 | properties_.removeProperty(propName); | |
1036 | } | |
1037 | ||
# | Line 1341 | Line 1039 | namespace OpenMD { | |
1039 | properties_.clearProperties(); | |
1040 | } | |
1041 | ||
1042 | < | std::vector<std::string> SimInfo::getPropertyNames() { |
1042 | > | vector<string> SimInfo::getPropertyNames() { |
1043 | return properties_.getPropertyNames(); | |
1044 | } | |
1045 | ||
1046 | < | std::vector<GenericData*> SimInfo::getProperties() { |
1046 | > | vector<GenericData*> SimInfo::getProperties() { |
1047 | return properties_.getProperties(); | |
1048 | } | |
1049 | ||
1050 | < | GenericData* SimInfo::getPropertyByName(const std::string& propName) { |
1050 | > | GenericData* SimInfo::getPropertyByName(const string& propName) { |
1051 | return properties_.getPropertyByName(propName); | |
1052 | } | |
1053 | ||
# | Line 1432 | Line 1130 | namespace OpenMD { | |
1130 | ||
1131 | } | |
1132 | ||
1133 | < | std::ostream& operator <<(std::ostream& o, SimInfo& info) { |
1133 | > | ostream& operator <<(ostream& o, SimInfo& info) { |
1134 | ||
1135 | return o; | |
1136 | } | |
# | Line 1475 | Line 1173 | namespace OpenMD { | |
1173 | ||
1174 | ||
1175 | [ Ixx -Ixy -Ixz ] | |
1176 | < | J =| -Iyx Iyy -Iyz | |
1176 | > | J =| -Iyx Iyy -Iyz | |
1177 | [ -Izx -Iyz Izz ] | |
1178 | */ | |
1179 | ||
# | Line 1582 | Line 1280 | namespace OpenMD { | |
1280 | return IOIndexToIntegrableObject.at(index); | |
1281 | } | |
1282 | ||
1283 | < | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1283 | > | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { |
1284 | IOIndexToIntegrableObject= v; | |
1285 | } | |
1286 | ||
# | Line 1624 | Line 1322 | namespace OpenMD { | |
1322 | return; | |
1323 | } | |
1324 | /* | |
1325 | < | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1325 | > | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { |
1326 | assert( v.size() == nAtoms_ + nRigidBodies_); | |
1327 | sdByGlobalIndex_ = v; | |
1328 | } | |
# | Line 1634 | Line 1332 | namespace OpenMD { | |
1332 | return sdByGlobalIndex_.at(index); | |
1333 | } | |
1334 | */ | |
1335 | + | int SimInfo::getNGlobalConstraints() { |
1336 | + | int nGlobalConstraints; |
1337 | + | #ifdef IS_MPI |
1338 | + | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
1339 | + | MPI_COMM_WORLD); |
1340 | + | #else |
1341 | + | nGlobalConstraints = nConstraints_; |
1342 | + | #endif |
1343 | + | return nGlobalConstraints; |
1344 | + | } |
1345 | + | |
1346 | }//end namespace OpenMD | |
1347 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |