Fixed incorrect specification of torsion-torsion terms in AMOEBA 2013 force field

plugins/amoeba/openmmapi/src/AmoebaTorsionTorsionForce.cpp

@@ -86,42 +86,55 @@ ForceImpl* AmoebaTorsionTorsionForce::createImpl() const {
 }
 
 AmoebaTorsionTorsionForce::TorsionTorsionGridInfo::TorsionTorsionGridInfo(const TorsionTorsionGrid& grid) {
-    _grid = 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;
+    _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());
+    _size[0]        = static_cast<int>(_grid.size());
+    _size[1]        = static_cast<int>(_grid[0].size());
     
-    if (grid[0][0].size() == 3) {
+    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])
+            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])
+            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];
+            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];
+                t[j] = _grid[j][i][2];
             if (periodic)
                 SplineFitter::createPeriodicSpline(x, t, deriv);
             else
@@ -136,7 +149,7 @@ AmoebaTorsionTorsionForce::TorsionTorsionGridInfo::TorsionTorsionGridInfo(const
         deriv.resize(ysize);
         for (int i = 0; i < xsize; i++) {
             for (int j = 0; j < ysize; j++)
-                t[j] = grid[i][j][2];
+                t[j] = _grid[i][j][2];
             if (periodic)
                 SplineFitter::createPeriodicSpline(y, t, deriv);
             else

wrappers/python/simtk/openmm/app/data/amoeba2013_gk.xml

 <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

@@ -2505,7 +2505,7 @@ class AmoebaTorsionTorsionGenerator:
                 gridRow.append(float(gridEntry.attrib['angle1']))
                 gridRow.append(float(gridEntry.attrib['angle2']))
                 gridRow.append(float(gridEntry.attrib['f']))
-                if 'fx' in gridEnergy.attrib:
+                if 'fx' in gridEntry.attrib:
                     gridRow.append(float(gridEntry.attrib['fx']))
                     gridRow.append(float(gridEntry.attrib['fy']))
                     gridRow.append(float(gridEntry.attrib['fxy']))