# | Line 36 | Line 36 | |
---|---|---|
36 | * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | |
37 | * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | |
38 | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). | |
39 | < | * [4] Vardeman & Gezelter, in progress (2009). |
39 | > | * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 | > | * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 | */ | |
42 | ||
43 | /** | |
# | Line 58 | Line 59 | |
59 | #include "utils/simError.h" | |
60 | #include "selection/SelectionManager.hpp" | |
61 | #include "io/ForceFieldOptions.hpp" | |
62 | < | #include "UseTheForce/ForceField.hpp" |
62 | > | #include "brains/ForceField.hpp" |
63 | #include "nonbonded/SwitchingFunction.hpp" | |
64 | #ifdef IS_MPI | |
65 | #include <mpi.h> | |
# | Line 71 | Line 72 | namespace OpenMD { | |
72 | forceField_(ff), simParams_(simParams), | |
73 | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | |
74 | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | |
75 | < | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), |
75 | > | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0), |
76 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), | |
77 | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | |
78 | < | nConstraints_(0), sman_(NULL), topologyDone_(false), |
78 | > | nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), |
79 | calcBoxDipole_(false), useAtomicVirial_(true) { | |
80 | ||
81 | MoleculeStamp* molStamp; | |
# | Line 87 | Line 88 | namespace OpenMD { | |
88 | ||
89 | vector<Component*> components = simParams->getComponents(); | |
90 | ||
91 | < | for (vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { |
91 | > | for (vector<Component*>::iterator i = components.begin(); |
92 | > | i !=components.end(); ++i) { |
93 | molStamp = (*i)->getMoleculeStamp(); | |
94 | nMolWithSameStamp = (*i)->getNMol(); | |
95 | ||
# | Line 224 | Line 226 | namespace OpenMD { | |
226 | ||
227 | ||
228 | void SimInfo::calcNdf() { | |
229 | < | int ndf_local; |
229 | > | int ndf_local, nfq_local; |
230 | MoleculeIterator i; | |
231 | vector<StuntDouble*>::iterator j; | |
232 | + | vector<Atom*>::iterator k; |
233 | + | |
234 | Molecule* mol; | |
235 | < | StuntDouble* integrableObject; |
235 | > | StuntDouble* sd; |
236 | > | Atom* atom; |
237 | ||
238 | ndf_local = 0; | |
239 | + | nfq_local = 0; |
240 | ||
241 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
236 | – | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
237 | – | integrableObject = mol->nextIntegrableObject(j)) { |
242 | ||
243 | + | for (sd = mol->beginIntegrableObject(j); sd != NULL; |
244 | + | sd = mol->nextIntegrableObject(j)) { |
245 | + | |
246 | ndf_local += 3; | |
247 | ||
248 | < | if (integrableObject->isDirectional()) { |
249 | < | if (integrableObject->isLinear()) { |
248 | > | if (sd->isDirectional()) { |
249 | > | if (sd->isLinear()) { |
250 | ndf_local += 2; | |
251 | } else { | |
252 | ndf_local += 3; | |
253 | } | |
254 | } | |
248 | – | |
255 | } | |
256 | + | |
257 | + | for (atom = mol->beginFluctuatingCharge(k); atom != NULL; |
258 | + | atom = mol->nextFluctuatingCharge(k)) { |
259 | + | if (atom->isFluctuatingCharge()) { |
260 | + | nfq_local++; |
261 | + | } |
262 | + | } |
263 | } | |
264 | ||
265 | + | ndfLocal_ = ndf_local; |
266 | + | |
267 | // n_constraints is local, so subtract them on each processor | |
268 | ndf_local -= nConstraints_; | |
269 | ||
270 | #ifdef IS_MPI | |
271 | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | |
272 | + | MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); |
273 | #else | |
274 | ndf_ = ndf_local; | |
275 | + | nGlobalFluctuatingCharges_ = nfq_local; |
276 | #endif | |
277 | ||
278 | // nZconstraints_ is global, as are the 3 COM translations for the | |
# | Line 298 | Line 315 | namespace OpenMD { | |
315 | MoleculeIterator i; | |
316 | vector<StuntDouble*>::iterator j; | |
317 | Molecule* mol; | |
318 | < | StuntDouble* integrableObject; |
318 | > | StuntDouble* sd; |
319 | ||
320 | // Raw degrees of freedom that we have to set | |
321 | ndfRaw_local = 0; | |
322 | ||
323 | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | |
307 | – | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
308 | – | integrableObject = mol->nextIntegrableObject(j)) { |
324 | ||
325 | + | for (sd = mol->beginIntegrableObject(j); sd != NULL; |
326 | + | sd = mol->nextIntegrableObject(j)) { |
327 | + | |
328 | ndfRaw_local += 3; | |
329 | ||
330 | < | if (integrableObject->isDirectional()) { |
331 | < | if (integrableObject->isLinear()) { |
330 | > | if (sd->isDirectional()) { |
331 | > | if (sd->isLinear()) { |
332 | ndfRaw_local += 2; | |
333 | } else { | |
334 | ndfRaw_local += 3; | |
# | Line 370 | Line 388 | namespace OpenMD { | |
388 | Molecule::RigidBodyIterator rbIter; | |
389 | RigidBody* rb; | |
390 | Molecule::IntegrableObjectIterator ii; | |
391 | < | StuntDouble* integrableObject; |
391 | > | StuntDouble* sd; |
392 | ||
393 | < | for (integrableObject = mol->beginIntegrableObject(ii); |
394 | < | integrableObject != NULL; |
377 | < | integrableObject = mol->nextIntegrableObject(ii)) { |
393 | > | for (sd = mol->beginIntegrableObject(ii); sd != NULL; |
394 | > | sd = mol->nextIntegrableObject(ii)) { |
395 | ||
396 | < | if (integrableObject->isRigidBody()) { |
397 | < | rb = static_cast<RigidBody*>(integrableObject); |
396 | > | if (sd->isRigidBody()) { |
397 | > | rb = static_cast<RigidBody*>(sd); |
398 | vector<Atom*> atoms = rb->getAtoms(); | |
399 | set<int> rigidAtoms; | |
400 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
# | Line 388 | Line 405 | namespace OpenMD { | |
405 | } | |
406 | } else { | |
407 | set<int> oneAtomSet; | |
408 | < | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
409 | < | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
408 | > | oneAtomSet.insert(sd->getGlobalIndex()); |
409 | > | atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet)); |
410 | } | |
411 | } | |
412 | ||
# | Line 523 | Line 540 | namespace OpenMD { | |
540 | Molecule::RigidBodyIterator rbIter; | |
541 | RigidBody* rb; | |
542 | Molecule::IntegrableObjectIterator ii; | |
543 | < | StuntDouble* integrableObject; |
543 | > | StuntDouble* sd; |
544 | ||
545 | < | for (integrableObject = mol->beginIntegrableObject(ii); |
546 | < | integrableObject != NULL; |
530 | < | integrableObject = mol->nextIntegrableObject(ii)) { |
545 | > | for (sd = mol->beginIntegrableObject(ii); sd != NULL; |
546 | > | sd = mol->nextIntegrableObject(ii)) { |
547 | ||
548 | < | if (integrableObject->isRigidBody()) { |
549 | < | rb = static_cast<RigidBody*>(integrableObject); |
548 | > | if (sd->isRigidBody()) { |
549 | > | rb = static_cast<RigidBody*>(sd); |
550 | vector<Atom*> atoms = rb->getAtoms(); | |
551 | set<int> rigidAtoms; | |
552 | for (int i = 0; i < static_cast<int>(atoms.size()); ++i) { | |
# | Line 541 | Line 557 | namespace OpenMD { | |
557 | } | |
558 | } else { | |
559 | set<int> oneAtomSet; | |
560 | < | oneAtomSet.insert(integrableObject->getGlobalIndex()); |
561 | < | atomGroups.insert(map<int, set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); |
560 | > | oneAtomSet.insert(sd->getGlobalIndex()); |
561 | > | atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet)); |
562 | } | |
563 | } | |
564 | ||
# | Line 766 | Line 782 | namespace OpenMD { | |
782 | ||
783 | void SimInfo::setupSimVariables() { | |
784 | useAtomicVirial_ = simParams_->getUseAtomicVirial(); | |
785 | < | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
785 | > | // we only call setAccumulateBoxDipole if the accumulateBoxDipole |
786 | > | // parameter is true |
787 | calcBoxDipole_ = false; | |
788 | if ( simParams_->haveAccumulateBoxDipole() ) | |
789 | if ( simParams_->getAccumulateBoxDipole() ) { | |
# | Line 776 | Line 793 | namespace OpenMD { | |
793 | set<AtomType*>::iterator i; | |
794 | set<AtomType*> atomTypes; | |
795 | atomTypes = getSimulatedAtomTypes(); | |
796 | < | int usesElectrostatic = 0; |
797 | < | int usesMetallic = 0; |
798 | < | int usesDirectional = 0; |
796 | > | bool usesElectrostatic = false; |
797 | > | bool usesMetallic = false; |
798 | > | bool usesDirectional = false; |
799 | > | bool usesFluctuatingCharges = false; |
800 | //loop over all of the atom types | |
801 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
802 | usesElectrostatic |= (*i)->isElectrostatic(); | |
803 | usesMetallic |= (*i)->isMetal(); | |
804 | usesDirectional |= (*i)->isDirectional(); | |
805 | < | } |
806 | < | |
807 | < | #ifdef IS_MPI |
808 | < | int temp; |
805 | > | usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); |
806 | > | } |
807 | > | |
808 | > | #ifdef IS_MPI |
809 | > | bool temp; |
810 | temp = usesDirectional; | |
811 | < | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
812 | < | |
811 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, |
812 | > | MPI::LOR); |
813 | > | |
814 | temp = usesMetallic; | |
815 | < | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
815 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, |
816 | > | MPI::LOR); |
817 | ||
818 | temp = usesElectrostatic; | |
819 | < | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
819 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, |
820 | > | MPI::LOR); |
821 | > | |
822 | > | temp = usesFluctuatingCharges; |
823 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, |
824 | > | MPI::LOR); |
825 | #else | |
826 | ||
827 | usesDirectionalAtoms_ = usesDirectional; | |
828 | usesMetallicAtoms_ = usesMetallic; | |
829 | usesElectrostaticAtoms_ = usesElectrostatic; | |
830 | + | usesFluctuatingCharges_ = usesFluctuatingCharges; |
831 | ||
832 | #endif | |
833 | ||
# | Line 954 | Line 981 | namespace OpenMD { | |
981 | ||
982 | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | |
983 | ||
984 | < | for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { |
984 | > | for (atom = mol->beginAtom(atomIter); atom != NULL; |
985 | > | atom = mol->nextAtom(atomIter)) { |
986 | atom->setSnapshotManager(sman_); | |
987 | } | |
988 | ||
989 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
989 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
990 | > | rb = mol->nextRigidBody(rbIter)) { |
991 | rb->setSnapshotManager(sman_); | |
992 | } | |
993 | ||
994 | < | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { |
994 | > | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; |
995 | > | cg = mol->nextCutoffGroup(cgIter)) { |
996 | cg->setSnapshotManager(sman_); | |
997 | } | |
998 | } | |
999 | ||
1000 | } | |
1001 | ||
972 | – | Vector3d SimInfo::getComVel(){ |
973 | – | SimInfo::MoleculeIterator i; |
974 | – | Molecule* mol; |
1002 | ||
976 | – | Vector3d comVel(0.0); |
977 | – | RealType totalMass = 0.0; |
978 | – | |
979 | – | |
980 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
981 | – | RealType mass = mol->getMass(); |
982 | – | totalMass += mass; |
983 | – | comVel += mass * mol->getComVel(); |
984 | – | } |
985 | – | |
986 | – | #ifdef IS_MPI |
987 | – | RealType tmpMass = totalMass; |
988 | – | Vector3d tmpComVel(comVel); |
989 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
990 | – | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
991 | – | #endif |
992 | – | |
993 | – | comVel /= totalMass; |
994 | – | |
995 | – | return comVel; |
996 | – | } |
997 | – | |
998 | – | Vector3d SimInfo::getCom(){ |
999 | – | SimInfo::MoleculeIterator i; |
1000 | – | Molecule* mol; |
1001 | – | |
1002 | – | Vector3d com(0.0); |
1003 | – | RealType totalMass = 0.0; |
1004 | – | |
1005 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1006 | – | RealType mass = mol->getMass(); |
1007 | – | totalMass += mass; |
1008 | – | com += mass * mol->getCom(); |
1009 | – | } |
1010 | – | |
1011 | – | #ifdef IS_MPI |
1012 | – | RealType tmpMass = totalMass; |
1013 | – | Vector3d tmpCom(com); |
1014 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1015 | – | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1016 | – | #endif |
1017 | – | |
1018 | – | com /= totalMass; |
1019 | – | |
1020 | – | return com; |
1021 | – | |
1022 | – | } |
1023 | – | |
1003 | ostream& operator <<(ostream& o, SimInfo& info) { | |
1004 | ||
1005 | return o; | |
1006 | } | |
1007 | ||
1008 | < | |
1030 | < | /* |
1031 | < | Returns center of mass and center of mass velocity in one function call. |
1032 | < | */ |
1033 | < | |
1034 | < | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
1035 | < | SimInfo::MoleculeIterator i; |
1036 | < | Molecule* mol; |
1037 | < | |
1038 | < | |
1039 | < | RealType totalMass = 0.0; |
1040 | < | |
1041 | < | |
1042 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1043 | < | RealType mass = mol->getMass(); |
1044 | < | totalMass += mass; |
1045 | < | com += mass * mol->getCom(); |
1046 | < | comVel += mass * mol->getComVel(); |
1047 | < | } |
1048 | < | |
1049 | < | #ifdef IS_MPI |
1050 | < | RealType tmpMass = totalMass; |
1051 | < | Vector3d tmpCom(com); |
1052 | < | Vector3d tmpComVel(comVel); |
1053 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1054 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1055 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1056 | < | #endif |
1057 | < | |
1058 | < | com /= totalMass; |
1059 | < | comVel /= totalMass; |
1060 | < | } |
1061 | < | |
1062 | < | /* |
1063 | < | Return intertia tensor for entire system and angular momentum Vector. |
1064 | < | |
1065 | < | |
1066 | < | [ Ixx -Ixy -Ixz ] |
1067 | < | J =| -Iyx Iyy -Iyz | |
1068 | < | [ -Izx -Iyz Izz ] |
1069 | < | */ |
1070 | < | |
1071 | < | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1072 | < | |
1073 | < | |
1074 | < | RealType xx = 0.0; |
1075 | < | RealType yy = 0.0; |
1076 | < | RealType zz = 0.0; |
1077 | < | RealType xy = 0.0; |
1078 | < | RealType xz = 0.0; |
1079 | < | RealType yz = 0.0; |
1080 | < | Vector3d com(0.0); |
1081 | < | Vector3d comVel(0.0); |
1082 | < | |
1083 | < | getComAll(com, comVel); |
1084 | < | |
1085 | < | SimInfo::MoleculeIterator i; |
1086 | < | Molecule* mol; |
1087 | < | |
1088 | < | Vector3d thisq(0.0); |
1089 | < | Vector3d thisv(0.0); |
1090 | < | |
1091 | < | RealType thisMass = 0.0; |
1092 | < | |
1093 | < | |
1094 | < | |
1095 | < | |
1096 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1097 | < | |
1098 | < | thisq = mol->getCom()-com; |
1099 | < | thisv = mol->getComVel()-comVel; |
1100 | < | thisMass = mol->getMass(); |
1101 | < | // Compute moment of intertia coefficients. |
1102 | < | xx += thisq[0]*thisq[0]*thisMass; |
1103 | < | yy += thisq[1]*thisq[1]*thisMass; |
1104 | < | zz += thisq[2]*thisq[2]*thisMass; |
1105 | < | |
1106 | < | // compute products of intertia |
1107 | < | xy += thisq[0]*thisq[1]*thisMass; |
1108 | < | xz += thisq[0]*thisq[2]*thisMass; |
1109 | < | yz += thisq[1]*thisq[2]*thisMass; |
1110 | < | |
1111 | < | angularMomentum += cross( thisq, thisv ) * thisMass; |
1112 | < | |
1113 | < | } |
1114 | < | |
1115 | < | |
1116 | < | inertiaTensor(0,0) = yy + zz; |
1117 | < | inertiaTensor(0,1) = -xy; |
1118 | < | inertiaTensor(0,2) = -xz; |
1119 | < | inertiaTensor(1,0) = -xy; |
1120 | < | inertiaTensor(1,1) = xx + zz; |
1121 | < | inertiaTensor(1,2) = -yz; |
1122 | < | inertiaTensor(2,0) = -xz; |
1123 | < | inertiaTensor(2,1) = -yz; |
1124 | < | inertiaTensor(2,2) = xx + yy; |
1125 | < | |
1126 | < | #ifdef IS_MPI |
1127 | < | Mat3x3d tmpI(inertiaTensor); |
1128 | < | Vector3d tmpAngMom; |
1129 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1130 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1131 | < | #endif |
1132 | < | |
1133 | < | return; |
1134 | < | } |
1135 | < | |
1136 | < | //Returns the angular momentum of the system |
1137 | < | Vector3d SimInfo::getAngularMomentum(){ |
1138 | < | |
1139 | < | Vector3d com(0.0); |
1140 | < | Vector3d comVel(0.0); |
1141 | < | Vector3d angularMomentum(0.0); |
1142 | < | |
1143 | < | getComAll(com,comVel); |
1144 | < | |
1145 | < | SimInfo::MoleculeIterator i; |
1146 | < | Molecule* mol; |
1147 | < | |
1148 | < | Vector3d thisr(0.0); |
1149 | < | Vector3d thisp(0.0); |
1150 | < | |
1151 | < | RealType thisMass; |
1152 | < | |
1153 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1154 | < | thisMass = mol->getMass(); |
1155 | < | thisr = mol->getCom()-com; |
1156 | < | thisp = (mol->getComVel()-comVel)*thisMass; |
1157 | < | |
1158 | < | angularMomentum += cross( thisr, thisp ); |
1159 | < | |
1160 | < | } |
1161 | < | |
1162 | < | #ifdef IS_MPI |
1163 | < | Vector3d tmpAngMom; |
1164 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1165 | < | #endif |
1166 | < | |
1167 | < | return angularMomentum; |
1168 | < | } |
1169 | < | |
1008 | > | |
1009 | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | |
1010 | < | return IOIndexToIntegrableObject.at(index); |
1010 | > | if (index >= IOIndexToIntegrableObject.size()) { |
1011 | > | sprintf(painCave.errMsg, |
1012 | > | "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n" |
1013 | > | "\tindex exceeds number of known objects!\n"); |
1014 | > | painCave.isFatal = 1; |
1015 | > | simError(); |
1016 | > | return NULL; |
1017 | > | } else |
1018 | > | return IOIndexToIntegrableObject.at(index); |
1019 | } | |
1020 | ||
1021 | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { | |
1022 | IOIndexToIntegrableObject= v; | |
1023 | } | |
1024 | ||
1178 | – | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1179 | – | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1180 | – | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1181 | – | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1182 | – | */ |
1183 | – | void SimInfo::getGyrationalVolume(RealType &volume){ |
1184 | – | Mat3x3d intTensor; |
1185 | – | RealType det; |
1186 | – | Vector3d dummyAngMom; |
1187 | – | RealType sysconstants; |
1188 | – | RealType geomCnst; |
1189 | – | |
1190 | – | geomCnst = 3.0/2.0; |
1191 | – | /* Get the inertial tensor and angular momentum for free*/ |
1192 | – | getInertiaTensor(intTensor,dummyAngMom); |
1193 | – | |
1194 | – | det = intTensor.determinant(); |
1195 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1196 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1197 | – | return; |
1198 | – | } |
1199 | – | |
1200 | – | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1201 | – | Mat3x3d intTensor; |
1202 | – | Vector3d dummyAngMom; |
1203 | – | RealType sysconstants; |
1204 | – | RealType geomCnst; |
1205 | – | |
1206 | – | geomCnst = 3.0/2.0; |
1207 | – | /* Get the inertial tensor and angular momentum for free*/ |
1208 | – | getInertiaTensor(intTensor,dummyAngMom); |
1209 | – | |
1210 | – | detI = intTensor.determinant(); |
1211 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1212 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1213 | – | return; |
1214 | – | } |
1215 | – | /* |
1216 | – | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { |
1217 | – | assert( v.size() == nAtoms_ + nRigidBodies_); |
1218 | – | sdByGlobalIndex_ = v; |
1219 | – | } |
1220 | – | |
1221 | – | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1222 | – | //assert(index < nAtoms_ + nRigidBodies_); |
1223 | – | return sdByGlobalIndex_.at(index); |
1224 | – | } |
1225 | – | */ |
1025 | int SimInfo::getNGlobalConstraints() { | |
1026 | int nGlobalConstraints; | |
1027 | #ifdef IS_MPI |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |