Skip to content

GitLab

Commit 14fb3330 authored by Lee-Ping's avatar Lee-Ping
Browse files

Merge branch 'amoeba' of github.com:peastman/openmm

parents c83f44dd bbe0ce70
Expand all Show whitespace changes
Inline Side-by-side
Showing with 2015 additions and 75 deletions
plugins/amoeba/openmmapi/include/openmm/AmoebaTorsionTorsionForce.h
View file @ 14fb3330
......@@ -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 @ 14fb3330
......@@ -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,102 @@ void AmoebaTorsionTorsionForce::setTorsionTorsionGrid(int index, const TorsionTo
ForceImpl* AmoebaTorsionTorsionForce::createImpl() const {
return new AmoebaTorsionTorsionForceImpl(*this);
}
AmoebaTorsionTorsionForce::TorsionTorsionGridInfo::TorsionTorsionGridInfo(const TorsionTorsionGrid& grid) {
if (grid[0][0][0] != grid[1][0][0])
_grid = grid;
else {
// We need to transpose the grid.
int xsize = grid[0].size();
int ysize = grid.size();
_grid.resize(xsize);
for (int i = 0; i < xsize; i++) {
_grid[i].resize(ysize);
for (int j = 0; j < ysize; j++)
_grid[i][j] = grid[j][i];
}
}
_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 < ysize; 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 @ 14fb3330
......@@ -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,16 +2557,14 @@ 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);
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(6);
grids[ii][jj][kk].resize(elementCount);
}
}
int index = 0;
......@@ -2574,18 +2572,17 @@ static double grid[4][625][6] = {
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++ ){
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/data/amoeba2013_gk.xml
View file @ 14fb3330
<ForceField>
<Info>
<Source>amoebapro13.prm</Source>
<DateGenerated>2014-03-11</DateGenerated>
<DateGenerated>2014-05-08</DateGenerated>
<Reference>Yue Shi, Zhen Xia, Jiajing Zhang, Robert Best, Chuanjie Wu, Jay W. Ponder, and Pengyu Ren. Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins. Journal of Chemical Theory and Computation, 9(9):4046–4063, 2013.</Reference>
</Info>
<AmoebaGeneralizedKirkwoodForce solventDielectric="78.3" soluteDielectric="1.0" includeCavityTerm="1" probeRadius="0.14" surfaceAreaFactor="-170.351730663">
......
wrappers/python/simtk/openmm/app/forcefield.py
View file @ 14fb3330
......@@ -2517,6 +2517,7 @@ class AmoebaTorsionTorsionGenerator:
gridRow.append(float(gridEntry.attrib['angle1']))
gridRow.append(float(gridEntry.attrib['angle2']))
gridRow.append(float(gridEntry.attrib['f']))
if 'fx' in gridEntry.attrib:
gridRow.append(float(gridEntry.attrib['fx']))
gridRow.append(float(gridEntry.attrib['fy']))
gridRow.append(float(gridEntry.attrib['fxy']))
......
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