# | Line 59 | Line 59 | |
---|---|---|
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 67 | 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 90 | 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), calcBoxDipole_(false) { |
95 | > | sman_(NULL), fortranInitialized_(false), calcBoxDipole_(false), |
96 | > | useAtomicVirial_(true) { |
97 | ||
98 | MoleculeStamp* molStamp; | |
99 | int nMolWithSameStamp; | |
# | Line 153 | Line 156 | namespace oopse { | |
156 | + nGlobalRigidBodies_; | |
157 | ||
158 | nGlobalMols_ = molStampIds_.size(); | |
156 | – | |
157 | – | #ifdef IS_MPI |
159 | molToProcMap_.resize(nGlobalMols_); | |
159 | – | #endif |
160 | – | |
160 | } | |
161 | ||
162 | SimInfo::~SimInfo() { | |
# | Line 664 | Line 663 | namespace oopse { | |
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; |
679 | > | useRF = 1; |
680 | } else if (myMethod == "SHIFTED_FORCE"){ | |
681 | useSF = 1; | |
682 | } else if (myMethod == "SHIFTED_POTENTIAL"){ | |
# | Line 687 | Line 688 | namespace oopse { | |
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 763 | Line 766 | namespace oopse { | |
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 | ||
# | Line 783 | Line 789 | namespace oopse { | |
789 | fInfo_.SIM_uses_SF = useSF; | |
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 1125 | Line 1162 | namespace oopse { | |
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 1442 | Line 1482 | namespace oopse { | |
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_); |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |