NAMD
Public Member Functions | List of all members
PmeKSpace Class Reference

#include <PmeKSpace.h>

Public Member Functions

 PmeKSpace (PmeGrid grid, int K2_start, int K2_end, int K3_start, int K3_end)
 
 ~PmeKSpace ()
 
double compute_energy (float q_arr[], const Lattice &lattice, double ewald, double virial[], int useCkLoop)
 
double compute_energy_orthogonal_helper (float q_arr[], const Lattice &lattice, double ewald, double virial[])
 
void compute_energy_orthogonal_subset (float q_arr[], double *recips, double partialVirial[], double *partialEnergy, int k1from, int k1to)
 

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 16 of file PmeKSpace.h.

Constructor & Destructor Documentation

◆ PmeKSpace()

PmeKSpace::PmeKSpace ( PmeGrid  grid,
int  K2_start,
int  K2_end,
int  K3_start,
int  K3_end 
)

Definition at line 63 of file PmeKSpace.C.

References compute_b_moduli(), PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, PmeGrid::order, and order.

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 }
compute_b_moduli
void compute_b_moduli(double *bm, int K, int order)
Definition: PmeKSpace.C:42
PmeGrid::K2
int K2
Definition: PmeBase.h:21
PmeGrid::K1
int K1
Definition: PmeBase.h:21
order
#define order
Definition: PmeRealSpace.C:235
PmeGrid::order
int order
Definition: PmeBase.h:23
PmeGrid::K3
int K3
Definition: PmeBase.h:21

◆ ~PmeKSpace()

PmeKSpace::~PmeKSpace ( )

Definition at line 125 of file PmeKSpace.C.

125  {
126  delete [] bm1;
127  delete [] bm2;
128  delete [] bm3;
129 
130  delete [] exp1;
131  delete [] exp2;
132  delete [] exp3;
133 }

Member Function Documentation

◆ compute_energy()

double PmeKSpace::compute_energy ( float  q_arr[],
const Lattice &  lattice,
double  ewald,
double  virial[],
int  useCkLoop 
)

Definition at line 321 of file PmeKSpace.C.

References Lattice::a(), Lattice::a_r(), Lattice::b(), Lattice::b_r(), Lattice::c(), Lattice::c_r(), compute_energy_orthogonal_helper(), PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, Vector::length(), M_PI, Lattice::orthogonal(), Vector::unit(), Lattice::volume(), Vector::x, Vector::y, and Vector::z.

Referenced by ComputePmeMgr::gridCalc2R(), and PmeXPencil::pme_kspace().

321  {
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 }
init_exp
static void init_exp(float *xp, int K, float recip, float kappa)
Definition: ComputeEwald.C:290
M_PI
#define M_PI
Definition: GoMolecule.C:40
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
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
Vector::length
NAMD_HOST_DEVICE BigReal length(void) const
Definition: Vector.h:202
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
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
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
Vector::y
BigReal y
Definition: Vector.h:74
Lattice::a
NAMD_HOST_DEVICE Vector a() const
Definition: Lattice.h:268
Vector::unit
NAMD_HOST_DEVICE Vector unit(void) const
Definition: Vector.h:215

◆ compute_energy_orthogonal_helper()

double PmeKSpace::compute_energy_orthogonal_helper ( float  q_arr[],
const Lattice &  lattice,
double  ewald,
double  virial[] 
)

Definition at line 240 of file PmeKSpace.C.

References Lattice::a_r(), ALLOCA, Lattice::b_r(), Lattice::c_r(), compute_energy_orthogonal_ckloop(), PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, M_PI, Lattice::volume(), Vector::x, Vector::y, and Vector::z.

Referenced by compute_energy().

240  {
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 }
init_exp
static void init_exp(float *xp, int K, float recip, float kappa)
Definition: ComputeEwald.C:290
M_PI
#define M_PI
Definition: GoMolecule.C:40
PmeGrid::K2
int K2
Definition: PmeBase.h:21
PmeGrid::K1
int K1
Definition: PmeBase.h:21
Vector::z
BigReal z
Definition: Vector.h:74
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
Lattice::c_r
NAMD_HOST_DEVICE Vector c_r() const
Definition: Lattice.h:286
PmeGrid::K3
int K3
Definition: PmeBase.h:21
ALLOCA
#define ALLOCA(TYPE, NAME, SIZE)
Definition: PmeKSpace.C:20
Vector::y
BigReal y
Definition: Vector.h:74
compute_energy_orthogonal_ckloop
static void compute_energy_orthogonal_ckloop(int first, int last, void *result, int paraNum, void *param)
Definition: PmeKSpace.C:214

◆ compute_energy_orthogonal_subset()

void PmeKSpace::compute_energy_orthogonal_subset ( float  q_arr[],
double *  recips,
double  partialVirial[],
double *  partialEnergy,
int  k1from,
int  k1to 
)

Definition at line 135 of file PmeKSpace.C.

References PmeGrid::K1, PmeGrid::K2, and PmeGrid::K3.

Referenced by compute_energy_orthogonal_ckloop().

135  {
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 }
PmeGrid::K2
int K2
Definition: PmeBase.h:21
PmeGrid::K1
int K1
Definition: PmeBase.h:21
PmeGrid::K3
int K3
Definition: PmeBase.h:21
The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.14