# | 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" | |
60 | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" | |
61 | #include "UseTheForce/doForces_interface.h" | |
62 | + | #include "UseTheForce/DarkSide/neighborLists_interface.h" |
63 | #include "UseTheForce/DarkSide/electrostatic_interface.h" | |
64 | #include "UseTheForce/DarkSide/switcheroo_interface.h" | |
65 | #include "utils/MemoryUtils.hpp" | |
66 | #include "utils/simError.h" | |
67 | #include "selection/SelectionManager.hpp" | |
68 | + | #include "io/ForceFieldOptions.hpp" |
69 | + | #include "UseTheForce/ForceField.hpp" |
70 | ||
71 | + | |
72 | #ifdef IS_MPI | |
73 | #include "UseTheForce/mpiComponentPlan.h" | |
74 | #include "UseTheForce/DarkSide/simParallel_interface.h" | |
# | Line 82 | Line 87 | namespace oopse { | |
87 | ||
88 | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | |
89 | forceField_(ff), simParams_(simParams), | |
90 | < | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
90 | > | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
91 | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | |
92 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
93 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), | |
94 | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), | |
95 | < | sman_(NULL), fortranInitialized_(false) { |
95 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false) { |
96 | ||
97 | MoleculeStamp* molStamp; | |
98 | int nMolWithSameStamp; | |
# | Line 288 | Line 293 | namespace oopse { | |
293 | ||
294 | } | |
295 | ||
296 | + | int SimInfo::getFdf() { |
297 | + | #ifdef IS_MPI |
298 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
299 | + | #else |
300 | + | fdf_ = fdf_local; |
301 | + | #endif |
302 | + | return fdf_; |
303 | + | } |
304 | + | |
305 | void SimInfo::calcNdfRaw() { | |
306 | int ndfRaw_local; | |
307 | ||
# | Line 589 | Line 603 | namespace oopse { | |
603 | /** @deprecate */ | |
604 | int isError = 0; | |
605 | ||
606 | + | setupCutoff(); |
607 | + | |
608 | setupElectrostaticSummationMethod( isError ); | |
609 | setupSwitchingFunction(); | |
610 | + | setupAccumulateBoxDipole(); |
611 | ||
612 | if(isError){ | |
613 | sprintf( painCave.errMsg, | |
# | Line 598 | Line 615 | namespace oopse { | |
615 | painCave.isFatal = 1; | |
616 | simError(); | |
617 | } | |
601 | – | |
602 | – | |
603 | – | setupCutoff(); |
618 | ||
619 | calcNdf(); | |
620 | calcNdfRaw(); | |
# | Line 650 | Line 664 | namespace oopse { | |
664 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
665 | int useRF; | |
666 | int useSF; | |
667 | + | int useSP; |
668 | + | int useBoxDipole; |
669 | std::string myMethod; | |
670 | ||
671 | // set the useRF logical | |
672 | useRF = 0; | |
673 | useSF = 0; | |
674 | + | useSP = 0; |
675 | ||
676 | ||
677 | if (simParams_->haveElectrostaticSummationMethod()) { | |
678 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
679 | toUpper(myMethod); | |
680 | < | if (myMethod == "REACTION_FIELD") { |
681 | < | useRF=1; |
682 | < | } else { |
683 | < | if (myMethod == "SHIFTED_FORCE") { |
684 | < | useSF = 1; |
685 | < | } |
680 | > | if (myMethod == "REACTION_FIELD"){ |
681 | > | useRF = 1; |
682 | > | } else if (myMethod == "SHIFTED_FORCE"){ |
683 | > | useSF = 1; |
684 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
685 | > | useSP = 1; |
686 | } | |
687 | } | |
688 | + | |
689 | + | if (simParams_->haveAccumulateBoxDipole()) |
690 | + | if (simParams_->getAccumulateBoxDipole()) |
691 | + | useBoxDipole = 1; |
692 | ||
693 | //loop over all of the atom types | |
694 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
# | Line 738 | Line 759 | namespace oopse { | |
759 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
760 | ||
761 | temp = useSF; | |
762 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
762 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
763 | ||
764 | + | temp = useSP; |
765 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
766 | + | |
767 | + | temp = useBoxDipole; |
768 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
769 | + | |
770 | #endif | |
771 | ||
772 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 757 | Line 784 | namespace oopse { | |
784 | fInfo_.SIM_uses_FLARB = useFLARB; | |
785 | fInfo_.SIM_uses_RF = useRF; | |
786 | fInfo_.SIM_uses_SF = useSF; | |
787 | < | |
788 | < | if( myMethod == "REACTION_FIELD") { |
762 | < | |
763 | < | if (simParams_->haveDielectric()) { |
764 | < | fInfo_.dielect = simParams_->getDielectric(); |
765 | < | } else { |
766 | < | sprintf(painCave.errMsg, |
767 | < | "SimSetup Error: No Dielectric constant was set.\n" |
768 | < | "\tYou are trying to use Reaction Field without" |
769 | < | "\tsetting a dielectric constant!\n"); |
770 | < | painCave.isFatal = 1; |
771 | < | simError(); |
772 | < | } |
773 | < | } |
774 | < | |
787 | > | fInfo_.SIM_uses_SP = useSP; |
788 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
789 | } | |
790 | ||
791 | void SimInfo::setupFortranSim() { | |
# | Line 788 | Line 802 | namespace oopse { | |
802 | } | |
803 | ||
804 | //calculate mass ratio of cutoff group | |
805 | < | std::vector<double> mfact; |
805 | > | std::vector<RealType> mfact; |
806 | SimInfo::MoleculeIterator mi; | |
807 | Molecule* mol; | |
808 | Molecule::CutoffGroupIterator ci; | |
809 | CutoffGroup* cg; | |
810 | Molecule::AtomIterator ai; | |
811 | Atom* atom; | |
812 | < | double totalMass; |
812 | > | RealType totalMass; |
813 | ||
814 | //to avoid memory reallocation, reserve enough space for mfact | |
815 | mfact.reserve(getNCutoffGroups()); | |
# | Line 856 | Line 870 | namespace oopse { | |
870 | "succesfully sent the simulation information to fortran.\n"); | |
871 | MPIcheckPoint(); | |
872 | #endif // is_mpi | |
873 | + | |
874 | + | // Setup number of neighbors in neighbor list if present |
875 | + | if (simParams_->haveNeighborListNeighbors()) { |
876 | + | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
877 | + | setNeighbors(&nlistNeighbors); |
878 | + | } |
879 | + | |
880 | + | |
881 | } | |
882 | ||
883 | ||
# | Line 920 | Line 942 | namespace oopse { | |
942 | ||
943 | void SimInfo::setupCutoff() { | |
944 | ||
945 | + | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
946 | + | |
947 | // Check the cutoff policy | |
948 | < | int cp = TRADITIONAL_CUTOFF_POLICY; |
949 | < | if (simParams_->haveCutoffPolicy()) { |
950 | < | std::string myPolicy = simParams_->getCutoffPolicy(); |
948 | > | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default |
949 | > | |
950 | > | std::string myPolicy; |
951 | > | if (forceFieldOptions_.haveCutoffPolicy()){ |
952 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); |
953 | > | }else if (simParams_->haveCutoffPolicy()) { |
954 | > | myPolicy = simParams_->getCutoffPolicy(); |
955 | > | } |
956 | > | |
957 | > | if (!myPolicy.empty()){ |
958 | toUpper(myPolicy); | |
959 | if (myPolicy == "MIX") { | |
960 | cp = MIX_CUTOFF_POLICY; | |
# | Line 946 | Line 977 | namespace oopse { | |
977 | notifyFortranCutoffPolicy(&cp); | |
978 | ||
979 | // Check the Skin Thickness for neighborlists | |
980 | < | double skin; |
980 | > | RealType skin; |
981 | if (simParams_->haveSkinThickness()) { | |
982 | skin = simParams_->getSkinThickness(); | |
983 | notifyFortranSkinThickness(&skin); | |
# | Line 958 | Line 989 | namespace oopse { | |
989 | if (simParams_->haveSwitchingRadius()) { | |
990 | rsw_ = simParams_->getSwitchingRadius(); | |
991 | } else { | |
992 | < | rsw_ = rcut_; |
992 | > | if (fInfo_.SIM_uses_Charges | |
993 | > | fInfo_.SIM_uses_Dipoles | |
994 | > | fInfo_.SIM_uses_RF) { |
995 | > | |
996 | > | rsw_ = 0.85 * rcut_; |
997 | > | sprintf(painCave.errMsg, |
998 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
999 | > | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" |
1000 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1001 | > | painCave.isFatal = 0; |
1002 | > | simError(); |
1003 | > | } else { |
1004 | > | rsw_ = rcut_; |
1005 | > | sprintf(painCave.errMsg, |
1006 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" |
1007 | > | "\tOOPSE will use the same value as the cutoffRadius.\n" |
1008 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); |
1009 | > | painCave.isFatal = 0; |
1010 | > | simError(); |
1011 | > | } |
1012 | } | |
1013 | + | |
1014 | notifyFortranCutoffs(&rcut_, &rsw_); | |
1015 | ||
1016 | } else { | |
# | Line 1016 | Line 1067 | namespace oopse { | |
1067 | int errorOut; | |
1068 | int esm = NONE; | |
1069 | int sm = UNDAMPED; | |
1070 | < | double alphaVal; |
1071 | < | double dielectric; |
1072 | < | |
1070 | > | RealType alphaVal; |
1071 | > | RealType dielectric; |
1072 | > | |
1073 | errorOut = isError; | |
1023 | – | alphaVal = simParams_->getDampingAlpha(); |
1024 | – | dielectric = simParams_->getDielectric(); |
1074 | ||
1075 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1076 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
# | Line 1038 | Line 1087 | namespace oopse { | |
1087 | if (myMethod == "SHIFTED_FORCE") { | |
1088 | esm = SHIFTED_FORCE; | |
1089 | } else { | |
1090 | < | if (myMethod == "REACTION_FIELD") { |
1090 | > | if (myMethod == "REACTION_FIELD") { |
1091 | esm = REACTION_FIELD; | |
1092 | + | dielectric = simParams_->getDielectric(); |
1093 | + | if (!simParams_->haveDielectric()) { |
1094 | + | // throw warning |
1095 | + | sprintf( painCave.errMsg, |
1096 | + | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1097 | + | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1098 | + | painCave.isFatal = 0; |
1099 | + | simError(); |
1100 | + | } |
1101 | } else { | |
1102 | // throw error | |
1103 | sprintf( painCave.errMsg, | |
# | Line 1066 | Line 1124 | namespace oopse { | |
1124 | if (myScreen == "DAMPED") { | |
1125 | sm = DAMPED; | |
1126 | if (!simParams_->haveDampingAlpha()) { | |
1127 | < | //throw error |
1127 | > | // first set a cutoff dependent alpha value |
1128 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1129 | > | alphaVal = 0.5125 - rcut_* 0.025; |
1130 | > | // for values rcut > 20.5, alpha is zero |
1131 | > | if (alphaVal < 0) alphaVal = 0; |
1132 | > | |
1133 | > | // throw warning |
1134 | sprintf( painCave.errMsg, | |
1135 | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | |
1136 | < | "\tA default value of %f (1/ang) will be used.\n", alphaVal); |
1136 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1137 | painCave.isFatal = 0; | |
1138 | simError(); | |
1139 | + | } else { |
1140 | + | alphaVal = simParams_->getDampingAlpha(); |
1141 | } | |
1142 | + | |
1143 | } else { | |
1144 | // throw error | |
1145 | sprintf( painCave.errMsg, | |
# | Line 1087 | Line 1154 | namespace oopse { | |
1154 | } | |
1155 | ||
1156 | // let's pass some summation method variables to fortran | |
1157 | < | setElectrostaticSumMethod( &esm ); |
1157 | > | setElectrostaticSummationMethod( &esm ); |
1158 | setFortranElectrostaticMethod( &esm ); | |
1159 | setScreeningMethod( &sm ); | |
1160 | setDampingAlpha( &alphaVal ); | |
# | Line 1121 | Line 1188 | namespace oopse { | |
1188 | ||
1189 | } | |
1190 | ||
1191 | + | void SimInfo::setupAccumulateBoxDipole() { |
1192 | + | |
1193 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1194 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1195 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1196 | + | setAccumulateBoxDipole(); |
1197 | + | calcBoxDipole_ = true; |
1198 | + | } |
1199 | + | |
1200 | + | } |
1201 | + | |
1202 | void SimInfo::addProperty(GenericData* genData) { | |
1203 | properties_.addProperty(genData); | |
1204 | } | |
# | Line 1177 | Line 1255 | namespace oopse { | |
1255 | Molecule* mol; | |
1256 | ||
1257 | Vector3d comVel(0.0); | |
1258 | < | double totalMass = 0.0; |
1258 | > | RealType totalMass = 0.0; |
1259 | ||
1260 | ||
1261 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1262 | < | double mass = mol->getMass(); |
1262 | > | RealType mass = mol->getMass(); |
1263 | totalMass += mass; | |
1264 | comVel += mass * mol->getComVel(); | |
1265 | } | |
1266 | ||
1267 | #ifdef IS_MPI | |
1268 | < | double tmpMass = totalMass; |
1268 | > | RealType tmpMass = totalMass; |
1269 | Vector3d tmpComVel(comVel); | |
1270 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1271 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1270 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1271 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1272 | #endif | |
1273 | ||
1274 | comVel /= totalMass; | |
# | Line 1203 | Line 1281 | namespace oopse { | |
1281 | Molecule* mol; | |
1282 | ||
1283 | Vector3d com(0.0); | |
1284 | < | double totalMass = 0.0; |
1284 | > | RealType totalMass = 0.0; |
1285 | ||
1286 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1287 | < | double mass = mol->getMass(); |
1287 | > | RealType mass = mol->getMass(); |
1288 | totalMass += mass; | |
1289 | com += mass * mol->getCom(); | |
1290 | } | |
1291 | ||
1292 | #ifdef IS_MPI | |
1293 | < | double tmpMass = totalMass; |
1293 | > | RealType tmpMass = totalMass; |
1294 | Vector3d tmpCom(com); | |
1295 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1296 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1295 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1296 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1297 | #endif | |
1298 | ||
1299 | com /= totalMass; | |
# | Line 1239 | Line 1317 | namespace oopse { | |
1317 | Molecule* mol; | |
1318 | ||
1319 | ||
1320 | < | double totalMass = 0.0; |
1320 | > | RealType totalMass = 0.0; |
1321 | ||
1322 | ||
1323 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1324 | < | double mass = mol->getMass(); |
1324 | > | RealType mass = mol->getMass(); |
1325 | totalMass += mass; | |
1326 | com += mass * mol->getCom(); | |
1327 | comVel += mass * mol->getComVel(); | |
1328 | } | |
1329 | ||
1330 | #ifdef IS_MPI | |
1331 | < | double tmpMass = totalMass; |
1331 | > | RealType tmpMass = totalMass; |
1332 | Vector3d tmpCom(com); | |
1333 | Vector3d tmpComVel(comVel); | |
1334 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1335 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1336 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1334 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1335 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1336 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1337 | #endif | |
1338 | ||
1339 | com /= totalMass; | |
# | Line 1274 | Line 1352 | namespace oopse { | |
1352 | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | |
1353 | ||
1354 | ||
1355 | < | double xx = 0.0; |
1356 | < | double yy = 0.0; |
1357 | < | double zz = 0.0; |
1358 | < | double xy = 0.0; |
1359 | < | double xz = 0.0; |
1360 | < | double yz = 0.0; |
1355 | > | RealType xx = 0.0; |
1356 | > | RealType yy = 0.0; |
1357 | > | RealType zz = 0.0; |
1358 | > | RealType xy = 0.0; |
1359 | > | RealType xz = 0.0; |
1360 | > | RealType yz = 0.0; |
1361 | Vector3d com(0.0); | |
1362 | Vector3d comVel(0.0); | |
1363 | ||
# | Line 1291 | Line 1369 | namespace oopse { | |
1369 | Vector3d thisq(0.0); | |
1370 | Vector3d thisv(0.0); | |
1371 | ||
1372 | < | double thisMass = 0.0; |
1372 | > | RealType thisMass = 0.0; |
1373 | ||
1374 | ||
1375 | ||
# | Line 1329 | Line 1407 | namespace oopse { | |
1407 | #ifdef IS_MPI | |
1408 | Mat3x3d tmpI(inertiaTensor); | |
1409 | Vector3d tmpAngMom; | |
1410 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1411 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1410 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1411 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1412 | #endif | |
1413 | ||
1414 | return; | |
# | Line 1351 | Line 1429 | namespace oopse { | |
1429 | Vector3d thisr(0.0); | |
1430 | Vector3d thisp(0.0); | |
1431 | ||
1432 | < | double thisMass; |
1432 | > | RealType thisMass; |
1433 | ||
1434 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1435 | thisMass = mol->getMass(); | |
# | Line 1364 | Line 1442 | namespace oopse { | |
1442 | ||
1443 | #ifdef IS_MPI | |
1444 | Vector3d tmpAngMom; | |
1445 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1445 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1446 | #endif | |
1447 | ||
1448 | return angularMomentum; | |
1449 | } | |
1450 | ||
1451 | < | |
1451 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { |
1452 | > | return IOIndexToIntegrableObject.at(index); |
1453 | > | } |
1454 | > | |
1455 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1456 | > | IOIndexToIntegrableObject= v; |
1457 | > | } |
1458 | > | |
1459 | > | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1460 | > | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1461 | > | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1462 | > | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1463 | > | */ |
1464 | > | void SimInfo::getGyrationalVolume(RealType &volume){ |
1465 | > | Mat3x3d intTensor; |
1466 | > | RealType det; |
1467 | > | Vector3d dummyAngMom; |
1468 | > | RealType sysconstants; |
1469 | > | RealType geomCnst; |
1470 | > | |
1471 | > | geomCnst = 3.0/2.0; |
1472 | > | /* Get the inertial tensor and angular momentum for free*/ |
1473 | > | getInertiaTensor(intTensor,dummyAngMom); |
1474 | > | |
1475 | > | det = intTensor.determinant(); |
1476 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1477 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1478 | > | return; |
1479 | > | } |
1480 | > | |
1481 | > | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1482 | > | Mat3x3d intTensor; |
1483 | > | Vector3d dummyAngMom; |
1484 | > | RealType sysconstants; |
1485 | > | RealType geomCnst; |
1486 | > | |
1487 | > | geomCnst = 3.0/2.0; |
1488 | > | /* Get the inertial tensor and angular momentum for free*/ |
1489 | > | getInertiaTensor(intTensor,dummyAngMom); |
1490 | > | |
1491 | > | detI = intTensor.determinant(); |
1492 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1493 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1494 | > | return; |
1495 | > | } |
1496 | > | /* |
1497 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1498 | > | assert( v.size() == nAtoms_ + nRigidBodies_); |
1499 | > | sdByGlobalIndex_ = v; |
1500 | > | } |
1501 | > | |
1502 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1503 | > | //assert(index < nAtoms_ + nRigidBodies_); |
1504 | > | return sdByGlobalIndex_.at(index); |
1505 | > | } |
1506 | > | */ |
1507 | }//end namespace oopse | |
1508 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |