LjPmeKSpace Class Reference

#include <LjPmeKSpace.h>

Public Member Functions
	LjPmeKSpace (LjPmeGrid grid)
	~LjPmeKSpace ()
double	compute_energy (float *q_arr, const Lattice &lattice, double ewald, double *virial)

Detailed Description

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 13 of file LjPmeKSpace.h.

Constructor & Destructor Documentation

LjPmeKSpace()

LjPmeKSpace::LjPmeKSpace

(

LjPmeGrid

grid

)

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 12 of file LjPmeKSpace.C.

References LjPmeGrid::K1, LjPmeGrid::K2, LjPmeGrid::K3, NAMD_die(), LjPmeGrid::order, and order.

  : 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::~LjPmeKSpace

(

)

Definition at line 46 of file LjPmeKSpace.C.

                          {
  delete [] bm1;
  delete [] bm2;
  delete [] bm3;
}

Member Function Documentation

compute_energy()

double LjPmeKSpace::compute_energy	(	float *	q_arr,
		const Lattice &	lattice,
		double	ewald,
		double *	virial
	)

Definition at line 52 of file LjPmeKSpace.C.

References Lattice::a_r(), Lattice::b_r(), Lattice::c_r(), LjPmeGrid::K1, LjPmeGrid::K2, LjPmeGrid::K3, M_PI, Lattice::orthogonal(), SQRT_PI, Lattice::volume(), Vector::x, Vector::y, and Vector::z.

Referenced by LjPmeMgr::gridCalculation().

                                                  {
  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;
}

---

The documentation for this class was generated from the following files:

• LjPmeKSpace.h
• LjPmeKSpace.C