ReferenceCMAPTorsionIxn.cpp

/* Portions copyright (c) 2010 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;

/**---------------------------------------------------------------------------------------

   Constructor

   --------------------------------------------------------------------------------------- */

ReferenceCMAPTorsionIxn::ReferenceCMAPTorsionIxn(const vector<vector<vector<RealOpenMM> > >& coeff,
        const vector<int>& torsionMaps, const vector<vector<int> >& torsionIndices) :
        coeff(coeff), torsionMaps(torsionMaps), torsionIndices(torsionIndices) {
}

/**---------------------------------------------------------------------------------------

   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(RealOpenMM** atomCoordinates, RealOpenMM** forces, RealOpenMM* 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, RealOpenMM** atomCoordinates, RealOpenMM** forces,
                     RealOpenMM* 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.

    RealOpenMM deltaA[3][ReferenceForce::LastDeltaRIndex];
    ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a1], deltaA[0]);
    ReferenceForce::getDeltaR(atomCoordinates[a2], atomCoordinates[a3], deltaA[1]);
    ReferenceForce::getDeltaR(atomCoordinates[a4], atomCoordinates[a3], deltaA[2]);
    RealOpenMM deltaB[3][ReferenceForce::LastDeltaRIndex];
    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]'

    RealOpenMM crossProductMemory[12];
    RealOpenMM* cpA[2];
    cpA[0] = crossProductMemory;
    cpA[1] = crossProductMemory + 3;
    RealOpenMM* cpB[2];
    cpB[0] = crossProductMemory + 6;
    cpB[1] = crossProductMemory + 9;

   // Compute the dihedral angles.

    RealOpenMM dotDihedral;
    RealOpenMM signOfAngle;
    RealOpenMM angleA =  getDihedralAngleBetweenThreeVectors(deltaA[0], deltaA[1], deltaA[2],
         cpA, &dotDihedral, deltaA[0], &signOfAngle, 1);
    RealOpenMM 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(double(coeff[map].size()));
    RealOpenMM delta = 2*M_PI/size;
    int s = (int) (angleA/delta);
    int t = (int) (angleB/delta);
    const vector<RealOpenMM>& c = coeff[map][s+size*t];
    RealOpenMM da = angleA/delta-s;
    RealOpenMM db = angleB/delta-t;

    // Evaluate the spline to determine the energy and gradients.

    RealOpenMM energy = 0;
    RealOpenMM dEdA = 0;
    RealOpenMM 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.

    RealOpenMM forceFactors[4];
    RealOpenMM normCross1 = DOT3(cpA[0], cpA[0]);
    RealOpenMM normBC = deltaA[1][ReferenceForce::RIndex];
    forceFactors[0] = (-dEdA*normBC)/normCross1;
    RealOpenMM 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++) {
        RealOpenMM f0 = forceFactors[0]*cpA[0][i];
        RealOpenMM f3 = forceFactors[3]*cpA[1][i];
        RealOpenMM 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++) {
        RealOpenMM f0 = forceFactors[0]*cpB[0][i];
        RealOpenMM f3 = forceFactors[3]*cpB[1][i];
        RealOpenMM 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.

   --------------------------------------------------------------------------------------- */

int ReferenceCMAPTorsionIxn::calculateBondIxn(int* atomIndices, RealOpenMM** atomCoordinates,
     RealOpenMM* parameters, RealOpenMM** forces, RealOpenMM* energyByBond, RealOpenMM* energyByAtom) const {
         return 0; // but it must return a value to avoid MSVC compile error
}