NAMD
LjPmeRealSpace.C
Go to the documentation of this file.
1 
7 #include <string.h>
8 #include "LjPmeRealSpace.h"
9 #include "Lattice.h"
10 #include "LjPmeBase.h"
11 
13  myGrid(grid), N(natoms) {
14  int order = myGrid.order;
15  M = new double[3*N*order];
16  dM = new double[3*N*order];
17 }
18 
20  delete [] M;
21  delete [] dM;
22 }
23 
24 
25 template <int order>
26 void LjPmeRealSpace::fill_b_spline(double *p) {
27  double fr[3];
28  double *Mi, *dMi;
29  int i, stride, index;
30 
31  stride = 3*order;
32  Mi = M; dMi = dM;
33  for (i=0; i<N; i++) {
34  index = 4*i;
35  fr[0] = p[index]; fr[0] -= (int)(fr[0]);
36  fr[1] = p[index+1]; fr[1] -= (int)(fr[1]);
37  fr[2] = p[index+2]; fr[2] -= (int)(fr[2]);
38  compute_LjPme_b_spline<order>(fr, Mi, dMi);
39  Mi += stride;
40  dMi += stride;
41  }
42 }
43 
44 void LjPmeRealSpace::fill_charges(float **q_arr, char *f_arr, char *fz_arr,
45  double *p) {
46  switch (myGrid.order) {
47  case 4:
48  fill_charges_order<4>(q_arr, f_arr, fz_arr, p);
49  break;
50  case 6:
51  fill_charges_order<6>(q_arr, f_arr, fz_arr, p);
52  break;
53  case 8:
54  fill_charges_order<8>(q_arr, f_arr, fz_arr, p);
55  break;
56  case 10:
57  fill_charges_order<10>(q_arr, f_arr, fz_arr, p);
58  break;
59  default: NAMD_die("unsupported LJPMEInterpOrder");
60  }
61 }
62 
63 template <int order>
64 void LjPmeRealSpace::fill_charges_order(float **q_arr, char *f_arr, char *fz_arr,
65  double *p) {
66 
67  int i, j, k, l;
68  int K1, K2, K3, dim2, dim3;
69  double *Mi;
70  Mi = M;
71  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2; dim3=myGrid.dim3;
72  int stride = 3*order;
73 
74  this->fill_b_spline<order>(p);
75 
76  for (i=0; i<N; i++) {
77  double q;
78  int u1, u2, u2i, u3i;
79  int index = 4*i;
80  u1 = (int)(p[index]);
81  u2i = (int)(p[index+1]);
82  u3i = (int)(p[index+2]);
83  q = p[index+3];
84  u1 -= order;
85  u2i -= order;
86  u3i -= order;
87  u3i += 1;
88  for (j=0; j<order; j++) {
89  double m1;
90  int ind1;
91  m1 = Mi[j]*q;
92  u1++;
93  ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2;
94  u2 = u2i;
95  for (k=0; k<order; k++) {
96  double m1m2;
97  int ind2;
98  m1m2 = m1*Mi[order+k];
99  u2++;
100  ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));
101  float *qline = q_arr[ind2];
102  if ( ! qline ) {
103  q_arr[ind2] = qline = new float[dim3];
104  memset( (void*) qline, 0, dim3 * sizeof(float) );
105  }
106  f_arr[ind2] = 1;
107  for (l=0; l<order; l++) {
108  double m3;
109  int ind;
110  m3 = Mi[2*order + l];
111  int u3 = u3i + l;
112  ind = u3 + (u3 < 0 ? K3 : 0);
113  qline[ind] += m1m2*m3;
114  }
115  }
116  }
117  Mi += stride;
118  for (l=0; l<order; l++) {
119  int u3 = u3i + l;
120  int ind = u3 + (u3 < 0 ? K3 : 0);
121  fz_arr[ind] = 1;
122  }
123  }
124 }
125 
126 void LjPmeRealSpace::compute_scaledForces(const float * const *q_arr, double *pos,
127  double *force, const Lattice &lattice) {
128  switch (myGrid.order) {
129  case 4:
130  compute_scaledForces_order<4>(q_arr, pos, force, lattice);
131  break;
132  case 6:
133  compute_scaledForces_order<6>(q_arr, pos, force, lattice);
134  break;
135  case 8:
136  compute_scaledForces_order<8>(q_arr, pos, force, lattice);
137  break;
138  case 10:
139  compute_scaledForces_order<10>(q_arr, pos, force, lattice);
140  break;
141  default: NAMD_die("unsupported LJPMEInterpOrder");
142  }
143 }
144 
145 template <int order>
146 void LjPmeRealSpace::compute_scaledForces_order(const float * const *q_arr, double *pos,
147  double *force, const Lattice &lattice) {
148  int i, j, k, l, stride;
149  double f1, f2, f3;
150  double *Mi, *dMi;
151  int K1, K2, K3, dim2;
152 
153  Vector recip1 = lattice.a_r();
154  Vector recip2 = lattice.b_r();
155  Vector recip3 = lattice.c_r();
156  double r1x = recip1.x;
157  double r1y = recip1.y;
158  double r1z = recip1.z;
159  double r2x = recip2.x;
160  double r2y = recip2.y;
161  double r2z = recip2.z;
162  double r3x = recip3.x;
163  double r3y = recip3.y;
164  double r3z = recip3.z;
165 
166  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3; dim2=myGrid.dim2;
167  stride=3*order;
168  Mi = M; dMi = dM;
169 
170  for (i=0; i<N; i++) {
171  double q;
172  int u1, u2, u2i, u3i;
173  f1=f2=f3=0.0;
174  u1 = (int)(pos[4*i]);
175  u2i = (int)(pos[4*i+1]);
176  u3i = (int)(pos[4*i+2]);
177  q = pos[4*i+3];
178  u1 -= order;
179  u2i -= order;
180  u3i -= order;
181  u3i += 1;
182  for (j=0; j<order; j++) {
183  double m1, d1;
184  int ind1;
185  m1=Mi[j]*q;
186  d1=K1*dMi[j]*q;
187  u1++;
188  ind1 = (u1 + (u1 < 0 ? K1 : 0))*dim2;
189  u2 = u2i;
190  for (k=0; k<order; k++) {
191  double m2, d2, m1m2, m1d2, d1m2;
192  int ind2;
193  m2=Mi[order+k];
194  d2=K2*dMi[order+k];
195  m1m2=m1*m2;
196  m1d2=m1*d2;
197  d1m2=d1*m2;
198  u2++;
199  ind2 = ind1 + (u2 + (u2 < 0 ? K2 : 0));
200  const float *qline = q_arr[ind2];
201  for (l=0; l<order; l++) {
202  double term, m3, d3;
203  int ind;
204  m3=Mi[2*order+l];
205  d3=K3*dMi[2*order+l];
206  int u3 = u3i + l;
207  ind = u3 + (u3 < 0 ? K3 : 0);
208  term = qline[ind];
209  f1 -= d1m2 * m3 * term;
210  f2 -= m1d2 * m3 * term;
211  f3 -= m1m2 * d3 * term;
212  }
213  }
214  }
215  Mi += stride;
216  dMi += stride;
217 
218  force[3*i] += f1*r1x + f2*r2x + f3*r3x;
219  force[3*i+1] += f1*r1y + f2*r2y + f3*r3y;
220  force[3*i+2] += f1*r1z + f2*r2z + f3*r3z;
221  }
222 }
223 
void compute_scaledForces(const float *const *q_arr, double *pos, double *force, const Lattice &lattice)
Definition: Vector.h:72
BigReal z
Definition: Vector.h:74
int dim2
Definition: LjPmeBase.h:22
int dim3
Definition: LjPmeBase.h:22
#define order
Definition: PmeRealSpace.C:235
BigReal x
Definition: Vector.h:74
NAMD_HOST_DEVICE Vector a_r() const
Definition: Lattice.h:284
NAMD_HOST_DEVICE Vector b_r() const
Definition: Lattice.h:285
void NAMD_die(const char *err_msg)
Definition: common.C:147
int K3
Definition: LjPmeBase.h:21
NAMD_HOST_DEVICE Vector c_r() const
Definition: Lattice.h:286
LjPmeRealSpace(LjPmeGrid grid, int natoms)
void fill_charges(float **q_arr, char *f_arr, char *fz_arr, double *p)
BigReal y
Definition: Vector.h:74
int K1
Definition: LjPmeBase.h:21
int K2
Definition: LjPmeBase.h:21
int order
Definition: LjPmeBase.h:23