Commit 682da7b2 authored by Peter Eastman's avatar Peter Eastman
Browse files

Added a separate scale factor for charge transfer

parent b53c0f10
...@@ -287,12 +287,13 @@ public: ...@@ -287,12 +287,13 @@ public:
* @param dipoleDipoleScale the factor by which to scale the Coulomb interaction between induced dipoles * @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 dispersionScale the factor by which to scale the dispersion interaction
* @param repulsionScale the factor by which to scale the Pauli repulsion * @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. * @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. * If false, an exception is thrown.
* @return the index of the exception that was added * @return the index of the exception that was added
*/ */
int addException(int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale, 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. * Get the scale factors for an interaction that should be calculated differently from others.
* *
...@@ -304,9 +305,10 @@ public: ...@@ -304,9 +305,10 @@ public:
* @param dipoleDipoleScale the factor by which to scale the Coulomb interaction between induced dipoles * @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 dispersionScale the factor by which to scale the dispersion interaction
* @param repulsionScale the factor by which to scale the Pauli repulsion * @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, 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. * 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 * If all scale factors are set to 0, this will cause the interaction to be completely omitted from
...@@ -320,9 +322,10 @@ public: ...@@ -320,9 +322,10 @@ public:
* @param dipoleDipoleScale the factor by which to scale the Coulomb interaction between induced dipoles * @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 dispersionScale the factor by which to scale the dispersion interaction
* @param repulsionScale the factor by which to scale the Pauli repulsion * @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, 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 * 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) * which is acceptable. This value is used to select the grid dimensions and separation (alpha)
...@@ -424,14 +427,14 @@ public: ...@@ -424,14 +427,14 @@ public:
class HippoNonbondedForce::ExceptionInfo { class HippoNonbondedForce::ExceptionInfo {
public: public:
int particle1, particle2; int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale; double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
ExceptionInfo() { ExceptionInfo() {
particle1 = particle2 = -1; 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), particle1(particle1), particle2(particle2), multipoleMultipoleScale(multipoleMultipoleScale), dipoleMultipoleScale(dipoleMultipoleScale),
dipoleDipoleScale(dipoleDipoleScale), dispersionScale(dispersionScale), repulsionScale(repulsionScale) { dipoleDipoleScale(dipoleDipoleScale), dispersionScale(dispersionScale), repulsionScale(repulsionScale), chargeTransferScale(chargeTransferScale) {
} }
}; };
......
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * 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 * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -173,7 +173,7 @@ void HippoNonbondedForce::setParticleParameters(int index, double charge, const ...@@ -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, 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)); map<pair<int, int>, int>::iterator iter = exceptionMap.find(pair<int, int>(particle1, particle2));
int newIndex; int newIndex;
if (iter == exceptionMap.end()) if (iter == exceptionMap.end())
...@@ -187,19 +187,19 @@ int HippoNonbondedForce::addException(int particle1, int particle2, double multi ...@@ -187,19 +187,19 @@ int HippoNonbondedForce::addException(int particle1, int particle2, double multi
msg << particle2; msg << particle2;
throw OpenMMException(msg.str()); 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; newIndex = iter->second;
exceptionMap.erase(iter->first); exceptionMap.erase(iter->first);
} }
else { 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; newIndex = exceptions.size()-1;
} }
exceptionMap[pair<int, int>(particle1, particle2)] = newIndex; exceptionMap[pair<int, int>(particle1, particle2)] = newIndex;
return newIndex; return newIndex;
} }
void HippoNonbondedForce::getExceptionParameters(int index, int& particle1, int& particle2, double& multipoleMultipoleScale, double& dipoleMultipoleScale, double& dipoleDipoleScale, 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); ASSERT_VALID_INDEX(index, exceptions);
particle1 = exceptions[index].particle1; particle1 = exceptions[index].particle1;
particle2 = exceptions[index].particle2; particle2 = exceptions[index].particle2;
...@@ -208,10 +208,11 @@ void HippoNonbondedForce::getExceptionParameters(int index, int& particle1, int& ...@@ -208,10 +208,11 @@ void HippoNonbondedForce::getExceptionParameters(int index, int& particle1, int&
dipoleDipoleScale = exceptions[index].dipoleDipoleScale; dipoleDipoleScale = exceptions[index].dipoleDipoleScale;
dispersionScale = exceptions[index].dispersionScale; dispersionScale = exceptions[index].dispersionScale;
repulsionScale = exceptions[index].repulsionScale; repulsionScale = exceptions[index].repulsionScale;
chargeTransferScale = exceptions[index].chargeTransferScale;
} }
void HippoNonbondedForce::setExceptionParameters(int index, int particle1, int particle2, double multipoleMultipoleScale, double dipoleMultipoleScale, double dipoleDipoleScale, 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); ASSERT_VALID_INDEX(index, exceptions);
exceptions[index].particle1 = particle1; exceptions[index].particle1 = particle1;
exceptions[index].particle2 = particle2; exceptions[index].particle2 = particle2;
...@@ -220,6 +221,7 @@ void HippoNonbondedForce::setExceptionParameters(int index, int particle1, int p ...@@ -220,6 +221,7 @@ void HippoNonbondedForce::setExceptionParameters(int index, int particle1, int p
exceptions[index].dipoleDipoleScale = dipoleDipoleScale; exceptions[index].dipoleDipoleScale = dipoleDipoleScale;
exceptions[index].dispersionScale = dispersionScale; exceptions[index].dispersionScale = dispersionScale;
exceptions[index].repulsionScale = repulsionScale; exceptions[index].repulsionScale = repulsionScale;
exceptions[index].chargeTransferScale = chargeTransferScale;
} }
void HippoNonbondedForce::getInducedDipoles(Context& context, vector<Vec3>& dipoles) { void HippoNonbondedForce::getInducedDipoles(Context& context, vector<Vec3>& dipoles) {
......
...@@ -2624,20 +2624,20 @@ public: ...@@ -2624,20 +2624,20 @@ public:
} }
void getParticlesInGroup(int index, vector<int>& particles) { void getParticlesInGroup(int index, vector<int>& particles) {
int particle1, particle2; int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale; double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
force.getExceptionParameters(index, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale); force.getExceptionParameters(index, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
particles.resize(2); particles.resize(2);
particles[0] = particle1; particles[0] = particle1;
particles[1] = particle2; particles[1] = particle2;
} }
bool areGroupsIdentical(int group1, int group2) { bool areGroupsIdentical(int group1, int group2) {
int particle1, particle2; int particle1, particle2;
double multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1; double multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1, chargeTransferScale1;
double multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2; double multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2, chargeTransferScale2;
force.getExceptionParameters(group1, particle1, particle2, multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1); force.getExceptionParameters(group1, particle1, particle2, multipoleMultipoleScale1, dipoleMultipoleScale1, dipoleDipoleScale1, dispersionScale1, repulsionScale1, chargeTransferScale1);
force.getExceptionParameters(group2, particle1, particle2, multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2); force.getExceptionParameters(group2, particle1, particle2, multipoleMultipoleScale2, dipoleMultipoleScale2, dipoleDipoleScale2, dispersionScale2, repulsionScale2, chargeTransferScale2);
return (multipoleMultipoleScale1 == multipoleMultipoleScale2 && dipoleMultipoleScale1 == dipoleMultipoleScale2 && return (multipoleMultipoleScale1 == multipoleMultipoleScale2 && dipoleMultipoleScale1 == dipoleMultipoleScale2 &&
dipoleDipoleScale1 == dipoleDipoleScale2 && dispersionScale1 == dispersionScale2 && repulsionScale1 == repulsionScale2); dipoleDipoleScale1 == dipoleDipoleScale2 && dispersionScale1 == dispersionScale2 && repulsionScale1 == repulsionScale2 && chargeTransferScale1 == chargeTransferScale2);
} }
private: private:
const HippoNonbondedForce& force; const HippoNonbondedForce& force;
...@@ -2778,27 +2778,28 @@ void CudaCalcHippoNonbondedForceKernel::initialize(const System& system, const H ...@@ -2778,27 +2778,28 @@ void CudaCalcHippoNonbondedForceKernel::initialize(const System& system, const H
// Record exceptions and exclusions. // Record exceptions and exclusions.
vector<double> exceptionScaleVec[5]; vector<double> exceptionScaleVec[6];
vector<int2> exceptionAtomsVec; vector<int2> exceptionAtomsVec;
for (int i = 0; i < force.getNumExceptions(); i++) { for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2; int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale; double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale); force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
exclusions[particle1].push_back(particle2); exclusions[particle1].push_back(particle2);
exclusions[particle2].push_back(particle1); 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)); exceptionAtomsVec.push_back(make_int2(particle1, particle2));
exceptionScaleVec[0].push_back(multipoleMultipoleScale); exceptionScaleVec[0].push_back(multipoleMultipoleScale);
exceptionScaleVec[1].push_back(dipoleMultipoleScale); exceptionScaleVec[1].push_back(dipoleMultipoleScale);
exceptionScaleVec[2].push_back(dipoleDipoleScale); exceptionScaleVec[2].push_back(dipoleDipoleScale);
exceptionScaleVec[3].push_back(dispersionScale); exceptionScaleVec[3].push_back(dispersionScale);
exceptionScaleVec[4].push_back(repulsionScale); exceptionScaleVec[4].push_back(repulsionScale);
exceptionScaleVec[5].push_back(chargeTransferScale);
} }
} }
if (exceptionAtomsVec.size() > 0) { if (exceptionAtomsVec.size() > 0) {
exceptionAtoms.initialize<int2>(cu, exceptionAtomsVec.size(), "exceptionAtoms"); exceptionAtoms.initialize<int2>(cu, exceptionAtomsVec.size(), "exceptionAtoms");
exceptionAtoms.upload(exceptionAtomsVec); 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].initialize(cu, exceptionAtomsVec.size(), elementSize, "exceptionScales");
exceptionScales[i].upload(exceptionScaleVec[i], true); exceptionScales[i].upload(exceptionScaleVec[i], true);
} }
...@@ -3288,6 +3289,7 @@ double CudaCalcHippoNonbondedForceKernel::execute(ContextImpl& context, bool inc ...@@ -3288,6 +3289,7 @@ double CudaCalcHippoNonbondedForceKernel::execute(ContextImpl& context, bool inc
computeExceptionsArgs.push_back(&exceptionScales[2].getDevicePointer()); computeExceptionsArgs.push_back(&exceptionScales[2].getDevicePointer());
computeExceptionsArgs.push_back(&exceptionScales[3].getDevicePointer()); computeExceptionsArgs.push_back(&exceptionScales[3].getDevicePointer());
computeExceptionsArgs.push_back(&exceptionScales[4].getDevicePointer()); computeExceptionsArgs.push_back(&exceptionScales[4].getDevicePointer());
computeExceptionsArgs.push_back(&exceptionScales[5].getDevicePointer());
computeExceptionsArgs.push_back(&coreCharge.getDevicePointer()); computeExceptionsArgs.push_back(&coreCharge.getDevicePointer());
computeExceptionsArgs.push_back(&valenceCharge.getDevicePointer()); computeExceptionsArgs.push_back(&valenceCharge.getDevicePointer());
computeExceptionsArgs.push_back(&alpha.getDevicePointer()); computeExceptionsArgs.push_back(&alpha.getDevicePointer());
...@@ -3697,26 +3699,27 @@ void CudaCalcHippoNonbondedForceKernel::copyParametersToContext(ContextImpl& con ...@@ -3697,26 +3699,27 @@ void CudaCalcHippoNonbondedForceKernel::copyParametersToContext(ContextImpl& con
// Record the per-exception parameters. // Record the per-exception parameters.
vector<double> exceptionScaleVec[5]; vector<double> exceptionScaleVec[6];
vector<int2> exceptionAtomsVec; vector<int2> exceptionAtomsVec;
for (int i = 0; i < force.getNumExceptions(); i++) { for (int i = 0; i < force.getNumExceptions(); i++) {
int particle1, particle2; int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale; double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale); force.getExceptionParameters(i, particle1, particle2, multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale);
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)); exceptionAtomsVec.push_back(make_int2(particle1, particle2));
exceptionScaleVec[0].push_back(multipoleMultipoleScale); exceptionScaleVec[0].push_back(multipoleMultipoleScale);
exceptionScaleVec[1].push_back(dipoleMultipoleScale); exceptionScaleVec[1].push_back(dipoleMultipoleScale);
exceptionScaleVec[2].push_back(dipoleDipoleScale); exceptionScaleVec[2].push_back(dipoleDipoleScale);
exceptionScaleVec[3].push_back(dispersionScale); exceptionScaleVec[3].push_back(dispersionScale);
exceptionScaleVec[4].push_back(repulsionScale); exceptionScaleVec[4].push_back(repulsionScale);
exceptionScaleVec[5].push_back(chargeTransferScale);
} }
} }
if (exceptionAtomsVec.size() > 0) { if (exceptionAtomsVec.size() > 0) {
if (!exceptionAtoms.isInitialized() || exceptionAtoms.getSize() != exceptionAtomsVec.size()) if (!exceptionAtoms.isInitialized() || exceptionAtoms.getSize() != exceptionAtomsVec.size())
throw OpenMMException("updateParametersInContext: The number of exceptions has changed"); throw OpenMMException("updateParametersInContext: The number of exceptions has changed");
exceptionAtoms.upload(exceptionAtomsVec); exceptionAtoms.upload(exceptionAtomsVec);
for (int i = 0; i < 5; i++) for (int i = 0; i < 6; i++)
exceptionScales[i].upload(exceptionScaleVec[i], true); exceptionScales[i].upload(exceptionScaleVec[i], true);
} }
else if (exceptionAtoms.isInitialized()) else if (exceptionAtoms.isInitialized())
......
...@@ -733,7 +733,7 @@ private: ...@@ -733,7 +733,7 @@ private:
CudaArray dpmeBsplineModuliX, dpmeBsplineModuliY, dpmeBsplineModuliZ; CudaArray dpmeBsplineModuliX, dpmeBsplineModuliY, dpmeBsplineModuliZ;
CudaArray pmePhi, pmePhidp, pmeCphi; CudaArray pmePhi, pmePhidp, pmeCphi;
CudaArray lastPositions; CudaArray lastPositions;
CudaArray exceptionScales[5]; CudaArray exceptionScales[6];
CudaArray exceptionAtoms; CudaArray exceptionAtoms;
CudaSort* sort; CudaSort* sort;
cufftHandle fftForward, fftBackward, dfftForward, dfftBackward; cufftHandle fftForward, fftBackward, dfftForward, dfftBackward;
......
...@@ -423,8 +423,8 @@ real bn5 = (9*bn4+alsq2n*exp2a)*rInv2; ...@@ -423,8 +423,8 @@ real bn5 = (9*bn4+alsq2n*exp2a)*rInv2;
} }
#endif #endif
#ifdef COMPUTING_EXCEPTIONS #ifdef COMPUTING_EXCEPTIONS
ctForce *= multipoleMultipoleScale; ctForce *= chargeTransferScale;
ctEnergy *= multipoleMultipoleScale; ctEnergy *= chargeTransferScale;
#endif #endif
tempEnergy += includeInteraction ? ctEnergy : 0; tempEnergy += includeInteraction ? ctEnergy : 0;
tempForce.z += includeInteraction ? ctForce*r : 0; tempForce.z += includeInteraction ? ctForce*r : 0;
......
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
extern "C" __global__ void computeNonbondedExceptions( extern "C" __global__ void computeNonbondedExceptions(
unsigned long long* __restrict__ forceBuffers, mixed* __restrict__ energyBuffer, unsigned long long* __restrict__ torqueBuffers, 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 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__ 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__ 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, 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( ...@@ -70,6 +70,7 @@ extern "C" __global__ void computeNonbondedExceptions(
real dipoleDipoleScale = ddScale[index]; real dipoleDipoleScale = ddScale[index];
real repulsionScale = repScale[index]; real repulsionScale = repScale[index];
real dispersionScale = dispScale[index]; real dispersionScale = dispScale[index];
real chargeTransferScale = ctScale[index];
real rInv = RSQRT(r2); real rInv = RSQRT(r2);
real r = r2*rInv; real r = r2*rInv;
real3 tempForce = make_real3(0); real3 tempForce = make_real3(0);
......
...@@ -69,9 +69,9 @@ void buildWaterSystem(System& system, int numWaters, HippoNonbondedForce* hippo) ...@@ -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, 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, 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); 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+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); 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); hippo->addException(3*mol+1, 3*mol+2, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0);
} }
} }
...@@ -1528,7 +1528,7 @@ void testChangingParameters() { ...@@ -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, 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, 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); 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); hippo->updateParametersInContext(context);
State state2 = context.getState(State::Energy | State::Forces); State state2 = context.getState(State::Energy | State::Forces);
ASSERT(fabs(state1.getPotentialEnergy()-state2.getPotentialEnergy()) > 1.0) ASSERT(fabs(state1.getPotentialEnergy()-state2.getPotentialEnergy()) > 1.0)
......
...@@ -64,7 +64,7 @@ AmoebaReferenceHippoNonbondedForce::AmoebaReferenceHippoNonbondedForce(const Hip ...@@ -64,7 +64,7 @@ AmoebaReferenceHippoNonbondedForce::AmoebaReferenceHippoNonbondedForce(const Hip
for (int i = 0; i < force.getNumExceptions(); i++) { for (int i = 0; i < force.getNumExceptions(); i++) {
Exception e; Exception e;
force.getExceptionParameters(i, e.particle1, e.particle2, e.multipoleMultipoleScale, e.dipoleMultipoleScale, 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.particle1, e.particle2)] = e;
exceptions[make_pair(e.particle2, e.particle1)] = e; exceptions[make_pair(e.particle2, e.particle1)] = e;
} }
...@@ -1099,8 +1099,8 @@ double AmoebaReferenceHippoNonbondedForce::calculateChargeTransferPairIxn(const ...@@ -1099,8 +1099,8 @@ double AmoebaReferenceHippoNonbondedForce::calculateChargeTransferPairIxn(const
} }
auto exception = exceptions.find(make_pair(particleI.index, particleK.index)); auto exception = exceptions.find(make_pair(particleI.index, particleK.index));
if (exception != exceptions.end()) { if (exception != exceptions.end()) {
energy *= exception->second.multipoleMultipoleScale; energy *= exception->second.chargeTransferScale;
dEnergydR *= exception->second.multipoleMultipoleScale; dEnergydR *= exception->second.chargeTransferScale;
} }
force[2] += dEnergydR; force[2] += dEnergydR;
return energy; return energy;
......
...@@ -194,7 +194,7 @@ protected: ...@@ -194,7 +194,7 @@ protected:
class Exception { class Exception {
public: public:
int particle1, particle2; int particle1, particle2;
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale; double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
}; };
unsigned int _numParticles; unsigned int _numParticles;
......
...@@ -67,9 +67,9 @@ void buildWaterSystem(System& system, int numWaters, HippoNonbondedForce* hippo) ...@@ -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, 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, 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); 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+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); 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); hippo->addException(3*mol+1, 3*mol+2, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0);
} }
} }
...@@ -1526,7 +1526,7 @@ void testChangingParameters() { ...@@ -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, 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, 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); 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); hippo->updateParametersInContext(context);
State state2 = context.getState(State::Energy | State::Forces); State state2 = context.getState(State::Energy | State::Forces);
ASSERT(fabs(state1.getPotentialEnergy()-state2.getPotentialEnergy()) > 1.0) ASSERT(fabs(state1.getPotentialEnergy()-state2.getPotentialEnergy()) > 1.0)
......
...@@ -100,9 +100,9 @@ void HippoNonbondedForceProxy::serialize(const void* object, SerializationNode& ...@@ -100,9 +100,9 @@ void HippoNonbondedForceProxy::serialize(const void* object, SerializationNode&
} }
SerializationNode& exceptions = node.createChildNode("Exceptions"); SerializationNode& exceptions = node.createChildNode("Exceptions");
for (int i = 0; i < force.getNumExceptions(); i++) { for (int i = 0; i < force.getNumExceptions(); i++) {
double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale; double multipoleMultipoleScale, dipoleMultipoleScale, dipoleDipoleScale, dispersionScale, repulsionScale, chargeTransferScale;
int p1, p2; 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"); SerializationNode& exception = exceptions.createChildNode("Exception");
exception.setIntProperty("p1", p1); exception.setIntProperty("p1", p1);
exception.setIntProperty("p2", p2); exception.setIntProperty("p2", p2);
...@@ -111,6 +111,7 @@ void HippoNonbondedForceProxy::serialize(const void* object, SerializationNode& ...@@ -111,6 +111,7 @@ void HippoNonbondedForceProxy::serialize(const void* object, SerializationNode&
exception.setDoubleProperty("ddScale", dipoleDipoleScale); exception.setDoubleProperty("ddScale", dipoleDipoleScale);
exception.setDoubleProperty("dispScale", dispersionScale); exception.setDoubleProperty("dispScale", dispersionScale);
exception.setDoubleProperty("repScale", repulsionScale); exception.setDoubleProperty("repScale", repulsionScale);
exception.setDoubleProperty("ctScale", chargeTransferScale);
} }
} }
...@@ -167,7 +168,8 @@ void* HippoNonbondedForceProxy::deserialize(const SerializationNode& node) const ...@@ -167,7 +168,8 @@ void* HippoNonbondedForceProxy::deserialize(const SerializationNode& node) const
const SerializationNode& exception = exceptions.getChildren()[i]; const SerializationNode& exception = exceptions.getChildren()[i];
force->addException(exception.getIntProperty("p1"), exception.getIntProperty("p2"), exception.getDoubleProperty("mmScale"), force->addException(exception.getIntProperty("p1"), exception.getIntProperty("p2"), exception.getDoubleProperty("mmScale"),
exception.getDoubleProperty("dmScale"), exception.getDoubleProperty("ddScale"), exception.getDoubleProperty("dmScale"), exception.getDoubleProperty("ddScale"),
exception.getDoubleProperty("dispScale"), exception.getDoubleProperty("repScale")); exception.getDoubleProperty("dispScale"), exception.getDoubleProperty("repScale"),
exception.getDoubleProperty("ctScale"));
} }
} }
catch (...) { catch (...) {
......
...@@ -67,7 +67,7 @@ void testSerialization() { ...@@ -67,7 +67,7 @@ void testSerialization() {
genrand_real2(sfmt), HippoNonbondedForce::Bisector, i+1, i+2, i+3); genrand_real2(sfmt), HippoNonbondedForce::Bisector, i+1, i+2, i+3);
} }
for (int i = 0; i < 2; i++) 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. // Serialize and then deserialize it.
...@@ -131,10 +131,10 @@ void testSerialization() { ...@@ -131,10 +131,10 @@ void testSerialization() {
for (int i = 0; i < force1.getNumExceptions(); i++) { for (int i = 0; i < force1.getNumExceptions(); i++) {
int p11, p21; int p11, p21;
int p12, p22; int p12, p22;
double mm1, dm1, dd1, disp1, rep1; double mm1, dm1, dd1, disp1, rep1, ct1;
double mm2, dm2, dd2, disp2, rep2; double mm2, dm2, dd2, disp2, rep2, ct2;
force1.getExceptionParameters(i, p11, p21, mm1, dm1, dd1, disp1, rep1); force1.getExceptionParameters(i, p11, p21, mm1, dm1, dd1, disp1, rep1, ct1);
force2.getExceptionParameters(i, p12, p22, mm2, dm2, dd2, disp2, rep2); force2.getExceptionParameters(i, p12, p22, mm2, dm2, dd2, disp2, rep2, ct2);
ASSERT_EQUAL(p11, p12); ASSERT_EQUAL(p11, p12);
ASSERT_EQUAL(p21, p22); ASSERT_EQUAL(p21, p22);
ASSERT_EQUAL(mm1, mm2); ASSERT_EQUAL(mm1, mm2);
...@@ -142,6 +142,7 @@ void testSerialization() { ...@@ -142,6 +142,7 @@ void testSerialization() {
ASSERT_EQUAL(dd1, dd2); ASSERT_EQUAL(dd1, dd2);
ASSERT_EQUAL(disp1, disp2); ASSERT_EQUAL(disp1, disp2);
ASSERT_EQUAL(rep1, rep2); ASSERT_EQUAL(rep1, rep2);
ASSERT_EQUAL(ct1, ct2);
} }
} }
......
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