NAMD
ComputeFmmSerial.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 "ComputeFmmSerial.h"
13 #include "ComputeFmmSerialMgr.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 <stdlib.h>
26 #include <stdio.h>
27 #include <errno.h>
28 
29 
30 #ifdef FMM_SOLVER
31 
32 // Calculate FMM:
33 //
34 // ufmmlap library from J. Huang - http://fastmultipole.org/Main/FMMSuite/
35 // (ufmmlap = Uniform FMM Laplace Solver)
36 //
37 // The positions must translated and scaled to be within the unit
38 // cube centered at the origin, i.e., [-1/2,1/2]^3. The potentials
39 // and forces must then also be scaled.
40 //
41 extern "C" void fmmlap_uni_(
42  int *natoms, // input: number of particles
43  double zat[][3], // input: particle positions, length natoms
44  double charge[], // input: particle charges, length natoms
45  double pot[], // output: particle potentials, length natoms
46  double field[][3], // output: particle forces, length natoms
47  int *nlev, // input: number of oct-tree levels
48  int *ier // output: error code
49  );
50 
51 #endif // FMM_SOLVER
52 
53 
56  float charge;
57  int id;
58 };
59 
61 
62 class FmmSerialCoordMsg : public CMessage_FmmSerialCoordMsg {
63 public:
65  int numAtoms;
68 };
69 
70 class FmmSerialForceMsg : public CMessage_FmmSerialForceMsg {
71 public:
73  BigReal virial[3][3];
75 };
76 
77 class ComputeFmmSerialMgr : public CBase_ComputeFmmSerialMgr {
78 public:
81 
82  void setCompute(ComputeFmmSerial *c) { fmmCompute = c; }
83 
86 
87 private:
88  CProxy_ComputeFmmSerialMgr fmmProxy;
89  ComputeFmmSerial *fmmCompute;
90 
91  int numSources;
92  int numArrived;
93  FmmSerialCoordMsg **coordMsgs;
94  int numAtoms;
95  ComputeFmmSerialAtom *coord;
96  FmmSerialForce *force;
97  FmmSerialForceMsg *oldmsg;
98 
99  int fmmsolver;
100  int nlevels;
101  double scaling;
102  Vector center;
103  double *chargebuffer;
104  double *epotbuffer;
105  double (*posbuffer)[3];
106  double (*forcebuffer)[3];
107 };
108 
110  fmmProxy(thisgroup), fmmCompute(0), numSources(0), numArrived(0),
111  coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0),
112  fmmsolver(0), nlevels(0), scaling(0), center(0),
113  chargebuffer(0), epotbuffer(0), posbuffer(0), forcebuffer(0)
114 {
115  CkpvAccess(BOCclass_group).computeFmmSerialMgr = thisgroup;
116 }
117 
119 {
120  for (int i=0; i < numSources; i++) delete coordMsgs[i];
121  delete [] coordMsgs;
122  delete [] coord;
123  delete [] force;
124  delete oldmsg;
125  if (chargebuffer) delete[] chargebuffer;
126  if (epotbuffer) delete[] epotbuffer;
127  if (posbuffer) delete[] posbuffer;
128  if (forcebuffer) delete[] forcebuffer;
129 }
130 
133 {
134  CProxy_ComputeFmmSerialMgr::ckLocalBranch(
135  CkpvAccess(BOCclass_group).computeFmmSerialMgr)->setCompute(this);
137 }
138 
140 {
141 }
142 
144 {
146 
147  // Skip computations if nothing to do.
148  if ( ! patchList[0].p->flags.doFullElectrostatics )
149  {
150  for (ap = ap.begin(); ap != ap.end(); ap++) {
151  CompAtom *x = (*ap).positionBox->open();
152  Results *r = (*ap).forceBox->open();
153  (*ap).positionBox->close(&x);
154  (*ap).forceBox->close(&r);
155  }
156  reduction->submit();
157  return;
158  }
159 
160  // allocate message
161  int numLocalAtoms = 0;
162  for (ap = ap.begin(); ap != ap.end(); ap++) {
163  numLocalAtoms += (*ap).p->getNumAtoms();
164  }
165 
166  FmmSerialCoordMsg *msg = new (numLocalAtoms, 0) FmmSerialCoordMsg;
167  msg->sourceNode = CkMyPe();
168  msg->numAtoms = numLocalAtoms;
169  msg->lattice = patchList[0].p->flags.lattice;
170  ComputeFmmSerialAtom *data_ptr = msg->coord;
171 
172  // get positions
173  for (ap = ap.begin(); ap != ap.end(); ap++) {
174  CompAtom *x = (*ap).positionBox->open();
175  CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
176  if ( patchList[0].p->flags.doMolly ) {
177  (*ap).positionBox->close(&x);
178  x = (*ap).avgPositionBox->open();
179  }
180  int numAtoms = (*ap).p->getNumAtoms();
181 
182  for(int i=0; i < numAtoms; i++)
183  {
184  data_ptr->position = x[i].position;
185  data_ptr->charge = x[i].charge;
186  data_ptr->id = xExt[i].id;
187  ++data_ptr;
188  }
189 
190  if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
191  else { (*ap).positionBox->close(&x); }
192  }
193 
194  CProxy_ComputeFmmSerialMgr fmmProxy(
195  CkpvAccess(BOCclass_group).computeFmmSerialMgr);
196  fmmProxy[0].recvCoord(msg);
197 }
198 
199 
201  int i;
202  if ( ! numSources ) {
203  numSources = (PatchMap::Object())->numNodesWithPatches();
204  coordMsgs = new FmmSerialCoordMsg*[numSources];
205  for ( i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
206  numArrived = 0;
207  numAtoms = Node::Object()->molecule->numAtoms;
208  coord = new ComputeFmmSerialAtom[numAtoms];
209  force = new FmmSerialForce[numAtoms];
210  }
211 
212  for ( i=0; i < msg->numAtoms; ++i ) {
213  coord[msg->coord[i].id] = msg->coord[i];
214  }
215 
216  coordMsgs[numArrived] = msg;
217  ++numArrived;
218 
219  if ( numArrived < numSources ) return;
220  numArrived = 0;
221 
222  // ALL DATA ARRIVED --- CALCULATE FORCES
223  Lattice lattice = msg->lattice;
225 
226  double energy = 0;
227  double virial[3][3] = { 0 }; // FMM does not calculate virial
228 
229  int rc = 0; // return code
230 
231  if ( ! fmmsolver ) {
232  //
233  // setup FMM solver
234  //
235 
236  // check boundary conditions
237  if (lattice.a_p() || lattice.b_p() || lattice.c_p()) {
238  NAMD_die("FMM solver requires non-periodic boundaries");
239  }
240 
241  // setup number of levels
242  if (simParams->FMMLevels > 0) {
243  // explicitly set in config file
244  nlevels = simParams->FMMLevels;
245  }
246  else {
247  // otherwise estimate number of levels as round(log_8(numAtoms))
248  nlevels = (int) floor(log((double)numAtoms) / log(8.) + 0.5);
249  }
250 
251  // find bounding cube length
252  if (numAtoms <= 0) {
253  NAMD_die("setting up FMM with no atoms");
254  }
255  Vector min, max;
256  min = max = coord[0].position;
257  for (i = 1; i < numAtoms; i++) {
258  Vector r = coord[i].position;
259  if (r.x < min.x) min.x = r.x;
260  else if (r.x > max.x) max.x = r.x;
261  if (r.y < min.y) min.y = r.y;
262  else if (r.y > max.y) max.y = r.y;
263  if (r.z < min.z) min.z = r.z;
264  else if (r.z > max.z) max.z = r.z;
265  }
266  double padding = simParams->FMMPadding;
267  if (padding <= 0) padding = 0.01; // pad by at least a delta margin
268  min -= padding;
269  max += padding;
270  double len = max.x - min.x;
271  if (len < max.y - min.y) len = max.y - min.y;
272  if (len < max.z - min.z) len = max.z - min.z;
273  scaling = 1.0 / len; // scale coordinates by length of the cube
274  center = 0.5*(min + max); // center of cube
275  iout << iINFO << "FMM scaling length set to " << len << " A\n" << endi;
276  iout << iINFO << "FMM center set to " << center << "\n" << endi;
277 
278  // allocate buffer space
279  chargebuffer = new double[numAtoms]; // double *
280  epotbuffer = new double[numAtoms]; // double *
281  posbuffer = new double[numAtoms][3]; // double(*)[3]
282  forcebuffer = new double[numAtoms][3]; // double(*)[3]
283  if (chargebuffer == 0 || epotbuffer == 0 ||
284  posbuffer == 0 || forcebuffer == 0) {
285  NAMD_die("can't allocate buffer space for FMM data");
286  }
287 
288  // scale the charges - these won't change
289  double celec = sqrt(COULOMB / simParams->dielectric);
290  for (i = 0; i < numAtoms; i++) {
291  chargebuffer[i] = celec * coord[i].charge;
292  }
293  fmmsolver = 1; // is setup
294  }
295 
296  // translate and scale positions into [-1/2,1/2]^3
297  for (i = 0; i < numAtoms; i++) {
298  posbuffer[i][0] = scaling*(coord[i].position.x - center.x);
299  posbuffer[i][1] = scaling*(coord[i].position.y - center.y);
300  posbuffer[i][2] = scaling*(coord[i].position.z - center.z);
301  }
302 
303  // call the FMM solver
304  int errcode = 0;
305 #ifdef FMM_SOLVER
306  fmmlap_uni_(&numAtoms, posbuffer, chargebuffer, epotbuffer, forcebuffer,
307  &nlevels, &errcode);
308 #else
309  NAMD_die("Must link to FMM library to use FMM\n");
310 #endif
311 
312  // scale force and potentials
313  for (i = 0; i < numAtoms; i++) {
314  double qs = chargebuffer[i] * scaling;
315  force[i].x = qs * scaling * forcebuffer[i][0];
316  force[i].y = qs * scaling * forcebuffer[i][1];
317  force[i].z = qs * scaling * forcebuffer[i][2];
318  energy += qs * epotbuffer[i];
319  }
320  energy *= 0.5;
321 
322  // distribute forces
323  for (int j=0; j < numSources; j++) {
324  FmmSerialCoordMsg *cmsg = coordMsgs[j];
325  coordMsgs[j] = 0;
326  FmmSerialForceMsg *fmsg = new (cmsg->numAtoms, 0) FmmSerialForceMsg;
327 
328  for (int i=0; i < cmsg->numAtoms; i++) {
329  fmsg->force[i] = force[cmsg->coord[i].id];
330  }
331 
332  if ( ! j ) { // set virial and energy only for first message
333  fmsg->energy = energy;
334  for (int k=0; k < 3; k++) {
335  for (int l=0; l < 3; l++) {
336  fmsg->virial[k][l] = virial[k][l];
337  }
338  }
339  }
340  else { // set other messages to zero, add into reduction only once
341  fmsg->energy = 0;
342  for (int k=0; k < 3; k++) {
343  for (int l=0; l < 3; l++) {
344  fmsg->virial[k][l] = 0;
345  }
346  }
347  }
348 
349  fmmProxy[cmsg->sourceNode].recvForce(fmsg);
350  delete cmsg;
351  }
352 }
353 
355  fmmCompute->saveResults(msg);
356  delete oldmsg;
357  oldmsg = msg;
358 }
359 
361 {
363 
364  FmmSerialForce *results_ptr = msg->force;
365 
366  // add in forces
367  for (ap = ap.begin(); ap != ap.end(); ap++) {
368  Results *r = (*ap).forceBox->open();
369  Force *f = r->f[Results::slow];
370  int numAtoms = (*ap).p->getNumAtoms();
371 
372  for(int i=0; i<numAtoms; ++i) {
373  f[i].x += results_ptr->x;
374  f[i].y += results_ptr->y;
375  f[i].z += results_ptr->z;
376  ++results_ptr;
377  }
378 
379  (*ap).forceBox->close(&r);
380  }
381 
382  reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += msg->energy;
383  reduction->item(REDUCTION_VIRIAL_SLOW_XX) += msg->virial[0][0];
384  reduction->item(REDUCTION_VIRIAL_SLOW_XY) += msg->virial[0][1];
385  reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += msg->virial[0][2];
386  reduction->item(REDUCTION_VIRIAL_SLOW_YX) += msg->virial[1][0];
387  reduction->item(REDUCTION_VIRIAL_SLOW_YY) += msg->virial[1][1];
388  reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += msg->virial[1][2];
389  reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += msg->virial[2][0];
390  reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += msg->virial[2][1];
391  reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += msg->virial[2][2];
392  reduction->submit();
393 }
394 
395 #include "ComputeFmmSerialMgr.def.h"
static Node * Object()
Definition: Node.h:86
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
virtual ~ComputeFmmSerial()
FmmSerialForce * force
int ComputeID
Definition: NamdTypes.h:183
static PatchMap * Object()
Definition: PatchMap.h:27
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
ComputeHomePatchList patchList
#define COULOMB
Definition: common.h:44
BigReal & item(int i)
Definition: ReductionMgr.h:312
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
void recvCoord(FmmSerialCoordMsg *)
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define iout
Definition: InfoStream.h:87
BigReal FMMPadding
Charge charge
Definition: NamdTypes.h:54
void setCompute(ComputeFmmSerial *c)
ComputeFmmSerialAtom * coord
ResizeArrayIter< T > end(void) const
Force * f[maxNumForces]
Definition: PatchTypes.h:67
ComputeFmmSerial(ComputeID c)
BigReal virial[3][3]
BigReal x
Definition: Vector.h:66
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
void saveResults(FmmSerialForceMsg *)
#define simParams
Definition: Output.C:127
void recvForce(FmmSerialForceMsg *)
BigReal y
Definition: Vector.h:66
Force FmmSerialForce
k< npairi;++k){TABENERGY(const int numtypes=simParams->tableNumTypes;const float table_spacing=simParams->tableSpacing;const int npertype=(int)(namdnearbyint(simParams->tableMaxDist/simParams->tableSpacing)+1);) int table_i=(r2iilist[2 *k] >> 14)+r2_delta_expc;const int j=pairlisti[k];#define p_j BigReal diffa=r2list[k]-r2_table[table_i];#define table_four_i TABENERGY(register const int tabtype=-1-(lj_pars->A< 0?lj_pars->A:0);) BigReal kqq=kq_i *p_j-> charge
BigReal dielectric
void submit(void)
Definition: ReductionMgr.h:323
infostream & endi(infostream &s)
Definition: InfoStream.C:38
int b_p() const
Definition: Lattice.h:274
gridSize x
int a_p() const
Definition: Lattice.h:273
Molecule * molecule
Definition: Node.h:176
ResizeArrayIter< T > begin(void) const
double BigReal
Definition: common.h:112
int c_p() const
Definition: Lattice.h:275