Skip to content

GitLab

Commit 7981628b authored by Peter Eastman's avatar Peter Eastman
Browse files

Merge branch 'peastman-amoeba'

parents cafd3d74 7eb5a10f
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 7321 additions and 89 deletions
docs/usersguide/application.rst
View file @ 7981628b
......@@ -654,13 +654,15 @@ For the main force field, OpenMM provides the following options:
================= ================================================================================
File Force Field
================= ================================================================================
amber96.xml AMBER96\ :cite:`Kollman1997`
amber96.xml AMBER96\ :cite:`Kollman1997`
amber99sb.xml AMBER99\ :cite:`Wang2000` with modified backbone torsions\ :cite:`Hornak2006`
amber99sbildn.xml AMBER99SB plus improved side chain torsions\ :cite:`Lindorff-Larsen2010`
amber99sbnmr.xml AMBER99SB with modifications to fit NMR data\ :cite:`Li2010`
amber03.xml AMBER03\ :cite:`Duan2003`
amber10.xml AMBER10
amoeba2009.xml AMOEBA\ :cite:`Ren2002`
amber99sbildn.xml AMBER99SB plus improved side chain torsions\ :cite:`Lindorff-Larsen2010`
amber99sbnmr.xml AMBER99SB with modifications to fit NMR data\ :cite:`Li2010`
amber03.xml AMBER03\ :cite:`Duan2003`
amber10.xml AMBER10
amoeba2009.xml AMOEBA 2009\ :cite:`Ren2002`. This force field is deprecated. It is
recommended to use AMOEBA 2013 instead.
amoeba2013.xml AMOEBA 2013\ :cite:`Shi2013`
================= ================================================================================
......@@ -692,15 +694,16 @@ the following files:
.. tabularcolumns:: |l|L|
================= ==============================================================================================
================= =================================================================================================
File Implicit Solvation Model
================= ==============================================================================================
amber96_obc.xml GBSA-OBC solvation model\ :cite:`Onufriev2004` for use with AMBER96 force field
amber99_obc.xml GBSA-OBC solvation model for use with AMBER99 force fields
amber03_obc.xml GBSA-OBC solvation model for use with AMBER03 force field
amber10_obc.xml GBSA-OBC solvation model for use with AMBER10 force field
amoeba2009_gk.xml Generalized Kirkwood solvation model\ :cite:`Schnieders2007` for use with AMOEBA force field
================= ==============================================================================================
================= =================================================================================================
amber96_obc.xml GBSA-OBC solvation model\ :cite:`Onufriev2004` for use with AMBER96 force field
amber99_obc.xml GBSA-OBC solvation model for use with AMBER99 force fields
amber03_obc.xml GBSA-OBC solvation model for use with AMBER03 force field
amber10_obc.xml GBSA-OBC solvation model for use with AMBER10 force field
amoeba2009_gk.xml Generalized Kirkwood solvation model\ :cite:`Schnieders2007` for use with AMOEBA 2009 force field
amoeba2013_gk.xml Generalized Kirkwood solvation model for use with AMOEBA 2013 force field
================= =================================================================================================
For example, to use the GBSA-OBC solvation model with the Amber99SB force field,
......
docs/usersguide/references.bib
View file @ 7981628b
......@@ -336,10 +336,20 @@
type = {Journal Article}
}
@article{Shi2013
author = {Shi, Yue and Xia, Zhen and Zhang, Jiajing and Best, Robert and Wu, Chuanjie and Ponder, Jay W. and Ren, Pengyu},
title = {Polarizable Atomic Multipole-Based AMOEBA Force Field for Proteins},
journal = {Journal of Chemical Theory and Computation},
volume = {9},
number = {9},
pages = {4046-4063},
year = {2013},
type = {Journal Article}
}
@article{Shirts2007
author = {Shirts, Michael R. and Mobley, David L. and Chodera, John D. and Pande, Vijay S.},
title = {Accurate and Efficient Corrections for Missing Dispersion Interactions in Molecular
Simulations},
title = {Accurate and Efficient Corrections for Missing Dispersion Interactions in Molecular Simulations},
journal = {Journal of Physical Chemistry B},
volume = {111},
pages = {13052-13063},
......
plugins/amoeba/openmmapi/include/openmm/AmoebaTorsionTorsionForce.h
View file @ 7981628b
......@@ -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 @ 7981628b
......@@ -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 @ 7981628b
......@@ -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/data/amoeba2013.xml 0 → 100644
View file @ 7981628b
This diff is collapsed.
wrappers/python/simtk/openmm/app/data/amoeba2013_gk.xml 0 → 100644
View file @ 7981628b
This diff is collapsed.
wrappers/python/simtk/openmm/app/forcefield.py
View file @ 7981628b
......@@ -2517,9 +2517,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 gridEntry.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