Skip to content

GitLab

Commit 9dd055eb authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Add unit test for NonbondedSoftcoreForce 

Cleanup code associated w/ NonbondedSoftcoreForce
parent ab60bff4
Hide whitespace changes
Inline Side-by-side
Showing with 337 additions and 721 deletions
plugins/freeEnergy/platforms/reference/src/ReferenceFreeEnergyKernels.cpp
View file @ 9dd055eb
......@@ -155,7 +155,6 @@ void ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel::initialize(const Syste
nonbondedMethod = CalcNonbondedSoftcoreForceKernel::NonbondedSoftcoreMethod(force.getNonbondedMethod());
nonbondedCutoff = (RealOpenMM) force.getCutoffDistance();
//softCoreLJLambda = (RealOpenMM) force.getSoftCoreLJLambda();
Vec3 boxVectors[3];
system.getDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2]);
......@@ -168,42 +167,6 @@ void ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel::initialize(const Syste
else
neighborList = new NeighborList();
#if 0
if (nonbondedMethod == Ewald || nonbondedMethod == PME) {
RealOpenMM ewaldErrorTol = (RealOpenMM) force.getEwaldErrorTolerance();
ewaldAlpha = (RealOpenMM) (std::sqrt(-std::log(ewaldErrorTol))/nonbondedCutoff);
RealOpenMM mx = periodicBoxSize[0]/nonbondedCutoff;
RealOpenMM my = periodicBoxSize[1]/nonbondedCutoff;
RealOpenMM mz = periodicBoxSize[2]/nonbondedCutoff;
RealOpenMM pi = (RealOpenMM) 3.1415926535897932385;
kmax[0] = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
kmax[1] = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
kmax[2] = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
if (kmax[0]%2 == 0)
kmax[0]++;
if (kmax[1]%2 == 0)
kmax[1]++;
if (kmax[2]%2 == 0)
kmax[2]++;
}
if (nonbondedMethod == Ewald || nonbondedMethod == PME) {
RealOpenMM ewaldErrorTol = (RealOpenMM) force.getEwaldErrorTolerance();
ewaldAlpha = (RealOpenMM) (std::sqrt(-std::log(ewaldErrorTol))/nonbondedCutoff);
RealOpenMM mx = periodicBoxSize[0]/nonbondedCutoff;
RealOpenMM my = periodicBoxSize[1]/nonbondedCutoff;
RealOpenMM mz = periodicBoxSize[2]/nonbondedCutoff;
RealOpenMM pi = (RealOpenMM) 3.1415926535897932385;
kmax[0] = (int)std::ceil(-(mx/pi)*std::log(ewaldErrorTol));
kmax[1] = (int)std::ceil(-(my/pi)*std::log(ewaldErrorTol));
kmax[2] = (int)std::ceil(-(mz/pi)*std::log(ewaldErrorTol));
if (kmax[0]%2 == 0)
kmax[0]++;
if (kmax[1]%2 == 0)
kmax[1]++;
if (kmax[2]%2 == 0)
kmax[2]++;
}
#endif
rfDielectric = (RealOpenMM)force.getReactionFieldDielectric();
}
......@@ -214,21 +177,18 @@ double ReferenceFreeEnergyCalcNonbondedSoftcoreForceKernel::execute(ContextImpl&
RealOpenMM energy = 0;
ReferenceFreeEnergyLJCoulombSoftcoreIxn clj;
// clj.setSoftCoreLJLambda( softCoreLJLambda );
bool periodic = (nonbondedMethod == CutoffPeriodic);
bool ewald = (nonbondedMethod == Ewald);
bool pme = (nonbondedMethod == PME);
if (nonbondedMethod != NoCutoff) {
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodicBoxSize, periodic || ewald || pme, nonbondedCutoff, 0.0);
computeNeighborListVoxelHash(*neighborList, numParticles, posData, exclusions, periodicBoxSize, periodic, nonbondedCutoff, 0.0);
clj.setUseCutoff(nonbondedCutoff, *neighborList, rfDielectric);
}
if (periodic || ewald || pme)
if (periodic)
clj.setPeriodic(periodicBoxSize);
if (ewald)
clj.setUseEwald(ewaldAlpha, kmax[0], kmax[1], kmax[2]);
if (pme)
clj.setUsePME(ewaldAlpha);
clj.calculatePairIxn(numParticles, posData, particleParamArray, exclusionArray, 0, forceData, 0, &energy);
ReferenceBondForce refBondForce;
ReferenceFreeEnergyLJCoulomb14Softcore nonbonded14;
if (nonbondedMethod == CutoffNonPeriodic || nonbondedMethod == CutoffPeriodic)
......
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulomb14Softcore.cpp
View file @ 9dd055eb
......@@ -100,13 +100,11 @@ ReferenceFreeEnergyLJCoulomb14Softcore::~ReferenceFreeEnergyLJCoulomb14Softcore(
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulomb14Softcore::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM q2, RealOpenMM epsfac,
RealOpenMM* parameters ) const {
RealOpenMM q2, RealOpenMM epsfac,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceFreeEnergyLJCoulomb14Softcore::getDerivedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
......@@ -218,7 +216,7 @@ void ReferenceFreeEnergyLJCoulomb14Softcore::calculateBondIxn( int* atomIndices,
*totalEnergy += energy;
}
/**---------------------------------------------------------------------------------------
/**---------------------------------------------------------------------------------------
Calculate LJ pair ixn between two atoms
......@@ -231,7 +229,7 @@ void ReferenceFreeEnergyLJCoulomb14Softcore::calculateBondIxn( int* atomIndices,
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulomb14Softcore::calculateOneLJ14Ixn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
......
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.cpp
View file @ 9dd055eb
......@@ -27,16 +27,8 @@
#include <complex>
#include "../SimTKUtilities/SimTKOpenMMCommon.h"
#include "../SimTKUtilities/SimTKOpenMMLog.h"
#include "../SimTKUtilities/SimTKOpenMMUtilities.h"
#include "ReferenceFreeEnergyLJCoulombSoftcoreIxn.h"
#include "ReferenceForce.h"
#include "PME.h"
// In case we're using some primitive version of Visual Studio this will
// make sure that erf() and erfc() are defined.
//#include "MSVC_erfc.h"
#include "openmm/internal/MSVC_erfc.h"
using std::vector;
using OpenMM::RealVec;
......@@ -47,14 +39,7 @@ using OpenMM::RealVec;
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulombSoftcoreIxn::ReferenceFreeEnergyLJCoulombSoftcoreIxn( ) : cutoff(false), periodic(false), ewald(false), pme(false), softCoreLJLambda(1.0) {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceFreeEnergyLJCoulombSoftcoreIxn::ReferenceFreeEnergyLJCoulombSoftcoreIxn";
// ---------------------------------------------------------------------------------------
ReferenceFreeEnergyLJCoulombSoftcoreIxn::ReferenceFreeEnergyLJCoulombSoftcoreIxn( ) : cutoff(false), periodic(false), ewald(false), pme(false) {
}
/**---------------------------------------------------------------------------------------
......@@ -64,13 +49,6 @@ ReferenceFreeEnergyLJCoulombSoftcoreIxn::ReferenceFreeEnergyLJCoulombSoftcoreIxn
--------------------------------------------------------------------------------------- */
ReferenceFreeEnergyLJCoulombSoftcoreIxn::~ReferenceFreeEnergyLJCoulombSoftcoreIxn( ){
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceFreeEnergyLJCoulombSoftcoreIxn::~ReferenceFreeEnergyLJCoulombSoftcoreIxn";
// ---------------------------------------------------------------------------------------
}
/**---------------------------------------------------------------------------------------
......@@ -85,11 +63,11 @@ ReferenceFreeEnergyLJCoulombSoftcoreIxn::~ReferenceFreeEnergyLJCoulombSoftcoreIx
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::setUseCutoff( RealOpenMM distance, const OpenMM::NeighborList& neighbors, RealOpenMM solventDielectric ) {
cutoff = true;
cutoff = true;
cutoffDistance = distance;
neighborList = &neighbors;
krf = pow(cutoffDistance, (RealOpenMM)-3.0)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
crf = (1.0f/cutoffDistance)*(3.0f*solventDielectric)/(2.0f*solventDielectric+1.0f);
neighborList = &neighbors;
krf = pow(cutoffDistance, (RealOpenMM)-3.0)*(solventDielectric-1.0f)/(2.0f*solventDielectric+1.0f);
crf = (1.0f/cutoffDistance)*(3.0f*solventDielectric)/(2.0f*solventDielectric+1.0f);
}
/**---------------------------------------------------------------------------------------
......@@ -105,60 +83,17 @@ ReferenceFreeEnergyLJCoulombSoftcoreIxn::~ReferenceFreeEnergyLJCoulombSoftcoreIx
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::setPeriodic( RealVec& boxSize ) {
assert(cutoff);
assert(boxSize[0] >= 2.0*cutoffDistance);
assert(boxSize[1] >= 2.0*cutoffDistance);
assert(boxSize[2] >= 2.0*cutoffDistance);
periodic = true;
periodic = true;
periodicBoxSize[0] = boxSize[0];
periodicBoxSize[1] = boxSize[1];
periodicBoxSize[2] = boxSize[2];
}
/**---------------------------------------------------------------------------------------
Set the force to use Ewald summation.
@param alpha the Ewald separation parameter
@param kmaxx the largest wave vector in the x direction
@param kmaxy the largest wave vector in the y direction
@param kmaxz the largest wave vector in the z direction
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::setUseEwald(RealOpenMM alpha, int kmaxx, int kmaxy, int kmaxz) {
alphaEwald = alpha;
numRx = kmaxx;
numRy = kmaxy;
numRz = kmaxz;
ewald = true;
}
/**---------------------------------------------------------------------------------------
Set the force to use Particle-Mesh Ewald (PME) summation.
@param alpha the Ewald separation parameter
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::setUsePME(RealOpenMM alpha) {
alphaEwald = alpha;
pme = true;
}
/**---------------------------------------------------------------------------------------
Set the soft core LJ lambda
@param lambda the soft core LJ lambda
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::setSoftCoreLJLambda(RealOpenMM lambda) {
softCoreLJLambda = lambda;
}
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ Coulomb ixn
......@@ -166,7 +101,7 @@ ReferenceFreeEnergyLJCoulombSoftcoreIxn::~ReferenceFreeEnergyLJCoulombSoftcoreIx
@param c6 c6
@param c12 c12
@param q1 q1 charge atom 1
@param epsfac epsfacSqrt ????????????
@param epsfac epsfacSqrt
@param parameters output parameters:
parameter[SigIndex] = 0.5*( (c12/c6)**1/6 ) (sigma/2)
parameter[EpsIndex] = sqrt(c6*c6/c12) (2*sqrt(epsilon))
......@@ -175,13 +110,11 @@ ReferenceFreeEnergyLJCoulombSoftcoreIxn::~ReferenceFreeEnergyLJCoulombSoftcoreIx
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const {
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const {
// ---------------------------------------------------------------------------------------
// static const char* methodName = "\nReferenceFreeEnergyLJCoulombSoftcoreIxn::getDerivedParameters";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
......@@ -207,440 +140,6 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::getDerivedParameters( RealOpenMM c
parameters[QIndex] = epsfacSqrt*q1;
}
/**---------------------------------------------------------------------------------------
Calculate Ewald ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateEwaldIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
#if 0
typedef std::complex<RealOpenMM> d_complex;
static const RealOpenMM epsilon = 1.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
static const RealOpenMM twelve = 12.0;
int kmax = (ewald ? std::max(numRx, std::max(numRy,numRz)) : 0);
RealOpenMM factorEwald = -1 / (4*alphaEwald*alphaEwald);
RealOpenMM SQRT_PI = sqrt(PI);
RealOpenMM TWO_PI = 2.0 * PI;
RealOpenMM recipCoeff = (RealOpenMM)(4*PI/(periodicBoxSize[0] * periodicBoxSize[1] * periodicBoxSize[2]) /epsilon);
RealOpenMM totalSelfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
RealOpenMM recipEnergy = 0.0;
RealOpenMM totalRecipEnergy = 0.0;
RealOpenMM vdwEnergy = 0.0;
// **************************************************************************************
// SELF ENERGY
// **************************************************************************************
for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
RealOpenMM selfEwaldEnergy = atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex] * alphaEwald/SQRT_PI;
totalSelfEwaldEnergy -= selfEwaldEnergy;
if( energyByAtom ){
energyByAtom[atomID] -= selfEwaldEnergy;
}
}
if( totalEnergy ){
*totalEnergy += totalSelfEwaldEnergy;
}
// **************************************************************************************
// RECIPROCAL SPACE EWALD ENERGY AND FORCES
// **************************************************************************************
// PME
if (pme) {
pme_t pmedata; /* abstract handle for PME data */
int ngrid[3];
RealOpenMM virial[3][3];
/* PME grid dimensions.
* We typically want to set this as the spacing rather than absolute dimensions, but
* to be able to reproduce results from other programs (e.g. Gromacs) we need to be
* able to set exact grid dimenisions occasionally.
*/
ngrid[0] = 16;
ngrid[1] = 16;
ngrid[2] = 16;
pme_init(&pmedata,alphaEwald,numberOfAtoms,ngrid,4,1);
pme_exec(pmedata,atomCoordinates,forces,atomParameters,periodicBoxSize,&recipEnergy,virial);
if( totalEnergy )
*totalEnergy += recipEnergy;
if( energyByAtom )
for(int n = 0; n < numberOfAtoms; n++)
energyByAtom[n] += recipEnergy;
pme_destroy(pmedata);
}
// Ewald method
else if (ewald) {
// setup reciprocal box
RealOpenMM recipBoxSize[3] = { TWO_PI / periodicBoxSize[0], TWO_PI / periodicBoxSize[1], TWO_PI / periodicBoxSize[2]};
// setup K-vectors
#define EIR(x, y, z) eir[(x)*numberOfAtoms*3+(y)*3+z]
vector<d_complex> eir(kmax*numberOfAtoms*3);
vector<d_complex> tab_xy(numberOfAtoms);
vector<d_complex> tab_qxyz(numberOfAtoms);
if (kmax < 1) {
std::stringstream message;
message << " kmax < 1 , Aborting" << std::endl;
SimTKOpenMMLog::printError( message );
}
for(int i = 0; (i < numberOfAtoms); i++) {
for(int m = 0; (m < 3); m++)
EIR(0, i, m) = d_complex(1,0);
for(int m=0; (m<3); m++)
EIR(1, i, m) = d_complex(cos(atomCoordinates[i][m]*recipBoxSize[m]),
sin(atomCoordinates[i][m]*recipBoxSize[m]));
for(int j=2; (j<kmax); j++)
for(int m=0; (m<3); m++)
EIR(j, i, m) = EIR(j-1, i, m) * EIR(1, i, m);
}
// calculate reciprocal space energy and forces
int lowry = 0;
int lowrz = 1;
for(int rx = 0; rx < numRx; rx++) {
RealOpenMM kx = rx * recipBoxSize[0];
for(int ry = lowry; ry < numRy; ry++) {
RealOpenMM ky = ry * recipBoxSize[1];
if(ry >= 0) {
for(int n = 0; n < numberOfAtoms; n++)
tab_xy[n] = EIR(rx, n, 0) * EIR(ry, n, 1);
}
else {
for(int n = 0; n < numberOfAtoms; n++)
tab_xy[n]= EIR(rx, n, 0) * conj (EIR(-ry, n, 1));
}
for (int rz = lowrz; rz < numRz; rz++) {
if( rz >= 0) {
for( int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * EIR(rz, n, 2));
}
else {
for( int n = 0; n < numberOfAtoms; n++)
tab_qxyz[n] = atomParameters[n][QIndex] * (tab_xy[n] * conj(EIR(-rz, n, 2)));
}
RealOpenMM cs = 0.0f;
RealOpenMM ss = 0.0f;
for( int n = 0; n < numberOfAtoms; n++) {
cs += tab_qxyz[n].real();
ss += tab_qxyz[n].imag();
}
RealOpenMM kz = rz * recipBoxSize[2];
RealOpenMM k2 = kx * kx + ky * ky + kz * kz;
RealOpenMM ak = exp(k2*factorEwald) / k2;
for(int n = 0; n < numberOfAtoms; n++) {
RealOpenMM force = ak * (cs * tab_qxyz[n].imag() - ss * tab_qxyz[n].real());
forces[n][0] += 2 * recipCoeff * force * kx ;
forces[n][1] += 2 * recipCoeff * force * ky ;
forces[n][2] += 2 * recipCoeff * force * kz ;
}
recipEnergy = recipCoeff * ak * ( cs * cs + ss * ss);
totalRecipEnergy += recipEnergy;
if( totalEnergy )
*totalEnergy += recipEnergy;
if( energyByAtom )
for(int n = 0; n < numberOfAtoms; n++)
energyByAtom[n] += recipEnergy;
lowrz = 1 - numRz;
}
lowry = 1 - numRy;
}
}
}
else {
std::stringstream message;
message << " Wrong method for Ewald summation, Aborting" << std::endl;
SimTKOpenMMLog::printError( message );
}
// **************************************************************************************
// SHORT-RANGE ENERGY AND FORCES
// **************************************************************************************
RealOpenMM totalVdwEnergy = 0.0f;
RealOpenMM totalRealSpaceEwaldEnergy = 0.0f;
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
int ii = pair.first;
int jj = pair.second;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
RealOpenMM sig = atomParameters[ii][SigIndex] + atomParameters[jj][SigIndex];
RealOpenMM sig2 = inverseR*sig;
sig2 *= sig2;
RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM eps = atomParameters[ii][EpsIndex]*atomParameters[jj][EpsIndex];
dEdR += eps*( twelve*sig6 - six )*sig6*inverseR*inverseR;
// accumulate forces
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force;
forces[jj][kk] -= force;
}
// accumulate energies
realSpaceEwaldEnergy = (RealOpenMM) (atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
vdwEnergy = eps*(sig6-one)*sig6;
totalVdwEnergy += vdwEnergy;
totalRealSpaceEwaldEnergy += realSpaceEwaldEnergy;
if( energyByAtom ){
energyByAtom[ii] += realSpaceEwaldEnergy + vdwEnergy;
energyByAtom[jj] += realSpaceEwaldEnergy + vdwEnergy;
}
}
if( totalEnergy )
*totalEnergy += totalRealSpaceEwaldEnergy + totalVdwEnergy;
// Now subtract off the exclusions, since they were implicitly included in the reciprocal space sum.
RealOpenMM totalExclusionEnergy = 0.0f;
for (int i = 0; i < numberOfAtoms; i++)
for (int j = 1; j <= exclusions[i][0]; j++)
if (exclusions[i][j] > i) {
int ii = i;
int jj = exclusions[i][j];
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erf(alphaR) - 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
// accumulate forces
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] -= force;
forces[jj][kk] += force;
}
// accumulate energies
realSpaceEwaldEnergy = (RealOpenMM) (atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erf(alphaR));
totalExclusionEnergy += realSpaceEwaldEnergy;
if( energyByAtom ){
energyByAtom[ii] -= realSpaceEwaldEnergy;
energyByAtom[jj] -= realSpaceEwaldEnergy;
}
}
if( totalEnergy )
*totalEnergy -= totalExclusionEnergy;
// ***********************************************************************
#endif
}
/**---------------------------------------------------------------------------------------
Calculate PME ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePMEIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy) const {
#if 0
RealOpenMM SQRT_PI = sqrt(PI);
static const RealOpenMM one = 1.0;
RealOpenMM selfEwaldEnergy = 0.0;
RealOpenMM realSpaceEwaldEnergy = 0.0;
RealOpenMM recipEnergy = 0.0;
// **************************************************************************************
// SELF ENERGY
// **************************************************************************************
for( int atomID = 0; atomID < numberOfAtoms; atomID++ ){
selfEwaldEnergy = selfEwaldEnergy + atomParameters[atomID][QIndex]*atomParameters[atomID][QIndex];
}
selfEwaldEnergy = selfEwaldEnergy * alphaEwald/SQRT_PI ;
// **************************************************************************************
// RECIPROCAL SPACE EWALD ENERGY AND FORCES
// **************************************************************************************
pme_t pmedata; /* abstract handle for PME data */
int ngrid[3];
RealOpenMM virial[3][3];
/* PME grid dimensions.
* We typically want to set this as the spacing rather than absolute dimensions, but
* to be able to reproduce results from other programs (e.g. Gromacs) we need to be
* able to set exact grid dimenisions occasionally.
*/
ngrid[0] = 16;
ngrid[1] = 16;
ngrid[2] = 16;
pme_init(&pmedata,alphaEwald,numberOfAtoms,ngrid,4,1);
pme_exec(pmedata,atomCoordinates,forces,atomParameters,periodicBoxSize,&recipEnergy,virial);
// **************************************************************************************
// SHORT-RANGE ENERGY AND FORCES
// **************************************************************************************
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
for( int atomID1 = 0; atomID1 < numberOfAtoms; atomID1++ ){
for( int atomID2 = atomID1 + 1; atomID2 < numberOfAtoms; atomID2++ ){
ReferenceForce::getDeltaRPeriodic( atomCoordinates[atomID2], atomCoordinates[atomID1], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[atomID1][QIndex]*atomParameters[atomID2][QIndex]*inverseR*erfc(alphaEwald*r));
}
}
// allocate and initialize exclusion array
vector<int> exclusionIndices(numberOfAtoms);
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
for( int jj = 1; jj <= exclusions[ii][0]; jj++ ){
exclusionIndices[exclusions[ii][jj]] = ii;
}
// loop over atom pairs
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
ReferenceForce::getDeltaRPeriodic( atomCoordinates[jj], atomCoordinates[ii], periodicBoxSize, deltaR[0] );
RealOpenMM r = deltaR[0][ReferenceForce::RIndex];
RealOpenMM r2 = deltaR[0][ReferenceForce::R2Index];
RealOpenMM inverseR = one/(deltaR[0][ReferenceForce::RIndex]);
RealOpenMM alphaR = alphaEwald * r;
realSpaceEwaldEnergy =
(RealOpenMM)(realSpaceEwaldEnergy + atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*inverseR*erfc(alphaR));
RealOpenMM dEdR = atomParameters[ii][QIndex] * atomParameters[jj][QIndex] * inverseR * inverseR * inverseR;
dEdR = (RealOpenMM)(dEdR * (erfc(alphaR) + 2 * alphaR * exp ( - alphaR * alphaR) / SQRT_PI ));
for( int kk = 0; kk < 3; kk++ ){
RealOpenMM force = dEdR*deltaR[0][kk];
forces[ii][kk] += force;
forces[jj][kk] -= force;
}
}
}
}
// ***********************************************************************
#endif
}
/**---------------------------------------------------------------------------------------
Calculate LJ Coulomb pair ixn
......@@ -654,32 +153,25 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePMEIxn( int numberOfAtoms
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtoms, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
if (ewald || pme)
return calculateEwaldIxn(numberOfAtoms, atomCoordinates, atomParameters, exclusions, fixedParameters, forces, energyByAtom, totalEnergy);
if (cutoff) {
for (int i = 0; i < (int) neighborList->size(); i++) {
OpenMM::AtomPair pair = (*neighborList)[i];
calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, totalEnergy);
}
}
else {
// allocate and initialize exclusion array
} else {
int* exclusionIndices = new int[numberOfAtoms];
for( int ii = 0; ii < numberOfAtoms; ii++ ){
exclusionIndices[ii] = -1;
}
// allocate and initialize exclusion array
std::vector<int> exclusionIndices( numberOfAtoms, -1 );
for( int ii = 0; ii < numberOfAtoms; ii++ ){
// set exclusions
......@@ -693,12 +185,11 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtom
for( int jj = ii+1; jj < numberOfAtoms; jj++ ){
if( exclusionIndices[jj] != ii ){
calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
calculateOneIxn(ii, jj, atomCoordinates, atomParameters, forces, totalEnergy);
}
}
}
delete[] exclusionIndices;
}
}
......@@ -711,18 +202,13 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculatePairIxn( int numberOfAtom
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, vector<RealVec>& atomCoordinates,
RealOpenMM** atomParameters, vector<RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const {
// ---------------------------------------------------------------------------------------
static const std::string methodName = "\nReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn";
RealOpenMM* totalEnergy ) const {
// ---------------------------------------------------------------------------------------
......@@ -738,10 +224,6 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, v
static const int threeI = 3;
// debug flag
static const int debug = -1;
static const int LastAtomIndex = 2;
RealOpenMM deltaR[2][ReferenceForce::LastDeltaRIndex];
......@@ -790,7 +272,7 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, v
// accumulate energies
if( totalEnergy || energyByAtom ) {
if( totalEnergy ) {
if (cutoff)
energy += atomParameters[ii][QIndex]*atomParameters[jj][QIndex]*(inverseR+krf*r2-crf);
else
......@@ -798,10 +280,6 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, v
if( totalEnergy )
*totalEnergy += energy;
if( energyByAtom ){
energyByAtom[ii] += energy;
energyByAtom[jj] += energy;
}
}
}
......@@ -818,16 +296,10 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneIxn( int ii, int jj, v
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneLJIxn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "\nReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneLJIxn";
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
// constants -- reduce Visual Studio warnings regarding conversions between float & double
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM six = 6.0;
......@@ -855,17 +327,11 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneLJIxn( RealOpenMM inve
--------------------------------------------------------------------------------------- */
void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneSoftCoreLJIxn( RealOpenMM r, RealOpenMM sig, RealOpenMM eps,
RealOpenMM lambda,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
RealOpenMM lambda,
RealOpenMM* dEdR, RealOpenMM* energy ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "\nReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneSoftCoreLJIxn";
// ---------------------------------------------------------------------------------------
// constants -- reduce Visual Studio warnings regarding conversions between float & double
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM four = 4.0;
......@@ -873,13 +339,6 @@ void ReferenceFreeEnergyLJCoulombSoftcoreIxn::calculateOneSoftCoreLJIxn( RealOpe
static const RealOpenMM twelve = 12.0;
static const RealOpenMM alphaLJ = 0.5;
#if 0
RealOpenMM dEdROrig = 0.0;
RealOpenMM E_Orig = 0.0;
static int maxPrint = 0;
calculateOneLJIxn( one/r, sig, eps, &dEdROrig, &E_Orig );
#endif
// soft-core LJ energy = lambda*4*eps*[ 1/{alphaLJ*(1-lambda) + (r/sig)**6}**2 - 1/{alphaLJ*(1-lambda) + (r/sig)**6} ]
eps *= lambda;
......@@ -895,10 +354,4 @@ calculateOneLJIxn( one/r, sig, eps, &dEdROrig, &E_Orig );
*energy += eps*softcoreLJInv*( softcoreLJInv - one );
#if 0
if( maxPrint++ < 5 ){
printf( "r=%14.6e sig=%14.6e eps=%14.6e lambda=%14.6e de[%14.6e %14.6e] e[%14.6e %14.6e] %14.6e %14.6e\n",
r, sig, eps/lambda, lambda, dEdROrig, *dEdR, E_Orig, *energy, softcoreLJInv, sig6 );
}
#endif
}
plugins/freeEnergy/platforms/reference/src/SimTKReference/ReferenceFreeEnergyLJCoulombSoftcoreIxn.h
View file @ 9dd055eb
......@@ -42,7 +42,6 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
RealOpenMM periodicBoxSize[3];
RealOpenMM cutoffDistance;
RealOpenMM krf, crf;
RealOpenMM softCoreLJLambda;
int numRx, numRy, numRz;
RealOpenMM alphaEwald;
......@@ -62,15 +61,13 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculateOneIxn( int atom1, int atom2, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
RealOpenMM* totalEnergy ) const;
public:
......@@ -116,41 +113,7 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
/**---------------------------------------------------------------------------------------
Set the force to use Ewald summation.
@param alpha the Ewald separation parameter
@param kmaxx the largest wave vector in the x direction
@param kmaxy the largest wave vector in the y direction
@param kmaxz the largest wave vector in the z direction
--------------------------------------------------------------------------------------- */
void setUseEwald(RealOpenMM alpha, int kmaxx, int kmaxy, int kmaxz);
/**---------------------------------------------------------------------------------------
Set the force to use Particle-Mesh Ewald (PME) summation.
@param alpha the Ewald separation parameter
--------------------------------------------------------------------------------------- */
void setUsePME(RealOpenMM alpha);
/**---------------------------------------------------------------------------------------
Set the soft core LJ lambda
@param lambda the soft core LJ lambda
--------------------------------------------------------------------------------------- */
void setSoftCoreLJLambda(RealOpenMM lambda);
/**---------------------------------------------------------------------------------------
Calculate parameters for LJ 1-4 ixn
Calculate parameters for ixn
@param c6 c6
@param c12 c12
......@@ -164,8 +127,7 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
--------------------------------------------------------------------------------------- */
void getDerivedParameters( RealOpenMM c6, RealOpenMM c12, RealOpenMM q1,
RealOpenMM epsfacSqrt,
RealOpenMM* parameters ) const;
RealOpenMM epsfacSqrt, RealOpenMM* parameters ) const;
/**---------------------------------------------------------------------------------------
......@@ -180,62 +142,16 @@ class ReferenceFreeEnergyLJCoulombSoftcoreIxn : public ReferencePairIxn {
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculatePairIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Ewald ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculateEwaldIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
/**---------------------------------------------------------------------------------------
Calculate PME ixn
@param numberOfAtoms number of atoms
@param atomCoordinates atom coordinates
@param atomParameters atom parameters (charges, c6, c12, ...) atomParameters[atomIndex][paramterIndex]
@param exclusions atom exclusion indices exclusions[atomIndex][atomToExcludeIndex]
exclusions[atomIndex][0] = number of exclusions
exclusions[atomIndex][1-no.] = atom indices of atoms to excluded from
interacting w/ atom atomIndex
@param fixedParameters non atom parameters (not currently used)
@param forces force array (forces added)
@param energyByAtom atom energy
@param totalEnergy total energy
--------------------------------------------------------------------------------------- */
void calculatePMEIxn( int numberOfAtoms, std::vector<OpenMM::RealVec>& atomCoordinates,
RealOpenMM** atomParameters, int** exclusions,
RealOpenMM* fixedParameters, std::vector<OpenMM::RealVec>& forces,
RealOpenMM* energyByAtom, RealOpenMM* totalEnergy ) const;
/**---------------------------------------------------------------------------------------
......@@ -250,7 +166,7 @@ private:
--------------------------------------------------------------------------------------- */
void calculateOneLJIxn( RealOpenMM inverseR, RealOpenMM sig, RealOpenMM eps,
RealOpenMM* dEdR, RealOpenMM* energy ) const;
RealOpenMM* dEdR, RealOpenMM* energy ) const;
/**---------------------------------------------------------------------------------------
......@@ -266,7 +182,7 @@ private:
--------------------------------------------------------------------------------------- */
void calculateOneSoftCoreLJIxn( RealOpenMM r, RealOpenMM sig, RealOpenMM eps,
RealOpenMM lambda, RealOpenMM* dEdR, RealOpenMM* energy ) const;
RealOpenMM lambda, RealOpenMM* dEdR, RealOpenMM* energy ) const;
};
......
plugins/freeEnergy/platforms/reference/tests/TestReferenceLJSoftcoreForce.cpp 0 → 100644
View file @ 9dd055eb
/* -------------------------------------------------------------------------- *
* OpenMM *
* -------------------------------------------------------------------------- *
* 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-2009 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 all the different force terms in the reference implementation of CustomGBForce.
*/
#include "../../../tests/AssertionUtilities.h"
#include "sfmt/SFMT.h"
#include "openmm/Context.h"
#include "openmm/CustomBondForce.h"
#include "openmm/CustomNonbondedForce.h"
#include "openmm/NonbondedForce.h"
#include "openmm/NonbondedSoftcoreForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include <iostream>
#include <vector>
#include <algorithm>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
static const int NoCutoff = 0;
static const int CutoffNonPeriodic = 1;
static const int CutoffPeriodic = 2;
void testNonbondedSoftcore( double lambda1, double lambda2, int nonbondedMethod ){
const int numMolecules = 70;
const int numParticles = numMolecules*2;
const double boxSize = 10.0;
const double reactionFieldDielectric = 80.0;
const double cutoffDistance = 0.4*boxSize;
// Create two systems: one with a NonbondedSoftcoreForce, and one using a CustomNonbondedForce to implement the same interaction.
System standardSystem;
System customSystem;
for (int i = 0; i < numParticles; i++) {
standardSystem.addParticle(1.0);
customSystem.addParticle(1.0);
}
standardSystem.setDefaultPeriodicBoxVectors(Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
customSystem.setDefaultPeriodicBoxVectors( Vec3(boxSize, 0, 0), Vec3(0, boxSize, 0), Vec3(0, 0, boxSize));
NonbondedSoftcoreForce* nonbondedSoftcoreForce = new NonbondedSoftcoreForce();
CustomNonbondedForce* customNonbonded;
CustomBondForce* customBond;
if( nonbondedMethod == NoCutoff ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::NoCutoff );
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customNonbonded->setNonbondedMethod( CustomNonbondedForce::NoCutoff );
customBond = new CustomBondForce("lambda*4*eps*(dem^2-dem)+138.935456*q/r;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
} else {
nonbondedSoftcoreForce->setCutoffDistance( cutoffDistance );
nonbondedSoftcoreForce->setReactionFieldDielectric( reactionFieldDielectric );
if( nonbondedMethod == CutoffNonPeriodic ){
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffNonPeriodic );
} else {
nonbondedSoftcoreForce->setNonbondedMethod( NonbondedSoftcoreForce::CutoffPeriodic );
}
customNonbonded = new CustomNonbondedForce("lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf);"
"q=q1*q2;"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda);"
"sigma=0.5*(sigma1+sigma2);"
"eps=sqrt(eps1*eps2);"
"lambda=min(lambda1,lambda2)");
customBond = new CustomBondForce("withinCutoff*(lambda*4*eps*(dem^2-dem)+138.935456*q*(1.0/r+(krf*r*r)-crf));"
"withinCutoff=step(cutoff-r);"
"dem=1.0/(soft+rsig);"
"rsig=(r/sigma)^6;"
"soft=0.5*(1.0-lambda)");
customNonbonded->setCutoffDistance( cutoffDistance );
if( nonbondedMethod == CutoffNonPeriodic ){
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffNonPeriodic );
} else {
customNonbonded->setNonbondedMethod( CustomNonbondedForce::CutoffPeriodic );
}
double eps2 = (reactionFieldDielectric - 1.0)/(2.0*reactionFieldDielectric+1.0);
double kValue = eps2/(cutoffDistance*cutoffDistance*cutoffDistance);
customNonbonded->addGlobalParameter("krf", kValue );
customBond->addGlobalParameter("krf", kValue );
double cValue = (1.0/cutoffDistance)*(3.0*reactionFieldDielectric)/(2.0*reactionFieldDielectric + 1.0);
customNonbonded->addGlobalParameter("crf", cValue );
customBond->addGlobalParameter("crf", cValue );
customBond->addGlobalParameter("cutoff", cutoffDistance );
}
customNonbonded->addPerParticleParameter("q");
customNonbonded->addPerParticleParameter("sigma");
customNonbonded->addPerParticleParameter("eps");
customNonbonded->addPerParticleParameter("lambda");
customBond->addPerBondParameter("q");
customBond->addPerBondParameter("sigma");
customBond->addPerBondParameter("eps");
customBond->addPerBondParameter("lambda");
vector<Vec3> positions(numParticles);
vector<Vec3> velocities(numParticles);
OpenMM_SFMT::SFMT sfmt;
init_gen_rand(0, sfmt);
vector<double> params(4);
// periodic boundary conditions not possible w/ CustomBond?
int includeExceptions = nonbondedMethod == CutoffPeriodic ? 0 : 1;
for (int i = 0; i < numMolecules; i++) {
if (i < numMolecules/2) {
double charge = 1.0;
nonbondedSoftcoreForce->addParticle( charge, 0.2, 0.5, lambda1);
nonbondedSoftcoreForce->addParticle(-charge, 0.1, 0.5, lambda1);
params[0] = charge;
params[1] = 0.2;
params[2] = 0.5;
params[3] = lambda1;
customNonbonded->addParticle(params);
params[0] = -charge;
params[1] = 0.1;
customNonbonded->addParticle(params);
if( includeExceptions && i && ((i%4) == 0) ){
vector<double> bondParams(4);
nonbondedSoftcoreForce->addException(i-4, i, charge*charge, 0.2, 0.5, false, lambda1);
customNonbonded->addExclusion( i-4,i);
bondParams[0] = charge*charge;
bondParams[1] = 0.2;
bondParams[2] = 0.5;
bondParams[3] = lambda1;
customBond->addBond(i-4,i, bondParams );
}
} else {
double charge = 1.2;
nonbondedSoftcoreForce->addParticle( charge, 0.2, 0.8, lambda2);
nonbondedSoftcoreForce->addParticle(-charge, 0.1, 0.8, lambda2);
params[0] = charge;
params[1] = 0.2;
params[2] = 0.8;
params[3] = lambda2;
customNonbonded->addParticle(params);
params[0] = -charge;
params[1] = 0.1;
customNonbonded->addParticle(params);
if( includeExceptions && i && ((i%4) == 0) ){
vector<double> bondParams(4);
nonbondedSoftcoreForce->addException(i-4, i, charge*charge, 0.2, 0.5, false, lambda1);
customNonbonded->addExclusion( i-4,i);
bondParams[0] = charge*charge;
bondParams[1] = 0.2;
bondParams[2] = 0.5;
bondParams[3] = lambda2;
customBond->addBond(i-4,i, bondParams );
}
}
positions[2*i] = Vec3(boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt), boxSize*genrand_real2(sfmt));
positions[2*i+1] = Vec3(positions[2*i][0]+1.0, positions[2*i][1], positions[2*i][2]);
velocities[2*i] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
velocities[2*i+1] = Vec3(genrand_real2(sfmt), genrand_real2(sfmt), genrand_real2(sfmt));
}
standardSystem.addForce(nonbondedSoftcoreForce);
customSystem.addForce(customNonbonded);
customSystem.addForce(customBond);
VerletIntegrator integrator1(0.01);
VerletIntegrator integrator2(0.01);
Context context1(standardSystem, integrator1, Platform::getPlatformByName( "Reference"));
context1.setPositions(positions);
context1.setVelocities(velocities);
State state1 = context1.getState(State::Forces | State::Energy);
Context context2(customSystem, integrator2, Platform::getPlatformByName( "Reference"));
context2.setPositions(positions);
context2.setVelocities(velocities);
State state2 = context2.getState(State::Forces | State::Energy);
// (void) fprintf( stderr, "%10.1f %10.1f %15.7e %15.7e\n", lambda1, lambda2, state1.getPotentialEnergy(), state2.getPotentialEnergy());
ASSERT_EQUAL_TOL(state1.getPotentialEnergy(), state2.getPotentialEnergy(), 1e-4);
for (int i = 0; i < numParticles; i++) {
ASSERT_EQUAL_VEC(state1.getForces()[i], state2.getForces()[i], 1e-4);
}
}
int main() {
try {
// test various combinations of lambdas and boundary conditions/cutoffs
testNonbondedSoftcore( 1.0, 1.0 , NoCutoff );
testNonbondedSoftcore( 1.0, 0.0 , NoCutoff );
testNonbondedSoftcore( 1.0, 0.5 , NoCutoff );
testNonbondedSoftcore( 0.0, 0.0 , NoCutoff );
testNonbondedSoftcore( 1.0, 1.0 , CutoffNonPeriodic );
testNonbondedSoftcore( 1.0, 0.0 , CutoffNonPeriodic );
testNonbondedSoftcore( 1.0, 0.5 , CutoffNonPeriodic );
testNonbondedSoftcore( 0.0, 0.0 , CutoffNonPeriodic );
testNonbondedSoftcore( 1.0, 1.0 , CutoffPeriodic );
testNonbondedSoftcore( 1.0, 0.0 , CutoffPeriodic );
testNonbondedSoftcore( 1.0, 0.5 , CutoffPeriodic );
testNonbondedSoftcore( 0.0, 0.0 , CutoffPeriodic );
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
cout << "Done" << endl;
return 0;
}
plugins/freeEnergy/platforms/reference/tests/TestReferenceOBCSoftcoreForce.cpp
View file @ 9dd055eb
......@@ -37,13 +37,12 @@
#include "sfmt/SFMT.h"
#include "openmm/Context.h"
#include "ReferencePlatform.h"
#include "openmm/CustomGBForce.h"
#include "openmm/GBSAOBCSoftcoreForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include <iostream>
#include <vector>
#include <algorithm>
......@@ -53,7 +52,7 @@ using namespace std;
const double TOL = 1e-5;
void testSoftcoreOBC( double lambda1, double lambda2 ){
void testOBCSoftcore( double lambda1, double lambda2 ){
const int numMolecules = 70;
const int numParticles = numMolecules*2;
......@@ -168,9 +167,9 @@ int main() {
// test various combinations of lambdas
testSoftcoreOBC( 1.0, 1.0 );
testSoftcoreOBC( 1.0, 0.0 );
testSoftcoreOBC( 1.0, 0.5 );
testOBCSoftcore( 1.0, 1.0 );
testOBCSoftcore( 1.0, 0.0 );
testOBCSoftcore( 1.0, 0.5 );
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
......
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