LjPmeRealSpace.C

Go to the documentation of this file.

#include <string.h>
#include "LjPmeRealSpace.h"
#include "Lattice.h"
#include "LjPmeBase.h"

LjPmeRealSpace::LjPmeRealSpace(LjPmeGrid grid, int natoms) : 
  myGrid(grid), N(natoms) {
  int order = myGrid.order;
  M = new double[3*N*order];
  dM = new double[3*N*order];
}

LjPmeRealSpace::~LjPmeRealSpace() {
  delete [] M;
  delete [] dM;
}


template <int order>
void LjPmeRealSpace::fill_b_spline(double *p) {
  double fr[3]; 
  double *Mi, *dMi;
  int i, stride, index;

  stride = 3*order;
  Mi = M; dMi = dM;
  for (i=0; i<N; i++) {
    index = 4*i;
    fr[0] = p[index];   fr[0] -= (int)(fr[0]);
    fr[1] = p[index+1]; fr[1] -= (int)(fr[1]);
    fr[2] = p[index+2]; fr[2] -= (int)(fr[2]);
    compute_LjPme_b_spline<order>(fr, Mi, dMi);
    Mi += stride;
    dMi += stride;
  }
}

void LjPmeRealSpace::fill_charges(float **q_arr, char *f_arr, char *fz_arr,
      double *p) {
  switch (myGrid.order) {
    case 4:
      fill_charges_order<4>(q_arr, f_arr, fz_arr, p);
      break;
    case 6:
      fill_charges_order<6>(q_arr, f_arr, fz_arr, p);
      break;
    case 8:
      fill_charges_order<8>(q_arr, f_arr, fz_arr, p);
      break;
    case 10:
      fill_charges_order<10>(q_arr, f_arr, fz_arr, p);
      break;
    default: NAMD_die("unsupported LJPMEInterpOrder");
  }
}

template <int order>
void LjPmeRealSpace::fill_charges_order(float **q_arr, char *f_arr, char *fz_arr,
      double *p) {
  
  int i, j, k, l;
  int K1, K2, K3, dim2, dim3;
  double *Mi;
  Mi = M;
  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2; dim3=myGrid.dim3;
  int stride = 3*order;

  this->fill_b_spline<order>(p);

  for (i=0; i<N; i++) {
    double q;
    int u1, u2, u2i, u3i;
    int index = 4*i;
    u1 = (int)(p[index]);
    u2i = (int)(p[index+1]);
    u3i = (int)(p[index+2]);
    q = p[index+3];
    u1 -= order;
    u2i -= order;
    u3i -= order;
    u3i += 1;
    for (j=0; j<order; j++) {
      double m1;
      int ind1;
      m1 = Mi[j]*q;
      u1++;
      ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2;
      u2 = u2i;
      for (k=0; k<order; k++) {
        double m1m2;
              int ind2;
        m1m2 = m1*Mi[order+k];
        u2++;
        ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));
        float *qline = q_arr[ind2];
        if ( ! qline ) {
          q_arr[ind2] = qline = new float[dim3];
          memset( (void*) qline, 0, dim3 * sizeof(float) );
        }
        f_arr[ind2] = 1;
        for (l=0; l<order; l++) {
          double m3;
          int ind;
          m3 = Mi[2*order + l];
          int u3 = u3i + l;
          ind = u3 + (u3 < 0 ? K3 : 0);
          qline[ind] += m1m2*m3; 
        }
      }
    }
    Mi += stride;
    for (l=0; l<order; l++) {
      int u3 = u3i + l;
      int ind = u3 + (u3 < 0 ? K3 : 0);
      fz_arr[ind] = 1;
    }
  }
}

void LjPmeRealSpace::compute_scaledForces(const float * const *q_arr, double *pos,
      double *force, const Lattice &lattice) {
  switch (myGrid.order) {
    case 4:
      compute_scaledForces_order<4>(q_arr, pos, force, lattice);
      break;
    case 6:
      compute_scaledForces_order<6>(q_arr, pos, force, lattice);
      break;
    case 8:
      compute_scaledForces_order<8>(q_arr, pos, force, lattice);
      break;
    case 10:
      compute_scaledForces_order<10>(q_arr, pos, force, lattice);
      break;
    default: NAMD_die("unsupported LJPMEInterpOrder");
  }
}

template <int order>
void LjPmeRealSpace::compute_scaledForces_order(const float * const *q_arr, double *pos,
      double *force, const Lattice &lattice) {
  int i, j, k, l, stride;
  double f1, f2, f3;
  double *Mi, *dMi;
  int K1, K2, K3, dim2;

  Vector recip1 = lattice.a_r();
  Vector recip2 = lattice.b_r();
  Vector recip3 = lattice.c_r();
  double r1x = recip1.x;
  double r1y = recip1.y;
  double r1z = recip1.z;
  double r2x = recip2.x;
  double r2y = recip2.y;
  double r2z = recip2.z;
  double r3x = recip3.x;
  double r3y = recip3.y;
  double r3z = recip3.z;

  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2;
  stride=3*order;
  Mi = M; dMi = dM;
 
  for (i=0; i<N; i++) {
    double q;
    int u1, u2, u2i, u3i;
    f1=f2=f3=0.0;
    u1 = (int)(pos[4*i]);
    u2i = (int)(pos[4*i+1]);
    u3i = (int)(pos[4*i+2]);
    q = pos[4*i+3];
    u1 -= order;
    u2i -= order;
    u3i -= order;
    u3i += 1;
    for (j=0; j<order; j++) {
      double m1, d1;
      int ind1;
      m1=Mi[j]*q;
      d1=K1*dMi[j]*q;
      u1++;
      ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2; 
      u2 = u2i;
      for (k=0; k<order; k++) {
        double m2, d2, m1m2, m1d2, d1m2;
        int ind2;
        m2=Mi[order+k];
        d2=K2*dMi[order+k];
        m1m2=m1*m2;
        m1d2=m1*d2;
        d1m2=d1*m2;
        u2++;
        ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));
        const float *qline = q_arr[ind2];
        for (l=0; l<order; l++) {
          double term, m3, d3;
          int ind;
          m3=Mi[2*order+l];
          d3=K3*dMi[2*order+l];
          int u3 = u3i + l;
          ind = u3 + (u3 < 0 ? K3 : 0);
          term = qline[ind];
          f1 -= d1m2 * m3 * term;
          f2 -= m1d2 * m3 * term;
          f3 -= m1m2 * d3 * term;
        }
      }
    }
    Mi += stride;
    dMi += stride;

    force[3*i]   += f1*r1x + f2*r2x + f3*r3x;
    force[3*i+1] += f1*r1y + f2*r2y + f3*r3y;
    force[3*i+2] += f1*r1z + f2*r2z + f3*r3z;
  }
}