NAMD
PmeKSpace.C
Go to the documentation of this file.
1 
7 #include "PmeKSpace.h"
8 #include <math.h>
9 #include <stdlib.h>
10 
11 #ifndef _NO_ALLOCA_H
12 #include <alloca.h>
13 #endif
14 
15 #include "PmeBase.inl"
16 #include "SimParameters.h"
17 #include "Node.h"
18 #include "ComputeMoaMgr.decl.h"
19 
20 #define ALLOCA(TYPE,NAME,SIZE) TYPE *NAME = (TYPE *) alloca((SIZE)*sizeof(TYPE))
21 
22 static void dftmod(double *bsp_mod, double *bsp_arr, int nfft) {
23  int j, k;
24  double twopi, arg, sum1, sum2;
25  double infft = 1.0/nfft;
26 /* Computes the modulus of the discrete fourier transform of bsp_arr, */
27 /* storing it into bsp_mod */
28  twopi = 2.0 * M_PI;
29 
30  for (k = 0; k <nfft; ++k) {
31  sum1 = 0.;
32  sum2 = 0.;
33  for (j = 0; j < nfft; ++j) {
34  arg = twopi * k * j * infft;
35  sum1 += bsp_arr[j] * cos(arg);
36  sum2 += bsp_arr[j] * sin(arg);
37  }
38  bsp_mod[k] = sum1*sum1 + sum2*sum2;
39  }
40 }
41 
42 void compute_b_moduli(double *bm, int K, int order) {
43  int i;
44  double fr[3];
45 
46  double *M = new double[3*order];
47  double *dM = new double[3*order];
48  double *scratch = new double[K];
49 
50  fr[0]=fr[1]=fr[2]=0.0;
51  compute_b_spline(fr,M,dM,order);
52  for (i=0; i<order; i++) bm[i] = M[i];
53  for (i=order; i<K; i++) bm[i] = 0.0;
54  dftmod(scratch, bm, K);
55  for (i=0; i<K; i++) bm[i] = 1.0/scratch[i];
56 
57 
58  delete [] scratch;
59  delete [] dM;
60  delete [] M;
61 }
62 
63 PmeKSpace::PmeKSpace(PmeGrid grid, int K2_start, int K2_end, int K3_start, int K3_end)
64  : myGrid(grid), k2_start(K2_start), k2_end(K2_end),
65  k3_start(K3_start), k3_end(K3_end) {
66  int K1, K2, K3, order;
67  K1=myGrid.K1; K2=myGrid.K2, K3=myGrid.K3; order=myGrid.order;
68 
69  bm1 = new double[K1];
70  bm2 = new double[K2];
71  bm3 = new double[K3];
72 
73  exp1 = new double[K1/2 + 1];
74  exp2 = new double[K2/2 + 1];
75  exp3 = new double[K3/2 + 1];
76 
77  compute_b_moduli(bm1, K1, order);
78  compute_b_moduli(bm2, K2, order);
79  compute_b_moduli(bm3, K3, order);
80 
81 }
82 
83 #ifdef OPENATOM_VERSION
84 PmeKSpace::PmeKSpace(PmeGrid grid, int K2_start, int K2_end, int K3_start, int K3_end, CProxy_ComputeMoaMgr moaProxy)
85  : myGrid(grid), k2_start(K2_start), k2_end(K2_end),
86  k3_start(K3_start), k3_end(K3_end) {
87  int K1, K2, K3, order;
88  K1=myGrid.K1; K2=myGrid.K2, K3=myGrid.K3; order=myGrid.order;
89 
90  SimParameters *simParams = Node::Object()->simParameters;
91 
92  bm1 = new double[K1];
93  bm2 = new double[K2];
94  bm3 = new double[K3];
95 
96  exp1 = new double[K1/2 + 1];
97  exp2 = new double[K2/2 + 1];
98  exp3 = new double[K3/2 + 1];
99 
100  compute_b_moduli(bm1, K1, order);
101  compute_b_moduli(bm2, K2, order);
102  compute_b_moduli(bm3, K3, order);
103 
104  if ( simParams->openatomOn ) {
105  CkPrintf("######################################################\n");
106  CkPrintf("Entering recvBmX loop on processor %d \n", CkMyPe() );
107  int i;
108  for ( i=0; i<=K1; ++i) {
109  CkPrintf("bm1 value pre-send %d = %d \n", i, bm1[i] );
110  }
111  for ( i=0; i<=K1; ++i) {
112  CkPrintf("bm1 reference pre-send %d = %d \n", i, &bm1[i] );
113  }
114  CkPrintf("bm1: %d \n", *bm1);
115  moaProxy[CkMyPe()].recvB(K2_start, K2_end, K3_start, K3_end, K1, K2, K3, bm1, bm2, bm3, order);
116 // qmmmcp->recvBm1(bm1, K1, order);
117 // qmmmcp->recvBm2(bm2, K2, order);
118 // qmmmcp->recvBm3(bm3, K3, order);
119  }
120 
121 
122 }
123 #endif // OPENATOM_VERSION
124 
125 PmeKSpace::~PmeKSpace() {
126  delete [] bm1;
127  delete [] bm2;
128  delete [] bm3;
129 
130  delete [] exp1;
131  delete [] exp2;
132  delete [] exp3;
133 }
134 
135 void PmeKSpace::compute_energy_orthogonal_subset(float q_arr[], double *recips, double *partialVirial, double *partialEnergy, int k1from, int k1to){
136 
137  double energy = 0.0;
138  double v0 = 0.;
139  double v1 = 0.;
140  double v2 = 0.;
141  double v3 = 0.;
142  double v4 = 0.;
143  double v5 = 0.;
144 
145  int k1, k2, k3;
146  int K1, K2, K3;
147  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3;
148 
149  double recipx = recips[0];
150  double recipy = recips[1];
151  double recipz = recips[2];
152 
153  int ind = k1from*(k2_end-k2_start)*(k3_end-k3_start)*2;
154 
155  for ( k1=k1from; k1<=k1to; ++k1 ) {
156  double m1, m11, b1, xp1;
157  b1 = bm1[k1];
158  int k1_s = k1<=K1/2 ? k1 : k1-K1;
159  m1 = k1_s*recipx;
160  m11 = m1*m1;
161  xp1 = i_pi_volume*exp1[abs(k1_s)];
162  for ( k2=k2_start; k2<k2_end; ++k2 ) {
163  double m2, m22, b1b2, xp2;
164  b1b2 = b1*bm2[k2];
165  int k2_s = k2<=K2/2 ? k2 : k2-K2;
166  m2 = k2_s*recipy;
167  m22 = m2*m2;
168  xp2 = exp2[abs(k2_s)]*xp1;
169 
170  k3 = k3_start;
171  if (k1==0 && k2==0 && k3==0) {
172  q_arr[ind++] = 0.0;
173  q_arr[ind++] = 0.0;
174  ++k3;
175  }
176  for (; k3<k3_end; ++k3 ) {
177  double m3, m33, xp3, msq, imsq, vir, fac;
178  double theta3, theta, q2, qr, qc, C;
179  theta3 = bm3[k3] *b1b2;
180  m3 = k3*recipz;
181  m33 = m3*m3;
182  xp3 = exp3[k3];
183  qr = q_arr[ind]; qc=q_arr[ind+1];
184  q2 = 2*(qr*qr + qc*qc)*theta3;
185  if ( (k3 == 0) || ( k3 == K3/2 && ! (K3 & 1) ) ) q2 *= 0.5;
186  msq = m11 + m22 + m33;
187  imsq = 1.0/msq;
188  C = xp2*xp3*imsq;
189  theta = theta3*C;
190  q_arr[ind] *= theta;
191  q_arr[ind+1] *= theta;
192  vir = -2*(piob+imsq);
193  fac = q2*C;
194  energy += fac;
195  v0 += fac*(1.0+vir*m11);
196  v1 += fac*vir*m1*m2;
197  v2 += fac*vir*m1*m3;
198  v3 += fac*(1.0+vir*m22);
199  v4 += fac*vir*m2*m3;
200  v5 += fac*(1.0+vir*m33);
201  ind += 2;
202  }
203  }
204  }
205 
206  *partialEnergy = 0.5*energy;
207  partialVirial[0] = 0.5*v0;
208  partialVirial[1] = 0.5*v1;
209  partialVirial[2] = 0.5*v2;
210  partialVirial[3] = 0.5*v3;
211  partialVirial[4] = 0.5*v4;
212  partialVirial[5] = 0.5*v5;
213 }
214 static inline void compute_energy_orthogonal_ckloop(int first, int last, void *result, int paraNum, void *param){
215  for ( int i = first; i <= last; ++i ) {
216  void **params = (void **)param;
217  PmeKSpace *kspace = (PmeKSpace *)params[0];
218  float *q_arr = (float *)params[1];
219  double *recips = (double *)params[2];
220  double *partialEnergy = (double *)params[3];
221  double *partialVirial = (double *)params[4];
222  int *unitDist = (int *)params[5];
223 
224  int unit = unitDist[0];
225  int remains = unitDist[1];
226  int k1from, k1to;
227  if(i<remains){
228  k1from = i*(unit+1);
229  k1to = k1from+unit;
230  }else{
231  k1from = remains*(unit+1)+(i-remains)*unit;
232  k1to = k1from+unit-1;
233  }
234  double *pEnergy = partialEnergy+i;
235  double *pVirial = partialVirial+i*6;
236  kspace->compute_energy_orthogonal_subset(q_arr, recips, pVirial, pEnergy, k1from, k1to);
237  }
238 }
239 
240 double PmeKSpace::compute_energy_orthogonal_helper(float *q_arr, const Lattice &lattice, double ewald, double *virial) {
241  double energy = 0.0;
242  double v0 = 0.;
243  double v1 = 0.;
244  double v2 = 0.;
245  double v3 = 0.;
246  double v4 = 0.;
247  double v5 = 0.;
248 
249  int n;
250  int k1, k2, k3, ind;
251  int K1, K2, K3;
252 
253  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3;
254 
255  i_pi_volume = 1.0/(M_PI * lattice.volume());
256  piob = M_PI/ewald;
257  piob *= piob;
258 
259 
260  double recipx = lattice.a_r().x;
261  double recipy = lattice.b_r().y;
262  double recipz = lattice.c_r().z;
263 
264  init_exp(exp1, K1, 0, K1, recipx);
265  init_exp(exp2, K2, k2_start, k2_end, recipy);
266  init_exp(exp3, K3, k3_start, k3_end, recipz);
267 
268  double recips[] = {recipx, recipy, recipz};
269  int NPARTS=CmiMyNodeSize(); //this controls the granularity of loop parallelism
270  int maxParts = ( K1 * ( k2_end - k2_start ) * ( k3_end - k3_start ) + 127 ) / 128;
271  if ( NPARTS > maxParts ) NPARTS = maxParts;
272  if ( NPARTS > K1 ) NPARTS = K1;
273  ALLOCA(double, partialEnergy, NPARTS);
274  ALLOCA(double, partialVirial, 6*NPARTS);
275  int unitDist[] = {K1/NPARTS, K1%NPARTS};
276 
277  //parallelize the following loop using CkLoop
278  void *params[] = {this, q_arr, recips, partialEnergy, partialVirial, unitDist};
279 
280 #if CMK_SMP && USE_CKLOOP
281  CkLoop_Parallelize(compute_energy_orthogonal_ckloop, 6, (void *)params, NPARTS, 0, NPARTS-1);
282 #endif
283 /*
284  //The transformed loop used to compute energy
285  int unit = K1/NPARTS;
286  int remains = K1%NPARTS;
287  for(int i=0; i<NPARTS; i++){
288  int k1from, k1to;
289  if(i<remains){
290  k1from = i*(unit+1);
291  k1to = k1from+unit;
292  }else{
293  k1from = remains*(unit+1)+(i-remains)*unit;
294  k1to = k1from+unit-1;
295  }
296  double *pEnergy = partialEnergy+i;
297  double *pVirial = partialVirial+i*6;
298  compute_energy_orthogonal_subset(q_arr, recips, pVirial, pEnergy, k1from, k1to);
299  }
300 */
301 
302  for(int i=0; i<NPARTS; i++){
303  v0 += partialVirial[i*6+0];
304  v1 += partialVirial[i*6+1];
305  v2 += partialVirial[i*6+2];
306  v3 += partialVirial[i*6+3];
307  v4 += partialVirial[i*6+4];
308  v5 += partialVirial[i*6+5];
309  energy += partialEnergy[i];
310  }
311 
312  virial[0] = v0;
313  virial[1] = v1;
314  virial[2] = v2;
315  virial[3] = v3;
316  virial[4] = v4;
317  virial[5] = v5;
318  return energy;
319 }
320 
321 double PmeKSpace::compute_energy(float *q_arr, const Lattice &lattice, double ewald, double *virial, int useCkLoop) {
322  double energy = 0.0;
323  double v0 = 0.;
324  double v1 = 0.;
325  double v2 = 0.;
326  double v3 = 0.;
327  double v4 = 0.;
328  double v5 = 0.;
329 
330  int n;
331  int k1, k2, k3, ind;
332  int K1, K2, K3;
333 
334  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3;
335 
336  i_pi_volume = 1.0/(M_PI * lattice.volume());
337  piob = M_PI/ewald;
338  piob *= piob;
339 
340  if ( lattice.orthogonal() ) {
341  // if ( 0 ) { // JCP FOR TESTING
342  //This branch is the usual call path.
343 #if CMK_SMP && USE_CKLOOP
344  if ( useCkLoop ) {
345  return compute_energy_orthogonal_helper(q_arr, lattice, ewald, virial);
346  }
347 #endif
348  double recipx = lattice.a_r().x;
349  double recipy = lattice.b_r().y;
350  double recipz = lattice.c_r().z;
351 
352  init_exp(exp1, K1, 0, K1, recipx);
353  init_exp(exp2, K2, k2_start, k2_end, recipy);
354  init_exp(exp3, K3, k3_start, k3_end, recipz);
355 
356  ind = 0;
357  for ( k1=0; k1<K1; ++k1 ) {
358  double m1, m11, b1, xp1;
359  b1 = bm1[k1];
360  int k1_s = k1<=K1/2 ? k1 : k1-K1;
361  m1 = k1_s*recipx;
362  m11 = m1*m1;
363  xp1 = i_pi_volume*exp1[abs(k1_s)];
364  for ( k2=k2_start; k2<k2_end; ++k2 ) {
365  double m2, m22, b1b2, xp2;
366  b1b2 = b1*bm2[k2];
367  int k2_s = k2<=K2/2 ? k2 : k2-K2;
368  m2 = k2_s*recipy;
369  m22 = m2*m2;
370  xp2 = exp2[abs(k2_s)]*xp1;
371  k3 = k3_start;
372  if ( k1==0 && k2==0 && k3==0 ) {
373  q_arr[ind++] = 0.0;
374  q_arr[ind++] = 0.0;
375  ++k3;
376  }
377  for ( ; k3<k3_end; ++k3 ) {
378  double m3, m33, xp3, msq, imsq, vir, fac;
379  double theta3, theta, q2, qr, qc, C;
380  theta3 = bm3[k3] *b1b2;
381  m3 = k3*recipz;
382  m33 = m3*m3;
383  xp3 = exp3[k3];
384  qr = q_arr[ind]; qc=q_arr[ind+1];
385  q2 = 2*(qr*qr + qc*qc)*theta3;
386  if ( (k3 == 0) || ( k3 == K3/2 && ! (K3 & 1) ) ) q2 *= 0.5;
387  msq = m11 + m22 + m33;
388  imsq = 1.0/msq;
389  C = xp2*xp3*imsq;
390  theta = theta3*C;
391  q_arr[ind] *= theta;
392  q_arr[ind+1] *= theta;
393  vir = -2*(piob+imsq);
394  fac = q2*C;
395  energy += fac;
396  v0 += fac*(1.0+vir*m11);
397  v1 += fac*vir*m1*m2;
398  v2 += fac*vir*m1*m3;
399  v3 += fac*(1.0+vir*m22);
400  v4 += fac*vir*m2*m3;
401  v5 += fac*(1.0+vir*m33);
402  ind += 2;
403  }
404  }
405  }
406 
407  } else if ( cross(lattice.a(),lattice.b()).unit() == lattice.c().unit() ) {
408  // } else if ( 0 ) { // JCP FOR TESTING
409  Vector recip1 = lattice.a_r();
410  Vector recip2 = lattice.b_r();
411  Vector recip3 = lattice.c_r();
412  double recip3_x = recip3.x;
413  double recip3_y = recip3.y;
414  double recip3_z = recip3.z;
415  init_exp(exp3, K3, k3_start, k3_end, recip3.length());
416 
417  ind = 0;
418  for ( k1=0; k1<K1; ++k1 ) {
419  double b1; Vector m1;
420  b1 = bm1[k1];
421  int k1_s = k1<=K1/2 ? k1 : k1-K1;
422  m1 = k1_s*recip1;
423  // xp1 = i_pi_volume*exp1[abs(k1_s)];
424  for ( k2=k2_start; k2<k2_end; ++k2 ) {
425  double xp2, b1b2, m2_x, m2_y, m2_z;
426  b1b2 = b1*bm2[k2];
427  int k2_s = k2<=K2/2 ? k2 : k2-K2;
428  m2_x = m1.x + k2_s*recip2.x;
429  m2_y = m1.y + k2_s*recip2.y;
430  m2_z = m1.z + k2_s*recip2.z;
431  // xp2 = exp2[abs(k2_s)]*xp1;
432  xp2 = i_pi_volume*exp(-piob*(m2_x*m2_x+m2_y*m2_y+m2_z*m2_z));
433  k3 = k3_start;
434  if ( k1==0 && k2==0 && k3==0 ) {
435  q_arr[ind++] = 0.0;
436  q_arr[ind++] = 0.0;
437  ++k3;
438  }
439  for ( ; k3<k3_end; ++k3 ) {
440  double xp3, msq, imsq, vir, fac;
441  double theta3, theta, q2, qr, qc, C;
442  double m_x, m_y, m_z;
443  theta3 = bm3[k3] *b1b2;
444  m_x = m2_x + k3*recip3_x;
445  m_y = m2_y + k3*recip3_y;
446  m_z = m2_z + k3*recip3_z;
447  msq = m_x*m_x + m_y*m_y + m_z*m_z;
448  xp3 = exp3[k3];
449  qr = q_arr[ind]; qc=q_arr[ind+1];
450  q2 = 2*(qr*qr + qc*qc)*theta3;
451  if ( (k3 == 0) || ( k3 == K3/2 && ! (K3 & 1) ) ) q2 *= 0.5;
452  imsq = 1.0/msq;
453  C = xp2*xp3*imsq;
454  theta = theta3*C;
455  q_arr[ind] *= theta;
456  q_arr[ind+1] *= theta;
457  vir = -2*(piob+imsq);
458  fac = q2*C;
459  energy += fac;
460  v0 += fac*(1.0+vir*m_x*m_x);
461  v1 += fac*vir*m_x*m_y;
462  v2 += fac*vir*m_x*m_z;
463  v3 += fac*(1.0+vir*m_y*m_y);
464  v4 += fac*vir*m_y*m_z;
465  v5 += fac*(1.0+vir*m_z*m_z);
466  ind += 2;
467  }
468  }
469  }
470 
471  } else {
472  Vector recip1 = lattice.a_r();
473  Vector recip2 = lattice.b_r();
474  Vector recip3 = lattice.c_r();
475  double recip3_x = recip3.x;
476  double recip3_y = recip3.y;
477  double recip3_z = recip3.z;
478 
479  ind = 0;
480  for ( k1=0; k1<K1; ++k1 ) {
481  double b1; Vector m1;
482  b1 = bm1[k1];
483  int k1_s = k1<=K1/2 ? k1 : k1-K1;
484  m1 = k1_s*recip1;
485  // xp1 = i_pi_volume*exp1[abs(k1_s)];
486  for ( k2=k2_start; k2<k2_end; ++k2 ) {
487  double b1b2, m2_x, m2_y, m2_z;
488  b1b2 = b1*bm2[k2];
489  int k2_s = k2<=K2/2 ? k2 : k2-K2;
490  m2_x = m1.x + k2_s*recip2.x;
491  m2_y = m1.y + k2_s*recip2.y;
492  m2_z = m1.z + k2_s*recip2.z;
493  // xp2 = exp2[abs(k2_s)]*xp1;
494  k3 = k3_start;
495  if ( k1==0 && k2==0 && k3==0 ) {
496  q_arr[ind++] = 0.0;
497  q_arr[ind++] = 0.0;
498  ++k3;
499  }
500  for ( ; k3<k3_end; ++k3 ) {
501  double xp3, msq, imsq, vir, fac;
502  double theta3, theta, q2, qr, qc, C;
503  double m_x, m_y, m_z;
504  theta3 = bm3[k3] *b1b2;
505  m_x = m2_x + k3*recip3_x;
506  m_y = m2_y + k3*recip3_y;
507  m_z = m2_z + k3*recip3_z;
508  msq = m_x*m_x + m_y*m_y + m_z*m_z;
509  // xp3 = exp3[k3];
510  xp3 = i_pi_volume*exp(-piob*msq);
511  qr = q_arr[ind]; qc=q_arr[ind+1];
512  q2 = 2*(qr*qr + qc*qc)*theta3;
513  if ( (k3 == 0) || ( k3 == K3/2 && ! (K3 & 1) ) ) q2 *= 0.5;
514  imsq = 1.0/msq;
515  C = xp3*imsq;
516  theta = theta3*C;
517  q_arr[ind] *= theta;
518  q_arr[ind+1] *= theta;
519  vir = -2*(piob+imsq);
520  fac = q2*C;
521  energy += fac;
522  v0 += fac*(1.0+vir*m_x*m_x);
523  v1 += fac*vir*m_x*m_y;
524  v2 += fac*vir*m_x*m_z;
525  v3 += fac*(1.0+vir*m_y*m_y);
526  v4 += fac*vir*m_y*m_z;
527  v5 += fac*(1.0+vir*m_z*m_z);
528  ind += 2;
529  }
530  }
531  }
532 
533  }
534 
535  virial[0] = 0.5 * v0;
536  virial[1] = 0.5 * v1;
537  virial[2] = 0.5 * v2;
538  virial[3] = 0.5 * v3;
539  virial[4] = 0.5 * v4;
540  virial[5] = 0.5 * v5;
541  return 0.5*energy;
542 }
543 
544 
545 double PmeKSpace::compute_energy_LJPME(float *q_arr, const Lattice &lattice, double ewald, double *virial, int useCkLoop) {
546  // reimplement calculation from LjPmeKSpace::compute_energy()
547  // using the indexing scheme here for distributed grid
548  // XXX Still need to optimize: one erfc per grid point, employ CkLoop
549  double energy = 0.0;
550  double v0 = 0.0;
551  double v1 = 0.0;
552  double v2 = 0.0;
553  double v3 = 0.0;
554  double v4 = 0.0;
555  double v5 = 0.0;
556 
557  int n;
558  int k1, k2, k3, ind;
559  int K1, K2, K3;
560  int maxK1, maxK2, maxK3;
561 
562  K1=myGrid.K1; K2=myGrid.K2; K3=myGrid.K3;
563  maxK1 = K1/2; maxK2 = K2/2; maxK3 = K3/2;
564 
565  // Constant needed to calculate the LJ-PME screening function
566  double box_factor = -M_PI*SQRT_PI/(3.0*lattice.volume());
567  double beta3 = ewald*ewald*ewald;
568  double fac1 = -2.0*ewald*M_PI*M_PI;
569  double fac2 = -M_PI*M_PI/(ewald*ewald);
570  double fac3 = 2.0*M_PI*M_PI*M_PI*SQRT_PI;
571  double fac4 = M_PI/ewald;
572 
573  if ( lattice.orthogonal() ) {
574  double recipx = lattice.a_r().x;
575  double recipy = lattice.b_r().y;
576  double recipz = lattice.c_r().z;
577  //printf("recip_x: %6.12f, recip_y: %6.12f, recip_z: %6.12f \n", recipx, recipy, recipz);
578 
579  ind = 0;
580  for ( k1=0; k1<K1; ++k1 ) {
581  double m1, m11, b1;
582  b1 = bm1[k1];
583  // Grid indices in the upper half correspond to negative frequencies!
584  int k1_s = k1<=maxK1 ? k1 : k1-K1;
585  m1 = k1_s*recipx;
586  m11 = m1*m1;
587  for ( k2=k2_start; k2<k2_end; ++k2 ) {
588  double m2, m22, b1b2;
589  b1b2 = b1*bm2[k2];
590  // Grid indices in the upper half correspond to negative frequencies!
591  int k2_s = k2<=maxK2 ? k2 : k2-K2;
592  m2 = k2_s*recipy;
593  m22 = m2*m2;
594  // We use symmetric in z dimension. Unlike electrostatic, we can have m == 0
595  for ( k3=k3_start; k3<k3_end; ++k3 ) {
596  double m3, m33, msq, mFreq, vir, tempE;
597  double b1b2b3, eTerm, q2, qr, qc;
598  // Inverse b-spline module (b1*b2*b3)
599  b1b2b3 = bm3[k3] *b1b2;
600  m3 = k3*recipz;
601  m33 = m3*m3;
602  // Calculate the square root of frequency
603  msq = m11 + m22 + m33;
604  mFreq = sqrt(msq);
605  qr = q_arr[ind]; qc=q_arr[ind+1];
606  // Using symmetary, multiply the structure factor by 2
607  q2 = 2.0*(qr*qr + qc*qc);
608  // except for (k3 == 0) || ( k3 == K3/2)
609  if ( (k3 == 0) || ( k3 == maxK3 && ! (K3 & 1) ) ) q2 *= 0.5;
610  // Calculate the energy term (f(x) / b_spline) for this frequency m
611  eTerm = b1b2b3*box_factor*((beta3 + fac1*msq)*exp(fac2*msq) + fac3*msq*mFreq*erfc(fac4*mFreq));
612  // Calculate the Long-range PME contribution ro dispersion energy
613  q_arr[ind] *= eTerm;
614  q_arr[ind+1] *= eTerm;
615  // Calculate the Energy
616  tempE = q2*eTerm;
617  energy += tempE;
618  // Calculate the virial using GROMACS formulation
619  vir = 3.0*box_factor*b1b2b3*(fac1*exp(fac2*msq)+fac3*mFreq*erfc(fac4*mFreq))*q2;
620  v0 += tempE + vir*m11;
621  v1 += vir*m1*m2;
622  v2 += vir*m1*m3;
623  v3 += tempE + vir*m22;
624  v4 += vir*m2*m3;
625  v5 += tempE + vir*m33;
626  ind += 2;
627  }
628  }
629  }
630  } else {
631  Vector recip1 = lattice.a_r();
632  Vector recip2 = lattice.b_r();
633  Vector recip3 = lattice.c_r();
634  double recip2_x = recip2.x;
635  double recip2_y = recip2.y;
636  double recip2_z = recip2.z;
637  double recip3_x = recip3.x;
638  double recip3_y = recip3.y;
639  double recip3_z = recip3.z;
640 
641  ind = 0;
642  for ( k1=0; k1<K1; ++k1 ) {
643  double b1; Vector m1;
644  b1 = bm1[k1];
645  // Grid indices in the upper half correspond to negative frequencies!
646  int k1_s = k1<=maxK1 ? k1 : k1-K1;
647  m1 = k1_s*recip1;
648  for ( k2=k2_start; k2<k2_end; ++k2 ) {
649  double b1b2, m2_x, m2_y, m2_z;
650  b1b2 = b1*bm2[k2];
651  // Grid indices in the upper half correspond to negative frequencies!
652  int k2_s = k2<=maxK2 ? k2 : k2-K2;
653  m2_x = m1.x + k2_s*recip2_x;
654  m2_y = m1.y + k2_s*recip2_y;
655  m2_z = m1.z + k2_s*recip2_z;
656  // We use symmetric in z dimension. Unlike electrostatic, we can have m == 0
657  for ( k3=k3_start; k3<k3_end; ++k3 ) {
658  double msq, mFreq, vir, tempE;
659  double b1b2b3, eTerm, q2, qr, qc;
660  double m_x, m_y, m_z;
661  // Inverse b-spline module (b1*b2*b3)
662  b1b2b3 = bm3[k3] *b1b2;
663  m_x = m2_x + k3*recip3_x;
664  m_y = m2_y + k3*recip3_y;
665  m_z = m2_z + k3*recip3_z;
666  // Calculate the square root of frequency
667  msq = m_x*m_x + m_y*m_y + m_z*m_z;
668  mFreq = sqrt(msq);
669  qr = q_arr[ind]; qc=q_arr[ind+1];
670  // Using symmetary, multiply the structure factor by 2
671  q2 = 2.0*(qr*qr + qc*qc);
672  // except for (k3 == 0) || ( k3 == K3/2)
673  if ( (k3 == 0) || ( k3 == maxK3 && ! (K3 & 1) ) ) q2 *= 0.5;
674  // Calculate the energy term (f(x) / b_spline) for this frequency m
675  eTerm = b1b2b3*box_factor*((beta3 + fac1*msq)*exp(fac2*msq) + fac3*msq*mFreq*erfc(fac4*mFreq));
676  // Calculate and write back convolution data to grid
677  q_arr[ind] *= eTerm;
678  q_arr[ind+1] *= eTerm;
679  // Calculate the Long-range PME contribution ro dispersion energy
680  tempE = q2*eTerm;
681  energy += tempE;
682  // Calculate the virial using GROMACS formulation
683  vir = 3.0*box_factor*b1b2b3*(fac1*exp(fac2*msq)+fac3*mFreq*erfc(fac4*mFreq))*q2;
684  v0 += tempE + vir*m_x*m_x;
685  v1 += vir*m_x*m_y;
686  v2 += vir*m_x*m_z;
687  v3 += tempE + vir*m_y*m_y;
688  v4 += vir*m_y*m_z;
689  v5 += tempE + vir*m_z*m_z;
690  ind += 2;
691  }
692  }
693  }
694  }
695 
696  virial[0] = 0.5 * v0;
697  virial[1] = 0.5 * v1;
698  virial[2] = 0.5 * v2;
699  virial[3] = 0.5 * v3;
700  virial[4] = 0.5 * v4;
701  virial[5] = 0.5 * v5;
702  return 0.5*energy;
703 }
704 
705 
706 void PmeKSpace::init_exp(double *xp, int K, int k_start, int k_end, double recip) {
707  int i;
708  double fac;
709  fac = -piob*recip*recip;
710  int i_start = k_start;
711  int i_end = k_end;
712  if ( k_start > K/2 ) {
713  i_start = K - k_end + 1;
714  i_end = K - k_start + 1;
715  } else if ( k_end > K/2 ) {
716  i_end = K/2 + 1;
717  i_start = K - k_end + 1;
718  if ( k_start < i_start ) i_start = k_start;
719  }
720 
721  for (i=i_start; i < i_end; i++)
722  xp[i] = exp(fac*i*i);
723 }
Node::Object
static Node * Object()
Definition: Node.h:86
PmeKSpace::compute_energy_orthogonal_subset
void compute_energy_orthogonal_subset(float q_arr[], double *recips, double partialVirial[], double *partialEnergy, int k1from, int k1to)
Definition: PmeKSpace.C:135
M_PI
#define M_PI
Definition: GoMolecule.C:40
compute_b_moduli
void compute_b_moduli(double *bm, int K, int order)
Definition: PmeKSpace.C:42
Lattice::c
NAMD_HOST_DEVICE Vector c() const
Definition: Lattice.h:270
Vector
Definition: Vector.h:72
PmeGrid::K2
int K2
Definition: PmeBase.h:21
Node::simParameters
SimParameters * simParameters
Definition: Node.h:181
PmeGrid::K1
int K1
Definition: PmeBase.h:21
Vector::z
BigReal z
Definition: Vector.h:74
Lattice::orthogonal
NAMD_HOST_DEVICE int orthogonal() const
Definition: Lattice.h:273
SQRT_PI
#define SQRT_PI
Definition: ComputeExt.C:30
dftmod
static void dftmod(double *bsp_mod, double *bsp_arr, int nfft)
Definition: PmeKSpace.C:22
order
#define order
Definition: PmeRealSpace.C:235
PmeKSpace::compute_energy
double compute_energy(float q_arr[], const Lattice &lattice, double ewald, double virial[], int useCkLoop)
Definition: PmeKSpace.C:321
PmeGrid::order
int order
Definition: PmeBase.h:23
Vector::length
NAMD_HOST_DEVICE BigReal length(void) const
Definition: Vector.h:202
PmeKSpace::compute_energy_LJPME
double compute_energy_LJPME(float q_arr[], const Lattice &lattice, double LJewald, double virial[], int useCkLoop)
Definition: PmeKSpace.C:545
Vector::x
BigReal x
Definition: Vector.h:74
Lattice::volume
NAMD_HOST_DEVICE BigReal volume(void) const
Definition: Lattice.h:293
Lattice::a_r
NAMD_HOST_DEVICE Vector a_r() const
Definition: Lattice.h:284
Lattice::b_r
NAMD_HOST_DEVICE Vector b_r() const
Definition: Lattice.h:285
PmeKSpace::~PmeKSpace
~PmeKSpace()
Definition: PmeKSpace.C:125
Lattice::c_r
NAMD_HOST_DEVICE Vector c_r() const
Definition: Lattice.h:286
Lattice::b
NAMD_HOST_DEVICE Vector b() const
Definition: Lattice.h:269
simParams
#define simParams
Definition: Output.C:131
PmeGrid::K3
int K3
Definition: PmeBase.h:21
PmeKSpace::compute_energy_orthogonal_helper
double compute_energy_orthogonal_helper(float q_arr[], const Lattice &lattice, double ewald, double virial[])
Definition: PmeKSpace.C:240
ALLOCA
#define ALLOCA(TYPE, NAME, SIZE)
Definition: PmeKSpace.C:20
Vector::y
BigReal y
Definition: Vector.h:74
Lattice::a
NAMD_HOST_DEVICE Vector a() const
Definition: Lattice.h:268
compute_b_spline
static void compute_b_spline(REAL *__restrict frac, REAL *M, REAL *dM, int order)
Definition: PmeBase.inl:86
compute_energy_orthogonal_ckloop
static void compute_energy_orthogonal_ckloop(int first, int last, void *result, int paraNum, void *param)
Definition: PmeKSpace.C:214
Vector::unit
NAMD_HOST_DEVICE Vector unit(void) const
Definition: Vector.h:215
PmeKSpace::PmeKSpace
PmeKSpace(PmeGrid grid, int K2_start, int K2_end, int K3_start, int K3_end)
Definition: PmeKSpace.C:63
Generated by   doxygen 1.8.14