NAMD
ComputeEwald.C
Go to the documentation of this file.
1 
7 /*
8  Forwards atoms to master node for force evaluation.
9 */
10 
11 #include "Node.h"
12 #include "PatchMap.h"
13 #include "PatchMap.inl"
14 #include "AtomMap.h"
15 #include "ComputeEwald.h"
16 #include "PatchMgr.h"
17 #include "Molecule.h"
18 #include "ReductionMgr.h"
19 #include "ComputeMgr.h"
20 #include "ComputeMgr.decl.h"
21 #include "ComputeNonbondedUtil.h"
22 #include "SimParameters.h"
23 #include "PmeBase.h"
24 #include <stdio.h>
25 
26 //#define DEBUGM
27 #define MIN_DEBUG_LEVEL 1
28 #include "Debug.h"
29 
31 public:
32  float x, y, z;
33  float cg;
34 };
35 
36 // table mapping atom types to array of structure factor combinations.
37 // For example, for ntypes=3, we have grids corresponding to atom types
38 // 0-0, 0-1, 0-2, 1-1, 1-2, 2-2, so type 0 gets added to grids 0,1,2, type 1
39 // to grids 1,3,4, and typ 2 to 2, 4, 5.
40 static int *generateAtomTypeTable(int ntypes) {
41  int *table = new int[ntypes*ntypes];
42  int ind = 0;
43  for (int i=0; i<ntypes; i++) {
44  for (int j=i; j<ntypes; j++) {
45  table[ntypes*i+j] = table[ntypes*j+i] = ind++;
46  }
47  }
48  return table;
49 }
50 
53 {
54  DebugM(3,"Constructing client\n");
55  comm = m;
57  kxmax = sp->pressureProfileEwaldX;
58  kymax = sp->pressureProfileEwaldY;
59  kzmax = sp->pressureProfileEwaldZ;
60 
61  ktot = (1+kxmax) * (2*kymax+1) * (2*kzmax+1);
63  pressureProfileSlabs = sp->pressureProfileSlabs;
64  numAtomTypes = sp->pressureProfileAtomTypes;
65  int nelements = 3*pressureProfileSlabs * (numAtomTypes*(numAtomTypes+1))/2;
66  pressureProfileData = new float[nelements];
67  reduction = ReductionMgr::Object()->willSubmit(
68  REDUCTIONS_PPROF_NONBONDED, nelements);
69 
70  // figure out who da masta be
71  numWorkingPes = (PatchMap::Object())->numNodesWithPatches();
72  masterNode = numWorkingPes - 1;
73 
74  recvCount = 0;
75  localAtoms = NULL;
76  localPartitions = NULL;
77 
78  expx = new float[kxmax+1];
79  expy = new float[kymax+1];
80  expz = new float[kzmax+1];
81 
82  if (CkMyPe() == masterNode) {
83  eiktotal = new float[2 * ktot * numAtomTypes];
84  memset(eiktotal, 0, 2 * ktot * numAtomTypes*sizeof(float));
85  } else {
86  eiktotal = NULL;
87  }
88  // space for exp(iky), k=-kymax, ..., kymax
89  eiky = new floatcomplex[2*kymax+1];
90  // space for exp(ikz), k=-kzmax, ..., kzmax
91  eikz = new floatcomplex[2*kzmax+1];
92  Qk = new float[3*ktot];
93 
94  gridsForAtomType = generateAtomTypeTable(numAtomTypes);
95 }
96 
98 {
99  delete reduction;
100  delete [] expx;
101  delete [] expy;
102  delete [] expz;
103  delete [] eiktotal;
104  delete [] eiky;
105  delete [] eikz;
106  delete [] pressureProfileData;
107  delete [] Qk;
108  delete [] gridsForAtomType;
109 
110  if (localAtoms) free(localAtoms);
111  if (localPartitions) free(localPartitions);
112 }
113 
116  // Skip computations if nothing to do.
117  if ( ! patchList[0].p->flags.doFullElectrostatics )
118  {
119  for (ap = ap.begin(); ap != ap.end(); ap++) {
120  CompAtom *x = (*ap).positionBox->open();
121  Results *r = (*ap).forceBox->open();
122  (*ap).positionBox->close(&x);
123  (*ap).forceBox->close(&r);
124  }
125  reduction->submit();
126  return;
127  }
128 
129 
130  lattice = patchList[0].p->lattice;
131  Vector o = lattice.origin();
132 
133  // recompute pressure profile cell parameters based on current lattice
134  pressureProfileThickness = lattice.c().z / pressureProfileSlabs;
135  pressureProfileMin = lattice.origin().z - 0.5*lattice.c().z;
136 
137  const BigReal coulomb_sqrt = sqrt( COULOMB * ComputeNonbondedUtil::scaling
139 
140  // get coordinates and store them
141  numLocalAtoms = 0;
142  for (ap = ap.begin(); ap != ap.end(); ap++) {
143  numLocalAtoms += (*ap).p->getNumAtoms();
144  }
145  localAtoms = (EwaldParticle *)realloc(localAtoms, numLocalAtoms*sizeof(EwaldParticle));
146  localPartitions = (int *)realloc(localPartitions, numLocalAtoms*sizeof(int));
147 
148  EwaldParticle *data_ptr = localAtoms;
149  int *part_ptr = localPartitions;
150 
151  for (ap = ap.begin(); ap != ap.end(); ap++) {
152  CompAtom *x = (*ap).positionBox->open();
153  // CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
154  Results *r = (*ap).forceBox->open();
155  if ( patchList[0].p->flags.doMolly ) {
156  (*ap).positionBox->close(&x);
157  x = (*ap).avgPositionBox->open();
158  }
159  int numAtoms = (*ap).p->getNumAtoms();
160 
161  for(int i=0; i<numAtoms; ++i) {
162  // wrap back to unit cell, centered on origin
163  Vector pos = x[i].position;
164  pos += lattice.wrap_delta(pos) - o;
165  *part_ptr++ = x[i].partition;
166  data_ptr->x = pos.x;
167  data_ptr->y = pos.y;
168  data_ptr->z = pos.z;
169  data_ptr->cg = coulomb_sqrt * x[i].charge;
170  ++data_ptr;
171  }
172  (*ap).positionBox->close(&x);
173  (*ap).forceBox->close(&r);
174  }
175 
176  // compute structure factor contribution from local atoms
177  // 2*ktot since charm++ uses float instead of floatcomplex
178  int msgsize = 2 * numAtomTypes * ktot;
179  ComputeEwaldMsg *msg = new (msgsize,0) ComputeEwaldMsg;
180  memset(msg->eik, 0, msgsize*sizeof(float));
181  compute_structurefactor(msg->eik);
182 
183  // send our partial sum
184  comm->sendComputeEwaldData(msg);
185 }
186 
188  // sum the data...
189  int nvecs = 2 * ktot * numAtomTypes;
190  for (int i=0; i<nvecs; i++) {
191  eiktotal[i] += msg->eik[i];
192  }
193  delete msg;
194  if (++recvCount == numWorkingPes) {
195  recvCount = 0;
196  int msgsize = 2 * ktot * numAtomTypes;
197  msg = new (msgsize,0) ComputeEwaldMsg;
198  memcpy(msg->eik, eiktotal, msgsize*sizeof(float));
199  memset(eiktotal, 0, msgsize*sizeof(float));
200  comm->sendComputeEwaldResults(msg);
201  }
202 }
203 
205  // receive total sum
206  computePprofile(msg->eik);
207  delete msg;
208  float scalefac = 1.0 / (2 * M_PI * lattice.volume());
209 
210  int nelements = 3*pressureProfileSlabs * (numAtomTypes*(numAtomTypes+1))/2;
211  for (int i=0; i<nelements; i++) {
212  reduction->item(i) += pressureProfileData[i] * scalefac;
213  }
214  reduction->submit();
215 }
216 
217 void ComputeEwald::compute_structurefactor(float *eik) {
218 
219  float recipx = lattice.a_r().x;
220  float recipy = lattice.b_r().y;
221  float recipz = lattice.c_r().z;
222  int i, j;
223  for (i=0; i<numLocalAtoms; i++) {
224  // compute exp(2 pi i r/L)
225  float krx = 2 * M_PI * localAtoms[i].x * recipx;
226  float kry = 2 * M_PI * localAtoms[i].y * recipy;
227  float krz = 2 * M_PI * localAtoms[i].z * recipz;
228  float cg = localAtoms[i].cg;
229  // sum structure factors for each atom type separately
230  const int offset = 2*ktot*localPartitions[i];
231 
232  floatcomplex eikx1(cos(krx), sin(krx));
233  floatcomplex eiky1(cos(kry), sin(kry));
234  floatcomplex eikz1(cos(krz), sin(krz));
235 
236  // store exp(2 pi i y j/Ly) for j = -kymax, ..., kymax
237  floatcomplex *ptr = eiky + kymax;
238  // n=0 case
239  ptr[0] = 1;
240  floatcomplex y(1,0);
241  for (j=1; j<= kymax; j++) {
242  y *= eiky1;
243  ptr[j] = y;
244  ptr[-j] = y.star();
245  }
246 
247  // store exp(2 pi i y j/Ly) for j = -kzmax, ..., kzmax
248  ptr = eikz + kzmax;
249  // n=0 case
250  ptr[0] = 1;
251  floatcomplex z(1,0);
252  for (j=1; j<= kzmax; j++) {
253  z *= eikz1;
254  ptr[j] = z;
255  ptr[-j] = z.star();
256  }
257 
258  // now loop over all k-vectors, computing S(k)
259  floatcomplex ex(cg);
260  int ind = offset; // index into eik
261  for (int kx=0; kx <= kxmax; kx++) {
262  for (int ky=0; ky <= 2*kymax; ky++) {
263  floatcomplex exy = eiky[ky];
264  exy *= ex;
265  const int max = 2*kzmax;
266  const float exyr = exy.r;
267  const float exyi = exy.i;
268  const float *eikzd = (const float *)eikz;
269 #pragma vector always
270  for (int kz=0; kz <= max; kz++) {
271  float ezr = *eikzd++;
272  float ezi = *eikzd++;
273  float eikr = ezr*exyr - ezi*exyi;
274  float eiki = ezr*exyi + ezi*exyr;
275 
276  // add this contribution to each grid to which the atom belongs
277  eik[ind ] += eikr;
278  eik[ind+1] += eiki;
279 
280  // next k vector
281  ind += 2;
282  }
283  }
284  ex *= eikx1;
285  }
286  }
287 }
288 
289 // compute exp(-k^2/4 kappa^2) for k=2pi*recip*n, n=0...K inclusive
290 static void init_exp(float *xp, int K, float recip, float kappa) {
291  float piob = M_PI / kappa;
292  piob *= piob;
293  float fac = -piob*recip*recip;
294  for (int i=0; i<= K; i++)
295  xp[i] = exp(fac*i*i);
296 }
297 
298 void ComputeEwald::computePprofile(const float *eik) const {
299  float recipx = lattice.a_r().x;
300  float recipy = lattice.b_r().y;
301  float recipz = lattice.c_r().z;
302 
303  init_exp(expx, kxmax, recipx, kappa);
304  init_exp(expy, kymax, recipy, kappa);
305  init_exp(expz, kzmax, recipz, kappa);
306 
307  //float energy = 0;
308  float piob = M_PI / kappa;
309  piob *= piob;
310 
311  // compute exp(-pi^2 m^2 / B^2)/m^2
312  int ind = 0;
313  for (int kx=0; kx <= kxmax; kx++) {
314  float m11 = recipx * kx;
315  m11 *= m11;
316  float xfac = expx[kx] * (kx ? 2 : 1);
317  for (int ky=-kymax; ky <= kymax; ky++) {
318  float m22 = recipy * ky;
319  m22 *= m22;
320  float xyfac = expy[abs(ky)] * xfac;
321  for (int kz=-kzmax; kz <= kzmax; kz++) {
322  float m33 = recipz * kz;
323  m33 *= m33;
324  float msq = m11 + m22 + m33;
325  float imsq = msq ? 1.0 / msq : 0;
326  float fac = expz[abs(kz)] * xyfac * imsq;
327 
328  float pfac = 2*(imsq + piob);
329  Qk[ind++] = fac*(1-pfac*m11);
330  Qk[ind++] = fac*(1-pfac*m22);
331  Qk[ind++] = fac*(1-pfac*m33);
332  }
333  }
334  }
335 
336  const int nslabs = pressureProfileSlabs;
337 
338  int nelements = 3*nslabs * (numAtomTypes*(numAtomTypes+1))/2;
339  memset(pressureProfileData, 0, nelements*sizeof(float));
340  int i, j;
341  for (i=0; i<numLocalAtoms; i++) {
342 
343  float krx = 2 * M_PI * localAtoms[i].x * recipx;
344  float kry = 2 * M_PI * localAtoms[i].y * recipy;
345  float krz = 2 * M_PI * localAtoms[i].z * recipz;
346  float cg = localAtoms[i].cg;
347  int atype = localPartitions[i];
348  const int *grids = gridsForAtomType+atype*numAtomTypes;
349 
350  // determine the slab where this particle is located
351  int slab = (int)floor((localAtoms[i].z - pressureProfileMin)/pressureProfileThickness);
352  if (slab < 0) slab += nslabs;
353  else if (slab >= nslabs) slab -= nslabs;
354  float *pprofptr = pressureProfileData + 3*slab;
355 
356  floatcomplex eikx1(cos(krx), sin(krx));
357  floatcomplex eiky1(cos(kry), sin(kry));
358  floatcomplex eikz1(cos(krz), sin(krz));
359 
360  // store exp(2 pi i y j/Ly) for j = -kymax, ..., kymax
361  floatcomplex *ptr = eiky + kymax;
362  // n=0 case
363  ptr[0] = 1;
364  floatcomplex y(1,0);
365  for (j=1; j<= kymax; j++) {
366  y *= eiky1;
367  ptr[j] = y;
368  ptr[-j] = y.star();
369  }
370 
371  // store exp(2 pi i y j/Ly) for j = -kzmax, ..., kzmax
372  ptr = eikz + kzmax;
373  // n=0 case
374  ptr[0] = 1;
375  floatcomplex z(1,0);
376  for (j=1; j<= kzmax; j++) {
377  z *= eikz1;
378  ptr[j] = z;
379  ptr[-j] = z.star();
380  }
381 
382  int ind = 0;
383  const float *Qkptr = Qk;
384  floatcomplex ex(cg);
385  const float *eikptr = eik;
386  for (int kx=0; kx <= kxmax; kx++) {
387  for (int ky=0; ky <= 2*kymax; ky++) {
388  floatcomplex exy = eiky[ky];
389  exy *= ex;
390  const int kzmax2 = 2*kzmax+1;
391  const float exyr = exy.r;
392  const float exyi = exy.i;
393  const float *eikzptr = (const float *)eikz;
394 #pragma vector always
395  for (int kz=0; kz < kzmax2; kz++) {
396  float ezr = *eikzptr++;
397  float ezi = *eikzptr++;
398  float exyzr = ezr * exyr - ezi * exyi;
399  float exyzi = ezr * exyi + ezi * exyr;
400 
401  // exyz holds exp(ikr) for this atom
402  // loop over atom types, adding contribution from each
403  // Re[exp(-ikr) * S(k)]
404  // add this contribution to each grid the atom belongs to
405  for (int igrid=0; igrid<numAtomTypes; igrid++) {
406  float eikr = eikptr[igrid*2*ktot ];
407  float eiki = eikptr[igrid*2*ktot+1];
408  int grid = grids[igrid];
409  int offset = 3*nslabs*grid;
410 
411  float E = exyzr * eikr + exyzi * eiki;
412  float vx = Qkptr[0] * E;
413  float vy = Qkptr[1] * E;
414  float vz = Qkptr[2] * E;
415  pprofptr[offset ] += vx;
416  pprofptr[offset+1] += vy;
417  pprofptr[offset+2] += vz;
418  }
419  // next k vector
420  eikptr += 2;
421  Qkptr += 3;
422  }
423  }
424  ex *= eikx1;
425  }
426  }
427 }
428 
static Node * Object()
Definition: Node.h:86
static void init_exp(float *xp, int K, float recip, float kappa)
Definition: ComputeEwald.C:290
unsigned char partition
Definition: NamdTypes.h:56
#define M_PI
Definition: GoMolecule.C:39
void recvData(ComputeEwaldMsg *)
Definition: ComputeEwald.C:187
void recvResults(ComputeEwaldMsg *)
Definition: ComputeEwald.C:204
int pressureProfileEwaldX
Vector a_r() const
Definition: Lattice.h:268
int ComputeID
Definition: NamdTypes.h:183
static int * generateAtomTypeTable(int ntypes)
Definition: ComputeEwald.C:40
static PatchMap * Object()
Definition: PatchMap.h:27
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
ComputeHomePatchList patchList
Vector c_r() const
Definition: Lattice.h:270
Vector wrap_delta(const Position &pos1) const
Definition: Lattice.h:206
#define COULOMB
Definition: common.h:44
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
ComputeEwald(ComputeID, ComputeMgr *)
Definition: ComputeEwald.C:51
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
int pressureProfileSlabs
Vector origin() const
Definition: Lattice.h:262
Vector b_r() const
Definition: Lattice.h:269
Charge charge
Definition: NamdTypes.h:54
ResizeArrayIter< T > end(void) const
void sendComputeEwaldData(ComputeEwaldMsg *)
Definition: ComputeMgr.C:1325
gridSize z
BigReal x
Definition: Vector.h:66
BigReal volume(void) const
Definition: Lattice.h:277
int pressureProfileEwaldY
int pressureProfileAtomTypes
BigReal y
Definition: Vector.h:66
gridSize y
void submit(void)
Definition: ReductionMgr.h:323
gridSize x
virtual ~ComputeEwald()
Definition: ComputeEwald.C:97
void sendComputeEwaldResults(ComputeEwaldMsg *)
Definition: ComputeMgr.C:1348
ResizeArrayIter< T > begin(void) const
Vector c() const
Definition: Lattice.h:254
double BigReal
Definition: common.h:112
int pressureProfileEwaldZ