# | Line 36 | Line 36 | |
---|---|---|
36 | * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | |
37 | * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | |
38 | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). | |
39 | < | * [4] Vardeman & Gezelter, in progress (2009). |
39 | > | * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 | > | * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 | */ | |
42 | ||
43 | /** | |
# | Line 54 | Line 55 | |
55 | #include "math/Vector3.hpp" | |
56 | #include "primitives/Molecule.hpp" | |
57 | #include "primitives/StuntDouble.hpp" | |
57 | – | #include "UseTheForce/doForces_interface.h" |
58 | – | #include "UseTheForce/DarkSide/neighborLists_interface.h" |
58 | #include "utils/MemoryUtils.hpp" | |
59 | #include "utils/simError.h" | |
60 | #include "selection/SelectionManager.hpp" | |
61 | #include "io/ForceFieldOptions.hpp" | |
62 | < | #include "UseTheForce/ForceField.hpp" |
62 | > | #include "brains/ForceField.hpp" |
63 | #include "nonbonded/SwitchingFunction.hpp" | |
65 | – | |
66 | – | |
64 | #ifdef IS_MPI | |
65 | < | #include "UseTheForce/mpiComponentPlan.h" |
66 | < | #include "UseTheForce/DarkSide/simParallel_interface.h" |
70 | < | #endif |
65 | > | #include <mpi.h> |
66 | > | #endif |
67 | ||
68 | using namespace std; | |
69 | namespace OpenMD { | |
# | Line 76 | Line 72 | namespace OpenMD { | |
72 | forceField_(ff), simParams_(simParams), | |
73 | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | |
74 | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | |
75 | < | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), |
75 | > | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0), |
76 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), | |
77 | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | |
78 | < | nConstraints_(0), sman_(NULL), fortranInitialized_(false), |
78 | > | nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), |
79 | calcBoxDipole_(false), useAtomicVirial_(true) { | |
80 | ||
81 | MoleculeStamp* molStamp; | |
# | Line 92 | Line 88 | namespace OpenMD { | |
88 | ||
89 | vector<Component*> components = simParams->getComponents(); | |
90 | ||
91 | < | for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
91 | > | for (vector<Component*>::iterator i = components.begin(); |
92 | > | i !=components.end(); ++i) { |
93 | molStamp = (*i)->getMoleculeStamp(); | |
94 | nMolWithSameStamp = (*i)->getNMol(); | |
95 | ||
# | Line 133 | Line 130 | namespace OpenMD { | |
130 | //equal to the total number of atoms minus number of atoms belong to | |
131 | //cutoff group defined in meta-data file plus the number of cutoff | |
132 | //groups defined in meta-data file | |
133 | + | |
134 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
135 | ||
136 | //every free atom (atom does not belong to rigid bodies) is an | |
# | Line 228 | Line 226 | namespace OpenMD { | |
226 | ||
227 | ||
228 | void SimInfo::calcNdf() { | |
229 | < | int ndf_local; |
229 | > | int ndf_local, nfq_local; |
230 | MoleculeIterator i; | |
231 | vector<StuntDouble*>::iterator j; | |
232 | + | vector<Atom*>::iterator k; |
233 | + | |
234 | Molecule* mol; | |
235 | < | StuntDouble* integrableObject; |
235 | > | StuntDouble* sd; |
236 | > | Atom* atom; |
237 | ||
238 | ndf_local = 0; | |
239 | + | nfq_local = 0; |
240 | ||
241 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
240 | – | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
241 | – | integrableObject = mol->nextIntegrableObject(j)) { |
242 | ||
243 | + | for (sd = mol->beginIntegrableObject(j); sd != NULL; |
244 | + | sd = mol->nextIntegrableObject(j)) { |
245 | + | |
246 | ndf_local += 3; | |
247 | ||
248 | < | if (integrableObject->isDirectional()) { |
249 | < | if (integrableObject->isLinear()) { |
248 | > | if (sd->isDirectional()) { |
249 | > | if (sd->isLinear()) { |
250 | ndf_local += 2; | |
251 | } else { | |
252 | ndf_local += 3; | |
253 | } | |
254 | } | |
252 | – | |
255 | } | |
256 | + | |
257 | + | for (atom = mol->beginFluctuatingCharge(k); atom != NULL; |
258 | + | atom = mol->nextFluctuatingCharge(k)) { |
259 | + | if (atom->isFluctuatingCharge()) { |
260 | + | nfq_local++; |
261 | + | } |
262 | + | } |
263 | } | |
264 | ||
265 | + | ndfLocal_ = ndf_local; |
266 | + | |
267 | // n_constraints is local, so subtract them on each processor | |
268 | ndf_local -= nConstraints_; | |
269 | ||
270 | #ifdef IS_MPI | |
271 | < | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
271 | > | MPI::COMM_WORLD.Allreduce(&ndf_local, &ndf_, 1, MPI::INT,MPI::SUM); |
272 | > | MPI::COMM_WORLD.Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1, |
273 | > | MPI::INT, MPI::SUM); |
274 | #else | |
275 | ndf_ = ndf_local; | |
276 | + | nGlobalFluctuatingCharges_ = nfq_local; |
277 | #endif | |
278 | ||
279 | // nZconstraints_ is global, as are the 3 COM translations for the | |
# | Line 270 | Line 284 | namespace OpenMD { | |
284 | ||
285 | int SimInfo::getFdf() { | |
286 | #ifdef IS_MPI | |
287 | < | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
287 | > | MPI::COMM_WORLD.Allreduce(&fdf_local, &fdf_, 1, MPI::INT, MPI::SUM); |
288 | #else | |
289 | fdf_ = fdf_local; | |
290 | #endif | |
291 | return fdf_; | |
292 | + | } |
293 | + | |
294 | + | unsigned int SimInfo::getNLocalCutoffGroups(){ |
295 | + | int nLocalCutoffAtoms = 0; |
296 | + | Molecule* mol; |
297 | + | MoleculeIterator mi; |
298 | + | CutoffGroup* cg; |
299 | + | Molecule::CutoffGroupIterator ci; |
300 | + | |
301 | + | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
302 | + | |
303 | + | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
304 | + | cg = mol->nextCutoffGroup(ci)) { |
305 | + | nLocalCutoffAtoms += cg->getNumAtom(); |
306 | + | |
307 | + | } |
308 | + | } |
309 | + | |
310 | + | return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; |
311 | } | |
312 | ||
313 | void SimInfo::calcNdfRaw() { | |
# | Line 283 | Line 316 | namespace OpenMD { | |
316 | MoleculeIterator i; | |
317 | vector<StuntDouble*>::iterator j; | |
318 | Molecule* mol; | |
319 | < | StuntDouble* integrableObject; |
319 | > | StuntDouble* sd; |
320 | ||
321 | // Raw degrees of freedom that we have to set | |
322 | ndfRaw_local = 0; | |
323 | ||
324 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
292 | – | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
293 | – | integrableObject = mol->nextIntegrableObject(j)) { |
325 | ||
326 | + | for (sd = mol->beginIntegrableObject(j); sd != NULL; |
327 | + | sd = mol->nextIntegrableObject(j)) { |
328 | + | |
329 | ndfRaw_local += 3; | |
330 | ||
331 | < | if (integrableObject->isDirectional()) { |
332 | < | if (integrableObject->isLinear()) { |
331 | > | if (sd->isDirectional()) { |
332 | > | if (sd->isLinear()) { |
333 | ndfRaw_local += 2; | |
334 | } else { | |
335 | ndfRaw_local += 3; | |
# | Line 306 | Line 340 | namespace OpenMD { | |
340 | } | |
341 | ||
342 | #ifdef IS_MPI | |
343 | < | MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
343 | > | MPI::COMM_WORLD.Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI::INT, MPI::SUM); |
344 | #else | |
345 | ndfRaw_ = ndfRaw_local; | |
346 | #endif | |
# | Line 319 | Line 353 | namespace OpenMD { | |
353 | ||
354 | ||
355 | #ifdef IS_MPI | |
356 | < | MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
356 | > | MPI::COMM_WORLD.Allreduce(&ndfTrans_local, &ndfTrans_, 1, |
357 | > | MPI::INT, MPI::SUM); |
358 | #else | |
359 | ndfTrans_ = ndfTrans_local; | |
360 | #endif | |
# | Line 355 | Line 390 | namespace OpenMD { | |
390 | Molecule::RigidBodyIterator rbIter; | |
391 | RigidBody* rb; | |
392 | Molecule::IntegrableObjectIterator ii; | |
393 | < | StuntDouble* integrableObject; |
393 | > | StuntDouble* sd; |
394 | ||
395 | < | for (integrableObject = mol->beginIntegrableObject(ii); |
396 | < | integrableObject != NULL; |
362 | < | integrableObject = mol->nextIntegrableObject(ii)) { |
395 | > | for (sd = mol->beginIntegrableObject(ii); sd != NULL; |
396 | > | sd = mol->nextIntegrableObject(ii)) { |
397 | ||
398 | < | if (integrableObject->isRigidBody()) { |
399 | < | rb = static_cast<RigidBody*>(integrableObject); |
398 | > | if (sd->isRigidBody()) { |
399 | > | rb = static_cast<RigidBody*>(sd); |
400 | vector<Atom*> atoms = rb->getAtoms(); | |
401 | set<int> rigidAtoms; | |
402 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
# | Line 373 | Line 407 | namespace OpenMD { | |
407 | } | |
408 | } else { | |
409 | set<int> oneAtomSet; | |
410 | < | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
411 | < | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
410 | > | oneAtomSet.insert(sd->getGlobalIndex()); |
411 | > | atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet)); |
412 | } | |
413 | } | |
414 | ||
# | Line 508 | Line 542 | namespace OpenMD { | |
542 | Molecule::RigidBodyIterator rbIter; | |
543 | RigidBody* rb; | |
544 | Molecule::IntegrableObjectIterator ii; | |
545 | < | StuntDouble* integrableObject; |
546 | < | |
547 | < | for (integrableObject = mol->beginIntegrableObject(ii); |
548 | < | integrableObject != NULL; |
515 | < | integrableObject = mol->nextIntegrableObject(ii)) { |
545 | > | StuntDouble* sd; |
546 | > | |
547 | > | for (sd = mol->beginIntegrableObject(ii); sd != NULL; |
548 | > | sd = mol->nextIntegrableObject(ii)) { |
549 | ||
550 | < | if (integrableObject->isRigidBody()) { |
551 | < | rb = static_cast<RigidBody*>(integrableObject); |
550 | > | if (sd->isRigidBody()) { |
551 | > | rb = static_cast<RigidBody*>(sd); |
552 | vector<Atom*> atoms = rb->getAtoms(); | |
553 | set<int> rigidAtoms; | |
554 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
# | Line 526 | Line 559 | namespace OpenMD { | |
559 | } | |
560 | } else { | |
561 | set<int> oneAtomSet; | |
562 | < | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
563 | < | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
562 | > | oneAtomSet.insert(sd->getGlobalIndex()); |
563 | > | atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet)); |
564 | } | |
565 | } | |
566 | ||
# | Line 657 | Line 690 | namespace OpenMD { | |
690 | /** | |
691 | * update | |
692 | * | |
693 | < | * Performs the global checks and variable settings after the objects have been |
694 | < | * created. |
693 | > | * Performs the global checks and variable settings after the |
694 | > | * objects have been created. |
695 | * | |
696 | */ | |
697 | < | void SimInfo::update() { |
665 | < | |
697 | > | void SimInfo::update() { |
698 | setupSimVariables(); | |
667 | – | setupCutoffs(); |
668 | – | setupSwitching(); |
669 | – | setupElectrostatics(); |
670 | – | setupNeighborlists(); |
671 | – | |
672 | – | #ifdef IS_MPI |
673 | – | setupFortranParallel(); |
674 | – | #endif |
675 | – | setupFortranSim(); |
676 | – | fortranInitialized_ = true; |
677 | – | |
699 | calcNdf(); | |
700 | calcNdfRaw(); | |
701 | calcNdfTrans(); | |
702 | } | |
703 | ||
704 | + | /** |
705 | + | * getSimulatedAtomTypes |
706 | + | * |
707 | + | * Returns an STL set of AtomType* that are actually present in this |
708 | + | * simulation. Must query all processors to assemble this information. |
709 | + | * |
710 | + | */ |
711 | set<AtomType*> SimInfo::getSimulatedAtomTypes() { | |
712 | SimInfo::MoleculeIterator mi; | |
713 | Molecule* mol; | |
# | Line 687 | Line 715 | namespace OpenMD { | |
715 | Atom* atom; | |
716 | set<AtomType*> atomTypes; | |
717 | ||
718 | < | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
719 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
718 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
719 | > | for(atom = mol->beginAtom(ai); atom != NULL; |
720 | > | atom = mol->nextAtom(ai)) { |
721 | atomTypes.insert(atom->getAtomType()); | |
722 | } | |
723 | } | |
724 | < | return atomTypes; |
725 | < | } |
724 | > | |
725 | > | #ifdef IS_MPI |
726 | ||
727 | < | /** |
728 | < | * setupCutoffs |
700 | < | * |
701 | < | * Sets the values of cutoffRadius and cutoffMethod |
702 | < | * |
703 | < | * cutoffRadius : realType |
704 | < | * If the cutoffRadius was explicitly set, use that value. |
705 | < | * If the cutoffRadius was not explicitly set: |
706 | < | * Are there electrostatic atoms? Use 12.0 Angstroms. |
707 | < | * No electrostatic atoms? Poll the atom types present in the |
708 | < | * simulation for suggested cutoff values (e.g. 2.5 * sigma). |
709 | < | * Use the maximum suggested value that was found. |
710 | < | * |
711 | < | * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) |
712 | < | * If cutoffMethod was explicitly set, use that choice. |
713 | < | * If cutoffMethod was not explicitly set, use SHIFTED_FORCE |
714 | < | */ |
715 | < | void SimInfo::setupCutoffs() { |
727 | > | // loop over the found atom types on this processor, and add their |
728 | > | // numerical idents to a vector: |
729 | ||
730 | < | if (simParams_->haveCutoffRadius()) { |
731 | < | cutoffRadius_ = simParams_->getCutoffRadius(); |
732 | < | } else { |
733 | < | if (usesElectrostaticAtoms_) { |
721 | < | sprintf(painCave.errMsg, |
722 | < | "SimInfo: No value was set for the cutoffRadius.\n" |
723 | < | "\tOpenMD will use a default value of 12.0 angstroms" |
724 | < | "\tfor the cutoffRadius.\n"); |
725 | < | painCave.isFatal = 0; |
726 | < | painCave.severity = OPENMD_INFO; |
727 | < | simError(); |
728 | < | cutoffRadius_ = 12.0; |
729 | < | } else { |
730 | < | RealType thisCut; |
731 | < | set<AtomType*>::iterator i; |
732 | < | set<AtomType*> atomTypes; |
733 | < | atomTypes = getSimulatedAtomTypes(); |
734 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
735 | < | thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); |
736 | < | cutoffRadius_ = max(thisCut, cutoffRadius_); |
737 | < | } |
738 | < | sprintf(painCave.errMsg, |
739 | < | "SimInfo: No value was set for the cutoffRadius.\n" |
740 | < | "\tOpenMD will use %lf angstroms.\n", |
741 | < | cutoffRadius_); |
742 | < | painCave.isFatal = 0; |
743 | < | painCave.severity = OPENMD_INFO; |
744 | < | simError(); |
745 | < | } |
746 | < | } |
730 | > | vector<int> foundTypes; |
731 | > | set<AtomType*>::iterator i; |
732 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) |
733 | > | foundTypes.push_back( (*i)->getIdent() ); |
734 | ||
735 | < | InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); |
735 | > | // count_local holds the number of found types on this processor |
736 | > | int count_local = foundTypes.size(); |
737 | ||
738 | < | map<string, CutoffMethod> stringToCutoffMethod; |
739 | < | stringToCutoffMethod["HARD"] = HARD; |
740 | < | stringToCutoffMethod["SWITCHING_FUNCTION"] = SWITCHING_FUNCTION; |
741 | < | stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; |
742 | < | stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; |
738 | > | int nproc = MPI::COMM_WORLD.Get_size(); |
739 | > | |
740 | > | // we need arrays to hold the counts and displacement vectors for |
741 | > | // all processors |
742 | > | vector<int> counts(nproc, 0); |
743 | > | vector<int> disps(nproc, 0); |
744 | > | |
745 | > | // fill the counts array |
746 | > | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], |
747 | > | 1, MPI::INT); |
748 | ||
749 | < | if (simParams_->haveCutoffMethod()) { |
750 | < | string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); |
751 | < | map<string, CutoffMethod>::iterator i; |
752 | < | i = stringToCutoffMethod.find(cutMeth); |
753 | < | if (i == stringToCutoffMethod.end()) { |
754 | < | sprintf(painCave.errMsg, |
762 | < | "SimInfo: Could not find chosen cutoffMethod %s\n" |
763 | < | "\tShould be one of: " |
764 | < | "HARD, SWITCHING_FUNCTION, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", |
765 | < | cutMeth.c_str()); |
766 | < | painCave.isFatal = 1; |
767 | < | painCave.severity = OPENMD_ERROR; |
768 | < | simError(); |
769 | < | } else { |
770 | < | cutoffMethod_ = i->second; |
771 | < | } |
772 | < | } else { |
773 | < | sprintf(painCave.errMsg, |
774 | < | "SimInfo: No value was set for the cutoffMethod.\n" |
775 | < | "\tOpenMD will use SHIFTED_FORCE.\n"); |
776 | < | painCave.isFatal = 0; |
777 | < | painCave.severity = OPENMD_INFO; |
778 | < | simError(); |
779 | < | cutoffMethod_ = SHIFTED_FORCE; |
749 | > | // use the processor counts to compute the displacement array |
750 | > | disps[0] = 0; |
751 | > | int totalCount = counts[0]; |
752 | > | for (int iproc = 1; iproc < nproc; iproc++) { |
753 | > | disps[iproc] = disps[iproc-1] + counts[iproc-1]; |
754 | > | totalCount += counts[iproc]; |
755 | } | |
756 | ||
757 | < | InteractionManager::Instance()->setCutoffMethod(cutoffMethod_); |
758 | < | } |
784 | < | |
785 | < | /** |
786 | < | * setupSwitching |
787 | < | * |
788 | < | * Sets the values of switchingRadius and |
789 | < | * If the switchingRadius was explicitly set, use that value (but check it) |
790 | < | * If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ |
791 | < | */ |
792 | < | void SimInfo::setupSwitching() { |
757 | > | // we need a (possibly redundant) set of all found types: |
758 | > | vector<int> ftGlobal(totalCount); |
759 | ||
760 | < | if (simParams_->haveSwitchingRadius()) { |
761 | < | switchingRadius_ = simParams_->getSwitchingRadius(); |
762 | < | if (switchingRadius_ > cutoffRadius_) { |
763 | < | sprintf(painCave.errMsg, |
798 | < | "SimInfo: switchingRadius (%f) is larger than cutoffRadius(%f)\n", |
799 | < | switchingRadius_, cutoffRadius_); |
800 | < | painCave.isFatal = 1; |
801 | < | painCave.severity = OPENMD_ERROR; |
802 | < | simError(); |
803 | < | } |
804 | < | } else { |
805 | < | switchingRadius_ = 0.85 * cutoffRadius_; |
806 | < | sprintf(painCave.errMsg, |
807 | < | "SimInfo: No value was set for the switchingRadius.\n" |
808 | < | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
809 | < | "\tswitchingRadius = %f. for this simulation\n", switchingRadius_); |
810 | < | painCave.isFatal = 0; |
811 | < | painCave.severity = OPENMD_WARNING; |
812 | < | simError(); |
813 | < | } |
814 | < | |
815 | < | InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); |
760 | > | // now spray out the foundTypes to all the other processors: |
761 | > | MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, |
762 | > | &ftGlobal[0], &counts[0], &disps[0], |
763 | > | MPI::INT); |
764 | ||
765 | < | SwitchingFunctionType ft; |
765 | > | vector<int>::iterator j; |
766 | > | |
767 | > | // foundIdents is a stl set, so inserting an already found ident |
768 | > | // will have no effect. |
769 | > | set<int> foundIdents; |
770 | > | |
771 | > | for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) |
772 | > | foundIdents.insert((*j)); |
773 | ||
774 | < | if (simParams_->haveSwitchingFunctionType()) { |
775 | < | string funcType = simParams_->getSwitchingFunctionType(); |
776 | < | toUpper(funcType); |
777 | < | if (funcType == "CUBIC") { |
778 | < | ft = cubic; |
779 | < | } else { |
780 | < | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
826 | < | ft = fifth_order_poly; |
827 | < | } else { |
828 | < | // throw error |
829 | < | sprintf( painCave.errMsg, |
830 | < | "SimInfo : Unknown switchingFunctionType. (Input file specified %s .)\n" |
831 | < | "\tswitchingFunctionType must be one of: " |
832 | < | "\"cubic\" or \"fifth_order_polynomial\".", |
833 | < | funcType.c_str() ); |
834 | < | painCave.isFatal = 1; |
835 | < | painCave.severity = OPENMD_ERROR; |
836 | < | simError(); |
837 | < | } |
838 | < | } |
839 | < | } |
774 | > | // now iterate over the foundIdents and get the actual atom types |
775 | > | // that correspond to these: |
776 | > | set<int>::iterator it; |
777 | > | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
778 | > | atomTypes.insert( forceField_->getAtomType((*it)) ); |
779 | > | |
780 | > | #endif |
781 | ||
782 | < | InteractionManager::Instance()->setSwitchingFunctionType(ft); |
782 | > | return atomTypes; |
783 | } | |
784 | ||
785 | < | /** |
786 | < | * setupSkinThickness |
787 | < | * |
788 | < | * If the skinThickness was explicitly set, use that value (but check it) |
789 | < | * If the skinThickness was not explicitly set: use 1.0 angstroms |
790 | < | */ |
791 | < | void SimInfo::setupSkinThickness() { |
792 | < | if (simParams_->haveSkinThickness()) { |
793 | < | skinThickness_ = simParams_->getSkinThickness(); |
794 | < | } else { |
854 | < | skinThickness_ = 1.0; |
855 | < | sprintf(painCave.errMsg, |
856 | < | "SimInfo Warning: No value was set for the skinThickness.\n" |
857 | < | "\tOpenMD will use a default value of %f Angstroms\n" |
858 | < | "\tfor this simulation\n", skinThickness_); |
859 | < | painCave.isFatal = 0; |
860 | < | simError(); |
861 | < | } |
862 | < | } |
863 | < | |
864 | < | void SimInfo::setupSimType() { |
785 | > | void SimInfo::setupSimVariables() { |
786 | > | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
787 | > | // we only call setAccumulateBoxDipole if the accumulateBoxDipole |
788 | > | // parameter is true |
789 | > | calcBoxDipole_ = false; |
790 | > | if ( simParams_->haveAccumulateBoxDipole() ) |
791 | > | if ( simParams_->getAccumulateBoxDipole() ) { |
792 | > | calcBoxDipole_ = true; |
793 | > | } |
794 | > | |
795 | set<AtomType*>::iterator i; | |
796 | set<AtomType*> atomTypes; | |
797 | < | atomTypes = getSimulatedAtomTypes(); |
798 | < | |
799 | < | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
800 | < | |
801 | < | int usesElectrostatic = 0; |
872 | < | int usesMetallic = 0; |
873 | < | int usesDirectional = 0; |
797 | > | atomTypes = getSimulatedAtomTypes(); |
798 | > | bool usesElectrostatic = false; |
799 | > | bool usesMetallic = false; |
800 | > | bool usesDirectional = false; |
801 | > | bool usesFluctuatingCharges = false; |
802 | //loop over all of the atom types | |
803 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
804 | usesElectrostatic |= (*i)->isElectrostatic(); | |
805 | usesMetallic |= (*i)->isMetal(); | |
806 | usesDirectional |= (*i)->isDirectional(); | |
807 | + | usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); |
808 | } | |
809 | ||
810 | < | #ifdef IS_MPI |
811 | < | int temp; |
810 | > | #ifdef IS_MPI |
811 | > | bool temp; |
812 | temp = usesDirectional; | |
813 | < | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
814 | < | |
813 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, |
814 | > | MPI::LOR); |
815 | > | |
816 | temp = usesMetallic; | |
817 | < | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
818 | < | |
817 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, |
818 | > | MPI::LOR); |
819 | > | |
820 | temp = usesElectrostatic; | |
821 | < | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
821 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, |
822 | > | MPI::LOR); |
823 | > | |
824 | > | temp = usesFluctuatingCharges; |
825 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, |
826 | > | MPI::LOR); |
827 | > | #else |
828 | > | |
829 | > | usesDirectionalAtoms_ = usesDirectional; |
830 | > | usesMetallicAtoms_ = usesMetallic; |
831 | > | usesElectrostaticAtoms_ = usesElectrostatic; |
832 | > | usesFluctuatingCharges_ = usesFluctuatingCharges; |
833 | > | |
834 | #endif | |
835 | < | fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; |
836 | < | fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; |
837 | < | fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; |
838 | < | fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; |
896 | < | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; |
897 | < | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; |
835 | > | |
836 | > | requiresPrepair_ = usesMetallicAtoms_ ? true : false; |
837 | > | requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; |
838 | > | requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; |
839 | } | |
840 | ||
900 | – | void SimInfo::setupFortranSim() { |
901 | – | int isError; |
902 | – | int nExclude, nOneTwo, nOneThree, nOneFour; |
903 | – | vector<int> fortranGlobalGroupMembership; |
904 | – | |
905 | – | notifyFortranSkinThickness(&skinThickness_); |
841 | ||
842 | < | int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
843 | < | int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
844 | < | notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); |
842 | > | vector<int> SimInfo::getGlobalAtomIndices() { |
843 | > | SimInfo::MoleculeIterator mi; |
844 | > | Molecule* mol; |
845 | > | Molecule::AtomIterator ai; |
846 | > | Atom* atom; |
847 | ||
848 | < | isError = 0; |
848 | > | vector<int> GlobalAtomIndices(getNAtoms(), 0); |
849 | > | |
850 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
851 | > | |
852 | > | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
853 | > | GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); |
854 | > | } |
855 | > | } |
856 | > | return GlobalAtomIndices; |
857 | > | } |
858 | ||
859 | < | //globalGroupMembership_ is filled by SimCreator |
860 | < | for (int i = 0; i < nGlobalAtoms_; i++) { |
861 | < | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
859 | > | |
860 | > | vector<int> SimInfo::getGlobalGroupIndices() { |
861 | > | SimInfo::MoleculeIterator mi; |
862 | > | Molecule* mol; |
863 | > | Molecule::CutoffGroupIterator ci; |
864 | > | CutoffGroup* cg; |
865 | > | |
866 | > | vector<int> GlobalGroupIndices; |
867 | > | |
868 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
869 | > | |
870 | > | //local index of cutoff group is trivial, it only depends on the |
871 | > | //order of travesing |
872 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
873 | > | cg = mol->nextCutoffGroup(ci)) { |
874 | > | GlobalGroupIndices.push_back(cg->getGlobalIndex()); |
875 | > | } |
876 | } | |
877 | + | return GlobalGroupIndices; |
878 | + | } |
879 | ||
880 | + | |
881 | + | void SimInfo::prepareTopology() { |
882 | + | |
883 | //calculate mass ratio of cutoff group | |
919 | – | vector<RealType> mfact; |
884 | SimInfo::MoleculeIterator mi; | |
885 | Molecule* mol; | |
886 | Molecule::CutoffGroupIterator ci; | |
# | Line 925 | Line 889 | namespace OpenMD { | |
889 | Atom* atom; | |
890 | RealType totalMass; | |
891 | ||
892 | < | //to avoid memory reallocation, reserve enough space for mfact |
893 | < | mfact.reserve(getNCutoffGroups()); |
892 | > | /** |
893 | > | * The mass factor is the relative mass of an atom to the total |
894 | > | * mass of the cutoff group it belongs to. By default, all atoms |
895 | > | * are their own cutoff groups, and therefore have mass factors of |
896 | > | * 1. We need some special handling for massless atoms, which |
897 | > | * will be treated as carrying the entire mass of the cutoff |
898 | > | * group. |
899 | > | */ |
900 | > | massFactors_.clear(); |
901 | > | massFactors_.resize(getNAtoms(), 1.0); |
902 | ||
903 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
904 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
904 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
905 | > | cg = mol->nextCutoffGroup(ci)) { |
906 | ||
907 | totalMass = cg->getMass(); | |
908 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
909 | // Check for massless groups - set mfact to 1 if true | |
910 | < | if (totalMass != 0) |
911 | < | mfact.push_back(atom->getMass()/totalMass); |
910 | > | if (totalMass != 0) |
911 | > | massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; |
912 | else | |
913 | < | mfact.push_back( 1.0 ); |
913 | > | massFactors_[atom->getLocalIndex()] = 1.0; |
914 | } | |
915 | } | |
916 | } | |
917 | ||
918 | < | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) |
946 | < | vector<int> identArray; |
918 | > | // Build the identArray_ |
919 | ||
920 | < | //to avoid memory reallocation, reserve enough space identArray |
921 | < | identArray.reserve(getNAtoms()); |
950 | < | |
920 | > | identArray_.clear(); |
921 | > | identArray_.reserve(getNAtoms()); |
922 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
923 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
924 | < | identArray.push_back(atom->getIdent()); |
924 | > | identArray_.push_back(atom->getIdent()); |
925 | } | |
926 | } | |
956 | – | |
957 | – | //fill molMembershipArray |
958 | – | //molMembershipArray is filled by SimCreator |
959 | – | vector<int> molMembershipArray(nGlobalAtoms_); |
960 | – | for (int i = 0; i < nGlobalAtoms_; i++) { |
961 | – | molMembershipArray[i] = globalMolMembership_[i] + 1; |
962 | – | } |
927 | ||
928 | < | //setup fortran simulation |
928 | > | //scan topology |
929 | ||
966 | – | nExclude = excludedInteractions_.getSize(); |
967 | – | nOneTwo = oneTwoInteractions_.getSize(); |
968 | – | nOneThree = oneThreeInteractions_.getSize(); |
969 | – | nOneFour = oneFourInteractions_.getSize(); |
970 | – | |
930 | int* excludeList = excludedInteractions_.getPairList(); | |
931 | int* oneTwoList = oneTwoInteractions_.getPairList(); | |
932 | int* oneThreeList = oneThreeInteractions_.getPairList(); | |
933 | int* oneFourList = oneFourInteractions_.getPairList(); | |
934 | ||
935 | < | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], |
977 | < | &nExclude, excludeList, |
978 | < | &nOneTwo, oneTwoList, |
979 | < | &nOneThree, oneThreeList, |
980 | < | &nOneFour, oneFourList, |
981 | < | &molMembershipArray[0], &mfact[0], &nCutoffGroups_, |
982 | < | &fortranGlobalGroupMembership[0], &isError); |
983 | < | |
984 | < | if( isError ){ |
985 | < | |
986 | < | sprintf( painCave.errMsg, |
987 | < | "There was an error setting the simulation information in fortran.\n" ); |
988 | < | painCave.isFatal = 1; |
989 | < | painCave.severity = OPENMD_ERROR; |
990 | < | simError(); |
991 | < | } |
992 | < | |
993 | < | |
994 | < | sprintf( checkPointMsg, |
995 | < | "succesfully sent the simulation information to fortran.\n"); |
996 | < | |
997 | < | errorCheckPoint(); |
998 | < | |
999 | < | // Setup number of neighbors in neighbor list if present |
1000 | < | if (simParams_->haveNeighborListNeighbors()) { |
1001 | < | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
1002 | < | setNeighbors(&nlistNeighbors); |
1003 | < | } |
1004 | < | |
1005 | < | |
935 | > | topologyDone_ = true; |
936 | } | |
937 | ||
1008 | – | |
1009 | – | void SimInfo::setupFortranParallel() { |
1010 | – | #ifdef IS_MPI |
1011 | – | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex |
1012 | – | vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
1013 | – | vector<int> localToGlobalCutoffGroupIndex; |
1014 | – | SimInfo::MoleculeIterator mi; |
1015 | – | Molecule::AtomIterator ai; |
1016 | – | Molecule::CutoffGroupIterator ci; |
1017 | – | Molecule* mol; |
1018 | – | Atom* atom; |
1019 | – | CutoffGroup* cg; |
1020 | – | mpiSimData parallelData; |
1021 | – | int isError; |
1022 | – | |
1023 | – | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
1024 | – | |
1025 | – | //local index(index in DataStorge) of atom is important |
1026 | – | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
1027 | – | localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; |
1028 | – | } |
1029 | – | |
1030 | – | //local index of cutoff group is trivial, it only depends on the order of travesing |
1031 | – | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
1032 | – | localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); |
1033 | – | } |
1034 | – | |
1035 | – | } |
1036 | – | |
1037 | – | //fill up mpiSimData struct |
1038 | – | parallelData.nMolGlobal = getNGlobalMolecules(); |
1039 | – | parallelData.nMolLocal = getNMolecules(); |
1040 | – | parallelData.nAtomsGlobal = getNGlobalAtoms(); |
1041 | – | parallelData.nAtomsLocal = getNAtoms(); |
1042 | – | parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); |
1043 | – | parallelData.nGroupsLocal = getNCutoffGroups(); |
1044 | – | parallelData.myNode = worldRank; |
1045 | – | MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); |
1046 | – | |
1047 | – | //pass mpiSimData struct and index arrays to fortran |
1048 | – | setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), |
1049 | – | &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), |
1050 | – | &localToGlobalCutoffGroupIndex[0], &isError); |
1051 | – | |
1052 | – | if (isError) { |
1053 | – | sprintf(painCave.errMsg, |
1054 | – | "mpiRefresh errror: fortran didn't like something we gave it.\n"); |
1055 | – | painCave.isFatal = 1; |
1056 | – | simError(); |
1057 | – | } |
1058 | – | |
1059 | – | sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
1060 | – | errorCheckPoint(); |
1061 | – | |
1062 | – | #endif |
1063 | – | } |
1064 | – | |
1065 | – | |
1066 | – | void SimInfo::setupSwitchingFunction() { |
1067 | – | |
1068 | – | } |
1069 | – | |
1070 | – | void SimInfo::setupAccumulateBoxDipole() { |
1071 | – | |
1072 | – | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1073 | – | if ( simParams_->haveAccumulateBoxDipole() ) |
1074 | – | if ( simParams_->getAccumulateBoxDipole() ) { |
1075 | – | calcBoxDipole_ = true; |
1076 | – | } |
1077 | – | |
1078 | – | } |
1079 | – | |
938 | void SimInfo::addProperty(GenericData* genData) { | |
939 | properties_.addProperty(genData); | |
940 | } | |
# | Line 1111 | Line 969 | namespace OpenMD { | |
969 | Molecule* mol; | |
970 | RigidBody* rb; | |
971 | Atom* atom; | |
972 | + | CutoffGroup* cg; |
973 | SimInfo::MoleculeIterator mi; | |
974 | Molecule::RigidBodyIterator rbIter; | |
975 | < | Molecule::AtomIterator atomIter;; |
975 | > | Molecule::AtomIterator atomIter; |
976 | > | Molecule::CutoffGroupIterator cgIter; |
977 | ||
978 | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
979 | ||
980 | < | for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { |
980 | > | for (atom = mol->beginAtom(atomIter); atom != NULL; |
981 | > | atom = mol->nextAtom(atomIter)) { |
982 | atom->setSnapshotManager(sman_); | |
983 | } | |
984 | ||
985 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
985 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
986 | > | rb = mol->nextRigidBody(rbIter)) { |
987 | rb->setSnapshotManager(sman_); | |
988 | } | |
989 | + | |
990 | + | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; |
991 | + | cg = mol->nextCutoffGroup(cgIter)) { |
992 | + | cg->setSnapshotManager(sman_); |
993 | + | } |
994 | } | |
995 | ||
996 | } | |
1130 | – | |
1131 | – | Vector3d SimInfo::getComVel(){ |
1132 | – | SimInfo::MoleculeIterator i; |
1133 | – | Molecule* mol; |
1134 | – | |
1135 | – | Vector3d comVel(0.0); |
1136 | – | RealType totalMass = 0.0; |
1137 | – | |
1138 | – | |
1139 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1140 | – | RealType mass = mol->getMass(); |
1141 | – | totalMass += mass; |
1142 | – | comVel += mass * mol->getComVel(); |
1143 | – | } |
1144 | – | |
1145 | – | #ifdef IS_MPI |
1146 | – | RealType tmpMass = totalMass; |
1147 | – | Vector3d tmpComVel(comVel); |
1148 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1149 | – | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1150 | – | #endif |
1151 | – | |
1152 | – | comVel /= totalMass; |
1153 | – | |
1154 | – | return comVel; |
1155 | – | } |
1156 | – | |
1157 | – | Vector3d SimInfo::getCom(){ |
1158 | – | SimInfo::MoleculeIterator i; |
1159 | – | Molecule* mol; |
1160 | – | |
1161 | – | Vector3d com(0.0); |
1162 | – | RealType totalMass = 0.0; |
1163 | – | |
1164 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1165 | – | RealType mass = mol->getMass(); |
1166 | – | totalMass += mass; |
1167 | – | com += mass * mol->getCom(); |
1168 | – | } |
997 | ||
1170 | – | #ifdef IS_MPI |
1171 | – | RealType tmpMass = totalMass; |
1172 | – | Vector3d tmpCom(com); |
1173 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1174 | – | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1175 | – | #endif |
998 | ||
1177 | – | com /= totalMass; |
1178 | – | |
1179 | – | return com; |
1180 | – | |
1181 | – | } |
1182 | – | |
999 | ostream& operator <<(ostream& o, SimInfo& info) { | |
1000 | ||
1001 | return o; | |
1002 | } | |
1003 | ||
1004 | < | |
1189 | < | /* |
1190 | < | Returns center of mass and center of mass velocity in one function call. |
1191 | < | */ |
1192 | < | |
1193 | < | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
1194 | < | SimInfo::MoleculeIterator i; |
1195 | < | Molecule* mol; |
1196 | < | |
1197 | < | |
1198 | < | RealType totalMass = 0.0; |
1199 | < | |
1200 | < | |
1201 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1202 | < | RealType mass = mol->getMass(); |
1203 | < | totalMass += mass; |
1204 | < | com += mass * mol->getCom(); |
1205 | < | comVel += mass * mol->getComVel(); |
1206 | < | } |
1207 | < | |
1208 | < | #ifdef IS_MPI |
1209 | < | RealType tmpMass = totalMass; |
1210 | < | Vector3d tmpCom(com); |
1211 | < | Vector3d tmpComVel(comVel); |
1212 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1213 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1214 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1215 | < | #endif |
1216 | < | |
1217 | < | com /= totalMass; |
1218 | < | comVel /= totalMass; |
1219 | < | } |
1220 | < | |
1221 | < | /* |
1222 | < | Return intertia tensor for entire system and angular momentum Vector. |
1223 | < | |
1224 | < | |
1225 | < | [ Ixx -Ixy -Ixz ] |
1226 | < | J =| -Iyx Iyy -Iyz | |
1227 | < | [ -Izx -Iyz Izz ] |
1228 | < | */ |
1229 | < | |
1230 | < | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1231 | < | |
1232 | < | |
1233 | < | RealType xx = 0.0; |
1234 | < | RealType yy = 0.0; |
1235 | < | RealType zz = 0.0; |
1236 | < | RealType xy = 0.0; |
1237 | < | RealType xz = 0.0; |
1238 | < | RealType yz = 0.0; |
1239 | < | Vector3d com(0.0); |
1240 | < | Vector3d comVel(0.0); |
1241 | < | |
1242 | < | getComAll(com, comVel); |
1243 | < | |
1244 | < | SimInfo::MoleculeIterator i; |
1245 | < | Molecule* mol; |
1246 | < | |
1247 | < | Vector3d thisq(0.0); |
1248 | < | Vector3d thisv(0.0); |
1249 | < | |
1250 | < | RealType thisMass = 0.0; |
1251 | < | |
1252 | < | |
1253 | < | |
1254 | < | |
1255 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1256 | < | |
1257 | < | thisq = mol->getCom()-com; |
1258 | < | thisv = mol->getComVel()-comVel; |
1259 | < | thisMass = mol->getMass(); |
1260 | < | // Compute moment of intertia coefficients. |
1261 | < | xx += thisq[0]*thisq[0]*thisMass; |
1262 | < | yy += thisq[1]*thisq[1]*thisMass; |
1263 | < | zz += thisq[2]*thisq[2]*thisMass; |
1264 | < | |
1265 | < | // compute products of intertia |
1266 | < | xy += thisq[0]*thisq[1]*thisMass; |
1267 | < | xz += thisq[0]*thisq[2]*thisMass; |
1268 | < | yz += thisq[1]*thisq[2]*thisMass; |
1269 | < | |
1270 | < | angularMomentum += cross( thisq, thisv ) * thisMass; |
1271 | < | |
1272 | < | } |
1273 | < | |
1274 | < | |
1275 | < | inertiaTensor(0,0) = yy + zz; |
1276 | < | inertiaTensor(0,1) = -xy; |
1277 | < | inertiaTensor(0,2) = -xz; |
1278 | < | inertiaTensor(1,0) = -xy; |
1279 | < | inertiaTensor(1,1) = xx + zz; |
1280 | < | inertiaTensor(1,2) = -yz; |
1281 | < | inertiaTensor(2,0) = -xz; |
1282 | < | inertiaTensor(2,1) = -yz; |
1283 | < | inertiaTensor(2,2) = xx + yy; |
1284 | < | |
1285 | < | #ifdef IS_MPI |
1286 | < | Mat3x3d tmpI(inertiaTensor); |
1287 | < | Vector3d tmpAngMom; |
1288 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1289 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1290 | < | #endif |
1291 | < | |
1292 | < | return; |
1293 | < | } |
1294 | < | |
1295 | < | //Returns the angular momentum of the system |
1296 | < | Vector3d SimInfo::getAngularMomentum(){ |
1297 | < | |
1298 | < | Vector3d com(0.0); |
1299 | < | Vector3d comVel(0.0); |
1300 | < | Vector3d angularMomentum(0.0); |
1301 | < | |
1302 | < | getComAll(com,comVel); |
1303 | < | |
1304 | < | SimInfo::MoleculeIterator i; |
1305 | < | Molecule* mol; |
1306 | < | |
1307 | < | Vector3d thisr(0.0); |
1308 | < | Vector3d thisp(0.0); |
1309 | < | |
1310 | < | RealType thisMass; |
1311 | < | |
1312 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1313 | < | thisMass = mol->getMass(); |
1314 | < | thisr = mol->getCom()-com; |
1315 | < | thisp = (mol->getComVel()-comVel)*thisMass; |
1316 | < | |
1317 | < | angularMomentum += cross( thisr, thisp ); |
1318 | < | |
1319 | < | } |
1320 | < | |
1321 | < | #ifdef IS_MPI |
1322 | < | Vector3d tmpAngMom; |
1323 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1324 | < | #endif |
1325 | < | |
1326 | < | return angularMomentum; |
1327 | < | } |
1328 | < | |
1004 | > | |
1005 | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | |
1006 | < | return IOIndexToIntegrableObject.at(index); |
1006 | > | if (index >= int(IOIndexToIntegrableObject.size())) { |
1007 | > | sprintf(painCave.errMsg, |
1008 | > | "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n" |
1009 | > | "\tindex exceeds number of known objects!\n"); |
1010 | > | painCave.isFatal = 1; |
1011 | > | simError(); |
1012 | > | return NULL; |
1013 | > | } else |
1014 | > | return IOIndexToIntegrableObject.at(index); |
1015 | } | |
1016 | ||
1017 | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { | |
1018 | IOIndexToIntegrableObject= v; | |
1019 | } | |
1020 | ||
1337 | – | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1338 | – | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1339 | – | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1340 | – | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1341 | – | */ |
1342 | – | void SimInfo::getGyrationalVolume(RealType &volume){ |
1343 | – | Mat3x3d intTensor; |
1344 | – | RealType det; |
1345 | – | Vector3d dummyAngMom; |
1346 | – | RealType sysconstants; |
1347 | – | RealType geomCnst; |
1348 | – | |
1349 | – | geomCnst = 3.0/2.0; |
1350 | – | /* Get the inertial tensor and angular momentum for free*/ |
1351 | – | getInertiaTensor(intTensor,dummyAngMom); |
1352 | – | |
1353 | – | det = intTensor.determinant(); |
1354 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1355 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1356 | – | return; |
1357 | – | } |
1358 | – | |
1359 | – | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1360 | – | Mat3x3d intTensor; |
1361 | – | Vector3d dummyAngMom; |
1362 | – | RealType sysconstants; |
1363 | – | RealType geomCnst; |
1364 | – | |
1365 | – | geomCnst = 3.0/2.0; |
1366 | – | /* Get the inertial tensor and angular momentum for free*/ |
1367 | – | getInertiaTensor(intTensor,dummyAngMom); |
1368 | – | |
1369 | – | detI = intTensor.determinant(); |
1370 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1371 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1372 | – | return; |
1373 | – | } |
1374 | – | /* |
1375 | – | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { |
1376 | – | assert( v.size() == nAtoms_ + nRigidBodies_); |
1377 | – | sdByGlobalIndex_ = v; |
1378 | – | } |
1379 | – | |
1380 | – | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1381 | – | //assert(index < nAtoms_ + nRigidBodies_); |
1382 | – | return sdByGlobalIndex_.at(index); |
1383 | – | } |
1384 | – | */ |
1021 | int SimInfo::getNGlobalConstraints() { | |
1022 | int nGlobalConstraints; | |
1023 | #ifdef IS_MPI | |
1024 | < | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
1025 | < | MPI_COMM_WORLD); |
1024 | > | MPI::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1, |
1025 | > | MPI::INT, MPI::SUM); |
1026 | #else | |
1027 | nGlobalConstraints = nConstraints_; | |
1028 | #endif |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |