LjPmeKSpace.C

Go to the documentation of this file.

#include <stdlib.h>
#include "LjPmeKSpace.h"
#include "Lattice.h"
#include "LjPmeBase.h"

LjPmeKSpace::LjPmeKSpace(LjPmeGrid grid) 
  : myGrid(grid) {
  int K1, K2, K3, order;
  K1=myGrid.K1; K2=myGrid.K2, K3=myGrid.K3; order=myGrid.order;

  bm1 = new double[K1];
  bm2 = new double[K2];
  bm3 = new double[K3];

  switch (myGrid.order) {
    case 4:
      compute_b_moduli<4>(bm1, K1);
      compute_b_moduli<4>(bm2, K2);
      compute_b_moduli<4>(bm3, K3);
      break;
    case 6:
      compute_b_moduli<6>(bm1, K1);
      compute_b_moduli<6>(bm2, K2);
      compute_b_moduli<6>(bm3, K3);
      break;
    case 8:
      compute_b_moduli<8>(bm1, K1);
      compute_b_moduli<8>(bm2, K2);
      compute_b_moduli<8>(bm3, K3);
      break;
    case 10:
      compute_b_moduli<10>(bm1, K1);
      compute_b_moduli<10>(bm2, K2);
      compute_b_moduli<10>(bm3, K3);
      break;
    default: NAMD_die("unsupported LJPMEInterpOrder");
  }
}

LjPmeKSpace::~LjPmeKSpace() {
  delete [] bm1;
  delete [] bm2;
  delete [] bm3;
}

double LjPmeKSpace::compute_energy(float *q_arr, const Lattice &lattice,
                    double ewald, double *virial) {
  double energy = 0.0;
  double v0 = 0.0;
  double v1 = 0.0;
  double v2 = 0.0;
  double v3 = 0.0;
  double v4 = 0.0;
  double v5 = 0.0;

  int n;
  int k1, k2, k3, ind;
  int K1, K2, K3;
  int maxK1, maxK2, maxK3;

  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3;
  maxK1 = K1/2; maxK2 = K2/2; maxK3 = K3/2;
  // Constant needed to calculate the LJ-PME screening function
  double box_factor = -M_PI*SQRT_PI/(3.0*lattice.volume());
  double beta3 = ewald*ewald*ewald;
  double fac1 = -2.0*ewald*M_PI*M_PI;
  double fac2 = -M_PI*M_PI/(ewald*ewald);
  double fac3 = 2.0*M_PI*M_PI*M_PI*SQRT_PI;
  double fac4 = M_PI/ewald;

        //printf("LJ-PME: Grid(%d, %d, %d) \n", K1, K2, K3);
  //printf("box_factor: %6.12f, beta3: %6.12f, fact1: %6.12f, fact2: %6.12f, fact3: %6.12f \n",
  //box_factor, beta3, fac1, fac2, fac3);

  if ( lattice.orthogonal() ) {
    double recipx = lattice.a_r().x;
    double recipy = lattice.b_r().y;
    double recipz = lattice.c_r().z;
    //printf("recip_x: %6.12f, recip_y: %6.12f, recip_z: %6.12f \n", recipx, recipy, recipz);

    ind = 0;
    for ( k1=0; k1<K1; ++k1 ) {
      double m1, m11, b1;
      b1 = bm1[k1];
      // Grid indices in the upper half correspond to negative frequencies!
      int k1_s = k1<=maxK1 ? k1 : k1-K1;
      m1 = k1_s*recipx;
      m11 = m1*m1;
      for ( k2=0; k2<K2; ++k2 ) {
        double m2, m22, b1b2;
        b1b2 = b1*bm2[k2];
        // Grid indices in the upper half correspond to negative frequencies!
        int k2_s = k2<=maxK2 ? k2 : k2-K2;
        m2 = k2_s*recipy;
        m22 = m2*m2;
        // We use symmetric in z dimension. Unlike electrostatic, we can have m == 0
        for ( k3=0; k3<=maxK3; ++k3 ) {
          double m3, m33, msq, mFreq, vir, tempE;
          double b1b2b3, eTerm, q2, qr, qc;
          // Inverse b-spline module (b1*b2*b3)
          b1b2b3 = bm3[k3] *b1b2;
          m3 = k3*recipz;
          m33 = m3*m3;
          // Calculate the square root of frequency
          msq = m11 + m22 + m33;
          mFreq = sqrt(msq);
          qr = q_arr[ind]; qc=q_arr[ind+1];
          // Using symmetary, multiply the structure factor by 2
          q2 = 2.0*(qr*qr + qc*qc);
          // except for (k3 == 0) || ( k3 == K3/2)
          if ( (k3 == 0) || ( k3 == maxK3 && ! (K3 & 1) ) ) q2 *= 0.5;
          // Calculate the energy term (f(x) / b_spline) for this frequency m
          eTerm = b1b2b3*box_factor*((beta3 + fac1*msq)*exp(fac2*msq) + fac3*msq*mFreq*erfc(fac4*mFreq));
          // Calculate the Long-range PME contribution ro dispersion energy
          q_arr[ind] *= eTerm;
          q_arr[ind+1] *= eTerm;
          // Calculate the Energy
          tempE = q2*eTerm;
          energy += tempE;
          // Calculate the virial using GROMACS formulation
          vir = 3.0*box_factor*b1b2b3*(fac1*exp(fac2*msq)+fac3*mFreq*erfc(fac4*mFreq))*q2;
          v0 += tempE + vir*m11;
          v1 += vir*m1*m2;
          v2 += vir*m1*m3;
          v3 += tempE + vir*m22;
          v4 += vir*m2*m3;
          v5 += tempE + vir*m33;

          ind += 2;
        }
      }
    }
  } else {
    Vector recip1 = lattice.a_r();
    Vector recip2 = lattice.b_r();
    Vector recip3 = lattice.c_r();
    double recip2_x = recip2.x;
    double recip2_y = recip2.y;
    double recip2_z = recip2.z;
    double recip3_x = recip3.x;
    double recip3_y = recip3.y;
    double recip3_z = recip3.z;

    ind = 0;
    for ( k1=0; k1<K1; ++k1 ) {
      double b1; Vector m1;
      b1 = bm1[k1];
      // Grid indices in the upper half correspond to negative frequencies!
      int k1_s = k1<=maxK1 ? k1 : k1-K1;
      m1 = k1_s*recip1;
      for ( k2=0; k2<K2; ++k2 ) {
        double b1b2, m2_x, m2_y, m2_z;
        b1b2 = b1*bm2[k2];
        // Grid indices in the upper half correspond to negative frequencies!
        int k2_s = k2<=maxK2 ? k2 : k2-K2;
        m2_x = m1.x + k2_s*recip2_x;
        m2_y = m1.y + k2_s*recip2_y;
        m2_z = m1.z + k2_s*recip2_z;
        // We use symmetric in z dimension. Unlike electrostatic, we can have m == 0
        for ( k3=0; k3<=maxK3; ++k3 ) {
          double xp3, msq, mFreq, vir, tempE;
          double b1b2b3, eTerm, q2, qr, qc;
          double m_x, m_y, m_z;
          // Inverse b-spline module (b1*b2*b3)
          b1b2b3 = bm3[k3] *b1b2;
          m_x = m2_x + k3*recip3_x;
          m_y = m2_y + k3*recip3_y;
          m_z = m2_z + k3*recip3_z;
          // Calculate the square root of frequency
          msq = m_x*m_x + m_y*m_y + m_z*m_z;
          mFreq = sqrt(msq);
          qr = q_arr[ind]; qc=q_arr[ind+1];
          // Using symmetary, multiply the structure factor by 2
          q2 = 2.0*(qr*qr + qc*qc);
          // except for (k3 == 0) || ( k3 == K3/2)
          if ( (k3 == 0) || ( k3 == maxK3 && ! (K3 & 1) ) ) q2 *= 0.5;
          // Calculate the energy term (f(x) / b_spline) for this frequency m
          eTerm = b1b2b3*box_factor*((beta3 + fac1*msq)*exp(fac2*msq) + fac3*msq*mFreq*erfc(fac4*mFreq));
          // Calculate and write back convolution data to grid
          q_arr[ind] *= eTerm;
          q_arr[ind+1] *= eTerm;
          // Calculate the Long-range PME contribution ro dispersion energy
          tempE = q2*eTerm;
          energy += tempE;
          // Calculate the virial using GROMACS formulation
          vir = 3.0*box_factor*b1b2b3*(fac1*exp(fac2*msq)+fac3*mFreq*erfc(fac4*mFreq))*q2;
          v0 += tempE + vir*m_x*m_x;
          v1 += vir*m_x*m_y;
          v2 += vir*m_x*m_z;
          v3 += tempE + vir*m_y*m_y;
          v4 += vir*m_y*m_z;
          v5 += tempE + vir*m_z*m_z;
          ind += 2;
        }
      }
    }

  }

  virial[0] = 0.5 * v0;
  virial[1] = 0.5 * v1;
  virial[2] = 0.5 * v2;
  virial[3] = 0.5 * v3;
  virial[4] = 0.5 * v4;
  virial[5] = 0.5 * v5;
  return 0.5*energy;
}

template <int order>
void LjPmeKSpace::compute_b_moduli(double *bm, int K) {
  int i;
  double fr[3];

  double *M = new double[3*order];
  double *dM = new double[3*order];
  double *scratch = new double[K];

  fr[0]=fr[1]=fr[2]=0.0;
  compute_LjPme_b_spline<order>(fr,M,dM);  
  for (i=0; i<order; i++) bm[i] = M[i];
  for (i=order; i<K; i++) bm[i] = 0.0;
  this->dftmod(scratch, bm, K);
  // Store the inverse bspline module
  for (i=0; i<K; i++) bm[i] = 1.0/scratch[i];

  delete [] scratch;
  delete [] dM;
  delete [] M;
}

void LjPmeKSpace::dftmod(double *bsp_mod, double *bsp_arr, int nfft) {
  int j, k;
  double twopi, arg, sum1, sum2;
  double infft = 1.0/nfft;
/* Computes the modulus of the discrete fourier transform of bsp_arr, */
/*  storing it into bsp_mod */
  twopi =  2.0 * M_PI;

  for (k = 0; k <nfft; ++k) {
    sum1 = 0.;
    sum2 = 0.;
    for (j = 0; j < nfft; ++j) {
      arg = twopi * k * j * infft;
      sum1 += bsp_arr[j] * cos(arg);
      sum2 += bsp_arr[j] * sin(arg);
    }
    bsp_mod[k] = sum1*sum1 + sum2*sum2;
  }
}