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 
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 
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 void PmeKSpace::init_exp(double *xp, int K, int k_start, int k_end, double recip) {
546  int i;
547  double fac;
548  fac = -piob*recip*recip;
549  int i_start = k_start;
550  int i_end = k_end;
551  if ( k_start > K/2 ) {
552  i_start = K - k_end + 1;
553  i_end = K - k_start + 1;
554  } else if ( k_end > K/2 ) {
555  i_end = K/2 + 1;
556  i_start = K - k_end + 1;
557  if ( k_start < i_start ) i_start = k_start;
558  }
559 
560  for (i=i_start; i < i_end; i++)
561  xp[i] = exp(fac*i*i);
562 }
static Node * Object()
Definition: Node.h:86
void compute_energy_orthogonal_subset(float q_arr[], double *recips, double partialVirial[], double *partialEnergy, int k1from, int k1to)
Definition: PmeKSpace.C:135
#define M_PI
Definition: GoMolecule.C:39
void compute_b_moduli(double *bm, int K, int order)
Definition: PmeKSpace.C:42
Vector a_r() const
Definition: Lattice.h:268
Definition: Vector.h:64
int K2
Definition: PmeBase.h:18
SimParameters * simParameters
Definition: Node.h:178
int K1
Definition: PmeBase.h:18
Vector c_r() const
Definition: Lattice.h:270
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
int orthogonal() const
Definition: Lattice.h:257
BigReal length(void) const
Definition: Vector.h:169
Vector b_r() const
Definition: Lattice.h:269
static void dftmod(double *bsp_mod, double *bsp_arr, int nfft)
Definition: PmeKSpace.C:22
#define order
Definition: PmeRealSpace.C:235
double compute_energy(float q_arr[], const Lattice &lattice, double ewald, double virial[], int useCkLoop)
Definition: PmeKSpace.C:321
int order
Definition: PmeBase.h:20
BigReal x
Definition: Vector.h:66
~PmeKSpace()
Definition: PmeKSpace.C:125
BigReal volume(void) const
Definition: Lattice.h:277
#define simParams
Definition: Output.C:127
int K3
Definition: PmeBase.h:18
double compute_energy_orthogonal_helper(float q_arr[], const Lattice &lattice, double ewald, double virial[])
Definition: PmeKSpace.C:240
#define ALLOCA(TYPE, NAME, SIZE)
Definition: PmeKSpace.C:20
BigReal y
Definition: Vector.h:66
Vector b() const
Definition: Lattice.h:253
static void compute_b_spline(REAL *__restrict frac, REAL *M, REAL *dM, int order)
Definition: PmeBase.inl:86
static void compute_energy_orthogonal_ckloop(int first, int last, void *result, int paraNum, void *param)
Definition: PmeKSpace.C:214
Vector a() const
Definition: Lattice.h:252
PmeKSpace(PmeGrid grid, int K2_start, int K2_end, int K3_start, int K3_end)
Definition: PmeKSpace.C:63
Vector unit(void) const
Definition: Vector.h:182
Vector c() const
Definition: Lattice.h:254