# | 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/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" |
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" | |
63 | < | #include "nonbonded/InteractionManager.hpp" |
68 | < | |
69 | < | |
63 | > | #include "nonbonded/SwitchingFunction.hpp" |
64 | #ifdef IS_MPI | |
65 | < | #include "UseTheForce/mpiComponentPlan.h" |
66 | < | #include "UseTheForce/DarkSide/simParallel_interface.h" |
73 | < | #endif |
65 | > | #include <mpi.h> |
66 | > | #endif |
67 | ||
68 | using namespace std; | |
69 | namespace OpenMD { | |
# | Line 79 | 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 136 | Line 129 | namespace OpenMD { | |
129 | //equal to the total number of atoms minus number of atoms belong to | |
130 | //cutoff group defined in meta-data file plus the number of cutoff | |
131 | //groups defined in meta-data file | |
132 | + | |
133 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
134 | ||
135 | //every free atom (atom does not belong to rigid bodies) is an | |
# | Line 231 | Line 225 | namespace OpenMD { | |
225 | ||
226 | ||
227 | void SimInfo::calcNdf() { | |
228 | < | int ndf_local; |
228 | > | int ndf_local, nfq_local; |
229 | MoleculeIterator i; | |
230 | vector<StuntDouble*>::iterator j; | |
231 | + | vector<Atom*>::iterator k; |
232 | + | |
233 | Molecule* mol; | |
234 | StuntDouble* integrableObject; | |
235 | + | Atom* atom; |
236 | ||
237 | ndf_local = 0; | |
238 | + | nfq_local = 0; |
239 | ||
240 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
241 | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; | |
# | Line 252 | Line 250 | namespace OpenMD { | |
250 | ndf_local += 3; | |
251 | } | |
252 | } | |
255 | – | |
253 | } | |
254 | + | for (atom = mol->beginFluctuatingCharge(k); atom != NULL; |
255 | + | atom = mol->nextFluctuatingCharge(k)) { |
256 | + | if (atom->isFluctuatingCharge()) { |
257 | + | nfq_local++; |
258 | + | } |
259 | + | } |
260 | } | |
261 | ||
262 | // n_constraints is local, so subtract them on each processor | |
# | Line 261 | Line 264 | namespace OpenMD { | |
264 | ||
265 | #ifdef IS_MPI | |
266 | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | |
267 | + | MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
268 | #else | |
269 | ndf_ = ndf_local; | |
270 | + | nGlobalFluctuatingCharges_ = nfq_local; |
271 | #endif | |
272 | ||
273 | // nZconstraints_ is global, as are the 3 COM translations for the | |
# | Line 279 | Line 284 | namespace OpenMD { | |
284 | #endif | |
285 | return fdf_; | |
286 | } | |
287 | + | |
288 | + | unsigned int SimInfo::getNLocalCutoffGroups(){ |
289 | + | int nLocalCutoffAtoms = 0; |
290 | + | Molecule* mol; |
291 | + | MoleculeIterator mi; |
292 | + | CutoffGroup* cg; |
293 | + | Molecule::CutoffGroupIterator ci; |
294 | ||
295 | + | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
296 | + | |
297 | + | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
298 | + | cg = mol->nextCutoffGroup(ci)) { |
299 | + | nLocalCutoffAtoms += cg->getNumAtom(); |
300 | + | |
301 | + | } |
302 | + | } |
303 | + | |
304 | + | return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; |
305 | + | } |
306 | + | |
307 | void SimInfo::calcNdfRaw() { | |
308 | int ndfRaw_local; | |
309 | ||
# | Line 656 | Line 680 | namespace OpenMD { | |
680 | molStampIds_.insert(molStampIds_.end(), nmol, curStampId); | |
681 | } | |
682 | ||
659 | – | void SimInfo::update() { |
660 | – | |
661 | – | setupSimType(); |
662 | – | setupCutoffRadius(); |
663 | – | setupSwitchingRadius(); |
664 | – | setupCutoffMethod(); |
665 | – | setupSkinThickness(); |
666 | – | setupSwitchingFunction(); |
667 | – | setupAccumulateBoxDipole(); |
683 | ||
684 | < | #ifdef IS_MPI |
685 | < | setupFortranParallel(); |
686 | < | #endif |
687 | < | setupFortranSim(); |
688 | < | fortranInitialized_ = true; |
689 | < | |
684 | > | /** |
685 | > | * update |
686 | > | * |
687 | > | * Performs the global checks and variable settings after the |
688 | > | * objects have been created. |
689 | > | * |
690 | > | */ |
691 | > | void SimInfo::update() { |
692 | > | setupSimVariables(); |
693 | calcNdf(); | |
694 | calcNdfRaw(); | |
695 | calcNdfTrans(); | |
696 | } | |
697 | ||
698 | + | /** |
699 | + | * getSimulatedAtomTypes |
700 | + | * |
701 | + | * Returns an STL set of AtomType* that are actually present in this |
702 | + | * simulation. Must query all processors to assemble this information. |
703 | + | * |
704 | + | */ |
705 | set<AtomType*> SimInfo::getSimulatedAtomTypes() { | |
706 | SimInfo::MoleculeIterator mi; | |
707 | Molecule* mol; | |
# | Line 684 | Line 709 | namespace OpenMD { | |
709 | Atom* atom; | |
710 | set<AtomType*> atomTypes; | |
711 | ||
712 | < | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
713 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
712 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
713 | > | for(atom = mol->beginAtom(ai); atom != NULL; |
714 | > | atom = mol->nextAtom(ai)) { |
715 | atomTypes.insert(atom->getAtomType()); | |
716 | } | |
717 | } | |
718 | < | return atomTypes; |
719 | < | } |
718 | > | |
719 | > | #ifdef IS_MPI |
720 | ||
721 | < | /** |
722 | < | * setupCutoffRadius |
697 | < | * |
698 | < | * If the cutoffRadius was explicitly set, use that value. |
699 | < | * If the cutoffRadius was not explicitly set: |
700 | < | * Are there electrostatic atoms? Use 12.0 Angstroms. |
701 | < | * No electrostatic atoms? Poll the atom types present in the |
702 | < | * simulation for suggested cutoff values (e.g. 2.5 * sigma). |
703 | < | * Use the maximum suggested value that was found. |
704 | < | */ |
705 | < | void SimInfo::setupCutoffRadius() { |
721 | > | // loop over the found atom types on this processor, and add their |
722 | > | // numerical idents to a vector: |
723 | ||
724 | < | if (simParams_->haveCutoffRadius()) { |
725 | < | cutoffRadius_ = simParams_->getCutoffRadius(); |
726 | < | } else { |
727 | < | if (usesElectrostaticAtoms_) { |
711 | < | sprintf(painCave.errMsg, |
712 | < | "SimInfo Warning: No value was set for the cutoffRadius.\n" |
713 | < | "\tOpenMD will use a default value of 12.0 angstroms" |
714 | < | "\tfor the cutoffRadius.\n"); |
715 | < | painCave.isFatal = 0; |
716 | < | simError(); |
717 | < | cutoffRadius_ = 12.0; |
718 | < | } else { |
719 | < | RealType thisCut; |
720 | < | set<AtomType*>::iterator i; |
721 | < | set<AtomType*> atomTypes; |
722 | < | atomTypes = getSimulatedAtomTypes(); |
723 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
724 | < | thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); |
725 | < | cutoffRadius_ = max(thisCut, cutoffRadius_); |
726 | < | } |
727 | < | sprintf(painCave.errMsg, |
728 | < | "SimInfo Warning: No value was set for the cutoffRadius.\n" |
729 | < | "\tOpenMD will use %lf angstroms.\n", |
730 | < | cutoffRadius_); |
731 | < | painCave.isFatal = 0; |
732 | < | simError(); |
733 | < | } |
734 | < | } |
724 | > | vector<int> foundTypes; |
725 | > | set<AtomType*>::iterator i; |
726 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) |
727 | > | foundTypes.push_back( (*i)->getIdent() ); |
728 | ||
729 | < | InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); |
730 | < | } |
729 | > | // count_local holds the number of found types on this processor |
730 | > | int count_local = foundTypes.size(); |
731 | > | |
732 | > | int nproc = MPI::COMM_WORLD.Get_size(); |
733 | > | |
734 | > | // we need arrays to hold the counts and displacement vectors for |
735 | > | // all processors |
736 | > | vector<int> counts(nproc, 0); |
737 | > | vector<int> disps(nproc, 0); |
738 | > | |
739 | > | // fill the counts array |
740 | > | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], |
741 | > | 1, MPI::INT); |
742 | ||
743 | < | /** |
744 | < | * setupSwitchingRadius |
745 | < | * |
746 | < | * If the switchingRadius was explicitly set, use that value (but check it) |
747 | < | * If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ |
748 | < | */ |
749 | < | void SimInfo::setupSwitchingRadius() { |
743 | > | // use the processor counts to compute the displacement array |
744 | > | disps[0] = 0; |
745 | > | int totalCount = counts[0]; |
746 | > | for (int iproc = 1; iproc < nproc; iproc++) { |
747 | > | disps[iproc] = disps[iproc-1] + counts[iproc-1]; |
748 | > | totalCount += counts[iproc]; |
749 | > | } |
750 | > | |
751 | > | // we need a (possibly redundant) set of all found types: |
752 | > | vector<int> ftGlobal(totalCount); |
753 | ||
754 | < | if (simParams_->haveSwitchingRadius()) { |
755 | < | switchingRadius_ = simParams_->getSwitchingRadius(); |
756 | < | if (switchingRadius_ > cutoffRadius_) { |
757 | < | sprintf(painCave.errMsg, |
751 | < | "SimInfo Error: switchingRadius (%f) is larger than cutoffRadius(%f)\n", |
752 | < | switchingRadius_, cutoffRadius_); |
753 | < | painCave.isFatal = 1; |
754 | < | simError(); |
754 | > | // now spray out the foundTypes to all the other processors: |
755 | > | MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, |
756 | > | &ftGlobal[0], &counts[0], &disps[0], |
757 | > | MPI::INT); |
758 | ||
759 | < | } |
757 | < | } else { |
758 | < | switchingRadius_ = 0.85 * cutoffRadius_; |
759 | < | sprintf(painCave.errMsg, |
760 | < | "SimInfo Warning: No value was set for the switchingRadius.\n" |
761 | < | "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" |
762 | < | "\tswitchingRadius = %f. for this simulation\n", switchingRadius_); |
763 | < | painCave.isFatal = 0; |
764 | < | simError(); |
765 | < | } |
766 | < | InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); |
767 | < | } |
759 | > | vector<int>::iterator j; |
760 | ||
761 | < | /** |
762 | < | * setupSkinThickness |
763 | < | * |
764 | < | * If the skinThickness was explicitly set, use that value (but check it) |
765 | < | * If the skinThickness was not explicitly set: use 1.0 angstroms |
766 | < | */ |
767 | < | void SimInfo::setupSkinThickness() { |
768 | < | if (simParams_->haveSkinThickness()) { |
769 | < | skinThickness_ = simParams_->getSkinThickness(); |
770 | < | } else { |
771 | < | skinThickness_ = 1.0; |
772 | < | sprintf(painCave.errMsg, |
773 | < | "SimInfo Warning: No value was set for the skinThickness.\n" |
774 | < | "\tOpenMD will use a default value of %f Angstroms\n" |
775 | < | "\tfor this simulation\n", skinThickness_); |
776 | < | painCave.isFatal = 0; |
785 | < | simError(); |
786 | < | } |
761 | > | // foundIdents is a stl set, so inserting an already found ident |
762 | > | // will have no effect. |
763 | > | set<int> foundIdents; |
764 | > | |
765 | > | for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) |
766 | > | foundIdents.insert((*j)); |
767 | > | |
768 | > | // now iterate over the foundIdents and get the actual atom types |
769 | > | // that correspond to these: |
770 | > | set<int>::iterator it; |
771 | > | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
772 | > | atomTypes.insert( forceField_->getAtomType((*it)) ); |
773 | > | |
774 | > | #endif |
775 | > | |
776 | > | return atomTypes; |
777 | } | |
778 | ||
779 | < | void SimInfo::setupSimType() { |
779 | > | void SimInfo::setupSimVariables() { |
780 | > | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
781 | > | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
782 | > | calcBoxDipole_ = false; |
783 | > | if ( simParams_->haveAccumulateBoxDipole() ) |
784 | > | if ( simParams_->getAccumulateBoxDipole() ) { |
785 | > | calcBoxDipole_ = true; |
786 | > | } |
787 | > | |
788 | set<AtomType*>::iterator i; | |
789 | set<AtomType*> atomTypes; | |
790 | < | atomTypes = getSimulatedAtomTypes(); |
793 | < | |
794 | < | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
795 | < | |
790 | > | atomTypes = getSimulatedAtomTypes(); |
791 | int usesElectrostatic = 0; | |
792 | int usesMetallic = 0; | |
793 | int usesDirectional = 0; | |
794 | + | int usesFluctuatingCharges = 0; |
795 | //loop over all of the atom types | |
796 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
797 | usesElectrostatic |= (*i)->isElectrostatic(); | |
798 | usesMetallic |= (*i)->isMetal(); | |
799 | usesDirectional |= (*i)->isDirectional(); | |
800 | + | usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); |
801 | } | |
802 | < | |
802 | > | |
803 | #ifdef IS_MPI | |
804 | int temp; | |
805 | temp = usesDirectional; | |
806 | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
807 | < | |
807 | > | |
808 | temp = usesMetallic; | |
809 | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
810 | < | |
810 | > | |
811 | temp = usesElectrostatic; | |
812 | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
813 | + | |
814 | + | temp = usesFluctuatingCharges; |
815 | + | MPI_Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
816 | + | #else |
817 | + | |
818 | + | usesDirectionalAtoms_ = usesDirectional; |
819 | + | usesMetallicAtoms_ = usesMetallic; |
820 | + | usesElectrostaticAtoms_ = usesElectrostatic; |
821 | + | usesFluctuatingCharges_ = usesFluctuatingCharges; |
822 | + | |
823 | #endif | |
824 | < | fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; |
825 | < | fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; |
826 | < | fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; |
827 | < | fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; |
821 | < | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; |
822 | < | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; |
824 | > | |
825 | > | requiresPrepair_ = usesMetallicAtoms_ ? true : false; |
826 | > | requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; |
827 | > | requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; |
828 | } | |
829 | ||
830 | < | void SimInfo::setupFortranSim() { |
831 | < | int isError; |
832 | < | int nExclude, nOneTwo, nOneThree, nOneFour; |
833 | < | vector<int> fortranGlobalGroupMembership; |
830 | > | |
831 | > | vector<int> SimInfo::getGlobalAtomIndices() { |
832 | > | SimInfo::MoleculeIterator mi; |
833 | > | Molecule* mol; |
834 | > | Molecule::AtomIterator ai; |
835 | > | Atom* atom; |
836 | > | |
837 | > | vector<int> GlobalAtomIndices(getNAtoms(), 0); |
838 | ||
839 | < | notifyFortranSkinThickness(&skinThickness_); |
839 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
840 | > | |
841 | > | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
842 | > | GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); |
843 | > | } |
844 | > | } |
845 | > | return GlobalAtomIndices; |
846 | > | } |
847 | ||
832 | – | int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
833 | – | int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
834 | – | notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); |
848 | ||
849 | < | isError = 0; |
849 | > | vector<int> SimInfo::getGlobalGroupIndices() { |
850 | > | SimInfo::MoleculeIterator mi; |
851 | > | Molecule* mol; |
852 | > | Molecule::CutoffGroupIterator ci; |
853 | > | CutoffGroup* cg; |
854 | ||
855 | < | //globalGroupMembership_ is filled by SimCreator |
856 | < | for (int i = 0; i < nGlobalAtoms_; i++) { |
857 | < | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
855 | > | vector<int> GlobalGroupIndices; |
856 | > | |
857 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
858 | > | |
859 | > | //local index of cutoff group is trivial, it only depends on the |
860 | > | //order of travesing |
861 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
862 | > | cg = mol->nextCutoffGroup(ci)) { |
863 | > | GlobalGroupIndices.push_back(cg->getGlobalIndex()); |
864 | > | } |
865 | } | |
866 | + | return GlobalGroupIndices; |
867 | + | } |
868 | ||
869 | + | |
870 | + | void SimInfo::prepareTopology() { |
871 | + | int nExclude, nOneTwo, nOneThree, nOneFour; |
872 | + | |
873 | //calculate mass ratio of cutoff group | |
844 | – | vector<RealType> mfact; |
874 | SimInfo::MoleculeIterator mi; | |
875 | Molecule* mol; | |
876 | Molecule::CutoffGroupIterator ci; | |
# | Line 850 | Line 879 | namespace OpenMD { | |
879 | Atom* atom; | |
880 | RealType totalMass; | |
881 | ||
882 | < | //to avoid memory reallocation, reserve enough space for mfact |
883 | < | mfact.reserve(getNCutoffGroups()); |
882 | > | /** |
883 | > | * The mass factor is the relative mass of an atom to the total |
884 | > | * mass of the cutoff group it belongs to. By default, all atoms |
885 | > | * are their own cutoff groups, and therefore have mass factors of |
886 | > | * 1. We need some special handling for massless atoms, which |
887 | > | * will be treated as carrying the entire mass of the cutoff |
888 | > | * group. |
889 | > | */ |
890 | > | massFactors_.clear(); |
891 | > | massFactors_.resize(getNAtoms(), 1.0); |
892 | ||
893 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
894 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
894 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
895 | > | cg = mol->nextCutoffGroup(ci)) { |
896 | ||
897 | totalMass = cg->getMass(); | |
898 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
899 | // Check for massless groups - set mfact to 1 if true | |
900 | < | if (totalMass != 0) |
901 | < | mfact.push_back(atom->getMass()/totalMass); |
900 | > | if (totalMass != 0) |
901 | > | massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; |
902 | else | |
903 | < | mfact.push_back( 1.0 ); |
903 | > | massFactors_[atom->getLocalIndex()] = 1.0; |
904 | } | |
905 | } | |
906 | } | |
907 | ||
908 | < | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) |
871 | < | vector<int> identArray; |
908 | > | // Build the identArray_ |
909 | ||
910 | < | //to avoid memory reallocation, reserve enough space identArray |
911 | < | identArray.reserve(getNAtoms()); |
875 | < | |
910 | > | identArray_.clear(); |
911 | > | identArray_.reserve(getNAtoms()); |
912 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
913 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
914 | < | identArray.push_back(atom->getIdent()); |
914 | > | identArray_.push_back(atom->getIdent()); |
915 | } | |
916 | } | |
881 | – | |
882 | – | //fill molMembershipArray |
883 | – | //molMembershipArray is filled by SimCreator |
884 | – | vector<int> molMembershipArray(nGlobalAtoms_); |
885 | – | for (int i = 0; i < nGlobalAtoms_; i++) { |
886 | – | molMembershipArray[i] = globalMolMembership_[i] + 1; |
887 | – | } |
917 | ||
918 | < | //setup fortran simulation |
918 | > | //scan topology |
919 | ||
920 | nExclude = excludedInteractions_.getSize(); | |
921 | nOneTwo = oneTwoInteractions_.getSize(); | |
# | Line 898 | Line 927 | namespace OpenMD { | |
927 | int* oneThreeList = oneThreeInteractions_.getPairList(); | |
928 | int* oneFourList = oneFourInteractions_.getPairList(); | |
929 | ||
930 | < | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], |
902 | < | &nExclude, excludeList, |
903 | < | &nOneTwo, oneTwoList, |
904 | < | &nOneThree, oneThreeList, |
905 | < | &nOneFour, oneFourList, |
906 | < | &molMembershipArray[0], &mfact[0], &nCutoffGroups_, |
907 | < | &fortranGlobalGroupMembership[0], &isError); |
908 | < | |
909 | < | if( isError ){ |
910 | < | |
911 | < | sprintf( painCave.errMsg, |
912 | < | "There was an error setting the simulation information in fortran.\n" ); |
913 | < | painCave.isFatal = 1; |
914 | < | painCave.severity = OPENMD_ERROR; |
915 | < | simError(); |
916 | < | } |
917 | < | |
918 | < | |
919 | < | sprintf( checkPointMsg, |
920 | < | "succesfully sent the simulation information to fortran.\n"); |
921 | < | |
922 | < | errorCheckPoint(); |
923 | < | |
924 | < | // Setup number of neighbors in neighbor list if present |
925 | < | if (simParams_->haveNeighborListNeighbors()) { |
926 | < | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
927 | < | setNeighbors(&nlistNeighbors); |
928 | < | } |
929 | < | |
930 | < | |
930 | > | topologyDone_ = true; |
931 | } | |
932 | ||
933 | – | |
934 | – | void SimInfo::setupFortranParallel() { |
935 | – | #ifdef IS_MPI |
936 | – | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex |
937 | – | vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
938 | – | vector<int> localToGlobalCutoffGroupIndex; |
939 | – | SimInfo::MoleculeIterator mi; |
940 | – | Molecule::AtomIterator ai; |
941 | – | Molecule::CutoffGroupIterator ci; |
942 | – | Molecule* mol; |
943 | – | Atom* atom; |
944 | – | CutoffGroup* cg; |
945 | – | mpiSimData parallelData; |
946 | – | int isError; |
947 | – | |
948 | – | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
949 | – | |
950 | – | //local index(index in DataStorge) of atom is important |
951 | – | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
952 | – | localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; |
953 | – | } |
954 | – | |
955 | – | //local index of cutoff group is trivial, it only depends on the order of travesing |
956 | – | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
957 | – | localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); |
958 | – | } |
959 | – | |
960 | – | } |
961 | – | |
962 | – | //fill up mpiSimData struct |
963 | – | parallelData.nMolGlobal = getNGlobalMolecules(); |
964 | – | parallelData.nMolLocal = getNMolecules(); |
965 | – | parallelData.nAtomsGlobal = getNGlobalAtoms(); |
966 | – | parallelData.nAtomsLocal = getNAtoms(); |
967 | – | parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); |
968 | – | parallelData.nGroupsLocal = getNCutoffGroups(); |
969 | – | parallelData.myNode = worldRank; |
970 | – | MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); |
971 | – | |
972 | – | //pass mpiSimData struct and index arrays to fortran |
973 | – | setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), |
974 | – | &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), |
975 | – | &localToGlobalCutoffGroupIndex[0], &isError); |
976 | – | |
977 | – | if (isError) { |
978 | – | sprintf(painCave.errMsg, |
979 | – | "mpiRefresh errror: fortran didn't like something we gave it.\n"); |
980 | – | painCave.isFatal = 1; |
981 | – | simError(); |
982 | – | } |
983 | – | |
984 | – | sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
985 | – | errorCheckPoint(); |
986 | – | |
987 | – | #endif |
988 | – | } |
989 | – | |
990 | – | |
991 | – | void SimInfo::setupSwitchingFunction() { |
992 | – | int ft = CUBIC; |
993 | – | |
994 | – | if (simParams_->haveSwitchingFunctionType()) { |
995 | – | string funcType = simParams_->getSwitchingFunctionType(); |
996 | – | toUpper(funcType); |
997 | – | if (funcType == "CUBIC") { |
998 | – | ft = CUBIC; |
999 | – | } else { |
1000 | – | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
1001 | – | ft = FIFTH_ORDER_POLY; |
1002 | – | } else { |
1003 | – | // throw error |
1004 | – | sprintf( painCave.errMsg, |
1005 | – | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n" |
1006 | – | "\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", |
1007 | – | funcType.c_str() ); |
1008 | – | painCave.isFatal = 1; |
1009 | – | simError(); |
1010 | – | } |
1011 | – | } |
1012 | – | } |
1013 | – | |
1014 | – | // send switching function notification to switcheroo |
1015 | – | setFunctionType(&ft); |
1016 | – | |
1017 | – | } |
1018 | – | |
1019 | – | void SimInfo::setupAccumulateBoxDipole() { |
1020 | – | |
1021 | – | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1022 | – | if ( simParams_->haveAccumulateBoxDipole() ) |
1023 | – | if ( simParams_->getAccumulateBoxDipole() ) { |
1024 | – | calcBoxDipole_ = true; |
1025 | – | } |
1026 | – | |
1027 | – | } |
1028 | – | |
933 | void SimInfo::addProperty(GenericData* genData) { | |
934 | properties_.addProperty(genData); | |
935 | } | |
# | Line 1060 | Line 964 | namespace OpenMD { | |
964 | Molecule* mol; | |
965 | RigidBody* rb; | |
966 | Atom* atom; | |
967 | + | CutoffGroup* cg; |
968 | SimInfo::MoleculeIterator mi; | |
969 | Molecule::RigidBodyIterator rbIter; | |
970 | < | Molecule::AtomIterator atomIter;; |
970 | > | Molecule::AtomIterator atomIter; |
971 | > | Molecule::CutoffGroupIterator cgIter; |
972 | ||
973 | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
974 | ||
# | Line 1072 | Line 978 | namespace OpenMD { | |
978 | ||
979 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
980 | rb->setSnapshotManager(sman_); | |
981 | + | } |
982 | + | |
983 | + | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { |
984 | + | cg->setSnapshotManager(sman_); |
985 | } | |
986 | } | |
987 | ||
# | Line 1301 | Line 1211 | namespace OpenMD { | |
1211 | ||
1212 | det = intTensor.determinant(); | |
1213 | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | |
1214 | < | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1214 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(det); |
1215 | return; | |
1216 | } | |
1217 | ||
# | Line 1317 | Line 1227 | namespace OpenMD { | |
1227 | ||
1228 | detI = intTensor.determinant(); | |
1229 | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | |
1230 | < | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1230 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(detI); |
1231 | return; | |
1232 | } | |
1233 | /* |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |