LjPmeCompute.C

Go to the documentation of this file.

#include "LjPmeCompute.h"
#include "LjPmeMgr.h"
#include "Molecule.h"
#include "SimParameters.h"


void LjPmeCompute::initialize(const double *pmeCoord, const double *parameter,
                const double *parameter14, const int *type,
                double *forceNonbond, double *forceSlow,
                const int &ljTableDim, const int &nAtoms) { 
        if (! initialized) {
                ljPmeCoord = pmeCoord; 
                ljParameter = parameter;
                ljParameter14 = parameter14; 
                atomType = type;
                ljForceNonbond = forceNonbond; 
                ljForceSlow = forceSlow;
                tableDim = ljTableDim;
                numAtoms = nAtoms;
                SimParameters *simParm = Node::Object()->simParameters;
                ljPmeMgr.initialize(simParm, numAtoms);
                initialized = true;
        } else {
                NAMD_die("LjPmeCompute has already been initialized!");
        }
}

void LjPmeCompute::computeLJpotential(const double &alphaLJ,
        const double &cutoff, const Lattice &lattice,
        double virialNonbonded[][3], double virialSlow[][3],
        double &energyNonbond, double &energySlow,
        bool doEnergy, bool doSlow) {
        if(!initialized) {
                NAMD_die("LjPmeCompute has been called without initialization!");
        }
        SimParameters *simParm = Node::Object()->simParameters;
        // Initialize the energy and virial values to zero
        energySlow = energyNonbond = 0.0;
        for(int i = 0; i < 3; ++i) {
                for(int j = 0; j < 3; ++j) {
                        virialNonbonded[i][j] = 0.0;
                        virialSlow[i][j] = 0.0;
                }
        }

        // Calculate the short range and long range Nonbonded LJ potential
        if(doSlow) {
                if (doEnergy) {
                    if (simParm->LJPMESerial) {
                        // Calculate the LJ nonBonded term
                        this->computeNonbonded<true, true>(alphaLJ, cutoff,
                                lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
                    }
                } else {
                    energySlow = 0.0;
                    if (simParm->LJPMESerial) {
                        this->computeNonbonded<true, false>(alphaLJ, cutoff,
                                lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
                    }
                }
                // Compute reciprocal force, virial, energy, and self energy
                ljPmeMgr.computeLongRange(ljPmeCoord, lattice, alphaLJ,
                        ljForceSlow, energySlow, virialSlow, doEnergy);
        } else {
                if (doEnergy) {
                        this->computeNonbonded<false, true>(alphaLJ, cutoff,
                                lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
                } else {
                        this->computeNonbonded<false, false>(alphaLJ, cutoff,
                                lattice, virialNonbonded, virialSlow, energyNonbond, energySlow);
                }
        }
}

template <bool doSlow, bool doEnergy>
void LjPmeCompute::computeNonbonded(const double &alphaLJ,
        const double &cutoff, const Lattice &lattice,
        double virialNonbonded[][3], double virialSlow[][3],
        double &energyNonbond, double &energySlow) {
        Molecule *molecule = Node::Object()->molecule;
        double dispersion = 0.0; double repulsion = 0.0;
        double nonBondEnergy = 0.0;
        double slowEnergy = 0.0;
        double r2 = 0.0; 
        double rcut2 = cutoff * cutoff; 
        double rcut2inv = 1.0 / rcut2; 
        double rcut6inv = rcut2inv * rcut2inv * rcut2inv;
        double rcut12inv = rcut6inv * rcut6inv;
        double aRc = alphaLJ * cutoff;
        double aRc2 = aRc * aRc;
        int exclusion = 0; 
        // Screening function at cutoff distance (1 - g(alpha*cutoff))
        double screen_rc = (1.0 - (1 + aRc2 + 0.5*aRc2*aRc2)*exp(-aRc2));
        Vector iCoord, jCoord, distance;
        for(int i = 0; i < numAtoms; ++i) {
                for(int j = i+1; j < numAtoms; ++j) {
                        double energy1 = 0;  // energy for one interaction
                        iCoord = Vector(ljPmeCoord[4*i], ljPmeCoord[4*i+1], ljPmeCoord[4*i+2]);
                        jCoord = Vector(ljPmeCoord[4*j], ljPmeCoord[4*j+1], ljPmeCoord[4*j+2]);
                        // Shortest vector from jCoords to iCoords (iCoord - jCoords)
                        distance = lattice.delta(iCoord, jCoord);
                        // Distance square
                        r2 = distance.length2();

                        if (r2 < rcut2) {

                        // Check exclusion flag; 0: no exclusion. 1: full exclusion. 2: 1-4 scaling
                        #ifdef MEM_OPT_VERSION
                        // For now we dont support MEM_OPT
                        exclusion = molecule->checkExclByIdx(i, i, j); // Not sure if this is correct!
                        #else
                        exclusion = molecule->checkexcl(i, j);
                        #endif

                        // Check if the pair of atoms have full exclusion (1)
                        if(exclusion == 1) {
                                // Remove the reciprocal contribution for excluded pairs
                                if(doSlow ) {
                                        double r2inv = 1.0/r2;
                                        double r6inv = r2inv * r2inv * r2inv;
                                        double r12inv = r6inv * r6inv;
                                        double aR2 = r2 * alphaLJ * alphaLJ;
                                        double aR4 = aR2 * aR2;
                                        // Get the C6 term (4*eps*sig^6) using geometric combining rule
                                        double c6 = ljPmeCoord[4*i+3] * ljPmeCoord[4*j+3];
                                        // Derivative of Screening function (1 - g'(alpha*r))
                                        double screen_dr = (1.0 - (1.0 + aR2 + 0.5*aR4 + aR4*aR2/6.0)*exp(-aR2));
                                        // Subtract LJ dispersion force that explicitly included in reciprocal sum
                                        // Check! I think this should be positive! OpenMM uses negative sign! 
                                        double slowForce = 6.0 * c6 * screen_dr * r6inv * r2inv;

                                        if(doEnergy) {
                                                // Screening function (1 - g(alpha*r))
                                                double screen_r = (1.0 - (1.0 + aR2 + 0.5*aR4)*exp(-aR2));
                                                // Subtract LJ dispersion energy that explicitly included in reciprocal sum
                                                // Since dispersion is negative value, subtraction turns to addition
                                                energy1 += c6 * screen_r * r6inv;
                                                slowEnergy += energy1;
                                        }

                                        // Store the force
                                        Vector force_r = slowForce * distance;
                                        ljForceSlow[3*i]   += force_r.x;
                                        ljForceSlow[3*i+1] += force_r.y;
                                        ljForceSlow[3*i+2] += force_r.z;
                                        ljForceSlow[3*j]   -= force_r.x;
                                        ljForceSlow[3*j+1] -= force_r.y;
                                        ljForceSlow[3*j+2] -= force_r.z;
                                        // Store virial
                                        virialSlow[0][0] += force_r.x * distance.x;
                                        virialSlow[0][1] += force_r.x * distance.y;
                                        virialSlow[0][2] += force_r.x * distance.z;
                                        virialSlow[1][0] += force_r.y * distance.x;
                                        virialSlow[1][1] += force_r.y * distance.y;
                                        virialSlow[1][2] += force_r.y * distance.z;
                                        virialSlow[2][0] += force_r.z * distance.x;
                                        virialSlow[2][1] += force_r.z * distance.y;
                                        virialSlow[2][2] += force_r.z * distance.z;
                                }

                        } else /* if(r2 < rcut2) */ { // Make sure the distance is within the cutoff if pairs are not excluded (0, 2)
                                int iType = atomType[i];
                                int jType = atomType[j];
                                int typeIndex = 2*(iType * tableDim + jType);
                                double r2inv = 1.0/r2;
                                double r6inv = r2inv * r2inv * r2inv;
                                double r12inv = r6inv * r6inv;
                                if (exclusion == 2) { // Use 1-4 parameters 
                                        repulsion = ljParameter14[typeIndex];
                                        dispersion = ljParameter14[typeIndex + 1];
                                } else { // No full and 1-4 exclusion
                                        repulsion = ljParameter[typeIndex];
                                        dispersion = ljParameter[typeIndex + 1];
                                }
                                // LJ energy
                                if(doEnergy) {
                                        // Add standard LJ energy term
                                        energy1 += (repulsion*r12inv - dispersion*r6inv);
                                        nonBondEnergy += energy1;
                                }
                                // LJ force
                                double nonBondForce = (12.0*repulsion*r12inv - 6.0*dispersion*r6inv)*r2inv;
                                
                                // LJ dispersion approximation contribution to force and energy
                                if(doSlow) {
                                        double aR2 = r2 * alphaLJ * alphaLJ;
                                        double aR4 = aR2 * aR2;
                                        // Get the C6 term (4*eps*sig^6) using geometric combining rule
                                        double c6 = ljPmeCoord[4*i+3] * ljPmeCoord[4*j+3];
                                        // Derivative of Screening function (1 - g'(alpha*r))
                                        double screen_dr = (1.0 - (1.0 + aR2 + 0.5*aR4 + aR4*aR2/6.0)*exp(-aR2));
                                        // Add LJ dispersion approximation contribution to force
                                        nonBondForce += 6.0 * c6 * screen_dr * r6inv * r2inv;

                                        if(doEnergy) {
                                                // Screening function (1 - g(alpha*r))
                                                double screen_r = (1.0 - (1.0 + aR2 + 0.5*aR4)*exp(-aR2));
                                                // Add LJ dispersion approximation contribution to energy
                                                double energySlow1 = c6 * screen_r * r6inv;
                                                /* Calculate the shifted energy */
                                                // LJ energy at cutoff distance
                                                double potentialShift = (repulsion*rcut12inv - dispersion*rcut6inv);
                                                // Add the LJ dispersion approximation energy at cutoff
                                                potentialShift += c6 * screen_rc * rcut6inv;
                                                // Subtract the potentialShift to have zero energy at cutoff
                                                energySlow1 -= potentialShift;
                                                energy1 += energySlow1;
                                                nonBondEnergy += energySlow1;
                                        }
                                }

                                // Store the force
                                Vector force_r = nonBondForce * distance;
                                ljForceNonbond[3*i]   += force_r.x;
                                ljForceNonbond[3*i+1] += force_r.y;
                                ljForceNonbond[3*i+2] += force_r.z;
                                ljForceNonbond[3*j]   -= force_r.x;
                                ljForceNonbond[3*j+1] -= force_r.y;
                                ljForceNonbond[3*j+2] -= force_r.z;

                                // Store virial
                                virialNonbonded[0][0] += force_r.x * distance.x;
                                virialNonbonded[0][1] += force_r.x * distance.y;
                                virialNonbonded[0][2] += force_r.x * distance.z;
                                virialNonbonded[1][0] += force_r.y * distance.x;
                                virialNonbonded[1][1] += force_r.y * distance.y;
                                virialNonbonded[1][2] += force_r.y * distance.z;
                                virialNonbonded[2][0] += force_r.z * distance.x;
                                virialNonbonded[2][1] += force_r.z * distance.y;
                                virialNonbonded[2][2] += force_r.z * distance.z;
                        }

                        }
                }
        }
        energyNonbond += nonBondEnergy;
        energySlow += slowEnergy;
}