# | Line 53 | Line 53 | |
---|---|---|
53 | #include "brains/SimInfo.hpp" | |
54 | #include "math/Vector3.hpp" | |
55 | #include "primitives/Molecule.hpp" | |
56 | + | #include "primitives/StuntDouble.hpp" |
57 | #include "UseTheForce/fCutoffPolicy.h" | |
58 | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" | |
59 | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" | |
# | Line 89 | Line 90 | namespace oopse { | |
90 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
91 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), | |
92 | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), | |
93 | < | sman_(NULL), fortranInitialized_(false) { |
93 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { |
94 | ||
95 | MoleculeStamp* molStamp; | |
96 | int nMolWithSameStamp; | |
# | Line 599 | Line 600 | namespace oopse { | |
600 | //setup fortran force field | |
601 | /** @deprecate */ | |
602 | int isError = 0; | |
603 | + | |
604 | + | setupCutoff(); |
605 | ||
606 | setupElectrostaticSummationMethod( isError ); | |
607 | setupSwitchingFunction(); | |
608 | + | setupAccumulateBoxDipole(); |
609 | ||
610 | if(isError){ | |
611 | sprintf( painCave.errMsg, | |
# | Line 609 | Line 613 | namespace oopse { | |
613 | painCave.isFatal = 1; | |
614 | simError(); | |
615 | } | |
612 | – | |
613 | – | |
614 | – | setupCutoff(); |
616 | ||
617 | calcNdf(); | |
618 | calcNdfRaw(); | |
# | Line 661 | 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 671 | 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 749 | 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 768 | 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 799 | 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 966 | 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 1056 | Line 1071 | namespace oopse { | |
1071 | int errorOut; | |
1072 | int esm = NONE; | |
1073 | int sm = UNDAMPED; | |
1074 | < | double alphaVal; |
1075 | < | double dielectric; |
1076 | < | |
1074 | > | RealType alphaVal; |
1075 | > | RealType dielectric; |
1076 | > | |
1077 | errorOut = isError; | |
1063 | – | alphaVal = simParams_->getDampingAlpha(); |
1078 | dielectric = simParams_->getDielectric(); | |
1079 | ||
1080 | if (simParams_->haveElectrostaticSummationMethod()) { | |
# | Line 1106 | Line 1120 | namespace oopse { | |
1120 | if (myScreen == "DAMPED") { | |
1121 | sm = DAMPED; | |
1122 | if (!simParams_->haveDampingAlpha()) { | |
1123 | < | //throw error |
1123 | > | // first set a cutoff dependent alpha value |
1124 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1125 | > | alphaVal = 0.5125 - rcut_* 0.025; |
1126 | > | // for values rcut > 20.5, alpha is zero |
1127 | > | if (alphaVal < 0) alphaVal = 0; |
1128 | > | |
1129 | > | // throw warning |
1130 | sprintf( painCave.errMsg, | |
1131 | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | |
1132 | < | "\tA default value of %f (1/ang) will be used.\n", alphaVal); |
1132 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1133 | painCave.isFatal = 0; | |
1134 | simError(); | |
1135 | } | |
# | Line 1161 | Line 1181 | namespace oopse { | |
1181 | ||
1182 | } | |
1183 | ||
1184 | + | void SimInfo::setupAccumulateBoxDipole() { |
1185 | + | |
1186 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1187 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1188 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1189 | + | setAccumulateBoxDipole(); |
1190 | + | calcBoxDipole_ = true; |
1191 | + | } |
1192 | + | |
1193 | + | } |
1194 | + | |
1195 | void SimInfo::addProperty(GenericData* genData) { | |
1196 | properties_.addProperty(genData); | |
1197 | } | |
# | Line 1217 | Line 1248 | namespace oopse { | |
1248 | Molecule* mol; | |
1249 | ||
1250 | Vector3d comVel(0.0); | |
1251 | < | double totalMass = 0.0; |
1251 | > | RealType totalMass = 0.0; |
1252 | ||
1253 | ||
1254 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1255 | < | double mass = mol->getMass(); |
1255 | > | RealType mass = mol->getMass(); |
1256 | totalMass += mass; | |
1257 | comVel += mass * mol->getComVel(); | |
1258 | } | |
1259 | ||
1260 | #ifdef IS_MPI | |
1261 | < | double tmpMass = totalMass; |
1261 | > | RealType tmpMass = totalMass; |
1262 | Vector3d tmpComVel(comVel); | |
1263 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1264 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1263 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1264 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1265 | #endif | |
1266 | ||
1267 | comVel /= totalMass; | |
# | Line 1243 | Line 1274 | namespace oopse { | |
1274 | Molecule* mol; | |
1275 | ||
1276 | Vector3d com(0.0); | |
1277 | < | double totalMass = 0.0; |
1277 | > | RealType totalMass = 0.0; |
1278 | ||
1279 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1280 | < | double mass = mol->getMass(); |
1280 | > | RealType mass = mol->getMass(); |
1281 | totalMass += mass; | |
1282 | com += mass * mol->getCom(); | |
1283 | } | |
1284 | ||
1285 | #ifdef IS_MPI | |
1286 | < | double tmpMass = totalMass; |
1286 | > | RealType tmpMass = totalMass; |
1287 | Vector3d tmpCom(com); | |
1288 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1289 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1288 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1289 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1290 | #endif | |
1291 | ||
1292 | com /= totalMass; | |
# | Line 1279 | Line 1310 | namespace oopse { | |
1310 | Molecule* mol; | |
1311 | ||
1312 | ||
1313 | < | double totalMass = 0.0; |
1313 | > | RealType totalMass = 0.0; |
1314 | ||
1315 | ||
1316 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1317 | < | double mass = mol->getMass(); |
1317 | > | RealType mass = mol->getMass(); |
1318 | totalMass += mass; | |
1319 | com += mass * mol->getCom(); | |
1320 | comVel += mass * mol->getComVel(); | |
1321 | } | |
1322 | ||
1323 | #ifdef IS_MPI | |
1324 | < | double tmpMass = totalMass; |
1324 | > | RealType tmpMass = totalMass; |
1325 | Vector3d tmpCom(com); | |
1326 | Vector3d tmpComVel(comVel); | |
1327 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1328 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1329 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1327 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1328 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1329 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1330 | #endif | |
1331 | ||
1332 | com /= totalMass; | |
# | Line 1314 | Line 1345 | namespace oopse { | |
1345 | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | |
1346 | ||
1347 | ||
1348 | < | double xx = 0.0; |
1349 | < | double yy = 0.0; |
1350 | < | double zz = 0.0; |
1351 | < | double xy = 0.0; |
1352 | < | double xz = 0.0; |
1353 | < | double yz = 0.0; |
1348 | > | RealType xx = 0.0; |
1349 | > | RealType yy = 0.0; |
1350 | > | RealType zz = 0.0; |
1351 | > | RealType xy = 0.0; |
1352 | > | RealType xz = 0.0; |
1353 | > | RealType yz = 0.0; |
1354 | Vector3d com(0.0); | |
1355 | Vector3d comVel(0.0); | |
1356 | ||
# | Line 1331 | Line 1362 | namespace oopse { | |
1362 | Vector3d thisq(0.0); | |
1363 | Vector3d thisv(0.0); | |
1364 | ||
1365 | < | double thisMass = 0.0; |
1365 | > | RealType thisMass = 0.0; |
1366 | ||
1367 | ||
1368 | ||
# | Line 1369 | Line 1400 | namespace oopse { | |
1400 | #ifdef IS_MPI | |
1401 | Mat3x3d tmpI(inertiaTensor); | |
1402 | Vector3d tmpAngMom; | |
1403 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1404 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1403 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1404 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1405 | #endif | |
1406 | ||
1407 | return; | |
# | Line 1391 | Line 1422 | namespace oopse { | |
1422 | Vector3d thisr(0.0); | |
1423 | Vector3d thisp(0.0); | |
1424 | ||
1425 | < | double thisMass; |
1425 | > | RealType thisMass; |
1426 | ||
1427 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1428 | thisMass = mol->getMass(); | |
# | Line 1404 | Line 1435 | namespace oopse { | |
1435 | ||
1436 | #ifdef IS_MPI | |
1437 | Vector3d tmpAngMom; | |
1438 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1438 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1439 | #endif | |
1440 | ||
1441 | return angularMomentum; | |
1442 | } | |
1443 | ||
1444 | < | |
1444 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { |
1445 | > | return IOIndexToIntegrableObject.at(index); |
1446 | > | } |
1447 | > | |
1448 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1449 | > | IOIndexToIntegrableObject= v; |
1450 | > | } |
1451 | > | |
1452 | > | /* |
1453 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1454 | > | assert( v.size() == nAtoms_ + nRigidBodies_); |
1455 | > | sdByGlobalIndex_ = v; |
1456 | > | } |
1457 | > | |
1458 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1459 | > | //assert(index < nAtoms_ + nRigidBodies_); |
1460 | > | return sdByGlobalIndex_.at(index); |
1461 | > | } |
1462 | > | */ |
1463 | }//end namespace oopse | |
1464 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |