Skip to content

GitLab

Commit 29654b80 authored by Mark Friedrichs's avatar Mark Friedrichs
Browse files

Added unit tests for Amoeba GeneralizedKirkwood and WcaDispersion forces on Cuda platform

parent e5b7f0b3
Show whitespace changes
Inline Side-by-side
Showing with 657 additions and 0 deletions
plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp 0 → 100644
View file @ 29654b80
/* -------------------------------------------------------------------------- *
* OpenMMAmoeba *
* -------------------------------------------------------------------------- *
* 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 Cuda implementation of CudaAmoebaMultipoleForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/System.h"
#include "openmm/AmoebaMultipoleForce.h"
#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;
// setup for 2 ammonia molecules
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){
// beginning of Multipole setup
System system;
AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
int numberOfParticles = 8;
amoebaMultipoleForce->setNonbondedMethod( AmoebaMultipoleForce::NoCutoff );
amoebaMultipoleForce->setPolarizationType( polarizationType );
amoebaMultipoleForce->setMutualInducedTargetEpsilon( 1.0e-06 );
amoebaMultipoleForce->setMutualInducedMaxIterations( 500 );
std::vector<double> nitrogenMolecularDipole(3);
std::vector<double> nitrogenMolecularQuadrupole(9);
nitrogenMolecularDipole[0] = 8.3832254e-03;
nitrogenMolecularDipole[1] = 0.0000000e+00;
nitrogenMolecularDipole[2] = 3.4232474e-03;
nitrogenMolecularQuadrupole[0] = -4.0406249e-04;
nitrogenMolecularQuadrupole[1] = 0.0000000e+00;
nitrogenMolecularQuadrupole[2] = -2.6883671e-04;
nitrogenMolecularQuadrupole[3] = 0.0000000e+00;
nitrogenMolecularQuadrupole[4] = 2.5463927e-04;
nitrogenMolecularQuadrupole[5] = 0.0000000e+00;
nitrogenMolecularQuadrupole[6] = -2.6883671e-04;
nitrogenMolecularQuadrupole[7] = 0.0000000e+00;
nitrogenMolecularQuadrupole[8] = 1.4942322e-04;
// first N
system.addParticle( 1.4007000e+01 );
amoebaMultipoleForce->addParticle( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 1, 2, 3, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 );
// 3 H attached to first N
std::vector<double> hydrogenMolecularDipole(3);
std::vector<double> hydrogenMolecularQuadrupole(9);
hydrogenMolecularDipole[0] = -1.7388763e-03;
hydrogenMolecularDipole[1] = 0.0000000e+00;
hydrogenMolecularDipole[2] = -4.6837475e-03;
hydrogenMolecularQuadrupole[0] = -4.4253841e-05;
hydrogenMolecularQuadrupole[1] = 0.0000000e+00;
hydrogenMolecularQuadrupole[2] = 1.5429571e-05;
hydrogenMolecularQuadrupole[3] = 0.0000000e+00;
hydrogenMolecularQuadrupole[4] = 4.1798924e-05;
hydrogenMolecularQuadrupole[5] = 0.0000000e+00;
hydrogenMolecularQuadrupole[6] = 1.5429571e-05;
hydrogenMolecularQuadrupole[7] = 0.0000000e+00;
hydrogenMolecularQuadrupole[8] = 2.4549167e-06;
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 2, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 3, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 0, 1, 2, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
// second N
system.addParticle( 1.4007000e+01 );
amoebaMultipoleForce->addParticle( -5.7960000e-01, nitrogenMolecularDipole, nitrogenMolecularQuadrupole, 2, 5, 6, 7, 3.9000000e-01, 3.1996314e-01, 1.0730000e-03 );
// 3 H attached to second N
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 6, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 7, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
amoebaMultipoleForce->addParticle( 1.9320000e-01, hydrogenMolecularDipole, hydrogenMolecularQuadrupole, 2, 4, 5, 6, 3.9000000e-01, 2.8135002e-01, 4.9600000e-04 );
// covalent maps
std::vector< int > covalentMap;
covalentMap.resize(0);
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 0, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 0 );
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 1, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 0 );
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 1 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 2, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 0 );
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 0 );
covalentMap.push_back( 1 );
covalentMap.push_back( 2 );
covalentMap.push_back( 3 );
amoebaMultipoleForce->setCovalentMap( 3, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 4, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 4 );
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 5, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 4 );
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 5 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 6, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 4 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(0), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(1), covalentMap );
covalentMap.resize(0);
covalentMap.push_back( 4 );
covalentMap.push_back( 5 );
covalentMap.push_back( 6 );
covalentMap.push_back( 7 );
amoebaMultipoleForce->setCovalentMap( 7, static_cast<OpenMM::AmoebaMultipoleForce::CovalentType>(4), covalentMap );
system.addForce(amoebaMultipoleForce);
// GK force
AmoebaGeneralizedKirkwoodForce* amoebaGeneralizedKirkwoodForce = new AmoebaGeneralizedKirkwoodForce();
amoebaGeneralizedKirkwoodForce->setSolventDielectric( 7.8300000e+01 );
amoebaGeneralizedKirkwoodForce->setSoluteDielectric( 1.0000000e+00 );
amoebaGeneralizedKirkwoodForce->setIncludeCavityTerm( includeCavityTerm );
// addParticle: charge, radius, scalingFactor
for( unsigned int ii = 0; ii < 2; ii++ ){
amoebaGeneralizedKirkwoodForce->addParticle( -5.7960000e-01, 1.5965000e-01, 6.9000000e-01 );
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
amoebaGeneralizedKirkwoodForce->addParticle( 1.9320000e-01, 1.2360000e-01, 6.9000000e-01 );
}
system.addForce(amoebaGeneralizedKirkwoodForce);
// 1-2 bonds needed
/*
AmoebaHarmonicBondForce* amoebaHarmonicBondForce = new AmoebaHarmonicBondForce();
// addBond: particle1, particle2, length, quadraticK
amoebaHarmonicBondForce->addBond( 0, 1, 0.0000000e+00, 0.0000000e+00 );
amoebaHarmonicBondForce->addBond( 0, 2, 0.0000000e+00, 0.0000000e+00 );
amoebaHarmonicBondForce->addBond( 0, 3, 0.0000000e+00, 0.0000000e+00 );
amoebaHarmonicBondForce->addBond( 4, 5, 0.0000000e+00, 0.0000000e+00 );
amoebaHarmonicBondForce->addBond( 4, 6, 0.0000000e+00, 0.0000000e+00 );
amoebaHarmonicBondForce->addBond( 4, 7, 0.0000000e+00, 0.0000000e+00 );
amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondCubic( -2.5500000e+01 );
amoebaHarmonicBondForce->setAmoebaGlobalHarmonicBondQuartic( 3.7931250e+02 );
system.addForce(amoebaHarmonicBondForce);
*/
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 );
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();
}
// compare forces and energies
static void compareForcesEnergy( std::string& testName, double expectedEnergy, double energy,
std::vector<Vec3>& expectedForces,
std::vector<Vec3>& forces, double tolerance, FILE* log ) {
//#define AMOEBA_DEBUG
#ifdef AMOEBA_DEBUG
if( log ){
double conversion = 1.0/4.184;
double energyAbsDiff = fabs( expectedEnergy - energy );
double energyRelDiff = 2.0*energyAbsDiff/( fabs( expectedEnergy ) + fabs( energy ) + 1.0e-08 );
(void) fprintf( log, "%s: expected energy=%14.7e %14.7e absDiff=%15.7e relDiff=%15.7e\n", testName.c_str(), conversion*expectedEnergy, conversion*energy,
conversion*energyAbsDiff, conversion*energyRelDiff );
if( conversion != 1.0 )conversion *= -0.1;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
double expectedNorm = sqrt( expectedForces[ii][0]*expectedForces[ii][0] +
expectedForces[ii][1]*expectedForces[ii][1] +
expectedForces[ii][2]*expectedForces[ii][2] );
double norm = sqrt( forces[ii][0]*forces[ii][0] + forces[ii][1]*forces[ii][1] + forces[ii][2]*forces[ii][2] );
double absDiff = fabs( norm - expectedNorm );
double relDiff = 2.0*absDiff/(fabs( norm ) + fabs( expectedNorm ) + 1.0e-08);
(void) fprintf( log, "%6u %15.7e %15.7e [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
conversion*absDiff, conversion*relDiff,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*forces[ii][2], conversion*expectedNorm, conversion*norm );
}
(void) fflush( log );
conversion = 1.0;
(void) fprintf( log, "\n%s: expected energy=%14.7e %14.7e no conversion\n", testName.c_str(), conversion*expectedEnergy, conversion*energy );
if( conversion != 1.0 )conversion = -1.0;
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,
conversion*expectedForces[ii][0], conversion*expectedForces[ii][1], conversion*expectedForces[ii][2],
conversion*forces[ii][0], conversion*forces[ii][1], conversion*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 GK direct polarization for system comprised of two ammonia molecules
static void testGeneralizedKirkwoodAmmoniaDirectPolarization( FILE* log ) {
std::string testName = "testGeneralizedKirkwoodAmmoniaDirectPolarization";
int numberOfParticles = 8;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Direct, 0, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -7.6636680e+01;
expectedForces[0] = Vec3( -6.9252994e+02, -8.9085133e+00, 9.6489739e+01 );
expectedForces[1] = Vec3( 1.5593797e+02, -6.0331931e+01, 1.5104507e+01 );
expectedForces[2] = Vec3( 1.5870088e+02, 6.1702809e+01, 6.7708985e+00 );
expectedForces[3] = Vec3( 1.4089885e+02, 7.5870617e+00, -1.1362294e+02 );
expectedForces[4] = Vec3( -1.8916205e+02, 2.1465549e-01, -4.3433152e+02 );
expectedForces[5] = Vec3( 1.0208290e+01, 6.2676753e+01, 1.4987953e+02 );
expectedForces[6] = Vec3( 4.0621859e+02, 1.8962203e-01, 1.3021956e+02 );
expectedForces[7] = Vec3( 9.7274235e+00, -6.3130458e+01, 1.4949024e+02 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
// test GK mutual polarization for system comprised of two ammonia molecules
static void testGeneralizedKirkwoodAmmoniaMutualPolarization( FILE* log ) {
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarization";
int numberOfParticles = 8;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Mutual, 0, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -7.8018875e+01;
expectedForces[0] = Vec3( -7.6820301e+02, -1.0102760e+01, 1.0094389e+02 );
expectedForces[1] = Vec3( 1.7037307e+02, -7.5621857e+01, 2.3320365e+01 );
expectedForces[2] = Vec3( 1.7353828e+02, 7.7199741e+01, 1.3965379e+01 );
expectedForces[3] = Vec3( 1.5045244e+02, 8.5784569e+00, -1.3377619e+02 );
expectedForces[4] = Vec3( -2.1811615e+02, -1.6818022e-01, -4.6103163e+02 );
expectedForces[5] = Vec3( 6.2091942e+00, 7.6748687e+01, 1.5883463e+02 );
expectedForces[6] = Vec3( 4.8035662e+02, 4.9704902e-01, 1.3948083e+02 );
expectedForces[7] = Vec3( 5.3895456e+00, -7.7131137e+01, 1.5826273e+02 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
// test GK mutual polarization for system comprised of two ammonia molecules
static void testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( FILE* log ) {
std::string testName = "testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm";
int numberOfParticles = 8;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Mutual, 1, forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
double expectedEnergy = -6.0434582e+01;
expectedForces[0] = Vec3( -7.8323218e+02, -1.0097644e+01, 1.0256890e+02 );
expectedForces[1] = Vec3( 1.7078480e+02, -7.1896701e+01, 2.0840172e+01 );
expectedForces[2] = Vec3( 1.7394089e+02, 7.3488594e+01, 1.1484648e+01 );
expectedForces[3] = Vec3( 1.5169364e+02, 8.5611299e+00, -1.2968050e+02 );
expectedForces[4] = Vec3( -2.1669693e+02, -1.5926823e-01, -4.6636274e+02 );
expectedForces[5] = Vec3( 8.7397444e+00, 7.3330990e+01, 1.6016898e+02 );
expectedForces[6] = Vec3( 4.8684950e+02, 4.8937161e-01, 1.4137061e+02 );
expectedForces[7] = Vec3( 7.9205382e+00, -7.3716473e+01, 1.5960993e+02 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
int main( int numberOfArguments, char* argv[] ) {
try {
std::cout << "TestCudaAmoebaMultipoleForce running test..." << std::endl;
registerAmoebaCudaKernelFactories();
FILE* log = NULL;
// test direct and mutual polarization cases and
// mutual polarization w/ the cavity term
testGeneralizedKirkwoodAmmoniaDirectPolarization( log );
testGeneralizedKirkwoodAmmoniaMutualPolarization( log );
testGeneralizedKirkwoodAmmoniaMutualPolarizationWithCavityTerm( 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 << "Done" << std::endl;
return 0;
}
plugins/amoeba/platforms/cuda/tests/TestCudaWcaDispersionForce.cpp 0 → 100644
View file @ 29654b80
/* -------------------------------------------------------------------------- *
* OpenMMAmoeba *
* -------------------------------------------------------------------------- *
* 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 Cuda implementation of CudaAmoebaWcaDispersionForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/System.h"
#include "openmm/AmoebaWcaDispersionForce.h"
#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;
void setupAndGetForcesEnergyWcaDispersionAmmonia( std::vector<Vec3>& forces, double& energy, FILE* log ){
// beginning of WcaDispersion setup
System system;
AmoebaWcaDispersionForce* amoebaWcaDispersionForce = new AmoebaWcaDispersionForce();;
int numberOfParticles = 8;
amoebaWcaDispersionForce->setEpso( 4.6024000e-01 );
amoebaWcaDispersionForce->setEpsh( 5.6484000e-02 );
amoebaWcaDispersionForce->setRmino( 1.7025000e-01 );
amoebaWcaDispersionForce->setRminh( 1.3275000e-01 );
amoebaWcaDispersionForce->setDispoff( 2.6000000e-02 );
amoebaWcaDispersionForce->setAwater( 3.3428000e+01 );
amoebaWcaDispersionForce->setSlevy( 1.0000000e+00 );
amoebaWcaDispersionForce->setShctd( 8.1000000e-01 );
// addParticle: radius, epsilon
for( unsigned int ii = 0; ii < 2; ii++ ){
system.addParticle( 1.4007000e+01 );
amoebaWcaDispersionForce->addParticle( 1.8550000e-01, 4.3932000e-01 );
system.addParticle( 1.0080000e+00 );
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 );
system.addParticle( 1.0080000e+00 );
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 );
system.addParticle( 1.0080000e+00 );
amoebaWcaDispersionForce->addParticle( 1.3500000e-01, 8.3680000e-02 );
}
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(amoebaWcaDispersionForce);
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, energy );
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 Wca dispersion
void testWcaDispersionAmmonia( FILE* log ) {
std::string testName = "testWcaDispersionAmmonia";
int numberOfParticles = 8;
std::vector<Vec3> forces;
double energy;
setupAndGetForcesEnergyWcaDispersionAmmonia( forces, energy, log );
std::vector<Vec3> expectedForces(numberOfParticles);
// TINKER-computed values
double expectedEnergy = -2.6981209e+01;
expectedForces[0] = Vec3( 4.7839388e+00, -7.3510133e-04, -5.0382764e-01 );
expectedForces[1] = Vec3( 1.4657758e+00, 1.2431003e+00, -6.7075886e-01 );
expectedForces[2] = Vec3( 1.4563936e+00, -1.2399917e+00, -6.7443841e-01 );
expectedForces[3] = Vec3( 2.1116744e+00, -2.7407512e-03, 1.3271245e+00 );
expectedForces[4] = Vec3( -4.7528440e+00, -1.5148066e-03, 1.2653813e+00 );
expectedForces[5] = Vec3( -1.1875619e+00, -1.2866678e+00, -3.9109060e-01 );
expectedForces[6] = Vec3( -2.6885679e+00, -4.3038639e-04, 3.3763583e-02 );
expectedForces[7] = Vec3( -1.1888087e+00, 1.2889802e+00, -3.8615387e-01 );
double tolerance = 1.0e-04;
compareForcesEnergy( testName, expectedEnergy, energy, expectedForces, forces, tolerance, log );
}
int main( int numberOfArguments, char* argv[] ) {
try {
std::cout << "TestCudaAmoebaWcaDispersionForce running test..." << std::endl;
registerAmoebaCudaKernelFactories();
FILE* log = NULL;
// test Wca dispersion force using two ammonia molecules
testWcaDispersionAmmonia( 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 << "Done" << 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