Skip to content

GitLab

Unverified Commit 9b523076 authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2301 from peastman/hippo2 

Improvements to HippoNonbondedForce
parents 43f619bc c60568b2
Hide whitespace changes
Inline Side-by-side
Showing with 342 additions and 137 deletions
plugins/amoeba/openmmapi/include/openmm/HippoNonbondedForce.h
View file @ 9b523076
......@@ -72,7 +72,7 @@ public:
PME = 1
};
enum ParticleAxisTypes { ZThenX = 0, Bisector = 1, ZBisect = 2, ThreeFold = 3, ZOnly = 4, NoAxisType = 5, LastAxisTypeIndex = 6 };
enum ParticleAxisTypes { ZThenX = 0, Bisector = 1, ZBisect = 2, ThreeFold = 3, ZOnly = 4, NoAxisType = 5 };
/**
* Create a HippoNonbondedForce.
......@@ -287,12 +287,13 @@ public:
* @param dipoleDipoleScale the factor by which to scale the Coulomb interaction between induced dipoles
* @param dispersionScale the factor by which to scale the dispersion interaction
* @param repulsionScale the factor by which to scale the Pauli repulsion
* @param chargeTransferScale the factor by which to scale the charge transfer interaction
* @param replace determines the behavior if there is already an exception for the same two particles. If true, the existing one is replaced.
* If false, an exception is thrown.
* @return the index of the exception that was added
*/
int addException(int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale,
double dispersionScale, double repulsionScale, bool replace = false);
double dispersionScale, double repulsionScale, double chargeTransferScale, bool replace = false);
/**
* Get the scale factors for an interaction that should be calculated differently from others.
*
......@@ -304,9 +305,10 @@ public:
* @param dipoleDipoleScale the factor by which to scale the Coulomb interaction between induced dipoles
* @param dispersionScale the factor by which to scale the dispersion interaction
* @param repulsionScale the factor by which to scale the Pauli repulsion
* @param chargeTransferScale the factor by which to scale the charge transfer interaction
*/
void getExceptionParameters(int index, int& particle1, int& particle2, double& multipoleMultipoleScale, double& dipoleMultipoleScale, double& dipoleDipoleScale,
double& dispersionScale, double& repulsionScale) const;
double& dispersionScale, double& repulsionScale, double& chargeTransferScale) const;
/**
* Set the scale factors for an interaction that should be calculated differently from others.
* If all scale factors are set to 0, this will cause the interaction to be completely omitted from
......@@ -320,9 +322,10 @@ public:
* @param dipoleDipoleScale the factor by which to scale the Coulomb interaction between induced dipoles
* @param dispersionScale the factor by which to scale the dispersion interaction
* @param repulsionScale the factor by which to scale the Pauli repulsion
* @param chargeTransferScale the factor by which to scale the charge transfer interaction
*/
void setExceptionParameters(int index, int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale,
double dispersionScale, double repulsionScale);
double dispersionScale, double repulsionScale, double chargeTransferScale);
/**
* Get the error tolerance for Ewald summation. This corresponds to the fractional error in the forces
* which is acceptable. This value is used to select the grid dimensions and separation (alpha)
......@@ -424,14 +427,14 @@ public:
class HippoNonbondedForce::ExceptionInfo {
public:
int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
ExceptionInfo() {
particle1 = particle2 = -1;
multipoleMultipoleScale = dipoleMultipoleScale = dipoleDipoleScale = dispersionScale = repulsionScale = 0.0;
multipoleMultipoleScale = dipoleMultipoleScale = dipoleDipoleScale = dispersionScale = repulsionScale = chargeTransferScale = 0.0;
}
ExceptionInfo(int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale, double dispersionScale, double repulsionScale) :
ExceptionInfo(int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale, double dispersionScale, double repulsionScale, double chargeTransferScale) :
particle1(particle1), particle2(particle2), multipoleMultipoleScale(multipoleMultipoleScale), dipoleMultipoleScale(dipoleMultipoleScale),
dipoleDipoleScale(dipoleDipoleScale), dispersionScale(dispersionScale), repulsionScale(repulsionScale) {
dipoleDipoleScale(dipoleDipoleScale), dispersionScale(dispersionScale), repulsionScale(repulsionScale), chargeTransferScale(chargeTransferScale) {
}
};
......
plugins/amoeba/openmmapi/src/HippoNonbondedForce.cpp
View file @ 9b523076
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2018 Stanford University and the Authors. *
* Portions copyright (c) 2008-2019 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
......@@ -173,7 +173,7 @@ void HippoNonbondedForce::setParticleParameters(int index, double charge, const
}
int HippoNonbondedForce::addException(int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale,
double dispersionScale, double repulsionScale, bool replace) {
double dispersionScale, double repulsionScale, double chargeTransferScale, bool replace) {
map<pair<int, int>, int>::iterator iter = exceptionMap.find(pair<int, int>(particle1, particle2));
int newIndex;
if (iter == exceptionMap.end())
......@@ -187,19 +187,19 @@ int HippoNonbondedForce::addException(int particle1, int particle2, double multi
msg << particle2;
throw OpenMMException(msg.str());
}
exceptions[iter->second] = ExceptionInfo(particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale);
exceptions[iter->second] = ExceptionInfo(particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
newIndex = iter->second;
exceptionMap.erase(iter->first);
}
else {
exceptions.push_back(ExceptionInfo(particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale));
exceptions.push_back(ExceptionInfo(particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale));
newIndex = exceptions.size()-1;
}
exceptionMap[pair<int, int>(particle1, particle2)] = newIndex;
return newIndex;
}
void HippoNonbondedForce::getExceptionParameters(int index, int& particle1, int& particle2, double& multipoleMultipoleScale, double& dipoleMultipoleScale, double& dipoleDipoleScale,
double& dispersionScale, double& repulsionScale) const {
double& dispersionScale, double& repulsionScale, double& chargeTransferScale) const {
ASSERT_VALID_INDEX(index, exceptions);
particle1 = exceptions[index].particle1;
particle2 = exceptions[index].particle2;
......@@ -208,10 +208,11 @@ void HippoNonbondedForce::getExceptionParameters(int index, int& particle1, int&
dipoleDipoleScale = exceptions[index].dipoleDipoleScale;
dispersionScale = exceptions[index].dispersionScale;
repulsionScale = exceptions[index].repulsionScale;
chargeTransferScale = exceptions[index].chargeTransferScale;
}
void HippoNonbondedForce::setExceptionParameters(int index, int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale,
double dispersionScale, double repulsionScale) {
double dispersionScale, double repulsionScale, double chargeTransferScale) {
ASSERT_VALID_INDEX(index, exceptions);
exceptions[index].particle1 = particle1;
exceptions[index].particle2 = particle2;
......@@ -220,6 +221,7 @@ void HippoNonbondedForce::setExceptionParameters(int index, int particle1, int p
exceptions[index].dipoleDipoleScale = dipoleDipoleScale;
exceptions[index].dispersionScale = dispersionScale;
exceptions[index].repulsionScale = repulsionScale;
exceptions[index].chargeTransferScale = chargeTransferScale;
}
void HippoNonbondedForce::getInducedDipoles(Context& context, vector<Vec3>& dipoles) {
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 9b523076
......@@ -2624,20 +2624,20 @@ public:
}
void getParticlesInGroup(int index, vector<int>& particles) {
int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale;
force.getExceptionParameters(index, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale);
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
force.getExceptionParameters(index, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
particles.resize(2);
particles[0] = particle1;
particles[1] = particle2;
}
bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2;
double multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1;
double multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2;
force.getExceptionParameters(group1, particle1, particle2, multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1);
force.getExceptionParameters(group2, particle1, particle2, multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2);
double multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1, chargeTransferScale1;
double multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2, chargeTransferScale2;
force.getExceptionParameters(group1, particle1, particle2, multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1, chargeTransferScale1);
force.getExceptionParameters(group2, particle1, particle2, multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2, chargeTransferScale2);
return (multipoleMultipoleScale1 == multipoleMultipoleScale2 && dipoleMultipoleScale1 == dipoleMultipoleScale2 &&
dipoleDipoleScale1 == dipoleDipoleScale2 && dispersionScale1 == dispersionScale2 && repulsionScale1 == repulsionScale2);
dipoleDipoleScale1 == dipoleDipoleScale2 && dispersionScale1 == dispersionScale2 && repulsionScale1 == repulsionScale2 && chargeTransferScale1 == chargeTransferScale2);
}
private:
const HippoNonbondedForce& force;
......@@ -2778,27 +2778,28 @@ void CudaCalcHippoNonbondedForceKernel::initialize(const System& system, const H
// Record exceptions and exclusions.
vector<double> exceptionScaleVec[5];
vector<double> exceptionScaleVec[6];
vector<int2> exceptionAtomsVec;
for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale;
force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale);
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
exclusions[particle1].push_back(particle2);
exclusions[particle2].push_back(particle1);
if (usePME || multipoleMultipoleScale != 0 || dipoleMultipoleScale != 0 || dipoleDipoleScale != 0 || dispersionScale != 0 || repulsionScale != 0) {
if (usePME || multipoleMultipoleScale != 0 || dipoleMultipoleScale != 0 || dipoleDipoleScale != 0 || dispersionScale != 0 || repulsionScale != 0 || chargeTransferScale != 0) {
exceptionAtomsVec.push_back(make_int2(particle1, particle2));
exceptionScaleVec[0].push_back(multipoleMultipoleScale);
exceptionScaleVec[1].push_back(dipoleMultipoleScale);
exceptionScaleVec[2].push_back(dipoleDipoleScale);
exceptionScaleVec[3].push_back(dispersionScale);
exceptionScaleVec[4].push_back(repulsionScale);
exceptionScaleVec[5].push_back(chargeTransferScale);
}
}
if (exceptionAtomsVec.size() > 0) {
exceptionAtoms.initialize<int2>(cu, exceptionAtomsVec.size(), "exceptionAtoms");
exceptionAtoms.upload(exceptionAtomsVec);
for (int i = 0; i < 5; i++) {
for (int i = 0; i < 6; i++) {
exceptionScales[i].initialize(cu, exceptionAtomsVec.size(), elementSize, "exceptionScales");
exceptionScales[i].upload(exceptionScaleVec[i], true);
}
......@@ -3288,6 +3289,7 @@ double CudaCalcHippoNonbondedForceKernel::execute(ContextImpl& context, bool inc
computeExceptionsArgs.push_back(&exceptionScales[2].getDevicePointer());
computeExceptionsArgs.push_back(&exceptionScales[3].getDevicePointer());
computeExceptionsArgs.push_back(&exceptionScales[4].getDevicePointer());
computeExceptionsArgs.push_back(&exceptionScales[5].getDevicePointer());
computeExceptionsArgs.push_back(&coreCharge.getDevicePointer());
computeExceptionsArgs.push_back(&valenceCharge.getDevicePointer());
computeExceptionsArgs.push_back(&alpha.getDevicePointer());
......@@ -3697,26 +3699,27 @@ void CudaCalcHippoNonbondedForceKernel::copyParametersToContext(ContextImpl& con
// Record the per-exception parameters.
vector<double> exceptionScaleVec[5];
vector<double> exceptionScaleVec[6];
vector<int2> exceptionAtomsVec;
for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale;
force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale);
if (usePME || multipoleMultipoleScale != 0 || dipoleMultipoleScale != 0 || dipoleDipoleScale != 0 || dispersionScale != 0 || repulsionScale != 0) {
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
if (usePME || multipoleMultipoleScale != 0 || dipoleMultipoleScale != 0 || dipoleDipoleScale != 0 || dispersionScale != 0 || repulsionScale != 0 || chargeTransferScale != 0) {
exceptionAtomsVec.push_back(make_int2(particle1, particle2));
exceptionScaleVec[0].push_back(multipoleMultipoleScale);
exceptionScaleVec[1].push_back(dipoleMultipoleScale);
exceptionScaleVec[2].push_back(dipoleDipoleScale);
exceptionScaleVec[3].push_back(dispersionScale);
exceptionScaleVec[4].push_back(repulsionScale);
exceptionScaleVec[5].push_back(chargeTransferScale);
}
}
if (exceptionAtomsVec.size() > 0) {
if (!exceptionAtoms.isInitialized() || exceptionAtoms.getSize() != exceptionAtomsVec.size())
throw OpenMMException("updateParametersInContext: The number of exceptions has changed");
exceptionAtoms.upload(exceptionAtomsVec);
for (int i = 0; i < 5; i++)
for (int i = 0; i < 6; i++)
exceptionScales[i].upload(exceptionScaleVec[i], true);
}
else if (exceptionAtoms.isInitialized())
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h
View file @ 9b523076
......@@ -733,7 +733,7 @@ private:
CudaArray dpmeBsplineModuliX, dpmeBsplineModuliY, dpmeBsplineModuliZ;
CudaArray pmePhi, pmePhidp, pmeCphi;
CudaArray lastPositions;
CudaArray exceptionScales[5];
CudaArray exceptionScales[6];
CudaArray exceptionAtoms;
CudaSort* sort;
cufftHandle fftForward, fftBackward, dfftForward, dfftBackward;
......
plugins/amoeba/platforms/cuda/src/kernels/hippoComputeField.cu
View file @ 9b523076
......@@ -32,16 +32,16 @@ __device__ void computeMutualFieldDampingFactors(real alphaI, real alphaJ, real
real aJ2 = alphaJ*alphaJ;
real A = aJ2/(aJ2-aI2);
real B = aI2/(aI2-aJ2);
real A2 = A*A;
real B2 = B*B;
fdamp3 = 1 - A2*(1 + arI + arI2*(one/2))*expARI -
B2*(1 + arJ + arJ2*(one/2))*expARJ -
2*A2*B*(1 + arI)*expARI -
2*B2*A*(1 + arJ)*expARJ;
fdamp5 = 1 - A2*(1 + arI + arI2*(one/2) + arI3*(one/6))*expARI -
B2*(1 + arJ + arJ2*(one/2) + arJ3*(one/6))*expARJ -
2*A2*B*(1 + arI + arI2*(one/3))*expARI -
2*B2*A*(1 + arJ + arJ2*(one/3))*expARJ;
real A2expARI = A*A*expARI;
real B2expARJ = B*B*expARJ;
fdamp3 = 1 - (1 + arI + arI2*(one/2))*A2expARI -
(1 + arJ + arJ2*(one/2))*B2expARJ -
(1 + arI)*2*B*A2expARI -
(1 + arJ)*2*A*B2expARJ;
fdamp5 = 1 - (1 + arI + arI2*(one/2) + arI3*(one/6))*A2expARI -
(1 + arJ + arJ2*(one/2) + arJ3*(one/6))*B2expARJ -
(1 + arI + arI2*(one/3))*2*B*A2expARI -
(1 + arJ + arJ2*(one/3))*2*A*B2expARJ;
}
}
......
plugins/amoeba/platforms/cuda/src/kernels/hippoFixedField.cu
View file @ 9b523076
......@@ -7,10 +7,20 @@ real invR7 = invR5*invR2;
// Calculate the error function damping terms.
real ralpha = PME_ALPHA*r;
real bn0 = erfc(ralpha)*invR;
real exp2a = EXP(-(ralpha*ralpha));
#ifdef USE_DOUBLE_PRECISION
const real erfcAlphaR = erfc(ralpha);
#else
// This approximation for erfc is from Abramowitz and Stegun (1964) p. 299. They cite the following as
// the original source: C. Hastings, Jr., Approximations for Digital Computers (1955). It has a maximum
// error of 1.5e-7.
const real t = RECIP(1.0f+0.3275911f*ralpha);
const real erfcAlphaR = (0.254829592f+(-0.284496736f+(1.421413741f+(-1.453152027f+1.061405429f*t)*t)*t)*t)*t*exp2a;
#endif
real bn0 = erfcAlphaR*invR;
real alsq2 = 2*PME_ALPHA*PME_ALPHA;
real alsq2n = 1/(SQRT_PI*PME_ALPHA);
real exp2a = EXP(-(ralpha*ralpha));
alsq2n *= alsq2;
real bn1 = (bn0+alsq2n*exp2a)*invR2;
alsq2n *= alsq2;
......
plugins/amoeba/platforms/cuda/src/kernels/hippoInteraction.cu
View file @ 9b523076
......@@ -168,7 +168,7 @@ real bn5 = (9*bn4+alsq2n*exp2a)*rInv2;
// Compute the force and torque.
real3 elecForce = -scale*(de*make_real3(0, 0, r) + term1*qiDipole1 + term2*qiDipole2 +
real3 elecForce = -scale*(make_real3(0, 0, de*r) + term1*qiDipole1 + term2*qiDipole2 +
term3*(diqkTemp-dkqiTemp) + term4*qi + term5*qk + term6*(qikTemp+qkiTemp));
real3 tI = scale*(-rr3ik*dikCross + term1*dirCross + term3*(dqik+dkqirCross) + term4*qirCross - term6*(qikrCross+qikCross));
real3 tK = scale*(rr3ik*dikCross + term2*dkrCross - term3*(dqik+diqkrCross) + term5*qkrCross - term6*(qkirCross-qikCross));
......@@ -423,8 +423,8 @@ real bn5 = (9*bn4+alsq2n*exp2a)*rInv2;
}
#endif
#ifdef COMPUTING_EXCEPTIONS
ctForce *= multipoleMultipoleScale;
ctEnergy *= multipoleMultipoleScale;
ctForce *= chargeTransferScale;
ctEnergy *= chargeTransferScale;
#endif
tempEnergy += includeInteraction ? ctEnergy : 0;
tempForce.z += includeInteraction ? ctForce*r : 0;
......
plugins/amoeba/platforms/cuda/src/kernels/hippoInteractionHeader.cu
View file @ 9b523076
......@@ -96,35 +96,35 @@ __device__ void computeOverlapDampingFactors(real alphaI, real alphaJ, real r,
real aJ2 = alphaJ*alphaJ;
real A = aJ2/(aJ2-aI2);
real B = aI2/(aI2-aJ2);
real A2 = A*A;
real B2 = B*B;
real A2expARI = A*A*expARI;
real B2expARJ = B*B*expARJ;
fdampJ1 = 1 - (1 + arJ*(one/2))*expARJ;
fdampJ3 = 1 - (1 + arJ + arJ2*(one/2))*expARJ;
fdampJ5 = 1 - (1 + arJ + arJ2*(one/2) + arJ3*(one/6))*expARJ;
fdampJ7 = 1 - (1 + arJ + arJ2*(one/2) + arJ3*(one/6) + arJ4*(one/30))*expARJ;
fdampJ9 = 1 - (1 + arJ + arJ2*(one/2) + arJ3*(one/6) + 4*arJ4*(one/105) + arJ5*(one/210))*expARJ;
fdampIJ1 = 1 - A2*(1 + 2*B + arI*(one/2))*expARI -
B2*(1 + 2*A + arJ*(one/2))*expARJ;
fdampIJ3 = 1 - A2*(1 + arI + arI2*(one/2))*expARI -
B2*(1 + arJ + arJ2*(one/2))*expARJ -
2*A2*B*(1 + arI)*expARI -
2*B2*A*(1 + arJ)*expARJ;
fdampIJ5 = 1 - A2*(1 + arI + arI2*(one/2) + arI3*(one/6))*expARI -
B2*(1 + arJ + arJ2*(one/2) + arJ3*(one/6))*expARJ -
2*A2*B*(1 + arI + arI2*(one/3))*expARI -
2*B2*A*(1 + arJ + arJ2*(one/3))*expARJ;
fdampIJ7 = 1 - A2*(1 + arI + arI2*(one/2) + arI3*(one/6) + arI4*(one/30))*expARI -
B2*(1 + arJ + arJ2*(one/2) + arJ3*(one/6) + arJ4*(one/30))*expARJ -
2*A2*B*(1 + arI + arI2*(two/5) + arI3*(one/15))*expARI -
2*B2*A*(1 + arJ + arJ2*(two/5) + arJ3*(one/15))*expARJ;
fdampIJ9 = 1 - A2*(1 + arI + arI2*(one/2) + arI3*(one/6) + arI4*4*(one/105) + arI5*(one/210))*expARI -
B2*(1 + arJ + arJ2*(one/2) + arJ3*(one/6) + arJ4*4*(one/105) + arJ5*(one/210))*expARJ -
2*A2*B*(1 + arI + arI2*(three/7) + arI3*(two/21) + arI4*(one/105))*expARI -
2*B2*A*(1 + arJ + arJ2*(three/7) + arJ3*(two/21) + arJ4*(one/105))*expARJ;
fdampIJ11 = 1 - A2*(1 + arI + arI2*(one/2) + arI3*(one/6) + arI4*(five/126) + arI5*(two/315) + arI6*(one/1890))*expARI -
B2*(1 + arJ + arJ2*(one/2) + arJ3*(one/6) + arJ4*(five/126) + arJ5*(two/315) + arJ6*(one/1890))*expARJ -
2*A2*B*(1 + arI + arI2*(four/9) + arI3*(one/9) + arI4*(one/63) + arI5*(one/945))*expARI -
2*B2*A*(1 + arJ + arJ2*(four/9) + arJ3*(one/9) + arJ4*(one/63) + arJ5*(one/945))*expARJ;
fdampIJ1 = 1 - (1 + 2*B + arI*(one/2))*A2expARI -
(1 + 2*A + arJ*(one/2))*B2expARJ;
fdampIJ3 = 1 - (1 + arI + arI2*(one/2))*A2expARI -
(1 + arJ + arJ2*(one/2))*B2expARJ -
(1 + arI)*2*B*A2expARI -
(1 + arJ)*2*A*B2expARJ;
fdampIJ5 = 1 - (1 + arI + arI2*(one/2) + arI3*(one/6))*A2expARI -
(1 + arJ + arJ2*(one/2) + arJ3*(one/6))*B2expARJ -
(1 + arI + arI2*(one/3))*2*B*A2expARI -
(1 + arJ + arJ2*(one/3))*2*A*B2expARJ;
fdampIJ7 = 1 - (1 + arI + arI2*(one/2) + arI3*(one/6) + arI4*(one/30))*A2expARI -
(1 + arJ + arJ2*(one/2) + arJ3*(one/6) + arJ4*(one/30))*B2expARJ -
(1 + arI + arI2*(two/5) + arI3*(one/15))*2*B*A2expARI -
(1 + arJ + arJ2*(two/5) + arJ3*(one/15))*2*A*B2expARJ;
fdampIJ9 = 1 - (1 + arI + arI2*(one/2) + arI3*(one/6) + arI4*4*(one/105) + arI5*(one/210))*A2expARI -
(1 + arJ + arJ2*(one/2) + arJ3*(one/6) + arJ4*4*(one/105) + arJ5*(one/210))*B2expARJ -
(1 + arI + arI2*(three/7) + arI3*(two/21) + arI4*(one/105))*2*B*A2expARI -
(1 + arJ + arJ2*(three/7) + arJ3*(two/21) + arJ4*(one/105))*2*A*B2expARJ;
fdampIJ11 = 1 - (1 + arI + arI2*(one/2) + arI3*(one/6) + arI4*(five/126) + arI5*(two/315) + arI6*(one/1890))*A2expARI -
(1 + arJ + arJ2*(one/2) + arJ3*(one/6) + arJ4*(five/126) + arJ5*(two/315) + arJ6*(one/1890))*B2expARJ -
(1 + arI + arI2*(four/9) + arI3*(one/9) + arI4*(one/63) + arI5*(one/945))*2*B*A2expARI -
(1 + arJ + arJ2*(four/9) + arJ3*(one/9) + arJ4*(one/63) + arJ5*(one/945))*2*A*B2expARJ;
}
}
......@@ -152,18 +152,18 @@ __device__ void computeDispersionDampingFactors(real alphaI, real alphaJ, real r
real aJ2 = alphaJ*alphaJ;
real A = aJ2/(aJ2-aI2);
real B = aI2/(aI2-aJ2);
real A2 = A*A;
real B2 = B*B;
fdamp3 = 1 - A2*(1 + arI + arI2*(one/2))*expARI -
B2*(1 + arJ + arJ2*(one/2))*expARJ -
2*A2*B*(1 + arI)*expARI -
2*B2*A*(1 + arJ)*expARJ;
fdamp5 = 1 - A2*(1 + arI + arI2*(one/2) + arI3*(one/6))*expARI -
B2*(1 + arJ + arJ2*(one/2) + arJ3*(one/6))*expARJ -
2*A2*B*(1 + arI + arI2*(one/3))*expARI -
2*B2*A*(1 + arJ + arJ2*(one/3))*expARJ;
ddamp = (A2*arI2*alphaI*expARI*(r*alphaI + 4*B - 1) +
(B2*arJ2*alphaJ*expARJ*(r*alphaJ + 4*A - 1)))*(one/4);
real A2expARI = A*A*expARI;
real B2expARJ = B*B*expARJ;
fdamp3 = 1 - (1 + arI + arI2*(one/2))*A2expARI -
(1 + arJ + arJ2*(one/2))*B2expARJ -
(1 + arI)*2*B*A2expARI -
(1 + arJ)*2*A*B2expARJ;
fdamp5 = 1 - (1 + arI + arI2*(one/2) + arI3*(one/6))*A2expARI -
(1 + arJ + arJ2*(one/2) + arJ3*(one/6))*B2expARJ -
(1 + arI + arI2*(one/3))*2*B*A2expARI -
(1 + arJ + arJ2*(one/3))*2*A*B2expARJ;
ddamp = (arI2*alphaI*A2expARI*(r*alphaI + 4*B - 1) +
(arJ2*alphaJ*B2expARJ*(r*alphaJ + 4*A - 1)))*(one/4);
}
fdamp = 1.5f*fdamp5 - 0.5f*fdamp3;
}
......@@ -210,21 +210,22 @@ __device__ void computeRepulsionDampingFactors(real pauliAlphaI, real pauliAlpha
real aJ2_3 = aJ2_2*aJ2;
real aJ2_4 = aJ2_2*aJ2_2;
real aJ2_5 = aJ2_3*aJ2_2;
real aJ2_6 = aJ2_3*aJ2_3;
real scale = 1/(aI2_2-aJ2_2);
real aI2aJ2expI = aI2*aJ2*expI;
real aI2aJ2expJ = aI2*aJ2*expJ;
pre = 8192*aI2_3*aJ2_3*(scale*scale*scale*scale);
real tmp = 4*aI2*aJ2*scale;
fexp = (arI2-tmp)*expJ + (arJ2+tmp)*expI;
fexp1 = (aI2*aJ2*r2 - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2*aJ2*r2 + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp2 = (aI2*aJ2*r2*(one/3) + aI2*aJ2_2*r3*(one/3) - (four/3)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2*aJ2*r2*(one/3) + aI2_2*aJ2*r3*(one/3) + (four/3)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp3 = (aI2*aJ2_3*r4*(one/15) + aI2*aJ2_2*r3*(one/5) + aI2*aJ2*r2*(one/5) - (four/15)*aI2*aJ2_4*r3*scale - (eight/5)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*scale*aI2*aJ2)*expJ +
(aI2_3*aJ2*r4*(one/15) + aI2_2*aJ2*r3*(one/5) + aI2*aJ2*r2*(one/5) + (four/15)*aI2_4*aJ2*r3*scale + (eight/5)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*scale*aI2*aJ2)*expI;
fexp4 = (aI2*aJ2_4*r5*(one/105) + (two/35)*aI2*aJ2_3*r4 + aI2*aJ2_2*r3*(one/7) + aI2*aJ2*r2*(one/7) - (four/105)*aI2*aJ2_5*r4*scale - (eight/21)*aI2*aJ2_4*r3*scale - (twelve/7)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2_4*aJ2*r5*(one/105) + (two/35)*aI2_3*aJ2*r4 + aI2_2*aJ2*r3*(one/7) + aI2*aJ2*r2*(one/7) + (four/105)*aI2_5*aJ2*r4*scale + (eight/21)*aI2_4*aJ2*r3*scale + (twelve/7)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp5 = (aI2*aJ2_5*r6*(one/945) + (two/189)*aI2*aJ2_4*r5 + aI2*aJ2_3*r4*(one/21) + aI2*aJ2_2*r3*(one/9) + aI2*aJ2*r2*(one/9) - (four/945)*aI2*aJ2_6*r5*scale - (four/63)*aI2*aJ2_5*r4*scale - (four/9)*aI2*aJ2_4*r3*scale - (sixteen/9)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2_5*aJ2*r6*(one/945) + (two/189)*aI2_4*aJ2*r5 + aI2_3*aJ2*r4*(one/21) + aI2_2*aJ2*r3*(one/9) + aI2*aJ2*r2*(one/9) + (four/945)*aI2_6*aJ2*r5*scale + (four/63)*aI2_5*aJ2*r4*scale + (four/9)*aI2_4*aJ2*r3*scale + (sixteen/9)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp1 = (r2 - (4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(r2 + (4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp2 = (r2*(one/3) + aJ2*r3*(one/3) - ((four/3)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(r2*(one/3) + aI2*r3*(one/3) + ((four/3)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp3 = (aJ2_2*r4*(one/15) + aJ2*r3*(one/5) + r2*(one/5) - ((four/15)*aJ2_3*r3 + (eight/5)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(aI2_2*r4*(one/15) + aI2*r3*(one/5) + r2*(one/5) + ((four/15)*aI2_3*r3 + (eight/5)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp4 = (aJ2_3*r5*(one/105) + (two/35)*aJ2_2*r4 + aJ2*r3*(one/7) + r2*(one/7) - ((four/105)*aJ2_4*r4 + (eight/21)*aJ2_3*r3 + (twelve/7)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(aI2_3*r5*(one/105) + (two/35)*aI2_2*r4 + aI2*r3*(one/7) + r2*(one/7) + ((four/105)*aI2_4*r4 + (eight/21)*aI2_3*r3 + (twelve/7)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp5 = (aJ2_4*r6*(one/945) + (two/189)*aJ2_3*r5 + aJ2_2*r4*(one/21) + aJ2*r3*(one/9) + r2*(one/9) - ((four/945)*aJ2_5*r5 + (four/63)*aJ2_4*r4 + (four/9)*aJ2_3*r3 + (sixteen/9)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(aI2_4*r6*(one/945) + (two/189)*aI2_3*r5 + aI2_2*r4*(one/21) + aI2*r3*(one/9) + r2*(one/9) + ((four/945)*aI2_5*r5 + (four/63)*aI2_4*r4 + (four/9)*aI2_3*r3 + (sixteen/9)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
}
fexp = fexp/r;
fexp1 = fexp1/r3;
......
plugins/amoeba/platforms/cuda/src/kernels/hippoMutualField.cu
View file @ 9b523076
......@@ -8,10 +8,20 @@ real invR2 = invR*invR;
real invR3 = invR*invR2;
#if USE_EWALD
real ralpha = PME_ALPHA*r;
real bn0 = erfc(ralpha)*invR;
real exp2a = EXP(-(ralpha*ralpha));
#ifdef USE_DOUBLE_PRECISION
const real erfcAlphaR = erfc(ralpha);
#else
// This approximation for erfc is from Abramowitz and Stegun (1964) p. 299. They cite the following as
// the original source: C. Hastings, Jr., Approximations for Digital Computers (1955). It has a maximum
// error of 1.5e-7.
const real t = RECIP(1.0f+0.3275911f*ralpha);
const real erfcAlphaR = (0.254829592f+(-0.284496736f+(1.421413741f+(-1.453152027f+1.061405429f*t)*t)*t)*t)*t*exp2a;
#endif
real bn0 = erfcAlphaR*invR;
real alsq2 = 2*PME_ALPHA*PME_ALPHA;
real alsq2n = 1/(SQRT_PI*PME_ALPHA);
real exp2a = EXP(-(ralpha*ralpha));
alsq2n *= alsq2;
real bn1 = (bn0+alsq2n*exp2a)*invR2;
alsq2n *= alsq2;
......
plugins/amoeba/platforms/cuda/src/kernels/hippoNonbondedExceptions.cu
View file @ 9b523076
......@@ -6,7 +6,7 @@
extern "C" __global__ void computeNonbondedExceptions(
unsigned long long* __restrict__ forceBuffers, mixed* __restrict__ energyBuffer, unsigned long long* __restrict__ torqueBuffers,
const real4* __restrict__ posq, const real3* __restrict__ extDipole, const int2* __restrict__ exceptionAtoms, const real* __restrict__ mmScale,
const real* __restrict__ dmScale, const real* __restrict__ ddScale, const real* __restrict__ dispScale, const real* __restrict__ repScale,
const real* __restrict__ dmScale, const real* __restrict__ ddScale, const real* __restrict__ dispScale, const real* __restrict__ repScale, const real* __restrict__ ctScale,
const real* __restrict__ coreCharge, const real* __restrict__ valenceCharge, const real* __restrict__ alpha, const real* __restrict__ epsilon,
const real* __restrict__ damping, const real* __restrict__ c6, const real* __restrict__ pauliK, const real* __restrict__ pauliQ,
const real* __restrict__ pauliAlpha, const real3* __restrict__ dipole, const real3* __restrict__ inducedDipole, const real* __restrict__ qXX,
......@@ -70,6 +70,7 @@ extern "C" __global__ void computeNonbondedExceptions(
real dipoleDipoleScale = ddScale[index];
real repulsionScale = repScale[index];
real dispersionScale = dispScale[index];
real chargeTransferScale = ctScale[index];
real rInv = RSQRT(r2);
real r = r2*rInv;
real3 tempForce = make_real3(0);
......
plugins/amoeba/platforms/cuda/tests/TestCudaHippoNonbondedForce.cpp
View file @ 9b523076
......@@ -69,9 +69,9 @@ void buildWaterSystem(System& system, int numWaters, HippoNonbondedForce* hippo)
hippo->addParticle(0.19140, {0.0, 0.0, ds*-0.20097}, {qs*0.03881, 0.0, 0.0, 0.0, qs*0.02214, 0.0, 0.0, 0.0, qs*-0.06095}, 1.0,
10*4.7909, 0.0, 10*3.5582, c6s*4.5670, ps*2.0037, -0.8086, 10*4.6450,
0.001*0.341, HippoNonbondedForce::ZThenX, 3*mol, 3*mol+1, -1);
hippo->addException(3*mol, 3*mol+1, 0.0, 0.0, 0.2, 0.0, 0.0);
hippo->addException(3*mol, 3*mol+2, 0.0, 0.0, 0.2, 0.0, 0.0);
hippo->addException(3*mol+1, 3*mol+2, 0.0, 0.0, 1.0, 0.0, 0.0);
hippo->addException(3*mol, 3*mol+1, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0);
hippo->addException(3*mol, 3*mol+2, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0);
hippo->addException(3*mol+1, 3*mol+2, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0);
}
}
......@@ -95,8 +95,8 @@ void checkForceEnergyConsistency(Context& context) {
for (int i = 0; i < numParticles; ++i) {
Vec3 p = state.getPositions()[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
positions2[i] = p-f*step;
positions3[i] = p+f*step;
}
context.setPositions(positions2);
State state2 = context.getState(State::Energy);
......@@ -136,7 +136,7 @@ void testWaterDimer() {
for (int i = 0; i < system.getNumParticles(); i++)
ASSERT_EQUAL_VEC(-expectedForces[i], state.getForces()[i], forceTol);
// Compare the induced dipoles to reference values computed with Tinker.
// Compare the permanent dipoles to reference values computed with Tinker.
vector<Vec3> expectedLabDipoles = {
Vec3(-1.3999971343167967e-3, 0.0, 2.5377493339976591e-3),
......@@ -1528,7 +1528,7 @@ void testChangingParameters() {
hippo->setParticleParameters(0, -0.2, {0.0, 0.0, 0.005}, {0.001, 0.0, 0.0, 0.0, -0.001, 0.0, 0.0, 0.0, 0.0}, 6.0,
20, 4.184*1326.0, 10*40.0, 0.03, 2.0, -2.4233, 10*4.3097,
0.001*0.795, HippoNonbondedForce::Bisector, 1, 2, -1);
hippo->setExceptionParameters(3, 3, 4, 0.1, 0.1, 0.3, 0.05, 0.05);
hippo->setExceptionParameters(3, 3, 4, 0.1, 0.1, 0.3, 0.05, 0.05, 0.15);
hippo->updateParametersInContext(context);
State state2 = context.getState(State::Energy | State::Forces);
ASSERT(fabs(state1.getPotentialEnergy()-state2.getPotentialEnergy()) > 1.0)
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHippoNonbondedForce.cpp
View file @ 9b523076
......@@ -64,7 +64,7 @@ AmoebaReferenceHippoNonbondedForce::AmoebaReferenceHippoNonbondedForce(const Hip
for (int i = 0; i < force.getNumExceptions(); i++) {
Exception e;
force.getExceptionParameters(i, e.particle1, e.particle2, e.multipoleMultipoleScale, e.dipoleMultipoleScale,
e.dipoleDipoleScale, e.dispersionScale, e.repulsionScale);
e.dipoleDipoleScale, e.dispersionScale, e.repulsionScale, e.chargeTransferScale);
exceptions[make_pair(e.particle1, e.particle2)] = e;
exceptions[make_pair(e.particle2, e.particle1)] = e;
}
......@@ -522,21 +522,22 @@ void AmoebaReferenceHippoNonbondedForce::computeRepulsionDampingFactors(const Mu
double aJ2_3 = aJ2_2*aJ2;
double aJ2_4 = aJ2_2*aJ2_2;
double aJ2_5 = aJ2_3*aJ2_2;
double aJ2_6 = aJ2_3*aJ2_3;
double scale = 1/(aI2_2-aJ2_2);
double aI2aJ2expI = aI2*aJ2*expI;
double aI2aJ2expJ = aI2*aJ2*expJ;
pre = 8192*aI2_3*aJ2_3*(scale*scale*scale*scale);
double tmp = 4*aI2*aJ2*scale;
fexp = (arI2-tmp)*expJ + (arJ2+tmp)*expI;
fexp1 = (aI2*aJ2*r2 - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2*aJ2*r2 + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp2 = (aI2*aJ2*r2*(1.0/3) + aI2*aJ2_2*r3*(1.0/3) - (4.0/3)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2*aJ2*r2*(1.0/3) + aI2_2*aJ2*r3*(1.0/3) + (4.0/3)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp3 = (aI2*aJ2_3*r4*(1.0/15) + aI2*aJ2_2*r3*(1.0/5) + aI2*aJ2*r2*(1.0/5) - (4.0/15)*aI2*aJ2_4*r3*scale - (8.0/5)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*scale*aI2*aJ2)*expJ +
(aI2_3*aJ2*r4*(1.0/15) + aI2_2*aJ2*r3*(1.0/5) + aI2*aJ2*r2*(1.0/5) + (4.0/15)*aI2_4*aJ2*r3*scale + (8.0/5)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*scale*aI2*aJ2)*expI;
fexp4 = (aI2*aJ2_4*r5*(1.0/105) + (2.0/35)*aI2*aJ2_3*r4 + aI2*aJ2_2*r3*(1.0/7) + aI2*aJ2*r2*(1.0/7) - (4.0/105)*aI2*aJ2_5*r4*scale - (8.0/21)*aI2*aJ2_4*r3*scale - (12.0/7)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2_4*aJ2*r5*(1.0/105) + (2.0/35)*aI2_3*aJ2*r4 + aI2_2*aJ2*r3*(1.0/7) + aI2*aJ2*r2*(1.0/7) + (4.0/105)*aI2_5*aJ2*r4*scale + (8.0/21)*aI2_4*aJ2*r3*scale + (12.0/7)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp5 = (aI2*aJ2_5*r6*(1.0/945) + (2.0/189)*aI2*aJ2_4*r5 + aI2*aJ2_3*r4*(1.0/21) + aI2*aJ2_2*r3*(1.0/9) + aI2*aJ2*r2*(1.0/9) - (4.0/945)*aI2*aJ2_6*r5*scale - (4.0/63)*aI2*aJ2_5*r4*scale - (4.0/9)*aI2*aJ2_4*r3*scale - (16.0/9)*aI2*aJ2_3*r2*scale - 4*aI2*aJ2_2*r*scale - 4*aI2*aJ2*scale)*expJ +
(aI2_5*aJ2*r6*(1.0/945) + (2.0/189)*aI2_4*aJ2*r5 + aI2_3*aJ2*r4*(1.0/21) + aI2_2*aJ2*r3*(1.0/9) + aI2*aJ2*r2*(1.0/9) + (4.0/945)*aI2_6*aJ2*r5*scale + (4.0/63)*aI2_5*aJ2*r4*scale + (4.0/9)*aI2_4*aJ2*r3*scale + (16.0/9)*aI2_3*aJ2*r2*scale + 4*aI2_2*aJ2*r*scale + 4*aI2*aJ2*scale)*expI;
fexp1 = (r2 - (4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(r2 + (4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp2 = (r2*(1.0/3) + aJ2*r3*(1.0/3) - ((4.0/3)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(r2*(1.0/3) + aI2*r3*(1.0/3) + ((4.0/3)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp3 = (aJ2_2*r4*(1.0/15) + aJ2*r3*(1.0/5) + r2*(1.0/5) - ((4.0/15)*aJ2_3*r3 + (8.0/5)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(aI2_2*r4*(1.0/15) + aI2*r3*(1.0/5) + r2*(1.0/5) + ((4.0/15)*aI2_3*r3 + (8.0/5)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp4 = (aJ2_3*r5*(1.0/105) + (2.0/35)*aJ2_2*r4 + aJ2*r3*(1.0/7) + r2*(1.0/7) - ((4.0/105)*aJ2_4*r4 + (8.0/21)*aJ2_3*r3 + (12.0/7)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(aI2_3*r5*(1.0/105) + (2.0/35)*aI2_2*r4 + aI2*r3*(1.0/7) + r2*(1.0/7) + ((4.0/105)*aI2_4*r4 + (8.0/21)*aI2_3*r3 + (12.0/7)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
fexp5 = (aJ2_4*r6*(1.0/945) + (2.0/189)*aJ2_3*r5 + aJ2_2*r4*(1.0/21) + aJ2*r3*(1.0/9) + r2*(1.0/9) - ((4.0/945)*aJ2_5*r5 + (4.0/63)*aJ2_4*r4 + (4.0/9)*aJ2_3*r3 + (16.0/9)*aJ2_2*r2 + 4*aJ2*r + 4)*scale)*aI2aJ2expJ +
(aI2_4*r6*(1.0/945) + (2.0/189)*aI2_3*r5 + aI2_2*r4*(1.0/21) + aI2*r3*(1.0/9) + r2*(1.0/9) + ((4.0/945)*aI2_5*r5 + (4.0/63)*aI2_4*r4 + (4.0/9)*aI2_3*r3 + (16.0/9)*aI2_2*r2 + 4*aI2*r + 4)*scale)*aI2aJ2expI;
}
fexp = fexp/r;
fexp1 = fexp1/r3;
......@@ -1098,8 +1099,8 @@ double AmoebaReferenceHippoNonbondedForce::calculateChargeTransferPairIxn(const
}
auto exception = exceptions.find(make_pair(particleI.index, particleK.index));
if (exception != exceptions.end()) {
energy *= exception->second.multipoleMultipoleScale;
dEnergydR *= exception->second.multipoleMultipoleScale;
energy *= exception->second.chargeTransferScale;
dEnergydR *= exception->second.chargeTransferScale;
}
force[2] += dEnergydR;
return energy;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceHippoNonbondedForce.h
View file @ 9b523076
......@@ -194,7 +194,7 @@ protected:
class Exception {
public:
int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
};
unsigned int _numParticles;
......
plugins/amoeba/platforms/reference/tests/TestReferenceHippoNonbondedForce.cpp
View file @ 9b523076
......@@ -67,9 +67,9 @@ void buildWaterSystem(System& system, int numWaters, HippoNonbondedForce* hippo)
hippo->addParticle(0.19140, {0.0, 0.0, ds*-0.20097}, {qs*0.03881, 0.0, 0.0, 0.0, qs*0.02214, 0.0, 0.0, 0.0, qs*-0.06095}, 1.0,
10*4.7909, 0.0, 10*3.5582, c6s*4.5670, ps*2.0037, -0.8086, 10*4.6450,
0.001*0.341, HippoNonbondedForce::ZThenX, 3*mol, 3*mol+1, -1);
hippo->addException(3*mol, 3*mol+1, 0.0, 0.0, 0.2, 0.0, 0.0);
hippo->addException(3*mol, 3*mol+2, 0.0, 0.0, 0.2, 0.0, 0.0);
hippo->addException(3*mol+1, 3*mol+2, 0.0, 0.0, 1.0, 0.0, 0.0);
hippo->addException(3*mol, 3*mol+1, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0);
hippo->addException(3*mol, 3*mol+2, 0.0, 0.0, 0.2, 0.0, 0.0, 0.0);
hippo->addException(3*mol+1, 3*mol+2, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0);
}
}
......@@ -93,8 +93,8 @@ void checkForceEnergyConsistency(Context& context) {
for (int i = 0; i < numParticles; ++i) {
Vec3 p = state.getPositions()[i];
Vec3 f = state.getForces()[i];
positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
positions2[i] = p-f*step;
positions3[i] = p+f*step;
}
context.setPositions(positions2);
State state2 = context.getState(State::Energy);
......@@ -134,7 +134,7 @@ void testWaterDimer() {
for (int i = 0; i < system.getNumParticles(); i++)
ASSERT_EQUAL_VEC(-expectedForces[i], state.getForces()[i], 1e-5);
// Compare the induced dipoles to reference values computed with Tinker.
// Compare the permanent dipoles to reference values computed with Tinker.
vector<Vec3> expectedLabDipoles = {
Vec3(-1.3999971343167967e-3, 0.0, 2.5377493339976591e-3),
......@@ -1526,7 +1526,7 @@ void testChangingParameters() {
hippo->setParticleParameters(0, -0.2, {0.0, 0.0, 0.005}, {0.001, 0.0, 0.0, 0.0, -0.001, 0.0, 0.0, 0.0, 0.0}, 6.0,
20, 4.184*1326.0, 10*40.0, 0.03, 2.0, -2.4233, 10*4.3097,
0.001*0.795, HippoNonbondedForce::Bisector, 1, 2, -1);
hippo->setExceptionParameters(3, 3, 4, 0.1, 0.1, 0.3, 0.05, 0.05);
hippo->setExceptionParameters(3, 3, 4, 0.1, 0.1, 0.3, 0.05, 0.05, 0.15);
hippo->updateParametersInContext(context);
State state2 = context.getState(State::Energy | State::Forces);
ASSERT(fabs(state1.getPotentialEnergy()-state2.getPotentialEnergy()) > 1.0)
......
plugins/amoeba/serialization/src/HippoNonbondedForceProxy.cpp
View file @ 9b523076
......@@ -62,7 +62,7 @@ void HippoNonbondedForceProxy::serialize(const void* object, SerializationNode&
node.setIntProperty("dpmeGridX", nx);
node.setIntProperty("dpmeGridY", ny);
node.setIntProperty("dpmeGridZ", nz);
SerializationNode& coefficients = node.createChildNode("extrapolationCoefficients");
SerializationNode& coefficients = node.createChildNode("ExtrapolationCoefficients");
vector<double> coeff = force.getExtrapolationCoefficients();
for (int i = 0; i < coeff.size(); i++) {
stringstream key;
......@@ -100,9 +100,9 @@ void HippoNonbondedForceProxy::serialize(const void* object, SerializationNode&
}
SerializationNode& exceptions = node.createChildNode("Exceptions");
for (int i = 0; i < force.getNumExceptions(); i++) {
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
int p1, p2;
force.getExceptionParameters(i, p1, p2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale);
force.getExceptionParameters(i, p1, p2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
SerializationNode& exception = exceptions.createChildNode("Exception");
exception.setIntProperty("p1", p1);
exception.setIntProperty("p2", p2);
......@@ -111,6 +111,7 @@ void HippoNonbondedForceProxy::serialize(const void* object, SerializationNode&
exception.setDoubleProperty("ddScale", dipoleDipoleScale);
exception.setDoubleProperty("dispScale", dispersionScale);
exception.setDoubleProperty("repScale", repulsionScale);
exception.setDoubleProperty("ctScale", chargeTransferScale);
}
}
......@@ -127,7 +128,7 @@ void* HippoNonbondedForceProxy::deserialize(const SerializationNode& node) const
force->setEwaldErrorTolerance(node.getDoubleProperty("ewaldErrorTolerance"));
force->setPMEParameters(node.getDoubleProperty("pmeAlpha"), node.getIntProperty("pmeGridX"), node.getIntProperty("pmeGridY"), node.getIntProperty("pmeGridZ"));
force->setDPMEParameters(node.getDoubleProperty("dpmeAlpha"), node.getIntProperty("dpmeGridX"), node.getIntProperty("dpmeGridY"), node.getIntProperty("dpmeGridZ"));
const SerializationNode& coefficients = node.getChildNode("extrapolationCoefficients");
const SerializationNode& coefficients = node.getChildNode("ExtrapolationCoefficients");
vector<double> coeff;
for (int i = 0; ; i++) {
stringstream key;
......@@ -167,7 +168,8 @@ void* HippoNonbondedForceProxy::deserialize(const SerializationNode& node) const
const SerializationNode& exception = exceptions.getChildren()[i];
force->addException(exception.getIntProperty("p1"), exception.getIntProperty("p2"), exception.getDoubleProperty("mmScale"),
exception.getDoubleProperty("dmScale"), exception.getDoubleProperty("ddScale"),
exception.getDoubleProperty("dispScale"), exception.getDoubleProperty("repScale"));
exception.getDoubleProperty("dispScale"), exception.getDoubleProperty("repScale"),
exception.getDoubleProperty("ctScale"));
}
}
catch (...) {
......
plugins/amoeba/serialization/tests/TestSerializeHippoNonbondedForce.cpp
View file @ 9b523076
......@@ -67,7 +67,7 @@ void testSerialization() {
genrand_real2(sfmt), HippoNonbondedForce::Bisector, i+1, i+2, i+3);
}
for (int i = 0; i < 2; i++)
force1.addException(i, i+1, genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
force1.addException(i, i+1, genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
// Serialize and then deserialize it.
......@@ -131,10 +131,10 @@ void testSerialization() {
for (int i = 0; i < force1.getNumExceptions(); i++) {
int p11, p21;
int p12, p22;
double mm1, dm1, dd1, disp1, rep1;
double mm2, dm2, dd2, disp2, rep2;
force1.getExceptionParameters(i, p11, p21, mm1, dm1, dd1, disp1, rep1);
force2.getExceptionParameters(i, p12, p22, mm2, dm2, dd2, disp2, rep2);
double mm1, dm1, dd1, disp1, rep1, ct1;
double mm2, dm2, dd2, disp2, rep2, ct2;
force1.getExceptionParameters(i, p11, p21, mm1, dm1, dd1, disp1, rep1, ct1);
force2.getExceptionParameters(i, p12, p22, mm2, dm2, dd2, disp2, rep2, ct2);
ASSERT_EQUAL(p11, p12);
ASSERT_EQUAL(p21, p22);
ASSERT_EQUAL(mm1, mm2);
......@@ -142,6 +142,7 @@ void testSerialization() {
ASSERT_EQUAL(dd1, dd2);
ASSERT_EQUAL(disp1, disp2);
ASSERT_EQUAL(rep1, rep2);
ASSERT_EQUAL(ct1, ct2);
}
}
......
wrappers/python/simtk/openmm/app/forcefield.py
View file @ 9b523076
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2018 Stanford University and the Authors.
Portions copyright (c) 2012-2019 Stanford University and the Authors.
Authors: Peter Eastman, Mark Friedrichs
Contributors:
......@@ -1334,6 +1334,70 @@ def _findBondsForExclusions(data, sys):
bondIndices.append((child, atom.index))
return bondIndices
def _findExclusions(bondIndices, maxSeparation, numAtoms):
"""Identify pairs of atoms in the same molecule separated by no more than maxSeparation bonds."""
bondedTo = [set() for i in range(numAtoms)]
for i, j in bondIndices:
bondedTo[i].add(j)
bondedTo[j].add(i)
# Identify all neighbors of each atom with each separation.
bondedWithSeparation = [bondedTo]
for i in range(maxSeparation-1):
lastBonds = bondedWithSeparation[-1]
newBonds = deepcopy(lastBonds)
for atom in range(numAtoms):
for a1 in lastBonds[atom]:
for a2 in bondedTo[a1]:
newBonds[atom].add(a2)
bondedWithSeparation.append(newBonds)
# Build the list of pairs.
pairs = []
for atom in range(numAtoms):
for otherAtom in bondedWithSeparation[-1][atom]:
if otherAtom > atom:
# Determine the minimum number of bonds between them.
sep = maxSeparation
for i in reversed(range(maxSeparation-1)):
if otherAtom in bondedWithSeparation[i][atom]:
sep -= 1
else:
break
pairs.append((atom, otherAtom, sep))
return pairs
def _findGroups(bondedTo):
"""Given bonds that connect atoms, identify the connected groups."""
atomGroup = [None]*len(bondedTo)
numGroups = 0
for i in range(len(bondedTo)):
if atomGroup[i] is None:
# Start a new group.
atomStack = [i]
neighborStack = [0]
group = numGroups
numGroups += 1
# Recursively tag all the bonded atoms.
while len(atomStack) > 0:
atom = atomStack[-1]
atomGroup[atom] = group
while neighborStack[-1] < len(bondedTo[atom]) and atomGroup[bondedTo[atom][neighborStack[-1]]] is not None:
neighborStack[-1] += 1
if neighborStack[-1] < len(bondedTo[atom]):
atomStack.append(bondedTo[atom][neighborStack[-1]])
neighborStack.append(0)
else:
atomStack.pop()
neighborStack.pop()
return atomGroup
def _countResidueAtoms(elements):
"""Count the number of atoms of each element in a residue."""
counts = {}
......@@ -1749,6 +1813,7 @@ def _matchImproper(data, torsion, generator):
break
return match
# The following classes are generators that know how to create Force subclasses and add them to a System that is being
# created. Each generator class must define two methods: 1) a static method that takes an etree Element and a ForceField,
# and returns the corresponding generator object; 2) a createForce() method that constructs the Force object and adds it
......@@ -5413,6 +5478,112 @@ parsers["AmoebaUreyBradleyForce"] = AmoebaUreyBradleyGenerator.parseElement
#=============================================================================================
## @private
class HippoNonbondedGenerator(object):
"""A HippoNonbondedGenerator constructs a HippoNonbondedForce."""
def __init__(self, forcefield, extrapCoeff):
self.ff = forcefield
self.extrapCoeff = extrapCoeff
self.exceptions = {}
@staticmethod
def parseElement(element, ff):
extrapCoeff = [float(c) for c in element.attrib['extrapolationCoefficients'].split(',')]
generator = HippoNonbondedGenerator(ff, extrapCoeff)
ff.registerGenerator(generator)
scaleNames = ('mmScale', 'dmScale', 'ddScale', 'dispScale', 'repScale', 'ctScale')
paramNames = ('charge', 'coreCharge', 'alpha', 'epsilon', 'damping', 'c6', 'pauliK', 'pauliQ', 'pauliAlpha', 'polarizability', 'axisType', 'd0', 'd1', 'd2', 'q11', 'q12', 'q13', 'q21', 'q22', 'q23', 'q31', 'q32', 'q33')
for ex in element.findall('Exception'):
separation = int(ex.attrib['separation'])
ingroup = ex.attrib['ingroup'].lower() == 'true'
key = (separation, ingroup)
if key in generator.exceptions:
raise ValueError('HippoNonbondedForce: multiple exceptions with separation=%d ingroup=%s' % (separation, ingroup))
generator.exceptions[key] = [float(ex.attrib[s]) for s in scaleNames]
generator.params = ForceField._AtomTypeParameters(ff, 'HippoNonbondedForce', 'Atom', paramNames)
generator.params.parseDefinitions(element)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
methodMap = {NoCutoff:mm.HippoNonbondedForce.NoCutoff,
PME:mm.HippoNonbondedForce.PME}
if nonbondedMethod not in methodMap:
raise ValueError('Illegal nonbonded method for HippoNonbondedForce')
# Build data structures we'll need for building local coordinate frames.
bondIndices = _findBondsForExclusions(data, sys)
pairs = _findExclusions(bondIndices, 2, len(data.atoms))
bonded12 = [set() for i in range(len(data.atoms))]
bonded13 = [set() for i in range(len(data.atoms))]
for atom1, atom2, sep in pairs:
if sep == 1:
bonded12[atom1].add(data.atoms[atom2])
bonded12[atom2].add(data.atoms[atom1])
else:
bonded13[atom1].add(data.atoms[atom2])
bonded13[atom2].add(data.atoms[atom1])
# Create the force.
force = mm.HippoNonbondedForce()
for atom in data.atoms:
values = self.params.getAtomParameters(atom, data)
params = [float(v) for v in values[:10]]
axisType = int(values[10])
dipole = [float(v) for v in values[11:14]]
quadrupole = [float(v) for v in values[14:23]]
extra = self.params.getExtraParameters(atom, data)
zAtom = self._findAxisAtom('zAtomType', extra, bonded12[atom.index], None, data, [])
xAtom = self._findAxisAtom('xAtomType', extra, bonded12[atom.index], bonded13[atom.index], data, [zAtom])
yAtom = self._findAxisAtom('yAtomType', extra, bonded12[atom.index], bonded13[atom.index], data, [zAtom, xAtom])
force.addParticle(params[0], dipole, quadrupole, *params[1:], axisType=axisType, multipoleAtomZ=zAtom, multipoleAtomX=xAtom, multipoleAtomY=yAtom)
force.setNonbondedMethod(methodMap[nonbondedMethod])
force.setExtrapolationCoefficients(self.extrapCoeff)
force.setCutoffDistance(nonbondedCutoff)
if args['switchDistance'] is not None:
force.setSwitchingDistance(args['switchDistance'])
if 'ewaldErrorTolerance' in args:
force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
sys.addForce(force)
def _findAxisAtom(self, paramName, params, bonded12, bonded13, data, exclude):
if paramName not in params:
return -1
atomType = params[paramName]
for atom in bonded12:
if data.atomType[atom] == atomType and atom.index not in exclude:
return atom.index
if bonded13 is not None:
for atom in bonded13:
if data.atomType[atom] == atomType and atom.index not in exclude:
return atom.index
raise ValueError('No bonded atom of type %s' % atomType)
def postprocessSystem(self, sys, data, args):
# Identify polarization groups.
bondIndices = _findBondsForExclusions(data, sys)
groupBondTypes = [self.params.getExtraParameters(atom, data)['groupTypes'].split(',') for atom in data.atoms]
groupBonds = [[] for i in range(len(data.atoms))]
for i,j in bondIndices:
if data.atomType[data.atoms[i]] in groupBondTypes[j]:
groupBonds[i].append(j)
groupBonds[j].append(i)
polarizationGroup = _findGroups(groupBonds)
# Create the exclusions.
maxSeparation = max(e[0] for e in self.exceptions)
hippo = [f for f in sys.getForces() if isinstance(f, mm.HippoNonbondedForce)][0]
pairs = _findExclusions(bondIndices, maxSeparation, hippo.getNumParticles())
for atom1, atom2, sep in pairs:
params = self.exceptions[(sep, polarizationGroup[atom1] == polarizationGroup[atom2])]
hippo.addException(atom1, atom2, *params)
parsers["HippoNonbondedForce"] = HippoNonbondedGenerator.parseElement
## @private
class DrudeGenerator(object):
"""A DrudeGenerator constructs a DrudeForce."""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment