Skip to content

GitLab

Commit 1db349e5 authored by Peter Eastman's avatar Peter Eastman
Browse files

Major cleanup of the AMOEBA API

parent a919f305
Show whitespace changes
Inline Side-by-side
Showing with 68 additions and 953 deletions
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp
View file @ 1db349e5
......@@ -36,7 +36,6 @@
#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"
......@@ -53,9 +52,11 @@
using namespace OpenMM;
const double TOL = 1e-4;
 
extern "C" void registerAmoebaCudaKernelFactories();
// setup for 2 ammonia molecules
 
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){
 
// beginning of Multipole setup
......@@ -316,7 +317,7 @@ static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::Amoeb
 
// setup for villin
 
static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
static void setupAndGetForcesEnergyMultipoleVillin( AmoebaMultipoleForce::PolarizationType polarizationType,
int includeCavityTerm, std::vector<Vec3>& forces, double& energy, FILE* log ){
 
// beginning of Multipole setup
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaHarmonicAngleForce.cpp
View file @ 1db349e5
......@@ -34,7 +34,6 @@
*/
#include "openmm/internal/AssertionUtilities.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
......@@ -44,6 +43,8 @@
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-5;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaHarmonicBondForce.cpp
View file @ 1db349e5
......@@ -35,7 +35,6 @@
#include "openmm/internal/AssertionUtilities.h"
#include "CudaPlatform.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
......@@ -45,6 +44,8 @@
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-5;
static void computeAmoebaHarmonicBondForce(int bondIndex, std::vector<Vec3>& positions, AmoebaHarmonicBondForce& amoebaHarmonicBondForce,
std::vector<Vec3>& forces, double* energy ) {
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaHarmonicInPlaneAngleForce.cpp
View file @ 1db349e5
......@@ -34,7 +34,6 @@
*/
#include "openmm/internal/AssertionUtilities.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
......@@ -44,6 +43,8 @@
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-5;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaMultipoleForce.cpp
View file @ 1db349e5
......@@ -36,7 +36,6 @@
#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"
......@@ -53,10 +52,12 @@
using namespace OpenMM;
const double TOL = 1e-4;
extern "C" void registerAmoebaCudaKernelFactories();
// setup for 2 ammonia molecules
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::AmoebaNonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
static void setupAndGetForcesEnergyMultipoleAmmonia( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::vector<Vec3>& forces, double& energy, FILE* log ){
// beginning of Multipole setup
......@@ -498,8 +499,8 @@ static void testMultipoleAmmoniaMutualPolarization( FILE* log ) {
// setup for box of 4 water molecules -- used to test PME
static void setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::AmoebaNonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
static void setupAndGetForcesEnergyMultipoleWater( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::vector<Vec3>& forces,
double& energy, FILE* log ){
......@@ -938,8 +939,8 @@ static void testQuadrupoleValidation( FILE* log ){
// this method does too much; I tried passing the context ptr back to
// the tests methods, but the tests would seg fault w/ a bad_alloc error
static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::AmoebaNonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
static void setupAndGetForcesEnergyMultipoleIonsAndWater( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::string testName,
std::vector<Vec3>& forces, double& energy, FILE* log ){
......@@ -1204,8 +1205,8 @@ static void testMultipoleIonsAndWaterPMEMutualPolarization( FILE* log ) {
// setup for box of 216 water molecules -- used to test PME
static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::AmoebaNonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::AmoebaPolarizationType polarizationType,
static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::NonbondedMethod nonbondedMethod,
AmoebaMultipoleForce::PolarizationType polarizationType,
double cutoff, int inputPmeGridDimension, std::string& testName,
std::vector<Vec3>& forces, double& energy,
std::vector< double >& outputMultipoleMoments,
......@@ -1993,8 +1994,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::Am
context.setPositions(positions);
if( testName == "testSystemMultipoleMoments" ){
Vec3 origin( 0.0, 0.0, 0.0 );
amoebaMultipoleForce->getSystemMultipoleMoments( origin, context, outputMultipoleMoments );
amoebaMultipoleForce->getSystemMultipoleMoments(context, outputMultipoleMoments );
} else if( testName == "testMultipoleGridPotential" ){
amoebaMultipoleForce->getElectrostaticPotential( inputGrid, context, outputGridPotential );
} else {
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaOutOfPlaneBendForce.cpp
View file @ 1db349e5
......@@ -34,7 +34,6 @@
*/
#include "openmm/internal/AssertionUtilities.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
......@@ -44,6 +43,8 @@
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-3;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaPME.cpp deleted 100644 → 0
View file @ a919f305
/* -------------------------------------------------------------------------- *
* 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-2010 Stanford University and the Authors. *
* Authors: Mark Friedrichs *
* 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 AmoebaMultipoleForce.
*/
#include "openmm/internal/AssertionUtilities.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "CudaPlatform.h"
#include "openmm/AmoebaMultipoleForce.h"
#include "openmm/System.h"
#include "openmm/VerletIntegrator.h"
#include "../src/SimTKUtilities/SimTKOpenMMRealType.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
using namespace std;
const double TOL = 1e-5;
void testPMEWater() {
Platform& platform = Platform::getPlatformByName("CUDA");
System system;
system.addParticle(16);
system.addParticle(1);
system.addParticle(1);
VerletIntegrator integrator(0.01);
AmoebaMultipoleForce* mp = new AmoebaMultipoleForce();
mp->setNonbondedMethod(AmoebaMultipoleForce::PME);
vector<double> dipole(3, 0.0);
dipole[2] = 7.556121361e-2;
vector<double> quadrupole(9, 0.0);
quadrupole[0] = 3.540307211e-2;
quadrupole[4] = -3.902570771e-2;
quadrupole[8] = 3.622635596e-3;
double damp = 9.707801995e-01*sqrt(0.1);
double polarity = 0.837*0.001;
mp->addParticle(-0.51966, dipole, quadrupole, 1, 1, 2, -1, 0.39, damp, polarity);
dipole[0] = -2.042094848e-2;
dipole[2] = -3.078753000e-2;
quadrupole[0] = -3.428482490e-3;
quadrupole[2] = -1.894859639e-4;
quadrupole[4] = -1.002408752e-2;
quadrupole[6] = -1.894859639e-4;
quadrupole[8] = 1.345257001e-2;
damp = 8.897068742e-01*sqrt(0.1);
polarity = 0.496*0.001;
mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 2, -1, 0.39, damp, polarity);
mp->addParticle(0.25983, dipole, quadrupole, 0, 0, 1, -1, 0.39, damp, polarity);
mp->setCutoffDistance(1.0);
std::vector<int> intVector;
intVector.push_back( 0 );
intVector.push_back( 1 );
intVector.push_back( 2 );
mp->setCovalentMap( 0, AmoebaMultipoleForce::PolarizationCovalent11, intVector );
mp->setCovalentMap( 1, AmoebaMultipoleForce::PolarizationCovalent11, intVector );
mp->setCovalentMap( 2, AmoebaMultipoleForce::PolarizationCovalent11, intVector );
intVector.resize(0); intVector.push_back( 1 ); intVector.push_back( 2 );
mp->setCovalentMap( 0, AmoebaMultipoleForce::Covalent12, intVector );
intVector.resize(0); intVector.push_back( 0 ); intVector.push_back( 2 );
mp->setCovalentMap( 1, AmoebaMultipoleForce::Covalent12, intVector );
intVector.resize(0); intVector.push_back( 0 ); intVector.push_back( 1 );
mp->setCovalentMap( 2, AmoebaMultipoleForce::Covalent12, intVector );
mp->setEwaldErrorTolerance(TOL);
system.setDefaultPeriodicBoxVectors(Vec3(2, 0, 0), Vec3(0, 2, 0), Vec3(0, 0, 2));
system.addForce(mp);
Context context(system, integrator, platform);
vector<Vec3> positions(3);
const double angle = 109.47*M_PI/180;
const double dOH = 0.1;
positions[0] = Vec3();
positions[1] = Vec3(dOH, 0, 0);
positions[2] = Vec3(dOH*std::cos(angle), dOH*std::sin(angle), 0);
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
const vector<Vec3>& forces = state.getForces();
#ifdef AMOEBA_DEBUG
(void) fprintf( stderr, "PME forces\n" );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( stderr, "%6u [%14.7e %14.7e %14.7e]\n", ii,
forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( stderr );
#endif
}
int main() {
try {
registerAmoebaCudaKernelFactories();
testPMEWater();
} catch(const exception& e) {
cout << "exception: " << e.what() << endl;
return 1;
}
std::cout << "Done" << std::endl;
return 0;
}
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaPiTorsionForce.cpp
View file @ 1db349e5
......@@ -34,7 +34,6 @@
*/
#include "openmm/internal/AssertionUtilities.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
......@@ -44,6 +43,8 @@
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-5;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaStretchBendForce.cpp
View file @ 1db349e5
......@@ -34,7 +34,6 @@
*/
#include "openmm/internal/AssertionUtilities.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
......@@ -44,9 +43,12 @@
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-5;
#define PI_M 3.141592653589
#define RADIAN 57.29577951308
const double DegreesToRadians = PI_M/180.0;
/* ---------------------------------------------------------------------------------------
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaTorsionTorsionForce.cpp
View file @ 1db349e5
......@@ -37,13 +37,15 @@
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "openmm/System.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/AmoebaTorsionTorsionForce.h"
#include "openmm/LangevinIntegrator.h"
#include <iostream>
#include <vector>
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-4;
TorsionTorsionGrid& getTorsionGrid( int gridIndex ) {
......
plugins/amoeba/platforms/cuda2/tests/TestCudaAmoebaVdwForce.cpp
View file @ 1db349e5
......@@ -36,7 +36,6 @@
#include "openmm/internal/AssertionUtilities.h"
#include "openmm/Context.h"
#include "OpenMMAmoeba.h"
#include "AmoebaTinkerParameterFile.h"
#include "openmm/System.h"
#include "openmm/AmoebaVdwForce.h"
#include "openmm/LangevinIntegrator.h"
......@@ -51,6 +50,9 @@
using namespace OpenMM;
extern "C" void registerAmoebaCudaKernelFactories();
const double TOL = 1e-4;
void testVdw( FILE* log ) {
......@@ -63,7 +65,6 @@ 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;
double mass, sigma, epsilon, reduction;
......@@ -94,7 +95,7 @@ void testVdw( FILE* log ) {
exclusions.push_back ( 5 );
}
system.addParticle(mass);
amoebaVdwForce->addParticle( indexIV, classIndex, sigma, epsilon, reduction );
amoebaVdwForce->addParticle( indexIV, sigma, epsilon, reduction );
amoebaVdwForce->setParticleExclusions( ii, exclusions );
}
LangevinIntegrator integrator(0.0, 0.1, 0.01);
......@@ -136,12 +137,11 @@ void testVdw( FILE* log ) {
}
for( int ii = 0; ii < amoebaVdwForce->getNumParticles(); ii++ ){
int indexIV;
int classIndex;
double sigma, epsilon, reduction;
amoebaVdwForce->getParticleParameters( ii, indexIV, classIndex, sigma, epsilon, reduction );
amoebaVdwForce->getParticleParameters( ii, indexIV, sigma, epsilon, reduction );
sigma *= AngstromToNm;
epsilon *= CalToJoule;
amoebaVdwForce->setParticleParameters( ii, indexIV, classIndex, sigma, epsilon, reduction );
amoebaVdwForce->setParticleParameters( ii, indexIV, sigma, epsilon, reduction );
}
platformName = "CUDA";
Context context(system, integrator, Platform::getPlatformByName( platformName ) );
......@@ -186,46 +186,44 @@ void setupAndGetForcesEnergyVdwAmmonia( const std::string& sigmaCombiningRule, c
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 );
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
amoebaVdwForce->setUseDispersionCorrection( 1 );
} else {
amoebaVdwForce->setPBC( 0 );
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::NoCutoff);
amoebaVdwForce->setUseDispersionCorrection( 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 );
amoebaVdwForce->addParticle( 0, 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 );
amoebaVdwForce->addParticle( 0, 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 );
amoebaVdwForce->addParticle( 0, 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 );
amoebaVdwForce->addParticle( 0, 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 );
amoebaVdwForce->addParticle( 4, 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 );
amoebaVdwForce->addParticle( 4, 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 );
amoebaVdwForce->addParticle( 4, 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 );
amoebaVdwForce->addParticle( 4, 1.3500000e-01, 8.3680000e-02, 9.1000000e-01 );
// ParticleExclusions
......@@ -463,7 +461,7 @@ void testVdwTaper( FILE* log ) {
std::string testName = "testVdwTaper";
int numberOfParticles = 8;
double boxDimension = -1.0;
double boxDimension = 50.0;
double cutoff = 0.25;
std::vector<Vec3> forces;
......@@ -532,33 +530,32 @@ void setupAndGetForcesEnergyVdwWater( const std::string& sigmaCombiningRule, con
int numberOfParticles = 648;
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 );
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::CutoffPeriodic);
amoebaVdwForce->setUseDispersionCorrection( includeVdwDispersionCorrection );
} else {
amoebaVdwForce->setPBC( 0 );
amoebaVdwForce->setNonbondedMethod(AmoebaVdwForce::NoCutoff);
amoebaVdwForce->setUseDispersionCorrection( 0 );
}
// addParticle: ivIndex, classIndex, radius, epsilon, reductionFactor
// addParticle: ivIndex, radius, epsilon, reductionFactor
int classIndex = 0;
for( unsigned int ii = 0; ii < numberOfParticles; ii += 3 ){
system.addParticle( 1.5995000e+01 );
amoebaVdwForce->addParticle( ii, 0, 1.7025000e-01, 4.6024000e-01, 0.0000000e+00 );
amoebaVdwForce->addParticle( ii, 1.7025000e-01, 4.6024000e-01, 0.0000000e+00 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( ii, 0, 1.3275000e-01, 5.6484000e-02, 9.1000000e-01 );
amoebaVdwForce->addParticle( ii, 1.3275000e-01, 5.6484000e-02, 9.1000000e-01 );
system.addParticle( 1.0080000e+00 );
amoebaVdwForce->addParticle( ii, 0, 1.3275000e-01, 5.6484000e-02, 9.1000000e-01 );
amoebaVdwForce->addParticle( ii, 1.3275000e-01, 5.6484000e-02, 9.1000000e-01 );
}
// exclusions
......
plugins/amoeba/platforms/cuda2/tests/TestCudaWcaDispersionForce.cpp
View file @ 1db349e5
......@@ -36,7 +36,6 @@
#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"
......@@ -53,6 +52,8 @@
using namespace OpenMM;
const double TOL = 1e-4;
extern "C" void registerAmoebaCudaKernelFactories();
void setupAndGetForcesEnergyWcaDispersionAmmonia( std::vector<Vec3>& forces, double& energy, FILE* log ){
// beginning of WcaDispersion setup
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp
View file @ 1db349e5
......@@ -53,10 +53,8 @@ extern "C" void initAmoebaReferenceKernels() {
AmoebaReferenceKernelFactory* factory = new AmoebaReferenceKernelFactory();
platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaUreyBradleyForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory);
platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(), factory);
......@@ -79,18 +77,12 @@ KernelImpl* AmoebaReferenceKernelFactory::createKernelImpl(std::string name, con
if (name == CalcAmoebaHarmonicBondForceKernel::Name())
return new ReferenceCalcAmoebaHarmonicBondForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaUreyBradleyForceKernel::Name())
return new ReferenceCalcAmoebaUreyBradleyForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaHarmonicAngleForceKernel::Name())
return new ReferenceCalcAmoebaHarmonicAngleForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaHarmonicInPlaneAngleForceKernel::Name())
return new ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaTorsionForceKernel::Name())
return new ReferenceCalcAmoebaTorsionForceKernel(name, platform, context.getSystem());
if (name == CalcAmoebaPiTorsionForceKernel::Name())
return new ReferenceCalcAmoebaPiTorsionForceKernel(name, platform, context.getSystem());
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.cpp
View file @ 1db349e5
......@@ -28,7 +28,6 @@
#include "AmoebaReferenceHarmonicBondForce.h"
#include "AmoebaReferenceHarmonicAngleForce.h"
#include "AmoebaReferenceHarmonicInPlaneAngleForce.h"
#include "AmoebaReferenceTorsionForce.h"
#include "AmoebaReferencePiTorsionForce.h"
#include "AmoebaReferenceStretchBendForce.h"
#include "AmoebaReferenceOutOfPlaneBendForce.h"
......@@ -38,7 +37,6 @@
#include "AmoebaReferenceWcaDispersionForce.h"
#include "openmm/internal/AmoebaTorsionTorsionForceImpl.h"
#include "openmm/internal/AmoebaWcaDispersionForceImpl.h"
#include "AmoebaReferenceUreyBradleyForce.h"
#include "ReferencePlatform.h"
#include "openmm/internal/ContextImpl.h"
#include "openmm/AmoebaMultipoleForce.h"
......@@ -111,41 +109,6 @@ double ReferenceCalcAmoebaHarmonicBondForceKernel::execute(ContextImpl& context,
// ***************************************************************************
ReferenceCalcAmoebaUreyBradleyForceKernel::ReferenceCalcAmoebaUreyBradleyForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaUreyBradleyForceKernel(name, platform), system(system) {
}
ReferenceCalcAmoebaUreyBradleyForceKernel::~ReferenceCalcAmoebaUreyBradleyForceKernel() {
}
void ReferenceCalcAmoebaUreyBradleyForceKernel::initialize(const System& system, const AmoebaUreyBradleyForce& force) {
numIxns = force.getNumInteractions();
for( int ii = 0; ii < numIxns; ii++) {
int particle1Index, particle2Index;
double lengthValue, kValue;
force.getUreyBradleyParameters(ii, particle1Index, particle2Index, lengthValue, kValue );
particle1.push_back( particle1Index );
particle2.push_back( particle2Index );
length.push_back( static_cast<RealOpenMM>( lengthValue ) );
kQuadratic.push_back( static_cast<RealOpenMM>( kValue ) );
}
globalUreyBradleyCubic = static_cast<RealOpenMM>(force.getAmoebaGlobalUreyBradleyCubic());
globalUreyBradleyQuartic = static_cast<RealOpenMM>(force.getAmoebaGlobalUreyBradleyQuartic());
}
double ReferenceCalcAmoebaUreyBradleyForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
AmoebaReferenceUreyBradleyForce amoebaReferenceUreyBradleyForce;
RealOpenMM energy = amoebaReferenceUreyBradleyForce.calculateForceAndEnergy( numIxns, posData, particle1, particle2, length, kQuadratic,
globalUreyBradleyCubic, globalUreyBradleyQuartic,
forceData );
return static_cast<double>(energy);
}
ReferenceCalcAmoebaHarmonicAngleForceKernel::ReferenceCalcAmoebaHarmonicAngleForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaHarmonicAngleForceKernel(name, platform), system(system) {
}
......@@ -220,55 +183,6 @@ double ReferenceCalcAmoebaHarmonicInPlaneAngleForceKernel::execute(ContextImpl&
return static_cast<double>(energy);
}
ReferenceCalcAmoebaTorsionForceKernel::ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaTorsionForceKernel(name, platform), system(system) {
}
ReferenceCalcAmoebaTorsionForceKernel::~ReferenceCalcAmoebaTorsionForceKernel() {
}
void ReferenceCalcAmoebaTorsionForceKernel::initialize(const System& system, const AmoebaTorsionForce& force) {
numTorsions = force.getNumTorsions();
torsionParameters1.resize( numTorsions );
torsionParameters2.resize( numTorsions );
torsionParameters3.resize( numTorsions );
for (int ii = 0; ii < numTorsions; ii++) {
int particle1Index, particle2Index, particle3Index, particle4Index;
std::vector<double> torsionParameter1;
std::vector<double> torsionParameter2;
std::vector<double> torsionParameter3;
std::vector<RealOpenMM> torsionParameters1F(3);
std::vector<RealOpenMM> torsionParameters2F(3);
std::vector<RealOpenMM> torsionParameters3F(3);
force.getTorsionParameters(ii, particle1Index, particle2Index, particle3Index, particle4Index, torsionParameter1, torsionParameter2, torsionParameter3 );
particle1.push_back( particle1Index );
particle2.push_back( particle2Index );
particle3.push_back( particle3Index );
particle4.push_back( particle4Index );
torsionParameters1[ii].resize( torsionParameter1.size() );
torsionParameters2[ii].resize( torsionParameter2.size() );
torsionParameters3[ii].resize( torsionParameter3.size() );
for ( unsigned int jj = 0; jj < torsionParameter1.size(); jj++) {
torsionParameters1[ii][jj] = static_cast<RealOpenMM>(torsionParameter1[jj]);
torsionParameters2[ii][jj] = static_cast<RealOpenMM>(torsionParameter2[jj]);
torsionParameters3[ii][jj] = static_cast<RealOpenMM>(torsionParameter3[jj]);
}
}
}
double ReferenceCalcAmoebaTorsionForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
AmoebaReferenceTorsionForce amoebaReferenceTorsionForce;
RealOpenMM energy = amoebaReferenceTorsionForce.calculateForceAndEnergy( numTorsions, posData, particle1, particle2, particle3, particle4,
torsionParameters1, torsionParameters2, torsionParameters3, forceData );
return static_cast<double>(energy);
}
ReferenceCalcAmoebaPiTorsionForceKernel::ReferenceCalcAmoebaPiTorsionForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaPiTorsionForceKernel(name, platform), system(system) {
}
......@@ -568,8 +482,7 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextI
return;
}
void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, const Vec3& origin,
std::vector< double >& outputMultipoleMonents){
void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMonents){
return;
}
......@@ -618,7 +531,7 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextI
ReferenceCalcAmoebaVdwForceKernel::ReferenceCalcAmoebaVdwForceKernel(std::string name, const Platform& platform, System& system) :
CalcAmoebaVdwForceKernel(name, platform), system(system) {
useNeighborList = 0;
useCutoff = 0;
usePBC = 0;
cutoff = 1.0e+10;
......@@ -634,7 +547,6 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
numParticles = system.getNumParticles();
indexIVs.resize( numParticles );
indexClasses.resize( numParticles );
allExclusions.resize( numParticles );
sigmas.resize( numParticles );
epsilons.resize( numParticles );
......@@ -642,26 +554,25 @@ void ReferenceCalcAmoebaVdwForceKernel::initialize(const System& system, const A
for( int ii = 0; ii < numParticles; ii++ ){
int indexIV, indexClass;
int indexIV;
double sigma, epsilon, reduction;
std::vector<int> exclusions;
force.getParticleParameters( ii, indexIV, indexClass, sigma, epsilon, reduction );
force.getParticleParameters( ii, indexIV, 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();
useNeighborList = force.getUseNeighborList();
usePBC = force.getPBC();
useCutoff = (force.getNonbondedMethod() != AmoebaVdwForce::NoCutoff);
usePBC = (force.getNonbondedMethod() == AmoebaVdwForce::CutoffPeriodic);
cutoff = force.getCutoff();
}
......@@ -670,7 +581,7 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
vector<RealVec>& posData = extractPositions(context);
vector<RealVec>& forceData = extractForces(context);
AmoebaReferenceVdwForce vdwForce( sigmaCombiningRule, epsilonCombiningRule );
if( useNeighborList ){
if( useCutoff ){
vdwForce.setCutoff( cutoff );
if( usePBC ){
vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::CutoffPeriodic);
......@@ -680,7 +591,7 @@ double ReferenceCalcAmoebaVdwForceKernel::execute(ContextImpl& context, bool inc
} else {
vdwForce.setNonbondedMethod( AmoebaReferenceVdwForce::NoCutoff );
}
RealOpenMM energy = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, indexClasses, sigmas, epsilons, reductions, allExclusions, forceData);
RealOpenMM energy = vdwForce.calculateForceAndEnergy( numParticles, posData, indexIVs, sigmas, epsilons, reductions, allExclusions, forceData);
return static_cast<double>(energy);
}
......
plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h
View file @ 1db349e5
......@@ -70,42 +70,6 @@ private:
System& system;
};
/**
* This kernel is invoked by AmoebaUreyBradleyForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcAmoebaUreyBradleyForceKernel : public CalcAmoebaUreyBradleyForceKernel {
public:
ReferenceCalcAmoebaUreyBradleyForceKernel(std::string name,
const Platform& platform,
System& system);
~ReferenceCalcAmoebaUreyBradleyForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the AmoebaUreyBradleyForce this kernel will be used for
*/
void initialize(const System& system, const AmoebaUreyBradleyForce& 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 numIxns;
std::vector<int> particle1;
std::vector<int> particle2;
std::vector<RealOpenMM> length;
std::vector<RealOpenMM> kQuadratic;
RealOpenMM globalUreyBradleyCubic;
RealOpenMM globalUreyBradleyQuartic;
System& system;
};
/**
* This kernel is invoked by AmoebaHarmonicAngleForce to calculate the forces acting on the system and the energy of the system.
*/
......@@ -181,41 +145,6 @@ private:
System& system;
};
/**
* This kernel is invoked by AmoebaTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
class ReferenceCalcAmoebaTorsionForceKernel : public CalcAmoebaTorsionForceKernel {
public:
ReferenceCalcAmoebaTorsionForceKernel(std::string name, const Platform& platform, System& system);
~ReferenceCalcAmoebaTorsionForceKernel();
/**
* Initialize the kernel.
*
* @param system the System this kernel will be applied to
* @param force the AmoebaTorsionForce this kernel will be used for
*/
void initialize(const System& system, const AmoebaTorsionForce& 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 numTorsions;
std::vector<int> particle1;
std::vector<int> particle2;
std::vector<int> particle3;
std::vector<int> particle4;
std::vector< std::vector<RealOpenMM> > torsionParameters1;
std::vector< std::vector<RealOpenMM> > torsionParameters2;
std::vector< std::vector<RealOpenMM> > torsionParameters3;
System& system;
};
/**
* This kernel is invoked by AmoebaPiTorsionForce to calculate the forces acting on the system and the energy of the system.
*/
......@@ -398,7 +327,6 @@ public:
/**
* Get the system multipole moments
*
* @param origin origin
* @param context context
* @param outputMultipoleMonents (charge,
dipole_x, dipole_y, dipole_z,
......@@ -406,7 +334,7 @@ public:
quadrupole_yx, quadrupole_yy, quadrupole_yz,
quadrupole_zx, quadrupole_zy, quadrupole_zz )
*/
void getSystemMultipoleMoments(ContextImpl& context, const Vec3& origin, std::vector< double >& outputMultipoleMonents);
void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMonents);
private:
int numMultipoles;
......@@ -456,11 +384,10 @@ public:
double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
private:
int numParticles;
int useNeighborList;
int useCutoff;
int usePBC;
double cutoff;
std::vector<int> indexIVs;
std::vector<int> indexClasses;
std::vector< std::vector<int> > allExclusions;
std::vector<RealOpenMM> sigmas;
std::vector<RealOpenMM> epsilons;
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionForce.cpp deleted 100644 → 0
View file @ a919f305
/* 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 "AmoebaReferenceTorsionForce.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param positionAtomD Cartesian coordinates of atom D
@param torsion1 vector of torsion params for first index (amplitude, phase, fold)
@param torsion2 vector of torsion params for second index (amplitude, phase, fold)
@param torsion3 vector of torsion params for third index (amplitude, phase, fold)
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceTorsionForce::calculateTorsionIxn( const RealVec& positionAtomA, const RealVec& positionAtomB,
const RealVec& positionAtomC, const RealVec& positionAtomD,
const std::vector<RealOpenMM>& torsionParameters1,
const std::vector<RealOpenMM>& torsionParameters2,
const std::vector<RealOpenMM>& torsionParameters3,
RealVec* forces ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceTorsionForce::calculateHarmonicForce";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM two = 2.0;
static const RealOpenMM three = 3.0;
// ---------------------------------------------------------------------------------------
enum { BA, CB, DC, CA, DB, LastDeltaIndex };
std::vector<RealOpenMM> deltaR[LastDeltaIndex];
AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomB, deltaR[BA] );
AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomC, deltaR[CB] );
AmoebaReferenceForce::loadDeltaR( positionAtomC, positionAtomD, deltaR[DC] );
AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomC, deltaR[CA] );
AmoebaReferenceForce::loadDeltaR( positionAtomB, positionAtomD, deltaR[DB] );
enum { Xt, Xu, Xtu, LastXtIndex };
std::vector<RealOpenMM> crossProducts[LastXtIndex];
for( unsigned int ii = 0; ii < LastXtIndex; ii++ ){
crossProducts[ii].resize(3);
}
AmoebaReferenceForce::getCrossProduct( deltaR[BA], deltaR[CB], crossProducts[Xt] );
AmoebaReferenceForce::getCrossProduct( deltaR[CB], deltaR[DC], crossProducts[Xu] );
AmoebaReferenceForce::getCrossProduct( crossProducts[Xt], crossProducts[Xu], crossProducts[Xtu] );
RealOpenMM rT2 = AmoebaReferenceForce::getNormSquared3( crossProducts[Xt] );
RealOpenMM rU2 = AmoebaReferenceForce::getNormSquared3( crossProducts[Xu] );
RealOpenMM rTrU = SQRT( rT2*rU2 );
if( rTrU <= zero ){
return zero;
}
RealOpenMM rCB = AmoebaReferenceForce::getNorm3( deltaR[CB] );
// ---------------------------------------------------------------------------------------
// cos(w), cos(2w), cos(3w), ...
// sin(w), sin(2w), sin(3w), ...
RealOpenMM cosine = AmoebaReferenceForce::getDotProduct3( crossProducts[Xt], crossProducts[Xu] );
cosine /= rTrU;
RealOpenMM sine = AmoebaReferenceForce::getDotProduct3( deltaR[CB], crossProducts[Xtu] );
sine /= (rCB*rTrU);
RealOpenMM cosine2 = cosine*cosine - sine*sine;
RealOpenMM sine2 = two*cosine*sine;
RealOpenMM cosine3 = cosine*cosine2 - sine*sine2;
RealOpenMM sine3 = cosine*sine2 + sine*cosine2;
/*
RealOpenMM cosine4 = cosine*cosine3 - sine*sine3;
RealOpenMM sine4 = cosine*sine3 + sine*cosine3;
RealOpenMM cosine5 = cosine*cosine4 - sine*sine4;
RealOpenMM sine5 = cosine*sine4 + sine*cosine4;
RealOpenMM cosine6 = cosine*cosine5 - sine*sine5;
RealOpenMM sine6 = cosine*sine5 + sine*cosine5;
*/
// ---------------------------------------------------------------------------------------
// dEdPhi prefactor
RealOpenMM dEdPhi = torsionParameters1[0]* (cosine* sin( torsionParameters1[1] ) - sine* cos( torsionParameters1[1] ) );
dEdPhi += torsionParameters2[0]*two* (cosine2*sin( torsionParameters2[1] ) - sine2*cos( torsionParameters2[1] ) );
dEdPhi += torsionParameters3[0]*three*(cosine3*sin( torsionParameters3[1] ) - sine3*cos( torsionParameters3[1] ) );
// ---------------------------------------------------------------------------------------
// dEdtu[0] = dEdT
// dEdtu[1] = dEdU
// tempVector[0] == dEdA: dEdT x CB
// tempVector[1] == dEdB: (CA x dEdT) + (dEdU x DC)
// tempVector[2] == dEdC: (dEdT x BA) + (DB x dEdU)
// tempVector[3] == dEdD: (dEdU x CB)
std::vector<RealOpenMM> dEdtu[2];
std::vector<RealOpenMM> tempVector[6];
for( unsigned int ii = 0; ii < 2; ii++ ){
dEdtu[ii].resize(3);
}
for( unsigned int ii = 0; ii < 6; ii++ ){
tempVector[ii].resize(3);
}
// dEdT & dEdU
AmoebaReferenceForce::getCrossProduct( crossProducts[Xt], deltaR[CB], tempVector[0] );
AmoebaReferenceForce::getCrossProduct( crossProducts[Xu], deltaR[CB], tempVector[1] );
RealOpenMM norm[2] = { dEdPhi/(rT2*rCB ), -dEdPhi/(rU2*rCB ) };
dEdtu[0][0] = norm[0]*tempVector[0][0];
dEdtu[0][1] = norm[0]*tempVector[0][1];
dEdtu[0][2] = norm[0]*tempVector[0][2];
dEdtu[1][0] = norm[1]*tempVector[1][0];
dEdtu[1][1] = norm[1]*tempVector[1][1];
dEdtu[1][2] = norm[1]*tempVector[1][2];
// dEdA
AmoebaReferenceForce::getCrossProduct( dEdtu[0], deltaR[CB], tempVector[0] );
// dEdB
AmoebaReferenceForce::getCrossProduct( deltaR[CA], dEdtu[0], tempVector[4] );
AmoebaReferenceForce::getCrossProduct( dEdtu[1], deltaR[DC], tempVector[1] );
// dEdC
AmoebaReferenceForce::getCrossProduct( dEdtu[0], deltaR[BA], tempVector[5] );
AmoebaReferenceForce::getCrossProduct( deltaR[DB], dEdtu[1], tempVector[2] );
tempVector[1][0] += tempVector[4][0];
tempVector[2][0] += tempVector[5][0];
tempVector[1][1] += tempVector[4][1];
tempVector[2][1] += tempVector[5][1];
tempVector[1][2] += tempVector[4][2];
tempVector[2][2] += tempVector[5][2];
// dEdD
AmoebaReferenceForce::getCrossProduct( dEdtu[1], deltaR[CB], tempVector[3] );
// ---------------------------------------------------------------------------------------
// forces
forces[0][0] = tempVector[0][0];
forces[0][1] = tempVector[0][1];
forces[0][2] = tempVector[0][2];
forces[1][0] = tempVector[1][0];
forces[1][1] = tempVector[1][1];
forces[1][2] = tempVector[1][2];
forces[2][0] = tempVector[2][0];
forces[2][1] = tempVector[2][1];
forces[2][2] = tempVector[2][2];
forces[3][0] = tempVector[3][0];
forces[3][1] = tempVector[3][1];
forces[3][2] = tempVector[3][2];
// ---------------------------------------------------------------------------------------
// calculate energy
RealOpenMM energy = torsionParameters1[0]*( one + cosine* cos( torsionParameters1[1] ) + sine *sin( torsionParameters1[1] ) );
energy += torsionParameters2[0]*( one + cosine2*cos( torsionParameters2[1] ) + sine2*sin( torsionParameters2[1] ) );
energy += torsionParameters3[0]*( one + cosine3*cos( torsionParameters3[1] ) + sine3*sin( torsionParameters3[1] ) );
return energy;
}
RealOpenMM AmoebaReferenceTorsionForce::calculateForceAndEnergy( int numTorsions, vector<RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector< std::vector<RealOpenMM> >& torsionParameters1,
const std::vector< std::vector<RealOpenMM> >& torsionParameters2,
const std::vector< std::vector<RealOpenMM> >& torsionParameters3,
vector<RealVec>& forceData ) const {
RealOpenMM energy = 0.0;
for (unsigned int ii = 0; ii < static_cast<unsigned int>(numTorsions); ii++) {
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
int particle3Index = particle3[ii];
int particle4Index = particle4[ii];
RealVec forces[4];
energy += calculateTorsionIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
torsionParameters1[ii], torsionParameters2[ii], torsionParameters3[ii], forces );
// accumulate forces
for( int jj = 0; jj < 3; jj++ ){
forceData[particle1Index][jj] -= forces[0][jj];
forceData[particle2Index][jj] -= forces[1][jj];
forceData[particle3Index][jj] -= forces[2][jj];
forceData[particle4Index][jj] -= forces[3][jj];
}
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceTorsionForce.h deleted 100644 → 0
View file @ a919f305
/* 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 __AmoebaReferenceTorsionForce_H__
#define __AmoebaReferenceTorsionForce_H__
#include "SimTKUtilities/RealVec.h"
#include <vector>
// ---------------------------------------------------------------------------------------
class AmoebaReferenceTorsionForce {
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceTorsionForce( ){};
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~AmoebaReferenceTorsionForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba torsion ixns (force and energy)
@param numTorsions number of torsions
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param particle3 particle 3 indices
@param particle4 particle 4 indices
@param torsionParameters1 first index torsion parameters (amplitude, phase, fold)
@param torsionParameters2 second index torsion parameters (amplitude, phase, fold)
@param torsionParameters3 third index torsion parameters (amplitude, phase, fold)
@param forces output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numTorsions, std::vector<OpenMM::RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<int>& particle3,
const std::vector<int>& particle4,
const std::vector< std::vector<RealOpenMM> >& torsionParameters1,
const std::vector< std::vector<RealOpenMM> >& torsionParameters2,
const std::vector< std::vector<RealOpenMM> >& torsionParameters3,
std::vector<OpenMM::RealVec>& forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Amoeba harmonic angle ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param positionAtomC Cartesian coordinates of atom C
@param positionAtomD Cartesian coordinates of atom D
@param torsion1 vector of torsion params for first index (amplitude, phase, fold)
@param torsion2 vector of torsion params for second index (amplitude, phase, fold)
@param torsion3 vector of torsion params for third index (amplitude, phase, fold)
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateTorsionIxn( const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
const OpenMM::RealVec& positionAtomC, const OpenMM::RealVec& positionAtomD,
const std::vector<RealOpenMM>& torsionParameters1,
const std::vector<RealOpenMM>& torsionParameters2,
const std::vector<RealOpenMM>& torsionParameters3,
OpenMM::RealVec* forces ) const;
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceTorsionForce___
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceUreyBradleyForce.cpp deleted 100644 → 0
View file @ a919f305
/* 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 "AmoebaReferenceUreyBradleyForce.h"
#include "AmoebaReferenceForce.h"
using std::vector;
using OpenMM::RealVec;
/**---------------------------------------------------------------------------------------
Calculate Amoeba UB ixn (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param idealLength ideal length
@param kForceConstant force constant
@param cubicK cubic force parameter
@param quarticK quartic force parameter
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM AmoebaReferenceUreyBradleyForce::calculateUreyBradleyIxn( const RealVec& positionAtomA, const RealVec& positionAtomB,
RealOpenMM idealLength, RealOpenMM kForceConstant,
RealOpenMM cubicK, RealOpenMM quarticK,
RealVec* forces ) const {
// ---------------------------------------------------------------------------------------
//static const std::string methodName = "AmoebaReferenceUreyBradleyForce::calculateHarmonicForce";
static const RealOpenMM zero = 0.0;
static const RealOpenMM one = 1.0;
static const RealOpenMM onePt5 = 1.5;
static const RealOpenMM two = 2.0;
// ---------------------------------------------------------------------------------------
// get deltaR, R2, and R between 2 atoms
std::vector<RealOpenMM> deltaR;
AmoebaReferenceForce::loadDeltaR( positionAtomA, positionAtomB, deltaR );
RealOpenMM r = AmoebaReferenceForce::getNorm3( deltaR );
// deltaIdeal = r - r_0
RealOpenMM deltaIdeal = r - idealLength;
RealOpenMM deltaIdeal2 = deltaIdeal*deltaIdeal;
RealOpenMM dEdR = (one + onePt5*cubicK*deltaIdeal + two*quarticK*deltaIdeal2);
dEdR *= two*kForceConstant*deltaIdeal;
dEdR = r > zero ? (dEdR/r) : zero;
forces[0][0] = dEdR*deltaR[0];
forces[0][1] = dEdR*deltaR[1];
forces[0][2] = dEdR*deltaR[2];
forces[1][0] = dEdR*deltaR[0];
forces[1][1] = dEdR*deltaR[1];
forces[1][2] = dEdR*deltaR[2];
RealOpenMM energy = kForceConstant*deltaIdeal2*( one + cubicK*deltaIdeal + quarticK*deltaIdeal2 );
return energy;
}
RealOpenMM AmoebaReferenceUreyBradleyForce::calculateForceAndEnergy( int numIxns,
vector<RealVec>& particlePositions,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<RealOpenMM>& length,
const std::vector<RealOpenMM>& kQuadratic,
RealOpenMM globalUreyBradleyCubic,
RealOpenMM globalUreyBradleyQuartic,
vector<RealVec>& forceData ) const {
RealOpenMM energy = 0.0;
for( int ii = 0; ii < numIxns; ii++ ){
int particle1Index = particle1[ii];
int particle2Index = particle2[ii];
RealOpenMM idealLength = length[ii];
RealOpenMM kForceConstant = kQuadratic[ii];
RealVec forces[2];
energy += calculateUreyBradleyIxn( particlePositions[particle1Index], particlePositions[particle2Index],
idealLength, kForceConstant, globalUreyBradleyCubic, globalUreyBradleyQuartic,
forces );
// accumulate forces
for( int jj = 0; jj < 3; jj++ ){
forceData[particle1Index][jj] += forces[0][jj];
forceData[particle2Index][jj] -= forces[1][jj];
}
}
return energy;
}
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceUreyBradleyForce.h deleted 100644 → 0
View file @ a919f305
/* 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 __AmoebaReferenceUreyBradleyForce_H__
#define __AmoebaReferenceUreyBradleyForce_H__
#include "SimTKUtilities/RealVec.h"
#include <vector>
// ---------------------------------------------------------------------------------------
class AmoebaReferenceUreyBradleyForce {
public:
/**---------------------------------------------------------------------------------------
Constructor
--------------------------------------------------------------------------------------- */
AmoebaReferenceUreyBradleyForce( ){};
/**---------------------------------------------------------------------------------------
Destructor
--------------------------------------------------------------------------------------- */
~AmoebaReferenceUreyBradleyForce( ){};
/**---------------------------------------------------------------------------------------
Calculate Amoeba Urey-Bradley ixns (force and energy)
@param numIxns number of ixns
@param posData particle positions
@param particle1 particle 1 indices
@param particle2 particle 2 indices
@param length ideal length
@param kForce force constant
@param cubic cubic force parameter
@param quartic quartic force parameter
@param forceDate output force vector
@return total energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateForceAndEnergy( int numIxns, std::vector<OpenMM::RealVec>& posData,
const std::vector<int>& particle1,
const std::vector<int>& particle2,
const std::vector<RealOpenMM>& length,
const std::vector<RealOpenMM>& kForce,
RealOpenMM cubic, RealOpenMM quartic,
std::vector<OpenMM::RealVec>& forceData ) const;
private:
/**---------------------------------------------------------------------------------------
Calculate Amoeba Urey-Bradley ixns (force and energy)
@param positionAtomA Cartesian coordinates of atom A
@param positionAtomB Cartesian coordinates of atom B
@param length ideal length
@param kForce force constant
@param cubic cubic force parameter
@param quartic quartic force parameter
@param forces force vector
@return energy
--------------------------------------------------------------------------------------- */
RealOpenMM calculateUreyBradleyIxn( const OpenMM::RealVec& positionAtomA, const OpenMM::RealVec& positionAtomB,
RealOpenMM idealLength, RealOpenMM kForce,
RealOpenMM cubic, RealOpenMM quartic,
OpenMM::RealVec* forces ) const;
};
// ---------------------------------------------------------------------------------------
#endif // _AmoebaReferenceUreyBradleyForce___
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp
View file @ 1db349e5
......@@ -210,7 +210,6 @@ RealOpenMM AmoebaReferenceVdwForce::calculatePairIxn( RealOpenMM combindedSigma,
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numParticles,
const vector<RealVec>& 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,
......@@ -293,7 +292,6 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyNoCutoff( int numPart
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyApplyCutoff( int numParticles,
const vector<RealVec>& 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,
......@@ -376,7 +374,6 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergyApplyCutoff( int numP
RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
const vector<RealVec>& 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,
......@@ -386,11 +383,11 @@ RealOpenMM AmoebaReferenceVdwForce::calculateForceAndEnergy( int numParticles,
if( getNonbondedMethod() == NoCutoff ){
return calculateForceAndEnergyNoCutoff( numParticles, particlePositions,
indexIVs, indexClasses, sigmas, epsilons,
indexIVs, sigmas, epsilons,
reductions, allExclusions, forces );
} else {
return calculateForceAndEnergyApplyCutoff( numParticles, particlePositions,
indexIVs, indexClasses, sigmas, epsilons,
indexIVs, sigmas, epsilons,
reductions, allExclusions, forces );
}
......
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