| 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) { | 
| 92 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(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 |  |  | 
| 602 |  |  | 
| 603 |  | setupElectrostaticSummationMethod( isError ); | 
| 604 |  | setupSwitchingFunction(); | 
| 605 | + | setupAccumulateBoxDipole(); | 
| 606 |  |  | 
| 607 |  | if(isError){ | 
| 608 |  | sprintf( painCave.errMsg, | 
| 662 |  | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | 
| 663 |  | int useRF; | 
| 664 |  | int useSF; | 
| 665 | + | int useSP; | 
| 666 | + | int useBoxDipole; | 
| 667 |  | std::string myMethod; | 
| 668 |  |  | 
| 669 |  | // set the useRF logical | 
| 674 |  | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 675 |  | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 676 |  | toUpper(myMethod); | 
| 677 | < | if (myMethod == "REACTION_FIELD") { | 
| 677 | > | if (myMethod == "REACTION_FIELD"){ | 
| 678 |  | useRF=1; | 
| 679 | < | } else { | 
| 680 | < | if (myMethod == "SHIFTED_FORCE") { | 
| 681 | < | useSF = 1; | 
| 682 | < | } | 
| 679 | > | } else if (myMethod == "SHIFTED_FORCE"){ | 
| 680 | > | useSF = 1; | 
| 681 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ | 
| 682 | > | useSP = 1; | 
| 683 |  | } | 
| 684 |  | } | 
| 685 | + |  | 
| 686 | + | if (simParams_->haveAccumulateBoxDipole()) | 
| 687 | + | if (simParams_->getAccumulateBoxDipole()) | 
| 688 | + | useBoxDipole = 1; | 
| 689 |  |  | 
| 690 |  | //loop over all of the atom types | 
| 691 |  | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 756 |  | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 757 |  |  | 
| 758 |  | temp = useSF; | 
| 759 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 759 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 760 | > |  | 
| 761 | > | temp = useSP; | 
| 762 | > | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 763 | > |  | 
| 764 | > | temp = useBoxDipole; | 
| 765 | > | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | 
| 766 |  |  | 
| 767 |  | #endif | 
| 768 |  |  | 
| 781 |  | fInfo_.SIM_uses_FLARB = useFLARB; | 
| 782 |  | fInfo_.SIM_uses_RF = useRF; | 
| 783 |  | fInfo_.SIM_uses_SF = useSF; | 
| 784 | + | fInfo_.SIM_uses_SP = useSP; | 
| 785 | + | fInfo_.SIM_uses_BoxDipole = useBoxDipole; | 
| 786 |  |  | 
| 787 |  | if( myMethod == "REACTION_FIELD") { | 
| 788 |  |  | 
| 814 |  | } | 
| 815 |  |  | 
| 816 |  | //calculate mass ratio of cutoff group | 
| 817 | < | std::vector<double> mfact; | 
| 817 | > | std::vector<RealType> mfact; | 
| 818 |  | SimInfo::MoleculeIterator mi; | 
| 819 |  | Molecule* mol; | 
| 820 |  | Molecule::CutoffGroupIterator ci; | 
| 821 |  | CutoffGroup* cg; | 
| 822 |  | Molecule::AtomIterator ai; | 
| 823 |  | Atom* atom; | 
| 824 | < | double totalMass; | 
| 824 | > | RealType totalMass; | 
| 825 |  |  | 
| 826 |  | //to avoid memory reallocation, reserve enough space for mfact | 
| 827 |  | mfact.reserve(getNCutoffGroups()); | 
| 944 |  |  | 
| 945 |  | #endif | 
| 946 |  |  | 
| 947 | < | double SimInfo::calcMaxCutoffRadius() { | 
| 947 | > | void SimInfo::setupCutoff() { | 
| 948 | > |  | 
| 949 | > | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | 
| 950 |  |  | 
| 951 | + | // Check the cutoff policy | 
| 952 | + | int cp =  TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | 
| 953 |  |  | 
| 954 | < | std::set<AtomType*> atomTypes; | 
| 955 | < | std::set<AtomType*>::iterator i; | 
| 956 | < | std::vector<double> cutoffRadius; | 
| 957 | < |  | 
| 958 | < | //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)); | 
| 954 | > | std::string myPolicy; | 
| 955 | > | if (forceFieldOptions_.haveCutoffPolicy()){ | 
| 956 | > | myPolicy = forceFieldOptions_.getCutoffPolicy(); | 
| 957 | > | }else if (simParams_->haveCutoffPolicy()) { | 
| 958 | > | myPolicy = simParams_->getCutoffPolicy(); | 
| 959 |  | } | 
| 960 |  |  | 
| 961 | < | 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(); | 
| 961 | > | if (!myPolicy.empty()){ | 
| 962 |  | toUpper(myPolicy); | 
| 963 |  | if (myPolicy == "MIX") { | 
| 964 |  | cp = MIX_CUTOFF_POLICY; | 
| 977 |  | } | 
| 978 |  | } | 
| 979 |  | } | 
| 980 | < | } | 
| 980 | > | } | 
| 981 | > | notifyFortranCutoffPolicy(&cp); | 
| 982 |  |  | 
| 983 | < |  | 
| 983 | > | // Check the Skin Thickness for neighborlists | 
| 984 | > | RealType skin; | 
| 985 |  | if (simParams_->haveSkinThickness()) { | 
| 986 | < | double skinThickness = simParams_->getSkinThickness(); | 
| 987 | < | } | 
| 986 | > | skin = simParams_->getSkinThickness(); | 
| 987 | > | notifyFortranSkinThickness(&skin); | 
| 988 | > | } | 
| 989 | > |  | 
| 990 | > | // Check if the cutoff was set explicitly: | 
| 991 | > | if (simParams_->haveCutoffRadius()) { | 
| 992 | > | rcut_ = simParams_->getCutoffRadius(); | 
| 993 | > | if (simParams_->haveSwitchingRadius()) { | 
| 994 | > | rsw_  = simParams_->getSwitchingRadius(); | 
| 995 | > | } else { | 
| 996 | > | if (fInfo_.SIM_uses_Charges | | 
| 997 | > | fInfo_.SIM_uses_Dipoles | | 
| 998 | > | fInfo_.SIM_uses_RF) { | 
| 999 | > |  | 
| 1000 | > | rsw_ = 0.85 * rcut_; | 
| 1001 | > | sprintf(painCave.errMsg, | 
| 1002 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 1003 | > | "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" | 
| 1004 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1005 | > | painCave.isFatal = 0; | 
| 1006 | > | simError(); | 
| 1007 | > | } else { | 
| 1008 | > | rsw_ = rcut_; | 
| 1009 | > | sprintf(painCave.errMsg, | 
| 1010 | > | "SimCreator Warning: No value was set for the switchingRadius.\n" | 
| 1011 | > | "\tOOPSE will use the same value as the cutoffRadius.\n" | 
| 1012 | > | "\tswitchingRadius = %f. for this simulation\n", rsw_); | 
| 1013 | > | painCave.isFatal = 0; | 
| 1014 | > | simError(); | 
| 1015 | > | } | 
| 1016 | > | } | 
| 1017 | > |  | 
| 1018 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1019 | > |  | 
| 1020 | > | } else { | 
| 1021 | > |  | 
| 1022 | > | // For electrostatic atoms, we'll assume a large safe value: | 
| 1023 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | 
| 1024 | > | sprintf(painCave.errMsg, | 
| 1025 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" | 
| 1026 | > | "\tOOPSE will use a default value of 15.0 angstroms" | 
| 1027 | > | "\tfor the cutoffRadius.\n"); | 
| 1028 | > | painCave.isFatal = 0; | 
| 1029 | > | simError(); | 
| 1030 | > | rcut_ = 15.0; | 
| 1031 | > |  | 
| 1032 | > | if (simParams_->haveElectrostaticSummationMethod()) { | 
| 1033 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | 
| 1034 | > | toUpper(myMethod); | 
| 1035 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { | 
| 1036 | > | if (simParams_->haveSwitchingRadius()){ | 
| 1037 | > | sprintf(painCave.errMsg, | 
| 1038 | > | "SimInfo Warning: A value was set for the switchingRadius\n" | 
| 1039 | > | "\teven though the electrostaticSummationMethod was\n" | 
| 1040 | > | "\tset to %s\n", myMethod.c_str()); | 
| 1041 | > | painCave.isFatal = 1; | 
| 1042 | > | simError(); | 
| 1043 | > | } | 
| 1044 | > | } | 
| 1045 | > | } | 
| 1046 | > |  | 
| 1047 | > | if (simParams_->haveSwitchingRadius()){ | 
| 1048 | > | rsw_ = simParams_->getSwitchingRadius(); | 
| 1049 | > | } else { | 
| 1050 | > | sprintf(painCave.errMsg, | 
| 1051 | > | "SimCreator Warning: No value was set for switchingRadius.\n" | 
| 1052 | > | "\tOOPSE will use a default value of\n" | 
| 1053 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); | 
| 1054 | > | painCave.isFatal = 0; | 
| 1055 | > | simError(); | 
| 1056 | > | rsw_ = 0.85 * rcut_; | 
| 1057 | > | } | 
| 1058 | > | notifyFortranCutoffs(&rcut_, &rsw_); | 
| 1059 | > | } else { | 
| 1060 | > | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | 
| 1061 | > | // We'll punt and let fortran figure out the cutoffs later. | 
| 1062 | > |  | 
| 1063 | > | notifyFortranYouAreOnYourOwn(); | 
| 1064 |  |  | 
| 1065 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); | 
| 1066 | < | // also send cutoff notification to electrostatics | 
| 1033 | < | setElectrostaticCutoffRadius(&rcut_, &rsw_); | 
| 1065 | > | } | 
| 1066 | > | } | 
| 1067 |  | } | 
| 1068 |  |  | 
| 1069 |  | void SimInfo::setupElectrostaticSummationMethod( int isError ) { | 
| 1071 |  | int errorOut; | 
| 1072 |  | int esm =  NONE; | 
| 1073 |  | int sm = UNDAMPED; | 
| 1074 | < | double alphaVal; | 
| 1075 | < | double dielectric; | 
| 1074 | > | RealType alphaVal; | 
| 1075 | > | RealType dielectric; | 
| 1076 |  |  | 
| 1077 |  | errorOut = isError; | 
| 1078 |  | alphaVal = simParams_->getDampingAlpha(); | 
| 1098 |  | } else { | 
| 1099 |  | // throw error | 
| 1100 |  | sprintf( painCave.errMsg, | 
| 1101 | < | "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() ); | 
| 1101 | > | "SimInfo error: Unknown electrostaticSummationMethod.\n" | 
| 1102 | > | "\t(Input file specified %s .)\n" | 
| 1103 | > | "\telectrostaticSummationMethod must be one of: \"none\",\n" | 
| 1104 | > | "\t\"shifted_potential\", \"shifted_force\", or \n" | 
| 1105 | > | "\t\"reaction_field\".\n", myMethod.c_str() ); | 
| 1106 |  | painCave.isFatal = 1; | 
| 1107 |  | simError(); | 
| 1108 |  | } | 
| 1123 |  | if (!simParams_->haveDampingAlpha()) { | 
| 1124 |  | //throw error | 
| 1125 |  | sprintf( painCave.errMsg, | 
| 1126 | < | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); | 
| 1126 | > | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | 
| 1127 | > | "\tA default value of %f (1/ang) will be used.\n", alphaVal); | 
| 1128 |  | painCave.isFatal = 0; | 
| 1129 |  | simError(); | 
| 1130 |  | } | 
| 1131 |  | } else { | 
| 1132 |  | // throw error | 
| 1133 |  | sprintf( painCave.errMsg, | 
| 1134 | < | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); | 
| 1134 | > | "SimInfo error: Unknown electrostaticScreeningMethod.\n" | 
| 1135 | > | "\t(Input file specified %s .)\n" | 
| 1136 | > | "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" | 
| 1137 | > | "or \"damped\".\n", myScreen.c_str() ); | 
| 1138 |  | painCave.isFatal = 1; | 
| 1139 |  | simError(); | 
| 1140 |  | } | 
| 1143 |  |  | 
| 1144 |  | // let's pass some summation method variables to fortran | 
| 1145 |  | setElectrostaticSummationMethod( &esm ); | 
| 1146 | + | setFortranElectrostaticMethod( &esm ); | 
| 1147 |  | setScreeningMethod( &sm ); | 
| 1148 |  | setDampingAlpha( &alphaVal ); | 
| 1149 |  | setReactionFieldDielectric( &dielectric ); | 
| 1150 | < | initFortranFF( &esm, &errorOut ); | 
| 1150 | > | initFortranFF( &errorOut ); | 
| 1151 |  | } | 
| 1152 |  |  | 
| 1153 |  | void SimInfo::setupSwitchingFunction() { | 
| 1176 |  |  | 
| 1177 |  | } | 
| 1178 |  |  | 
| 1179 | + | void SimInfo::setupAccumulateBoxDipole() { | 
| 1180 | + |  | 
| 1181 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 1182 | + | if ( simParams_->haveAccumulateBoxDipole() ) | 
| 1183 | + | if ( simParams_->getAccumulateBoxDipole() ) { | 
| 1184 | + | setAccumulateBoxDipole(); | 
| 1185 | + | calcBoxDipole_ = true; | 
| 1186 | + | } | 
| 1187 | + |  | 
| 1188 | + | } | 
| 1189 | + |  | 
| 1190 |  | void SimInfo::addProperty(GenericData* genData) { | 
| 1191 |  | properties_.addProperty(genData); | 
| 1192 |  | } | 
| 1243 |  | Molecule* mol; | 
| 1244 |  |  | 
| 1245 |  | Vector3d comVel(0.0); | 
| 1246 | < | double totalMass = 0.0; | 
| 1246 | > | RealType totalMass = 0.0; | 
| 1247 |  |  | 
| 1248 |  |  | 
| 1249 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1250 | < | double mass = mol->getMass(); | 
| 1250 | > | RealType mass = mol->getMass(); | 
| 1251 |  | totalMass += mass; | 
| 1252 |  | comVel += mass * mol->getComVel(); | 
| 1253 |  | } | 
| 1254 |  |  | 
| 1255 |  | #ifdef IS_MPI | 
| 1256 | < | double tmpMass = totalMass; | 
| 1256 | > | RealType tmpMass = totalMass; | 
| 1257 |  | Vector3d tmpComVel(comVel); | 
| 1258 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1259 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1258 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1259 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1260 |  | #endif | 
| 1261 |  |  | 
| 1262 |  | comVel /= totalMass; | 
| 1269 |  | Molecule* mol; | 
| 1270 |  |  | 
| 1271 |  | Vector3d com(0.0); | 
| 1272 | < | double totalMass = 0.0; | 
| 1272 | > | RealType totalMass = 0.0; | 
| 1273 |  |  | 
| 1274 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1275 | < | double mass = mol->getMass(); | 
| 1275 | > | RealType mass = mol->getMass(); | 
| 1276 |  | totalMass += mass; | 
| 1277 |  | com += mass * mol->getCom(); | 
| 1278 |  | } | 
| 1279 |  |  | 
| 1280 |  | #ifdef IS_MPI | 
| 1281 | < | double tmpMass = totalMass; | 
| 1281 | > | RealType tmpMass = totalMass; | 
| 1282 |  | Vector3d tmpCom(com); | 
| 1283 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1284 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1283 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1284 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1285 |  | #endif | 
| 1286 |  |  | 
| 1287 |  | com /= totalMass; | 
| 1305 |  | Molecule* mol; | 
| 1306 |  |  | 
| 1307 |  |  | 
| 1308 | < | double totalMass = 0.0; | 
| 1308 | > | RealType totalMass = 0.0; | 
| 1309 |  |  | 
| 1310 |  |  | 
| 1311 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1312 | < | double mass = mol->getMass(); | 
| 1312 | > | RealType mass = mol->getMass(); | 
| 1313 |  | totalMass += mass; | 
| 1314 |  | com += mass * mol->getCom(); | 
| 1315 |  | comVel += mass * mol->getComVel(); | 
| 1316 |  | } | 
| 1317 |  |  | 
| 1318 |  | #ifdef IS_MPI | 
| 1319 | < | double tmpMass = totalMass; | 
| 1319 | > | RealType tmpMass = totalMass; | 
| 1320 |  | Vector3d tmpCom(com); | 
| 1321 |  | Vector3d tmpComVel(comVel); | 
| 1322 | < | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1323 | < | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1324 | < | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1322 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1323 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1324 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1325 |  | #endif | 
| 1326 |  |  | 
| 1327 |  | com /= totalMass; | 
| 1340 |  | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ | 
| 1341 |  |  | 
| 1342 |  |  | 
| 1343 | < | double xx = 0.0; | 
| 1344 | < | double yy = 0.0; | 
| 1345 | < | double zz = 0.0; | 
| 1346 | < | double xy = 0.0; | 
| 1347 | < | double xz = 0.0; | 
| 1348 | < | double yz = 0.0; | 
| 1343 | > | RealType xx = 0.0; | 
| 1344 | > | RealType yy = 0.0; | 
| 1345 | > | RealType zz = 0.0; | 
| 1346 | > | RealType xy = 0.0; | 
| 1347 | > | RealType xz = 0.0; | 
| 1348 | > | RealType yz = 0.0; | 
| 1349 |  | Vector3d com(0.0); | 
| 1350 |  | Vector3d comVel(0.0); | 
| 1351 |  |  | 
| 1357 |  | Vector3d thisq(0.0); | 
| 1358 |  | Vector3d thisv(0.0); | 
| 1359 |  |  | 
| 1360 | < | double thisMass = 0.0; | 
| 1360 | > | RealType thisMass = 0.0; | 
| 1361 |  |  | 
| 1362 |  |  | 
| 1363 |  |  | 
| 1395 |  | #ifdef IS_MPI | 
| 1396 |  | Mat3x3d tmpI(inertiaTensor); | 
| 1397 |  | Vector3d tmpAngMom; | 
| 1398 | < | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1399 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1398 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1399 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1400 |  | #endif | 
| 1401 |  |  | 
| 1402 |  | return; | 
| 1417 |  | Vector3d thisr(0.0); | 
| 1418 |  | Vector3d thisp(0.0); | 
| 1419 |  |  | 
| 1420 | < | double thisMass; | 
| 1420 | > | RealType thisMass; | 
| 1421 |  |  | 
| 1422 |  | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { | 
| 1423 |  | thisMass = mol->getMass(); | 
| 1430 |  |  | 
| 1431 |  | #ifdef IS_MPI | 
| 1432 |  | Vector3d tmpAngMom; | 
| 1433 | < | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); | 
| 1433 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); | 
| 1434 |  | #endif | 
| 1435 |  |  | 
| 1436 |  | return angularMomentum; |