| 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" | 
| 85 |  |  | 
| 86 |  | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | 
| 87 |  | forceField_(ff), simParams_(simParams), | 
| 88 | < | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 88 | > | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 89 |  | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | 
| 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; | 
| 291 |  |  | 
| 292 |  | } | 
| 293 |  |  | 
| 294 | + | int SimInfo::getFdf() { | 
| 295 | + | #ifdef IS_MPI | 
| 296 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 297 | + | #else | 
| 298 | + | fdf_ = fdf_local; | 
| 299 | + | #endif | 
| 300 | + | return fdf_; | 
| 301 | + | } | 
| 302 | + |  | 
| 303 |  | void SimInfo::calcNdfRaw() { | 
| 304 |  | int ndfRaw_local; | 
| 305 |  |  | 
| 603 |  |  | 
| 604 |  | setupElectrostaticSummationMethod( isError ); | 
| 605 |  | setupSwitchingFunction(); | 
| 606 | + | setupAccumulateBoxDipole(); | 
| 607 |  |  | 
| 608 |  | if(isError){ | 
| 609 |  | sprintf( painCave.errMsg, | 
| 663 |  | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | 
| 664 |  | int useRF; | 
| 665 |  | int useSF; | 
| 666 | + | int useSP; | 
| 667 | + | int useBoxDipole; | 
| 668 |  | std::string myMethod; | 
| 669 |  |  | 
| 670 |  | // set the useRF logical | 
| 675 |  | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 676 |  | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 677 |  | toUpper(myMethod); | 
| 678 | < | if (myMethod == "REACTION_FIELD") { | 
| 678 | > | if (myMethod == "REACTION_FIELD"){ | 
| 679 |  | useRF=1; | 
| 680 | < | } else { | 
| 681 | < | if (myMethod == "SHIFTED_FORCE") { | 
| 682 | < | useSF = 1; | 
| 683 | < | } | 
| 680 | > | } else if (myMethod == "SHIFTED_FORCE"){ | 
| 681 | > | useSF = 1; | 
| 682 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ | 
| 683 | > | useSP = 1; | 
| 684 |  | } | 
| 685 |  | } | 
| 686 | + |  | 
| 687 | + | if (simParams_->haveAccumulateBoxDipole()) | 
| 688 | + | if (simParams_->getAccumulateBoxDipole()) | 
| 689 | + | useBoxDipole = 1; | 
| 690 |  |  | 
| 691 |  | //loop over all of the atom types | 
| 692 |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 757 |  | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 758 |  |  | 
| 759 |  | temp = useSF; | 
| 760 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 760 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 761 |  |  | 
| 762 | + | temp = useSP; | 
| 763 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 764 | + |  | 
| 765 | + | temp = useBoxDipole; | 
| 766 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 767 | + |  | 
| 768 |  | #endif | 
| 769 |  |  | 
| 770 |  | fInfo_.SIM_uses_PBC = usePBC; | 
| 782 |  | fInfo_.SIM_uses_FLARB = useFLARB; | 
| 783 |  | fInfo_.SIM_uses_RF = useRF; | 
| 784 |  | fInfo_.SIM_uses_SF = useSF; | 
| 785 | + | fInfo_.SIM_uses_SP = useSP; | 
| 786 | + | fInfo_.SIM_uses_BoxDipole = useBoxDipole; | 
| 787 |  |  | 
| 788 |  | if( myMethod == "REACTION_FIELD") { | 
| 789 |  |  | 
| 815 |  | } | 
| 816 |  |  | 
| 817 |  | //calculate mass ratio of cutoff group | 
| 818 | < | std::vector<double> mfact; | 
| 818 | > | std::vector<RealType> mfact; | 
| 819 |  | SimInfo::MoleculeIterator mi; | 
| 820 |  | Molecule* mol; | 
| 821 |  | Molecule::CutoffGroupIterator ci; | 
| 822 |  | CutoffGroup* cg; | 
| 823 |  | Molecule::AtomIterator ai; | 
| 824 |  | Atom* atom; | 
| 825 | < | double totalMass; | 
| 825 | > | RealType totalMass; | 
| 826 |  |  | 
| 827 |  | //to avoid memory reallocation, reserve enough space for mfact | 
| 828 |  | mfact.reserve(getNCutoffGroups()); | 
| 982 |  | notifyFortranCutoffPolicy(&cp); | 
| 983 |  |  | 
| 984 |  | // Check the Skin Thickness for neighborlists | 
| 985 | < | double skin; | 
| 985 | > | RealType skin; | 
| 986 |  | if (simParams_->haveSkinThickness()) { | 
| 987 |  | skin = simParams_->getSkinThickness(); | 
| 988 |  | notifyFortranSkinThickness(&skin); | 
| 1072 |  | int errorOut; | 
| 1073 |  | int esm =  NONE; | 
| 1074 |  | int sm = UNDAMPED; | 
| 1075 | < | double alphaVal; | 
| 1076 | < | double dielectric; | 
| 1075 | > | RealType alphaVal; | 
| 1076 | > | RealType dielectric; | 
| 1077 |  |  | 
| 1078 |  | errorOut = isError; | 
| 1079 |  | alphaVal = simParams_->getDampingAlpha(); | 
| 1174 |  |  | 
| 1175 |  | // send switching function notification to switcheroo | 
| 1176 |  | setFunctionType(&ft); | 
| 1177 | + |  | 
| 1178 | + | } | 
| 1179 |  |  | 
| 1180 | + | void SimInfo::setupAccumulateBoxDipole() { | 
| 1181 | + |  | 
| 1182 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 1183 | + | if ( simParams_->haveAccumulateBoxDipole() ) | 
| 1184 | + | if ( simParams_->getAccumulateBoxDipole() ) { | 
| 1185 | + | setAccumulateBoxDipole(); | 
| 1186 | + | calcBoxDipole_ = true; | 
| 1187 | + | } | 
| 1188 | + |  | 
| 1189 |  | } | 
| 1190 |  |  | 
| 1191 |  | void SimInfo::addProperty(GenericData* genData) { | 
| 1244 |  | Molecule* mol; | 
| 1245 |  |  | 
| 1246 |  | Vector3d comVel(0.0); | 
| 1247 | < | double totalMass = 0.0; | 
| 1247 | > | RealType totalMass = 0.0; | 
| 1248 |  |  | 
| 1249 |  |  | 
| 1250 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1251 | < | double mass = mol->getMass(); | 
| 1251 | > | RealType mass = mol->getMass(); | 
| 1252 |  | totalMass += mass; | 
| 1253 |  | comVel += mass * mol->getComVel(); | 
| 1254 |  | } | 
| 1255 |  |  | 
| 1256 |  | #ifdef IS_MPI | 
| 1257 | < | double tmpMass = totalMass; | 
| 1257 | > | RealType tmpMass = totalMass; | 
| 1258 |  | Vector3d tmpComVel(comVel); | 
| 1259 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1260 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1259 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1260 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1261 |  | #endif | 
| 1262 |  |  | 
| 1263 |  | comVel /= totalMass; | 
| 1270 |  | Molecule* mol; | 
| 1271 |  |  | 
| 1272 |  | Vector3d com(0.0); | 
| 1273 | < | double totalMass = 0.0; | 
| 1273 | > | RealType totalMass = 0.0; | 
| 1274 |  |  | 
| 1275 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1276 | < | double mass = mol->getMass(); | 
| 1276 | > | RealType mass = mol->getMass(); | 
| 1277 |  | totalMass += mass; | 
| 1278 |  | com += mass * mol->getCom(); | 
| 1279 |  | } | 
| 1280 |  |  | 
| 1281 |  | #ifdef IS_MPI | 
| 1282 | < | double tmpMass = totalMass; | 
| 1282 | > | RealType tmpMass = totalMass; | 
| 1283 |  | Vector3d tmpCom(com); | 
| 1284 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1285 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1284 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1285 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1286 |  | #endif | 
| 1287 |  |  | 
| 1288 |  | com /= totalMass; | 
| 1306 |  | Molecule* mol; | 
| 1307 |  |  | 
| 1308 |  |  | 
| 1309 | < | double totalMass = 0.0; | 
| 1309 | > | RealType totalMass = 0.0; | 
| 1310 |  |  | 
| 1311 |  |  | 
| 1312 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1313 | < | double mass = mol->getMass(); | 
| 1313 | > | RealType mass = mol->getMass(); | 
| 1314 |  | totalMass += mass; | 
| 1315 |  | com += mass * mol->getCom(); | 
| 1316 |  | comVel += mass * mol->getComVel(); | 
| 1317 |  | } | 
| 1318 |  |  | 
| 1319 |  | #ifdef IS_MPI | 
| 1320 | < | double tmpMass = totalMass; | 
| 1320 | > | RealType tmpMass = totalMass; | 
| 1321 |  | Vector3d tmpCom(com); | 
| 1322 |  | Vector3d tmpComVel(comVel); | 
| 1323 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1324 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1325 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1323 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1324 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1325 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1326 |  | #endif | 
| 1327 |  |  | 
| 1328 |  | com /= totalMass; | 
| 1341 |  | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | 
| 1342 |  |  | 
| 1343 |  |  | 
| 1344 | < | double xx = 0.0; | 
| 1345 | < | double yy = 0.0; | 
| 1346 | < | double zz = 0.0; | 
| 1347 | < | double xy = 0.0; | 
| 1348 | < | double xz = 0.0; | 
| 1349 | < | double yz = 0.0; | 
| 1344 | > | RealType xx = 0.0; | 
| 1345 | > | RealType yy = 0.0; | 
| 1346 | > | RealType zz = 0.0; | 
| 1347 | > | RealType xy = 0.0; | 
| 1348 | > | RealType xz = 0.0; | 
| 1349 | > | RealType yz = 0.0; | 
| 1350 |  | Vector3d com(0.0); | 
| 1351 |  | Vector3d comVel(0.0); | 
| 1352 |  |  | 
| 1358 |  | Vector3d thisq(0.0); | 
| 1359 |  | Vector3d thisv(0.0); | 
| 1360 |  |  | 
| 1361 | < | double thisMass = 0.0; | 
| 1361 | > | RealType thisMass = 0.0; | 
| 1362 |  |  | 
| 1363 |  |  | 
| 1364 |  |  | 
| 1396 |  | #ifdef IS_MPI | 
| 1397 |  | Mat3x3d tmpI(inertiaTensor); | 
| 1398 |  | Vector3d tmpAngMom; | 
| 1399 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1400 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1399 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1400 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1401 |  | #endif | 
| 1402 |  |  | 
| 1403 |  | return; | 
| 1418 |  | Vector3d thisr(0.0); | 
| 1419 |  | Vector3d thisp(0.0); | 
| 1420 |  |  | 
| 1421 | < | double thisMass; | 
| 1421 | > | RealType thisMass; | 
| 1422 |  |  | 
| 1423 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1424 |  | thisMass = mol->getMass(); | 
| 1431 |  |  | 
| 1432 |  | #ifdef IS_MPI | 
| 1433 |  | Vector3d tmpAngMom; | 
| 1434 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1434 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1435 |  | #endif | 
| 1436 |  |  | 
| 1437 |  | return angularMomentum; | 
| 1438 |  | } | 
| 1439 |  |  | 
| 1440 | < |  | 
| 1440 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | 
| 1441 | > | return IOIndexToIntegrableObject.at(index); | 
| 1442 | > | } | 
| 1443 | > |  | 
| 1444 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { | 
| 1445 | > | IOIndexToIntegrableObject= v; | 
| 1446 | > | } | 
| 1447 | > |  | 
| 1448 | > | /* | 
| 1449 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { | 
| 1450 | > | assert( v.size() == nAtoms_ + nRigidBodies_); | 
| 1451 | > | sdByGlobalIndex_ = v; | 
| 1452 | > | } | 
| 1453 | > |  | 
| 1454 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { | 
| 1455 | > | //assert(index < nAtoms_ + nRigidBodies_); | 
| 1456 | > | return sdByGlobalIndex_.at(index); | 
| 1457 | > | } | 
| 1458 | > | */ | 
| 1459 |  | }//end namespace oopse | 
| 1460 |  |  |