text: 
  /* Portions copyright (c) 2010-2016 Stanford University and Simbios.
   * Contributors: Peter Eastman
   *
   * Permission is hereby granted, free of charge, to any person obtaining
   * a copy of this software and associated documentation files (the
   * "Software"), to deal in the Software without restriction, including
   * without limitation the rights to use, copy, modify, merge, publish,
   * distribute, sublicense, and/or sell copies of the Software, and to
   * permit persons to whom the Software is furnished to do so, subject
   * to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included
   * in all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
   * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
   * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
   * IN NO EVENT SHALL THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE
   * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
   * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
   * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
   */
  
  #include "ReferenceCMAPTorsionIxn.h"
  #include "ReferenceForce.h"
  
  using std::vector;
  using namespace OpenMM;
  
  /**---------------------------------------------------------------------------------------
  
     Constructor
  
     --------------------------------------------------------------------------------------- */
  
  ReferenceCMAPTorsionIxn::ReferenceCMAPTorsionIxn(const vector<vector<vector<double> > >& coeff,
          const vector<int>& torsionMaps, const vector<vector<int> >& torsionIndices) :
          coeff(coeff), torsionMaps(torsionMaps), torsionIndices(torsionIndices), usePeriodic(false) {
  }
  
  void ReferenceCMAPTorsionIxn::setPeriodic(OpenMM::Vec3* vectors) {
      usePeriodic = true;
      boxVectors[0] = vectors[0];
      boxVectors[1] = vectors[1];
      boxVectors[2] = vectors[2];
  }
  
  /**---------------------------------------------------------------------------------------
  
     Calculate torsion interaction
  
     @param atomIndices      bond indices
     @param atomCoordinates  atom coordinates
     @param parameters       parameters
     @param forces           force array (forces added)
     @param energyByBond     bond energy
     @param energy           atom energy
  
     --------------------------------------------------------------------------------------- */
  
  void ReferenceCMAPTorsionIxn::calculateIxn(vector<Vec3>& atomCoordinates, vector<Vec3>& forces, double* totalEnergy) const {
      for (unsigned int i = 0; i < torsionMaps.size(); i++)
          calculateOneIxn(i, atomCoordinates, forces, totalEnergy);
  }
  
  /**---------------------------------------------------------------------------------------
  
     Calculate the interaction due to a single torsion pair
  
     @param index            the index of the torsion
     @param atomCoordinates  atom coordinates
     @param forces           force array (forces added)
     @param totalEnergy      total energy
  
       --------------------------------------------------------------------------------------- */
  
  void ReferenceCMAPTorsionIxn::calculateOneIxn(int index, vector<Vec3>& atomCoordinates, vector<Vec3>& forces,
                       double* totalEnergy) const {
      int map = torsionMaps[index];
      int a1 = torsionIndices[index][0];
      int a2 = torsionIndices[index][1];
      int a3 = torsionIndices[index][2];
      int a4 = torsionIndices[index][3];
      int b1 = torsionIndices[index][4];
      int b2 = torsionIndices[index][5];
      int b3 = torsionIndices[index][6];
      int b4 = torsionIndices[index][7];
  
      // Compute deltas between the various atoms involved.
  
      double deltaA[3][ReferenceForce::LastDeltaRIndex];
      double deltaB[3][ReferenceForce::LastDeltaRIndex];
      if (usePeriodic) {
          ReferenceForce::getDeltaRPeriodic(atomCoordinates[a2], atomCoordinates[a1], boxVectors, deltaA[0]);
          ReferenceForce::getDeltaRPeriodic(atomCoordinates[a2], atomCoordinates[a3], boxVectors, deltaA[1]);
          ReferenceForce::getDeltaRPeriodic(atomCoordinates[a4], atomCoordinates[a3], boxVectors, deltaA[2]);
          ReferenceForce::getDeltaRPeriodic(atomCoordinates[b2], atomCoordinates[b1], boxVectors, deltaB[0]);
          ReferenceForce::getDeltaRPeriodic(atomCoordinates[b2], atomCoordinates[b3], boxVectors, deltaB[1]);
          ReferenceForce::getDeltaRPeriodic(atomCoordinates[b4], atomCoordinates[b3], boxVectors, deltaB[2]);
      }
      else {
          ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a1], deltaA[0]);
          ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a3], deltaA[1]);
          ReferenceForce::getDeltaR(atomCoordinates[a4], atomCoordinates[a3], deltaA[2]);
          ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b1], deltaB[0]);
          ReferenceForce::getDeltaR(atomCoordinates[b2], atomCoordinates[b3], deltaB[1]);
          ReferenceForce::getDeltaR(atomCoordinates[b4], atomCoordinates[b3], deltaB[2]);
      }
  
      // Visual Studio complains if crossProduct declared as 'crossProduct[2][3]'
  
      double crossProductMemory[12];
      double* cpA[2];
      cpA[0] = crossProductMemory;
      cpA[1] = crossProductMemory + 3;
      double* cpB[2];
      cpB[0] = crossProductMemory + 6;
      cpB[1] = crossProductMemory + 9;
  
     // Compute the dihedral angles.
  
      double dotDihedral;
      double signOfAngle;
      double angleA =  getDihedralAngleBetweenThreeVectors(deltaA[0], deltaA[1], deltaA[2],
           cpA, &dotDihedral, deltaA[0], &signOfAngle, 1);
      double angleB =  getDihedralAngleBetweenThreeVectors(deltaB[0], deltaB[1], deltaB[2],
           cpB, &dotDihedral, deltaB[0], &signOfAngle, 1);
      angleA = fmod(angleA+2.0*M_PI, 2.0*M_PI);
      angleB = fmod(angleB+2.0*M_PI, 2.0*M_PI);
  
      // Identify which patch this is in.
  
      int size = (int) sqrt(coeff[map].size());
      double delta = 2*M_PI/size;
      int s = (int) (angleA/delta);
      int t = (int) (angleB/delta);
      const vector<double>& c = coeff[map][s+size*t];
      double da = angleA/delta-s;
      double db = angleB/delta-t;
  
      // Evaluate the spline to determine the energy and gradients.
  
      double energy = 0;
      double dEdA = 0;
      double dEdB = 0;
      for (int i = 3; i >= 0; i--) {
          energy = da*energy + ((c[i*4+3]*db + c[i*4+2])*db + c[i*4+1])*db + c[i*4+0];
          dEdA = db*dEdA + (3.0*c[i+3*4]*da + 2.0*c[i+2*4])*da + c[i+1*4];
          dEdB = da*dEdB + (3.0*c[i*4+3]*db + 2.0*c[i*4+2])*db + c[i*4+1];
      }
      dEdA /= delta;
      dEdB /= delta;
      if (totalEnergy != NULL)
          *totalEnergy += energy;
  
      // Apply the force to the first torsion.
  
      double forceFactors[4];
      double normCross1 = DOT3(cpA[0], cpA[0]);
      double normBC = deltaA[1][ReferenceForce::RIndex];
      forceFactors[0] = (-dEdA*normBC)/normCross1;
      double normCross2 = DOT3(cpA[1], cpA[1]);
      forceFactors[3] = (dEdA*normBC)/normCross2;
      forceFactors[1] = DOT3(deltaA[0], deltaA[1]);
      forceFactors[1] /= deltaA[1][ReferenceForce::R2Index];
      forceFactors[2] = DOT3(deltaA[2], deltaA[1]);
      forceFactors[2] /= deltaA[1][ReferenceForce::R2Index];
      for (int i = 0; i < 3; i++) {
          double f0 = forceFactors[0]*cpA[0][i];
          double f3 = forceFactors[3]*cpA[1][i];
          double s = forceFactors[1]*f0 - forceFactors[2]*f3;
          forces[a1][i] += f0;
          forces[a2][i] -= f0-s;
          forces[a3][i] -= f3+s;
          forces[a4][i] += f3;
      }
  
      // Apply the force to the second torsion.
  
      normCross1 = DOT3(cpB[0], cpB[0]);
      normBC = deltaB[1][ReferenceForce::RIndex];
      forceFactors[0] = (-dEdB*normBC)/normCross1;
      normCross2 = DOT3(cpB[1], cpB[1]);
      forceFactors[3] = (dEdB*normBC)/normCross2;
      forceFactors[1] = DOT3(deltaB[0], deltaB[1]);
      forceFactors[1] /= deltaB[1][ReferenceForce::R2Index];
      forceFactors[2] = DOT3(deltaB[2], deltaB[1]);
      forceFactors[2] /= deltaB[1][ReferenceForce::R2Index];
      for (int i = 0; i < 3; i++) {
          double f0 = forceFactors[0]*cpB[0][i];
          double f3 = forceFactors[3]*cpB[1][i];
          double s = forceFactors[1]*f0 - forceFactors[2]*f3;
          forces[b1][i] += f0;
          forces[b2][i] -= f0-s;
          forces[b3][i] -= f3+s;
          forces[b4][i] += f3;
      }
  }
  
  // ---------------------------------------------------------------------------------------
  
  /**---------------------------------------------------------------------------------------
  
     This is present only because we must define it to subclass ReferenceBondIxn.  It is never called.
  
     --------------------------------------------------------------------------------------- */
  
  void ReferenceCMAPTorsionIxn::calculateBondIxn(vector<int>& atomIndices, vector<Vec3>& atomCoordinates,
          vector<double>& parameters, vector<Vec3>& forces, double* totalEnergy, double* energyParamDerivs) {
  }
  