| 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" |
| 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" |
| 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; |
| 156 |
|
+ nGlobalRigidBodies_; |
| 157 |
|
|
| 158 |
|
nGlobalMols_ = molStampIds_.size(); |
| 156 |
– |
|
| 157 |
– |
#ifdef IS_MPI |
| 159 |
|
molToProcMap_.resize(nGlobalMols_); |
| 159 |
– |
#endif |
| 160 |
– |
|
| 160 |
|
} |
| 161 |
|
|
| 162 |
|
SimInfo::~SimInfo() { |
| 663 |
|
int useSF; |
| 664 |
|
int useSP; |
| 665 |
|
int useBoxDipole; |
| 666 |
+ |
|
| 667 |
|
std::string myMethod; |
| 668 |
|
|
| 669 |
|
// set the useRF logical |
| 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(); |
| 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; |
| 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() { |
| 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 |
|
|
| 874 |
– |
#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; |
| 937 |
|
} |
| 938 |
|
|
| 939 |
|
sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
| 940 |
< |
MPIcheckPoint(); |
| 940 |
> |
errorCheckPoint(); |
| 941 |
|
|
| 942 |
< |
|
| 942 |
> |
#endif |
| 943 |
|
} |
| 944 |
|
|
| 931 |
– |
#endif |
| 932 |
– |
|
| 945 |
|
void SimInfo::setupCutoff() { |
| 946 |
|
|
| 947 |
|
ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); |
| 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(); |
| 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 |
|
|
| 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" |
| 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. |
| 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_); |