# | Line 54 | Line 54 | |
---|---|---|
54 | #include "math/Vector3.hpp" | |
55 | #include "primitives/Molecule.hpp" | |
56 | #include "primitives/StuntDouble.hpp" | |
57 | – | #include "UseTheForce/doForces_interface.h" |
58 | – | #include "UseTheForce/DarkSide/neighborLists_interface.h" |
57 | #include "utils/MemoryUtils.hpp" | |
58 | #include "utils/simError.h" | |
59 | #include "selection/SelectionManager.hpp" | |
60 | #include "io/ForceFieldOptions.hpp" | |
61 | #include "UseTheForce/ForceField.hpp" | |
62 | #include "nonbonded/SwitchingFunction.hpp" | |
65 | – | |
66 | – | |
63 | #ifdef IS_MPI | |
64 | < | #include "UseTheForce/mpiComponentPlan.h" |
65 | < | #include "UseTheForce/DarkSide/simParallel_interface.h" |
70 | < | #endif |
64 | > | #include <mpi.h> |
65 | > | #endif |
66 | ||
67 | using namespace std; | |
68 | namespace OpenMD { | |
# | Line 79 | Line 74 | namespace OpenMD { | |
74 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
75 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), | |
76 | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | |
77 | < | nConstraints_(0), sman_(NULL), fortranInitialized_(false), |
77 | > | nConstraints_(0), sman_(NULL), topologyDone_(false), |
78 | calcBoxDipole_(false), useAtomicVirial_(true) { | |
79 | ||
80 | MoleculeStamp* molStamp; | |
# | Line 133 | Line 128 | namespace OpenMD { | |
128 | //equal to the total number of atoms minus number of atoms belong to | |
129 | //cutoff group defined in meta-data file plus the number of cutoff | |
130 | //groups defined in meta-data file | |
131 | + | |
132 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
133 | ||
134 | //every free atom (atom does not belong to rigid bodies) is an | |
# | Line 276 | Line 272 | namespace OpenMD { | |
272 | #endif | |
273 | return fdf_; | |
274 | } | |
275 | + | |
276 | + | unsigned int SimInfo::getNLocalCutoffGroups(){ |
277 | + | int nLocalCutoffAtoms = 0; |
278 | + | Molecule* mol; |
279 | + | MoleculeIterator mi; |
280 | + | CutoffGroup* cg; |
281 | + | Molecule::CutoffGroupIterator ci; |
282 | ||
283 | + | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
284 | + | |
285 | + | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
286 | + | cg = mol->nextCutoffGroup(ci)) { |
287 | + | nLocalCutoffAtoms += cg->getNumAtom(); |
288 | + | |
289 | + | } |
290 | + | } |
291 | + | |
292 | + | return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; |
293 | + | } |
294 | + | |
295 | void SimInfo::calcNdfRaw() { | |
296 | int ndfRaw_local; | |
297 | ||
# | Line 657 | Line 672 | namespace OpenMD { | |
672 | /** | |
673 | * update | |
674 | * | |
675 | < | * Performs the global checks and variable settings after the objects have been |
676 | < | * created. |
675 | > | * Performs the global checks and variable settings after the |
676 | > | * objects have been created. |
677 | * | |
678 | */ | |
679 | < | void SimInfo::update() { |
665 | < | |
679 | > | void SimInfo::update() { |
680 | setupSimVariables(); | |
667 | – | setupCutoffs(); |
668 | – | setupSwitching(); |
669 | – | setupElectrostatics(); |
670 | – | setupNeighborlists(); |
671 | – | |
672 | – | #ifdef IS_MPI |
673 | – | setupFortranParallel(); |
674 | – | #endif |
675 | – | setupFortranSim(); |
676 | – | fortranInitialized_ = true; |
677 | – | |
681 | calcNdf(); | |
682 | calcNdfRaw(); | |
683 | calcNdfTrans(); | |
684 | } | |
685 | ||
686 | + | /** |
687 | + | * getSimulatedAtomTypes |
688 | + | * |
689 | + | * Returns an STL set of AtomType* that are actually present in this |
690 | + | * simulation. Must query all processors to assemble this information. |
691 | + | * |
692 | + | */ |
693 | set<AtomType*> SimInfo::getSimulatedAtomTypes() { | |
694 | SimInfo::MoleculeIterator mi; | |
695 | Molecule* mol; | |
# | Line 687 | Line 697 | namespace OpenMD { | |
697 | Atom* atom; | |
698 | set<AtomType*> atomTypes; | |
699 | ||
700 | < | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
701 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
700 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
701 | > | for(atom = mol->beginAtom(ai); atom != NULL; |
702 | > | atom = mol->nextAtom(ai)) { |
703 | atomTypes.insert(atom->getAtomType()); | |
704 | } | |
705 | } | |
706 | < | return atomTypes; |
707 | < | } |
706 | > | |
707 | > | #ifdef IS_MPI |
708 | ||
709 | < | /** |
710 | < | * 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() { |
709 | > | // loop over the found atom types on this processor, and add their |
710 | > | // numerical idents to a vector: |
711 | ||
712 | < | if (simParams_->haveCutoffRadius()) { |
713 | < | cutoffRadius_ = simParams_->getCutoffRadius(); |
714 | < | } else { |
715 | < | 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 | < | } |
712 | > | vector<int> foundTypes; |
713 | > | set<AtomType*>::iterator i; |
714 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) |
715 | > | foundTypes.push_back( (*i)->getIdent() ); |
716 | ||
717 | < | InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); |
717 | > | // count_local holds the number of found types on this processor |
718 | > | int count_local = foundTypes.size(); |
719 | ||
720 | < | map<string, CutoffMethod> stringToCutoffMethod; |
721 | < | stringToCutoffMethod["HARD"] = HARD; |
722 | < | stringToCutoffMethod["SWITCHING_FUNCTION"] = SWITCHING_FUNCTION; |
723 | < | stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL; |
724 | < | stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; |
720 | > | int nproc = MPI::COMM_WORLD.Get_size(); |
721 | > | |
722 | > | // we need arrays to hold the counts and displacement vectors for |
723 | > | // all processors |
724 | > | vector<int> counts(nproc, 0); |
725 | > | vector<int> disps(nproc, 0); |
726 | > | |
727 | > | // fill the counts array |
728 | > | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], |
729 | > | 1, MPI::INT); |
730 | ||
731 | < | if (simParams_->haveCutoffMethod()) { |
732 | < | string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); |
733 | < | map<string, CutoffMethod>::iterator i; |
734 | < | i = stringToCutoffMethod.find(cutMeth); |
735 | < | if (i == stringToCutoffMethod.end()) { |
736 | < | 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; |
731 | > | // use the processor counts to compute the displacement array |
732 | > | disps[0] = 0; |
733 | > | int totalCount = counts[0]; |
734 | > | for (int iproc = 1; iproc < nproc; iproc++) { |
735 | > | disps[iproc] = disps[iproc-1] + counts[iproc-1]; |
736 | > | totalCount += counts[iproc]; |
737 | } | |
738 | ||
739 | < | InteractionManager::Instance()->setCutoffMethod(cutoffMethod_); |
740 | < | } |
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() { |
739 | > | // we need a (possibly redundant) set of all found types: |
740 | > | vector<int> ftGlobal(totalCount); |
741 | ||
742 | < | if (simParams_->haveSwitchingRadius()) { |
743 | < | switchingRadius_ = simParams_->getSwitchingRadius(); |
744 | < | if (switchingRadius_ > cutoffRadius_) { |
745 | < | 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_); |
742 | > | // now spray out the foundTypes to all the other processors: |
743 | > | MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, |
744 | > | &ftGlobal[0], &counts[0], &disps[0], |
745 | > | MPI::INT); |
746 | ||
747 | < | SwitchingFunctionType ft; |
818 | < | |
819 | < | if (simParams_->haveSwitchingFunctionType()) { |
820 | < | string funcType = simParams_->getSwitchingFunctionType(); |
821 | < | toUpper(funcType); |
822 | < | if (funcType == "CUBIC") { |
823 | < | ft = cubic; |
824 | < | } else { |
825 | < | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
826 | < | ft = fifth_order_poly; |
827 | < | } else { |
828 | < | // throw error |
829 | < | sprintf( painCave.errMsg, |
830 | < | "SimInfo : Unknown switchingFunctionType. (Input file specified %s .)\n" |
831 | < | "\tswitchingFunctionType must be one of: " |
832 | < | "\"cubic\" or \"fifth_order_polynomial\".", |
833 | < | funcType.c_str() ); |
834 | < | painCave.isFatal = 1; |
835 | < | painCave.severity = OPENMD_ERROR; |
836 | < | simError(); |
837 | < | } |
838 | < | } |
839 | < | } |
747 | > | vector<int>::iterator j; |
748 | ||
749 | < | InteractionManager::Instance()->setSwitchingFunctionType(ft); |
750 | < | } |
749 | > | // foundIdents is a stl set, so inserting an already found ident |
750 | > | // will have no effect. |
751 | > | set<int> foundIdents; |
752 | ||
753 | < | /** |
754 | < | * setupSkinThickness |
755 | < | * |
756 | < | * If the skinThickness was explicitly set, use that value (but check it) |
757 | < | * If the skinThickness was not explicitly set: use 1.0 angstroms |
758 | < | */ |
759 | < | void SimInfo::setupSkinThickness() { |
760 | < | if (simParams_->haveSkinThickness()) { |
761 | < | skinThickness_ = simParams_->getSkinThickness(); |
762 | < | } else { |
763 | < | skinThickness_ = 1.0; |
764 | < | 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 | < | } |
753 | > | for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) |
754 | > | foundIdents.insert((*j)); |
755 | > | |
756 | > | // now iterate over the foundIdents and get the actual atom types |
757 | > | // that correspond to these: |
758 | > | set<int>::iterator it; |
759 | > | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
760 | > | atomTypes.insert( forceField_->getAtomType((*it)) ); |
761 | > | |
762 | > | #endif |
763 | > | |
764 | > | return atomTypes; |
765 | } | |
766 | ||
767 | < | void SimInfo::setupSimType() { |
767 | > | void SimInfo::setupSimVariables() { |
768 | > | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
769 | > | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
770 | > | calcBoxDipole_ = false; |
771 | > | if ( simParams_->haveAccumulateBoxDipole() ) |
772 | > | if ( simParams_->getAccumulateBoxDipole() ) { |
773 | > | calcBoxDipole_ = true; |
774 | > | } |
775 | > | |
776 | set<AtomType*>::iterator i; | |
777 | set<AtomType*> atomTypes; | |
778 | < | atomTypes = getSimulatedAtomTypes(); |
868 | < | |
869 | < | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
870 | < | |
778 | > | atomTypes = getSimulatedAtomTypes(); |
779 | int usesElectrostatic = 0; | |
780 | int usesMetallic = 0; | |
781 | int usesDirectional = 0; | |
# | Line 877 | Line 785 | namespace OpenMD { | |
785 | usesMetallic |= (*i)->isMetal(); | |
786 | usesDirectional |= (*i)->isDirectional(); | |
787 | } | |
788 | < | |
788 | > | |
789 | #ifdef IS_MPI | |
790 | int temp; | |
791 | temp = usesDirectional; | |
792 | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
793 | < | |
793 | > | |
794 | temp = usesMetallic; | |
795 | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
796 | < | |
796 | > | |
797 | temp = usesElectrostatic; | |
798 | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
799 | + | #else |
800 | + | |
801 | + | usesDirectionalAtoms_ = usesDirectional; |
802 | + | usesMetallicAtoms_ = usesMetallic; |
803 | + | usesElectrostaticAtoms_ = usesElectrostatic; |
804 | + | |
805 | #endif | |
806 | < | fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; |
807 | < | fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; |
808 | < | fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; |
809 | < | fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; |
896 | < | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; |
897 | < | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; |
806 | > | |
807 | > | requiresPrepair_ = usesMetallicAtoms_ ? true : false; |
808 | > | requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; |
809 | > | requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; |
810 | } | |
811 | ||
812 | < | void SimInfo::setupFortranSim() { |
813 | < | int isError; |
814 | < | int nExclude, nOneTwo, nOneThree, nOneFour; |
815 | < | vector<int> fortranGlobalGroupMembership; |
812 | > | |
813 | > | vector<int> SimInfo::getGlobalAtomIndices() { |
814 | > | SimInfo::MoleculeIterator mi; |
815 | > | Molecule* mol; |
816 | > | Molecule::AtomIterator ai; |
817 | > | Atom* atom; |
818 | > | |
819 | > | vector<int> GlobalAtomIndices(getNAtoms(), 0); |
820 | ||
821 | < | notifyFortranSkinThickness(&skinThickness_); |
821 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
822 | > | |
823 | > | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
824 | > | GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); |
825 | > | } |
826 | > | } |
827 | > | return GlobalAtomIndices; |
828 | > | } |
829 | ||
907 | – | int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; |
908 | – | int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; |
909 | – | notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); |
830 | ||
831 | < | isError = 0; |
831 | > | vector<int> SimInfo::getGlobalGroupIndices() { |
832 | > | SimInfo::MoleculeIterator mi; |
833 | > | Molecule* mol; |
834 | > | Molecule::CutoffGroupIterator ci; |
835 | > | CutoffGroup* cg; |
836 | ||
837 | < | //globalGroupMembership_ is filled by SimCreator |
838 | < | for (int i = 0; i < nGlobalAtoms_; i++) { |
839 | < | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
837 | > | vector<int> GlobalGroupIndices; |
838 | > | |
839 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
840 | > | |
841 | > | //local index of cutoff group is trivial, it only depends on the |
842 | > | //order of travesing |
843 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
844 | > | cg = mol->nextCutoffGroup(ci)) { |
845 | > | GlobalGroupIndices.push_back(cg->getGlobalIndex()); |
846 | > | } |
847 | } | |
848 | + | return GlobalGroupIndices; |
849 | + | } |
850 | ||
851 | + | |
852 | + | void SimInfo::prepareTopology() { |
853 | + | int nExclude, nOneTwo, nOneThree, nOneFour; |
854 | + | |
855 | //calculate mass ratio of cutoff group | |
919 | – | vector<RealType> mfact; |
856 | SimInfo::MoleculeIterator mi; | |
857 | Molecule* mol; | |
858 | Molecule::CutoffGroupIterator ci; | |
# | Line 925 | Line 861 | namespace OpenMD { | |
861 | Atom* atom; | |
862 | RealType totalMass; | |
863 | ||
864 | < | //to avoid memory reallocation, reserve enough space for mfact |
865 | < | mfact.reserve(getNCutoffGroups()); |
864 | > | /** |
865 | > | * The mass factor is the relative mass of an atom to the total |
866 | > | * mass of the cutoff group it belongs to. By default, all atoms |
867 | > | * are their own cutoff groups, and therefore have mass factors of |
868 | > | * 1. We need some special handling for massless atoms, which |
869 | > | * will be treated as carrying the entire mass of the cutoff |
870 | > | * group. |
871 | > | */ |
872 | > | massFactors_.clear(); |
873 | > | massFactors_.resize(getNAtoms(), 1.0); |
874 | ||
875 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
876 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
876 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
877 | > | cg = mol->nextCutoffGroup(ci)) { |
878 | ||
879 | totalMass = cg->getMass(); | |
880 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
881 | // Check for massless groups - set mfact to 1 if true | |
882 | < | if (totalMass != 0) |
883 | < | mfact.push_back(atom->getMass()/totalMass); |
882 | > | if (totalMass != 0) |
883 | > | massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; |
884 | else | |
885 | < | mfact.push_back( 1.0 ); |
885 | > | massFactors_[atom->getLocalIndex()] = 1.0; |
886 | } | |
887 | } | |
888 | } | |
889 | ||
890 | < | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) |
946 | < | vector<int> identArray; |
890 | > | // Build the identArray_ |
891 | ||
892 | < | //to avoid memory reallocation, reserve enough space identArray |
893 | < | identArray.reserve(getNAtoms()); |
950 | < | |
892 | > | identArray_.clear(); |
893 | > | identArray_.reserve(getNAtoms()); |
894 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
895 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
896 | < | identArray.push_back(atom->getIdent()); |
896 | > | identArray_.push_back(atom->getIdent()); |
897 | } | |
898 | } | |
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 | – | } |
899 | ||
900 | < | //setup fortran simulation |
900 | > | //scan topology |
901 | ||
902 | nExclude = excludedInteractions_.getSize(); | |
903 | nOneTwo = oneTwoInteractions_.getSize(); | |
# | Line 973 | Line 909 | namespace OpenMD { | |
909 | int* oneThreeList = oneThreeInteractions_.getPairList(); | |
910 | int* oneFourList = oneFourInteractions_.getPairList(); | |
911 | ||
912 | < | 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 | < | |
1006 | < | } |
1007 | < | |
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 |
912 | > | topologyDone_ = true; |
913 | } | |
914 | ||
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 | – | |
915 | void SimInfo::addProperty(GenericData* genData) { | |
916 | properties_.addProperty(genData); | |
917 | } | |
# | Line 1111 | Line 946 | namespace OpenMD { | |
946 | Molecule* mol; | |
947 | RigidBody* rb; | |
948 | Atom* atom; | |
949 | + | CutoffGroup* cg; |
950 | SimInfo::MoleculeIterator mi; | |
951 | Molecule::RigidBodyIterator rbIter; | |
952 | < | Molecule::AtomIterator atomIter;; |
952 | > | Molecule::AtomIterator atomIter; |
953 | > | Molecule::CutoffGroupIterator cgIter; |
954 | ||
955 | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
956 | ||
# | Line 1124 | Line 961 | namespace OpenMD { | |
961 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
962 | rb->setSnapshotManager(sman_); | |
963 | } | |
964 | + | |
965 | + | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { |
966 | + | cg->setSnapshotManager(sman_); |
967 | + | } |
968 | } | |
969 | ||
970 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |