| 59 |  | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" | 
| 60 |  | #include "UseTheForce/doForces_interface.h" | 
| 61 |  | #include "UseTheForce/DarkSide/electrostatic_interface.h" | 
| 62 | – | #include "UseTheForce/notifyCutoffs_interface.h" | 
| 62 |  | #include "UseTheForce/DarkSide/switcheroo_interface.h" | 
| 63 |  | #include "utils/MemoryUtils.hpp" | 
| 64 |  | #include "utils/simError.h" | 
| 65 |  | #include "selection/SelectionManager.hpp" | 
| 66 | + | #include "io/ForceFieldOptions.hpp" | 
| 67 | + | #include "UseTheForce/ForceField.hpp" | 
| 68 |  |  | 
| 69 |  | #ifdef IS_MPI | 
| 70 |  | #include "UseTheForce/mpiComponentPlan.h" | 
| 82 |  | return result; | 
| 83 |  | } | 
| 84 |  |  | 
| 85 | < | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, | 
| 86 | < | ForceField* ff, Globals* simParams) : | 
| 87 | < | stamps_(stamps), forceField_(ff), simParams_(simParams), | 
| 87 | < | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 85 | > | SimInfo::SimInfo(ForceField* ff, Globals* simParams) : | 
| 86 | > | forceField_(ff), simParams_(simParams), | 
| 87 | > | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 88 |  | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | 
| 89 |  | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 90 |  | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0), | 
| 91 |  | nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0), | 
| 92 |  | sman_(NULL), fortranInitialized_(false) { | 
| 93 |  |  | 
| 94 | – |  | 
| 95 | – | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; | 
| 94 |  | MoleculeStamp* molStamp; | 
| 95 |  | int nMolWithSameStamp; | 
| 96 |  | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | 
| 98 |  | CutoffGroupStamp* cgStamp; | 
| 99 |  | RigidBodyStamp* rbStamp; | 
| 100 |  | int nRigidAtoms = 0; | 
| 101 | < |  | 
| 102 | < | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { | 
| 103 | < | molStamp = i->first; | 
| 104 | < | nMolWithSameStamp = i->second; | 
| 101 | > | std::vector<Component*> components = simParams->getComponents(); | 
| 102 | > |  | 
| 103 | > | for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) { | 
| 104 | > | molStamp = (*i)->getMoleculeStamp(); | 
| 105 | > | nMolWithSameStamp = (*i)->getNMol(); | 
| 106 |  |  | 
| 107 |  | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 108 |  |  | 
| 109 |  | //calculate atoms in molecules | 
| 110 |  | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 111 |  |  | 
| 113 | – |  | 
| 112 |  | //calculate atoms in cutoff groups | 
| 113 |  | int nAtomsInGroups = 0; | 
| 114 |  | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); | 
| 115 |  |  | 
| 116 |  | for (int j=0; j < nCutoffGroupsInStamp; j++) { | 
| 117 | < | cgStamp = molStamp->getCutoffGroup(j); | 
| 117 | > | cgStamp = molStamp->getCutoffGroupStamp(j); | 
| 118 |  | nAtomsInGroups += cgStamp->getNMembers(); | 
| 119 |  | } | 
| 120 |  |  | 
| 127 |  | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); | 
| 128 |  |  | 
| 129 |  | for (int j=0; j < nRigidBodiesInStamp; j++) { | 
| 130 | < | rbStamp = molStamp->getRigidBody(j); | 
| 130 | > | rbStamp = molStamp->getRigidBodyStamp(j); | 
| 131 |  | nAtomsInRigidBodies += rbStamp->getNMembers(); | 
| 132 |  | } | 
| 133 |  |  | 
| 166 |  | } | 
| 167 |  | molecules_.clear(); | 
| 168 |  |  | 
| 171 | – | delete stamps_; | 
| 169 |  | delete sman_; | 
| 170 |  | delete simParams_; | 
| 171 |  | delete forceField_; | 
| 272 |  | } | 
| 273 |  | } | 
| 274 |  |  | 
| 275 | < | }//end for (integrableObject) | 
| 276 | < | }// end for (mol) | 
| 275 | > | } | 
| 276 | > | } | 
| 277 |  |  | 
| 278 |  | // n_constraints is local, so subtract them on each processor | 
| 279 |  | ndf_local -= nConstraints_; | 
| 290 |  |  | 
| 291 |  | } | 
| 292 |  |  | 
| 293 | + | int SimInfo::getFdf() { | 
| 294 | + | #ifdef IS_MPI | 
| 295 | + | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 296 | + | #else | 
| 297 | + | fdf_ = fdf_local; | 
| 298 | + | #endif | 
| 299 | + | return fdf_; | 
| 300 | + | } | 
| 301 | + |  | 
| 302 |  | void SimInfo::calcNdfRaw() { | 
| 303 |  | int ndfRaw_local; | 
| 304 |  |  | 
| 799 |  | } | 
| 800 |  |  | 
| 801 |  | //calculate mass ratio of cutoff group | 
| 802 | < | std::vector<double> mfact; | 
| 802 | > | std::vector<RealType> mfact; | 
| 803 |  | SimInfo::MoleculeIterator mi; | 
| 804 |  | Molecule* mol; | 
| 805 |  | Molecule::CutoffGroupIterator ci; | 
| 806 |  | CutoffGroup* cg; | 
| 807 |  | Molecule::AtomIterator ai; | 
| 808 |  | Atom* atom; | 
| 809 | < | double totalMass; | 
| 809 | > | RealType totalMass; | 
| 810 |  |  | 
| 811 |  | //to avoid memory reallocation, reserve enough space for mfact | 
| 812 |  | mfact.reserve(getNCutoffGroups()); | 
| 929 |  |  | 
| 930 |  | #endif | 
| 931 |  |  | 
| 932 | < | double SimInfo::calcMaxCutoffRadius() { | 
| 932 | > | void SimInfo::setupCutoff() { | 
| 933 | > |  | 
| 934 | > | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | 
| 935 |  |  | 
| 936 | + | // Check the cutoff policy | 
| 937 | + | int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | 
| 938 |  |  | 
| 939 | < | std::set<AtomType*> atomTypes; | 
| 940 | < | std::set<AtomType*>::iterator i; | 
| 941 | < | std::vector<double> cutoffRadius; | 
| 942 | < |  | 
| 943 | < | //get the unique atom types | 
| 934 | < | atomTypes = getUniqueAtomTypes(); | 
| 935 | < |  | 
| 936 | < | //query the max cutoff radius among these atom types | 
| 937 | < | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 938 | < | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); | 
| 939 | > | std::string myPolicy; | 
| 940 | > | if (forceFieldOptions_.haveCutoffPolicy()){ | 
| 941 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); | 
| 942 | > | }else if (simParams_->haveCutoffPolicy()) { | 
| 943 | > | myPolicy = simParams_->getCutoffPolicy(); | 
| 944 |  | } | 
| 945 |  |  | 
| 946 | < | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); | 
| 942 | < | #ifdef IS_MPI | 
| 943 | < | //pick the max cutoff radius among the processors | 
| 944 | < | #endif | 
| 945 | < |  | 
| 946 | < | return maxCutoffRadius; | 
| 947 | < | } | 
| 948 | < |  | 
| 949 | < | void SimInfo::getCutoff(double& rcut, double& rsw) { | 
| 950 | < |  | 
| 951 | < | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 952 | < |  | 
| 953 | < | if (!simParams_->haveCutoffRadius()){ | 
| 954 | < | sprintf(painCave.errMsg, | 
| 955 | < | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 956 | < | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 957 | < | "\tfor the cutoffRadius.\n"); | 
| 958 | < | painCave.isFatal = 0; | 
| 959 | < | simError(); | 
| 960 | < | rcut = 15.0; | 
| 961 | < | } else{ | 
| 962 | < | rcut = simParams_->getCutoffRadius(); | 
| 963 | < | } | 
| 964 | < |  | 
| 965 | < | if (!simParams_->haveSwitchingRadius()){ | 
| 966 | < | sprintf(painCave.errMsg, | 
| 967 | < | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 968 | < | "\tOOPSE will use a default value of\n" | 
| 969 | < | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 970 | < | painCave.isFatal = 0; | 
| 971 | < | simError(); | 
| 972 | < | rsw = 0.85 * rcut; | 
| 973 | < | } else{ | 
| 974 | < | rsw = simParams_->getSwitchingRadius(); | 
| 975 | < | } | 
| 976 | < |  | 
| 977 | < | } else { | 
| 978 | < | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | 
| 979 | < | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | 
| 980 | < |  | 
| 981 | < | if (simParams_->haveCutoffRadius()) { | 
| 982 | < | rcut = simParams_->getCutoffRadius(); | 
| 983 | < | } else { | 
| 984 | < | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | 
| 985 | < | rcut = calcMaxCutoffRadius(); | 
| 986 | < | } | 
| 987 | < |  | 
| 988 | < | if (simParams_->haveSwitchingRadius()) { | 
| 989 | < | rsw  = simParams_->getSwitchingRadius(); | 
| 990 | < | } else { | 
| 991 | < | rsw = rcut; | 
| 992 | < | } | 
| 993 | < |  | 
| 994 | < | } | 
| 995 | < | } | 
| 996 | < |  | 
| 997 | < | void SimInfo::setupCutoff() { | 
| 998 | < | getCutoff(rcut_, rsw_); | 
| 999 | < | double rnblist = rcut_ + 1; // skin of neighbor list | 
| 1000 | < |  | 
| 1001 | < | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | 
| 1002 | < |  | 
| 1003 | < | int cp =  TRADITIONAL_CUTOFF_POLICY; | 
| 1004 | < | if (simParams_->haveCutoffPolicy()) { | 
| 1005 | < | std::string myPolicy = simParams_->getCutoffPolicy(); | 
| 946 | > | if (!myPolicy.empty()){ | 
| 947 |  | toUpper(myPolicy); | 
| 948 |  | if (myPolicy == "MIX") { | 
| 949 |  | cp = MIX_CUTOFF_POLICY; | 
| 962 |  | } | 
| 963 |  | } | 
| 964 |  | } | 
| 965 | < | } | 
| 965 | > | } | 
| 966 | > | notifyFortranCutoffPolicy(&cp); | 
| 967 |  |  | 
| 968 | < |  | 
| 968 | > | // Check the Skin Thickness for neighborlists | 
| 969 | > | RealType skin; | 
| 970 |  | if (simParams_->haveSkinThickness()) { | 
| 971 | < | double skinThickness = simParams_->getSkinThickness(); | 
| 972 | < | } | 
| 971 | > | skin = simParams_->getSkinThickness(); | 
| 972 | > | notifyFortranSkinThickness(&skin); | 
| 973 | > | } | 
| 974 | > |  | 
| 975 | > | // Check if the cutoff was set explicitly: | 
| 976 | > | if (simParams_->haveCutoffRadius()) { | 
| 977 | > | rcut_ = simParams_->getCutoffRadius(); | 
| 978 | > | if (simParams_->haveSwitchingRadius()) { | 
| 979 | > | rsw_  = simParams_->getSwitchingRadius(); | 
| 980 | > | } else { | 
| 981 | > | if (fInfo_.SIM_uses_Charges | | 
| 982 | > | fInfo_.SIM_uses_Dipoles | | 
| 983 | > | fInfo_.SIM_uses_RF) { | 
| 984 | > |  | 
| 985 | > | rsw_ = 0.85 * rcut_; | 
| 986 | > | sprintf(painCave.errMsg, | 
| 987 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 988 | > | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" | 
| 989 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 990 | > | painCave.isFatal = 0; | 
| 991 | > | simError(); | 
| 992 | > | } else { | 
| 993 | > | rsw_ = rcut_; | 
| 994 | > | sprintf(painCave.errMsg, | 
| 995 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 996 | > | "\tOOPSE will use the same value as the cutoffRadius.\n" | 
| 997 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 998 | > | painCave.isFatal = 0; | 
| 999 | > | simError(); | 
| 1000 | > | } | 
| 1001 | > | } | 
| 1002 | > |  | 
| 1003 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1004 | > |  | 
| 1005 | > | } else { | 
| 1006 | > |  | 
| 1007 | > | // For electrostatic atoms, we'll assume a large safe value: | 
| 1008 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 1009 | > | sprintf(painCave.errMsg, | 
| 1010 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 1011 | > | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 1012 | > | "\tfor the cutoffRadius.\n"); | 
| 1013 | > | painCave.isFatal = 0; | 
| 1014 | > | simError(); | 
| 1015 | > | rcut_ = 15.0; | 
| 1016 | > |  | 
| 1017 | > | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1018 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1019 | > | toUpper(myMethod); | 
| 1020 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 1021 | > | if (simParams_->haveSwitchingRadius()){ | 
| 1022 | > | sprintf(painCave.errMsg, | 
| 1023 | > | "SimInfo Warning: A value was set for the switchingRadius\n" | 
| 1024 | > | "\teven though the electrostaticSummationMethod was\n" | 
| 1025 | > | "\tset to %s\n", myMethod.c_str()); | 
| 1026 | > | painCave.isFatal = 1; | 
| 1027 | > | simError(); | 
| 1028 | > | } | 
| 1029 | > | } | 
| 1030 | > | } | 
| 1031 | > |  | 
| 1032 | > | if (simParams_->haveSwitchingRadius()){ | 
| 1033 | > | rsw_ = simParams_->getSwitchingRadius(); | 
| 1034 | > | } else { | 
| 1035 | > | sprintf(painCave.errMsg, | 
| 1036 | > | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 1037 | > | "\tOOPSE will use a default value of\n" | 
| 1038 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 1039 | > | painCave.isFatal = 0; | 
| 1040 | > | simError(); | 
| 1041 | > | rsw_ = 0.85 * rcut_; | 
| 1042 | > | } | 
| 1043 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1044 | > | } else { | 
| 1045 | > | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | 
| 1046 | > | // We'll punt and let fortran figure out the cutoffs later. | 
| 1047 | > |  | 
| 1048 | > | notifyFortranYouAreOnYourOwn(); | 
| 1049 |  |  | 
| 1050 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | 
| 1051 | < | // also send cutoff notification to electrostatics | 
| 1033 | < | setElectrostaticCutoffRadius(&rcut_, &rsw_); | 
| 1050 | > | } | 
| 1051 | > | } | 
| 1052 |  | } | 
| 1053 |  |  | 
| 1054 |  | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | 
| 1056 |  | int errorOut; | 
| 1057 |  | int esm =  NONE; | 
| 1058 |  | int sm = UNDAMPED; | 
| 1059 | < | double alphaVal; | 
| 1060 | < | double dielectric; | 
| 1059 | > | RealType alphaVal; | 
| 1060 | > | RealType dielectric; | 
| 1061 |  |  | 
| 1062 |  | errorOut = isError; | 
| 1063 |  | alphaVal = simParams_->getDampingAlpha(); | 
| 1083 |  | } else { | 
| 1084 |  | // throw error | 
| 1085 |  | sprintf( painCave.errMsg, | 
| 1086 | < | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() ); | 
| 1086 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" | 
| 1087 | > | "\t(Input file specified %s .)\n" | 
| 1088 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" | 
| 1089 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" | 
| 1090 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); | 
| 1091 |  | painCave.isFatal = 1; | 
| 1092 |  | simError(); | 
| 1093 |  | } | 
| 1108 |  | if (!simParams_->haveDampingAlpha()) { | 
| 1109 |  | //throw error | 
| 1110 |  | sprintf( painCave.errMsg, | 
| 1111 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); | 
| 1111 | > | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1112 | > | "\tA default value of %f (1/ang) will be used.\n", alphaVal); | 
| 1113 |  | painCave.isFatal = 0; | 
| 1114 |  | simError(); | 
| 1115 |  | } | 
| 1116 |  | } else { | 
| 1117 |  | // throw error | 
| 1118 |  | sprintf( painCave.errMsg, | 
| 1119 | < | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); | 
| 1119 | > | "SimInfo error: Unknown electrostaticScreeningMethod.\n" | 
| 1120 | > | "\t(Input file specified %s .)\n" | 
| 1121 | > | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" | 
| 1122 | > | "or \"damped\".\n", myScreen.c_str() ); | 
| 1123 |  | painCave.isFatal = 1; | 
| 1124 |  | simError(); | 
| 1125 |  | } | 
| 1128 |  |  | 
| 1129 |  | // let's pass some summation method variables to fortran | 
| 1130 |  | setElectrostaticSummationMethod( &esm ); | 
| 1131 | + | setFortranElectrostaticMethod( &esm ); | 
| 1132 |  | setScreeningMethod( &sm ); | 
| 1133 |  | setDampingAlpha( &alphaVal ); | 
| 1134 |  | setReactionFieldDielectric( &dielectric ); | 
| 1135 | < | initFortranFF( &esm, &errorOut ); | 
| 1135 | > | initFortranFF( &errorOut ); | 
| 1136 |  | } | 
| 1137 |  |  | 
| 1138 |  | void SimInfo::setupSwitchingFunction() { | 
| 1217 |  | Molecule* mol; | 
| 1218 |  |  | 
| 1219 |  | Vector3d comVel(0.0); | 
| 1220 | < | double totalMass = 0.0; | 
| 1220 | > | RealType totalMass = 0.0; | 
| 1221 |  |  | 
| 1222 |  |  | 
| 1223 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1224 | < | double mass = mol->getMass(); | 
| 1224 | > | RealType mass = mol->getMass(); | 
| 1225 |  | totalMass += mass; | 
| 1226 |  | comVel += mass * mol->getComVel(); | 
| 1227 |  | } | 
| 1228 |  |  | 
| 1229 |  | #ifdef IS_MPI | 
| 1230 | < | double tmpMass = totalMass; | 
| 1230 | > | RealType tmpMass = totalMass; | 
| 1231 |  | Vector3d tmpComVel(comVel); | 
| 1232 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1233 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1232 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1233 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1234 |  | #endif | 
| 1235 |  |  | 
| 1236 |  | comVel /= totalMass; | 
| 1243 |  | Molecule* mol; | 
| 1244 |  |  | 
| 1245 |  | Vector3d com(0.0); | 
| 1246 | < | double totalMass = 0.0; | 
| 1246 | > | RealType totalMass = 0.0; | 
| 1247 |  |  | 
| 1248 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1249 | < | double mass = mol->getMass(); | 
| 1249 | > | RealType mass = mol->getMass(); | 
| 1250 |  | totalMass += mass; | 
| 1251 |  | com += mass * mol->getCom(); | 
| 1252 |  | } | 
| 1253 |  |  | 
| 1254 |  | #ifdef IS_MPI | 
| 1255 | < | double tmpMass = totalMass; | 
| 1255 | > | RealType tmpMass = totalMass; | 
| 1256 |  | Vector3d tmpCom(com); | 
| 1257 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1258 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1257 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1258 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1259 |  | #endif | 
| 1260 |  |  | 
| 1261 |  | com /= totalMass; | 
| 1279 |  | Molecule* mol; | 
| 1280 |  |  | 
| 1281 |  |  | 
| 1282 | < | double totalMass = 0.0; | 
| 1282 | > | RealType totalMass = 0.0; | 
| 1283 |  |  | 
| 1284 |  |  | 
| 1285 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1286 | < | double mass = mol->getMass(); | 
| 1286 | > | RealType mass = mol->getMass(); | 
| 1287 |  | totalMass += mass; | 
| 1288 |  | com += mass * mol->getCom(); | 
| 1289 |  | comVel += mass * mol->getComVel(); | 
| 1290 |  | } | 
| 1291 |  |  | 
| 1292 |  | #ifdef IS_MPI | 
| 1293 | < | double tmpMass = totalMass; | 
| 1293 | > | RealType tmpMass = totalMass; | 
| 1294 |  | Vector3d tmpCom(com); | 
| 1295 |  | Vector3d tmpComVel(comVel); | 
| 1296 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1297 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1298 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1296 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1297 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1298 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1299 |  | #endif | 
| 1300 |  |  | 
| 1301 |  | com /= totalMass; | 
| 1314 |  | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | 
| 1315 |  |  | 
| 1316 |  |  | 
| 1317 | < | double xx = 0.0; | 
| 1318 | < | double yy = 0.0; | 
| 1319 | < | double zz = 0.0; | 
| 1320 | < | double xy = 0.0; | 
| 1321 | < | double xz = 0.0; | 
| 1322 | < | double yz = 0.0; | 
| 1317 | > | RealType xx = 0.0; | 
| 1318 | > | RealType yy = 0.0; | 
| 1319 | > | RealType zz = 0.0; | 
| 1320 | > | RealType xy = 0.0; | 
| 1321 | > | RealType xz = 0.0; | 
| 1322 | > | RealType yz = 0.0; | 
| 1323 |  | Vector3d com(0.0); | 
| 1324 |  | Vector3d comVel(0.0); | 
| 1325 |  |  | 
| 1331 |  | Vector3d thisq(0.0); | 
| 1332 |  | Vector3d thisv(0.0); | 
| 1333 |  |  | 
| 1334 | < | double thisMass = 0.0; | 
| 1334 | > | RealType thisMass = 0.0; | 
| 1335 |  |  | 
| 1336 |  |  | 
| 1337 |  |  | 
| 1369 |  | #ifdef IS_MPI | 
| 1370 |  | Mat3x3d tmpI(inertiaTensor); | 
| 1371 |  | Vector3d tmpAngMom; | 
| 1372 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1373 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1372 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1373 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1374 |  | #endif | 
| 1375 |  |  | 
| 1376 |  | return; | 
| 1391 |  | Vector3d thisr(0.0); | 
| 1392 |  | Vector3d thisp(0.0); | 
| 1393 |  |  | 
| 1394 | < | double thisMass; | 
| 1394 | > | RealType thisMass; | 
| 1395 |  |  | 
| 1396 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1397 |  | thisMass = mol->getMass(); | 
| 1404 |  |  | 
| 1405 |  | #ifdef IS_MPI | 
| 1406 |  | Vector3d tmpAngMom; | 
| 1407 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1407 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1408 |  | #endif | 
| 1409 |  |  | 
| 1410 |  | return angularMomentum; |