NAMD
ComputeMsmSerial.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "Node.h"
9 #include "PatchMap.h"
10 #include "PatchMap.inl"
11 #include "AtomMap.h"
12 #include "ComputeMsmSerial.h"
13 #include "ComputeMsmSerialMgr.decl.h"
14 #include "PatchMgr.h"
15 #include "Molecule.h"
16 #include "ReductionMgr.h"
17 #include "ComputeMgr.h"
18 #include "ComputeMgr.decl.h"
19 // #define DEBUGM
20 #define MIN_DEBUG_LEVEL 3
21 #include "Debug.h"
22 #include "SimParameters.h"
23 #include "WorkDistrib.h"
24 #include "varsizemsg.h"
25 #include "msm.h"
26 #include <stdlib.h>
27 #include <stdio.h>
28 #include <errno.h>
29 
30 
33  float charge;
34  int id;
35 };
36 
38 
39 class MsmSerialCoordMsg : public CMessage_MsmSerialCoordMsg {
40 public:
42  int numAtoms;
45 };
46 
47 class MsmSerialForceMsg : public CMessage_MsmSerialForceMsg {
48 public:
50  BigReal virial[3][3];
52 };
53 
54 class ComputeMsmSerialMgr : public CBase_ComputeMsmSerialMgr {
55 public:
58 
59  void setCompute(ComputeMsmSerial *c) { msmCompute = c; }
60 
63 
64 private:
65  CProxy_ComputeMsmSerialMgr msmProxy;
66  ComputeMsmSerial *msmCompute;
67 
68  int numSources;
69  int numArrived;
70  MsmSerialCoordMsg **coordMsgs;
71  int numAtoms;
72  ComputeMsmSerialAtom *coord;
73  MsmSerialForce *force;
74  MsmSerialForceMsg *oldmsg;
75 
76  NL_Msm *msmsolver;
77  double *msmcoord;
78  double *msmforce;
79 };
80 
82  msmProxy(thisgroup), msmCompute(0), numSources(0), numArrived(0),
83  coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0),
84  msmsolver(0), msmcoord(0), msmforce(0)
85 {
86  CkpvAccess(BOCclass_group).computeMsmSerialMgr = thisgroup;
87 }
88 
90 {
91  for (int i=0; i < numSources; i++) delete coordMsgs[i];
92  delete [] coordMsgs;
93  delete [] coord;
94  delete [] force;
95  delete oldmsg;
96  if (msmsolver) NL_msm_destroy(msmsolver);
97  if (msmcoord) delete[] msmcoord;
98  if (msmforce) delete[] msmforce;
99 }
100 
103 {
104  CProxy_ComputeMsmSerialMgr::ckLocalBranch(
105  CkpvAccess(BOCclass_group).computeMsmSerialMgr)->setCompute(this);
107 }
108 
110 {
111 }
112 
114 {
116 
117  // Skip computations if nothing to do.
118  if ( ! patchList[0].p->flags.doFullElectrostatics )
119  {
120  for (ap = ap.begin(); ap != ap.end(); ap++) {
121  CompAtom *x = (*ap).positionBox->open();
122  Results *r = (*ap).forceBox->open();
123  (*ap).positionBox->close(&x);
124  (*ap).forceBox->close(&r);
125  }
126  reduction->submit();
127  return;
128  }
129 
130  // allocate message
131  int numLocalAtoms = 0;
132  for (ap = ap.begin(); ap != ap.end(); ap++) {
133  numLocalAtoms += (*ap).p->getNumAtoms();
134  }
135 
136  MsmSerialCoordMsg *msg = new (numLocalAtoms, 0) MsmSerialCoordMsg;
137  msg->sourceNode = CkMyPe();
138  msg->numAtoms = numLocalAtoms;
139  msg->lattice = patchList[0].p->flags.lattice;
140  ComputeMsmSerialAtom *data_ptr = msg->coord;
141 
142  // get positions
143  for (ap = ap.begin(); ap != ap.end(); ap++) {
144  CompAtom *x = (*ap).positionBox->open();
145  CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
146  if ( patchList[0].p->flags.doMolly ) {
147  (*ap).positionBox->close(&x);
148  x = (*ap).avgPositionBox->open();
149  }
150  int numAtoms = (*ap).p->getNumAtoms();
151 
152  for(int i=0; i < numAtoms; i++)
153  {
154  data_ptr->position = x[i].position;
155  data_ptr->charge = x[i].charge;
156  data_ptr->id = xExt[i].id;
157  ++data_ptr;
158  }
159 
160  if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
161  else { (*ap).positionBox->close(&x); }
162  }
163 
164  CProxy_ComputeMsmSerialMgr msmProxy(
165  CkpvAccess(BOCclass_group).computeMsmSerialMgr);
166  msmProxy[0].recvCoord(msg);
167 }
168 
169 
171  Vector &u, Vector &v, Vector &w, Vector &c,
172  Vector &ru, Vector &rv, Vector &rw,
173  int isperiodic, int numatoms, const ComputeMsmSerialAtom *coord,
174  double padding, double gridspacing) {
175  const double NL_ZERO_TOLERANCE = 1e-4;
176  double xlen = 1, inv_xlen = 1; /* don't rescale periodic directions */
177  double ylen = 1, inv_ylen = 1;
178  double zlen = 1, inv_zlen = 1;
179  double ulen_1 = u.rlength();
180  double vlen_1 = v.rlength();
181  double wlen_1 = w.rlength();
182  double rpadx = padding * ulen_1; /* padding distance in recip space */
183  double rpady = padding * vlen_1;
184  double rpadz = padding * wlen_1;
185  double rhx = gridspacing * ulen_1; /* grid spacing in recip space */
186  double rhy = gridspacing * vlen_1;
187  double rhz = gridspacing * wlen_1;
188  Vector s, d;
189  Vector rc = 0; // don't move center along periodic directions
190  double xmin, xmax, ymin, ymax, zmin, zmax;
191  int is_periodic_x = (isperiodic & NL_MSM_PERIODIC_VEC1);
192  int is_periodic_y = (isperiodic & NL_MSM_PERIODIC_VEC2);
193  int is_periodic_z = (isperiodic & NL_MSM_PERIODIC_VEC3);
194  int i;
195 
196  /* affine linear transformation of coordinates to reciprocal space,
197  * where periodic vector directions map into [-0.5, 0.5) */
198  //printf("*** center=%.4f %.4f %.4f\n", c.x, c.y, c.z);
199  d = coord[0].position - c;
200  s.x = ru * d; // dot product
201  s.y = rv * d;
202  s.z = rw * d;
203  xmin = xmax = s.x;
204  ymin = ymax = s.y;
205  zmin = zmax = s.z;
206  for (i = 1; i < numatoms; i++) {
207  d = coord[i].position - c;
208  s.x = ru * d; // dot product
209  s.y = rv * d;
210  s.z = rw * d;
211  if (s.x < xmin) xmin = s.x;
212  else if (s.x > xmax) xmax = s.x;
213  if (s.y < ymin) ymin = s.y;
214  else if (s.y > ymax) ymax = s.y;
215  if (s.z < zmin) zmin = s.z;
216  else if (s.z > zmax) zmax = s.z;
217  }
218 #if 0
219  printf("*** xmin=%.4f xmax=%.4f\n", xmin, xmax);
220  printf("*** ymin=%.4f ymax=%.4f\n", ymin, ymax);
221  printf("*** zmin=%.4f zmax=%.4f\n", zmin, zmax);
222 #endif
223 
224  if ( ! is_periodic_x) {
225  xmax += rpadx; /* pad the edges */
226  xmin -= rpadx;
227  if (rhx > 0) { /* restrict center to rhx lattice points */
228  double mupper = ceil(xmax / (2*rhx));
229  double mlower = floor(xmin / (2*rhx));
230  xmax = 2*rhx*mupper;
231  xmin = 2*rhx*mlower;
232  }
233  rc.x = 0.5*(xmin + xmax);
234  xlen = xmax - xmin;
235  if (xlen < NL_ZERO_TOLERANCE) {
236  xlen = 1; /* leave scaling unchanged */
237  }
238  else {
239  inv_xlen = 1/xlen;
240  }
241  }
242 
243  if ( ! is_periodic_y) {
244  ymax += rpady; /* pad the edges */
245  ymin -= rpady;
246  if (rhy > 0) { /* restrict center to rhy lattice points */
247  double mupper = ceil(ymax / (2*rhy));
248  double mlower = floor(ymin / (2*rhy));
249  ymax = 2*rhy*mupper;
250  ymin = 2*rhy*mlower;
251  }
252  rc.y = 0.5*(ymin + ymax);
253  ylen = ymax - ymin;
254  if (ylen < NL_ZERO_TOLERANCE) {
255  ylen = 1; /* leave scaling unchanged */
256  }
257  else {
258  inv_ylen = 1/ylen;
259  }
260  }
261 
262  if ( ! is_periodic_z) {
263  zmax += rpadz; /* pad the edges */
264  zmin -= rpadz;
265  if (rhz > 0) { /* restrict center to rhz lattice points */
266  double mupper = ceil(zmax / (2*rhz));
267  double mlower = floor(zmin / (2*rhz));
268  zmax = 2*rhz*mupper;
269  zmin = 2*rhz*mlower;
270  }
271  rc.z = 0.5*(zmin + zmax);
272  zlen = zmax - zmin;
273  if (zlen < NL_ZERO_TOLERANCE) {
274  zlen = 1; /* leave scaling unchanged */
275  }
276  else {
277  inv_zlen = 1/zlen;
278  }
279  }
280 
281 #if 0
282  printf("xmin=%g xmax=%g\n", xmin, xmax);
283  printf("ymin=%g ymax=%g\n", ymin, ymax);
284  printf("zmin=%g zmax=%g\n", zmin, zmax);
285  printf("xlen=%g ylen=%g zlen=%g\n", xlen, ylen, zlen);
286  printf("rc_x=%g rc_y=%g rc_z=%g\n", rc.x, rc.y, rc.z);
287 #endif
288 
289  /* transform new center back to real space, then rescale basis vectors */
290  c.x = (u.x*rc.x + v.x*rc.y + w.x*rc.z) + c.x;
291  c.y = (u.y*rc.x + v.y*rc.y + w.y*rc.z) + c.y;
292  c.z = (u.z*rc.x + v.z*rc.y + w.z*rc.z) + c.z;
293 
294 #if 0
295  printf("c_x=%g c_y=%g c_z=%g\n", c.x, c.y, c.z);
296 #endif
297 
298  u *= xlen;
299  v *= ylen;
300  w *= zlen;
301 
302  ru *= inv_xlen;
303  rv *= inv_ylen;
304  rw *= inv_zlen;
305 }
306 
307 
309  if ( ! numSources ) {
310  numSources = (PatchMap::Object())->numNodesWithPatches();
311  coordMsgs = new MsmSerialCoordMsg*[numSources];
312  for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
313  numArrived = 0;
314  numAtoms = Node::Object()->molecule->numAtoms;
315  coord = new ComputeMsmSerialAtom[numAtoms];
316  force = new MsmSerialForce[numAtoms];
317  }
318 
319  int i;
320  for ( i=0; i < msg->numAtoms; ++i ) {
321  coord[msg->coord[i].id] = msg->coord[i];
322  }
323 
324  coordMsgs[numArrived] = msg;
325  ++numArrived;
326 
327  if ( numArrived < numSources ) return;
328  numArrived = 0;
329 
330  // ALL DATA ARRIVED --- CALCULATE FORCES
331  Lattice lattice = msg->lattice;
333 
334  double energy = 0;
335  double virial[3][3];
336 
337  int rc = 0; // return code
338 
339  if ( ! msmsolver ) {
340  //
341  // setup MSM solver
342  //
343  msmsolver = NL_msm_create();
344  if ( ! msmsolver ) NAMD_die("unable to create MSM solver");
345  double dielectric = simParams->dielectric;
346  double cutoff = simParams->cutoff;
347  double gridspacing = simParams->MSMGridSpacing;
348  double padding = simParams->MSMPadding;
349  int approx = simParams->MSMApprox;
350  int split = simParams->MSMSplit;
351  int nlevels = simParams->MSMLevels;
352  int msmflags = 0;
353  msmflags |= (lattice.a_p() ? NL_MSM_PERIODIC_VEC1 : 0);
354  msmflags |= (lattice.b_p() ? NL_MSM_PERIODIC_VEC2 : 0);
355  msmflags |= (lattice.c_p() ? NL_MSM_PERIODIC_VEC3 : 0);
356  msmflags |= NL_MSM_COMPUTE_LONG_RANGE; // compute only long-range part
357  //msmflags |= NL_MSM_COMPUTE_ALL;
358  //printf("msmflags = %x\n", msmflags);
359  rc = NL_msm_configure(msmsolver, gridspacing, approx, split, nlevels);
360  if (rc) NAMD_die("unable to configure MSM solver");
361  Vector u=lattice.a(), v=lattice.b(), w=lattice.c(), c=lattice.origin();
362  Vector ru=lattice.a_r(), rv=lattice.b_r(), rw=lattice.c_r();
363  if ((msmflags & NL_MSM_PERIODIC_ALL) != NL_MSM_PERIODIC_ALL) {
364  // called only if there is some non-periodic boundary
365  int isperiodic = (msmflags & NL_MSM_PERIODIC_ALL);
366  //printf("calling rescale\n");
367  rescale_nonperiodic_cell(u, v, w, c, ru, rv, rw,
368  isperiodic, numAtoms, coord, padding, gridspacing);
369  }
370  double vec1[3], vec2[3], vec3[3], center[3];
371  vec1[0] = u.x;
372  vec1[1] = u.y;
373  vec1[2] = u.z;
374  vec2[0] = v.x;
375  vec2[1] = v.y;
376  vec2[2] = v.z;
377  vec3[0] = w.x;
378  vec3[1] = w.y;
379  vec3[2] = w.z;
380  center[0] = c.x;
381  center[1] = c.y;
382  center[2] = c.z;
383 #if 0
384  printf("dielectric = %g\n", dielectric);
385  printf("vec1 = %g %g %g\n", vec1[0], vec1[1], vec1[2]);
386  printf("vec2 = %g %g %g\n", vec2[0], vec2[1], vec2[2]);
387  printf("vec3 = %g %g %g\n", vec3[0], vec3[1], vec3[2]);
388  printf("center = %g %g %g\n", center[0], center[1], center[2]);
389  printf("cutoff = %g\n", cutoff);
390  printf("numatoms = %d\n", numAtoms);
391 #endif
392  rc = NL_msm_setup(msmsolver, cutoff, vec1, vec2, vec3, center, msmflags);
393  if (rc) NAMD_die("unable to set up MSM solver");
394  msmcoord = new double[4*numAtoms];
395  msmforce = new double[3*numAtoms];
396  if (msmcoord==0 || msmforce==0) NAMD_die("can't allocate MSM atom buffers");
397  // scale charges - these won't change
398  double celec = sqrt(COULOMB / simParams->dielectric);
399  for (i = 0; i < numAtoms; i++) {
400  msmcoord[4*i+3] = celec * coord[i].charge;
401  }
402  }
403 
404  // evaluate long-range MSM forces
405  for (i = 0; i < numAtoms; i++) {
406  msmcoord[4*i ] = coord[i].position.x;
407  msmcoord[4*i+1] = coord[i].position.y;
408  msmcoord[4*i+2] = coord[i].position.z;
409  }
410  for (i = 0; i < numAtoms; i++) {
411  msmforce[3*i ] = 0;
412  msmforce[3*i+1] = 0;
413  msmforce[3*i+2] = 0;
414  }
415  rc = NL_msm_compute_force(msmsolver, msmforce, &energy, msmcoord, numAtoms);
416  if (rc) NAMD_die("error evaluating MSM forces");
417  for (i = 0; i < numAtoms; i++) {
418  force[i].x = msmforce[3*i ];
419  force[i].y = msmforce[3*i+1];
420  force[i].z = msmforce[3*i+2];
421  }
422  // MSM does not yet calculate virial
423  for (int k=0; k < 3; k++) {
424  for (int l=0; l < 3; l++) {
425  virial[k][l] = 0;
426  }
427  }
428 
429  // distribute forces
430  for (int j=0; j < numSources; j++) {
431  MsmSerialCoordMsg *cmsg = coordMsgs[j];
432  coordMsgs[j] = 0;
433  MsmSerialForceMsg *fmsg = new (cmsg->numAtoms, 0) MsmSerialForceMsg;
434 
435  for (int i=0; i < cmsg->numAtoms; i++) {
436  fmsg->force[i] = force[cmsg->coord[i].id];
437  }
438 
439  if ( ! j ) { // set virial and energy only for first message
440  fmsg->energy = energy;
441  for (int k=0; k < 3; k++) {
442  for (int l=0; l < 3; l++) {
443  fmsg->virial[k][l] = virial[k][l];
444  }
445  }
446  }
447  else { // set other messages to zero, add into reduction only once
448  fmsg->energy = 0;
449  for (int k=0; k < 3; k++) {
450  for (int l=0; l < 3; l++) {
451  fmsg->virial[k][l] = 0;
452  }
453  }
454  }
455 
456  msmProxy[cmsg->sourceNode].recvForce(fmsg);
457  delete cmsg;
458  }
459 }
460 
462  msmCompute->saveResults(msg);
463  delete oldmsg;
464  oldmsg = msg;
465 }
466 
468 {
470 
471  MsmSerialForce *results_ptr = msg->force;
472 
473  // add in forces
474  for (ap = ap.begin(); ap != ap.end(); ap++) {
475  Results *r = (*ap).forceBox->open();
476  Force *f = r->f[Results::slow];
477  int numAtoms = (*ap).p->getNumAtoms();
478 
479  for(int i=0; i<numAtoms; ++i) {
480  f[i].x += results_ptr->x;
481  f[i].y += results_ptr->y;
482  f[i].z += results_ptr->z;
483  ++results_ptr;
484  }
485 
486  (*ap).forceBox->close(&r);
487  }
488 
489  reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += msg->energy;
490  reduction->item(REDUCTION_VIRIAL_SLOW_XX) += msg->virial[0][0];
491  reduction->item(REDUCTION_VIRIAL_SLOW_XY) += msg->virial[0][1];
492  reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += msg->virial[0][2];
493  reduction->item(REDUCTION_VIRIAL_SLOW_YX) += msg->virial[1][0];
494  reduction->item(REDUCTION_VIRIAL_SLOW_YY) += msg->virial[1][1];
495  reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += msg->virial[1][2];
496  reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += msg->virial[2][0];
497  reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += msg->virial[2][1];
498  reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += msg->virial[2][2];
499  reduction->submit();
500 }
501 
502 #include "ComputeMsmSerialMgr.def.h"
503 
static Node * Object()
Definition: Node.h:86
void setCompute(ComputeMsmSerial *c)
Vector a_r() const
Definition: Lattice.h:268
BigReal virial[3][3]
int ComputeID
Definition: NamdTypes.h:183
static PatchMap * Object()
Definition: PatchMap.h:27
Definition: Vector.h:64
Force MsmSerialForce
SimParameters * simParameters
Definition: Node.h:178
ComputeHomePatchList patchList
Vector c_r() const
Definition: Lattice.h:270
static int numatoms
Definition: ScriptTcl.C:64
#define COULOMB
Definition: common.h:44
BigReal & item(int i)
Definition: ReductionMgr.h:312
BigReal z
Definition: Vector.h:66
ComputeMsmSerialAtom * coord
Position position
Definition: NamdTypes.h:53
NL_Msm * NL_msm_create(void)
Definition: msm.c:6
int NL_msm_configure(NL_Msm *pm, double gridspacing, int approx, int split, int nlevels)
Definition: msm.c:36
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
void recvForce(MsmSerialForceMsg *)
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
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
MsmSerialForce * force
BigReal rlength(void)
Definition: Vector.h:177
void saveResults(MsmSerialForceMsg *)
Force * f[maxNumForces]
Definition: PatchTypes.h:67
int NL_msm_compute_force(NL_Msm *pm, double *felec, double *uelec, const double *atom, int natoms)
Definition: msm.c:56
BigReal x
Definition: Vector.h:66
virtual ~ComputeMsmSerial()
int numAtoms
Definition: Molecule.h:556
BigReal MSMGridSpacing
void NAMD_die(const char *err_msg)
Definition: common.C:83
std::vector< std::string > split(const std::string &text, std::string delimiter)
Definition: MoleculeQM.C:73
void recvCoord(MsmSerialCoordMsg *)
ComputeMsmSerial(ComputeID c)
#define simParams
Definition: Output.C:127
static void rescale_nonperiodic_cell(Vector &u, Vector &v, Vector &w, Vector &c, Vector &ru, Vector &rv, Vector &rw, int isperiodic, int numatoms, const ComputeMsmSerialAtom *coord, double padding, double gridspacing)
BigReal y
Definition: Vector.h:66
Vector b() const
Definition: Lattice.h:253
BigReal dielectric
void submit(void)
Definition: ReductionMgr.h:323
void NL_msm_destroy(NL_Msm *pm)
Definition: msm.c:28
int b_p() const
Definition: Lattice.h:274
BigReal MSMPadding
gridSize x
int a_p() const
Definition: Lattice.h:273
Molecule * molecule
Definition: Node.h:176
Vector a() const
Definition: Lattice.h:252
int NL_msm_setup(NL_Msm *msm, double cutoff, double cellvec1[3], double cellvec2[3], double cellvec3[3], double cellcenter[3], int msmflags)
Definition: msm_setup.c:62
ResizeArrayIter< T > begin(void) const
Vector c() const
Definition: Lattice.h:254
double BigReal
Definition: common.h:112
int c_p() const
Definition: Lattice.h:275