Skip to content

GitLab

Commit 1f2cbefe authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Reference VdwForce and Test

parent c7387612
Show whitespace changes
Inline Side-by-side
Showing with 837 additions and 95 deletions
plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp
View file @ 1f2cbefe
......@@ -93,10 +93,12 @@ void AmoebaVdwForce::setParticleExclusions( int particleIndex, std::vector< int
void AmoebaVdwForce::getParticleExclusions( int particleIndex, std::vector< int >& outputExclusions ) const {
if( particleIndex < exclusions.size() ){
outputExclusions.resize( exclusions[particleIndex].size() );
for( unsigned int ii = 0; ii < exclusions[particleIndex].size(); ii++ ){
outputExclusions[ii] = exclusions[particleIndex][ii];
}
}
}
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp
View file @ 1f2cbefe
......@@ -46,6 +46,7 @@ extern "C" void OPENMM_AMOEBA_REFERENCE_EXPORT registerKernelFactories() {
platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(), factory);
/*
platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory);
......@@ -84,6 +85,9 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
if (name == CalcAmoebaTorsionTorsionForceKernel::Name())
return new ReferenceCalcAmoebaTorsionTorsionForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaVdwForceKernel::Name())
return new ReferenceCalcAmoebaVdwForceKernel(name, platform, context.getSystem());
/*
if (name == CalcAmoebaMultipoleForceKernel::Name())
......@@ -92,9 +96,6 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
if (name == CalcAmoebaGeneralizedKirkwoodForceKernel::Name())
return new ReferenceCalcAmoebaGeneralizedKirkwoodForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaVdwForceKernel::Name())
return new ReferenceCalcAmoebaVdwForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaWcaDispersionForceKernel::Name())
return new ReferenceCalcAmoebaWcaDispersionForceKernel(name, platform, context.getSystem());
*/
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 1f2cbefe
......@@ -33,6 +33,7 @@
#include "AmoebaReferenceStretchBendForce.h"
#include "AmoebaReferenceOutOfPlaneBendForce.h"
#include "AmoebaReferenceTorsionTorsionForce.h"
#include "AmoebaReferenceVdwForce.h"
#include "ReferencePlatform.h"
#include "openmm/internal/ContextImpl.h"
//#include "internal/AmoebaMultipoleForceImpl.h"
......@@ -731,67 +732,58 @@ double ReferenceCalcAmoebaTorsionTorsionForceKernel::execute(ContextImpl& contex
// // handled in computeAmoebaMultipoleForce()
// return 0.0;
//}
//
//static void computeAmoebaVdwForce( AmoebaReferenceData& data ) {
//
// amoebaGpuContext gpu = data.getAmoebaGpu();
// data.initializeGpu();
//
// // Vdw14_7F
//
// kCalculateAmoebaVdw14_7Forces(gpu);
//}
//
//ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system) :
// CalcAmoebaVdwForceKernel(name, platform), system(system) {
// data.incrementKernelCount();
//}
//
//ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
// data.decrementKernelCount();
//}
//
//void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
//
// // per-particle parameters
//
// int numParticles = system.getNumParticles();
//
// std::vector<int> indexIVs(numParticles);
// std::vector<int> indexClasses(numParticles);
// std::vector< std::vector<int> > allExclusions(numParticles);
// std::vector<RealOpenMM> sigmas(numParticles);
// std::vector<RealOpenMM> epsilons(numParticles);
// std::vector<RealOpenMM> reductions(numParticles);
// for( int ii = 0; ii < numParticles; ii++ ){
//
// int indexIV, indexClass;
// double sigma, epsilon, reduction;
// std::vector<int> exclusions;
//
// force.getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
// force.getParticleExclusions( ii, exclusions );
// for( unsigned int jj = 0; jj < exclusions.size(); jj++ ){
// allExclusions[ii].push_back( exclusions[jj] );
// }
//
// indexIVs[ii] = indexIV;
// indexClasses[ii] = indexClass;
// sigmas[ii] = static_cast<RealOpenMM>( sigma );
// epsilons[ii] = static_cast<RealOpenMM>( epsilon );
// reductions[ii] = static_cast<RealOpenMM>( reduction );
// }
//
// gpuSetAmoebaVdwParameters( data.getAmoebaGpu(), indexIVs, indexClasses, sigmas, epsilons, reductions,
// force.getSigmaCombiningRule(), force.getEpsilonCombiningRule(),
// allExclusions );
//}
//
//double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
// computeAmoebaVdwForce( data );
// return 0.0;
//}
//
ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaVdwForceKernel(name, platform), system(system) {
}
ReferenceCalcAmoebaVdwForceKernel::~ReferenceCalcAmoebaVdwForceKernel() {
}
void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const AmoebaVdwForce& force) {
// per-particle parameters
numParticles = system.getNumParticles();
indexIVs.resize( numParticles );
indexClasses.resize( numParticles );
allExclusions.resize( numParticles );
sigmas.resize( numParticles );
epsilons.resize( numParticles );
reductions.resize( numParticles );
for( int ii = 0; ii < numParticles; ii++ ){
int indexIV, indexClass;
double sigma, epsilon, reduction;
std::vector<int> exclusions;
force.getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
force.getParticleExclusions( ii, exclusions );
for( unsigned int jj = 0; jj < exclusions.size(); jj++ ){
allExclusions[ii].push_back( exclusions[jj] );
}
indexIVs[ii] = indexIV;
indexClasses[ii] = indexClass;
sigmas[ii] = static_cast<RealOpenMM>( sigma );
epsilons[ii] = static_cast<RealOpenMM>( epsilon );
reductions[ii] = static_cast<RealOpenMM>( reduction );
}
sigmaCombiningRule = force.getSigmaCombiningRule();
epsilonCombiningRule = force.getEpsilonCombiningRule();
}
double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
RealOpenMM** posData = extractPositions(context);
RealOpenMM** forceData = extractForces(context);
AmoebaReferenceVdwForce vdwForce( sigmaCombiningRule, epsilonCombiningRule, AmoebaReferenceVdwForce::NoCutoff );
RealOpenMM energy = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, indexClasses, sigmas, epsilons, reductions, allExclusions, forceData);
return static_cast<double>(energy);
}
///* -------------------------------------------------------------------------- *
// * AmoebaWcaDispersion *
// * -------------------------------------------------------------------------- */
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 1f2cbefe
......@@ -379,34 +379,43 @@ private:
// private:
// System& system;
// };
//
// /**
// * This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
// */
// class ReferenceCalcAmoebaVdwForceKernel : public CalcAmoebaVdwForceKernel {
// public:
// ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system);
// ~ReferenceCalcAmoebaVdwForceKernel();
// /**
// * Initialize the kernel.
// *
// * @param system the System this kernel will be applied to
// * @param force the AmoebaMultipoleForce this kernel will be used for
// */
// void initialize(const System& system, const AmoebaVdwForce& force);
// /**
// * Execute the kernel to calculate the forces and/or energy.
// *
// * @param context the context in which to execute this kernel
// * @param includeForces true if forces should be calculated
// * @param includeEnergy true if the energy should be calculated
// * @return the potential energy due to the force
// */
// double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
// private:
// System& system;
// };
//
/**
* This kernel is invoked to calculate the vdw forces acting on the system and the energy of the system.
*/
class ReferenceCalcAmoebaVdwForceKernel : public CalcAmoebaVdwForceKernel {
public:
ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system);
~ReferenceCalcAmoebaVdwForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the AmoebaMultipoleForce this kernel will be used for
*/
void initialize(const System& system, const AmoebaVdwForce& force);
/**
* Execute the kernel to calculate the forces and/or energy.
*
* @param context the context in which to execute this kernel
* @param includeForces true if forces should be calculated
* @param includeEnergy true if the energy should be calculated
* @return the potential energy due to the force
*/
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
int numParticles;
std::vector<int> indexIVs;
std::vector<int> indexClasses;
std::vector< std::vector<int> > allExclusions;
std::vector<RealOpenMM> sigmas;
std::vector<RealOpenMM> epsilons;
std::vector<RealOpenMM> reductions;
std::string sigmaCombiningRule;
std::string epsilonCombiningRule;
System& system;
};
// /**
// * This kernel is invoked to calculate the WCA dispersion forces acting on the system and the energy of the system.
// */
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp 0 → 100644
View file @ 1f2cbefe
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* 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.
*/
#include "AmoebaReferenceForce.h"
#include "AmoebaReferenceVdwForce.h"
#include <vector>
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( ) : _nonbondedMethod(NoCutoff) {
setSigmaCombiningRule( "ARITHMETIC" );
setEpsilonCombiningRule( "GEOMETRIC" );
}
AmoebaReferenceVdwForce::AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, NonbondedMethod nonbondedMethod ) :
_nonbondedMethod(nonbondedMethod) {
setSigmaCombiningRule( sigmaCombiningRule );
setEpsilonCombiningRule( epsilonCombiningRule );
}
AmoebaReferenceVdwForce::NonbondedMethod AmoebaReferenceVdwForce::getNonbondedMethod( void ) const {
return _nonbondedMethod;
}
void AmoebaReferenceVdwForce::setNonbondedMethod( AmoebaReferenceVdwForce::NonbondedMethod nonbondedMethod ){
_nonbondedMethod = nonbondedMethod;
}
void AmoebaReferenceVdwForce::setSigmaCombiningRule( const std::string& sigmaCombiningRule ){
_sigmaCombiningRule = sigmaCombiningRule;
// convert to upper case and set combining function
std::transform( _sigmaCombiningRule.begin(), _sigmaCombiningRule.end(), _sigmaCombiningRule.begin(), (int(*)(int)) std::toupper);
if( _sigmaCombiningRule == "GEOMETRIC" ){
_combineSigmas = &AmoebaReferenceVdwForce::geometricSigmaCombiningRule;
} else if( _sigmaCombiningRule == "CUBIC-MEAN" ){
_combineSigmas = &AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule;
} else {
_combineSigmas = &AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule;
}
}
std::string AmoebaReferenceVdwForce::getSigmaCombiningRule( void ) const {
return _sigmaCombiningRule;
}
RealOpenMM AmoebaReferenceVdwForce::arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const {
return (sigmaI + sigmaJ);
}
RealOpenMM AmoebaReferenceVdwForce::geometricSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const {
return 2.0*SQRT(sigmaI*sigmaJ);
}
RealOpenMM AmoebaReferenceVdwForce::cubicMeanSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const {
const RealOpenMM zero = 0.0;
RealOpenMM sigmaI2 = sigmaI*sigmaI;
RealOpenMM sigmaJ2 = sigmaJ*sigmaJ;
return sigmaI != zero && sigmaJ != 0.0 ? 2.0*(sigmaI2*sigmaI + sigmaJ2*sigmaJ)/(sigmaI2 + sigmaJ2) : zero;
}
void AmoebaReferenceVdwForce::setEpsilonCombiningRule( const std::string& epsilonCombiningRule ){
_epsilonCombiningRule = epsilonCombiningRule;
std::transform( _epsilonCombiningRule.begin(), _epsilonCombiningRule.end(), _epsilonCombiningRule.begin(), (int(*)(int)) std::toupper);
// convert to upper case and set combining function
if( _epsilonCombiningRule == "ARITHMETIC" ){
_combineEpsilons = &AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule;
} else if( _epsilonCombiningRule == "HARMONIC" ){
_combineEpsilons = &AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule;
} else if( _epsilonCombiningRule == "HHG" ){
_combineEpsilons = &AmoebaReferenceVdwForce::hhgEpsilonCombiningRule;
} else {
_combineEpsilons = &AmoebaReferenceVdwForce::geometricEpsilonCombiningRule;
}
}
std::string AmoebaReferenceVdwForce::getEpsilonCombiningRule( void ) const {
return _epsilonCombiningRule;
}
RealOpenMM AmoebaReferenceVdwForce::arithmeticEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
return 0.5*(epsilonI + epsilonJ);
}
RealOpenMM AmoebaReferenceVdwForce::geometricEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
return 2.0*SQRT(epsilonI*epsilonJ);
}
RealOpenMM AmoebaReferenceVdwForce::harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
return (epsilonI != 0.0 && epsilonJ != 0.0) ? 2.0*(epsilonI*epsilonJ)/(epsilonI + epsilonJ) : 0.0;
}
RealOpenMM AmoebaReferenceVdwForce::hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const {
RealOpenMM denominator = SQRT(epsilonI) + SQRT(epsilonJ);
return (epsilonI != 0.0 && epsilonJ != 0.0) ? 4.0*(epsilonI*epsilonJ)/(denominator*denominator) : 0.0;
}
void AmoebaReferenceVdwForce::addReducedForce( unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign,
Vec3& force, RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = 1.0;
// ---------------------------------------------------------------------------------------
forces[particleI][0] += sign*force[0]*reduction;
forces[particleI][1] += sign*force[1]*reduction;
forces[particleI][2] += sign*force[2]*reduction;
forces[particleIV][0] += sign*force[0]*(one - reduction);
forces[particleIV][1] += sign*force[1]*(one - reduction);
forces[particleIV][2] += sign*force[2]*(one - reduction);
}
RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma, RealOpenMM combindedEpsilon,
const Vec3& particleIPosition,
const Vec3& particleJPosition,
Vec3& force ) const {
// ---------------------------------------------------------------------------------------
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM seven = 7.0;
static const RealOpenMM dhal = 0.07;
static const RealOpenMM ghal = 0.12;
// ---------------------------------------------------------------------------------------
RealOpenMM xr = particleIPosition[0] - particleJPosition[0];
RealOpenMM yr = particleIPosition[1] - particleJPosition[1];
RealOpenMM zr = particleIPosition[2] - particleJPosition[2];
RealOpenMM r_ij_2 = xr*xr + yr*yr + zr*zr;
RealOpenMM sigma_7 = combindedSigma*combindedSigma*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_7 = r_ij_6*r_ij;
RealOpenMM rho = r_ij_7 + ghal*sigma_7;
RealOpenMM tau = (dhal + one)/(r_ij + dhal*combindedSigma);
RealOpenMM tau_7 = tau*tau*tau;
tau_7 = tau_7*tau_7*tau;
RealOpenMM dtau = tau/(dhal + one);
RealOpenMM ratio = (sigma_7/rho);
RealOpenMM gtau = combindedEpsilon*tau_7*r_ij_6*(ghal+one)*ratio*ratio;
RealOpenMM energy = combindedEpsilon*tau_7*sigma_7*( (ghal+one)*sigma_7/rho - two);
RealOpenMM dEdR = -seven*(dtau*energy + gtau);
dEdR /= r_ij;
force[0] = dEdR*xr;
force[1] = dEdR*yr;
force[2] = dEdR*zr;
return energy;
}
RealOpenMM AmoebaReferenceVdwForce::calculateNoCutoffForceAndEnergy( int numParticles,
RealOpenMM** const particlePositions,
const std::vector<int>& indexIVs,
const std::vector<int>& indexClasses,
const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& allExclusions,
RealOpenMM** forces ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceVdwForce::calculateNoCutoffForceAndEnergy";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
// set reduced coordinates
std::vector<Vec3> reducedPositions;
reducedPositions.resize(numParticles);
for( unsigned int ii = 0; ii < 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
// (1) initialize exclusion vector
RealOpenMM energy = zero;
std::vector<unsigned int> exclusions(numParticles, 0);
for( unsigned int ii = 0; ii < numParticles; ii++ ){
RealOpenMM sigmaI = sigmas[ii];
RealOpenMM epsilonI = epsilons[ii];
for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){
exclusions[allExclusions[ii][jj]] = 1;
}
for( unsigned int jj = ii+1; jj < numParticles; jj++ ){
if( exclusions[jj] == 0 ){
RealOpenMM combindedSigma = (this->*_combineSigmas)(sigmaI, sigmas[jj] );
RealOpenMM combindedEpsilon = (this->*_combineEpsilons)(epsilonI, epsilons[jj] );
Vec3 force;
energy += calculatePairIxn( combindedSigma, combindedEpsilon,
reducedPositions[ii], reducedPositions[jj],
force );
if( indexIVs[ii] == ii ){
forces[ii][0] -= force[0];
forces[ii][1] -= force[1];
forces[ii][2] -= force[2];
} else {
addReducedForce( ii, indexIVs[ii], reductions[ii], -one, force, forces );
}
if( indexIVs[jj] == jj ){
forces[jj][0] += force[0];
forces[jj][1] += force[1];
forces[jj][2] += force[2];
} else {
addReducedForce( jj, indexIVs[jj], reductions[jj], one, force, forces );
}
}
}
for( unsigned int jj = 0; jj < allExclusions[ii].size(); jj++ ){
exclusions[allExclusions[ii][jj]] = 0;
}
}
return energy;
}
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
RealOpenMM** const particlePositions,
const std::vector<int>& indexIVs,
const std::vector<int>& indexClasses,
const std::vector<RealOpenMM>& sigmas,
const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& allExclusions,
RealOpenMM** forces ) const {
if( getNonbondedMethod() == NoCutoff || 1 ){
return calculateNoCutoffForceAndEnergy( numParticles, particlePositions,
indexIVs, indexClasses, sigmas, epsilons,
reductions, allExclusions, forces );
}
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.h 0 → 100644
View file @ 1f2cbefe
/* Portions copyright (c) 2006 Stanford University and Simbios.
* Contributors: Pande Group
*
* 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.
*/
#ifndef __AmoebaReferenceVdwForce_H__
#define __AmoebaReferenceVdwForce_H__
#include "SimTKUtilities/SimTKOpenMMRealType.h"
#include "openmm/Vec3.h"
#include <string>
using namespace OpenMM;
class AmoebaReferenceVdwForce;
typedef RealOpenMM (AmoebaReferenceVdwForce::*CombiningFunction)( RealOpenMM x, RealOpenMM y) const;
// ---------------------------------------------------------------------------------------
class AmoebaReferenceVdwForce {
public:
/**
* This is an enumeration of the different methods that may be used for handling long range Vdw forces.
*/
enum NonbondedMethod {
/**
* No cutoff is applied to the interactions. The full set of N^2 interactions is computed exactly.
* This necessarily means that periodic boundary conditions cannot be used. This is the default.
*/
NoCutoff = 0,
/**
* Interactions beyond the cutoff distance are ignored.
*/
CutoffNonPeriodic = 1,
/**
* Periodic boundary conditions are used, so that each particle interacts only with the nearest periodic copy of
* each other particle. Interactions beyond the cutoff distance are ignored.
*/
CutoffPeriodic = 2,
};
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceVdwForce( );
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceVdwForce( const std::string& sigmaCombiningRule,
const std::string& epsilonCombiningRule,
NonbondedMethod nonbondedMethod );
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~AmoebaReferenceVdwForce( ){};
/**---------------------------------------------------------------------------------------
Get nonbonded method
@return nonbonded method
--------------------------------------------------------------------------------------- */
NonbondedMethod getNonbondedMethod( void ) const;
/**---------------------------------------------------------------------------------------
Set nonbonded method
@param nonbonded method
--------------------------------------------------------------------------------------- */
void setNonbondedMethod( NonbondedMethod nonbondedMethod );
/**---------------------------------------------------------------------------------------
Set sigma combining rule
@param sigmaCombiningRule rule: GEOMETRIC, CUBIC-MEAN, ARITHMETIC (default)
--------------------------------------------------------------------------------------- */
void setSigmaCombiningRule( const std::string& sigmaCombiningRule );
/**---------------------------------------------------------------------------------------
Get sigma combining rule
@return sigmaCombiningRule
--------------------------------------------------------------------------------------- */
std::string getSigmaCombiningRule( void ) const;
/**---------------------------------------------------------------------------------------
Set epsilon combining rule
@param epsilonCombiningRule rule: GEOMETRIC, CUBIC-MEAN, ARITHMETIC (default)
--------------------------------------------------------------------------------------- */
void setEpsilonCombiningRule( const std::string& epsilonCombiningRule );
/**---------------------------------------------------------------------------------------
Get epsilon combining rule
@return epsilonCombiningRule
--------------------------------------------------------------------------------------- */
std::string getEpsilonCombiningRule( void ) const;
/**---------------------------------------------------------------------------------------
Calculate Amoeba Hal vdw ixns
@param numParticles number of particles
@param particlePositions Cartesian coordinates of particles
@param indexIVs position index for associated reducing particle
@param indexClasses class index for combining sigmas/epsilons (not currently used)
@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
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numParticles, RealOpenMM** const particlePositions,
const std::vector<int>& indexIVs,
const std::vector<int>& indexClasses,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& vdwExclusions,
RealOpenMM** forces ) const;
private:
std::string _sigmaCombiningRule;
std::string _epsilonCombiningRule;
NonbondedMethod _nonbondedMethod;
CombiningFunction _combineSigmas;
RealOpenMM arithmeticSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
RealOpenMM geometricSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
RealOpenMM cubicMeanSigmaCombiningRule( RealOpenMM sigmaI, RealOpenMM sigmaJ ) const;
CombiningFunction _combineEpsilons;
RealOpenMM arithmeticEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM geometricEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM harmonicEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
RealOpenMM hhgEpsilonCombiningRule( RealOpenMM epsilonI, RealOpenMM epsilonJ ) const;
/**---------------------------------------------------------------------------------------
Add reduced forces to force vector
@param particleI index of particleI
@param particleIV index of particleIV
@param reduction reduction factor
@param sign +1 or -1 for add/sutracting forces
@param force force vector to add
@param forces force vector for particles
--------------------------------------------------------------------------------------- */
void addReducedForce( unsigned int particleI, unsigned int particleIV,
RealOpenMM reduction, RealOpenMM sign,
Vec3& force, RealOpenMM** forces ) const;
/**---------------------------------------------------------------------------------------
Calculate pair ixn
@param combindedSigma combined sigmas
@param combindedEpsilon combined epsilons
@param particleIPosition particle I position
@param particleJPosition particle J position
@param force output force
@return energy for ixn
--------------------------------------------------------------------------------------- */
RealOpenMM calculatePairIxn( RealOpenMM combindedSigma, RealOpenMM combindedEpsilon,
const Vec3& particleIPosition, const Vec3& particleJPosition,
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 indexClasses class index for combining sigmas/epsilons (not currently used)
@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 calculateNoCutoffForceAndEnergy( int numParticles, RealOpenMM** const particlePositions,
const std::vector<int>& indexIVs,
const std::vector<int>& indexClasses,
const std::vector<RealOpenMM>& sigmas, const std::vector<RealOpenMM>& epsilons,
const std::vector<RealOpenMM>& reductions,
const std::vector< std::vector<int> >& vdwExclusions,
RealOpenMM** forces ) const;
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceVdwForce___
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaVdwForce.cpp 0 → 100644
View file @ 1f2cbefe
/* -------------------------------------------------------------------------- *
* 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 Stanford University and the Authors. *
* Authors: Peter Eastman *
* Contributors: *
* *
* Permission is hereby granted, free of charge, to any person obtaining a *
* copy of this software and associated documentation files (the "Software"), *
* to deal in the Software without restriction, including without limitation *
* the rights to use, copy, modify, merge, publish, distribute, sublicense, *
* and/or sell copies of the Software, and to permit persons to whom the *
* Software is furnished to do so, subject to the following conditions: *
* *
* The above copyright notice and this permission notice shall be included in *
* all copies or substantial portions of the Software. *
* *
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, *
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL *
* THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, *
* DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR *
* OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
* USE OR OTHER DEALINGS IN THE SOFTWARE. *
* -------------------------------------------------------------------------- */
/**
* This tests the Reference implementation of ReferenceAmoebaVdwForce.
*/
#include "../../../tests/AssertionUtilities.h"
//#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "AmoebaOpenMM.h"
#include "openmm/System.h"
#include "AmoebaVdwForce.h"
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
const double TOL = 1e-4;
void testVdw( FILE* log ) {
System system;
int numberOfParticles = 6;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();
std::string sigmaCombiningRule = std::string("CUBIC-MEAN");
amoebaVdwForce->setSigmaCombiningRule( sigmaCombiningRule );
std::string epsilonCombiningRule = std::string("HHG");
amoebaVdwForce->setEpsilonCombiningRule( epsilonCombiningRule );
for( int ii = 0; ii < numberOfParticles; ii++ ){
int indexIV, indexClass;
double mass, sigma, epsilon, reduction;
std::vector< int > exclusions;
if( ii == 0 || ii == 3 ){
mass = 16.0;
indexIV = ii;
indexClass = 70;
sigma = 1.70250E+00;
epsilon = 1.10000E-01;
reduction = 0.0;
} else {
mass = 1.0;
indexIV = ii < 3 ? 0 : 3;
indexClass = 71;
sigma = 1.32750E+00;
epsilon = 1.35000E-02;
reduction = 0.91;
}
// exclusions
if( ii < 3 ){
exclusions.push_back ( 0 );
exclusions.push_back ( 1 );
exclusions.push_back ( 2 );
} else {
exclusions.push_back ( 3 );
exclusions.push_back ( 4 );
exclusions.push_back ( 5 );
}
system.addParticle(mass);
amoebaVdwForce->addParticle( indexIV, indexClass, sigma, epsilon, reduction );
amoebaVdwForce->setParticleExclusions( ii, exclusions );
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
std::vector<Vec3> positions(numberOfParticles);
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy;
positions[0] = Vec3( -0.254893450E+02, -0.876646600E+01, 0.174761600E+01 );
positions[1] = Vec3( -0.263489690E+02, -0.907798000E+01, 0.205385100E+01 );
positions[2] = Vec3( -0.252491680E+02, -0.949411200E+01, 0.115017600E+01 );
positions[3] = Vec3( 0.172827200E+01, 0.195873090E+02, 0.100059800E+01 );
positions[4] = Vec3( 0.129370700E+01, 0.190112810E+02, 0.169576300E+01 );
positions[5] = Vec3( 0.256122300E+01, 0.191601930E+02, 0.854382000E+00 );
double offset = 27.0;
for( int ii = 0; ii < 3; ii++ ){
positions[ii][0] += offset;
positions[ii][1] += offset;
}
expectedForces[0] = Vec3( -0.729561040E+03, 0.425828484E+04, -0.769114213E+03 );
expectedForces[1] = Vec3( 0.181000041E+02, 0.328216639E+02, -0.126210511E+02 );
expectedForces[2] = Vec3( -0.943743014E+00, 0.199728310E+02, 0.884567842E+00 );
expectedForces[3] = Vec3( 0.615734500E+01, -0.747350431E+03, 0.264726489E+03 );
expectedForces[4] = Vec3( 0.735772031E+03, -0.353310112E+04, 0.490066356E+03 );
expectedForces[5] = Vec3( -0.295245970E+02, -0.306277797E+02, 0.260578506E+02 );
expectedEnergy = 0.740688488E+03;
system.addForce(amoebaVdwForce);
std::string platformName = "Reference";
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
std::vector<Vec3> forces = state.getForces();
const double conversion = -1.0;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
forces[ii][0] *= conversion;
forces[ii][1] *= conversion;
forces[ii][2] *= conversion;
}
if( log ){
(void) fprintf( log, "computeAmoebaVdwForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2], forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
}
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
}
int main( int numberOfArguments, char* argv[] ) {
try {
std::cout << "TestReferenceAmoebaVdwForce running test..." << std::endl;
Platform::loadPluginsFromDirectory( Platform::getDefaultPluginsDirectory() );
FILE* log = NULL;
//FILE* log = stderr;
//FILE* log = fopen( "AmoebaVdwForce1.log", "w" );;
testVdw( log );
if( log && log != stderr )
(void) fclose( log );
}
catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
return 1;
}
std::cout << "PASS - Test succeeded." << std::endl;
return 0;
}
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