# | Line 53 | Line 53 | |
---|---|---|
53 | #include "brains/SimInfo.hpp" | |
54 | #include "math/Vector3.hpp" | |
55 | #include "primitives/Molecule.hpp" | |
56 | + | #include "primitives/StuntDouble.hpp" |
57 | #include "UseTheForce/fCutoffPolicy.h" | |
58 | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" | |
59 | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" | |
60 | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" | |
61 | #include "UseTheForce/doForces_interface.h" | |
62 | + | #include "UseTheForce/DarkSide/neighborLists_interface.h" |
63 | #include "UseTheForce/DarkSide/electrostatic_interface.h" | |
64 | #include "UseTheForce/DarkSide/switcheroo_interface.h" | |
65 | #include "utils/MemoryUtils.hpp" | |
# | Line 66 | Line 68 | |
68 | #include "io/ForceFieldOptions.hpp" | |
69 | #include "UseTheForce/ForceField.hpp" | |
70 | ||
71 | + | |
72 | #ifdef IS_MPI | |
73 | #include "UseTheForce/mpiComponentPlan.h" | |
74 | #include "UseTheForce/DarkSide/simParallel_interface.h" | |
# | Line 89 | Line 92 | namespace oopse { | |
92 | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | |
93 | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), | |
94 | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), | |
95 | < | sman_(NULL), fortranInitialized_(false) { |
95 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), |
96 | > | useAtomicVirial_(true) { |
97 | ||
98 | MoleculeStamp* molStamp; | |
99 | int nMolWithSameStamp; | |
# | Line 152 | Line 156 | namespace oopse { | |
156 | + nGlobalRigidBodies_; | |
157 | ||
158 | nGlobalMols_ = molStampIds_.size(); | |
155 | – | |
156 | – | #ifdef IS_MPI |
159 | molToProcMap_.resize(nGlobalMols_); | |
158 | – | #endif |
159 | – | |
160 | } | |
161 | ||
162 | SimInfo::~SimInfo() { | |
# | Line 600 | Line 600 | namespace oopse { | |
600 | /** @deprecate */ | |
601 | int isError = 0; | |
602 | ||
603 | + | setupCutoff(); |
604 | + | |
605 | setupElectrostaticSummationMethod( isError ); | |
606 | setupSwitchingFunction(); | |
607 | + | setupAccumulateBoxDipole(); |
608 | ||
609 | if(isError){ | |
610 | sprintf( painCave.errMsg, | |
# | Line 609 | Line 612 | namespace oopse { | |
612 | painCave.isFatal = 1; | |
613 | simError(); | |
614 | } | |
612 | – | |
613 | – | |
614 | – | setupCutoff(); |
615 | ||
616 | calcNdf(); | |
617 | calcNdfRaw(); | |
# | Line 661 | Line 661 | namespace oopse { | |
661 | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); | |
662 | int useRF; | |
663 | int useSF; | |
664 | + | int useSP; |
665 | + | int useBoxDipole; |
666 | + | |
667 | std::string myMethod; | |
668 | ||
669 | // set the useRF logical | |
670 | useRF = 0; | |
671 | useSF = 0; | |
672 | + | useSP = 0; |
673 | ||
674 | ||
675 | if (simParams_->haveElectrostaticSummationMethod()) { | |
676 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
677 | toUpper(myMethod); | |
678 | < | if (myMethod == "REACTION_FIELD") { |
679 | < | useRF=1; |
680 | < | } else { |
681 | < | if (myMethod == "SHIFTED_FORCE") { |
682 | < | useSF = 1; |
683 | < | } |
678 | > | if (myMethod == "REACTION_FIELD"){ |
679 | > | useRF = 1; |
680 | > | } else if (myMethod == "SHIFTED_FORCE"){ |
681 | > | useSF = 1; |
682 | > | } else if (myMethod == "SHIFTED_POTENTIAL"){ |
683 | > | useSP = 1; |
684 | } | |
685 | } | |
686 | + | |
687 | + | if (simParams_->haveAccumulateBoxDipole()) |
688 | + | if (simParams_->getAccumulateBoxDipole()) |
689 | + | useBoxDipole = 1; |
690 | ||
691 | + | useAtomicVirial_ = simParams_->getUseAtomicVirial(); |
692 | + | |
693 | //loop over all of the atom types | |
694 | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | |
695 | useLennardJones |= (*i)->isLennardJones(); | |
# | Line 749 | Line 759 | namespace oopse { | |
759 | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); | |
760 | ||
761 | temp = useSF; | |
762 | < | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
762 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
763 | ||
764 | + | temp = useSP; |
765 | + | MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
766 | + | |
767 | + | temp = useBoxDipole; |
768 | + | MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
769 | + | |
770 | + | temp = useAtomicVirial_; |
771 | + | MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
772 | + | |
773 | #endif | |
774 | ||
775 | fInfo_.SIM_uses_PBC = usePBC; | |
# | Line 768 | Line 787 | namespace oopse { | |
787 | fInfo_.SIM_uses_FLARB = useFLARB; | |
788 | fInfo_.SIM_uses_RF = useRF; | |
789 | fInfo_.SIM_uses_SF = useSF; | |
790 | < | |
791 | < | if( myMethod == "REACTION_FIELD") { |
792 | < | |
774 | < | if (simParams_->haveDielectric()) { |
775 | < | fInfo_.dielect = simParams_->getDielectric(); |
776 | < | } else { |
777 | < | sprintf(painCave.errMsg, |
778 | < | "SimSetup Error: No Dielectric constant was set.\n" |
779 | < | "\tYou are trying to use Reaction Field without" |
780 | < | "\tsetting a dielectric constant!\n"); |
781 | < | painCave.isFatal = 1; |
782 | < | simError(); |
783 | < | } |
784 | < | } |
785 | < | |
790 | > | fInfo_.SIM_uses_SP = useSP; |
791 | > | fInfo_.SIM_uses_BoxDipole = useBoxDipole; |
792 | > | fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_; |
793 | } | |
794 | ||
795 | void SimInfo::setupFortranSim() { | |
# | Line 849 | Line 856 | namespace oopse { | |
856 | int nGlobalExcludes = 0; | |
857 | int* globalExcludes = NULL; | |
858 | int* excludeList = exclude_.getExcludeList(); | |
859 | < | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList , |
860 | < | &nGlobalExcludes, globalExcludes, &molMembershipArray[0], |
861 | < | &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); |
862 | < | |
859 | > | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], |
860 | > | &nExclude, excludeList , &nGlobalExcludes, globalExcludes, |
861 | > | &molMembershipArray[0], &mfact[0], &nCutoffGroups_, |
862 | > | &fortranGlobalGroupMembership[0], &isError); |
863 | > | |
864 | if( isError ){ | |
865 | < | |
865 | > | |
866 | sprintf( painCave.errMsg, | |
867 | "There was an error setting the simulation information in fortran.\n" ); | |
868 | painCave.isFatal = 1; | |
869 | painCave.severity = OOPSE_ERROR; | |
870 | simError(); | |
871 | } | |
872 | < | |
873 | < | #ifdef IS_MPI |
872 | > | |
873 | > | |
874 | sprintf( checkPointMsg, | |
875 | "succesfully sent the simulation information to fortran.\n"); | |
876 | < | MPIcheckPoint(); |
877 | < | #endif // is_mpi |
876 | > | |
877 | > | errorCheckPoint(); |
878 | > | |
879 | > | // Setup number of neighbors in neighbor list if present |
880 | > | if (simParams_->haveNeighborListNeighbors()) { |
881 | > | int nlistNeighbors = simParams_->getNeighborListNeighbors(); |
882 | > | setNeighbors(&nlistNeighbors); |
883 | > | } |
884 | > | |
885 | > | |
886 | } | |
887 | ||
888 | ||
873 | – | #ifdef IS_MPI |
889 | void SimInfo::setupFortranParallel() { | |
890 | < | |
890 | > | #ifdef IS_MPI |
891 | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex | |
892 | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); | |
893 | std::vector<int> localToGlobalCutoffGroupIndex; | |
# | Line 922 | Line 937 | namespace oopse { | |
937 | } | |
938 | ||
939 | sprintf(checkPointMsg, " mpiRefresh successful.\n"); | |
940 | < | MPIcheckPoint(); |
940 | > | errorCheckPoint(); |
941 | ||
942 | < | |
942 | > | #endif |
943 | } | |
944 | ||
930 | – | #endif |
931 | – | |
945 | void SimInfo::setupCutoff() { | |
946 | ||
947 | ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); | |
# | Line 936 | Line 949 | namespace oopse { | |
949 | // Check the cutoff policy | |
950 | int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default | |
951 | ||
952 | + | // Set LJ shifting bools to false |
953 | + | ljsp_ = false; |
954 | + | ljsf_ = false; |
955 | + | |
956 | std::string myPolicy; | |
957 | if (forceFieldOptions_.haveCutoffPolicy()){ | |
958 | myPolicy = forceFieldOptions_.getCutoffPolicy(); | |
# | Line 999 | Line 1016 | namespace oopse { | |
1016 | simError(); | |
1017 | } | |
1018 | } | |
1019 | < | |
1020 | < | notifyFortranCutoffs(&rcut_, &rsw_); |
1019 | > | |
1020 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
1021 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
1022 | > | toUpper(myMethod); |
1023 | > | |
1024 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
1025 | > | ljsp_ = true; |
1026 | > | } else if (myMethod == "SHIFTED_FORCE") { |
1027 | > | ljsf_ = true; |
1028 | > | } |
1029 | > | } |
1030 | > | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1031 | ||
1032 | } else { | |
1033 | ||
# | Line 1017 | Line 1044 | namespace oopse { | |
1044 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1045 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
1046 | toUpper(myMethod); | |
1047 | < | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1047 | > | |
1048 | > | // For the time being, we're tethering the LJ shifted behavior to the |
1049 | > | // electrostaticSummationMethod keyword options |
1050 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
1051 | > | ljsp_ = true; |
1052 | > | } else if (myMethod == "SHIFTED_FORCE") { |
1053 | > | ljsf_ = true; |
1054 | > | } |
1055 | > | if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { |
1056 | if (simParams_->haveSwitchingRadius()){ | |
1057 | sprintf(painCave.errMsg, | |
1058 | "SimInfo Warning: A value was set for the switchingRadius\n" | |
# | Line 1040 | Line 1075 | namespace oopse { | |
1075 | simError(); | |
1076 | rsw_ = 0.85 * rcut_; | |
1077 | } | |
1078 | < | notifyFortranCutoffs(&rcut_, &rsw_); |
1078 | > | |
1079 | > | notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_); |
1080 | > | |
1081 | } else { | |
1082 | // We didn't set rcut explicitly, and we don't have electrostatic atoms, so | |
1083 | // We'll punt and let fortran figure out the cutoffs later. | |
# | Line 1058 | Line 1095 | namespace oopse { | |
1095 | int sm = UNDAMPED; | |
1096 | RealType alphaVal; | |
1097 | RealType dielectric; | |
1098 | < | |
1098 | > | |
1099 | errorOut = isError; | |
1063 | – | alphaVal = simParams_->getDampingAlpha(); |
1064 | – | dielectric = simParams_->getDielectric(); |
1100 | ||
1101 | if (simParams_->haveElectrostaticSummationMethod()) { | |
1102 | std::string myMethod = simParams_->getElectrostaticSummationMethod(); | |
# | Line 1078 | Line 1113 | namespace oopse { | |
1113 | if (myMethod == "SHIFTED_FORCE") { | |
1114 | esm = SHIFTED_FORCE; | |
1115 | } else { | |
1116 | < | if (myMethod == "REACTION_FIELD") { |
1116 | > | if (myMethod == "REACTION_FIELD") { |
1117 | esm = REACTION_FIELD; | |
1118 | + | dielectric = simParams_->getDielectric(); |
1119 | + | if (!simParams_->haveDielectric()) { |
1120 | + | // throw warning |
1121 | + | sprintf( painCave.errMsg, |
1122 | + | "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n" |
1123 | + | "\tA default value of %f will be used for the dielectric.\n", dielectric); |
1124 | + | painCave.isFatal = 0; |
1125 | + | simError(); |
1126 | + | } |
1127 | } else { | |
1128 | // throw error | |
1129 | sprintf( painCave.errMsg, | |
# | Line 1106 | Line 1150 | namespace oopse { | |
1150 | if (myScreen == "DAMPED") { | |
1151 | sm = DAMPED; | |
1152 | if (!simParams_->haveDampingAlpha()) { | |
1153 | < | //throw error |
1153 | > | // first set a cutoff dependent alpha value |
1154 | > | // we assume alpha depends linearly with rcut from 0 to 20.5 ang |
1155 | > | alphaVal = 0.5125 - rcut_* 0.025; |
1156 | > | // for values rcut > 20.5, alpha is zero |
1157 | > | if (alphaVal < 0) alphaVal = 0; |
1158 | > | |
1159 | > | // throw warning |
1160 | sprintf( painCave.errMsg, | |
1161 | "SimInfo warning: dampingAlpha was not specified in the input file.\n" | |
1162 | < | "\tA default value of %f (1/ang) will be used.\n", alphaVal); |
1162 | > | "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_); |
1163 | painCave.isFatal = 0; | |
1164 | simError(); | |
1165 | + | } else { |
1166 | + | alphaVal = simParams_->getDampingAlpha(); |
1167 | } | |
1168 | + | |
1169 | } else { | |
1170 | // throw error | |
1171 | sprintf( painCave.errMsg, | |
# | Line 1158 | Line 1211 | namespace oopse { | |
1211 | ||
1212 | // send switching function notification to switcheroo | |
1213 | setFunctionType(&ft); | |
1214 | + | |
1215 | + | } |
1216 | + | |
1217 | + | void SimInfo::setupAccumulateBoxDipole() { |
1218 | + | |
1219 | + | // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true |
1220 | + | if ( simParams_->haveAccumulateBoxDipole() ) |
1221 | + | if ( simParams_->getAccumulateBoxDipole() ) { |
1222 | + | setAccumulateBoxDipole(); |
1223 | + | calcBoxDipole_ = true; |
1224 | + | } |
1225 | ||
1226 | } | |
1227 | ||
# | Line 1410 | Line 1474 | namespace oopse { | |
1474 | return angularMomentum; | |
1475 | } | |
1476 | ||
1477 | < | |
1477 | > | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { |
1478 | > | return IOIndexToIntegrableObject.at(index); |
1479 | > | } |
1480 | > | |
1481 | > | void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) { |
1482 | > | IOIndexToIntegrableObject= v; |
1483 | > | } |
1484 | > | |
1485 | > | /* Returns the Volume of the simulation based on a ellipsoid with semi-axes |
1486 | > | based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 |
1487 | > | where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to |
1488 | > | V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. |
1489 | > | */ |
1490 | > | void SimInfo::getGyrationalVolume(RealType &volume){ |
1491 | > | Mat3x3d intTensor; |
1492 | > | RealType det; |
1493 | > | Vector3d dummyAngMom; |
1494 | > | RealType sysconstants; |
1495 | > | RealType geomCnst; |
1496 | > | |
1497 | > | geomCnst = 3.0/2.0; |
1498 | > | /* Get the inertial tensor and angular momentum for free*/ |
1499 | > | getInertiaTensor(intTensor,dummyAngMom); |
1500 | > | |
1501 | > | det = intTensor.determinant(); |
1502 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1503 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det); |
1504 | > | return; |
1505 | > | } |
1506 | > | |
1507 | > | void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ |
1508 | > | Mat3x3d intTensor; |
1509 | > | Vector3d dummyAngMom; |
1510 | > | RealType sysconstants; |
1511 | > | RealType geomCnst; |
1512 | > | |
1513 | > | geomCnst = 3.0/2.0; |
1514 | > | /* Get the inertial tensor and angular momentum for free*/ |
1515 | > | getInertiaTensor(intTensor,dummyAngMom); |
1516 | > | |
1517 | > | detI = intTensor.determinant(); |
1518 | > | sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; |
1519 | > | volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI); |
1520 | > | return; |
1521 | > | } |
1522 | > | /* |
1523 | > | void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) { |
1524 | > | assert( v.size() == nAtoms_ + nRigidBodies_); |
1525 | > | sdByGlobalIndex_ = v; |
1526 | > | } |
1527 | > | |
1528 | > | StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { |
1529 | > | //assert(index < nAtoms_ + nRigidBodies_); |
1530 | > | return sdByGlobalIndex_.at(index); |
1531 | > | } |
1532 | > | */ |
1533 | }//end namespace oopse | |
1534 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |