Skip to content

GitLab

Commit 25e3852a authored by Peter Eastman's avatar Peter Eastman
Browse files

AmoebaTorsionTorsionForce can calculate derivatives of the grid points automatically

parent 075e56ba
Hide whitespace changes
Inline Side-by-side
Showing with 125 additions and 73 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaTorsionTorsionForce.h
View file @ 25e3852a
......@@ -126,13 +126,15 @@ public:
* Set the torsion-torsion grid at the specified index
*
* @param index the index of the torsion-torsion for which to get parameters
* @param grid grid
* @param grid either 3 or 6 values may be specified per grid point. If the derivatives
* are omitted, they are calculated automatically by fitting a 2D spline to
* the energies.
* grid[x][y][0] = x value
* grid[x][y][1] = y value
* grid[x][y][2] = function value
* grid[x][y][3] = dfdx value
* grid[x][y][4] = dfdy value
* grid[x][y][5] = dfd(xy) value
* grid[x][y][2] = energy
* grid[x][y][3] = dEdx value
* grid[x][y][4] = dEdy value
* grid[x][y][5] = dEd(xy) value
*/
void setTorsionTorsionGrid(int index, const std::vector<std::vector<std::vector<double> > >& grid);
......@@ -181,29 +183,7 @@ public:
_spacing[0] = _spacing[1] = 1.0;
}
TorsionTorsionGridInfo(const TorsionTorsionGrid& grid) {
_grid.resize(grid.size());
for(unsigned int kk = 0; kk < grid.size(); kk++){
_grid[kk].resize(grid[kk].size());
for(unsigned int jj = 0; jj < grid[kk].size(); jj++){
_grid[kk][jj].resize(grid[kk][jj].size());
for(unsigned int ii = 0; ii < grid[kk][jj].size(); ii++){
_grid[kk][jj][ii] = grid[kk][jj][ii];
}
}
}
_startValues[0] = _grid[0][0][0];
_startValues[1] = _grid[0][0][1];
_spacing[0] = static_cast<double>(_grid.size()-1)/360.0;
_spacing[1] = static_cast<double>(grid.size()-1)/360.0;
_size[0] = static_cast<int>(grid.size());
_size[1] = static_cast<int>(grid[0].size());
}
TorsionTorsionGridInfo(const TorsionTorsionGrid& grid);
const TorsionTorsionGrid& getTorsionTorsionGrid(void) const {
return _grid;
......
plugins/amoeba/openmmapi/src/AmoebaTorsionTorsionForce.cpp
View file @ 25e3852a
......@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. *
* *
* Portions copyright (c) 2008-2009 Stanford University and the Authors. *
* Portions copyright (c) 2008-2014 Stanford University and the Authors. *
* Authors: *
* Contributors: *
* *
......@@ -33,8 +33,10 @@
#include "openmm/OpenMMException.h"
#include "openmm/AmoebaTorsionTorsionForce.h"
#include "openmm/internal/AmoebaTorsionTorsionForceImpl.h"
#include "openmm/internal/SplineFitter.h"
using namespace OpenMM;
using namespace std;
AmoebaTorsionTorsionForce::AmoebaTorsionTorsionForce() {
}
......@@ -82,3 +84,89 @@ void AmoebaTorsionTorsionForce::setTorsionTorsionGrid(int index, const TorsionTo
ForceImpl* AmoebaTorsionTorsionForce::createImpl() const {
return new AmoebaTorsionTorsionForceImpl(*this);
}
AmoebaTorsionTorsionForce::TorsionTorsionGridInfo::TorsionTorsionGridInfo(const TorsionTorsionGrid& grid) {
_grid = grid;
_startValues[0] = _grid[0][0][0];
_startValues[1] = _grid[0][0][1];
_spacing[0] = static_cast<double>(_grid.size()-1)/360.0;
_spacing[1] = static_cast<double>(grid.size()-1)/360.0;
_size[0] = static_cast<int>(grid.size());
_size[1] = static_cast<int>(grid[0].size());
if (grid[0][0].size() == 3) {
// We need to compute the derivatives ourselves. First determine if the grid is periodic.
int xsize = _size[0];
int ysize = _size[1];
bool periodic = true;
for (int i = 0; i < xsize; i++)
if (grid[i][0][2] != grid[i][ysize-1][2])
periodic = false;
for (int i = 0; i < ysize; i++)
if (grid[0][i][2] != grid[xsize-1][i][2])
periodic = false;
// Compute derivatives with respect to the first angle.
vector<double> x(xsize), y(ysize);
for (int i = 0; i < xsize; i++)
x[i] = grid[i][0][0];
for (int i = 0; i < xsize; i++)
y[i] = grid[0][i][1];
vector<double> d1(xsize*ysize), d2(xsize*ysize), d12(xsize*ysize);
vector<double> t(xsize), deriv(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < xsize; j++)
t[j] = grid[j][i][2];
if (periodic)
SplineFitter::createPeriodicSpline(x, t, deriv);
else
SplineFitter::createNaturalSpline(x, t, deriv);
for (int j = 0; j < xsize; j++)
d1[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
}
// Compute derivatives with respect to the second angle.
t.resize(ysize);
deriv.resize(ysize);
for (int i = 0; i < xsize; i++) {
for (int j = 0; j < ysize; j++)
t[j] = grid[i][j][2];
if (periodic)
SplineFitter::createPeriodicSpline(y, t, deriv);
else
SplineFitter::createNaturalSpline(y, t, deriv);
for (int j = 0; j < ysize; j++)
d2[i+xsize*j] = SplineFitter::evaluateSplineDerivative(y, t, deriv, y[j]);
}
// Compute cross derivatives.
t.resize(xsize);
deriv.resize(xsize);
for (int i = 0; i < ysize; i++) {
for (int j = 0; j < xsize; j++)
t[j] = d2[j+xsize*i];
if (periodic)
SplineFitter::createPeriodicSpline(x, t, deriv);
else
SplineFitter::createNaturalSpline(x, t, deriv);
for (int j = 0; j < xsize; j++)
d12[j+xsize*i] = SplineFitter::evaluateSplineDerivative(x, t, deriv, x[j]);
}
// Add the derivatives to the grid.
for (int i = 0; i < xsize; i++)
for (int j = 0; j < ysize; j++) {
_grid[i][j].push_back(d1[i+xsize*j]);
_grid[i][j].push_back(d2[i+xsize*j]);
_grid[i][j].push_back(d12[i+xsize*j]);
}
}
}
plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaTorsionTorsionForce.cpp
View file @ 25e3852a
......@@ -49,7 +49,7 @@ extern "C" OPENMM_EXPORT void registerAmoebaReferenceKernelFactories();
const double TOL = 1e-4;
TorsionTorsionGrid& getTorsionGrid( int gridIndex ) {
TorsionTorsionGrid getTorsionGrid(int gridIndex, bool includeDerivs) {
static double grid[4][625][6] = {
{
......@@ -2557,35 +2557,32 @@ static double grid[4][625][6] = {
{ 165.0000, 180.0000, -0.182999000E+01, 0.377952854E-01, 0.233583295E-01, -0.109828932E-02 },
{ 180.0000, 180.0000, -0.146854000E+01, 0.491175487E-02, 0.195601580E-02, -0.163177030E-02 } } };
// static std::vector< std::vector< std::vector<double> > > TorsionTorsionGrid
static std::vector<TorsionTorsionGrid> grids;
if( grids.size() == 0 ){
grids.resize(4);
for( int ii = 0; ii < 4; ii++ ){
grids[ii].resize( 25 );
for( int jj = 0; jj < 25; jj++ ){
grids[ii][jj].resize(25);
for( int kk = 0; kk < 25; kk++ ){
grids[ii][jj][kk].resize(6);
}
int elementCount = (includeDerivs ? 6 : 3);
std::vector<TorsionTorsionGrid> grids(4);
for( int ii = 0; ii < 4; ii++ ){
grids[ii].resize( 25 );
for( int jj = 0; jj < 25; jj++ ){
grids[ii][jj].resize(25);
for( int kk = 0; kk < 25; kk++ ){
grids[ii][jj][kk].resize(elementCount);
}
int index = 0;
for( int jj = 0; jj < 25; jj++ ){
for( int kk = 0; kk < 25; kk++ ){
int jjIndex = static_cast<int>(((grid[ii][index][0] + 180.0)/15.0)+1.0e-05);
int kkIndex = static_cast<int>(((grid[ii][index][1] + 180.0)/15.0)+1.0e-05);
for( int ll = 0; ll < 6; ll++ ){
grids[ii][kk][jj][ll] = grid[ii][index][ll];
}
index++;
}
int index = 0;
for( int jj = 0; jj < 25; jj++ ){
for( int kk = 0; kk < 25; kk++ ){
int jjIndex = static_cast<int>(((grid[ii][index][0] + 180.0)/15.0)+1.0e-05);
int kkIndex = static_cast<int>(((grid[ii][index][1] + 180.0)/15.0)+1.0e-05);
for( int ll = 0; ll < elementCount; ll++ ){
grids[ii][kk][jj][ll] = grid[ii][index][ll];
}
index++;
}
}
}
return grids[gridIndex];
}
void testTorsionTorsion( FILE* log, int systemId ) {
void testTorsionTorsion(int systemId, bool includeDerivs) {
System system;
int numberOfParticles = 6;
......@@ -2645,11 +2642,11 @@ void testTorsionTorsion( FILE* log, int systemId ) {
expectedEnergy = -3.372536909E+00;
}
amoebaTorsionTorsionForce->addTorsionTorsion( 0, 1, 2, 3, 4, chiralCheckAtomIndex, 0 );
amoebaTorsionTorsionForce->setTorsionTorsionGrid( 0, getTorsionGrid( gridIndex ) );
amoebaTorsionTorsionForce->addTorsionTorsion(0, 1, 2, 3, 4, chiralCheckAtomIndex, 0);
amoebaTorsionTorsionForce->setTorsionTorsionGrid(0, getTorsionGrid(gridIndex, includeDerivs));
system.addForce(amoebaTorsionTorsionForce);
Context context(system, integrator, Platform::getPlatformByName( "Reference"));
Context context(system, integrator, Platform::getPlatformByName("Reference"));
context.setPositions(positions);
State state = context.getState(State::Forces | State::Energy);
......@@ -2662,18 +2659,6 @@ void testTorsionTorsion( FILE* log, int systemId ) {
forces[ii][2] *= conversion;
}
#ifdef AMOEBA_DEBUG
if( log ){
(void) fprintf( log, "computeAmoebaTorsionTorsionForces: expected energy=%14.7e %14.7e\n", expectedEnergy, state.getPotentialEnergy() );
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
(void) fprintf( log, "%6u [%14.7e %14.7e %14.7e] [%14.7e %14.7e %14.7e]\n", ii,
expectedForces[ii][0], expectedForces[ii][1], expectedForces[ii][2],
forces[ii][0], forces[ii][1], forces[ii][2] );
}
(void) fflush( log );
}
#endif
double tolerance = 1.0e-03;
for( unsigned int ii = 0; ii < forces.size(); ii++ ){
ASSERT_EQUAL_VEC( expectedForces[ii], forces[ii], tolerance );
......@@ -2687,10 +2672,8 @@ int main( int numberOfArguments, char* argv[] ) {
try {
std::cout << "TestReferenceAmoebaTorsionTorsionForce running test..." << std::endl;
registerAmoebaReferenceKernelFactories();
//registerAmoebaCudaKernelFactories();
FILE* log = NULL;
testTorsionTorsion( log, 1 );
testTorsionTorsion(1, true);
testTorsionTorsion(1, false);
} catch(const std::exception& e) {
std::cout << "exception: " << e.what() << std::endl;
std::cout << "FAIL - ERROR. Test failed." << std::endl;
......
wrappers/python/simtk/openmm/app/forcefield.py
View file @ 25e3852a
......@@ -2505,9 +2505,10 @@ class AmoebaTorsionTorsionGenerator:
gridRow.append(float(gridEntry.attrib['angle1']))
gridRow.append(float(gridEntry.attrib['angle2']))
gridRow.append(float(gridEntry.attrib['f']))
gridRow.append(float(gridEntry.attrib['fx']))
gridRow.append(float(gridEntry.attrib['fy']))
gridRow.append(float(gridEntry.attrib['fxy']))
if 'fx' in gridEnergy.attrib:
gridRow.append(float(gridEntry.attrib['fx']))
gridRow.append(float(gridEntry.attrib['fy']))
gridRow.append(float(gridEntry.attrib['fxy']))
gridCol.append(gridRow)
gridColIndex += 1
......
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