# | 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 58 | Line 59 | |
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 | + | #ifdef IS_MPI |
65 | + | #include <mpi.h> |
66 | + | #endif |
67 | ||
68 | using namespace std; | |
69 | namespace OpenMD { | |
# | Line 68 | 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), topologyDone_(false), |
78 | > | nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), |
79 | calcBoxDipole_(false), useAtomicVirial_(true) { | |
80 | ||
81 | MoleculeStamp* molStamp; | |
# | Line 125 | 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 | |
128 | – | std::cerr << "nGA = " << nGlobalAtoms_ << "\n"; |
129 | – | std::cerr << "nCA = " << nCutoffAtoms << "\n"; |
130 | – | std::cerr << "nG = " << nGroups << "\n"; |
132 | ||
133 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
133 | – | |
134 | – | std::cerr << "nGCG = " << nGlobalCutoffGroups_ << "\n"; |
134 | ||
135 | //every free atom (atom does not belong to rigid bodies) is an | |
136 | //integrable object therefore the total number of integrable objects | |
# | 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 | + | cerr << "ndfLocal_ = " << ndfLocal_ << "\n"; |
264 | + | |
265 | // n_constraints is local, so subtract them on each processor | |
266 | ndf_local -= nConstraints_; | |
267 | ||
268 | #ifdef IS_MPI | |
269 | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | |
270 | + | MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
271 | #else | |
272 | ndf_ = ndf_local; | |
273 | + | nGlobalFluctuatingCharges_ = nfq_local; |
274 | #endif | |
275 | ||
276 | // nZconstraints_ is global, as are the 3 COM translations for the | |
# | Line 274 | Line 287 | namespace OpenMD { | |
287 | #endif | |
288 | return fdf_; | |
289 | } | |
290 | + | |
291 | + | unsigned int SimInfo::getNLocalCutoffGroups(){ |
292 | + | int nLocalCutoffAtoms = 0; |
293 | + | Molecule* mol; |
294 | + | MoleculeIterator mi; |
295 | + | CutoffGroup* cg; |
296 | + | Molecule::CutoffGroupIterator ci; |
297 | ||
298 | + | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
299 | + | |
300 | + | for (cg = mol->beginCutoffGroup(ci); cg != NULL; |
301 | + | cg = mol->nextCutoffGroup(ci)) { |
302 | + | nLocalCutoffAtoms += cg->getNumAtom(); |
303 | + | |
304 | + | } |
305 | + | } |
306 | + | |
307 | + | return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; |
308 | + | } |
309 | + | |
310 | void SimInfo::calcNdfRaw() { | |
311 | int ndfRaw_local; | |
312 | ||
# | Line 680 | Line 712 | namespace OpenMD { | |
712 | Atom* atom; | |
713 | set<AtomType*> atomTypes; | |
714 | ||
715 | < | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
716 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
715 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
716 | > | for(atom = mol->beginAtom(ai); atom != NULL; |
717 | > | atom = mol->nextAtom(ai)) { |
718 | atomTypes.insert(atom->getAtomType()); | |
719 | } | |
720 | } | |
721 | < | |
721 | > | |
722 | #ifdef IS_MPI | |
723 | ||
724 | // loop over the found atom types on this processor, and add their | |
725 | // numerical idents to a vector: | |
726 | < | |
726 | > | |
727 | vector<int> foundTypes; | |
728 | set<AtomType*>::iterator i; | |
729 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) | |
# | Line 699 | Line 732 | namespace OpenMD { | |
732 | // count_local holds the number of found types on this processor | |
733 | int count_local = foundTypes.size(); | |
734 | ||
735 | < | // 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); |
735 | > | int nproc = MPI::COMM_WORLD.Get_size(); |
736 | ||
737 | < | // create a vector to hold the globally found types, and resize it: |
738 | < | vector<int> ftGlobal; |
739 | < | ftGlobal.resize(count); |
740 | < | vector<int> counts; |
737 | > | // we need arrays to hold the counts and displacement vectors for |
738 | > | // all processors |
739 | > | vector<int> counts(nproc, 0); |
740 | > | vector<int> disps(nproc, 0); |
741 | ||
742 | < | int nproc = MPI::COMM_WORLD.Get_size(); |
743 | < | counts.resize(nproc); |
744 | < | vector<int> disps; |
745 | < | disps.resize(nproc); |
742 | > | // fill the counts array |
743 | > | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], |
744 | > | 1, MPI::INT); |
745 | > | |
746 | > | // use the processor counts to compute the displacement array |
747 | > | disps[0] = 0; |
748 | > | int totalCount = counts[0]; |
749 | > | for (int iproc = 1; iproc < nproc; iproc++) { |
750 | > | disps[iproc] = disps[iproc-1] + counts[iproc-1]; |
751 | > | totalCount += counts[iproc]; |
752 | > | } |
753 | ||
754 | < | // now spray out the foundTypes to all the other processors: |
754 | > | // we need a (possibly redundant) set of all found types: |
755 | > | vector<int> ftGlobal(totalCount); |
756 | ||
757 | + | // now spray out the foundTypes to all the other processors: |
758 | MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, | |
759 | < | &ftGlobal[0], &counts[0], &disps[0], MPI::INT); |
759 | > | &ftGlobal[0], &counts[0], &disps[0], |
760 | > | MPI::INT); |
761 | ||
762 | + | vector<int>::iterator j; |
763 | + | |
764 | // foundIdents is a stl set, so inserting an already found ident | |
765 | // will have no effect. | |
766 | set<int> foundIdents; | |
767 | < | vector<int>::iterator j; |
767 | > | |
768 | for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) | |
769 | foundIdents.insert((*j)); | |
770 | ||
771 | // now iterate over the foundIdents and get the actual atom types | |
772 | // that correspond to these: | |
773 | set<int>::iterator it; | |
774 | < | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
774 | > | for (it = foundIdents.begin(); it != foundIdents.end(); ++it) |
775 | atomTypes.insert( forceField_->getAtomType((*it)) ); | |
776 | ||
777 | #endif | |
778 | < | |
778 | > | |
779 | return atomTypes; | |
780 | } | |
781 | ||
# | Line 745 | Line 787 | namespace OpenMD { | |
787 | if ( simParams_->getAccumulateBoxDipole() ) { | |
788 | calcBoxDipole_ = true; | |
789 | } | |
790 | < | |
790 | > | |
791 | set<AtomType*>::iterator i; | |
792 | set<AtomType*> atomTypes; | |
793 | atomTypes = getSimulatedAtomTypes(); | |
794 | int usesElectrostatic = 0; | |
795 | int usesMetallic = 0; | |
796 | int usesDirectional = 0; | |
797 | + | int usesFluctuatingCharges = 0; |
798 | //loop over all of the atom types | |
799 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
800 | usesElectrostatic |= (*i)->isElectrostatic(); | |
801 | usesMetallic |= (*i)->isMetal(); | |
802 | usesDirectional |= (*i)->isDirectional(); | |
803 | + | usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); |
804 | } | |
805 | < | |
805 | > | |
806 | #ifdef IS_MPI | |
807 | int temp; | |
808 | temp = usesDirectional; | |
809 | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
810 | < | |
810 | > | |
811 | temp = usesMetallic; | |
812 | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
813 | < | |
813 | > | |
814 | temp = usesElectrostatic; | |
815 | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
816 | + | |
817 | + | temp = usesFluctuatingCharges; |
818 | + | MPI_Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
819 | + | #else |
820 | + | |
821 | + | usesDirectionalAtoms_ = usesDirectional; |
822 | + | usesMetallicAtoms_ = usesMetallic; |
823 | + | usesElectrostaticAtoms_ = usesElectrostatic; |
824 | + | usesFluctuatingCharges_ = usesFluctuatingCharges; |
825 | + | |
826 | #endif | |
827 | + | |
828 | + | requiresPrepair_ = usesMetallicAtoms_ ? true : false; |
829 | + | requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; |
830 | + | requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; |
831 | } | |
832 | ||
833 | ||
# | Line 824 | Line 882 | namespace OpenMD { | |
882 | Atom* atom; | |
883 | RealType totalMass; | |
884 | ||
885 | < | //to avoid memory reallocation, reserve enough space for massFactors_ |
885 | > | /** |
886 | > | * The mass factor is the relative mass of an atom to the total |
887 | > | * mass of the cutoff group it belongs to. By default, all atoms |
888 | > | * are their own cutoff groups, and therefore have mass factors of |
889 | > | * 1. We need some special handling for massless atoms, which |
890 | > | * will be treated as carrying the entire mass of the cutoff |
891 | > | * group. |
892 | > | */ |
893 | massFactors_.clear(); | |
894 | < | massFactors_.reserve(getNCutoffGroups()); |
894 | > | massFactors_.resize(getNAtoms(), 1.0); |
895 | ||
896 | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
897 | for (cg = mol->beginCutoffGroup(ci); cg != NULL; | |
# | Line 835 | Line 900 | namespace OpenMD { | |
900 | totalMass = cg->getMass(); | |
901 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
902 | // Check for massless groups - set mfact to 1 if true | |
903 | < | if (totalMass != 0) |
904 | < | massFactors_.push_back(atom->getMass()/totalMass); |
903 | > | if (totalMass != 0) |
904 | > | massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; |
905 | else | |
906 | < | massFactors_.push_back( 1.0 ); |
906 | > | massFactors_[atom->getLocalIndex()] = 1.0; |
907 | } | |
908 | } | |
909 | } | |
# | Line 865 | Line 930 | namespace OpenMD { | |
930 | int* oneThreeList = oneThreeInteractions_.getPairList(); | |
931 | int* oneFourList = oneFourInteractions_.getPairList(); | |
932 | ||
868 | – | //setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray_[0], |
869 | – | // &nExclude, excludeList, |
870 | – | // &nOneTwo, oneTwoList, |
871 | – | // &nOneThree, oneThreeList, |
872 | – | // &nOneFour, oneFourList, |
873 | – | // &molMembershipArray[0], &mfact[0], &nCutoffGroups_, |
874 | – | // &fortranGlobalGroupMembership[0], &isError); |
875 | – | |
933 | topologyDone_ = true; | |
934 | } | |
935 | ||
# | Line 1157 | Line 1214 | namespace OpenMD { | |
1214 | ||
1215 | det = intTensor.determinant(); | |
1216 | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | |
1217 | < | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1217 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(det); |
1218 | return; | |
1219 | } | |
1220 | ||
# | Line 1173 | Line 1230 | namespace OpenMD { | |
1230 | ||
1231 | detI = intTensor.determinant(); | |
1232 | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; | |
1233 | < | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1233 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(detI); |
1234 | return; | |
1235 | } | |
1236 | /* |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |