Skip to content

GitLab

Commit 51f82c02 authored by Peter Eastman's avatar Peter Eastman
Browse files

Merge branch 'genpt' of https://github.com/andysim/openmm into 

extrapolated
parents cab0faf8 2913b686
Hide whitespace changes
Inline Side-by-side
Showing with 937 additions and 17 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h
View file @ 51f82c02
......@@ -79,7 +79,13 @@ public:
/**
* Direct polarization
*/
Direct = 1
Direct = 1,
/**
* Optimized perturbation theory
*/
OPT = 2
};
enum MultipoleAxisTypes { ZThenX = 0, Bisector = 1, ZBisect = 2, ThreeFold = 3, ZOnly = 4, NoAxisType = 5, LastAxisTypeIndex = 6 };
......@@ -298,6 +304,21 @@ public:
*/
void setMutualInducedTargetEpsilon(double inputMutualInducedTargetEpsilon);
/**
* Set the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the perturbation
* theory algorithm for induced dipoles.
*
* @param optCoefficients a vector whose mth entry specifies the coefficient for mu_m
*
*/
void setOPTCoefficients(const std::vector<double> &OPTFullCoefficientsIn);
/**
* Get the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the perturbation
* theory algorithm for induced dipoles.
*/
const std::vector<double>& getOPTCoefficients() const;
/**
* 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)
......@@ -384,6 +405,8 @@ private:
int pmeBSplineOrder;
std::vector<int> pmeGridDimension;
int mutualInducedMaxIterations;
std::vector<double> OPTFullCoefficients;
double mutualInducedTargetEpsilon;
double scalingDistanceCutoff;
double electricConstant;
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp
View file @ 51f82c02
......@@ -61,6 +61,19 @@ void AmoebaMultipoleForce::setPolarizationType(AmoebaMultipoleForce::Polarizatio
polarizationType = type;
}
void AmoebaMultipoleForce::setOPTCoefficients(const std::vector<double> &OPTFullCoefficientsIn)
{
size_t maxPTOrder = OPTFullCoefficientsIn.size();
OPTFullCoefficients.resize(maxPTOrder);
std::copy(OPTFullCoefficientsIn.begin(), OPTFullCoefficientsIn.end(), OPTFullCoefficients.begin());
}
const std::vector<double> & AmoebaMultipoleForce::getOPTCoefficients() const
{
return OPTFullCoefficients;
}
double AmoebaMultipoleForce::getCutoffDistance() const {
return cutoffDistance;
}
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 51f82c02
......@@ -565,6 +565,8 @@ void ReferenceCalcAmoebaMultipoleForceKernel::initialize(const System& system, c
if (polarizationType == AmoebaMultipoleForce::Mutual) {
mutualInducedMaxIterations = force.getMutualInducedMaxIterations();
mutualInducedTargetEpsilon = force.getMutualInducedTargetEpsilon();
} else if (polarizationType == AmoebaMultipoleForce::OPT) {
OPTFullCoefficients = force.getOPTCoefficients();
}
// PME
......@@ -667,6 +669,9 @@ AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmo
amoebaReferenceMultipoleForce->setMaximumMutualInducedDipoleIterations(mutualInducedMaxIterations);
} else if (polarizationType == AmoebaMultipoleForce::Direct) {
amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Direct);
} else if (polarizationType == AmoebaMultipoleForce::OPT) {
amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::OPT);
amoebaReferenceMultipoleForce->setOPTCoefficients(OPTFullCoefficients);
} else {
throw OpenMMException("Polarization type not recognzied.");
}
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 51f82c02
......@@ -432,6 +432,7 @@ private:
int mutualInducedMaxIterations;
RealOpenMM mutualInducedTargetEpsilon;
std::vector<double> OPTFullCoefficients;
bool usePme;
RealOpenMM alphaEwald;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @ 51f82c02
......@@ -25,7 +25,6 @@
#include "AmoebaReferenceMultipoleForce.h"
#include "jama_svd.h"
#include <algorithm>
// In case we're using some primitive version of Visual Studio this will
// make sure that erf() and erfc() are defined.
#include "openmm/internal/MSVC_erfc.h"
......@@ -162,6 +161,19 @@ RealOpenMM AmoebaReferenceMultipoleForce::getMutualInducedDipoleTargetEpsilon()
return _mutualInducedDipoleTargetEpsilon;
}
void AmoebaReferenceMultipoleForce::setOPTCoefficients(const std::vector<RealOpenMM> &OPTFullCoefficients)
{
_maxPTOrder = OPTFullCoefficients.size(); // This accounts for the zero-based counting; actual highest order is 1 less
_OPTFullCoefficients.resize(_maxPTOrder);
_OPTPartCoefficients.resize(_maxPTOrder);
std::copy(OPTFullCoefficients.begin(), OPTFullCoefficients.end(), _OPTFullCoefficients.begin());
for(int i = 0; i < _maxPTOrder; ++i){
_OPTPartCoefficients[i] = 0.0;
for(int j = i; j < _maxPTOrder; ++j)
_OPTPartCoefficients[i] += _OPTFullCoefficients[j];
}
}
void AmoebaReferenceMultipoleForce::setMutualInducedDipoleTargetEpsilon(RealOpenMM mutualInducedDipoleTargetEpsilon)
{
_mutualInducedDipoleTargetEpsilon = mutualInducedDipoleTargetEpsilon;
......@@ -801,6 +813,9 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
RealVec deltaR = particleJ.position - particleI.position;
RealOpenMM r = SQRT(deltaR.dot(deltaR));
vector<RealOpenMM> rrI(2);
// If we're using the OPT algorithm, we need to compute the field gradient, so ask for one more rrI value.
if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT)
rrI.push_back(0.0);
getAndScaleInverseRs(particleI.dampingFactor, particleJ.dampingFactor,
particleI.thole, particleJ.thole, r, rrI);
......@@ -811,6 +826,50 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipolePairIxns(const Multipo
for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) {
calculateInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, rr3, rr5, deltaR,
*updateInducedDipoleFields[ii].inducedDipoles, updateInducedDipoleFields[ii].inducedDipoleField);
if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
// Compute and store the field gradient for later use.
RealOpenMM dx = deltaR[0];
RealOpenMM dy = deltaR[1];
RealOpenMM dz = deltaR[2];
OpenMM::RealVec &dipolesI = (*updateInducedDipoleFields[ii].inducedDipoles)[particleI.particleIndex];
RealOpenMM xDipole = dipolesI[0];
RealOpenMM yDipole = dipolesI[1];
RealOpenMM zDipole = dipolesI[2];
RealOpenMM muDotR = xDipole*dx + yDipole*dy + zDipole*dz;
RealOpenMM Exx = muDotR*dx*dx*rrI[2] - (2.0*xDipole*dx + muDotR)*rrI[1];
RealOpenMM Eyy = muDotR*dy*dy*rrI[2] - (2.0*yDipole*dy + muDotR)*rrI[1];
RealOpenMM Ezz = muDotR*dz*dz*rrI[2] - (2.0*zDipole*dz + muDotR)*rrI[1];
RealOpenMM Exy = muDotR*dx*dy*rrI[2] - (xDipole*dy + yDipole*dx)*rrI[1];
RealOpenMM Exz = muDotR*dx*dz*rrI[2] - (xDipole*dz + zDipole*dx)*rrI[1];
RealOpenMM Eyz = muDotR*dy*dz*rrI[2] - (yDipole*dz + zDipole*dy)*rrI[1];
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][0] -= Exx;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][1] -= Eyy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][2] -= Ezz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][3] -= Exy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][4] -= Exz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][5] -= Eyz;
OpenMM::RealVec &dipolesJ = (*updateInducedDipoleFields[ii].inducedDipoles)[particleJ.particleIndex];
xDipole = dipolesJ[0];
yDipole = dipolesJ[1];
zDipole = dipolesJ[2];
muDotR = xDipole*dx + yDipole*dy + zDipole*dz;
Exx = muDotR*dx*dx*rrI[2] - (2.0*xDipole*dx + muDotR)*rrI[1];
Eyy = muDotR*dy*dy*rrI[2] - (2.0*yDipole*dy + muDotR)*rrI[1];
Ezz = muDotR*dz*dz*rrI[2] - (2.0*zDipole*dz + muDotR)*rrI[1];
Exy = muDotR*dx*dy*rrI[2] - (xDipole*dy + yDipole*dx)*rrI[1];
Exz = muDotR*dx*dz*rrI[2] - (xDipole*dz + zDipole*dx)*rrI[1];
Eyz = muDotR*dy*dz*rrI[2] - (yDipole*dz + zDipole*dy)*rrI[1];
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][0] += Exx;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][1] += Eyy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][2] += Ezz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][3] += Exy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][4] += Exz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][5] += Eyz;
}
}
}
......@@ -899,6 +958,86 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesBySOR(const vector<Mult
setMutualInducedDipoleIterations(iteration);
}
void AmoebaReferenceMultipoleForce::convergeInduceDipolesByOPT(const vector<MultipoleParticleData>& particleData, vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleField) {
if (_OPTFullCoefficients.empty()){
std::stringstream message;
message << "An OPT calcultion was requested, but setOPTCoefficients() was not called.";
throw OpenMMException(message.str());
}
_ptDipoleD.clear();
_ptDipoleP.clear();
UpdateInducedDipoleFieldStruct& fieldD = updateInducedDipoleField[0];
UpdateInducedDipoleFieldStruct& fieldP = updateInducedDipoleField[1];
// Start by storing the direct dipoles as PT0
vector<RealVec> thisDipoleD;
vector<RealVec> thisDipoleP;
for(int atom = 0; atom < _numParticles; ++atom){
thisDipoleD.push_back((*fieldD.inducedDipoles)[atom]);
thisDipoleP.push_back((*fieldP.inducedDipoles)[atom]);
}
_ptDipoleD.push_back(thisDipoleD);
_ptDipoleP.push_back(thisDipoleP);
// Make sure there is some storage available for the field derivatives
std::vector<RealOpenMM> zeros(6, 0.0);
fieldD.inducedDipoleFieldGradient.resize(_numParticles);
fieldP.inducedDipoleFieldGradient.resize(_numParticles);
// Recursively apply alpha.Tau to the µ_(n) components to generate µ_(n+1), and store the result
for(int order = 1; order < _maxPTOrder; ++order){
std::fill(fieldD.inducedDipoleFieldGradient.begin(), fieldD.inducedDipoleFieldGradient.end(), zeros);
std::fill(fieldP.inducedDipoleFieldGradient.begin(), fieldP.inducedDipoleFieldGradient.end(), zeros);
calculateInducedDipoleFields(particleData, updateInducedDipoleField);
vector<RealVec> thisDipoleD;
vector<RealVec> thisDipoleP;
for(int atom = 0; atom < _numParticles; ++atom){
(*fieldD.inducedDipoles)[atom] = fieldD.inducedDipoleField[atom] * particleData[atom].polarity;
(*fieldP.inducedDipoles)[atom] = fieldP.inducedDipoleField[atom] * particleData[atom].polarity;
thisDipoleD.push_back((*fieldD.inducedDipoles)[atom]);
thisDipoleP.push_back((*fieldP.inducedDipoles)[atom]);
}
_ptDipoleD.push_back(thisDipoleD);
_ptDipoleP.push_back(thisDipoleP);
vector<RealOpenMM> fieldGradD(6*_numParticles, 0.0);
vector<RealOpenMM> fieldGradP(6*_numParticles, 0.0);
for(int atom = 0; atom < _numParticles; ++atom){
for(int component = 0; component < 6; ++component){
fieldGradD[6*atom + component] = fieldD.inducedDipoleFieldGradient[atom][component];
fieldGradP[6*atom + component] = fieldP.inducedDipoleFieldGradient[atom][component];
}
}
_ptDipoleFieldGradientD.push_back(fieldGradD);
_ptDipoleFieldGradientP.push_back(fieldGradP);
}
// Take a linear combination of the µ_(n) components to form the total dipole
RealVec zeroVec(0.0, 0.0, 0.0);
std::fill(_inducedDipole.begin(), _inducedDipole.end(), zeroVec);
std::fill(_inducedDipolePolar.begin(), _inducedDipolePolar.end(), zeroVec);
for(int order = 0; order < _maxPTOrder; ++order){
for(int atom = 0; atom < _numParticles; ++atom){
_inducedDipole[atom] += _ptDipoleD[order][atom] * _OPTPartCoefficients[order];
_inducedDipolePolar[atom] += _ptDipoleP[order][atom] * _OPTPartCoefficients[order];
}
}
// Copy the combined dipoles over to compute the field
for(int order = 0; order < _maxPTOrder; ++order){
for(int atom = 0; atom < _numParticles; ++atom){
(*fieldD.inducedDipoles)[atom] = _inducedDipole[atom];
(*fieldP.inducedDipoles)[atom] = _inducedDipolePolar[atom];
}
}
calculateInducedDipoleFields(particleData, updateInducedDipoleField);
setMutualInducedDipoleConverged(true);
}
void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<MultipoleParticleData>& particleData, vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleField) {
int numFields = updateInducedDipoleField.size();
vector<vector<vector<RealVec> > > prevDipoles(numFields);
......@@ -964,6 +1103,7 @@ void AmoebaReferenceMultipoleForce::convergeInduceDipolesByDIIS(const vector<Mul
}
}
}
}
void AmoebaReferenceMultipoleForce::computeDIISCoefficients(const vector<vector<RealVec> >& prevErrors, vector<RealOpenMM>& coefficients) const {
......@@ -1037,15 +1177,19 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const vector<Multipo
// UpdateInducedDipoleFieldStruct contains induced dipole, fixed multipole fields and fields
// due to other induced dipoles at each site
if (getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual){
convergeInduceDipolesByDIIS(particleData, updateInducedDipoleField);
}if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
convergeInduceDipolesByOPT(particleData, updateInducedDipoleField);
}
convergeInduceDipolesByDIIS(particleData, updateInducedDipoleField);
}
RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostaticPairIxn(const MultipoleParticleData& particleI,
const MultipoleParticleData& particleK,
const vector<RealOpenMM>& scalingFactors,
vector<RealVec>& forces,
vector<RealVec>& torque) const
vector<RealVec>& torque) const
{
unsigned int iIndex = particleI.particleIndex;
unsigned int kIndex = particleK.particleIndex;
......@@ -1675,6 +1819,9 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Mu
vector<RealVec>& torques,
vector<RealVec>& forces)
{
const int deriv1[] = {1, 4, 7, 8, 10, 15, 17, 13, 14, 19};
const int deriv2[] = {2, 7, 5, 9, 13, 11, 18, 15, 19, 16};
const int deriv3[] = {3, 8, 9, 6, 14, 16, 12, 19, 17, 18};
RealOpenMM energy = 0.0;
vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX);
......@@ -1683,7 +1830,6 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Mu
}
// main loop over particle pairs
for (unsigned int ii = 0; ii < particleData.size(); ii++) {
for (unsigned int jj = ii+1; jj < particleData.size(); jj++) {
......@@ -1700,6 +1846,36 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateElectrostatic(const vector<Mu
}
}
}
if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
RealOpenMM prefac = (_electric/_dielectric);
for (int i = 0; i < _numParticles; i++) {
// Compute the µ(m) T µ(n) force contributions here
for(int l = 0; l < _maxPTOrder-1; ++l) {
for(int m = 0; m < _maxPTOrder-1-l; ++m) {
RealOpenMM p = _OPTPartCoefficients[l+m+1];
if(std::fabs(p) < 1e-6) continue;
forces[i][0] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+0]
+ _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+3]
+ _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+4]);
forces[i][1] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+3]
+ _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+1]
+ _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+5]);
forces[i][2] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+4]
+ _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+5]
+ _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+2]);
forces[i][0] += 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+0]
+ _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+3]
+ _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+4]);
forces[i][1] += 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+3]
+ _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+1]
+ _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+5]);
forces[i][2] += 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+4]
+ _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+5]
+ _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+2]);
}
}
}
}
return energy;
}
......@@ -5821,13 +5997,11 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::computeReciprocalSpaceInducedDipole
f[1] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j2];
f[2] += (inducedDipole[k]+inducedDipolePolar[k])*_phi[20*i+j3];
if (polarizationType == AmoebaReferenceMultipoleForce::Mutual)
{
if (polarizationType == AmoebaReferenceMultipoleForce::Mutual) {
f[0] += (inducedDipole[k]*_phip[10*i+j1] + inducedDipolePolar[k]*_phid[10*i+j1]);
f[1] += (inducedDipole[k]*_phip[10*i+j2] + inducedDipolePolar[k]*_phid[10*i+j2]);
f[2] += (inducedDipole[k]*_phip[10*i+j3] + inducedDipolePolar[k]*_phid[10*i+j3]);
}
}
for (int k = 0; k < 10; k++) {
......@@ -5880,6 +6054,7 @@ void AmoebaReferencePmeMultipoleForce::recordInducedDipoleField(vector<RealVec>&
fieldPolar[i][1] -= _phip[10*i+1]*fracToCart[1][0] + _phip[10*i+2]*fracToCart[1][1] + _phip[10*i+3]*fracToCart[1][2];
fieldPolar[i][2] -= _phip[10*i+1]*fracToCart[2][0] + _phip[10*i+2]*fracToCart[2][1] + _phip[10*i+3]*fracToCart[2][2];
}
}
void AmoebaReferencePmeMultipoleForce::calculateReciprocalSpaceInducedDipoleField(vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields)
......@@ -5916,7 +6091,69 @@ void AmoebaReferencePmeMultipoleForce::calculateInducedDipoleFields(const vector
calculateReciprocalSpaceInducedDipoleField(updateInducedDipoleFields);
// self ixn
if(getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
// While we have the reciprocal space (fractional coordinate) field gradient available, add it to the real space
// terms computed above, after transforming to Cartesian coordinates. This allows real and reciprocal space
// dipole response force contributions to be computed together.
RealVec fracToCart[3];
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++)
fracToCart[i][j] = _pmeGridDimensions[j]*_recipBoxVectors[i][j];
for (int i = 0; i < _numParticles; i++) {
RealOpenMM EmatD[3][3] = {
{ _phid[10*i+4], _phid[10*i+7], _phid[10*i+8] },
{ _phid[10*i+7], _phid[10*i+5], _phid[10*i+9] },
{ _phid[10*i+8], _phid[10*i+9], _phid[10*i+6] }
};
RealOpenMM Exx = 0.0, Eyy = 0.0, Ezz = 0.0, Exy = 0.0, Exz = 0.0, Eyz = 0.0;
for(int k = 0; k < 3; ++k){
for(int l = 0; l < 3; ++l){
Exx += fracToCart[0][k] * EmatD[k][l] * fracToCart[0][l];
Eyy += fracToCart[1][k] * EmatD[k][l] * fracToCart[1][l];
Ezz += fracToCart[2][k] * EmatD[k][l] * fracToCart[2][l];
Exy += fracToCart[0][k] * EmatD[k][l] * fracToCart[1][l];
Exz += fracToCart[0][k] * EmatD[k][l] * fracToCart[2][l];
Eyz += fracToCart[1][k] * EmatD[k][l] * fracToCart[2][l];
}
}
updateInducedDipoleFields[0].inducedDipoleFieldGradient[i][0] -= Exx;
updateInducedDipoleFields[0].inducedDipoleFieldGradient[i][1] -= Eyy;
updateInducedDipoleFields[0].inducedDipoleFieldGradient[i][2] -= Ezz;
updateInducedDipoleFields[0].inducedDipoleFieldGradient[i][3] -= Exy;
updateInducedDipoleFields[0].inducedDipoleFieldGradient[i][4] -= Exz;
updateInducedDipoleFields[0].inducedDipoleFieldGradient[i][5] -= Eyz;
RealOpenMM EmatP[3][3] = {
{ _phip[10*i+4], _phip[10*i+7], _phip[10*i+8] },
{ _phip[10*i+7], _phip[10*i+5], _phip[10*i+9] },
{ _phip[10*i+8], _phip[10*i+9], _phip[10*i+6] }
};
Exx = 0.0; Eyy = 0.0; Ezz = 0.0; Exy = 0.0; Exz = 0.0; Eyz = 0.0;
for(int k = 0; k < 3; ++k){
for(int l = 0; l < 3; ++l){
Exx += fracToCart[0][k] * EmatP[k][l] * fracToCart[0][l];
Eyy += fracToCart[1][k] * EmatP[k][l] * fracToCart[1][l];
Ezz += fracToCart[2][k] * EmatP[k][l] * fracToCart[2][l];
Exy += fracToCart[0][k] * EmatP[k][l] * fracToCart[1][l];
Exz += fracToCart[0][k] * EmatP[k][l] * fracToCart[2][l];
Eyz += fracToCart[1][k] * EmatP[k][l] * fracToCart[2][l];
}
}
updateInducedDipoleFields[1].inducedDipoleFieldGradient[i][0] -= Exx;
updateInducedDipoleFields[1].inducedDipoleFieldGradient[i][1] -= Eyy;
updateInducedDipoleFields[1].inducedDipoleFieldGradient[i][2] -= Ezz;
updateInducedDipoleFields[1].inducedDipoleFieldGradient[i][3] -= Exy;
updateInducedDipoleFields[1].inducedDipoleFieldGradient[i][4] -= Exz;
updateInducedDipoleFields[1].inducedDipoleFieldGradient[i][5] -= Eyz;
}
}
// self ixn
RealOpenMM term = (4.0/3.0)*(_alphaEwald*_alphaEwald*_alphaEwald)/SQRT_PI;
for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) {
......@@ -5980,10 +6217,15 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
alsq2n *= alsq2;
RealOpenMM bn2 = (3.0*bn1+alsq2n*exp2a)/r2;
alsq2n *= alsq2;
RealOpenMM bn3 = (5.0*bn2+alsq2n*exp2a)/r2;
// compute the error function scaled and unscaled terms
RealOpenMM scale3 = 1.0;
RealOpenMM scale5 = 1.0;
RealOpenMM scale7 = 1.0;
RealOpenMM damp = particleI.dampingFactor*particleJ.dampingFactor;
if (damp != 0.0) {
......@@ -5996,24 +6238,73 @@ void AmoebaReferencePmeMultipoleForce::calculateDirectInducedDipolePairIxns(cons
RealOpenMM expdamp = expf(damp);
scale3 = 1.0 - expdamp;
scale5 = 1.0 - expdamp*(1.0-damp);
scale7 = 1.0 - (1.0 - damp + (0.6*damp*damp))*expdamp;
}
}
RealOpenMM dsc3 = uscale*scale3;
RealOpenMM dsc5 = uscale*scale5;
RealOpenMM dsc7 = uscale*scale7;
RealOpenMM r3 = (r*r2);
RealOpenMM r5 = (r3*r2);
RealOpenMM r7 = (r5*r2);
RealOpenMM rr3 = (1.0-dsc3)/r3;
RealOpenMM rr5 = 3.0*(1.0-dsc5)/r5;
RealOpenMM rr7 = 15.0*(1.0-dsc7)/r7;
RealOpenMM preFactor1 = rr3 - bn1;
RealOpenMM preFactor2 = bn2 - rr5;
RealOpenMM preFactor3 = bn3 - rr7;
for (unsigned int ii = 0; ii < updateInducedDipoleFields.size(); ii++) {
calculateDirectInducedDipolePairIxn(particleI.particleIndex, particleJ.particleIndex, preFactor1, preFactor2, deltaR,
*updateInducedDipoleFields[ii].inducedDipoles,
updateInducedDipoleFields[ii].inducedDipoleField);
}
if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
// Compute and store the field gradient for later use.
RealOpenMM dx = deltaR[0];
RealOpenMM dy = deltaR[1];
RealOpenMM dz = deltaR[2];
OpenMM::RealVec &dipolesI = (*updateInducedDipoleFields[ii].inducedDipoles)[particleI.particleIndex];
RealOpenMM xDipole = dipolesI[0];
RealOpenMM yDipole = dipolesI[1];
RealOpenMM zDipole = dipolesI[2];
RealOpenMM muDotR = xDipole*dx + yDipole*dy + zDipole*dz;
RealOpenMM Exx = muDotR*dx*dx*preFactor3 - (2.0*xDipole*dx + muDotR)*preFactor2;
RealOpenMM Eyy = muDotR*dy*dy*preFactor3 - (2.0*yDipole*dy + muDotR)*preFactor2;
RealOpenMM Ezz = muDotR*dz*dz*preFactor3 - (2.0*zDipole*dz + muDotR)*preFactor2;
RealOpenMM Exy = muDotR*dx*dy*preFactor3 - (xDipole*dy + yDipole*dx)*preFactor2;
RealOpenMM Exz = muDotR*dx*dz*preFactor3 - (xDipole*dz + zDipole*dx)*preFactor2;
RealOpenMM Eyz = muDotR*dy*dz*preFactor3 - (yDipole*dz + zDipole*dy)*preFactor2;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][0] -= Exx;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][1] -= Eyy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][2] -= Ezz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][3] -= Exy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][4] -= Exz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleJ.particleIndex][5] -= Eyz;
OpenMM::RealVec &dipolesJ = (*updateInducedDipoleFields[ii].inducedDipoles)[particleJ.particleIndex];
xDipole = dipolesJ[0];
yDipole = dipolesJ[1];
zDipole = dipolesJ[2];
muDotR = xDipole*dx + yDipole*dy + zDipole*dz;
Exx = muDotR*dx*dx*preFactor3 - (2.0*xDipole*dx + muDotR)*preFactor2;
Eyy = muDotR*dy*dy*preFactor3 - (2.0*yDipole*dy + muDotR)*preFactor2;
Ezz = muDotR*dz*dz*preFactor3 - (2.0*zDipole*dz + muDotR)*preFactor2;
Exy = muDotR*dx*dy*preFactor3 - (xDipole*dy + yDipole*dx)*preFactor2;
Exz = muDotR*dx*dz*preFactor3 - (xDipole*dz + zDipole*dx)*preFactor2;
Eyz = muDotR*dy*dz*preFactor3 - (yDipole*dz + zDipole*dy)*preFactor2;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][0] += Exx;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][1] += Eyy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][2] += Ezz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][3] += Exy;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][4] += Exz;
updateInducedDipoleFields[ii].inducedDipoleFieldGradient[particleI.particleIndex][5] += Eyz;
}
}
}
RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeSelfEnergy(const vector<MultipoleParticleData>& particleData) const
......@@ -6429,7 +6720,7 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPair
RealOpenMM iEJY = qiUinpJ[0]*Vjip[1] + qiUindJ[0]*Vjid[1] - qiUinpJ[1]*Vjip[0] - qiUindJ[1]*Vjid[0];
// Add in the induced-induced terms, if needed.
if(getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual){
if(getPolarizationType() == AmoebaReferenceMultipoleForce::Mutual) {
// Uind-Uind terms (m=0)
RealOpenMM eCoef = -fourThirds*rInvVec[3]*(3.0*(uScale*thole_d0 + bVec[3]) + alphaRVec[3]*X);
RealOpenMM dCoef = rInvVec[4]*(6.0*(uScale*dthole_d0 + bVec[3]) + 4.0*alphaRVec[5]*X);
......@@ -6478,7 +6769,6 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculatePmeDirectElectrostaticPair
RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector<MultipoleParticleData>& particleData,
vector<RealVec>& torques, vector<RealVec>& forces)
{
RealOpenMM energy = 0.0;
vector<RealOpenMM> scaleFactors(LAST_SCALE_TYPE_INDEX);
for (unsigned int kk = 0; kk < scaleFactors.size(); kk++) {
......@@ -6504,10 +6794,43 @@ RealOpenMM AmoebaReferencePmeMultipoleForce::calculateElectrostatic(const vector
}
}
// The polarization energy
calculatePmeSelfTorque(particleData, torques);
energy += computeReciprocalSpaceInducedDipoleForceAndEnergy(getPolarizationType(), particleData, forces, torques);
energy += computeReciprocalSpaceFixedMultipoleForceAndEnergy(particleData, forces, torques);
energy += calculatePmeSelfEnergy(particleData);
// Now that both the direct and reciprocal space contributions have been added, we can compute the dipole
// response contributions to the forces, if we're using the OPT polarization algorithm.
if (getPolarizationType() == AmoebaReferenceMultipoleForce::OPT){
RealOpenMM prefac = (_electric/_dielectric);
for (int i = 0; i < _numParticles; i++) {
// Compute the µ(m) T µ(n) force contributions here
for(int l = 0; l < _maxPTOrder-1; ++l) {
for(int m = 0; m < _maxPTOrder-1-l; ++m) {
RealOpenMM p = _OPTPartCoefficients[l+m+1];
if(std::fabs(p) < 1e-6) continue;
forces[i][0] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+0]
+ _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+3]
+ _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+4]);
forces[i][1] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+3]
+ _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+1]
+ _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+5]);
forces[i][2] += 0.5*p*prefac*(_ptDipoleD[l][i][0]*_ptDipoleFieldGradientP[m][6*i+4]
+ _ptDipoleD[l][i][1]*_ptDipoleFieldGradientP[m][6*i+5]
+ _ptDipoleD[l][i][2]*_ptDipoleFieldGradientP[m][6*i+2]);
forces[i][0] += 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+0]
+ _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+3]
+ _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+4]);
forces[i][1] += 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+3]
+ _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+1]
+ _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+5]);
forces[i][2] += 0.5*p*prefac*(_ptDipoleP[l][i][0]*_ptDipoleFieldGradientD[m][6*i+4]
+ _ptDipoleP[l][i][1]*_ptDipoleFieldGradientD[m][6*i+5]
+ _ptDipoleP[l][i][2]*_ptDipoleFieldGradientD[m][6*i+2]);
}
}
}
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.h
View file @ 51f82c02
......@@ -345,11 +345,16 @@ public:
*/
Mutual = 0,
/**
/**
* Direct polarization
*/
Direct = 1
};
Direct = 1,
/**
* Optimized perturbation theory
*/
OPT = 2
};
/**
* Constructor
......@@ -422,6 +427,15 @@ public:
*/
RealOpenMM getMutualInducedDipoleEpsilon() const;
/**
* Set the coefficients for the µ_0, µ_1, µ_2, µ_n terms in the pertubation
* theory algorithm for induced dipoles
*
* @param optCoefficients a vector whose mth entry specifies the coefficient for µ_m
*
*/
void setOPTCoefficients(const std::vector<RealOpenMM> &OPTFullCoefficients);
/**
* Set the target epsilon for converging mutual induced dipoles.
*
......@@ -628,6 +642,7 @@ protected:
std::vector<OpenMM::RealVec>* fixedMultipoleField;
std::vector<OpenMM::RealVec>* inducedDipoles;
std::vector<OpenMM::RealVec> inducedDipoleField;
std::vector<std::vector<RealOpenMM> > inducedDipoleFieldGradient;
};
unsigned int _numParticles;
......@@ -651,10 +666,19 @@ protected:
std::vector<RealVec> _fixedMultipoleFieldPolar;
std::vector<RealVec> _inducedDipole;
std::vector<RealVec> _inducedDipolePolar;
std::vector< std::vector<RealVec> > _ptDipoleP;
std::vector< std::vector<RealVec> > _ptDipoleD;
std::vector<std::vector<RealOpenMM> > _ptDipoleFieldGradientP;
std::vector<std::vector<RealOpenMM> > _ptDipoleFieldGradientD;
int _mutualInducedDipoleConverged;
int _mutualInducedDipoleIterations;
int _maximumMutualInducedDipoleIterations;
int _maxPTOrder;
std::vector<RealOpenMM> _OPTFullCoefficients;
std::vector<RealOpenMM> _OPTPartCoefficients;
RealOpenMM _mutualInducedDipoleEpsilon;
RealOpenMM _mutualInducedDipoleTargetEpsilon;
RealOpenMM _polarSOR;
......@@ -904,7 +928,7 @@ protected:
/**
* Calculate fields due induced dipoles at each site.
*
*
* @param particleI positions and parameters (charge, labFrame dipoles, quadrupoles, ...) for particle I
* @param particleJ positions and parameters (charge, labFrame dipoles, quadrupoles, ...) for particle J
* @param updateInducedDipoleFields vector of UpdateInducedDipoleFieldStruct containing input induced dipoles and output fields
......@@ -920,6 +944,14 @@ protected:
*/
virtual void calculateInducedDipoleFields(const std::vector<MultipoleParticleData>& particleData,
std::vector<UpdateInducedDipoleFieldStruct>& updateInducedDipoleFields);
/**
* Calculated induced dipoles using Optimized Perturbation Theory.
*
* @param particleData vector of particle positions and parameters (charge, labFrame dipoles, quadrupoles, ...)
* @param updateInducedDipoleFields vector of UpdateInducedDipoleFieldStruct containing input induced dipoles and output fields
*/
void convergeInduceDipolesByOPT(const std::vector<MultipoleParticleData>& particleData,
std::vector<UpdateInducedDipoleFieldStruct>& calculateInducedDipoleField);
/**
* Converge induced dipoles.
*
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaPTPolarization.cpp 0 → 100644
View file @ 51f82c02
/* -------------------------------------------------------------------------- *
* OpenMMAmoeba *
* -------------------------------------------------------------------------- *
* This is part of the OpenMM molecular simulation toolkit originating from *
* 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) 2008-2015 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the Reference implementation of the PT polarization algorithms in ReferenceAmoebaMultipoleForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "openmm/AmoebaMultipoleForce.h"
#include "openmm/LangevinIntegrator.h"
#include "openmm/Vec3.h"
#include "openmm/Units.h"
#include <iostream>
#include <iomanip>
#include <vector>
#include <stdlib.h>
#include <stdio.h>
#define ASSERT_EQUAL_TOL_MOD(expected, found, tol, testname) {double _scale_ = std::abs(expected) > 1.0 ? std::abs(expected) : 1.0; if (!(std::abs((expected)-(found))/_scale_ <= (tol))) {std::stringstream details; details << testname << " Expected "<<(expected)<<", found "<<(found); throwException(__FILE__, __LINE__, details.str());}};
#define ASSERT_EQUAL_VEC_MOD(expected, found, tol,testname) {ASSERT_EQUAL_TOL_MOD((expected)[0], (found)[0], (tol),(testname)); ASSERT_EQUAL_TOL_MOD((expected)[1], (found)[1], (tol),(testname)); ASSERT_EQUAL_TOL_MOD((expected)[2], (found)[2], (tol),(testname));};
using namespace OpenMM;
using namespace std;
extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
const double TOL = 1e-4;
// print the energy and forces out, in AKMA units, to allow comparison with TINKER
static void printEnergyAndForces(double energy, vector<Vec3> &forces){
size_t natoms = forces.size();
double sf = 1.0;
std::cout << "Energy (SI):" << std::setw(20) << std::setprecision(10) << energy << std::endl;
std::cout << "Forces (SI):" << std::endl;
for(int i = 0; i < natoms; ++i){
std::cout << i+1 << "\t" << std::setw(20) << std::setprecision(10) << forces[i][0]*sf <<
std::setw(20) << std::setprecision(10) << forces[i][1]*sf <<
std::setw(20) << std::setprecision(10) << forces[i][2]*sf << std::endl;
}
sf = -OpenMM::KcalPerKJ/10.0;
std::cout << "Energy (AKMA):" << std::setw(20) << std::setprecision(10) << energy*OpenMM::KcalPerKJ << std::endl;
std::cout << "Forces (AKMA):" << std::endl;
for(int i = 0; i < natoms; ++i){
std::cout << i+1 << "\t" << std::setw(20) << std::setprecision(10) << forces[i][0]*sf <<
std::setw(20) << std::setprecision(10) << forces[i][1]*sf <<
std::setw(20) << std::setprecision(10) << forces[i][2]*sf << std::endl;
}
}
// compare forces and energies
static void compareForcesEnergy(std::string& testName, double expectedEnergy, double energy,
const std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance) {
for (unsigned int ii = 0; ii < forces.size(); ii++) {
ASSERT_EQUAL_VEC_MOD(expectedForces[ii], forces[ii], tolerance, testName);
}
ASSERT_EQUAL_TOL_MOD(expectedEnergy, energy, tolerance, testName);
}
// compare relative differences in force norms and energies
static void compareForceNormsEnergy(std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
const std::vector<Vec3>& forces, double tolerance) {
for (unsigned int ii = 0; ii < forces.size(); ii++) {
double expectedNorm = sqrt(expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2]);
double norm = sqrt(forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2]);
double absDiff = fabs(norm - expectedNorm);
double relDiff = 2.0*absDiff/(fabs(norm) + fabs(expectedNorm) + 1.0e-08);
if (relDiff > tolerance && absDiff > 0.001) {
std::stringstream details;
details << testName << "Relative difference in norms " << relDiff << " larger than allowed tolerance at particle=" << ii;
details << ": norms=" << norm << " expected norm=" << expectedNorm;
throwException(__FILE__, __LINE__, details.str());
}
}
double energyAbsDiff = fabs(expectedEnergy - energy);
double energyRelDiff = 2.0*energyAbsDiff/(fabs(expectedEnergy) + fabs(energy) + 1.0e-08);
if (energyRelDiff > tolerance) {
std::stringstream details;
details << testName << "Relative difference in energies " << energyRelDiff << " larger than allowed tolerance.";
details << "Energies=" << energy << " expected energy=" << expectedEnergy;
throwException(__FILE__, __LINE__, details.str());
}
}
vector<Vec3> setupWaterDimer(System& system, AmoebaMultipoleForce* amoebaMultipoleForce, bool use_pol_groups) {
const int NATOMS = 6;
const char* atom_types[NATOMS] = {"O", "H1", "H2", "O", "H1", "H2"};
const double coords[NATOMS][3] = {
{ 2.000000, 2.000000, 2.000000},
{ 2.500000, 2.000000, 3.000000},
{ 1.500000, 2.000000, 3.000000},
{ 0.000000, 0.000000, 0.000000},
{ 0.500000, 0.000000, 1.000000},
{ -0.500000, 0.000000, 1.000000}
};
std::map < std::string, double > tholemap;
std::map < std::string, double > polarmap;
std::map < std::string, double > chargemap;
std::map < std::string, std::vector<double> > dipolemap;
std::map < std::string, std::vector<double> > quadrupolemap;
std::map < std::string, AmoebaMultipoleForce::MultipoleAxisTypes > axesmap;
std::map < std::string, std::vector<int> > anchormap;
std::map < std::string, double > massmap;
std::map < std::string, std::vector<int> > polgrpmap;
std::map < std::string, std::vector<int> > cov12map;
std::map < std::string, std::vector<int> > cov13map;
axesmap["O"] = AmoebaMultipoleForce::Bisector;
axesmap["H1"] = AmoebaMultipoleForce::ZThenX;
axesmap["H2"] = AmoebaMultipoleForce::ZThenX;
chargemap["O"] = -0.51966;
chargemap["H1"] = 0.25983;
chargemap["H2"] = 0.25983;
int oanc[3] = {1, 2, 0};
int h1anc[3] = {-1, 1, 0};
int h2anc[3] = {-2, -1, 0};
std::vector<int> oancv(&oanc[0], &oanc[3]);
std::vector<int> h1ancv(&h1anc[0], &h1anc[3]);
std::vector<int> h2ancv(&h2anc[0], &h2anc[3]);
anchormap["O"] = oancv;
anchormap["H1"] = h1ancv;
anchormap["H2"] = h2ancv;
double od[3] = {0.0, 0.0, 0.00755612136146};
double hd[3] = {-0.00204209484795, 0.0, -0.00307875299958};
std::vector<double> odv(&od[0], &od[3]);
std::vector<double> hdv(&hd[0], &hd[3]);
dipolemap["O"] = odv;
dipolemap["H1"] = hdv;
dipolemap["H2"] = hdv;
double oq[9] = {0.000354030721139, 0.0, 0.0,
0.0, -0.000390257077096, 0.0,
0.0, 0.0, 3.62263559571e-05};
double hq[9] = {-3.42848248983e-05, 0.0, -1.89485963908e-06,
0.0, -0.000100240875193, 0.0,
-1.89485963908e-06, 0.0, 0.000134525700091};
std::vector<double> oqv(&oq[0], &oq[9]);
std::vector<double> hqv(&hq[0], &hq[9]);
quadrupolemap["O"] = oqv;
quadrupolemap["H1"] = hqv;
quadrupolemap["H2"] = hqv;
polarmap["O"] = 0.3069876538;
polarmap["H1"] = 0.2813500172;
polarmap["H2"] = 0.2813500172;
polarmap["O"] = 0.000837;
polarmap["H1"] = 0.000496;
polarmap["H2"] = 0.000496;
tholemap["O"] = 0.3900;
tholemap["H1"] = 0.3900;
tholemap["H2"] = 0.3900;
massmap["O"] = 15.999;
massmap["H1"] = 1.0080000;
massmap["H2"] = 1.0080000;
int opg[3] = {0,1,2};
int h1pg[3] = {-1,0,1};
int h2pg[3] = {-2,-1,0};
std::vector<int> opgv(&opg[0], &opg[3]);
std::vector<int> h1pgv(&h1pg[0], &h1pg[3]);
std::vector<int> h2pgv(&h2pg[0], &h2pg[3]);
if(!use_pol_groups){
opgv.clear();
h1pgv.clear();
h2pgv.clear();
}
polgrpmap["O"] = opgv;
polgrpmap["H1"] = h1pgv;
polgrpmap["H2"] = h2pgv;
int cov12o[2] = {1,2};
int cov12h1[1] = {-1};
int cov12h2[1] = {-2};
std::vector<int> cov12ov(&cov12o[0], &cov12o[2]);
std::vector<int> cov12h1v(&cov12h1[0], &cov12h1[1]);
std::vector<int> cov12h2v(&cov12h2[0], &cov12h2[1]);
cov12map["O"] = cov12ov;
cov12map["H1"] = cov12h1v;
cov12map["H2"] = cov12h2v;
int cov13h1[1] = {1};
int cov13h2[1] = {-1};
std::vector<int> cov13h1v(&cov13h1[0], &cov13h1[1]);
std::vector<int> cov13h2v(&cov13h2[0], &cov13h2[1]);
cov13map["O"] = std::vector<int>();
cov13map["H1"] = cov13h1v;
cov13map["H2"] = cov13h2v;
std::vector<Vec3> positions(NATOMS);
for(int atom = 0; atom < NATOMS; ++atom){
const char* element = atom_types[atom];
double damp = polarmap[element];
double alpha = pow(damp, 1.0/6.0);
int atomz = atom + anchormap[element][0];
int atomx = atom + anchormap[element][1];
int atomy = anchormap[element][2]==0 ? -1 : atom + anchormap[element][2];
amoebaMultipoleForce->addMultipole(chargemap[element], dipolemap[element], quadrupolemap[element],
axesmap[element], atomz, atomx, atomy, tholemap[element], alpha, damp);
system.addParticle(massmap[element]);
double offset =0.0;
positions[atom] = Vec3(coords[atom][0]+offset, coords[atom][1]+offset, coords[atom][2]+offset)*OpenMM::NmPerAngstrom;
// Polarization groups
std::vector<int> tmppol;
std::vector<int>& polgrps = polgrpmap[element];
for(int i=0; i < polgrps.size(); ++i)
tmppol.push_back(polgrps[i]+atom);
if(!tmppol.empty())
amoebaMultipoleForce->setCovalentMap(atom, AmoebaMultipoleForce::PolarizationCovalent11, tmppol);
// 1-2 covalent groups
std::vector<int> tmp12;
std::vector<int>& cov12s = cov12map[element];
for(int i=0; i < cov12s.size(); ++i)
tmp12.push_back(cov12s[i]+atom);
if(!tmp12.empty())
amoebaMultipoleForce->setCovalentMap(atom, AmoebaMultipoleForce::Covalent12, tmp12);
// 1-3 covalent groups
std::vector<int> tmp13;
std::vector<int>& cov13s = cov13map[element];
for(int i=0; i < cov13s.size(); ++i)
tmp13.push_back(cov13s[i]+atom);
if(!tmp13.empty())
amoebaMultipoleForce->setCovalentMap(atom, AmoebaMultipoleForce::Covalent13, tmp13);
}
system.addForce(amoebaMultipoleForce);
return positions;
}
static void check_finite_differences(vector<Vec3> analytic_forces, Context &context, vector<Vec3> positions)
{
// Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
double norm = 0.0;
for (int i = 0; i < (int) analytic_forces.size(); ++i)
norm += analytic_forces[i].dot(analytic_forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
vector<Vec3> positions2(analytic_forces.size()), positions3(analytic_forces.size());
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = analytic_forces[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);
}
context.setPositions(positions2);
State state2 = context.getState(State::Energy);
context.setPositions(positions3);
State state3 = context.getState(State::Energy);
ASSERT_EQUAL_TOL(norm, (state2.getPotentialEnergy()-state3.getPotentialEnergy())/stepSize, 1e-4);
}
static void testWaterDimerExPTPolarizationTriclinicPME() {
std::string testName = "testWaterDimerExPTPolarizationTriclinicPME";
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
vector<Vec3> coords = setupWaterDimer(system, amoebaMultipoleForce, true);
system.setDefaultPeriodicBoxVectors(Vec3(2.0, 0.0, 0.0),
Vec3(0.2, 2.0, 0.0),
Vec3(0.1, 0.5, 2.0));
amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::PME);
amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
std::vector<double> coefs;
coefs.push_back(0.0); // The mu_0 coefficient
coefs.push_back(-0.3); // The mu_1 coefficient
coefs.push_back(0.0); // The mu_2 coefficient
coefs.push_back(1.3); // The mu_3 coefficient
amoebaMultipoleForce->setOPTCoefficients(coefs);
amoebaMultipoleForce->setCutoffDistance(9.0*OpenMM::NmPerAngstrom);
amoebaMultipoleForce->setAEwald(4);
amoebaMultipoleForce->setEwaldErrorTolerance(1.0e-06);
std::vector<int> pmeGridDimension(3);
pmeGridDimension[0] = pmeGridDimension[1] = pmeGridDimension[2] = 64;
amoebaMultipoleForce->setPmeGridDimensions(pmeGridDimension);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
context.setPositions(coords);
OpenMM::State state = context.getState(State::Forces | State::Energy);
std::vector<Vec3> forces = state.getForces();
double energy = state.getPotentialEnergy();
// printEnergyAndForces(energy, forces);
double expectedEnergy = -1.945797427;
std::vector<Vec3> expectedForces(forces.size());
expectedForces[0] = Vec3( -131.1099603, -187.2725558, 36.94657685);
expectedForces[1] = Vec3( 38.6397841, 2.410997985, 8.008437937);
expectedForces[2] = Vec3( 38.69034185, 117.5018257, 32.43097836);
expectedForces[3] = Vec3( -117.3212339, -102.3366145, -30.50621066);
expectedForces[4] = Vec3( 124.8343077, 169.7729804, -24.10742414);
expectedForces[5] = Vec3( 46.26244074, -0.07194110719, -22.77727325);
double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
check_finite_differences(forces, context, coords);
}
static void testWaterDimerExPTPolarizationTriclinicPMENoPolGroups() {
std::string testName = "testWaterDimerExPTPolarizationTriclinicPMENoPolGroups";
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
vector<Vec3> coords = setupWaterDimer(system, amoebaMultipoleForce, false);
system.setDefaultPeriodicBoxVectors(Vec3(2.0, 0.0, 0.0),
Vec3(0.2, 2.0, 0.0),
Vec3(0.1, 0.5, 2.0));
amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::PME);
amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
std::vector<double> coefs;
coefs.push_back(0.0); // The mu_0 coefficient
coefs.push_back(-0.3); // The mu_1 coefficient
coefs.push_back(0.0); // The mu_2 coefficient
coefs.push_back(1.3); // The mu_3 coefficient
amoebaMultipoleForce->setOPTCoefficients(coefs);
amoebaMultipoleForce->setCutoffDistance(9.0*OpenMM::NmPerAngstrom);
amoebaMultipoleForce->setAEwald(4);
amoebaMultipoleForce->setEwaldErrorTolerance(1.0e-06);
std::vector<int> pmeGridDimension(3);
pmeGridDimension[0] = pmeGridDimension[1] = pmeGridDimension[2] = 64;
amoebaMultipoleForce->setPmeGridDimensions(pmeGridDimension);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
context.setPositions(coords);
OpenMM::State state = context.getState(State::Forces | State::Energy);
std::vector<Vec3> forces = state.getForces();
double energy = state.getPotentialEnergy();
// printEnergyAndForces(energy, forces);
double expectedEnergy = -1.840068409;
std::vector<Vec3> expectedForces(forces.size());
expectedForces[0] = Vec3( -69.85154559, -104.2092334, 3.586495334);
expectedForces[1] = Vec3( 19.50350452, -14.5844519, 9.400418341);
expectedForces[2] = Vec3( 16.75641493, 75.15006506, 19.14553199);
expectedForces[3] = Vec3( -67.24268213, -47.39994175, -18.81277222);
expectedForces[4] = Vec3( 75.15808251, 110.6109313, 4.355432435);
expectedForces[5] = Vec3( 25.67255306, -19.56378113, -17.68217953);
double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
check_finite_differences(forces, context, coords);
}
static void testWaterDimerExPTPolarizationNoCutoff() {
std::string testName = "testWaterDimerExPTPolarizationNoCutoff";
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
vector<Vec3> coords = setupWaterDimer(system, amoebaMultipoleForce, true);
amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
std::vector<double> coefs;
coefs.push_back(0.0); // The mu_0 coefficient
coefs.push_back(-0.3); // The mu_1 coefficient
coefs.push_back(0.0); // The mu_2 coefficient
coefs.push_back(1.3); // The mu_3 coefficient
amoebaMultipoleForce->setOPTCoefficients(coefs);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
context.setPositions(coords);
OpenMM::State state = context.getState(State::Forces | State::Energy);
std::vector<Vec3> forces = state.getForces();
double energy = state.getPotentialEnergy();
// printEnergyAndForces(energy, forces);
double expectedEnergy = -1.399194432;
std::vector<Vec3> expectedForces(forces.size());
expectedForces[0] = Vec3( -130.7294487, -186.3287444, 41.40628056);
expectedForces[1] = Vec3( 38.90143386, 2.140957908, 5.564712102);
expectedForces[2] = Vec3( 38.32881448, 117.0462626, 29.90093041);
expectedForces[3] = Vec3( -117.1147396, -101.6981494, -25.55733439);
expectedForces[4] = Vec3( 124.7421318, 169.1571359, -26.38724373);
expectedForces[5] = Vec3( 45.87180816, -0.3174626947, -24.92734495);
double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
check_finite_differences(forces, context, coords);
}
static void testWaterDimerExPTPolarizationNoCutoffNoPolGroups() {
std::string testName = "testWaterDimerExPTPolarizationNoCutoffNoPolGroups";
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
vector<Vec3> coords = setupWaterDimer(system, amoebaMultipoleForce, false);
amoebaMultipoleForce->setNonbondedMethod(AmoebaMultipoleForce::NoCutoff);
amoebaMultipoleForce->setPolarizationType(AmoebaMultipoleForce::OPT);
std::vector<double> coefs;
coefs.push_back(0.0); // The mu_0 coefficient
coefs.push_back(-0.3); // The mu_1 coefficient
coefs.push_back(0.0); // The mu_2 coefficient
coefs.push_back(1.3); // The mu_3 coefficient
amoebaMultipoleForce->setOPTCoefficients(coefs);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
context.setPositions(coords);
OpenMM::State state = context.getState(State::Forces | State::Energy);
std::vector<Vec3> forces = state.getForces();
double energy = state.getPotentialEnergy();
// printEnergyAndForces(energy, forces);
double expectedEnergy = -1.56926564;
std::vector<Vec3> expectedForces(forces.size());
expectedForces[0] = Vec3( -69.623843, -103.7006124, 6.162774255);
expectedForces[1] = Vec3( 19.54326912, -14.69441322, 8.014369439);
expectedForces[2] = Vec3( 16.65441143, 74.88100242, 17.70364405);
expectedForces[3] = Vec3( -67.10049929, -47.08900953, -16.01092086);
expectedForces[4] = Vec3( 74.98800293, 110.2649458, 3.020145768);
expectedForces[5] = Vec3( 25.53865881, -19.66191302, -18.89001266);
double tolerance = 1.0e-04;
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
check_finite_differences(forces, context, coords);
}
int main(int numberOfArguments, char* argv[]) {
try {
std::cout << "TestReferenceAmoebaPTPolarization running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
/*
* Water dimer energy / force tests under various conditions.
*/
// PME, triclinic
testWaterDimerExPTPolarizationTriclinicPME();
// PME, triclinic, no polarization groups
testWaterDimerExPTPolarizationTriclinicPMENoPolGroups();
// No cutoff
testWaterDimerExPTPolarizationNoCutoff();
// No cutoff, no polarization groups
testWaterDimerExPTPolarizationNoCutoffNoPolGroups();
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
}
wrappers/python/src/swig_doxygen/swigInputConfig.py
View file @ 51f82c02
......@@ -247,6 +247,7 @@ UNITS = {
("AmoebaMultipoleForce", "getPmeBSplineOrder") : ( None,()),
("AmoebaMultipoleForce", "getMutualInducedMaxIterations") : ( None, ()),
("AmoebaMultipoleForce", "getMutualInducedTargetEpsilon") : ( None, ()),
("AmoebaMultipoleForce", "getOPTCoefficients") : ( None, ()),
("AmoebaMultipoleForce", "getEwaldErrorTolerance") : ( None, ()),
("AmoebaMultipoleForce", "getPmeGridDimensions") : ( None,()),
......
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