# | Line 63 | Line 63 | |
---|---|---|
63 | #include "utils/MemoryUtils.hpp" | |
64 | #include "utils/simError.h" | |
65 | #include "selection/SelectionManager.hpp" | |
66 | + | #include "io/ForceFieldOptions.hpp" |
67 | + | #include "UseTheForce/ForceField.hpp" |
68 | ||
69 | #ifdef IS_MPI | |
70 | #include "UseTheForce/mpiComponentPlan.h" | |
# | Line 80 | Line 82 | namespace oopse { | |
82 | return result; | |
83 | } | |
84 | ||
85 | < | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, |
86 | < | ForceField* ff, Globals* simParams) : |
87 | < | stamps_(stamps), forceField_(ff), simParams_(simParams), |
86 | < | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
85 | > | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : |
86 | > | forceField_(ff), simParams_(simParams), |
87 | > | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
88 | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | |
89 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
90 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), | |
91 | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), | |
92 | < | sman_(NULL), fortranInitialized_(false) { |
92 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { |
93 | ||
93 | – | |
94 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; |
94 | MoleculeStamp* molStamp; | |
95 | int nMolWithSameStamp; | |
96 | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | |
# | Line 99 | Line 98 | namespace oopse { | |
98 | CutoffGroupStamp* cgStamp; | |
99 | RigidBodyStamp* rbStamp; | |
100 | int nRigidAtoms = 0; | |
101 | < | |
102 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { |
103 | < | molStamp = i->first; |
104 | < | nMolWithSameStamp = i->second; |
101 | > | std::vector<Component*> components = simParams->getComponents(); |
102 | > | |
103 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
104 | > | molStamp = (*i)->getMoleculeStamp(); |
105 | > | nMolWithSameStamp = (*i)->getNMol(); |
106 | ||
107 | addMoleculeStamp(molStamp, nMolWithSameStamp); | |
108 | ||
109 | //calculate atoms in molecules | |
110 | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | |
111 | – | |
111 | ||
112 | //calculate atoms in cutoff groups | |
113 | int nAtomsInGroups = 0; | |
114 | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | |
115 | ||
116 | for (int j=0; j < nCutoffGroupsInStamp; j++) { | |
117 | < | cgStamp = molStamp->getCutoffGroup(j); |
117 | > | cgStamp = molStamp->getCutoffGroupStamp(j); |
118 | nAtomsInGroups += cgStamp->getNMembers(); | |
119 | } | |
120 | ||
# | Line 128 | Line 127 | namespace oopse { | |
127 | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | |
128 | ||
129 | for (int j=0; j < nRigidBodiesInStamp; j++) { | |
130 | < | rbStamp = molStamp->getRigidBody(j); |
130 | > | rbStamp = molStamp->getRigidBodyStamp(j); |
131 | nAtomsInRigidBodies += rbStamp->getNMembers(); | |
132 | } | |
133 | ||
# | Line 167 | Line 166 | namespace oopse { | |
166 | } | |
167 | molecules_.clear(); | |
168 | ||
170 | – | delete stamps_; |
169 | delete sman_; | |
170 | delete simParams_; | |
171 | delete forceField_; | |
# | Line 274 | Line 272 | namespace oopse { | |
272 | } | |
273 | } | |
274 | ||
275 | < | }//end for (integrableObject) |
276 | < | }// end for (mol) |
275 | > | } |
276 | > | } |
277 | ||
278 | // n_constraints is local, so subtract them on each processor | |
279 | ndf_local -= nConstraints_; | |
# | Line 292 | Line 290 | namespace oopse { | |
290 | ||
291 | } | |
292 | ||
293 | + | int SimInfo::getFdf() { |
294 | + | #ifdef IS_MPI |
295 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
296 | + | #else |
297 | + | fdf_ = fdf_local; |
298 | + | #endif |
299 | + | return fdf_; |
300 | + | } |
301 | + | |
302 | void SimInfo::calcNdfRaw() { | |
303 | int ndfRaw_local; | |
304 | ||
# | Line 595 | Line 602 | namespace oopse { | |
602 | ||
603 | setupElectrostaticSummationMethod( isError ); | |
604 | setupSwitchingFunction(); | |
605 | + | setupAccumulateBoxDipole(); |
606 | ||
607 | if(isError){ | |
608 | sprintf( painCave.errMsg, | |
# | Line 654 | Line 662 | namespace oopse { | |
662 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
663 | int useRF; | |
664 | int useSF; | |
665 | + | int useSP; |
666 | + | int useBoxDipole; |
667 | std::string myMethod; | |
668 | ||
669 | // set the useRF logical | |
# | Line 664 | Line 674 | namespace oopse { | |
674 | if (simParams_->haveElectrostaticSummationMethod()) { | |
675 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
676 | toUpper(myMethod); | |
677 | < | if (myMethod == "REACTION_FIELD") { |
677 | > | if (myMethod == "REACTION_FIELD"){ |
678 | useRF=1; | |
679 | < | } else { |
680 | < | if (myMethod == "SHIFTED_FORCE") { |
681 | < | useSF = 1; |
682 | < | } |
679 | > | } else if (myMethod == "SHIFTED_FORCE"){ |
680 | > | useSF = 1; |
681 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
682 | > | useSP = 1; |
683 | } | |
684 | } | |
685 | + | |
686 | + | if (simParams_->haveAccumulateBoxDipole()) |
687 | + | if (simParams_->getAccumulateBoxDipole()) |
688 | + | useBoxDipole = 1; |
689 | ||
690 | //loop over all of the atom types | |
691 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
# | Line 742 | Line 756 | namespace oopse { | |
756 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
757 | ||
758 | temp = useSF; | |
759 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
759 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
760 | ||
761 | + | temp = useSP; |
762 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
763 | + | |
764 | + | temp = useBoxDipole; |
765 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
766 | + | |
767 | #endif | |
768 | ||
769 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 761 | Line 781 | namespace oopse { | |
781 | fInfo_.SIM_uses_FLARB = useFLARB; | |
782 | fInfo_.SIM_uses_RF = useRF; | |
783 | fInfo_.SIM_uses_SF = useSF; | |
784 | + | fInfo_.SIM_uses_SP = useSP; |
785 | + | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
786 | ||
787 | if( myMethod == "REACTION_FIELD") { | |
788 | ||
# | Line 792 | Line 814 | namespace oopse { | |
814 | } | |
815 | ||
816 | //calculate mass ratio of cutoff group | |
817 | < | std::vector<double> mfact; |
817 | > | std::vector<RealType> mfact; |
818 | SimInfo::MoleculeIterator mi; | |
819 | Molecule* mol; | |
820 | Molecule::CutoffGroupIterator ci; | |
821 | CutoffGroup* cg; | |
822 | Molecule::AtomIterator ai; | |
823 | Atom* atom; | |
824 | < | double totalMass; |
824 | > | RealType totalMass; |
825 | ||
826 | //to avoid memory reallocation, reserve enough space for mfact | |
827 | mfact.reserve(getNCutoffGroups()); | |
# | Line 924 | Line 946 | namespace oopse { | |
946 | ||
947 | void SimInfo::setupCutoff() { | |
948 | ||
949 | + | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
950 | + | |
951 | // Check the cutoff policy | |
952 | < | int cp = TRADITIONAL_CUTOFF_POLICY; |
953 | < | if (simParams_->haveCutoffPolicy()) { |
954 | < | std::string myPolicy = simParams_->getCutoffPolicy(); |
952 | > | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default |
953 | > | |
954 | > | std::string myPolicy; |
955 | > | if (forceFieldOptions_.haveCutoffPolicy()){ |
956 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); |
957 | > | }else if (simParams_->haveCutoffPolicy()) { |
958 | > | myPolicy = simParams_->getCutoffPolicy(); |
959 | > | } |
960 | > | |
961 | > | if (!myPolicy.empty()){ |
962 | toUpper(myPolicy); | |
963 | if (myPolicy == "MIX") { | |
964 | cp = MIX_CUTOFF_POLICY; | |
# | Line 950 | Line 981 | namespace oopse { | |
981 | notifyFortranCutoffPolicy(&cp); | |
982 | ||
983 | // Check the Skin Thickness for neighborlists | |
984 | < | double skin; |
984 | > | RealType skin; |
985 | if (simParams_->haveSkinThickness()) { | |
986 | skin = simParams_->getSkinThickness(); | |
987 | notifyFortranSkinThickness(&skin); | |
# | Line 962 | Line 993 | namespace oopse { | |
993 | if (simParams_->haveSwitchingRadius()) { | |
994 | rsw_ = simParams_->getSwitchingRadius(); | |
995 | } else { | |
996 | < | rsw_ = rcut_; |
996 | > | if (fInfo_.SIM_uses_Charges | |
997 | > | fInfo_.SIM_uses_Dipoles | |
998 | > | fInfo_.SIM_uses_RF) { |
999 | > | |
1000 | > | rsw_ = 0.85 * rcut_; |
1001 | > | sprintf(painCave.errMsg, |
1002 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1003 | > | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" |
1004 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1005 | > | painCave.isFatal = 0; |
1006 | > | simError(); |
1007 | > | } else { |
1008 | > | rsw_ = rcut_; |
1009 | > | sprintf(painCave.errMsg, |
1010 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1011 | > | "\tOOPSE will use the same value as the cutoffRadius.\n" |
1012 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1013 | > | painCave.isFatal = 0; |
1014 | > | simError(); |
1015 | > | } |
1016 | } | |
1017 | + | |
1018 | notifyFortranCutoffs(&rcut_, &rsw_); | |
1019 | ||
1020 | } else { | |
# | Line 1020 | Line 1071 | namespace oopse { | |
1071 | int errorOut; | |
1072 | int esm = NONE; | |
1073 | int sm = UNDAMPED; | |
1074 | < | double alphaVal; |
1075 | < | double dielectric; |
1074 | > | RealType alphaVal; |
1075 | > | RealType dielectric; |
1076 | ||
1077 | errorOut = isError; | |
1078 | alphaVal = simParams_->getDampingAlpha(); | |
# | Line 1092 | Line 1143 | namespace oopse { | |
1143 | ||
1144 | // let's pass some summation method variables to fortran | |
1145 | setElectrostaticSummationMethod( &esm ); | |
1146 | < | notifyFortranElectrostaticMethod( &esm ); |
1146 | > | setFortranElectrostaticMethod( &esm ); |
1147 | setScreeningMethod( &sm ); | |
1148 | setDampingAlpha( &alphaVal ); | |
1149 | setReactionFieldDielectric( &dielectric ); | |
# | Line 1125 | Line 1176 | namespace oopse { | |
1176 | ||
1177 | } | |
1178 | ||
1179 | + | void SimInfo::setupAccumulateBoxDipole() { |
1180 | + | |
1181 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1182 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1183 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1184 | + | setAccumulateBoxDipole(); |
1185 | + | calcBoxDipole_ = true; |
1186 | + | } |
1187 | + | |
1188 | + | } |
1189 | + | |
1190 | void SimInfo::addProperty(GenericData* genData) { | |
1191 | properties_.addProperty(genData); | |
1192 | } | |
# | Line 1181 | Line 1243 | namespace oopse { | |
1243 | Molecule* mol; | |
1244 | ||
1245 | Vector3d comVel(0.0); | |
1246 | < | double totalMass = 0.0; |
1246 | > | RealType totalMass = 0.0; |
1247 | ||
1248 | ||
1249 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1250 | < | double mass = mol->getMass(); |
1250 | > | RealType mass = mol->getMass(); |
1251 | totalMass += mass; | |
1252 | comVel += mass * mol->getComVel(); | |
1253 | } | |
1254 | ||
1255 | #ifdef IS_MPI | |
1256 | < | double tmpMass = totalMass; |
1256 | > | RealType tmpMass = totalMass; |
1257 | Vector3d tmpComVel(comVel); | |
1258 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1259 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1258 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1259 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1260 | #endif | |
1261 | ||
1262 | comVel /= totalMass; | |
# | Line 1207 | Line 1269 | namespace oopse { | |
1269 | Molecule* mol; | |
1270 | ||
1271 | Vector3d com(0.0); | |
1272 | < | double totalMass = 0.0; |
1272 | > | RealType totalMass = 0.0; |
1273 | ||
1274 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1275 | < | double mass = mol->getMass(); |
1275 | > | RealType mass = mol->getMass(); |
1276 | totalMass += mass; | |
1277 | com += mass * mol->getCom(); | |
1278 | } | |
1279 | ||
1280 | #ifdef IS_MPI | |
1281 | < | double tmpMass = totalMass; |
1281 | > | RealType tmpMass = totalMass; |
1282 | Vector3d tmpCom(com); | |
1283 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1284 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1283 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1284 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1285 | #endif | |
1286 | ||
1287 | com /= totalMass; | |
# | Line 1243 | Line 1305 | namespace oopse { | |
1305 | Molecule* mol; | |
1306 | ||
1307 | ||
1308 | < | double totalMass = 0.0; |
1308 | > | RealType totalMass = 0.0; |
1309 | ||
1310 | ||
1311 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1312 | < | double mass = mol->getMass(); |
1312 | > | RealType mass = mol->getMass(); |
1313 | totalMass += mass; | |
1314 | com += mass * mol->getCom(); | |
1315 | comVel += mass * mol->getComVel(); | |
1316 | } | |
1317 | ||
1318 | #ifdef IS_MPI | |
1319 | < | double tmpMass = totalMass; |
1319 | > | RealType tmpMass = totalMass; |
1320 | Vector3d tmpCom(com); | |
1321 | Vector3d tmpComVel(comVel); | |
1322 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1323 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1324 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1322 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1323 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1324 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1325 | #endif | |
1326 | ||
1327 | com /= totalMass; | |
# | Line 1278 | Line 1340 | namespace oopse { | |
1340 | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | |
1341 | ||
1342 | ||
1343 | < | double xx = 0.0; |
1344 | < | double yy = 0.0; |
1345 | < | double zz = 0.0; |
1346 | < | double xy = 0.0; |
1347 | < | double xz = 0.0; |
1348 | < | double yz = 0.0; |
1343 | > | RealType xx = 0.0; |
1344 | > | RealType yy = 0.0; |
1345 | > | RealType zz = 0.0; |
1346 | > | RealType xy = 0.0; |
1347 | > | RealType xz = 0.0; |
1348 | > | RealType yz = 0.0; |
1349 | Vector3d com(0.0); | |
1350 | Vector3d comVel(0.0); | |
1351 | ||
# | Line 1295 | Line 1357 | namespace oopse { | |
1357 | Vector3d thisq(0.0); | |
1358 | Vector3d thisv(0.0); | |
1359 | ||
1360 | < | double thisMass = 0.0; |
1360 | > | RealType thisMass = 0.0; |
1361 | ||
1362 | ||
1363 | ||
# | Line 1333 | Line 1395 | namespace oopse { | |
1395 | #ifdef IS_MPI | |
1396 | Mat3x3d tmpI(inertiaTensor); | |
1397 | Vector3d tmpAngMom; | |
1398 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1399 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1398 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1399 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1400 | #endif | |
1401 | ||
1402 | return; | |
# | Line 1355 | Line 1417 | namespace oopse { | |
1417 | Vector3d thisr(0.0); | |
1418 | Vector3d thisp(0.0); | |
1419 | ||
1420 | < | double thisMass; |
1420 | > | RealType thisMass; |
1421 | ||
1422 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1423 | thisMass = mol->getMass(); | |
# | Line 1368 | Line 1430 | namespace oopse { | |
1430 | ||
1431 | #ifdef IS_MPI | |
1432 | Vector3d tmpAngMom; | |
1433 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1433 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1434 | #endif | |
1435 | ||
1436 | return angularMomentum; |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |