Skip to content

GitLab

Commit 71de4b1b authored by peastman's avatar peastman
Browse files

Merge pull request #1337 from peastman/extrapolated 

AMOEBA extrapolated polarization
parents 83408400 49cbc791
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 2384 additions and 236 deletions
docs-source/usersguide/application.rst
View file @ 71de4b1b
......@@ -922,29 +922,33 @@ If :code:`vdwCutoff` is not specified, then the value of
Specifying the Polarization Method
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
OpenMM allows the setting of several other parameters particular to the AMOEBA
force field. The :code:`mutualInducedTargetEpsilon` option allows you to
specify the accuracy to which the induced dipoles are calculated at each time
step; the default value is 0.01. The :code:`polarization` setting
determines whether the calculation of the induced dipoles is continued until the
dipoles are self-consistent to within the tolerance specified by
:code:`mutualInducedTargetEpsilon` or whether a quick estimate of the induced
dipoles is used instead. The first option corresponds to the
:code:`polarization='mutual'` setting and is the default; the quick estimate
option is given by :code:`polarization='direct'` and in this case,
:code:`mutualInducedTargetEpsilon` is ignored, if provided. Simulations using
:code:`polarization='direct'` will be significantly faster than those with
:code:`polarization='mutual'`\ , but less accurate. Examples using the two
options are given below:
::
system = forcefield.createSystem(nonbondedMethod=PME,
nonbondedCutoff=1*nanometer,ewaldErrorTolerance=0.00001,
vdwCutoff=1.2*nanometer, mutualInducedTargetEpsilon=0.01)
system = forcefield.createSystem(nonbondedMethod=PME,
nonbondedCutoff=1*nanometer,ewaldErrorTolerance=0.00001,
vdwCutoff=1.2*nanometer, polarization ='direct')
When using the AMOEBA force field, OpenMM allows the induced dipoles to be
calculated in any of three different ways. The slowest but potentially most
accurate method is to iterate the calculation until the dipoles converge to a
specified tolerance. To select this, specify :code:`polarization='mutual'`.
Use the :code:`mutualInducedTargetEpsilon` option to select the tolerance; for
most situations, a value of 0.00001 works well. Alternatively you can specify
:code:`polarization='extrapolated'`. This uses an analytic approximation
:cite:`Simmonett2015` to estimate what the fully converged dipoles will be without
actually continuing the calculation to convergence. In many cases this can be
significantly faster with only a small loss in accuracy. Finally, you can
specify :code:`polarization='direct'` to use the direct polarization
approximation, in which induced dipoles depend only on the fixed multipoles, not
on other induced dipoles. This is even faster, but it produces very different
forces from mutual polarization, so it should only be used with force fields
that have been specifically parameterized for use with this approximation.
Here are examples of using each method:
::
system = forcefield.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
vdwCutoff=1.2*nanometer, polarization='mutual', mutualInducedTargetEpsilon=0.00001)
system = forcefield.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
vdwCutoff=1.2*nanometer, polarization='extrapolated')
system = forcefield.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
vdwCutoff=1.2*nanometer, polarization='direct')
Implicit Solvent and Solute Dielectrics
......
docs-source/usersguide/references.bib
View file @ 71de4b1b
......@@ -389,6 +389,17 @@
type = {Journal Article}
}
@article{Simmonett2015
author = {Simmonett, Andrew C. and Pickard, Frank C. and Shao, Yihan and Cheatham, Thomas E. and Brooks, Bernard R.},
title = {Efficient treatment of induced dipoles},
journal = {Journal of Chemical Physics},
year = {2015},
volume = {143},
number = {7},
pages = {074115},
type = {Journal Article}
}
@article{Sindhikara2009,
author = {Sindhikara, Daniel J. and Kim, Seonah and Voter,
Arthur F. and Roitberg, Adrian E.},
......
examples/benchmark.py
View file @ 71de4b1b
......@@ -36,22 +36,18 @@ def runOneTest(testName, options):
if amoeba:
constraints = None
epsilon = float(options.epsilon)
if epsilon == 0:
polarization = 'direct'
else:
polarization = 'mutual'
if explicit:
ff = app.ForceField('amoeba2009.xml')
pdb = app.PDBFile('5dfr_solv-cube_equil.pdb')
cutoff = 0.7*unit.nanometers
vdwCutoff = 0.9*unit.nanometers
system = ff.createSystem(pdb.topology, nonbondedMethod=app.PME, nonbondedCutoff=cutoff, vdwCutoff=vdwCutoff, constraints=constraints, ewaldErrorTolerance=0.00075, mutualInducedTargetEpsilon=epsilon, polarization=polarization)
system = ff.createSystem(pdb.topology, nonbondedMethod=app.PME, nonbondedCutoff=cutoff, vdwCutoff=vdwCutoff, constraints=constraints, ewaldErrorTolerance=0.00075, mutualInducedTargetEpsilon=epsilon, polarization=options.polarization)
else:
ff = app.ForceField('amoeba2009.xml', 'amoeba2009_gk.xml')
pdb = app.PDBFile('5dfr_minimized.pdb')
cutoff = 2.0*unit.nanometers
vdwCutoff = 1.2*unit.nanometers
system = ff.createSystem(pdb.topology, nonbondedMethod=app.NoCutoff, constraints=constraints, mutualInducedTargetEpsilon=epsilon, polarization=polarization)
system = ff.createSystem(pdb.topology, nonbondedMethod=app.NoCutoff, constraints=constraints, mutualInducedTargetEpsilon=epsilon, polarization=options.polarization)
for f in system.getForces():
if isinstance(f, mm.AmoebaMultipoleForce) or isinstance(f, mm.AmoebaVdwForce) or isinstance(f, mm.AmoebaGeneralizedKirkwoodForce) or isinstance(f, mm.AmoebaWcaDispersionForce):
f.setForceGroup(1)
......@@ -127,7 +123,8 @@ parser.add_option('--platform', dest='platform', choices=platformNames, help='na
parser.add_option('--test', dest='test', choices=('gbsa', 'rf', 'pme', 'amoebagk', 'amoebapme'), help='the test to perform: gbsa, rf, pme, amoebagk, or amoebapme [default: all]')
parser.add_option('--pme-cutoff', default='0.9', dest='cutoff', type='float', help='direct space cutoff for PME in nm [default: 0.9]')
parser.add_option('--seconds', default='60', dest='seconds', type='float', help='target simulation length in seconds [default: 60]')
parser.add_option('--mutual-epsilon', default='1e-4', dest='epsilon', type='float', help='mutual induced epsilon for AMOEBA [default: 1e-4]')
parser.add_option('--polarization', default='mutual', dest='polarization', choices=('direct', 'extrapolated', 'mutual'), help='the polarization method for AMOEBA: direct, extrapolated, or mutual [default: mutual]')
parser.add_option('--mutual-epsilon', default='1e-5', dest='epsilon', type='float', help='mutual induced epsilon for AMOEBA [default: 1e-5]')
parser.add_option('--heavy-hydrogens', action='store_true', default=False, dest='heavy', help='repartition mass to allow a larger time step')
parser.add_option('--device', default=None, dest='device', help='device index for CUDA or OpenCL')
parser.add_option('--precision', default='single', dest='precision', choices=('single', 'mixed', 'double'), help='precision mode for CUDA or OpenCL: single, mixed, or double [default: single]')
......
plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h
View file @ 71de4b1b
......@@ -9,7 +9,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Authors: Mark Friedrichs, Peter Eastman *
* Contributors: *
* *
......@@ -72,14 +72,25 @@ public:
enum PolarizationType {
/**
* Mutual polarization
* Full mutually induced polarization. The dipoles are iterated until the converge to the accuracy specified
* by getMutualInducedTargetEpsilon().
*/
Mutual = 0,
/**
* Direct polarization
* Direct polarization approximation. The induced dipoles depend only on the fixed multipoles, not on other
* induced dipoles.
*/
Direct = 1
Direct = 1,
/**
* Extrapolated perturbation theory approximation. The dipoles are iterated a few times, and then an analytic
* approximation is used to extrapolate to the fully converged values. Call setExtrapolationCoefficients()
* to set the coefficients used for the extrapolation. The default coefficients used in this release are
* [0, -0.3, 0, 1.3], but be aware that those may change in a future release.
*/
Extrapolated = 2
};
enum MultipoleAxisTypes { ZThenX = 0, Bisector = 1, ZBisect = 2, ThreeFold = 3, ZOnly = 4, NoAxisType = 5, LastAxisTypeIndex = 6 };
......@@ -298,6 +309,23 @@ public:
*/
void setMutualInducedTargetEpsilon(double inputMutualInducedTargetEpsilon);
/**
* Set the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the extrapolation
* algorithm for induced dipoles.
*
* @param coefficients a vector whose mth entry specifies the coefficient for mu_m. The length of this
* vector determines how many iterations are performed.
*
*/
void setExtrapolationCoefficients(const std::vector<double> &coefficients);
/**
* Get the coefficients for the mu_0, mu_1, mu_2, ..., mu_n terms in the extrapolation
* algorithm for induced dipoles. In this release, the default values for the coefficients are
* [0, -0.3, 0, 1.3], but be aware that those may change in a future release.
*/
const std::vector<double>& getExtrapolationCoefficients() 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 +412,8 @@ private:
int pmeBSplineOrder;
std::vector<int> pmeGridDimension;
int mutualInducedMaxIterations;
std::vector<double> extrapolationCoefficients;
double mutualInducedTargetEpsilon;
double scalingDistanceCutoff;
double electricConstant;
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp
View file @ 71de4b1b
......@@ -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-2009 Stanford University and the Authors. *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Authors: *
* Contributors: *
* *
......@@ -43,6 +43,10 @@ AmoebaMultipoleForce::AmoebaMultipoleForce() : nonbondedMethod(NoCutoff), polari
mutualInducedTargetEpsilon(1.0e-02), scalingDistanceCutoff(100.0), electricConstant(138.9354558456), aewald(0.0) {
pmeGridDimension.resize(3);
pmeGridDimension[0] = pmeGridDimension[1] = pmeGridDimension[2];
extrapolationCoefficients.push_back(0.0);
extrapolationCoefficients.push_back(-0.3);
extrapolationCoefficients.push_back(0.0);
extrapolationCoefficients.push_back(1.3);
}
AmoebaMultipoleForce::NonbondedMethod AmoebaMultipoleForce::getNonbondedMethod() const {
......@@ -61,6 +65,14 @@ void AmoebaMultipoleForce::setPolarizationType(AmoebaMultipoleForce::Polarizatio
polarizationType = type;
}
void AmoebaMultipoleForce::setExtrapolationCoefficients(const std::vector<double> &coefficients) {
extrapolationCoefficients = coefficients;
}
const std::vector<double> & AmoebaMultipoleForce::getExtrapolationCoefficients() const {
return extrapolationCoefficients;
}
double AmoebaMultipoleForce::getCutoffDistance() const {
return cutoffDistance;
}
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h
View file @ 71de4b1b
......@@ -385,14 +385,17 @@ private:
const char* getSortKey() const {return "value.y";}
};
void initializeScaleFactors();
void computeInducedField(void** recipBoxVectorPointer);
bool iterateDipolesByDIIS(int iteration);
void computeExtrapolatedDipoles(void** recipBoxVectorPointer);
void ensureMultipolesValid(ContextImpl& context);
template <class T, class T4, class M4> void computeSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
int numMultipoles, maxInducedIterations;
int numMultipoles, maxInducedIterations, maxExtrapolationOrder;
int fixedFieldThreads, inducedFieldThreads, electrostaticsThreads;
int gridSizeX, gridSizeY, gridSizeZ;
double alpha, inducedEpsilon;
bool usePME, hasQuadrupoles, hasInitializedScaleFactors, hasInitializedFFT, multipolesAreValid;
AmoebaMultipoleForce::PolarizationType polarizationType;
CudaContext& cu;
const System& system;
std::vector<int3> covalentFlagValues;
......@@ -422,6 +425,18 @@ private:
CudaArray* prevErrors;
CudaArray* diisMatrix;
CudaArray* diisCoefficients;
CudaArray* extrapolatedDipole;
CudaArray* extrapolatedDipolePolar;
CudaArray* extrapolatedDipoleGk;
CudaArray* extrapolatedDipoleGkPolar;
CudaArray* inducedDipoleFieldGradient;
CudaArray* inducedDipoleFieldGradientPolar;
CudaArray* inducedDipoleFieldGradientGk;
CudaArray* inducedDipoleFieldGradientGkPolar;
CudaArray* extrapolatedDipoleFieldGradient;
CudaArray* extrapolatedDipoleFieldGradientPolar;
CudaArray* extrapolatedDipoleFieldGradientGk;
CudaArray* extrapolatedDipoleFieldGradientGkPolar;
CudaArray* polarizability;
CudaArray* covalentFlags;
CudaArray* polarizationGroupFlags;
......@@ -444,6 +459,7 @@ private:
CUfunction pmeGridIndexKernel, pmeSpreadFixedMultipolesKernel, pmeSpreadInducedDipolesKernel, pmeFinishSpreadChargeKernel, pmeConvolutionKernel;
CUfunction pmeFixedPotentialKernel, pmeInducedPotentialKernel, pmeFixedForceKernel, pmeInducedForceKernel, pmeRecordInducedFieldDipolesKernel, computePotentialKernel;
CUfunction recordDIISDipolesKernel, buildMatrixKernel;
CUfunction initExtrapolatedKernel, iterateExtrapolatedKernel, computeExtrapolatedKernel, addExtrapolatedGradientKernel;
CUfunction pmeTransformMultipolesKernel, pmeTransformPotentialKernel;
CudaCalcAmoebaGeneralizedKirkwoodForceKernel* gkKernel;
static const int PmeOrder = 5;
......@@ -512,6 +528,7 @@ private:
const System& system;
bool includeSurfaceArea, hasInitializedKernels;
int computeBornSumThreads, gkForceThreads, chainRuleThreads, ediffThreads;
AmoebaMultipoleForce::PolarizationType polarizationType;
std::map<std::string, std::string> defines;
CudaArray* params;
CudaArray* bornSum;
......
plugins/amoeba/platforms/cuda/src/kernels/gkEDiffPairForce.cu
View file @ 71de4b1b
......@@ -365,7 +365,22 @@ __device__ void computeOneEDiffInteractionT3(AtomData4& atom1, volatile AtomData
// correction to convert mutual to direct polarization force
#ifdef DIRECT_POLARIZATION
#ifdef MUTUAL_POLARIZATION
real findmp1 = uscale*(scip2*ddsc3_1 - ddsc5_1*(sci3*scip4+scip3*sci4));
real findmp2 = uscale*(scip2*ddsc3_2 - ddsc5_2*(sci3*scip4+scip3*sci4));
real findmp3 = uscale*(scip2*ddsc3_3 - ddsc5_3*(sci3*scip4+scip3*sci4));
ftm2i1 -= 0.5f*findmp1;
ftm2i2 -= 0.5f*findmp2;
ftm2i3 -= 0.5f*findmp3;
real sci3X = sci3 - sci3Y;
real sci4X = sci4 - sci4Y;
real scip3X = scip3 - scip3Y;
real scip4X = scip4 - scip4Y;
ftm2i1 += 0.5f*uscale*(-ddsc5_1*(sci3X*scip4X+scip3X*sci4X));
ftm2i2 += 0.5f*uscale*(-ddsc5_2*(sci3X*scip4X+scip3X*sci4X));
ftm2i3 += 0.5f*uscale*(-ddsc5_3*(sci3X*scip4X+scip3X*sci4X));
#else
real gfd = (scip2*scale3i - 5*rr2*(scip3*sci4+sci3*scip4)*scale5i);
real fdir1 = gfd*xr + scale5i* (sci4*atom1.inducedDipolePolarS.x+scip4*atom1.inducedDipoleS.x + sci3*atom2.inducedDipolePolarS.x+scip3*atom2.inducedDipoleS.x);
real fdir2 = gfd*yr + scale5i* (sci4*atom1.inducedDipolePolarS.y+scip4*atom1.inducedDipoleS.y + sci3*atom2.inducedDipolePolarS.y+scip3*atom2.inducedDipoleS.y);
......@@ -385,21 +400,6 @@ __device__ void computeOneEDiffInteractionT3(AtomData4& atom1, volatile AtomData
ftm2i1 += 0.5f*fdir1;
ftm2i2 += 0.5f*fdir2;
ftm2i3 += 0.5f*fdir3;
#else
real findmp1 = uscale*(scip2*ddsc3_1 - ddsc5_1*(sci3*scip4+scip3*sci4));
real findmp2 = uscale*(scip2*ddsc3_2 - ddsc5_2*(sci3*scip4+scip3*sci4));
real findmp3 = uscale*(scip2*ddsc3_3 - ddsc5_3*(sci3*scip4+scip3*sci4));
ftm2i1 -= 0.5f*findmp1;
ftm2i2 -= 0.5f*findmp2;
ftm2i3 -= 0.5f*findmp3;
real sci3X = sci3 - sci3Y;
real sci4X = sci4 - sci4Y;
real scip3X = scip3 - scip3Y;
real scip4X = scip4 - scip4Y;
ftm2i1 += 0.5f*uscale*(-ddsc5_1*(sci3X*scip4X+scip3X*sci4X));
ftm2i2 += 0.5f*uscale*(-ddsc5_2*(sci3X*scip4X+scip3X*sci4X));
ftm2i3 += 0.5f*uscale*(-ddsc5_3*(sci3X*scip4X+scip3X*sci4X));
#endif
#endif
......
plugins/amoeba/platforms/cuda/src/kernels/multipoleElectrostatics.cu
View file @ 71de4b1b
......@@ -327,7 +327,7 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, bool has
real iEIY = qiUinpI.x*Vijp[1] + qiUindI.x*Vijd[1] - qiUinpI.y*Vijp[0] - qiUindI.y*Vijd[0];
real iEJY = qiUinpJ.x*Vjip[1] + qiUindJ.x*Vjid[1] - qiUinpJ.y*Vjip[0] - qiUindJ.y*Vjid[0];
#ifdef USE_MUTUAL_POLARIZATION
#ifdef MUTUAL_POLARIZATION
// Uind-Uind terms (m=0)
real eCoef = -4*rInvVec[3]*thole_d0;
real dCoef = 6*rInvVec[4]*dthole_d0;
......
plugins/amoeba/platforms/cuda/src/kernels/multipolePme.cu
View file @ 71de4b1b
......@@ -643,8 +643,8 @@ extern "C" __global__ void computeInducedPotentialFromGrid(const real2* __restri
// Process the atoms in spatially sorted order. This improves cache performance when loading
// the grid values.
for (int i = blockIdx.x*blockDim.x+threadIdx.x; i < NUM_ATOMS; i += blockDim.x*gridDim.x) {
int m = pmeAtomGridIndex[i].x;
for (int atom = blockIdx.x*blockDim.x+threadIdx.x; atom < NUM_ATOMS; atom += blockDim.x*gridDim.x) {
int m = pmeAtomGridIndex[atom].x;
real4 pos = posq[m];
pos -= periodicBoxVecZ*floor(pos.z*recipBoxVecZ.z+0.5f);
pos -= periodicBoxVecY*floor(pos.y*recipBoxVecY.z+0.5f);
......@@ -1051,7 +1051,7 @@ extern "C" __global__ void computeInducedDipoleForceAndEnergy(real4* __restrict_
f.x += (inducedDipole[k]+inducedDipolePolar[k])*phi[i+NUM_ATOMS*j1];
f.y += (inducedDipole[k]+inducedDipolePolar[k])*phi[i+NUM_ATOMS*j2];
f.z += (inducedDipole[k]+inducedDipolePolar[k])*phi[i+NUM_ATOMS*j3];
#ifndef DIRECT_POLARIZATION
#ifdef MUTUAL_POLARIZATION
f.x += (inducedDipole[k]*phip[i+NUM_ATOMS*j1] + inducedDipolePolar[k]*phid[i+NUM_ATOMS*j1]);
f.y += (inducedDipole[k]*phip[i+NUM_ATOMS*j2] + inducedDipolePolar[k]*phid[i+NUM_ATOMS*j2]);
f.z += (inducedDipole[k]*phip[i+NUM_ATOMS*j3] + inducedDipolePolar[k]*phid[i+NUM_ATOMS*j3]);
......@@ -1073,8 +1073,12 @@ extern "C" __global__ void computeInducedDipoleForceAndEnergy(real4* __restrict_
energyBuffer[blockIdx.x*blockDim.x+threadIdx.x] += 0.25f*EPSILON_FACTOR*energy;
}
extern "C" __global__ void recordInducedFieldDipoles(const real* __restrict__ phid, real* const __restrict__ phip,
long long* __restrict__ inducedField, long long* __restrict__ inducedFieldPolar, real3 recipBoxVecX, real3 recipBoxVecY, real3 recipBoxVecZ) {
extern "C" __global__ void recordInducedFieldDipoles(const real* __restrict__ phid, real* const __restrict__ phip, long long* __restrict__ inducedField,
long long* __restrict__ inducedFieldPolar, const real* __restrict__ inducedDipole, const real* __restrict__ inducedDipolePolar,
#ifdef EXTRAPOLATED_POLARIZATION
unsigned long long* __restrict__ fieldGradient, unsigned long long* __restrict__ fieldGradientPolar,
#endif
real3 recipBoxVecX, real3 recipBoxVecY, real3 recipBoxVecZ) {
__shared__ real fracToCart[3][3];
if (threadIdx.x == 0) {
fracToCart[0][0] = GRID_SIZE_X*recipBoxVecX.x;
......@@ -1088,12 +1092,62 @@ extern "C" __global__ void recordInducedFieldDipoles(const real* __restrict__ ph
fracToCart[2][2] = GRID_SIZE_Z*recipBoxVecZ.z;
}
__syncthreads();
real selfDipoleScale = (4/(real) 3)*(EWALD_ALPHA*EWALD_ALPHA*EWALD_ALPHA)/SQRT_PI;
for (int i = blockIdx.x*blockDim.x+threadIdx.x; i < NUM_ATOMS; i += blockDim.x*gridDim.x) {
inducedField[i] -= (long long) (0x100000000*(phid[i+NUM_ATOMS]*fracToCart[0][0] + phid[i+NUM_ATOMS*2]*fracToCart[0][1] + phid[i+NUM_ATOMS*3]*fracToCart[0][2]));
inducedField[i+PADDED_NUM_ATOMS] -= (long long) (0x100000000*(phid[i+NUM_ATOMS]*fracToCart[1][0] + phid[i+NUM_ATOMS*2]*fracToCart[1][1] + phid[i+NUM_ATOMS*3]*fracToCart[1][2]));
inducedField[i+PADDED_NUM_ATOMS*2] -= (long long) (0x100000000*(phid[i+NUM_ATOMS]*fracToCart[2][0] + phid[i+NUM_ATOMS*2]*fracToCart[2][1] + phid[i+NUM_ATOMS*3]*fracToCart[2][2]));
inducedFieldPolar[i] -= (long long) (0x100000000*(phip[i+NUM_ATOMS]*fracToCart[0][0] + phip[i+NUM_ATOMS*2]*fracToCart[0][1] + phip[i+NUM_ATOMS*3]*fracToCart[0][2]));
inducedFieldPolar[i+PADDED_NUM_ATOMS] -= (long long) (0x100000000*(phip[i+NUM_ATOMS]*fracToCart[1][0] + phip[i+NUM_ATOMS*2]*fracToCart[1][1] + phip[i+NUM_ATOMS*3]*fracToCart[1][2]));
inducedFieldPolar[i+PADDED_NUM_ATOMS*2] -= (long long) (0x100000000*(phip[i+NUM_ATOMS]*fracToCart[2][0] + phip[i+NUM_ATOMS*2]*fracToCart[2][1] + phip[i+NUM_ATOMS*3]*fracToCart[2][2]));
inducedField[i] -= (long long) (0x100000000*(phid[i+NUM_ATOMS]*fracToCart[0][0] + phid[i+NUM_ATOMS*2]*fracToCart[0][1] + phid[i+NUM_ATOMS*3]*fracToCart[0][2] - selfDipoleScale*inducedDipole[3*i]));
inducedField[i+PADDED_NUM_ATOMS] -= (long long) (0x100000000*(phid[i+NUM_ATOMS]*fracToCart[1][0] + phid[i+NUM_ATOMS*2]*fracToCart[1][1] + phid[i+NUM_ATOMS*3]*fracToCart[1][2] - selfDipoleScale*inducedDipole[3*i+1]));
inducedField[i+PADDED_NUM_ATOMS*2] -= (long long) (0x100000000*(phid[i+NUM_ATOMS]*fracToCart[2][0] + phid[i+NUM_ATOMS*2]*fracToCart[2][1] + phid[i+NUM_ATOMS*3]*fracToCart[2][2] - selfDipoleScale*inducedDipole[3*i+2]));
inducedFieldPolar[i] -= (long long) (0x100000000*(phip[i+NUM_ATOMS]*fracToCart[0][0] + phip[i+NUM_ATOMS*2]*fracToCart[0][1] + phip[i+NUM_ATOMS*3]*fracToCart[0][2] - selfDipoleScale*inducedDipolePolar[3*i]));
inducedFieldPolar[i+PADDED_NUM_ATOMS] -= (long long) (0x100000000*(phip[i+NUM_ATOMS]*fracToCart[1][0] + phip[i+NUM_ATOMS*2]*fracToCart[1][1] + phip[i+NUM_ATOMS*3]*fracToCart[1][2] - selfDipoleScale*inducedDipolePolar[3*i+1]));
inducedFieldPolar[i+PADDED_NUM_ATOMS*2] -= (long long) (0x100000000*(phip[i+NUM_ATOMS]*fracToCart[2][0] + phip[i+NUM_ATOMS*2]*fracToCart[2][1] + phip[i+NUM_ATOMS*3]*fracToCart[2][2] - selfDipoleScale*inducedDipolePolar[3*i+2]));
#ifdef EXTRAPOLATED_POLARIZATION
// Compute and store the field gradients for later use.
real EmatD[3][3] = {
{phid[i+NUM_ATOMS*4], phid[i+NUM_ATOMS*7], phid[i+NUM_ATOMS*8]},
{phid[i+NUM_ATOMS*7], phid[i+NUM_ATOMS*5], phid[i+NUM_ATOMS*9]},
{phid[i+NUM_ATOMS*8], phid[i+NUM_ATOMS*9], phid[i+NUM_ATOMS*6]}
};
real Exx = 0, Eyy = 0, Ezz = 0, Exy = 0, Exz = 0, Eyz = 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];
}
}
atomicAdd(&fieldGradient[6*i+0], static_cast<unsigned long long>((long long) (-Exx*0x100000000)));
atomicAdd(&fieldGradient[6*i+1], static_cast<unsigned long long>((long long) (-Eyy*0x100000000)));
atomicAdd(&fieldGradient[6*i+2], static_cast<unsigned long long>((long long) (-Ezz*0x100000000)));
atomicAdd(&fieldGradient[6*i+3], static_cast<unsigned long long>((long long) (-Exy*0x100000000)));
atomicAdd(&fieldGradient[6*i+4], static_cast<unsigned long long>((long long) (-Exz*0x100000000)));
atomicAdd(&fieldGradient[6*i+5], static_cast<unsigned long long>((long long) (-Eyz*0x100000000)));
real EmatP[3][3] = {
{phip[i+NUM_ATOMS*4], phip[i+NUM_ATOMS*7], phip[i+NUM_ATOMS*8]},
{phip[i+NUM_ATOMS*7], phip[i+NUM_ATOMS*5], phip[i+NUM_ATOMS*9]},
{phip[i+NUM_ATOMS*8], phip[i+NUM_ATOMS*9], phip[i+NUM_ATOMS*6]}
};
Exx = 0; Eyy = 0; Ezz = 0; Exy = 0; Exz = 0; Eyz = 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];
}
}
atomicAdd(&fieldGradientPolar[6*i+0], static_cast<unsigned long long>((long long) (-Exx*0x100000000)));
atomicAdd(&fieldGradientPolar[6*i+1], static_cast<unsigned long long>((long long) (-Eyy*0x100000000)));
atomicAdd(&fieldGradientPolar[6*i+2], static_cast<unsigned long long>((long long) (-Ezz*0x100000000)));
atomicAdd(&fieldGradientPolar[6*i+3], static_cast<unsigned long long>((long long) (-Exy*0x100000000)));
atomicAdd(&fieldGradientPolar[6*i+4], static_cast<unsigned long long>((long long) (-Exz*0x100000000)));
atomicAdd(&fieldGradientPolar[6*i+5], static_cast<unsigned long long>((long long) (-Eyz*0x100000000)));
#endif
}
}
plugins/amoeba/platforms/cuda/src/kernels/pmeMultipoleElectrostatics.cu
View file @ 71de4b1b
......@@ -365,7 +365,7 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, bool has
real iEIY = qiUinpI.x*Vijp[1] + qiUindI.x*Vijd[1] - qiUinpI.y*Vijp[0] - qiUindI.y*Vijd[0];
real iEJY = qiUinpJ.x*Vjip[1] + qiUindJ.x*Vjid[1] - qiUinpJ.y*Vjip[0] - qiUindJ.y*Vjid[0];
#ifdef USE_MUTUAL_POLARIZATION
#ifdef MUTUAL_POLARIZATION
// Uind-Uind terms (m=0)
real eCoef = -fourThirds*rInvVec[3]*(3*(thole_d0 + bVec[3]) + alphaRVec[3]*X);
real dCoef = rInvVec[4]*(6*(dthole_d0 + bVec[3]) + 4*alphaRVec[5]*X);
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
View file @ 71de4b1b
......@@ -50,10 +50,36 @@
 
 
using namespace OpenMM;
using namespace std;
const double TOL = 1e-4;
 
extern "C" void registerAmoebaCudaKernelFactories();
 
static void checkFiniteDifferences(vector<Vec3> 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) forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
vector<Vec3> positions2(forces.size()), positions3(forces.size());
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = 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(state3.getPotentialEnergy()+norm*stepSize, state2.getPotentialEnergy(), 1e-4);
}
// setup for 2 ammonia molecules
 
static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce,
......@@ -289,6 +315,10 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
energy = state.getPotentialEnergy();
// Check that the forces and energy are consistent.
checkFiniteDifferences(forces, context, positions);
}
 
// setup for villin
......@@ -6972,6 +7002,10 @@ static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Polariz
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
energy = state.getPotentialEnergy();
// Check that the forces and energy are consistent.
checkFiniteDifferences(forces, context, positions);
}
 
// compare forces and energies
......@@ -7049,6 +7083,25 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization() {
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
}
 
static void testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization() {
std::string testName = "testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization";
int numberOfParticles = 8;
std::vector<Vec3> forces;
double energy;
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("CUDA"));
// We don't have reference values for this case, but at least check that force and energy are consistent.
getForcesEnergyMultipoleAmmonia(context, forces, energy);
}
// test GK mutual polarization for system comprised of two ammonia molecules
 
static void testGeneralizedKirkwoodAmmoniaMutualPolarization() {
......@@ -7765,6 +7818,19 @@ static void testGeneralizedKirkwoodVillinDirectPolarization() {
compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
}
 
static void testGeneralizedKirkwoodVillinExtrapolatedPolarization() {
std::string testName = "testGeneralizedKirkwoodVillinExtrapolatedPolarization";
int numberOfParticles = 596;
std::vector<Vec3> forces;
double energy;
// We don't have reference values for this case, but at least check that force and energy are consistent.
setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Extrapolated, 0, forces, energy);
}
// test GK mutual polarization for villin system
 
static void testGeneralizedKirkwoodVillinMutualPolarization() {
......@@ -8401,9 +8467,10 @@ int main(int argc, char* argv[]) {
 
testGeneralizedKirkwoodAmmoniaDirectPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarization();
testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm();
testGeneralizedKirkwoodVillinDirectPolarization();
testGeneralizedKirkwoodVillinExtrapolatedPolarization();
testGeneralizedKirkwoodVillinMutualPolarization();
 
} catch(const std::exception& e) {
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 71de4b1b
......@@ -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::Extrapolated) {
extrapolationCoefficients = force.getExtrapolationCoefficients();
}
// PME
......@@ -667,6 +669,9 @@ AmoebaReferenceMultipoleForce* ReferenceCalcAmoebaMultipoleForceKernel::setupAmo
amoebaReferenceMultipoleForce->setMaximumMutualInducedDipoleIterations(mutualInducedMaxIterations);
} else if (polarizationType == AmoebaMultipoleForce::Direct) {
amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Direct);
} else if (polarizationType == AmoebaMultipoleForce::Extrapolated) {
amoebaReferenceMultipoleForce->setPolarizationType(AmoebaReferenceMultipoleForce::Extrapolated);
amoebaReferenceMultipoleForce->setExtrapolationCoefficients(extrapolationCoefficients);
} else {
throw OpenMMException("Polarization type not recognzied.");
}
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 71de4b1b
......@@ -432,6 +432,7 @@ private:
int mutualInducedMaxIterations;
RealOpenMM mutualInducedTargetEpsilon;
std::vector<double> extrapolationCoefficients;
bool usePme;
RealOpenMM alphaEwald;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.h
View file @ 71de4b1b
......@@ -345,11 +345,16 @@ public:
*/
Mutual = 0,
/**
/**
* Direct polarization
*/
Direct = 1
};
Direct = 1,
/**
* Extrapolated perturbation theory
*/
Extrapolated = 2
};
/**
* Constructor
......@@ -422,6 +427,15 @@ public:
*/
RealOpenMM getMutualInducedDipoleEpsilon() const;
/**
* Set the coefficients for the µ_0, µ_1, µ_2, µ_n terms in the extrapolation
* theory algorithm for induced dipoles
*
* @param optCoefficients a vector whose mth entry specifies the coefficient for µ_m
*
*/
void setExtrapolationCoefficients(const std::vector<RealOpenMM> &coefficients);
/**
* Set the target epsilon for converging mutual induced dipoles.
*
......@@ -624,10 +638,13 @@ protected:
* Helper class used in calculating induced dipoles
*/
struct UpdateInducedDipoleFieldStruct {
UpdateInducedDipoleFieldStruct(std::vector<OpenMM::RealVec>& inputFixed_E_Field, std::vector<OpenMM::RealVec>& inputInducedDipoles);
UpdateInducedDipoleFieldStruct(std::vector<OpenMM::RealVec>& inputFixed_E_Field, std::vector<OpenMM::RealVec>& inputInducedDipoles, std::vector<std::vector<RealVec> >& extrapolatedDipoles, std::vector<std::vector<RealOpenMM> >& extrapolatedDipoleFieldGradient);
std::vector<OpenMM::RealVec>* fixedMultipoleField;
std::vector<OpenMM::RealVec>* inducedDipoles;
std::vector<std::vector<RealVec> >* extrapolatedDipoles;
std::vector<std::vector<RealOpenMM> >* extrapolatedDipoleFieldGradient;
std::vector<OpenMM::RealVec> inducedDipoleField;
std::vector<std::vector<RealOpenMM> > inducedDipoleFieldGradient;
};
unsigned int _numParticles;
......@@ -651,10 +668,17 @@ 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> _extrapolationCoefficients;
std::vector<RealOpenMM> _extPartCoefficients;
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 extrapolated 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 convergeInduceDipolesByExtrapolation(const std::vector<MultipoleParticleData>& particleData,
std::vector<UpdateInducedDipoleFieldStruct>& calculateInducedDipoleField);
/**
* Converge induced dipoles.
*
......@@ -1191,6 +1223,10 @@ private:
std::vector<RealVec> _gkField;
std::vector<RealVec> _inducedDipoleS;
std::vector<RealVec> _inducedDipolePolarS;
std::vector<std::vector<RealVec> > _ptDipolePS;
std::vector<std::vector<RealVec> > _ptDipoleDS;
std::vector<std::vector<RealOpenMM> > _ptDipoleFieldGradientPS;
std::vector<std::vector<RealOpenMM> > _ptDipoleFieldGradientDS;
int _includeCavityTerm;
RealOpenMM _probeRadius;
......
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp
View file @ 71de4b1b
......@@ -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-2012 Stanford University and the Authors. *
* Portions copyright (c) 2008-2015 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* Contributors: *
* *
......@@ -50,11 +50,36 @@
 
 
using namespace OpenMM;
using namespace std;
 
extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
 
const double TOL = 1e-4;
 
static void checkFiniteDifferences(vector<Vec3> 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) forces.size(); ++i)
norm += forces[i].dot(forces[i]);
norm = std::sqrt(norm);
const double stepSize = 1e-3;
double step = 0.5*stepSize/norm;
vector<Vec3> positions2(forces.size()), positions3(forces.size());
for (int i = 0; i < (int) positions.size(); ++i) {
Vec3 p = positions[i];
Vec3 f = 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);
}
// setup for 2 ammonia molecules
 
static void setupMultipoleAmmonia(System& system, AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce,
......@@ -291,6 +316,10 @@ static void getForcesEnergyMultipoleAmmonia(Context& context, std::vector<Vec3>&
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
energy = state.getPotentialEnergy();
// Check that the forces and energy are consistent.
checkFiniteDifferences(forces, context, positions);
}
 
// setup for villin
......@@ -6974,6 +7003,10 @@ static void setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Polariz
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
energy = state.getPotentialEnergy();
// Check that the forces and energy are consistent.
checkFiniteDifferences(forces, context, positions);
}
 
// compare forces and energies
......@@ -7052,6 +7085,25 @@ static void testGeneralizedKirkwoodAmmoniaDirectPolarization() {
compareForcesEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
}
 
static void testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization() {
std::string testName = "testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization";
int numberOfParticles = 8;
std::vector<Vec3> forces;
double energy;
System system;
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
setupMultipoleAmmonia(system, amoebaGeneralizedKirkwoodForce, AmoebaMultipoleForce::Direct, 0);
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName("Reference"));
// We don't have reference values for this case, but at least check that force and energy are consistent.
getForcesEnergyMultipoleAmmonia(context, forces, energy);
}
// test GK mutual polarization for system comprised of two ammonia molecules
 
static void testGeneralizedKirkwoodAmmoniaMutualPolarization() {
......@@ -7770,6 +7822,19 @@ static void testGeneralizedKirkwoodVillinDirectPolarization() {
compareForceNormsEnergy(testName, expectedEnergy, energy, expectedForces, forces, tolerance);
}
 
static void testGeneralizedKirkwoodVillinExtrapolatedPolarization() {
std::string testName = "testGeneralizedKirkwoodVillinExtrapolatedPolarization";
int numberOfParticles = 596;
std::vector<Vec3> forces;
double energy;
// We don't have reference values for this case, but at least check that force and energy are consistent.
setupAndGetForcesEnergyMultipoleVillin(AmoebaMultipoleForce::Extrapolated, 0, forces, energy);
}
// test GK mutual polarization for villin system
 
static void testGeneralizedKirkwoodVillinMutualPolarization() {
......@@ -8405,8 +8470,10 @@ int main(int numberOfArguments, char* argv[]) {
 
testGeneralizedKirkwoodAmmoniaMutualPolarization();
testGeneralizedKirkwoodAmmoniaDirectPolarization();
testGeneralizedKirkwoodAmmoniaExtrapolatedPolarization();
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm();
testGeneralizedKirkwoodVillinDirectPolarization();
testGeneralizedKirkwoodVillinExtrapolatedPolarization();
testGeneralizedKirkwoodVillinMutualPolarization();
 
}
......
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