Skip to content

GitLab

Commit 7492dc48 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Modified Reference platform of AmoebaVdwForce to handle PBC, cutoffs, tapering... 

Modified Reference platform of AmoebaVdwForce to handle PBC, cutoffs, tapering and dispersion correction
parent c6ebf6e2
Expand all Show whitespace changes
Inline Side-by-side
Showing with 1702 additions and 186 deletions
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 7492dc48
...@@ -41,6 +41,7 @@ ...@@ -41,6 +41,7 @@
#include "openmm/internal/ContextImpl.h" #include "openmm/internal/ContextImpl.h"
#include "openmm/AmoebaMultipoleForce.h" #include "openmm/AmoebaMultipoleForce.h"
#include "openmm/internal/AmoebaMultipoleForceImpl.h" #include "openmm/internal/AmoebaMultipoleForceImpl.h"
#include "openmm/internal/AmoebaVdwForceImpl.h"
#include <cmath> #include <cmath>
#ifdef _MSC_VER #ifdef _MSC_VER
...@@ -534,10 +535,14 @@ ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string ...@@ -534,10 +535,14 @@ ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string
useCutoff = 0; useCutoff = 0;
usePBC = 0; usePBC = 0;
cutoff = 1.0e+10; cutoff = 1.0e+10;
neighborList = NULL;
} }
ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() { ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
if( neighborList ){
delete neighborList;
}
} }
void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) { void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
...@@ -561,7 +566,7 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A ...@@ -561,7 +566,7 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
force.getParticleParameters( ii, indexIV, sigma, epsilon, reduction ); force.getParticleParameters( ii, indexIV, sigma, epsilon, reduction );
force.getParticleExclusions( ii, exclusions ); force.getParticleExclusions( ii, exclusions );
for( unsigned int jj = 0; jj < exclusions.size(); jj++ ){ for( unsigned int jj = 0; jj < exclusions.size(); jj++ ){
allExclusions[ii].push_back( exclusions[jj] ); allExclusions[ii].insert( exclusions[jj] );
} }
indexIVs[ii] = indexIV; indexIVs[ii] = indexIV;
...@@ -574,6 +579,9 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A ...@@ -574,6 +579,9 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
useCutoff = (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff); useCutoff = (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff);
usePBC = (force.getNonbondedMethod() == AmoebaVdwForce::CutoffPeriodic); usePBC = (force.getNonbondedMethod() == AmoebaVdwForce::CutoffPeriodic);
cutoff = force.getCutoff(); cutoff = force.getCutoff();
neighborList = useCutoff ? new NeighborList() : NULL;
dispersionCoefficient = force.getUseDispersionCorrection() ? AmoebaVdwForceImpl::calcDispersionCorrection(system, force) : 0.0;
} }
double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) { double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
...@@ -581,17 +589,27 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc ...@@ -581,17 +589,27 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
vector<RealVec>& posData = extractPositions(context); vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context); vector<RealVec>& forceData = extractForces(context);
AmoebaReferenceVdwForce vdwForce( sigmaCombiningRule, epsilonCombiningRule ); AmoebaReferenceVdwForce vdwForce( sigmaCombiningRule, epsilonCombiningRule );
RealOpenMM energy;
if( useCutoff ){ if( useCutoff ){
vdwForce.setCutoff( cutoff ); vdwForce.setCutoff( cutoff );
computeNeighborListVoxelHash( *neighborList, numParticles, posData, allExclusions, extractBoxSize(context), usePBC, cutoff, 0.0);
if( usePBC ){ if( usePBC ){
vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffPeriodic); vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffPeriodic);
RealVec& box = extractBoxSize(context);
double minAllowedSize = 1.999999*cutoff;
if (box[0] < minAllowedSize || box[1] < minAllowedSize || box[2] < minAllowedSize){
throw OpenMMException("The periodic box size has decreased to less than twice the cutoff.");
}
vdwForce.setPeriodicBox(box);
energy = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, sigmas, epsilons, reductions, *neighborList, forceData);
energy += dispersionCoefficient/(box[0]*box[1]*box[2]);
} else { } else {
vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffNonPeriodic); vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffNonPeriodic);
} }
} else { } else {
vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::NoCutoff ); vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::NoCutoff );
energy = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, sigmas, epsilons, reductions, allExclusions, forceData);
} }
RealOpenMM energy = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, sigmas, epsilons, reductions, allExclusions, forceData);
return static_cast<double>(energy); return static_cast<double>(energy);
} }
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 7492dc48
...@@ -370,7 +370,7 @@ public: ...@@ -370,7 +370,7 @@ public:
* Initialize the kernel. * Initialize the kernel.
* *
* @param system the System this kernel will be applied to * @param system the System this kernel will be applied to
* @param force the AmoebaMultipoleForce this kernel will be used for * @param force the AmoebaVdwForce this kernel will be used for
*/ */
void initialize(const System& system, const AmoebaVdwForce& force); void initialize(const System& system, const AmoebaVdwForce& force);
/** /**
...@@ -387,8 +387,9 @@ private: ...@@ -387,8 +387,9 @@ private:
int useCutoff; int useCutoff;
int usePBC; int usePBC;
double cutoff; double cutoff;
double dispersionCoefficient;
std::vector<int> indexIVs; std::vector<int> indexIVs;
std::vector< std::vector<int> > allExclusions; std::vector< std::set<int> > allExclusions;
std::vector<RealOpenMM> sigmas; std::vector<RealOpenMM> sigmas;
std::vector<RealOpenMM> epsilons; std::vector<RealOpenMM> epsilons;
std::vector<RealOpenMM> reductions; std::vector<RealOpenMM> reductions;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @ 7492dc48
...@@ -30,17 +30,21 @@ ...@@ -30,17 +30,21 @@
using std::vector; using std::vector;
using OpenMM::RealVec; using OpenMM::RealVec;
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( ) : _nonbondedMethod(NoCutoff) { AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( ) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
_cutoff = 1.0e+10; setTaperCoefficients( _cutoff );
setSigmaCombiningRule( "ARITHMETIC" ); setSigmaCombiningRule( "ARITHMETIC" );
setEpsilonCombiningRule( "GEOMETRIC" ); setEpsilonCombiningRule( "GEOMETRIC" );
_periodicBoxDimensions = RealVec( 0.0, 0.0, 0.0 );
} }
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule ) : _nonbondedMethod(NoCutoff) {
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule ) : _nonbondedMethod(NoCutoff), _cutoff(1.0e+10), _taperCutoffFactor(0.9) {
setTaperCoefficients( _cutoff );
setSigmaCombiningRule( sigmaCombiningRule ); setSigmaCombiningRule( sigmaCombiningRule );
setEpsilonCombiningRule( epsilonCombiningRule ); setEpsilonCombiningRule( epsilonCombiningRule );
_periodicBoxDimensions = RealVec( 0.0, 0.0, 0.0 );
} }
AmoebaReferenceVdwForce::NonbondedMethod AmoebaReferenceVdwForce::getNonbondedMethod( void ) const { AmoebaReferenceVdwForce::NonbondedMethod AmoebaReferenceVdwForce::getNonbondedMethod( void ) const {
...@@ -51,14 +55,36 @@ void AmoebaReferenceVdwForce::setNonbondedMethod( AmoebaReferenceVdwForce::Nonbo ...@@ -51,14 +55,36 @@ void AmoebaReferenceVdwForce::setNonbondedMethod( AmoebaReferenceVdwForce::Nonbo
_nonbondedMethod = nonbondedMethod; _nonbondedMethod = nonbondedMethod;
} }
void AmoebaReferenceVdwForce::setTaperCoefficients( double cutoff ){
_taperCutoff = cutoff*_taperCutoffFactor;
if( _taperCutoff != cutoff ){
_taperCoefficients[C3] = 10.0/pow(_taperCutoff - cutoff, 3.0);
_taperCoefficients[C4] = 15.0/pow(_taperCutoff - cutoff, 4.0);
_taperCoefficients[C5] = 6.0/pow(_taperCutoff - cutoff, 5.0);
} else {
_taperCoefficients[C3] = 0.0;
_taperCoefficients[C4] = 0.0;
_taperCoefficients[C5] = 0.0;
}
}
void AmoebaReferenceVdwForce::setCutoff( double cutoff ){ void AmoebaReferenceVdwForce::setCutoff( double cutoff ){
_cutoff = cutoff; _cutoff = cutoff;
setTaperCoefficients( _cutoff );
} }
double AmoebaReferenceVdwForce::getCutoff( void ) const { double AmoebaReferenceVdwForce::getCutoff( void ) const {
return _cutoff; return _cutoff;
} }
void AmoebaReferenceVdwForce::setPeriodicBox( const RealVec& box ){
_periodicBoxDimensions = box;
}
RealVec AmoebaReferenceVdwForce::getPeriodicBox( void ) const {
return _periodicBoxDimensions;
}
void AmoebaReferenceVdwForce::setSigmaCombiningRule( const std::string& sigmaCombiningRule ){ void AmoebaReferenceVdwForce::setSigmaCombiningRule( const std::string& sigmaCombiningRule ){
_sigmaCombiningRule = sigmaCombiningRule; _sigmaCombiningRule = sigmaCombiningRule;
...@@ -174,12 +200,18 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma, ...@@ -174,12 +200,18 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma,
RealOpenMM xr = particleIPosition[0] - particleJPosition[0]; RealOpenMM xr = particleIPosition[0] - particleJPosition[0];
RealOpenMM yr = particleIPosition[1] - particleJPosition[1]; RealOpenMM yr = particleIPosition[1] - particleJPosition[1];
RealOpenMM zr = particleIPosition[2] - particleJPosition[2]; RealOpenMM zr = particleIPosition[2] - particleJPosition[2];
RealOpenMM r_ij_2 = xr*xr + yr*yr + zr*zr;
if( _nonbondedMethod == CutoffPeriodic ){
xr -= static_cast<RealOpenMM>((floor(xr/_periodicBoxDimensions[0]+0.5)*_periodicBoxDimensions[0]));
yr -= static_cast<RealOpenMM>((floor(yr/_periodicBoxDimensions[1]+0.5)*_periodicBoxDimensions[1]));
zr -= static_cast<RealOpenMM>((floor(zr/_periodicBoxDimensions[2]+0.5)*_periodicBoxDimensions[2]));
}
RealOpenMM r_ij_2 = xr*xr + yr*yr + zr*zr;
RealOpenMM r_ij = SQRT(r_ij_2);
RealOpenMM sigma_7 = combindedSigma*combindedSigma*combindedSigma; RealOpenMM sigma_7 = combindedSigma*combindedSigma*combindedSigma;
sigma_7 = sigma_7*sigma_7*combindedSigma; sigma_7 = sigma_7*sigma_7*combindedSigma;
RealOpenMM r_ij = SQRT(r_ij_2);
RealOpenMM r_ij_6 = r_ij_2*r_ij_2*r_ij_2; RealOpenMM r_ij_6 = r_ij_2*r_ij_2*r_ij_2;
RealOpenMM r_ij_7 = r_ij_6*r_ij; RealOpenMM r_ij_7 = r_ij_6*r_ij;
...@@ -197,6 +229,17 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma, ...@@ -197,6 +229,17 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma,
RealOpenMM energy = combindedEpsilon*tau_7*sigma_7*( (ghal+one)*sigma_7/rho - two); RealOpenMM energy = combindedEpsilon*tau_7*sigma_7*( (ghal+one)*sigma_7/rho - two);
RealOpenMM dEdR = -seven*(dtau*energy + gtau); RealOpenMM dEdR = -seven*(dtau*energy + gtau);
// tapering
if( (_nonbondedMethod == CutoffNonPeriodic || _nonbondedMethod == CutoffPeriodic) && r_ij > _taperCutoff ){
RealOpenMM delta = r_ij - _taperCutoff;
RealOpenMM taper = 1.0 + delta*delta*delta*(_taperCoefficients[C3] + delta*(_taperCoefficients[C4] + delta*_taperCoefficients[C5]));
RealOpenMM dtaper = delta*delta*(3.0*_taperCoefficients[C3] + delta*(4.0*_taperCoefficients[C4] + delta*5.0*_taperCoefficients[C5]));
dEdR = energy*dtaper + dEdR*taper;
energy *= taper;
}
dEdR /= r_ij; dEdR /= r_ij;
force[0] = dEdR*xr; force[0] = dEdR*xr;
...@@ -207,13 +250,34 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma, ...@@ -207,13 +250,34 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma,
} }
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numParticles, void AmoebaReferenceVdwForce::setReducedPositions( int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& reductions,
std::vector<Vec3>& reducedPositions ) const {
static const RealOpenMM zero = 0.0;
reducedPositions.resize(numParticles);
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
if( reductions[ii] != zero ){
int reductionIndex = indexIVs[ii];
reducedPositions[ii] = Vec3( reductions[ii]*( particlePositions[ii][0] - particlePositions[reductionIndex][0] ) + particlePositions[reductionIndex][0],
reductions[ii]*( particlePositions[ii][1] - particlePositions[reductionIndex][1] ) + particlePositions[reductionIndex][1],
reductions[ii]*( particlePositions[ii][2] - particlePositions[reductionIndex][2] ) + particlePositions[reductionIndex][2] );
} else {
reducedPositions[ii] = Vec3( particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2] );
}
}
}
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
const vector<RealVec>& particlePositions, const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions, const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& allExclusions, const std::vector< std::set<int> >& allExclusions,
vector<RealVec>& forces ) const { vector<RealVec>& forces ) const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -227,20 +291,16 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numPart ...@@ -227,20 +291,16 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numPart
// set reduced coordinates // set reduced coordinates
std::vector<Vec3> reducedPositions; std::vector<Vec3> reducedPositions;
reducedPositions.resize(numParticles); setReducedPositions( numParticles, particlePositions, indexIVs, reductions, reducedPositions );
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
if( reductions[ii] != zero ){ // loop over all particle pairs
int reductionIndex = indexIVs[ii];
reducedPositions[ii] = Vec3( reductions[ii]*( particlePositions[ii][0] - particlePositions[reductionIndex][0] ) + particlePositions[reductionIndex][0],
reductions[ii]*( particlePositions[ii][1] - particlePositions[reductionIndex][1] ) + particlePositions[reductionIndex][1],
reductions[ii]*( particlePositions[ii][2] - particlePositions[reductionIndex][2] ) + particlePositions[reductionIndex][2] );
} else {
reducedPositions[ii] = Vec3( particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2] );
}
}
// loop over all ixns
// (1) initialize exclusion vector // (1) initialize exclusion vector
// (2) calculate pair ixn, if not excluded
// (3) accumulate forces: if particle is a site where interaction position != particle position,
// then call addReducedForce() to apportion force to particle and its covalent partner
// based on reduction factor
// (4) reset exclusion vector
RealOpenMM energy = zero; RealOpenMM energy = zero;
std::vector<unsigned int> exclusions(numParticles, 0); std::vector<unsigned int> exclusions(numParticles, 0);
...@@ -248,8 +308,8 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numPart ...@@ -248,8 +308,8 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numPart
RealOpenMM sigmaI = sigmas[ii]; RealOpenMM sigmaI = sigmas[ii];
RealOpenMM epsilonI = epsilons[ii]; RealOpenMM epsilonI = epsilons[ii];
for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){ for( std::set<int>::iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++ ){
exclusions[allExclusions[ii][jj]] = 1; exclusions[*jj] = 1;
} }
for( unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++ ){ for( unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++ ){
...@@ -281,21 +341,21 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numPart ...@@ -281,21 +341,21 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numPart
} }
} }
for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){ for( std::set<int>::iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++ ){
exclusions[allExclusions[ii][jj]] = 0; exclusions[*jj] = 0;
} }
} }
return energy; return energy;
} }
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyApplyCutoff( int numParticles, RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
const vector<RealVec>& particlePositions, const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions, const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& allExclusions, const NeighborList& neighborList,
vector<RealVec>& forces ) const { vector<RealVec>& forces ) const {
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
...@@ -309,96 +369,44 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyApplyCutoff( int numP ...@@ -309,96 +369,44 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyApplyCutoff( int numP
// set reduced coordinates // set reduced coordinates
std::vector<Vec3> reducedPositions; std::vector<Vec3> reducedPositions;
reducedPositions.resize(numParticles); setReducedPositions( numParticles, particlePositions, indexIVs, reductions, reducedPositions );
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
if( reductions[ii] != zero ){
int reductionIndex = indexIVs[ii];
reducedPositions[ii] = Vec3( reductions[ii]*( particlePositions[ii][0] - particlePositions[reductionIndex][0] ) + particlePositions[reductionIndex][0],
reductions[ii]*( particlePositions[ii][1] - particlePositions[reductionIndex][1] ) + particlePositions[reductionIndex][1],
reductions[ii]*( particlePositions[ii][2] - particlePositions[reductionIndex][2] ) + particlePositions[reductionIndex][2] );
} else {
reducedPositions[ii] = Vec3( particlePositions[ii][0], particlePositions[ii][1], particlePositions[ii][2] );
}
}
// loop over all ixns // loop over neighbor list
// (1) initialize exclusion vector // (1) calculate pair vdw ixn
// (2) accumulate forces: if particle is a site where interaction position != particle position,
// then call addReducedForce() to apportion force to particle and its covalent partner
// based on reduction factor
RealOpenMM energy = zero; RealOpenMM energy = zero;
std::vector<unsigned int> exclusions(numParticles, 0); for( unsigned int ii = 0; ii < neighborList.size(); ii++ ){
for( unsigned int ii = 0; ii < static_cast<unsigned int>(numParticles); ii++ ){
RealOpenMM sigmaI = sigmas[ii]; OpenMM::AtomPair pair = neighborList[ii];
RealOpenMM epsilonI = epsilons[ii]; int siteI = pair.first;
for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){ int siteJ = pair.second;
exclusions[allExclusions[ii][jj]] = 1;
}
for( unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++ ){ RealOpenMM combindedSigma = (this->*_combineSigmas)( sigmas[siteI], sigmas[siteJ] );
if( exclusions[jj] == 0 ){ RealOpenMM combindedEpsilon = (this->*_combineEpsilons)( epsilons[siteI], epsilons[siteJ] );
RealOpenMM combindedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj] );
RealOpenMM combindedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj] );
Vec3 force; Vec3 force;
energy += calculatePairIxn( combindedSigma, combindedEpsilon, energy += calculatePairIxn( combindedSigma, combindedEpsilon,
reducedPositions[ii], reducedPositions[jj], reducedPositions[siteI], reducedPositions[siteJ], force );
force );
if( indexIVs[ii] == ii ){ if( indexIVs[siteI] == siteI ){
forces[ii][0] -= force[0]; forces[siteI][0] -= force[0];
forces[ii][1] -= force[1]; forces[siteI][1] -= force[1];
forces[ii][2] -= force[2]; forces[siteI][2] -= force[2];
} else { } else {
addReducedForce( ii, indexIVs[ii], reductions[ii], -one, force, forces ); addReducedForce( siteI, indexIVs[siteI], reductions[siteI], -one, force, forces );
} }
if( indexIVs[jj] == jj ){ if( indexIVs[siteJ] == siteJ ){
forces[jj][0] += force[0]; forces[siteJ][0] += force[0];
forces[jj][1] += force[1]; forces[siteJ][1] += force[1];
forces[jj][2] += force[2]; forces[siteJ][2] += force[2];
} else { } else {
addReducedForce( jj, indexIVs[jj], reductions[jj], one, force, forces ); addReducedForce( siteJ, indexIVs[siteJ], reductions[siteJ], one, force, forces );
}
}
} }
for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){
exclusions[allExclusions[ii][jj]] = 0;
}
} }
return energy; return energy;
} }
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
const vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& allExclusions,
vector<RealVec>& forces ) const {
if( getNonbondedMethod() == NoCutoff ){
return calculateForceAndEnergyNoCutoff( numParticles, particlePositions,
indexIVs, sigmas, epsilons,
reductions, allExclusions, forces );
} else {
return calculateForceAndEnergyApplyCutoff( numParticles, particlePositions,
indexIVs, sigmas, epsilons,
reductions, allExclusions, forces );
}
}
/*
double cutoff = force.getCutoff();
double taperCutoff = cutoff*0.9;
replacements["CUTOFF_DISTANCE"] = cu.doubleToString(force.getCutoff());
replacements["TAPER_CUTOFF"] = cu.doubleToString(taperCutoff);
replacements["TAPER_C3"] = cu.doubleToString(10/pow(taperCutoff-cutoff, 3.0));
replacements["TAPER_C4"] = cu.doubleToString(15/pow(taperCutoff-cutoff, 4.0));
replacements["TAPER_C5"] = cu.doubleToString(6/pow(taperCutoff-cutoff, 5.0));
*/
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h
View file @ 7492dc48
...@@ -27,6 +27,7 @@ ...@@ -27,6 +27,7 @@
#include "SimTKUtilities/RealVec.h" #include "SimTKUtilities/RealVec.h"
#include "openmm/Vec3.h" #include "openmm/Vec3.h"
#include "SimTKReference/ReferenceNeighborList.h"
#include <string> #include <string>
#include <vector> #include <vector>
...@@ -169,6 +170,26 @@ public: ...@@ -169,6 +170,26 @@ public:
std::string getEpsilonCombiningRule( void ) const; std::string getEpsilonCombiningRule( void ) const;
/**---------------------------------------------------------------------------------------
Set box dimensions
@param box box dimensions
--------------------------------------------------------------------------------------- */
void setPeriodicBox( const RealVec& box );
/**---------------------------------------------------------------------------------------
Get box dimensions
@return box dimensions
--------------------------------------------------------------------------------------- */
RealVec getPeriodicBox( void ) const;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
Calculate Amoeba Hal vdw ixns Calculate Amoeba Hal vdw ixns
...@@ -190,16 +211,49 @@ public: ...@@ -190,16 +211,49 @@ public:
const std::vector<int>& indexIVs, const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons, const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions, const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& vdwExclusions, const std::vector< std::set<int> >& vdwExclusions,
std::vector<OpenMM::RealVec>& forces ) const;
/**---------------------------------------------------------------------------------------
Calculate Vdw ixn using neighbor list
@param numParticles number of particles
@param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle
@param sigmas particle sigmas
@param epsilons particle epsilons
@param reductions particle reduction factors
@param neighborList neighbor list
@param forces add forces to this vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const NeighborList& neighborList,
std::vector<OpenMM::RealVec>& forces ) const; std::vector<OpenMM::RealVec>& forces ) const;
private: private:
// taper coefficient indices
static const int C3=0;
static const int C4=1;
static const int C5=2;
std::string _sigmaCombiningRule; std::string _sigmaCombiningRule;
std::string _epsilonCombiningRule; std::string _epsilonCombiningRule;
NonbondedMethod _nonbondedMethod; NonbondedMethod _nonbondedMethod;
double _cutoff; double _cutoff;
double _taperCutoffFactor;
double _taperCutoff;
RealOpenMM _taperCoefficients[3];
RealVec _periodicBoxDimensions;
CombiningFunction _combineSigmas; CombiningFunction _combineSigmas;
RealOpenMM arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const; RealOpenMM arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
RealOpenMM geometricSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const; RealOpenMM geometricSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
...@@ -211,6 +265,27 @@ private: ...@@ -211,6 +265,27 @@ private:
RealOpenMM harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const; RealOpenMM harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const; RealOpenMM hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
/**---------------------------------------------------------------------------------------
Set reduced positions: position used to calculate vdw interaction is moved towards
covalent partner
@param numParticles number of particles
@param particlePositions current particle positions
@param indexIVs particle index of covalent partner
@param reductions fraction of bond length to move particle interacting site;
reductions[i] = zero,
if interacting position == particle position
@param reducedPositions output: modfied or original position depending on whether
reduction factor is nonzero
--------------------------------------------------------------------------------------- */
void setReducedPositions( int numParticles, const std::vector<RealVec>& particlePositions,
const std::vector<int>& indexIVs, const std::vector<RealOpenMM>& reductions,
std::vector<Vec3>& reducedPositions ) const;
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
Add reduced forces to force vector Add reduced forces to force vector
...@@ -228,6 +303,16 @@ private: ...@@ -228,6 +303,16 @@ private:
RealOpenMM reduction, RealOpenMM sign, RealOpenMM reduction, RealOpenMM sign,
Vec3& force, std::vector<OpenMM::RealVec>& forces ) const; Vec3& force, std::vector<OpenMM::RealVec>& forces ) const;
/**---------------------------------------------------------------------------------------
Set taper coefficients
@param cutoff cutoff
--------------------------------------------------------------------------------------- */
void setTaperCoefficients( double cutoff );
/**--------------------------------------------------------------------------------------- /**---------------------------------------------------------------------------------------
Calculate pair ixn Calculate pair ixn
...@@ -246,54 +331,6 @@ private: ...@@ -246,54 +331,6 @@ private:
const Vec3& particleIPosition, const Vec3& particleJPosition, const Vec3& particleIPosition, const Vec3& particleJPosition,
Vec3& force ) const; Vec3& force ) const;
/**---------------------------------------------------------------------------------------
Calculate Vdw ixns w/ no cutoff
@param numParticles number of particles
@param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle
@param sigmas particle sigmas
@param epsilons particle epsilons
@param reductions particle reduction factors
@param vdwExclusions particle exclusions
@param forces add forces to this vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergyNoCutoff( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& vdwExclusions,
std::vector<OpenMM::RealVec>& forces ) const;
/**---------------------------------------------------------------------------------------
Calculate Vdw ixns w/ cutoff
@param numParticles number of particles
@param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle
@param sigmas particle sigmas
@param epsilons particle epsilons
@param reductions particle reduction factors
@param vdwExclusions particle exclusions
@param forces add forces to this vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergyApplyCutoff( int numParticles, const std::vector<OpenMM::RealVec>& particlePositions,
const std::vector<int>& indexIVs,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& vdwExclusions,
std::vector<OpenMM::RealVec>& forces ) const;
}; };
// --------------------------------------------------------------------------------------- // ---------------------------------------------------------------------------------------
......
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