# | Line 60 | Line 60 | |
---|---|---|
60 | #include "io/ForceFieldOptions.hpp" | |
61 | #include "UseTheForce/ForceField.hpp" | |
62 | #include "nonbonded/SwitchingFunction.hpp" | |
63 | + | #ifdef IS_MPI |
64 | + | #include <mpi.h> |
65 | + | #endif |
66 | ||
67 | using namespace std; | |
68 | namespace OpenMD { | |
# | Line 71 | 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 125 | 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 | |
128 | – | std::cerr << "nGA = " << nGlobalAtoms_ << "\n"; |
129 | – | std::cerr << "nCA = " << nCutoffAtoms << "\n"; |
130 | – | std::cerr << "nG = " << nGroups << "\n"; |
131 | ||
132 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
133 | – | |
134 | – | std::cerr << "nGCG = " << nGlobalCutoffGroups_ << "\n"; |
133 | ||
134 | //every free atom (atom does not belong to rigid bodies) is an | |
135 | //integrable object therefore the total number of integrable objects | |
# | Line 273 | Line 271 | namespace OpenMD { | |
271 | fdf_ = fdf_local; | |
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() { | |
# | Line 680 | 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 | < | |
706 | > | |
707 | #ifdef IS_MPI | |
708 | ||
709 | // loop over the found atom types on this processor, and add their | |
710 | // numerical idents to a vector: | |
711 | < | |
711 | > | |
712 | vector<int> foundTypes; | |
713 | set<AtomType*>::iterator i; | |
714 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) | |
# | Line 699 | Line 717 | namespace OpenMD { | |
717 | // count_local holds the number of found types on this processor | |
718 | int count_local = foundTypes.size(); | |
719 | ||
720 | < | // count holds the total number of found types on all processors |
703 | < | // (some will be redundant with the ones found locally): |
704 | < | int count; |
705 | < | MPI::COMM_WORLD.Allreduce(&count_local, &count, 1, MPI::INT, MPI::SUM); |
720 | > | int nproc = MPI::COMM_WORLD.Get_size(); |
721 | ||
722 | < | // create a vector to hold the globally found types, and resize it: |
723 | < | vector<int> ftGlobal; |
724 | < | ftGlobal.resize(count); |
725 | < | vector<int> counts; |
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 | < | int nproc = MPI::COMM_WORLD.Get_size(); |
728 | < | counts.resize(nproc); |
729 | < | vector<int> disps; |
730 | < | disps.resize(nproc); |
727 | > | // fill the counts array |
728 | > | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], |
729 | > | 1, MPI::INT); |
730 | > | |
731 | > | // use the processor counts to compute the displacement array |
732 | > | disps[0] = 0; |
733 | > | int totalCount = counts[0]; |
734 | > | for (int iproc = 1; iproc < nproc; iproc++) { |
735 | > | disps[iproc] = disps[iproc-1] + counts[iproc-1]; |
736 | > | totalCount += counts[iproc]; |
737 | > | } |
738 | ||
739 | < | // now spray out the foundTypes to all the other processors: |
739 | > | // we need a (possibly redundant) set of all found types: |
740 | > | vector<int> ftGlobal(totalCount); |
741 | ||
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], MPI::INT); |
744 | > | &ftGlobal[0], &counts[0], &disps[0], |
745 | > | MPI::INT); |
746 | ||
747 | + | vector<int>::iterator j; |
748 | + | |
749 | // foundIdents is a stl set, so inserting an already found ident | |
750 | // will have no effect. | |
751 | set<int> foundIdents; | |
752 | < | vector<int>::iterator j; |
752 | > | |
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) |
759 | > | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
760 | atomTypes.insert( forceField_->getAtomType((*it)) ); | |
761 | ||
762 | #endif | |
763 | < | |
763 | > | |
764 | return atomTypes; | |
765 | } | |
766 | ||
# | Line 745 | Line 772 | namespace OpenMD { | |
772 | if ( simParams_->getAccumulateBoxDipole() ) { | |
773 | calcBoxDipole_ = true; | |
774 | } | |
775 | < | |
775 | > | |
776 | set<AtomType*>::iterator i; | |
777 | set<AtomType*> atomTypes; | |
778 | atomTypes = getSimulatedAtomTypes(); | |
# | Line 758 | 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_; |
777 | < | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; |
778 | < | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; |
806 | > | |
807 | > | requiresPrepair_ = usesMetallicAtoms_ ? true : false; |
808 | > | requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; |
809 | > | requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; |
810 | } | |
811 | ||
812 | ||
# | Line 818 | Line 849 | namespace OpenMD { | |
849 | } | |
850 | ||
851 | ||
852 | < | void SimInfo::setupFortran() { |
822 | < | int isError; |
852 | > | void SimInfo::prepareTopology() { |
853 | int nExclude, nOneTwo, nOneThree, nOneFour; | |
824 | – | vector<int> fortranGlobalGroupMembership; |
825 | – | |
826 | – | isError = 0; |
854 | ||
828 | – | //globalGroupMembership_ is filled by SimCreator |
829 | – | for (int i = 0; i < nGlobalAtoms_; i++) { |
830 | – | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
831 | – | } |
832 | – | |
855 | //calculate mass ratio of cutoff group | |
834 | – | vector<RealType> mfact; |
856 | SimInfo::MoleculeIterator mi; | |
857 | Molecule* mol; | |
858 | Molecule::CutoffGroupIterator ci; | |
# | Line 840 | 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 | } | |
# | Line 866 | Line 896 | namespace OpenMD { | |
896 | identArray_.push_back(atom->getIdent()); | |
897 | } | |
898 | } | |
869 | – | |
870 | – | //fill molMembershipArray |
871 | – | //molMembershipArray is filled by SimCreator |
872 | – | vector<int> molMembershipArray(nGlobalAtoms_); |
873 | – | for (int i = 0; i < nGlobalAtoms_; i++) { |
874 | – | molMembershipArray[i] = globalMolMembership_[i] + 1; |
875 | – | } |
899 | ||
900 | < | //setup fortran simulation |
900 | > | //scan topology |
901 | ||
902 | nExclude = excludedInteractions_.getSize(); | |
903 | nOneTwo = oneTwoInteractions_.getSize(); | |
# | Line 886 | Line 909 | namespace OpenMD { | |
909 | int* oneThreeList = oneThreeInteractions_.getPairList(); | |
910 | int* oneFourList = oneFourInteractions_.getPairList(); | |
911 | ||
912 | < | //setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray_[0], |
890 | < | // &nExclude, excludeList, |
891 | < | // &nOneTwo, oneTwoList, |
892 | < | // &nOneThree, oneThreeList, |
893 | < | // &nOneFour, oneFourList, |
894 | < | // &molMembershipArray[0], &mfact[0], &nCutoffGroups_, |
895 | < | // &fortranGlobalGroupMembership[0], &isError); |
896 | < | |
897 | < | // if( isError ){ |
898 | < | // |
899 | < | // sprintf( painCave.errMsg, |
900 | < | // "There was an error setting the simulation information in fortran.\n" ); |
901 | < | // painCave.isFatal = 1; |
902 | < | // painCave.severity = OPENMD_ERROR; |
903 | < | // simError(); |
904 | < | //} |
905 | < | |
906 | < | |
907 | < | // sprintf( checkPointMsg, |
908 | < | // "succesfully sent the simulation information to fortran.\n"); |
909 | < | |
910 | < | // errorCheckPoint(); |
911 | < | |
912 | < | // Setup number of neighbors in neighbor list if present |
913 | < | //if (simParams_->haveNeighborListNeighbors()) { |
914 | < | // int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
915 | < | // setNeighbors(&nlistNeighbors); |
916 | < | //} |
917 | < | |
918 | < | #ifdef IS_MPI |
919 | < | // mpiSimData parallelData; |
920 | < | |
921 | < | //fill up mpiSimData struct |
922 | < | // parallelData.nMolGlobal = getNGlobalMolecules(); |
923 | < | // parallelData.nMolLocal = getNMolecules(); |
924 | < | // parallelData.nAtomsGlobal = getNGlobalAtoms(); |
925 | < | // parallelData.nAtomsLocal = getNAtoms(); |
926 | < | // parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); |
927 | < | // parallelData.nGroupsLocal = getNCutoffGroups(); |
928 | < | // parallelData.myNode = worldRank; |
929 | < | // MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); |
930 | < | |
931 | < | //pass mpiSimData struct and index arrays to fortran |
932 | < | //setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), |
933 | < | // &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), |
934 | < | // &localToGlobalCutoffGroupIndex[0], &isError); |
935 | < | |
936 | < | // if (isError) { |
937 | < | // sprintf(painCave.errMsg, |
938 | < | // "mpiRefresh errror: fortran didn't like something we gave it.\n"); |
939 | < | // painCave.isFatal = 1; |
940 | < | // simError(); |
941 | < | // } |
942 | < | |
943 | < | // sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
944 | < | // errorCheckPoint(); |
945 | < | #endif |
946 | < | |
947 | < | // initFortranFF(&isError); |
948 | < | // if (isError) { |
949 | < | // sprintf(painCave.errMsg, |
950 | < | // "initFortranFF errror: fortran didn't like something we gave it.\n"); |
951 | < | // painCave.isFatal = 1; |
952 | < | // simError(); |
953 | < | // } |
954 | < | // fortranInitialized_ = true; |
912 | > | topologyDone_ = true; |
913 | } | |
914 | ||
915 | void SimInfo::addProperty(GenericData* genData) { |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |