# | 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/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" | |
64 | – | |
64 | #ifdef IS_MPI | |
65 | < | #include "UseTheForce/mpiComponentPlan.h" |
66 | < | #include "UseTheForce/DarkSide/simParallel_interface.h" |
68 | < | #endif |
65 | > | #include <mpi.h> |
66 | > | #endif |
67 | ||
68 | using namespace std; | |
69 | namespace OpenMD { | |
# | Line 74 | 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 131 | 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 226 | 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 247 | Line 250 | namespace OpenMD { | |
250 | ndf_local += 3; | |
251 | } | |
252 | } | |
250 | – | |
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 | + | ndfLocal_ = ndf_local; |
263 | + | |
264 | // n_constraints is local, so subtract them on each processor | |
265 | ndf_local -= nConstraints_; | |
266 | ||
267 | #ifdef IS_MPI | |
268 | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | |
269 | + | MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
270 | #else | |
271 | ndf_ = ndf_local; | |
272 | + | nGlobalFluctuatingCharges_ = nfq_local; |
273 | #endif | |
274 | ||
275 | // nZconstraints_ is global, as are the 3 COM translations for the | |
# | Line 274 | Line 286 | namespace OpenMD { | |
286 | #endif | |
287 | return fdf_; | |
288 | } | |
289 | + | |
290 | + | unsigned int SimInfo::getNLocalCutoffGroups(){ |
291 | + | int nLocalCutoffAtoms = 0; |
292 | + | Molecule* mol; |
293 | + | MoleculeIterator mi; |
294 | + | CutoffGroup* cg; |
295 | + | Molecule::CutoffGroupIterator ci; |
296 | + | |
297 | + | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
298 | + | |
299 | + | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
300 | + | cg = mol->nextCutoffGroup(ci)) { |
301 | + | nLocalCutoffAtoms += cg->getNumAtom(); |
302 | + | |
303 | + | } |
304 | + | } |
305 | + | |
306 | + | return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; |
307 | + | } |
308 | ||
309 | void SimInfo::calcNdfRaw() { | |
310 | int ndfRaw_local; | |
# | Line 680 | Line 711 | namespace OpenMD { | |
711 | Atom* atom; | |
712 | set<AtomType*> atomTypes; | |
713 | ||
714 | < | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
715 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
714 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
715 | > | for(atom = mol->beginAtom(ai); atom != NULL; |
716 | > | atom = mol->nextAtom(ai)) { |
717 | atomTypes.insert(atom->getAtomType()); | |
718 | } | |
719 | } | |
720 | < | |
720 | > | |
721 | #ifdef IS_MPI | |
722 | ||
723 | // loop over the found atom types on this processor, and add their | |
724 | // numerical idents to a vector: | |
725 | < | |
725 | > | |
726 | vector<int> foundTypes; | |
727 | set<AtomType*>::iterator i; | |
728 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) | |
# | Line 699 | Line 731 | namespace OpenMD { | |
731 | // count_local holds the number of found types on this processor | |
732 | int count_local = foundTypes.size(); | |
733 | ||
702 | – | // 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); |
706 | – | |
707 | – | // create a vector to hold the globally found types, and resize it: |
708 | – | vector<int> ftGlobal; |
709 | – | ftGlobal.resize(count); |
710 | – | vector<int> counts; |
711 | – | |
734 | int nproc = MPI::COMM_WORLD.Get_size(); | |
735 | < | counts.resize(nproc); |
736 | < | vector<int> disps; |
737 | < | disps.resize(nproc); |
735 | > | |
736 | > | // we need arrays to hold the counts and displacement vectors for |
737 | > | // all processors |
738 | > | vector<int> counts(nproc, 0); |
739 | > | vector<int> disps(nproc, 0); |
740 | ||
741 | < | // now spray out the foundTypes to all the other processors: |
741 | > | // fill the counts array |
742 | > | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], |
743 | > | 1, MPI::INT); |
744 | > | |
745 | > | // use the processor counts to compute the displacement array |
746 | > | disps[0] = 0; |
747 | > | int totalCount = counts[0]; |
748 | > | for (int iproc = 1; iproc < nproc; iproc++) { |
749 | > | disps[iproc] = disps[iproc-1] + counts[iproc-1]; |
750 | > | totalCount += counts[iproc]; |
751 | > | } |
752 | > | |
753 | > | // we need a (possibly redundant) set of all found types: |
754 | > | vector<int> ftGlobal(totalCount); |
755 | ||
756 | + | // now spray out the foundTypes to all the other processors: |
757 | MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, | |
758 | < | &ftGlobal[0], &counts[0], &disps[0], MPI::INT); |
758 | > | &ftGlobal[0], &counts[0], &disps[0], |
759 | > | MPI::INT); |
760 | ||
761 | + | vector<int>::iterator j; |
762 | + | |
763 | // foundIdents is a stl set, so inserting an already found ident | |
764 | // will have no effect. | |
765 | set<int> foundIdents; | |
766 | < | vector<int>::iterator j; |
766 | > | |
767 | for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) | |
768 | foundIdents.insert((*j)); | |
769 | ||
770 | // now iterate over the foundIdents and get the actual atom types | |
771 | // that correspond to these: | |
772 | set<int>::iterator it; | |
773 | < | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
773 | > | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
774 | atomTypes.insert( forceField_->getAtomType((*it)) ); | |
775 | ||
776 | #endif | |
777 | < | |
777 | > | |
778 | return atomTypes; | |
779 | } | |
780 | ||
# | Line 745 | Line 786 | namespace OpenMD { | |
786 | if ( simParams_->getAccumulateBoxDipole() ) { | |
787 | calcBoxDipole_ = true; | |
788 | } | |
789 | < | |
789 | > | |
790 | set<AtomType*>::iterator i; | |
791 | set<AtomType*> atomTypes; | |
792 | atomTypes = getSimulatedAtomTypes(); | |
793 | < | int usesElectrostatic = 0; |
794 | < | int usesMetallic = 0; |
795 | < | int usesDirectional = 0; |
793 | > | bool usesElectrostatic = false; |
794 | > | bool usesMetallic = false; |
795 | > | bool usesDirectional = false; |
796 | > | bool usesFluctuatingCharges = false; |
797 | //loop over all of the atom types | |
798 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
799 | usesElectrostatic |= (*i)->isElectrostatic(); | |
800 | usesMetallic |= (*i)->isMetal(); | |
801 | usesDirectional |= (*i)->isDirectional(); | |
802 | + | usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); |
803 | } | |
804 | ||
805 | < | #ifdef IS_MPI |
806 | < | int temp; |
805 | > | #ifdef IS_MPI |
806 | > | bool temp; |
807 | temp = usesDirectional; | |
808 | < | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
809 | < | |
808 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, |
809 | > | MPI::LOR); |
810 | > | |
811 | temp = usesMetallic; | |
812 | < | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
813 | < | |
812 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, |
813 | > | MPI::LOR); |
814 | > | |
815 | temp = usesElectrostatic; | |
816 | < | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
816 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, |
817 | > | MPI::LOR); |
818 | > | |
819 | > | temp = usesFluctuatingCharges; |
820 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, |
821 | > | MPI::LOR); |
822 | > | #else |
823 | > | |
824 | > | usesDirectionalAtoms_ = usesDirectional; |
825 | > | usesMetallicAtoms_ = usesMetallic; |
826 | > | usesElectrostaticAtoms_ = usesElectrostatic; |
827 | > | usesFluctuatingCharges_ = usesFluctuatingCharges; |
828 | > | |
829 | #endif | |
830 | < | fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_; |
831 | < | fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; |
832 | < | fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; |
833 | < | fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; |
777 | < | fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; |
778 | < | fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; |
830 | > | |
831 | > | requiresPrepair_ = usesMetallicAtoms_ ? true : false; |
832 | > | requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; |
833 | > | requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; |
834 | } | |
835 | ||
836 | < | void SimInfo::setupFortran() { |
837 | < | int isError; |
838 | < | int nExclude, nOneTwo, nOneThree, nOneFour; |
839 | < | vector<int> fortranGlobalGroupMembership; |
836 | > | |
837 | > | vector<int> SimInfo::getGlobalAtomIndices() { |
838 | > | SimInfo::MoleculeIterator mi; |
839 | > | Molecule* mol; |
840 | > | Molecule::AtomIterator ai; |
841 | > | Atom* atom; |
842 | > | |
843 | > | vector<int> GlobalAtomIndices(getNAtoms(), 0); |
844 | ||
845 | < | isError = 0; |
845 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
846 | > | |
847 | > | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
848 | > | GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); |
849 | > | } |
850 | > | } |
851 | > | return GlobalAtomIndices; |
852 | > | } |
853 | ||
854 | < | //globalGroupMembership_ is filled by SimCreator |
855 | < | for (int i = 0; i < nGlobalAtoms_; i++) { |
856 | < | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
854 | > | |
855 | > | vector<int> SimInfo::getGlobalGroupIndices() { |
856 | > | SimInfo::MoleculeIterator mi; |
857 | > | Molecule* mol; |
858 | > | Molecule::CutoffGroupIterator ci; |
859 | > | CutoffGroup* cg; |
860 | > | |
861 | > | vector<int> GlobalGroupIndices; |
862 | > | |
863 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
864 | > | |
865 | > | //local index of cutoff group is trivial, it only depends on the |
866 | > | //order of travesing |
867 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
868 | > | cg = mol->nextCutoffGroup(ci)) { |
869 | > | GlobalGroupIndices.push_back(cg->getGlobalIndex()); |
870 | > | } |
871 | } | |
872 | + | return GlobalGroupIndices; |
873 | + | } |
874 | ||
875 | + | |
876 | + | void SimInfo::prepareTopology() { |
877 | + | int nExclude, nOneTwo, nOneThree, nOneFour; |
878 | + | |
879 | //calculate mass ratio of cutoff group | |
794 | – | vector<RealType> mfact; |
880 | SimInfo::MoleculeIterator mi; | |
881 | Molecule* mol; | |
882 | Molecule::CutoffGroupIterator ci; | |
# | Line 800 | Line 885 | namespace OpenMD { | |
885 | Atom* atom; | |
886 | RealType totalMass; | |
887 | ||
888 | < | //to avoid memory reallocation, reserve enough space for mfact |
889 | < | mfact.reserve(getNCutoffGroups()); |
888 | > | /** |
889 | > | * The mass factor is the relative mass of an atom to the total |
890 | > | * mass of the cutoff group it belongs to. By default, all atoms |
891 | > | * are their own cutoff groups, and therefore have mass factors of |
892 | > | * 1. We need some special handling for massless atoms, which |
893 | > | * will be treated as carrying the entire mass of the cutoff |
894 | > | * group. |
895 | > | */ |
896 | > | massFactors_.clear(); |
897 | > | massFactors_.resize(getNAtoms(), 1.0); |
898 | ||
899 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
900 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
900 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
901 | > | cg = mol->nextCutoffGroup(ci)) { |
902 | ||
903 | totalMass = cg->getMass(); | |
904 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
905 | // Check for massless groups - set mfact to 1 if true | |
906 | < | if (totalMass != 0) |
907 | < | mfact.push_back(atom->getMass()/totalMass); |
906 | > | if (totalMass != 0) |
907 | > | massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; |
908 | else | |
909 | < | mfact.push_back( 1.0 ); |
909 | > | massFactors_[atom->getLocalIndex()] = 1.0; |
910 | } | |
911 | } | |
912 | } | |
913 | ||
914 | < | //fill ident array of local atoms (it is actually ident of |
821 | < | //AtomType, it is so confusing !!!) |
822 | < | vector<int> identArray; |
914 | > | // Build the identArray_ |
915 | ||
916 | < | //to avoid memory reallocation, reserve enough space identArray |
917 | < | identArray.reserve(getNAtoms()); |
826 | < | |
916 | > | identArray_.clear(); |
917 | > | identArray_.reserve(getNAtoms()); |
918 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
919 | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | |
920 | < | identArray.push_back(atom->getIdent()); |
920 | > | identArray_.push_back(atom->getIdent()); |
921 | } | |
922 | } | |
832 | – | |
833 | – | //fill molMembershipArray |
834 | – | //molMembershipArray is filled by SimCreator |
835 | – | vector<int> molMembershipArray(nGlobalAtoms_); |
836 | – | for (int i = 0; i < nGlobalAtoms_; i++) { |
837 | – | molMembershipArray[i] = globalMolMembership_[i] + 1; |
838 | – | } |
923 | ||
924 | < | //setup fortran simulation |
924 | > | //scan topology |
925 | ||
926 | nExclude = excludedInteractions_.getSize(); | |
927 | nOneTwo = oneTwoInteractions_.getSize(); | |
# | Line 849 | Line 933 | namespace OpenMD { | |
933 | int* oneThreeList = oneThreeInteractions_.getPairList(); | |
934 | int* oneFourList = oneFourInteractions_.getPairList(); | |
935 | ||
936 | < | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], |
853 | < | &nExclude, excludeList, |
854 | < | &nOneTwo, oneTwoList, |
855 | < | &nOneThree, oneThreeList, |
856 | < | &nOneFour, oneFourList, |
857 | < | &molMembershipArray[0], &mfact[0], &nCutoffGroups_, |
858 | < | &fortranGlobalGroupMembership[0], &isError); |
859 | < | |
860 | < | if( isError ){ |
861 | < | |
862 | < | sprintf( painCave.errMsg, |
863 | < | "There was an error setting the simulation information in fortran.\n" ); |
864 | < | painCave.isFatal = 1; |
865 | < | painCave.severity = OPENMD_ERROR; |
866 | < | simError(); |
867 | < | } |
868 | < | |
869 | < | |
870 | < | sprintf( checkPointMsg, |
871 | < | "succesfully sent the simulation information to fortran.\n"); |
872 | < | |
873 | < | errorCheckPoint(); |
874 | < | |
875 | < | // Setup number of neighbors in neighbor list if present |
876 | < | if (simParams_->haveNeighborListNeighbors()) { |
877 | < | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
878 | < | setNeighbors(&nlistNeighbors); |
879 | < | } |
880 | < | |
881 | < | #ifdef IS_MPI |
882 | < | //SimInfo is responsible for creating localToGlobalAtomIndex and |
883 | < | //localToGlobalGroupIndex |
884 | < | vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
885 | < | vector<int> localToGlobalCutoffGroupIndex; |
886 | < | mpiSimData parallelData; |
887 | < | |
888 | < | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
889 | < | |
890 | < | //local index(index in DataStorge) of atom is important |
891 | < | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
892 | < | localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; |
893 | < | } |
894 | < | |
895 | < | //local index of cutoff group is trivial, it only depends on the order of travesing |
896 | < | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
897 | < | localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); |
898 | < | } |
899 | < | |
900 | < | } |
901 | < | |
902 | < | //fill up mpiSimData struct |
903 | < | parallelData.nMolGlobal = getNGlobalMolecules(); |
904 | < | parallelData.nMolLocal = getNMolecules(); |
905 | < | parallelData.nAtomsGlobal = getNGlobalAtoms(); |
906 | < | parallelData.nAtomsLocal = getNAtoms(); |
907 | < | parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); |
908 | < | parallelData.nGroupsLocal = getNCutoffGroups(); |
909 | < | parallelData.myNode = worldRank; |
910 | < | MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); |
911 | < | |
912 | < | //pass mpiSimData struct and index arrays to fortran |
913 | < | setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), |
914 | < | &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), |
915 | < | &localToGlobalCutoffGroupIndex[0], &isError); |
916 | < | |
917 | < | if (isError) { |
918 | < | sprintf(painCave.errMsg, |
919 | < | "mpiRefresh errror: fortran didn't like something we gave it.\n"); |
920 | < | painCave.isFatal = 1; |
921 | < | simError(); |
922 | < | } |
923 | < | |
924 | < | sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
925 | < | errorCheckPoint(); |
926 | < | #endif |
927 | < | fortranInitialized_ = true; |
936 | > | topologyDone_ = true; |
937 | } | |
938 | ||
939 | void SimInfo::addProperty(GenericData* genData) { | |
# | Line 961 | Line 970 | namespace OpenMD { | |
970 | Molecule* mol; | |
971 | RigidBody* rb; | |
972 | Atom* atom; | |
973 | + | CutoffGroup* cg; |
974 | SimInfo::MoleculeIterator mi; | |
975 | Molecule::RigidBodyIterator rbIter; | |
976 | < | Molecule::AtomIterator atomIter;; |
976 | > | Molecule::AtomIterator atomIter; |
977 | > | Molecule::CutoffGroupIterator cgIter; |
978 | ||
979 | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
980 | ||
# | Line 974 | Line 985 | namespace OpenMD { | |
985 | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | |
986 | rb->setSnapshotManager(sman_); | |
987 | } | |
988 | + | |
989 | + | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { |
990 | + | cg->setSnapshotManager(sman_); |
991 | + | } |
992 | } | |
993 | ||
994 | } | |
995 | ||
981 | – | Vector3d SimInfo::getComVel(){ |
982 | – | SimInfo::MoleculeIterator i; |
983 | – | Molecule* mol; |
996 | ||
985 | – | Vector3d comVel(0.0); |
986 | – | RealType totalMass = 0.0; |
987 | – | |
988 | – | |
989 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
990 | – | RealType mass = mol->getMass(); |
991 | – | totalMass += mass; |
992 | – | comVel += mass * mol->getComVel(); |
993 | – | } |
994 | – | |
995 | – | #ifdef IS_MPI |
996 | – | RealType tmpMass = totalMass; |
997 | – | Vector3d tmpComVel(comVel); |
998 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
999 | – | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1000 | – | #endif |
1001 | – | |
1002 | – | comVel /= totalMass; |
1003 | – | |
1004 | – | return comVel; |
1005 | – | } |
1006 | – | |
1007 | – | Vector3d SimInfo::getCom(){ |
1008 | – | SimInfo::MoleculeIterator i; |
1009 | – | Molecule* mol; |
1010 | – | |
1011 | – | Vector3d com(0.0); |
1012 | – | RealType totalMass = 0.0; |
1013 | – | |
1014 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1015 | – | RealType mass = mol->getMass(); |
1016 | – | totalMass += mass; |
1017 | – | com += mass * mol->getCom(); |
1018 | – | } |
1019 | – | |
1020 | – | #ifdef IS_MPI |
1021 | – | RealType tmpMass = totalMass; |
1022 | – | Vector3d tmpCom(com); |
1023 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1024 | – | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1025 | – | #endif |
1026 | – | |
1027 | – | com /= totalMass; |
1028 | – | |
1029 | – | return com; |
1030 | – | |
1031 | – | } |
1032 | – | |
997 | ostream& operator <<(ostream& o, SimInfo& info) { | |
998 | ||
999 | return o; | |
1000 | } | |
1001 | ||
1002 | < | |
1039 | < | /* |
1040 | < | Returns center of mass and center of mass velocity in one function call. |
1041 | < | */ |
1042 | < | |
1043 | < | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
1044 | < | SimInfo::MoleculeIterator i; |
1045 | < | Molecule* mol; |
1046 | < | |
1047 | < | |
1048 | < | RealType totalMass = 0.0; |
1049 | < | |
1050 | < | |
1051 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1052 | < | RealType mass = mol->getMass(); |
1053 | < | totalMass += mass; |
1054 | < | com += mass * mol->getCom(); |
1055 | < | comVel += mass * mol->getComVel(); |
1056 | < | } |
1057 | < | |
1058 | < | #ifdef IS_MPI |
1059 | < | RealType tmpMass = totalMass; |
1060 | < | Vector3d tmpCom(com); |
1061 | < | Vector3d tmpComVel(comVel); |
1062 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1063 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1064 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1065 | < | #endif |
1066 | < | |
1067 | < | com /= totalMass; |
1068 | < | comVel /= totalMass; |
1069 | < | } |
1070 | < | |
1071 | < | /* |
1072 | < | Return intertia tensor for entire system and angular momentum Vector. |
1073 | < | |
1074 | < | |
1075 | < | [ Ixx -Ixy -Ixz ] |
1076 | < | J =| -Iyx Iyy -Iyz | |
1077 | < | [ -Izx -Iyz Izz ] |
1078 | < | */ |
1079 | < | |
1080 | < | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1081 | < | |
1082 | < | |
1083 | < | RealType xx = 0.0; |
1084 | < | RealType yy = 0.0; |
1085 | < | RealType zz = 0.0; |
1086 | < | RealType xy = 0.0; |
1087 | < | RealType xz = 0.0; |
1088 | < | RealType yz = 0.0; |
1089 | < | Vector3d com(0.0); |
1090 | < | Vector3d comVel(0.0); |
1091 | < | |
1092 | < | getComAll(com, comVel); |
1093 | < | |
1094 | < | SimInfo::MoleculeIterator i; |
1095 | < | Molecule* mol; |
1096 | < | |
1097 | < | Vector3d thisq(0.0); |
1098 | < | Vector3d thisv(0.0); |
1099 | < | |
1100 | < | RealType thisMass = 0.0; |
1101 | < | |
1102 | < | |
1103 | < | |
1104 | < | |
1105 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1106 | < | |
1107 | < | thisq = mol->getCom()-com; |
1108 | < | thisv = mol->getComVel()-comVel; |
1109 | < | thisMass = mol->getMass(); |
1110 | < | // Compute moment of intertia coefficients. |
1111 | < | xx += thisq[0]*thisq[0]*thisMass; |
1112 | < | yy += thisq[1]*thisq[1]*thisMass; |
1113 | < | zz += thisq[2]*thisq[2]*thisMass; |
1114 | < | |
1115 | < | // compute products of intertia |
1116 | < | xy += thisq[0]*thisq[1]*thisMass; |
1117 | < | xz += thisq[0]*thisq[2]*thisMass; |
1118 | < | yz += thisq[1]*thisq[2]*thisMass; |
1119 | < | |
1120 | < | angularMomentum += cross( thisq, thisv ) * thisMass; |
1121 | < | |
1122 | < | } |
1123 | < | |
1124 | < | |
1125 | < | inertiaTensor(0,0) = yy + zz; |
1126 | < | inertiaTensor(0,1) = -xy; |
1127 | < | inertiaTensor(0,2) = -xz; |
1128 | < | inertiaTensor(1,0) = -xy; |
1129 | < | inertiaTensor(1,1) = xx + zz; |
1130 | < | inertiaTensor(1,2) = -yz; |
1131 | < | inertiaTensor(2,0) = -xz; |
1132 | < | inertiaTensor(2,1) = -yz; |
1133 | < | inertiaTensor(2,2) = xx + yy; |
1134 | < | |
1135 | < | #ifdef IS_MPI |
1136 | < | Mat3x3d tmpI(inertiaTensor); |
1137 | < | Vector3d tmpAngMom; |
1138 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1139 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1140 | < | #endif |
1141 | < | |
1142 | < | return; |
1143 | < | } |
1144 | < | |
1145 | < | //Returns the angular momentum of the system |
1146 | < | Vector3d SimInfo::getAngularMomentum(){ |
1147 | < | |
1148 | < | Vector3d com(0.0); |
1149 | < | Vector3d comVel(0.0); |
1150 | < | Vector3d angularMomentum(0.0); |
1151 | < | |
1152 | < | getComAll(com,comVel); |
1153 | < | |
1154 | < | SimInfo::MoleculeIterator i; |
1155 | < | Molecule* mol; |
1156 | < | |
1157 | < | Vector3d thisr(0.0); |
1158 | < | Vector3d thisp(0.0); |
1159 | < | |
1160 | < | RealType thisMass; |
1161 | < | |
1162 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1163 | < | thisMass = mol->getMass(); |
1164 | < | thisr = mol->getCom()-com; |
1165 | < | thisp = (mol->getComVel()-comVel)*thisMass; |
1166 | < | |
1167 | < | angularMomentum += cross( thisr, thisp ); |
1168 | < | |
1169 | < | } |
1170 | < | |
1171 | < | #ifdef IS_MPI |
1172 | < | Vector3d tmpAngMom; |
1173 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1174 | < | #endif |
1175 | < | |
1176 | < | return angularMomentum; |
1177 | < | } |
1178 | < | |
1002 | > | |
1003 | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | |
1004 | return IOIndexToIntegrableObject.at(index); | |
1005 | } | |
# | Line 1183 | Line 1007 | namespace OpenMD { | |
1007 | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { | |
1008 | IOIndexToIntegrableObject= v; | |
1009 | } | |
1186 | – | |
1187 | – | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1188 | – | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1189 | – | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1190 | – | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1191 | – | */ |
1192 | – | void SimInfo::getGyrationalVolume(RealType &volume){ |
1193 | – | Mat3x3d intTensor; |
1194 | – | RealType det; |
1195 | – | Vector3d dummyAngMom; |
1196 | – | RealType sysconstants; |
1197 | – | RealType geomCnst; |
1198 | – | |
1199 | – | geomCnst = 3.0/2.0; |
1200 | – | /* Get the inertial tensor and angular momentum for free*/ |
1201 | – | getInertiaTensor(intTensor,dummyAngMom); |
1202 | – | |
1203 | – | det = intTensor.determinant(); |
1204 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1205 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1206 | – | return; |
1207 | – | } |
1208 | – | |
1209 | – | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1210 | – | Mat3x3d intTensor; |
1211 | – | Vector3d dummyAngMom; |
1212 | – | RealType sysconstants; |
1213 | – | RealType geomCnst; |
1214 | – | |
1215 | – | geomCnst = 3.0/2.0; |
1216 | – | /* Get the inertial tensor and angular momentum for free*/ |
1217 | – | getInertiaTensor(intTensor,dummyAngMom); |
1218 | – | |
1219 | – | detI = intTensor.determinant(); |
1220 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1221 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1222 | – | return; |
1223 | – | } |
1010 | /* | |
1011 | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { | |
1012 | assert( v.size() == nAtoms_ + nRigidBodies_); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |