# | Line 790 | Line 790 | namespace OpenMD { | |
---|---|---|
790 | set<AtomType*>::iterator i; | |
791 | set<AtomType*> atomTypes; | |
792 | atomTypes = getSimulatedAtomTypes(); | |
793 | < | int usesElectrostatic = 0; |
794 | < | int usesMetallic = 0; |
795 | < | int usesDirectional = 0; |
796 | < | int usesFluctuatingCharges = 0; |
793 | > | bool usesElectrostatic = false; |
794 | > | bool usesMetallic = false; |
795 | > | bool usesDirectional = false; |
796 | > | bool usesFluctuatingCharges = false; |
797 | //loop over all of the atom types | |
798 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
799 | usesElectrostatic |= (*i)->isElectrostatic(); | |
# | Line 801 | Line 801 | namespace OpenMD { | |
801 | usesDirectional |= (*i)->isDirectional(); | |
802 | usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); | |
803 | } | |
804 | < | |
805 | < | #ifdef IS_MPI |
806 | < | int temp; |
804 | > | |
805 | > | #ifdef IS_MPI |
806 | > | bool temp; |
807 | temp = usesDirectional; | |
808 | < | MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
809 | < | |
808 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, |
809 | > | MPI::LOR); |
810 | > | |
811 | temp = usesMetallic; | |
812 | < | MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
812 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, |
813 | > | MPI::LOR); |
814 | ||
815 | temp = usesElectrostatic; | |
816 | < | MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
816 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, |
817 | > | MPI::LOR); |
818 | ||
819 | temp = usesFluctuatingCharges; | |
820 | < | MPI_Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
820 | > | MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, |
821 | > | MPI::LOR); |
822 | #else | |
823 | ||
824 | usesDirectionalAtoms_ = usesDirectional; | |
# | Line 986 | Line 990 | namespace OpenMD { | |
990 | cg->setSnapshotManager(sman_); | |
991 | } | |
992 | } | |
989 | – | |
990 | – | } |
991 | – | |
992 | – | Vector3d SimInfo::getComVel(){ |
993 | – | SimInfo::MoleculeIterator i; |
994 | – | Molecule* mol; |
995 | – | |
996 | – | Vector3d comVel(0.0); |
997 | – | RealType totalMass = 0.0; |
993 | ||
999 | – | |
1000 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1001 | – | RealType mass = mol->getMass(); |
1002 | – | totalMass += mass; |
1003 | – | comVel += mass * mol->getComVel(); |
1004 | – | } |
1005 | – | |
1006 | – | #ifdef IS_MPI |
1007 | – | RealType tmpMass = totalMass; |
1008 | – | Vector3d tmpComVel(comVel); |
1009 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1010 | – | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1011 | – | #endif |
1012 | – | |
1013 | – | comVel /= totalMass; |
1014 | – | |
1015 | – | return comVel; |
994 | } | |
995 | ||
1018 | – | Vector3d SimInfo::getCom(){ |
1019 | – | SimInfo::MoleculeIterator i; |
1020 | – | Molecule* mol; |
996 | ||
1022 | – | Vector3d com(0.0); |
1023 | – | RealType totalMass = 0.0; |
1024 | – | |
1025 | – | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1026 | – | RealType mass = mol->getMass(); |
1027 | – | totalMass += mass; |
1028 | – | com += mass * mol->getCom(); |
1029 | – | } |
1030 | – | |
1031 | – | #ifdef IS_MPI |
1032 | – | RealType tmpMass = totalMass; |
1033 | – | Vector3d tmpCom(com); |
1034 | – | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1035 | – | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1036 | – | #endif |
1037 | – | |
1038 | – | com /= totalMass; |
1039 | – | |
1040 | – | return com; |
1041 | – | |
1042 | – | } |
1043 | – | |
997 | ostream& operator <<(ostream& o, SimInfo& info) { | |
998 | ||
999 | return o; | |
1000 | } | |
1001 | ||
1002 | < | |
1050 | < | /* |
1051 | < | Returns center of mass and center of mass velocity in one function call. |
1052 | < | */ |
1053 | < | |
1054 | < | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
1055 | < | SimInfo::MoleculeIterator i; |
1056 | < | Molecule* mol; |
1057 | < | |
1058 | < | |
1059 | < | RealType totalMass = 0.0; |
1060 | < | |
1061 | < | |
1062 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1063 | < | RealType mass = mol->getMass(); |
1064 | < | totalMass += mass; |
1065 | < | com += mass * mol->getCom(); |
1066 | < | comVel += mass * mol->getComVel(); |
1067 | < | } |
1068 | < | |
1069 | < | #ifdef IS_MPI |
1070 | < | RealType tmpMass = totalMass; |
1071 | < | Vector3d tmpCom(com); |
1072 | < | Vector3d tmpComVel(comVel); |
1073 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1074 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1075 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1076 | < | #endif |
1077 | < | |
1078 | < | com /= totalMass; |
1079 | < | comVel /= totalMass; |
1080 | < | } |
1081 | < | |
1082 | < | /* |
1083 | < | Return intertia tensor for entire system and angular momentum Vector. |
1084 | < | |
1085 | < | |
1086 | < | [ Ixx -Ixy -Ixz ] |
1087 | < | J =| -Iyx Iyy -Iyz | |
1088 | < | [ -Izx -Iyz Izz ] |
1089 | < | */ |
1090 | < | |
1091 | < | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1092 | < | |
1093 | < | |
1094 | < | RealType xx = 0.0; |
1095 | < | RealType yy = 0.0; |
1096 | < | RealType zz = 0.0; |
1097 | < | RealType xy = 0.0; |
1098 | < | RealType xz = 0.0; |
1099 | < | RealType yz = 0.0; |
1100 | < | Vector3d com(0.0); |
1101 | < | Vector3d comVel(0.0); |
1102 | < | |
1103 | < | getComAll(com, comVel); |
1104 | < | |
1105 | < | SimInfo::MoleculeIterator i; |
1106 | < | Molecule* mol; |
1107 | < | |
1108 | < | Vector3d thisq(0.0); |
1109 | < | Vector3d thisv(0.0); |
1110 | < | |
1111 | < | RealType thisMass = 0.0; |
1112 | < | |
1113 | < | |
1114 | < | |
1115 | < | |
1116 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1117 | < | |
1118 | < | thisq = mol->getCom()-com; |
1119 | < | thisv = mol->getComVel()-comVel; |
1120 | < | thisMass = mol->getMass(); |
1121 | < | // Compute moment of intertia coefficients. |
1122 | < | xx += thisq[0]*thisq[0]*thisMass; |
1123 | < | yy += thisq[1]*thisq[1]*thisMass; |
1124 | < | zz += thisq[2]*thisq[2]*thisMass; |
1125 | < | |
1126 | < | // compute products of intertia |
1127 | < | xy += thisq[0]*thisq[1]*thisMass; |
1128 | < | xz += thisq[0]*thisq[2]*thisMass; |
1129 | < | yz += thisq[1]*thisq[2]*thisMass; |
1130 | < | |
1131 | < | angularMomentum += cross( thisq, thisv ) * thisMass; |
1132 | < | |
1133 | < | } |
1134 | < | |
1135 | < | |
1136 | < | inertiaTensor(0,0) = yy + zz; |
1137 | < | inertiaTensor(0,1) = -xy; |
1138 | < | inertiaTensor(0,2) = -xz; |
1139 | < | inertiaTensor(1,0) = -xy; |
1140 | < | inertiaTensor(1,1) = xx + zz; |
1141 | < | inertiaTensor(1,2) = -yz; |
1142 | < | inertiaTensor(2,0) = -xz; |
1143 | < | inertiaTensor(2,1) = -yz; |
1144 | < | inertiaTensor(2,2) = xx + yy; |
1145 | < | |
1146 | < | #ifdef IS_MPI |
1147 | < | Mat3x3d tmpI(inertiaTensor); |
1148 | < | Vector3d tmpAngMom; |
1149 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1150 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1151 | < | #endif |
1152 | < | |
1153 | < | return; |
1154 | < | } |
1155 | < | |
1156 | < | //Returns the angular momentum of the system |
1157 | < | Vector3d SimInfo::getAngularMomentum(){ |
1158 | < | |
1159 | < | Vector3d com(0.0); |
1160 | < | Vector3d comVel(0.0); |
1161 | < | Vector3d angularMomentum(0.0); |
1162 | < | |
1163 | < | getComAll(com,comVel); |
1164 | < | |
1165 | < | SimInfo::MoleculeIterator i; |
1166 | < | Molecule* mol; |
1167 | < | |
1168 | < | Vector3d thisr(0.0); |
1169 | < | Vector3d thisp(0.0); |
1170 | < | |
1171 | < | RealType thisMass; |
1172 | < | |
1173 | < | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1174 | < | thisMass = mol->getMass(); |
1175 | < | thisr = mol->getCom()-com; |
1176 | < | thisp = (mol->getComVel()-comVel)*thisMass; |
1177 | < | |
1178 | < | angularMomentum += cross( thisr, thisp ); |
1179 | < | |
1180 | < | } |
1181 | < | |
1182 | < | #ifdef IS_MPI |
1183 | < | Vector3d tmpAngMom; |
1184 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); |
1185 | < | #endif |
1186 | < | |
1187 | < | return angularMomentum; |
1188 | < | } |
1189 | < | |
1002 | > | |
1003 | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | |
1004 | return IOIndexToIntegrableObject.at(index); | |
1005 | } | |
# | Line 1194 | Line 1007 | namespace OpenMD { | |
1007 | void SimInfo::setIOIndexToIntegrableObject(const vector<StuntDouble*>& v) { | |
1008 | IOIndexToIntegrableObject= v; | |
1009 | } | |
1197 | – | |
1198 | – | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1199 | – | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1200 | – | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1201 | – | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1202 | – | */ |
1203 | – | void SimInfo::getGyrationalVolume(RealType &volume){ |
1204 | – | Mat3x3d intTensor; |
1205 | – | RealType det; |
1206 | – | Vector3d dummyAngMom; |
1207 | – | RealType sysconstants; |
1208 | – | RealType geomCnst; |
1209 | – | |
1210 | – | geomCnst = 3.0/2.0; |
1211 | – | /* Get the inertial tensor and angular momentum for free*/ |
1212 | – | getInertiaTensor(intTensor,dummyAngMom); |
1213 | – | |
1214 | – | det = intTensor.determinant(); |
1215 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1216 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(det); |
1217 | – | return; |
1218 | – | } |
1219 | – | |
1220 | – | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1221 | – | Mat3x3d intTensor; |
1222 | – | Vector3d dummyAngMom; |
1223 | – | RealType sysconstants; |
1224 | – | RealType geomCnst; |
1225 | – | |
1226 | – | geomCnst = 3.0/2.0; |
1227 | – | /* Get the inertial tensor and angular momentum for free*/ |
1228 | – | getInertiaTensor(intTensor,dummyAngMom); |
1229 | – | |
1230 | – | detI = intTensor.determinant(); |
1231 | – | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1232 | – | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(detI); |
1233 | – | return; |
1234 | – | } |
1010 | /* | |
1011 | void SimInfo::setStuntDoubleFromGlobalIndex(vector<StuntDouble*> v) { | |
1012 | assert( v.size() == nAtoms_ + nRigidBodies_); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |