Skip to content

GitLab

Commit f0612350 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Added new tests to TestCudaAmoebaVdwForce; fixed bug in AmoebaVdw tapering;... 

Added new tests to TestCudaAmoebaVdwForce; fixed bug in AmoebaVdw tapering; non-default AmoebaVdw rules for calculating combinded epsilon 
value changed to agree w/ TINKER
parent d0426ba9
Hide whitespace changes
Inline Side-by-side
Showing with 415 additions and 49 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaHarmonicAngleForce.h
View file @ f0612350
......@@ -122,7 +122,7 @@ public:
* @param particle1 the index of the first particle connected by the angle
* @param particle2 the index of the second particle connected by the angle
* @param particle3 the index of the third particle connected by the angle
* @param length the angle measured in radians
* @param length the angle measured in degrees
* @param quadratic k the quadratic harmonic force constant for the angle, measured in kJ/mol/radian^2
* @return the index of the angle that was added
*/
......@@ -135,7 +135,7 @@ public:
* @param particle1 the index of the first particle connected by the angle
* @param particle2 the index of the second particle connected by the angle
* @param particle3 the index of the third particle connected by the angle
* @param length the equilibrium angle, measured in radians
* @param length the equilibrium angle, measured in degress
* @param quadratic k the quadratic harmonic force constant for the angle, measured in kJ/mol/radian^2
*/
void getAngleParameters(int index, int& particle1, int& particle2, int& particle3, double& length, double& quadraticK ) const;
......@@ -147,7 +147,7 @@ public:
* @param particle1 the index of the first particle connected by the angle
* @param particle2 the index of the second particle connected by the angle
* @param particle3 the index of the third particle connected by the angle
* @param length the equilibrium angle, measured in radians
* @param length the equilibrium angle, measured in degrees
* @param quadratic k the quadratic harmonic force constant for the angle, measured in kJ/mol/radian^2
*/
void setAngleParameters(int index, int particle1, int particle2, int particle3, double length, double quadraticK );
......
plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp
View file @ f0612350
......@@ -81,11 +81,6 @@ int AmoebaMultipoleForce::getPmeBSplineOrder( void ) const {
return pmeBSplineOrder;
}
/*
void AmoebaMultipoleForce::setPmeBSplineOrder(int inputBSplineOrder) {
pmeBSplineOrder = inputBSplineOrder;
} */
void AmoebaMultipoleForce::getPmeGridDimensions( std::vector<int>& gridDimension ) const {
if( gridDimension.size() < 3 ){
gridDimension.resize(3);
......@@ -107,14 +102,6 @@ void AmoebaMultipoleForce::setPmeGridDimensions( const std::vector<int>& gridDim
pmeGridDimension[2] = gridDimension[2];
return;
}
/*
AmoebaMultipoleForce::MutualInducedIterationMethod AmoebaMultipoleForce::getMutualInducedIterationMethod( void ) const {
return mutualInducedIterationMethod;
}
void AmoebaMultipoleForce::setMutualInducedIterationMethod( AmoebaMultipoleForce::MutualInducedIterationMethod inputMutualInducedIterationMethod ) {
mutualInducedIterationMethod = inputMutualInducedIterationMethod;
} */
int AmoebaMultipoleForce::getMutualInducedMaxIterations( void ) const {
return mutualInducedMaxIterations;
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ f0612350
......@@ -801,7 +801,8 @@ static void computeAmoebaMultipoleForce( AmoebaCudaData& data ) {
data.incrementMultipoleForceCount();
if( 0 && data.getLog() ){
(void) fprintf( data.getLog(), "In computeAmoebaMultipoleForce hasAmoebaGeneralizedKirkwood=%d\n", data.getHasAmoebaGeneralizedKirkwood() );
(void) fprintf( data.getLog(), "In computeAmoebaMultipoleForce hasAmoebaGeneralizedKirkwood=%d\n",
data.getHasAmoebaGeneralizedKirkwood() );
(void) fflush( data.getLog());
}
......@@ -929,7 +930,6 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
std::vector<int> minCovalentIndices(numMultipoles);
std::vector<int> minCovalentPolarizationIndices(numMultipoles);
//float scalingDistanceCutoff = static_cast<float>(force.getScalingDistanceCutoff());
float scalingDistanceCutoff = 50.0f;
std::vector<AmoebaMultipoleForce::CovalentType> covalentList;
......@@ -1007,7 +1007,7 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
}
int polarizationType = static_cast<int>(force.getPolarizationType());
int nonbondedMethod = static_cast<int>(force.getNonbondedMethod());
int nonbondedMethod = static_cast<int>(force.getNonbondedMethod());
if( nonbondedMethod != 0 && nonbondedMethod != 1 ){
throw OpenMMException("AmoebaMultipoleForce nonbonded method not recognized.\n");
}
......@@ -1170,6 +1170,7 @@ static void computeAmoebaVdwForce( AmoebaCudaData& data ) {
// Vdw14_7F
kCalculateAmoebaVdw14_7Forces(gpu, data.getUseVdwNeighborList());
}
/* -------------------------------------------------------------------------- *
......@@ -1181,11 +1182,12 @@ public:
ForceInfo(const AmoebaVdwForce& force) : force(force) {
}
bool areParticlesIdentical(int particle1, int particle2) {
int iv1, iv2, class1, class2;
int iv1, iv2;
int classIndex1, classIndex2;
double sigma1, sigma2, epsilon1, epsilon2, reduction1, reduction2;
force.getParticleParameters(particle1, iv1, class1, sigma1, epsilon1, reduction1);
force.getParticleParameters(particle2, iv2, class2, sigma2, epsilon2, reduction2);
return (class1 == class2 && sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2);
force.getParticleParameters(particle1, iv1, classIndex1, sigma1, epsilon1, reduction1);
force.getParticleParameters(particle2, iv2, classIndex2, sigma2, epsilon2, reduction2);
return (sigma1 == sigma2 && epsilon1 == epsilon2 && reduction1 == reduction2);
}
private:
const AmoebaVdwForce& force;
......@@ -1207,31 +1209,30 @@ void CudaCalcAmoebaVdwForceKernel::initialize(const System& system, const Amoeba
int numParticles = system.getNumParticles();
std::vector<int> indexIVs(numParticles);
std::vector<int> indexClasses(numParticles);
std::vector< std::vector<int> > allExclusions(numParticles);
std::vector<float> sigmas(numParticles);
std::vector<float> epsilons(numParticles);
std::vector<float> reductions(numParticles);
for( int ii = 0; ii < numParticles; ii++ ){
int indexIV, indexClass;
int indexIV;
int classIndex;
double sigma, epsilon, reduction;
std::vector<int> exclusions;
force.getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
force.getParticleParameters( ii, indexIV, classIndex, 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<float>( sigma );
epsilons[ii] = static_cast<float>( epsilon );
reductions[ii] = static_cast<float>( reduction );
}
gpuSetAmoebaVdwParameters( data.getAmoebaGpu(), indexIVs, indexClasses, sigmas, epsilons, reductions,
gpuSetAmoebaVdwParameters( data.getAmoebaGpu(), indexIVs, sigmas, epsilons, reductions,
force.getSigmaCombiningRule(), force.getEpsilonCombiningRule(),
allExclusions, force.getPBC(), static_cast<float>(force.getCutoff()) );
data.getAmoebaGpu()->gpuContext->forces.push_back(new ForceInfo(force));
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaCudaGpu.cpp
View file @ f0612350
......@@ -1886,7 +1886,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
} else if( nonbondedMethod == 1 ){
amoebaGpu->multipoleNonbondedMethod = AMOEBA_PARTICLE_MESH_EWALD;
} else {
throw OpenMM::OpenMMException("multipoleNonbondedMethod not recognized.\n" );
throw OpenMM::OpenMMException("MultipoleNonbondedMethod not recognized.\n" );
}
amoebaGpu->amoebaSim.sqrtPi = std::sqrt( 3.14159265358f );
......@@ -2030,7 +2030,7 @@ void gpuSetAmoebaMultipoleParameters(amoebaGpuContext amoebaGpu, const std::vect
// polarization covalent info
const int minCovalentPolarizationIndex = minCovalentPolarizationIndices[ii];
const int minCovalentPolarizationIndex = minCovalentPolarizationIndices[ii];
amoebaGpu->polarizationDegree[particlesOffset] = minCovalentPolarizationIndex;
for( unsigned int jj = 4; jj < covalentInfo.size(); jj++ ){
......@@ -2463,7 +2463,6 @@ static void lookupVdwTaper( float r, double vdwTaperCut, double delta,
extern "C"
void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
const std::vector<int>& indexIVs,
const std::vector<int>& indexClasses,
const std::vector<float>& sigmas,
const std::vector<float>& epsilons,
const std::vector<float>& reductions,
......@@ -2484,7 +2483,7 @@ void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
amoebaGpu->amoebaSim.vdwUsePBC = usePBC;
amoebaGpu->amoebaSim.vdwCutoff = cutoff;
amoebaGpu->amoebaSim.vdwCutoff2 = cutoff*cutoff;
double vdwTaper = 0.90f;
double vdwTaper = 0.90f;
if( vdwTaper < 1.0 ){
double vdwTaperCoefficients[6];
double vdwCut = cutoff;
......
plugins/amoeba/platforms/cuda/src/kernels/amoebaGpuTypes.h
View file @ f0612350
......@@ -310,7 +310,6 @@ void gpuSetAmoebaGrycukParameters( amoebaGpuContext amoebaGpu , float innerDiele
extern "C"
void gpuSetAmoebaVdwParameters( amoebaGpuContext amoebaGpu,
const std::vector<int>& indexIVs,
const std::vector<int>& indexClasses,
const std::vector<float>& sigmas,
const std::vector<float>& epsilons,
const std::vector<float>& reductions,
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaElectrostatic.cu
View file @ f0612350
......@@ -110,7 +110,7 @@ struct ElectrostaticParticle {
#define one 1.0f
__device__ void calculateElectrostaticPairIxnOrig_kernel( ElectrostaticParticle& atomI, ElectrostaticParticle& atomJ,
float* scalingFactors, float4* outputForce, float4 outputTorque[2]){
float* scalingFactors, float4* outputForce, float4 outputTorque[2]){
float deltaR[3];
......
plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaVdw14_7.cu
View file @ f0612350
......@@ -91,6 +91,10 @@ __device__ void loadVdw14_7Shared( struct Vdw14_7Particle* sA, unsigned int atom
__device__ void getVdw14_7CombindedSigmaEpsilon_kernel( int sigmaCombiningRule, float iSigma, float jSigma, float* combindedSigma,
int epsilonCombiningRule, float iEpsilon, float jEpsilon, float* combindedEpsilon )
{
// ARITHMETIC = 1
// GEOMETRIC = 2
// CUBIC-MEAN = 3
if( sigmaCombiningRule == 1 ){
*combindedSigma = iSigma + jSigma;
} else if( sigmaCombiningRule == 2 ){
......@@ -101,14 +105,17 @@ __device__ void getVdw14_7CombindedSigmaEpsilon_kernel( int sigmaCombiningRule,
*combindedSigma = 2.0f*( iSigma2*iSigma + jSigma2*jSigma )/( iSigma2 + jSigma2 );
}
// ARITHMETIC = 1
// GEOMETRIC = 2
// HARMONIC = 3
// HHG = 4
if( epsilonCombiningRule == 1 ){
*combindedEpsilon = iEpsilon + jEpsilon;
*combindedEpsilon = 0.5f*(iEpsilon + jEpsilon);
} else if( epsilonCombiningRule == 2 ){
*combindedEpsilon = 2.0f*sqrtf( iEpsilon*jEpsilon );
*combindedEpsilon = sqrtf( iEpsilon*jEpsilon );
} else if( epsilonCombiningRule == 3 ){
float iEpsilon2 = iEpsilon*iEpsilon;
float jEpsilon2 = jEpsilon*jEpsilon;
*combindedEpsilon = 2.0f*( iEpsilon2*iEpsilon + jEpsilon2*jEpsilon )/( iEpsilon2 + jEpsilon2 );
*combindedEpsilon = 2.0f*( iEpsilon*jEpsilon )/( iEpsilon + jEpsilon );
} else {
float epsilonS = sqrtf( iEpsilon ) + sqrtf( jEpsilon );
*combindedEpsilon = 4.0f*( iEpsilon*jEpsilon )/( epsilonS*epsilonS );
......@@ -215,7 +222,7 @@ __device__ void calculateVdw14_7PairIxn_kernel( float combindedSigma, float c
float gTau = combindedEpsilon*tau7*r6*gammaHal*tmp*tmp;
*energy = combindedEpsilon*combindedSigma7*tau7*( (combindedSigma7*gammaHal*rhoInverse) - 2.0f);
float deltaE = (-7.0f*(dTau*(*energy) + gTau))*rI;
float deltaE = (-7.0f*(dTau*(*energy) + gTau));
if( r > cAmoebaSim.vdwTaperCutoff ){
......@@ -226,11 +233,12 @@ __device__ void calculateVdw14_7PairIxn_kernel( float combindedSigma, float c
*energy *= taper;
}
deltaE *= rI;
force[0] *= deltaE;
force[1] *= deltaE;
force[2] *= deltaE;
}
// perform reduction of force on H's and add to heavy atom partner
......
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaVdwForce.cpp
View file @ f0612350
......@@ -42,6 +42,13 @@
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
#include <stdlib.h>
#include <stdio.h>
#define ASSERT_EQUAL_TOL_MOD(expected, found, tol, testname) {double _scale_ = std::abs(expected) > 1.0 ? std::abs(expected) : 1.0; if (!(std::abs((expected)-(found))/_scale_ <= (tol))) {std::stringstream details; details << testname << " Expected "<<(expected)<<", found "<<(found); throwException(__FILE__, __LINE__, details.str());}};
#define ASSERT_EQUAL_VEC_MOD(expected, found, tol,testname) {ASSERT_EQUAL_TOL_MOD((expected)[0], (found)[0], (tol),(testname)); ASSERT_EQUAL_TOL_MOD((expected)[1], (found)[1], (tol),(testname)); ASSERT_EQUAL_TOL_MOD((expected)[2], (found)[2], (tol),(testname));};
using namespace OpenMM;
const double TOL = 1e-4;
......@@ -56,22 +63,20 @@ void testVdw( FILE* log ) {
std::string epsilonCombiningRule = std::string("HHG");
amoebaVdwForce->setEpsilonCombiningRule( epsilonCombiningRule );
int classIndex = 0;
for( int ii = 0; ii < numberOfParticles; ii++ ){
int indexIV, indexClass;
int indexIV;
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;
......@@ -89,7 +94,7 @@ void testVdw( FILE* log ) {
exclusions.push_back ( 5 );
}
system.addParticle(mass);
amoebaVdwForce->addParticle( indexIV, indexClass, sigma, epsilon, reduction );
amoebaVdwForce->addParticle( indexIV, classIndex, sigma, epsilon, reduction );
amoebaVdwForce->setParticleExclusions( ii, exclusions );
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
......@@ -130,12 +135,13 @@ void testVdw( FILE* log ) {
positions[ii][2] *= AngstromToNm;
}
for( int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++ ){
int indexIV, indexClass;
int indexIV;
int classIndex;
double sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
amoebaVdwForce->getParticleParameters( ii, indexIV, classIndex, sigma, epsilon, reduction );
sigma *= AngstromToNm;
epsilon *= CalToJoule;
amoebaVdwForce->setParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
amoebaVdwForce->setParticleParameters( ii, indexIV, classIndex, sigma, epsilon, reduction );
}
platformName = "Cuda";
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
......@@ -170,6 +176,342 @@ void testVdw( FILE* log ) {
ASSERT_EQUAL_TOL( expectedEnergy, state.getPotentialEnergy(), tolerance );
}
void setupAndGetForcesEnergyVdwAmmonia( const std::string& sigmaCombiningRule, const std::string& epsilonCombiningRule, double cutoff,
double boxDimension, std::vector<Vec3>& forces, double& energy, FILE* log ){
// beginning of Vdw setup
System system;
AmoebaVdwForce* amoebaVdwForce = new AmoebaVdwForce();;
int numberOfParticles = 8;
amoebaVdwForce->setSigmaCombiningRule( sigmaCombiningRule );
amoebaVdwForce->setEpsilonCombiningRule( epsilonCombiningRule );
amoebaVdwForce->setUseNeighborList( 1 );
amoebaVdwForce->setCutoff( cutoff );
if( boxDimension > 0.0 ){
Vec3 a( boxDimension, 0.0, 0.0 );
Vec3 b( 0.0, boxDimension, 0.0 );
Vec3 c( 0.0, 0.0, boxDimension );
system.setDefaultPeriodicBoxVectors( a, b, c );
amoebaVdwForce->setPBC( 1 );
} else {
amoebaVdwForce->setPBC( 0 );
}
// addParticle: ivIndex, radius, epsilon, reductionFactor
int classIndex = 0;
system.addParticle( 1.4007000e+01 );
amoebaVdwForce->addParticle( 0, classIndex, 1.8550000e-01, 4.3932000e-01, 0.0000000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( 0, classIndex, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( 0, classIndex, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( 0, classIndex, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01 );
system.addParticle( 1.4007000e+01 );
amoebaVdwForce->addParticle( 4, classIndex, 1.8550000e-01, 4.3932000e-01, 0.0000000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( 4, classIndex, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( 4, classIndex, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( 4, classIndex, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01 );
// ParticleExclusions
std::vector< int > exclusions;
exclusions.resize(0);
exclusions.push_back( 0 );
exclusions.push_back( 1 );
exclusions.push_back( 2 );
exclusions.push_back( 3 );
amoebaVdwForce->setParticleExclusions( 0, exclusions );
exclusions.resize(0);
exclusions.push_back( 1 );
exclusions.push_back( 0 );
exclusions.push_back( 2 );
exclusions.push_back( 3 );
amoebaVdwForce->setParticleExclusions( 1, exclusions );
exclusions.resize(0);
exclusions.push_back( 2 );
exclusions.push_back( 0 );
exclusions.push_back( 1 );
exclusions.push_back( 3 );
amoebaVdwForce->setParticleExclusions( 2, exclusions );
exclusions.resize(0);
exclusions.push_back( 3 );
exclusions.push_back( 0 );
exclusions.push_back( 1 );
exclusions.push_back( 2 );
amoebaVdwForce->setParticleExclusions( 3, exclusions );
exclusions.resize(0);
exclusions.push_back( 4 );
exclusions.push_back( 5 );
exclusions.push_back( 6 );
exclusions.push_back( 7 );
amoebaVdwForce->setParticleExclusions( 4, exclusions );
exclusions.resize(0);
exclusions.push_back( 5 );
exclusions.push_back( 4 );
exclusions.push_back( 6 );
exclusions.push_back( 7 );
amoebaVdwForce->setParticleExclusions( 5, exclusions );
exclusions.resize(0);
exclusions.push_back( 6 );
exclusions.push_back( 4 );
exclusions.push_back( 5 );
exclusions.push_back( 7 );
amoebaVdwForce->setParticleExclusions( 6, exclusions );
exclusions.resize(0);
exclusions.push_back( 7 );
exclusions.push_back( 4 );
exclusions.push_back( 5 );
exclusions.push_back( 6 );
amoebaVdwForce->setParticleExclusions( 7, exclusions );
// end of Vdw setup
std::vector<Vec3> positions(numberOfParticles);
positions[0] = Vec3( 1.5927280e-01, 1.7000000e-06, 1.6491000e-03 );
positions[1] = Vec3( 2.0805540e-01, -8.1258800e-02, 3.7282500e-02 );
positions[2] = Vec3( 2.0843610e-01, 8.0953200e-02, 3.7462200e-02 );
positions[3] = Vec3( 1.7280780e-01, 2.0730000e-04, -9.8741700e-02 );
positions[4] = Vec3( -1.6743680e-01, 1.5900000e-05, -6.6149000e-03 );
positions[5] = Vec3( -2.0428260e-01, 8.1071500e-02, 4.1343900e-02 );
positions[6] = Vec3( -6.7308300e-02, 1.2800000e-05, 1.0623300e-02 );
positions[7] = Vec3( -2.0426290e-01, -8.1231400e-02, 4.1033500e-02 );
system.addForce(amoebaVdwForce);
std::string platformName;
platformName = "Cuda";
LangevinIntegrator integrator(0.0, 0.1, 0.01);
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
forces = state.getForces();
energy = state.getPotentialEnergy();
}
void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) {
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e\n", testName.c_str(), 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 );
}
#endif
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC_MOD( expectedForces[ii], forces[ii], tolerance, testName );
}
ASSERT_EQUAL_TOL_MOD( expectedEnergy, energy, tolerance, testName );
}
// test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG
void testVdwAmmoniaCubicMeanHhg( FILE* log ) {
std::string testName = "testVdwAmmoniaCubicMeanHhg";
int numberOfParticles = 8;
double boxDimension = -1.0;
double cutoff = 9000000.0;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.8012258e+00;
expectedForces[0] = Vec3( 2.9265247e+02, -1.4507808e-02, -6.9562123e+00 );
expectedForces[1] = Vec3( -2.2451693e+00, 4.8143073e-01, -2.0041494e-01 );
expectedForces[2] = Vec3( -2.2440698e+00, -4.7905450e-01, -2.0125284e-01 );
expectedForces[3] = Vec3( -1.0840394e+00, -5.8531253e-04, 2.6934135e-01 );
expectedForces[4] = Vec3( -5.6305662e+01, 1.4733908e-03, -1.8083306e-01 );
expectedForces[5] = Vec3( 1.6750145e+00, -3.2448374e-01, -1.8030914e-01 );
expectedForces[6] = Vec3( -2.3412420e+02, 1.0754069e-02, 7.6287492e+00 );
expectedForces[7] = Vec3( 1.6756544e+00, 3.2497316e-01, -1.7906832e-01 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
void testVdwAmmoniaArithmeticArithmetic( FILE* log ) {
std::string testName = "testVdwAmmoniaArithmeticArithmetic";
int numberOfParticles = 8;
double boxDimension = -1.0;
double cutoff = 9000000.0;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia( "ARITHMETIC", "ARITHMETIC", cutoff, boxDimension, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.2252403e+00;
expectedForces[0] = Vec3( 3.0603839e+02, -1.5550310e-02, -7.2661707e+00 );
expectedForces[1] = Vec3( -2.7801357e+00, 5.8805051e-01, -2.5907269e-01 );
expectedForces[2] = Vec3( -2.7753968e+00, -5.8440732e-01, -2.5969111e-01 );
expectedForces[3] = Vec3( -2.2496416e+00, -1.1797440e-03, 5.5501757e-01 );
expectedForces[4] = Vec3( -5.5077629e+01, 8.3417114e-04, -3.3668921e-01 );
expectedForces[5] = Vec3( 2.3752452e+00, -4.6788669e-01, -2.4907764e-01 );
expectedForces[6] = Vec3( -2.4790697e+02, 1.1419770e-02, 8.0629999e+00 );
expectedForces[7] = Vec3( 2.3761408e+00, 4.6871961e-01, -2.4731607e-01 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
// test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric
void testVdwAmmoniaGeometricGeometric( FILE* log ) {
std::string testName = "testVdwAmmoniaGeometricGeometric";
int numberOfParticles = 8;
double boxDimension = -1.0;
double cutoff = 9000000.0;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia( "GEOMETRIC", "GEOMETRIC", cutoff, boxDimension, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 2.5249914e+00;
expectedForces[0] = Vec3( 2.1169631e+02, -1.0710925e-02, -4.3728025e+00 );
expectedForces[1] = Vec3( -2.2585621e+00, 4.8409995e-01, -2.0188344e-01 );
expectedForces[2] = Vec3( -2.2551351e+00, -4.8124855e-01, -2.0246986e-01 );
expectedForces[3] = Vec3( -1.7178028e+00, -9.0851787e-04, 4.2466975e-01 );
expectedForces[4] = Vec3( -4.8302147e+01, 9.6603376e-04, -5.7972068e-01 );
expectedForces[5] = Vec3( 1.8100634e+00, -3.5214093e-01, -1.9357207e-01 );
expectedForces[6] = Vec3( -1.6078365e+02, 7.2117601e-03, 5.3180261e+00 );
expectedForces[7] = Vec3( 1.8109211e+00, 3.5273117e-01, -1.9224723e-01 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
void testVdwAmmoniaCubicMeanHarmonic( FILE* log ) {
std::string testName = "testVdwAmmoniaCubicMeanHarmonic";
int numberOfParticles = 8;
double boxDimension = -1.0;
double cutoff = 9000000.0;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HARMONIC", cutoff, boxDimension, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 4.1369069e+00;
expectedForces[0] = Vec3( 2.5854436e+02, -1.2779529e-02, -5.9041148e+00 );
expectedForces[1] = Vec3( -2.0832419e+00, 4.4915831e-01, -1.8266000e-01 );
expectedForces[2] = Vec3( -2.0823991e+00, -4.4699804e-01, -1.8347141e-01 );
expectedForces[3] = Vec3( -9.5914714e-01, -5.2162026e-04, 2.3873165e-01 );
expectedForces[4] = Vec3( -5.3724787e+01, 1.4838241e-03, -2.8089191e-01 );
expectedForces[5] = Vec3( 1.5074325e+00, -2.9016397e-01, -1.6385118e-01 );
expectedForces[6] = Vec3( -2.0271029e+02, 9.2367947e-03, 6.6389988e+00 );
expectedForces[7] = Vec3( 1.5080748e+00, 2.9058422e-01, -1.6274118e-01 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
// test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on
// particle 4 due to reduction applied to NH
// the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region
void testVdwTaper( FILE* log ) {
std::string testName = "testVdwTaper";
int numberOfParticles = 8;
double boxDimension = -1.0;
double cutoff = 0.25;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 3.5478444e+00;
expectedForces[0] = Vec3( 5.6710779e+02, -2.7391004e-02, -1.7867730e+01 );
expectedForces[1] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 );
expectedForces[2] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 );
expectedForces[3] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 );
expectedForces[4] = Vec3( -5.1039701e+01, 2.4651903e-03, 1.6080957e+00 );
expectedForces[5] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 );
expectedForces[6] = Vec3( -5.1606809e+02, 2.4925813e-02, 1.6259634e+01 );
expectedForces[7] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
// test PBC
void testVdwPBC( FILE* log ) {
std::string testName = "testVdwPBC";
int numberOfParticles = 8;
double boxDimension = 0.6;
double cutoff = 0.25;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyVdwAmmonia( "CUBIC-MEAN", "HHG", cutoff, boxDimension, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = 8.4385405e+00;
expectedForces[0] = Vec3( 5.1453069e+02, 4.9751912e-01, -1.2759570e+01 );
expectedForces[1] = Vec3( -2.5622586e+02, -4.6524265e+01, 2.4281465e+01 );
expectedForces[2] = Vec3( -2.7538705e+02, 5.1831690e+01, 2.7367710e+01 );
expectedForces[3] = Vec3( -0.0000000e+00, -0.0000000e+00, -0.0000000e+00 );
expectedForces[4] = Vec3( 3.0883034e+02, -5.8876974e+00, -5.8286122e+01 );
expectedForces[5] = Vec3( 1.1319359e+02, -3.2047069e-01, 1.6181231e+00 );
expectedForces[6] = Vec3( -5.1606809e+02, 2.4925813e-02, 1.6259634e+01 );
expectedForces[7] = Vec3( 1.1112638e+02, 3.7829857e-01, 1.5187587e+00 );
// tolerance is higher here due to interpolation used in setting tapering coefficients;
// if tapering turned off, then absolute difference < 2.0e-05
double tolerance = 5.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
int main( int numberOfArguments, char* argv[] ) {
......@@ -178,7 +520,37 @@ int main( int numberOfArguments, char* argv[] ) {
registerAmoebaCudaKernelFactories();
FILE* log = NULL;
testVdw( log );
// tests using two ammonia molecules
// test VDW w/ sigmaRule=CubicMean and epsilonRule=HHG
testVdwAmmoniaCubicMeanHhg( log );
// test VDW w/ sigmaRule=Arithmetic and epsilonRule=Arithmetic
testVdwAmmoniaArithmeticArithmetic( log );
// test VDW w/ sigmaRule=Geometric and epsilonRule=Geometric
testVdwAmmoniaGeometricGeometric( log );
// test VDW w/ sigmaRule=CubicMean and epsilonRule=Harmonic
testVdwAmmoniaCubicMeanHarmonic( log );
// test w/ cutoff=0.25 nm; single ixn between two particles (0 and 6); force nonzero on
// particle 4 due to reduction applied to NH
// the distance between 0 and 6 is ~ 0.235 so the ixn is in the tapered region
testVdwTaper( log );
// test PBC
testVdwPBC( log );
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::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