# | Line 52 | Line 52 | |
---|---|---|
52 | #include "brains/SimInfo.hpp" | |
53 | #include "math/Vector3.hpp" | |
54 | #include "primitives/Molecule.hpp" | |
55 | + | #include "UseTheForce/fCutoffPolicy.h" |
56 | + | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
57 | + | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" |
58 | + | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" |
59 | #include "UseTheForce/doForces_interface.h" | |
60 | + | #include "UseTheForce/DarkSide/electrostatic_interface.h" |
61 | #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" | |
# | Line 80 | Line 86 | namespace oopse { | |
86 | MoleculeStamp* molStamp; | |
87 | int nMolWithSameStamp; | |
88 | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups | |
89 | < | int nGroups = 0; //total cutoff groups defined in meta-data file |
89 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
90 | CutoffGroupStamp* cgStamp; | |
91 | RigidBodyStamp* rbStamp; | |
92 | int nRigidAtoms = 0; | |
# | Line 105 | Line 111 | namespace oopse { | |
111 | } | |
112 | ||
113 | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; | |
114 | + | |
115 | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; | |
116 | ||
117 | //calculate atoms in rigid bodies | |
# | Line 121 | Line 128 | namespace oopse { | |
128 | ||
129 | } | |
130 | ||
131 | < | //every free atom (atom does not belong to cutoff groups) is a cutoff group |
132 | < | //therefore the total number of cutoff groups in the system is equal to |
133 | < | //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data |
134 | < | //file plus the number of cutoff groups defined in meta-data file |
131 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff |
132 | > | //group therefore the total number of cutoff groups in the system is |
133 | > | //equal to the total number of atoms minus number of atoms belong to |
134 | > | //cutoff group defined in meta-data file plus the number of cutoff |
135 | > | //groups defined in meta-data file |
136 | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | |
129 | – | |
130 | – | //every free atom (atom does not belong to rigid bodies) is an integrable object |
131 | – | //therefore the total number of integrable objects in the system is equal to |
132 | – | //the total number of atoms minus number of atoms belong to rigid body defined in meta-data |
133 | – | //file plus the number of rigid bodies defined in meta-data file |
134 | – | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; |
137 | ||
138 | + | //every free atom (atom does not belong to rigid bodies) is an |
139 | + | //integrable object therefore the total number of integrable objects |
140 | + | //in the system is equal to the total number of atoms minus number of |
141 | + | //atoms belong to rigid body defined in meta-data file plus the number |
142 | + | //of rigid bodies defined in meta-data file |
143 | + | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
144 | + | + nGlobalRigidBodies_; |
145 | + | |
146 | nGlobalMols_ = molStampIds_.size(); | |
147 | ||
148 | #ifdef IS_MPI | |
# | Line 462 | Line 472 | namespace oopse { | |
472 | //setup fortran force field | |
473 | /** @deprecate */ | |
474 | int isError = 0; | |
475 | < | initFortranFF( &fInfo_.SIM_uses_RF , &isError ); |
475 | > | |
476 | > | setupElectrostaticSummationMethod( isError ); |
477 | > | setupSwitchingFunction(); |
478 | > | |
479 | if(isError){ | |
480 | sprintf( painCave.errMsg, | |
481 | "ForceField error: There was an error initializing the forceField in fortran.\n" ); | |
# | Line 517 | Line 530 | namespace oopse { | |
530 | int useDirectionalAtom = 0; | |
531 | int useElectrostatics = 0; | |
532 | //usePBC and useRF are from simParams | |
533 | < | int usePBC = simParams_->getPBC(); |
534 | < | int useRF = simParams_->getUseRF(); |
533 | > | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
534 | > | int useRF; |
535 | > | int useSF; |
536 | > | std::string myMethod; |
537 | > | |
538 | > | // set the useRF logical |
539 | > | useRF = 0; |
540 | > | useSF = 0; |
541 | > | |
542 | > | |
543 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
544 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
545 | > | toUpper(myMethod); |
546 | > | if (myMethod == "REACTION_FIELD") { |
547 | > | useRF=1; |
548 | > | } else { |
549 | > | if (myMethod == "SHIFTED_FORCE") { |
550 | > | useSF = 1; |
551 | > | } |
552 | > | } |
553 | > | } |
554 | ||
555 | //loop over all of the atom types | |
556 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
# | Line 583 | Line 615 | namespace oopse { | |
615 | ||
616 | temp = useRF; | |
617 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
618 | < | |
618 | > | |
619 | > | temp = useSF; |
620 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
621 | > | |
622 | #endif | |
623 | ||
624 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 599 | Line 634 | namespace oopse { | |
634 | fInfo_.SIM_uses_Shapes = useShape; | |
635 | fInfo_.SIM_uses_FLARB = useFLARB; | |
636 | fInfo_.SIM_uses_RF = useRF; | |
637 | + | fInfo_.SIM_uses_SF = useSF; |
638 | ||
639 | < | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { |
640 | < | |
639 | > | if( myMethod == "REACTION_FIELD") { |
640 | > | |
641 | if (simParams_->haveDielectric()) { | |
642 | fInfo_.dielect = simParams_->getDielectric(); | |
643 | } else { | |
# | Line 611 | Line 647 | namespace oopse { | |
647 | "\tsetting a dielectric constant!\n"); | |
648 | painCave.isFatal = 1; | |
649 | simError(); | |
650 | < | } |
615 | < | |
616 | < | } else { |
617 | < | fInfo_.dielect = 0.0; |
650 | > | } |
651 | } | |
652 | ||
653 | } | |
# | Line 650 | Line 683 | namespace oopse { | |
683 | ||
684 | totalMass = cg->getMass(); | |
685 | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | |
686 | < | mfact.push_back(atom->getMass()/totalMass); |
686 | > | // Check for massless groups - set mfact to 1 if true |
687 | > | if (totalMass != 0) |
688 | > | mfact.push_back(atom->getMass()/totalMass); |
689 | > | else |
690 | > | mfact.push_back( 1.0 ); |
691 | } | |
692 | ||
693 | } | |
# | Line 786 | Line 823 | namespace oopse { | |
823 | ||
824 | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { | |
825 | ||
826 | < | if (!simParams_->haveRcut()){ |
826 | > | if (!simParams_->haveCutoffRadius()){ |
827 | sprintf(painCave.errMsg, | |
828 | "SimCreator Warning: No value was set for the cutoffRadius.\n" | |
829 | "\tOOPSE will use a default value of 15.0 angstroms" | |
# | Line 795 | Line 832 | namespace oopse { | |
832 | simError(); | |
833 | rcut = 15.0; | |
834 | } else{ | |
835 | < | rcut = simParams_->getRcut(); |
835 | > | rcut = simParams_->getCutoffRadius(); |
836 | } | |
837 | ||
838 | < | if (!simParams_->haveRsw()){ |
838 | > | if (!simParams_->haveSwitchingRadius()){ |
839 | sprintf(painCave.errMsg, | |
840 | "SimCreator Warning: No value was set for switchingRadius.\n" | |
841 | "\tOOPSE will use a default value of\n" | |
842 | < | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
842 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
843 | painCave.isFatal = 0; | |
844 | simError(); | |
845 | < | rsw = 0.95 * rcut; |
845 | > | rsw = 0.85 * rcut; |
846 | } else{ | |
847 | < | rsw = simParams_->getRsw(); |
847 | > | rsw = simParams_->getSwitchingRadius(); |
848 | } | |
849 | ||
850 | } else { | |
851 | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in | |
852 | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used | |
853 | ||
854 | < | if (simParams_->haveRcut()) { |
855 | < | rcut = simParams_->getRcut(); |
854 | > | if (simParams_->haveCutoffRadius()) { |
855 | > | rcut = simParams_->getCutoffRadius(); |
856 | } else { | |
857 | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system | |
858 | rcut = calcMaxCutoffRadius(); | |
859 | } | |
860 | ||
861 | < | if (simParams_->haveRsw()) { |
862 | < | rsw = simParams_->getRsw(); |
861 | > | if (simParams_->haveSwitchingRadius()) { |
862 | > | rsw = simParams_->getSwitchingRadius(); |
863 | } else { | |
864 | rsw = rcut; | |
865 | } | |
# | Line 830 | Line 867 | namespace oopse { | |
867 | } | |
868 | } | |
869 | ||
870 | < | void SimInfo::setupCutoff() { |
870 | > | void SimInfo::setupCutoff() { |
871 | getCutoff(rcut_, rsw_); | |
872 | double rnblist = rcut_ + 1; // skin of neighbor list | |
873 | ||
874 | //Pass these cutoff radius etc. to fortran. This function should be called once and only once | |
875 | < | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); |
875 | > | |
876 | > | int cp = TRADITIONAL_CUTOFF_POLICY; |
877 | > | if (simParams_->haveCutoffPolicy()) { |
878 | > | std::string myPolicy = simParams_->getCutoffPolicy(); |
879 | > | toUpper(myPolicy); |
880 | > | if (myPolicy == "MIX") { |
881 | > | cp = MIX_CUTOFF_POLICY; |
882 | > | } else { |
883 | > | if (myPolicy == "MAX") { |
884 | > | cp = MAX_CUTOFF_POLICY; |
885 | > | } else { |
886 | > | if (myPolicy == "TRADITIONAL") { |
887 | > | cp = TRADITIONAL_CUTOFF_POLICY; |
888 | > | } else { |
889 | > | // throw error |
890 | > | sprintf( painCave.errMsg, |
891 | > | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
892 | > | painCave.isFatal = 1; |
893 | > | simError(); |
894 | > | } |
895 | > | } |
896 | > | } |
897 | > | } |
898 | > | |
899 | > | |
900 | > | if (simParams_->haveSkinThickness()) { |
901 | > | double skinThickness = simParams_->getSkinThickness(); |
902 | > | } |
903 | > | |
904 | > | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
905 | > | // also send cutoff notification to electrostatics |
906 | > | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
907 | } | |
908 | ||
909 | + | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
910 | + | |
911 | + | int errorOut; |
912 | + | int esm = NONE; |
913 | + | int sm = UNDAMPED; |
914 | + | double alphaVal; |
915 | + | double dielectric; |
916 | + | |
917 | + | errorOut = isError; |
918 | + | alphaVal = simParams_->getDampingAlpha(); |
919 | + | dielectric = simParams_->getDielectric(); |
920 | + | |
921 | + | if (simParams_->haveElectrostaticSummationMethod()) { |
922 | + | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
923 | + | toUpper(myMethod); |
924 | + | if (myMethod == "NONE") { |
925 | + | esm = NONE; |
926 | + | } else { |
927 | + | if (myMethod == "SWITCHING_FUNCTION") { |
928 | + | esm = SWITCHING_FUNCTION; |
929 | + | } else { |
930 | + | if (myMethod == "SHIFTED_POTENTIAL") { |
931 | + | esm = SHIFTED_POTENTIAL; |
932 | + | } else { |
933 | + | if (myMethod == "SHIFTED_FORCE") { |
934 | + | esm = SHIFTED_FORCE; |
935 | + | } else { |
936 | + | if (myMethod == "REACTION_FIELD") { |
937 | + | esm = REACTION_FIELD; |
938 | + | } else { |
939 | + | // throw error |
940 | + | sprintf( painCave.errMsg, |
941 | + | "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() ); |
942 | + | painCave.isFatal = 1; |
943 | + | simError(); |
944 | + | } |
945 | + | } |
946 | + | } |
947 | + | } |
948 | + | } |
949 | + | } |
950 | + | |
951 | + | if (simParams_->haveElectrostaticScreeningMethod()) { |
952 | + | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
953 | + | toUpper(myScreen); |
954 | + | if (myScreen == "UNDAMPED") { |
955 | + | sm = UNDAMPED; |
956 | + | } else { |
957 | + | if (myScreen == "DAMPED") { |
958 | + | sm = DAMPED; |
959 | + | if (!simParams_->haveDampingAlpha()) { |
960 | + | //throw error |
961 | + | sprintf( painCave.errMsg, |
962 | + | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); |
963 | + | painCave.isFatal = 0; |
964 | + | simError(); |
965 | + | } |
966 | + | } else { |
967 | + | // throw error |
968 | + | sprintf( painCave.errMsg, |
969 | + | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
970 | + | painCave.isFatal = 1; |
971 | + | simError(); |
972 | + | } |
973 | + | } |
974 | + | } |
975 | + | |
976 | + | // let's pass some summation method variables to fortran |
977 | + | setElectrostaticSummationMethod( &esm ); |
978 | + | setScreeningMethod( &sm ); |
979 | + | setDampingAlpha( &alphaVal ); |
980 | + | setReactionFieldDielectric( &dielectric ); |
981 | + | initFortranFF( &esm, &errorOut ); |
982 | + | } |
983 | + | |
984 | + | void SimInfo::setupSwitchingFunction() { |
985 | + | int ft = CUBIC; |
986 | + | |
987 | + | if (simParams_->haveSwitchingFunctionType()) { |
988 | + | std::string funcType = simParams_->getSwitchingFunctionType(); |
989 | + | toUpper(funcType); |
990 | + | if (funcType == "CUBIC") { |
991 | + | ft = CUBIC; |
992 | + | } else { |
993 | + | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
994 | + | ft = FIFTH_ORDER_POLY; |
995 | + | } else { |
996 | + | // throw error |
997 | + | sprintf( painCave.errMsg, |
998 | + | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
999 | + | painCave.isFatal = 1; |
1000 | + | simError(); |
1001 | + | } |
1002 | + | } |
1003 | + | } |
1004 | + | |
1005 | + | // send switching function notification to switcheroo |
1006 | + | setFunctionType(&ft); |
1007 | + | |
1008 | + | } |
1009 | + | |
1010 | void SimInfo::addProperty(GenericData* genData) { | |
1011 | properties_.addProperty(genData); | |
1012 | } | |
# | Line 945 | Line 1114 | namespace oopse { | |
1114 | ||
1115 | return o; | |
1116 | } | |
1117 | + | |
1118 | + | |
1119 | + | /* |
1120 | + | Returns center of mass and center of mass velocity in one function call. |
1121 | + | */ |
1122 | + | |
1123 | + | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
1124 | + | SimInfo::MoleculeIterator i; |
1125 | + | Molecule* mol; |
1126 | + | |
1127 | + | |
1128 | + | double totalMass = 0.0; |
1129 | + | |
1130 | ||
1131 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1132 | + | double mass = mol->getMass(); |
1133 | + | totalMass += mass; |
1134 | + | com += mass * mol->getCom(); |
1135 | + | comVel += mass * mol->getComVel(); |
1136 | + | } |
1137 | + | |
1138 | + | #ifdef IS_MPI |
1139 | + | double tmpMass = totalMass; |
1140 | + | Vector3d tmpCom(com); |
1141 | + | Vector3d tmpComVel(comVel); |
1142 | + | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1143 | + | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1144 | + | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1145 | + | #endif |
1146 | + | |
1147 | + | com /= totalMass; |
1148 | + | comVel /= totalMass; |
1149 | + | } |
1150 | + | |
1151 | + | /* |
1152 | + | Return intertia tensor for entire system and angular momentum Vector. |
1153 | + | |
1154 | + | |
1155 | + | [ Ixx -Ixy -Ixz ] |
1156 | + | J =| -Iyx Iyy -Iyz | |
1157 | + | [ -Izx -Iyz Izz ] |
1158 | + | */ |
1159 | + | |
1160 | + | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1161 | + | |
1162 | + | |
1163 | + | double xx = 0.0; |
1164 | + | double yy = 0.0; |
1165 | + | double zz = 0.0; |
1166 | + | double xy = 0.0; |
1167 | + | double xz = 0.0; |
1168 | + | double yz = 0.0; |
1169 | + | Vector3d com(0.0); |
1170 | + | Vector3d comVel(0.0); |
1171 | + | |
1172 | + | getComAll(com, comVel); |
1173 | + | |
1174 | + | SimInfo::MoleculeIterator i; |
1175 | + | Molecule* mol; |
1176 | + | |
1177 | + | Vector3d thisq(0.0); |
1178 | + | Vector3d thisv(0.0); |
1179 | + | |
1180 | + | double thisMass = 0.0; |
1181 | + | |
1182 | + | |
1183 | + | |
1184 | + | |
1185 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1186 | + | |
1187 | + | thisq = mol->getCom()-com; |
1188 | + | thisv = mol->getComVel()-comVel; |
1189 | + | thisMass = mol->getMass(); |
1190 | + | // Compute moment of intertia coefficients. |
1191 | + | xx += thisq[0]*thisq[0]*thisMass; |
1192 | + | yy += thisq[1]*thisq[1]*thisMass; |
1193 | + | zz += thisq[2]*thisq[2]*thisMass; |
1194 | + | |
1195 | + | // compute products of intertia |
1196 | + | xy += thisq[0]*thisq[1]*thisMass; |
1197 | + | xz += thisq[0]*thisq[2]*thisMass; |
1198 | + | yz += thisq[1]*thisq[2]*thisMass; |
1199 | + | |
1200 | + | angularMomentum += cross( thisq, thisv ) * thisMass; |
1201 | + | |
1202 | + | } |
1203 | + | |
1204 | + | |
1205 | + | inertiaTensor(0,0) = yy + zz; |
1206 | + | inertiaTensor(0,1) = -xy; |
1207 | + | inertiaTensor(0,2) = -xz; |
1208 | + | inertiaTensor(1,0) = -xy; |
1209 | + | inertiaTensor(1,1) = xx + zz; |
1210 | + | inertiaTensor(1,2) = -yz; |
1211 | + | inertiaTensor(2,0) = -xz; |
1212 | + | inertiaTensor(2,1) = -yz; |
1213 | + | inertiaTensor(2,2) = xx + yy; |
1214 | + | |
1215 | + | #ifdef IS_MPI |
1216 | + | Mat3x3d tmpI(inertiaTensor); |
1217 | + | Vector3d tmpAngMom; |
1218 | + | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1219 | + | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1220 | + | #endif |
1221 | + | |
1222 | + | return; |
1223 | + | } |
1224 | + | |
1225 | + | //Returns the angular momentum of the system |
1226 | + | Vector3d SimInfo::getAngularMomentum(){ |
1227 | + | |
1228 | + | Vector3d com(0.0); |
1229 | + | Vector3d comVel(0.0); |
1230 | + | Vector3d angularMomentum(0.0); |
1231 | + | |
1232 | + | getComAll(com,comVel); |
1233 | + | |
1234 | + | SimInfo::MoleculeIterator i; |
1235 | + | Molecule* mol; |
1236 | + | |
1237 | + | Vector3d thisr(0.0); |
1238 | + | Vector3d thisp(0.0); |
1239 | + | |
1240 | + | double thisMass; |
1241 | + | |
1242 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1243 | + | thisMass = mol->getMass(); |
1244 | + | thisr = mol->getCom()-com; |
1245 | + | thisp = (mol->getComVel()-comVel)*thisMass; |
1246 | + | |
1247 | + | angularMomentum += cross( thisr, thisp ); |
1248 | + | |
1249 | + | } |
1250 | + | |
1251 | + | #ifdef IS_MPI |
1252 | + | Vector3d tmpAngMom; |
1253 | + | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1254 | + | #endif |
1255 | + | |
1256 | + | return angularMomentum; |
1257 | + | } |
1258 | + | |
1259 | + | |
1260 | }//end namespace oopse | |
1261 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |