# | 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" | |
# | Line 66 | Line 68 | |
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 89 | Line 92 | namespace oopse { | |
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 | > | useAtomicVirial_(true) { |
97 | ||
98 | MoleculeStamp* molStamp; | |
99 | int nMolWithSameStamp; | |
# | Line 599 | Line 603 | namespace oopse { | |
603 | //setup fortran force field | |
604 | /** @deprecate */ | |
605 | int isError = 0; | |
606 | + | |
607 | + | setupCutoff(); |
608 | ||
609 | setupElectrostaticSummationMethod( isError ); | |
610 | setupSwitchingFunction(); | |
611 | + | setupAccumulateBoxDipole(); |
612 | ||
613 | if(isError){ | |
614 | sprintf( painCave.errMsg, | |
# | Line 609 | Line 616 | namespace oopse { | |
616 | painCave.isFatal = 1; | |
617 | simError(); | |
618 | } | |
612 | – | |
613 | – | |
614 | – | setupCutoff(); |
619 | ||
620 | calcNdf(); | |
621 | calcNdfRaw(); | |
# | Line 661 | Line 665 | namespace oopse { | |
665 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
666 | int useRF; | |
667 | int useSF; | |
668 | + | int useSP; |
669 | + | int useBoxDipole; |
670 | + | |
671 | std::string myMethod; | |
672 | ||
673 | // set the useRF logical | |
674 | useRF = 0; | |
675 | useSF = 0; | |
676 | + | useSP = 0; |
677 | ||
678 | ||
679 | if (simParams_->haveElectrostaticSummationMethod()) { | |
680 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
681 | toUpper(myMethod); | |
682 | < | if (myMethod == "REACTION_FIELD") { |
683 | < | useRF=1; |
684 | < | } else { |
685 | < | if (myMethod == "SHIFTED_FORCE") { |
686 | < | useSF = 1; |
687 | < | } |
682 | > | if (myMethod == "REACTION_FIELD"){ |
683 | > | useRF = 1; |
684 | > | } else if (myMethod == "SHIFTED_FORCE"){ |
685 | > | useSF = 1; |
686 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
687 | > | useSP = 1; |
688 | } | |
689 | } | |
690 | + | |
691 | + | if (simParams_->haveAccumulateBoxDipole()) |
692 | + | if (simParams_->getAccumulateBoxDipole()) |
693 | + | useBoxDipole = 1; |
694 | ||
695 | + | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
696 | + | |
697 | //loop over all of the atom types | |
698 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
699 | useLennardJones |= (*i)->isLennardJones(); | |
# | Line 749 | Line 763 | namespace oopse { | |
763 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
764 | ||
765 | temp = useSF; | |
766 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
766 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
767 | ||
768 | + | temp = useSP; |
769 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
770 | + | |
771 | + | temp = useBoxDipole; |
772 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
773 | + | |
774 | + | temp = useAtomicVirial_; |
775 | + | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
776 | + | |
777 | #endif | |
778 | ||
779 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 768 | Line 791 | namespace oopse { | |
791 | fInfo_.SIM_uses_FLARB = useFLARB; | |
792 | fInfo_.SIM_uses_RF = useRF; | |
793 | fInfo_.SIM_uses_SF = useSF; | |
794 | < | |
795 | < | if( myMethod == "REACTION_FIELD") { |
796 | < | |
774 | < | if (simParams_->haveDielectric()) { |
775 | < | fInfo_.dielect = simParams_->getDielectric(); |
776 | < | } else { |
777 | < | sprintf(painCave.errMsg, |
778 | < | "SimSetup Error: No Dielectric constant was set.\n" |
779 | < | "\tYou are trying to use Reaction Field without" |
780 | < | "\tsetting a dielectric constant!\n"); |
781 | < | painCave.isFatal = 1; |
782 | < | simError(); |
783 | < | } |
784 | < | } |
785 | < | |
794 | > | fInfo_.SIM_uses_SP = useSP; |
795 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
796 | > | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
797 | } | |
798 | ||
799 | void SimInfo::setupFortranSim() { | |
# | Line 799 | Line 810 | namespace oopse { | |
810 | } | |
811 | ||
812 | //calculate mass ratio of cutoff group | |
813 | < | std::vector<double> mfact; |
813 | > | std::vector<RealType> mfact; |
814 | SimInfo::MoleculeIterator mi; | |
815 | Molecule* mol; | |
816 | Molecule::CutoffGroupIterator ci; | |
817 | CutoffGroup* cg; | |
818 | Molecule::AtomIterator ai; | |
819 | Atom* atom; | |
820 | < | double totalMass; |
820 | > | RealType totalMass; |
821 | ||
822 | //to avoid memory reallocation, reserve enough space for mfact | |
823 | mfact.reserve(getNCutoffGroups()); | |
# | Line 867 | Line 878 | namespace oopse { | |
878 | "succesfully sent the simulation information to fortran.\n"); | |
879 | MPIcheckPoint(); | |
880 | #endif // is_mpi | |
881 | + | |
882 | + | // Setup number of neighbors in neighbor list if present |
883 | + | if (simParams_->haveNeighborListNeighbors()) { |
884 | + | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
885 | + | setNeighbors(&nlistNeighbors); |
886 | + | } |
887 | + | |
888 | + | |
889 | } | |
890 | ||
891 | ||
# | Line 936 | Line 955 | namespace oopse { | |
955 | // Check the cutoff policy | |
956 | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | |
957 | ||
958 | + | // Set LJ shifting bools to false |
959 | + | ljsp_ = false; |
960 | + | ljsf_ = false; |
961 | + | |
962 | std::string myPolicy; | |
963 | if (forceFieldOptions_.haveCutoffPolicy()){ | |
964 | myPolicy = forceFieldOptions_.getCutoffPolicy(); | |
# | Line 966 | Line 989 | namespace oopse { | |
989 | notifyFortranCutoffPolicy(&cp); | |
990 | ||
991 | // Check the Skin Thickness for neighborlists | |
992 | < | double skin; |
992 | > | RealType skin; |
993 | if (simParams_->haveSkinThickness()) { | |
994 | skin = simParams_->getSkinThickness(); | |
995 | notifyFortranSkinThickness(&skin); | |
# | Line 999 | Line 1022 | namespace oopse { | |
1022 | simError(); | |
1023 | } | |
1024 | } | |
1025 | + | |
1026 | + | if (simParams_->haveElectrostaticSummationMethod()) { |
1027 | + | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1028 | + | toUpper(myMethod); |
1029 | + | |
1030 | + | if (myMethod == "SHIFTED_POTENTIAL") { |
1031 | + | ljsp_ = true; |
1032 | + | } else if (myMethod == "SHIFTED_FORCE") { |
1033 | + | ljsf_ = true; |
1034 | + | } |
1035 | + | } |
1036 | + | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1037 | ||
1003 | – | notifyFortranCutoffs(&rcut_, &rsw_); |
1004 | – | |
1038 | } else { | |
1039 | ||
1040 | // For electrostatic atoms, we'll assume a large safe value: | |
# | Line 1017 | Line 1050 | namespace oopse { | |
1050 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1051 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
1052 | toUpper(myMethod); | |
1053 | < | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1053 | > | |
1054 | > | // For the time being, we're tethering the LJ shifted behavior to the |
1055 | > | // electrostaticSummationMethod keyword options |
1056 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
1057 | > | ljsp_ = true; |
1058 | > | } else if (myMethod == "SHIFTED_FORCE") { |
1059 | > | ljsf_ = true; |
1060 | > | } |
1061 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1062 | if (simParams_->haveSwitchingRadius()){ | |
1063 | sprintf(painCave.errMsg, | |
1064 | "SimInfo Warning: A value was set for the switchingRadius\n" | |
# | Line 1040 | Line 1081 | namespace oopse { | |
1081 | simError(); | |
1082 | rsw_ = 0.85 * rcut_; | |
1083 | } | |
1084 | < | notifyFortranCutoffs(&rcut_, &rsw_); |
1084 | > | |
1085 | > | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1086 | > | |
1087 | } else { | |
1088 | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | |
1089 | // We'll punt and let fortran figure out the cutoffs later. | |
# | Line 1056 | Line 1099 | namespace oopse { | |
1099 | int errorOut; | |
1100 | int esm = NONE; | |
1101 | int sm = UNDAMPED; | |
1102 | < | double alphaVal; |
1103 | < | double dielectric; |
1104 | < | |
1102 | > | RealType alphaVal; |
1103 | > | RealType dielectric; |
1104 | > | |
1105 | errorOut = isError; | |
1063 | – | alphaVal = simParams_->getDampingAlpha(); |
1064 | – | dielectric = simParams_->getDielectric(); |
1106 | ||
1107 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1108 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
# | Line 1078 | Line 1119 | namespace oopse { | |
1119 | if (myMethod == "SHIFTED_FORCE") { | |
1120 | esm = SHIFTED_FORCE; | |
1121 | } else { | |
1122 | < | if (myMethod == "REACTION_FIELD") { |
1122 | > | if (myMethod == "REACTION_FIELD") { |
1123 | esm = REACTION_FIELD; | |
1124 | + | dielectric = simParams_->getDielectric(); |
1125 | + | if (!simParams_->haveDielectric()) { |
1126 | + | // throw warning |
1127 | + | sprintf( painCave.errMsg, |
1128 | + | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1129 | + | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1130 | + | painCave.isFatal = 0; |
1131 | + | simError(); |
1132 | + | } |
1133 | } else { | |
1134 | // throw error | |
1135 | sprintf( painCave.errMsg, | |
# | Line 1106 | Line 1156 | namespace oopse { | |
1156 | if (myScreen == "DAMPED") { | |
1157 | sm = DAMPED; | |
1158 | if (!simParams_->haveDampingAlpha()) { | |
1159 | < | //throw error |
1159 | > | // first set a cutoff dependent alpha value |
1160 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1161 | > | alphaVal = 0.5125 - rcut_* 0.025; |
1162 | > | // for values rcut > 20.5, alpha is zero |
1163 | > | if (alphaVal < 0) alphaVal = 0; |
1164 | > | |
1165 | > | // throw warning |
1166 | sprintf( painCave.errMsg, | |
1167 | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | |
1168 | < | "\tA default value of %f (1/ang) will be used.\n", alphaVal); |
1168 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1169 | painCave.isFatal = 0; | |
1170 | simError(); | |
1171 | + | } else { |
1172 | + | alphaVal = simParams_->getDampingAlpha(); |
1173 | } | |
1174 | + | |
1175 | } else { | |
1176 | // throw error | |
1177 | sprintf( painCave.errMsg, | |
# | Line 1161 | Line 1220 | namespace oopse { | |
1220 | ||
1221 | } | |
1222 | ||
1223 | + | void SimInfo::setupAccumulateBoxDipole() { |
1224 | + | |
1225 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1226 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1227 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1228 | + | setAccumulateBoxDipole(); |
1229 | + | calcBoxDipole_ = true; |
1230 | + | } |
1231 | + | |
1232 | + | } |
1233 | + | |
1234 | void SimInfo::addProperty(GenericData* genData) { | |
1235 | properties_.addProperty(genData); | |
1236 | } | |
# | Line 1217 | Line 1287 | namespace oopse { | |
1287 | Molecule* mol; | |
1288 | ||
1289 | Vector3d comVel(0.0); | |
1290 | < | double totalMass = 0.0; |
1290 | > | RealType totalMass = 0.0; |
1291 | ||
1292 | ||
1293 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1294 | < | double mass = mol->getMass(); |
1294 | > | RealType mass = mol->getMass(); |
1295 | totalMass += mass; | |
1296 | comVel += mass * mol->getComVel(); | |
1297 | } | |
1298 | ||
1299 | #ifdef IS_MPI | |
1300 | < | double tmpMass = totalMass; |
1300 | > | RealType tmpMass = totalMass; |
1301 | Vector3d tmpComVel(comVel); | |
1302 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1303 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1302 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1303 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1304 | #endif | |
1305 | ||
1306 | comVel /= totalMass; | |
# | Line 1243 | Line 1313 | namespace oopse { | |
1313 | Molecule* mol; | |
1314 | ||
1315 | Vector3d com(0.0); | |
1316 | < | double totalMass = 0.0; |
1316 | > | RealType totalMass = 0.0; |
1317 | ||
1318 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1319 | < | double mass = mol->getMass(); |
1319 | > | RealType mass = mol->getMass(); |
1320 | totalMass += mass; | |
1321 | com += mass * mol->getCom(); | |
1322 | } | |
1323 | ||
1324 | #ifdef IS_MPI | |
1325 | < | double tmpMass = totalMass; |
1325 | > | RealType tmpMass = totalMass; |
1326 | Vector3d tmpCom(com); | |
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); |
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 | #endif | |
1330 | ||
1331 | com /= totalMass; | |
# | Line 1279 | Line 1349 | namespace oopse { | |
1349 | Molecule* mol; | |
1350 | ||
1351 | ||
1352 | < | double totalMass = 0.0; |
1352 | > | RealType totalMass = 0.0; |
1353 | ||
1354 | ||
1355 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1356 | < | double mass = mol->getMass(); |
1356 | > | RealType mass = mol->getMass(); |
1357 | totalMass += mass; | |
1358 | com += mass * mol->getCom(); | |
1359 | comVel += mass * mol->getComVel(); | |
1360 | } | |
1361 | ||
1362 | #ifdef IS_MPI | |
1363 | < | double tmpMass = totalMass; |
1363 | > | RealType tmpMass = totalMass; |
1364 | Vector3d tmpCom(com); | |
1365 | Vector3d tmpComVel(comVel); | |
1366 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1367 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1368 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1366 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1367 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1368 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1369 | #endif | |
1370 | ||
1371 | com /= totalMass; | |
# | Line 1314 | Line 1384 | namespace oopse { | |
1384 | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | |
1385 | ||
1386 | ||
1387 | < | double xx = 0.0; |
1388 | < | double yy = 0.0; |
1389 | < | double zz = 0.0; |
1390 | < | double xy = 0.0; |
1391 | < | double xz = 0.0; |
1392 | < | double yz = 0.0; |
1387 | > | RealType xx = 0.0; |
1388 | > | RealType yy = 0.0; |
1389 | > | RealType zz = 0.0; |
1390 | > | RealType xy = 0.0; |
1391 | > | RealType xz = 0.0; |
1392 | > | RealType yz = 0.0; |
1393 | Vector3d com(0.0); | |
1394 | Vector3d comVel(0.0); | |
1395 | ||
# | Line 1331 | Line 1401 | namespace oopse { | |
1401 | Vector3d thisq(0.0); | |
1402 | Vector3d thisv(0.0); | |
1403 | ||
1404 | < | double thisMass = 0.0; |
1404 | > | RealType thisMass = 0.0; |
1405 | ||
1406 | ||
1407 | ||
# | Line 1369 | Line 1439 | namespace oopse { | |
1439 | #ifdef IS_MPI | |
1440 | Mat3x3d tmpI(inertiaTensor); | |
1441 | Vector3d tmpAngMom; | |
1442 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1443 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1442 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1443 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1444 | #endif | |
1445 | ||
1446 | return; | |
# | Line 1391 | Line 1461 | namespace oopse { | |
1461 | Vector3d thisr(0.0); | |
1462 | Vector3d thisp(0.0); | |
1463 | ||
1464 | < | double thisMass; |
1464 | > | RealType thisMass; |
1465 | ||
1466 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
1467 | thisMass = mol->getMass(); | |
# | Line 1404 | Line 1474 | namespace oopse { | |
1474 | ||
1475 | #ifdef IS_MPI | |
1476 | Vector3d tmpAngMom; | |
1477 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1477 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1478 | #endif | |
1479 | ||
1480 | return angularMomentum; | |
1481 | } | |
1482 | ||
1483 | < | |
1483 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { |
1484 | > | return IOIndexToIntegrableObject.at(index); |
1485 | > | } |
1486 | > | |
1487 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1488 | > | IOIndexToIntegrableObject= v; |
1489 | > | } |
1490 | > | |
1491 | > | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1492 | > | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1493 | > | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1494 | > | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1495 | > | */ |
1496 | > | void SimInfo::getGyrationalVolume(RealType &volume){ |
1497 | > | Mat3x3d intTensor; |
1498 | > | RealType det; |
1499 | > | Vector3d dummyAngMom; |
1500 | > | RealType sysconstants; |
1501 | > | RealType geomCnst; |
1502 | > | |
1503 | > | geomCnst = 3.0/2.0; |
1504 | > | /* Get the inertial tensor and angular momentum for free*/ |
1505 | > | getInertiaTensor(intTensor,dummyAngMom); |
1506 | > | |
1507 | > | det = intTensor.determinant(); |
1508 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1509 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1510 | > | return; |
1511 | > | } |
1512 | > | |
1513 | > | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1514 | > | Mat3x3d intTensor; |
1515 | > | Vector3d dummyAngMom; |
1516 | > | RealType sysconstants; |
1517 | > | RealType geomCnst; |
1518 | > | |
1519 | > | geomCnst = 3.0/2.0; |
1520 | > | /* Get the inertial tensor and angular momentum for free*/ |
1521 | > | getInertiaTensor(intTensor,dummyAngMom); |
1522 | > | |
1523 | > | detI = intTensor.determinant(); |
1524 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1525 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1526 | > | return; |
1527 | > | } |
1528 | > | /* |
1529 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1530 | > | assert( v.size() == nAtoms_ + nRigidBodies_); |
1531 | > | sdByGlobalIndex_ = v; |
1532 | > | } |
1533 | > | |
1534 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1535 | > | //assert(index < nAtoms_ + nRigidBodies_); |
1536 | > | return sdByGlobalIndex_.at(index); |
1537 | > | } |
1538 | > | */ |
1539 | }//end namespace oopse | |
1540 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |