Merge pull request #455 from peastman/master

Implemented AMOEBA 2013 force field

docs/usersguide/application.rst

@@ -660,7 +660,9 @@ amber99sbildn.xml  AMBER99SB plus improved side chain torsions\ :cite:`Lindorff-
 amber99sbnmr.xml   AMBER99SB with modifications to fit NMR data\ :cite:`Li2010`
 amber03.xml        AMBER03\ :cite:`Duan2003`
 amber10.xml        AMBER10
-amoeba2009.xml     AMOEBA\ :cite:`Ren2002`                                                       
+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
-=================  ==============================================================================================
+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

@@ -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

@@ -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

@@ -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

@@ -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/forcefield.py

@@ -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']))