NAMD
ComputeExt.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 "ComputeExt.h"
13 #include "ComputeExtMgr.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 #ifndef SQRT_PI
30 #define SQRT_PI 1.7724538509055160273 /* mathematica 15 digits*/
31 #endif
32 
35  float charge;
36  int id;
37 };
38 
39 class ExtCoordMsg : public CMessage_ExtCoordMsg {
40 public:
42  int numAtoms;
45 };
46 
47 class ExtForceMsg : public CMessage_ExtForceMsg {
48 public:
50  BigReal virial[3][3];
52 };
53 
54 class ComputeExtMgr : public CBase_ComputeExtMgr {
55 public:
56  ComputeExtMgr();
58 
59  void setCompute(ComputeExt *c) { extCompute = c; }
60 
61  void recvCoord(ExtCoordMsg *);
62  void recvForce(ExtForceMsg *);
63 
64 private:
65  CProxy_ComputeExtMgr extProxy;
66  ComputeExt *extCompute;
67 
68  int numSources;
69  int numArrived;
70  ExtCoordMsg **coordMsgs;
71  int numAtoms;
72  ComputeExtAtom *coord;
73  ExtForce *force;
74  ExtForceMsg *oldmsg;
75 };
76 
78  extProxy(thisgroup), extCompute(0), numSources(0), numArrived(0),
79  coordMsgs(0), coord(0), force(0), oldmsg(0), numAtoms(0) {
80  CkpvAccess(BOCclass_group).computeExtMgr = thisgroup;
81 }
82 
84  for ( int i=0; i<numSources; ++i ) { delete coordMsgs[i]; }
85  delete [] coordMsgs;
86  delete [] coord;
87  delete [] force;
88  delete oldmsg;
89 }
90 
93 {
94  CProxy_ComputeExtMgr::ckLocalBranch(
95  CkpvAccess(BOCclass_group).computeExtMgr)->setCompute(this);
96 
98 
99 }
100 
102 {
103 }
104 
106 {
108 
109 #if 0
110  // Skip computations if nothing to do.
111  if ( ! patchList[0].p->flags.doFullElectrostatics )
112  {
113  for (ap = ap.begin(); ap != ap.end(); ap++) {
114  CompAtom *x = (*ap).positionBox->open();
115  Results *r = (*ap).forceBox->open();
116  (*ap).positionBox->close(&x);
117  (*ap).forceBox->close(&r);
118  }
119  reduction->submit();
120  return;
121  }
122 #endif
123 
124  // allocate message
125  int numLocalAtoms = 0;
126  for (ap = ap.begin(); ap != ap.end(); ap++) {
127  numLocalAtoms += (*ap).p->getNumAtoms();
128  }
129 
130  ExtCoordMsg *msg = new (numLocalAtoms, 0) ExtCoordMsg;
131  msg->sourceNode = CkMyPe();
132  msg->numAtoms = numLocalAtoms;
133  msg->lattice = patchList[0].p->flags.lattice;
134  ComputeExtAtom *data_ptr = msg->coord;
135 
136  // get positions
137  for (ap = ap.begin(); ap != ap.end(); ap++) {
138  CompAtom *x = (*ap).positionBox->open();
139  CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
140 #if 0
141  if ( patchList[0].p->flags.doMolly ) {
142  (*ap).positionBox->close(&x);
143  x = (*ap).avgPositionBox->open();
144  }
145 #endif
146  int numAtoms = (*ap).p->getNumAtoms();
147 
148  for(int i=0; i<numAtoms; ++i)
149  {
150  data_ptr->position = x[i].position;
151  data_ptr->charge = x[i].charge;
152  data_ptr->id = xExt[i].id;
153  ++data_ptr;
154  }
155 
156 #if 0
157  if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
158  else { (*ap).positionBox->close(&x); }
159 #endif
160  (*ap).positionBox->close(&x);
161  }
162 
163  CProxy_ComputeExtMgr extProxy(CkpvAccess(BOCclass_group).computeExtMgr);
164  extProxy[0].recvCoord(msg);
165 
166 }
167 
169  if ( ! numSources ) {
170  numSources = (PatchMap::Object())->numNodesWithPatches();
171  coordMsgs = new ExtCoordMsg*[numSources];
172  for ( int i=0; i<numSources; ++i ) { coordMsgs[i] = 0; }
173  numArrived = 0;
174  numAtoms = Node::Object()->molecule->numAtoms;
175  coord = new ComputeExtAtom[numAtoms];
176  force = new ExtForce[numAtoms];
177  }
178 
179  int i;
180  for ( i=0; i < msg->numAtoms; ++i ) {
181  coord[msg->coord[i].id] = msg->coord[i];
182  }
183 
184  coordMsgs[numArrived] = msg;
185  ++numArrived;
186 
187  if ( numArrived < numSources ) return;
188  numArrived = 0;
189 
190  // ALL DATA ARRIVED --- CALCULATE FORCES
191  Lattice lattice = msg->lattice;
193  FILE *file;
194  int iret;
195 
196  // write coordinates to file
197  //iout << "writing to file " << simParams->extCoordFilename << "\n" << endi;
198  file = fopen(simParams->extCoordFilename,"w");
199  if ( ! file ) { NAMD_die(strerror(errno)); }
200  for ( i=0; i<numAtoms; ++i ) {
201  int id = coord[i].id + 1;
202  double charge = coord[i].charge;
203  double x = coord[i].position.x;
204  double y = coord[i].position.y;
205  double z = coord[i].position.z;
206  iret = fprintf(file,"%d %f %f %f %f\n",id,charge,x,y,z);
207  if ( iret < 0 ) { NAMD_die(strerror(errno)); }
208  }
209  // write periodic cell lattice (0 0 0 if non-periodic)
210  Vector a = lattice.a(); if ( ! lattice.a_p() ) a = Vector(0,0,0);
211  Vector b = lattice.b(); if ( ! lattice.b_p() ) b = Vector(0,0,0);
212  Vector c = lattice.c(); if ( ! lattice.c_p() ) c = Vector(0,0,0);
213  iret = fprintf(file,"%f %f %f\n%f %f %f\n%f %f %f\n",
214  a.x, a.y, a.z, b.x, b.y, b.z, c.x, c.y, c.z);
215  if ( iret < 0 ) { NAMD_die(strerror(errno)); }
216  fclose(file);
217 
218  // run user-specified command
219  //iout << "running command " << simParams->extForcesCommand << "\n" << endi;
220  iret = system(simParams->extForcesCommand);
221  if ( iret == -1 ) { NAMD_die(strerror(errno)); }
222  if ( iret ) { NAMD_die("Error running command for external forces."); }
223 
224  // remove coordinate file
225  iret = remove(simParams->extCoordFilename);
226  if ( iret ) { NAMD_die(strerror(errno)); }
227 
228  // read forces from file (overwrite positions)
229  //iout << "reading from file " << simParams->extForceFilename << "\n" << endi;
230  file = fopen(simParams->extForceFilename,"r");
231  if ( ! file ) { NAMD_die(strerror(errno)); }
232  for ( i=0; i<numAtoms; ++i ) {
233  int id, replace;
234  double x, y, z;
235  iret = fscanf(file,"%d %d %lf %lf %lf\n", &id, &replace, &x, &y, &z);
236  if ( iret != 5 ) { NAMD_die("Error reading external forces file."); }
237  if ( id != i + 1 ) { NAMD_die("Atom ID error in external forces file."); }
238  force[i].force.x = x; force[i].force.y = y; force[i].force.z = z;
239  force[i].replace = replace;
240  }
241  // read energy and virial if they are present
242  // virial used by NAMD is -'ve of normal convention, so reverse it!
243  // virial[i][j] in file should be sum of -1 * f_i * r_j
244  double energy;
245  double virial[3][3];
246  iret = fscanf(file,"%lf\n", &energy);
247  if ( iret != 1 ) {
248  energy = 0;
249  for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] = 0;
250  } else {
251  iret = fscanf(file,"%lf %lf %lf\n%lf %lf %lf\n%lf %lf %lf\n",
252  &virial[0][0], &virial[0][1], &virial[0][2],
253  &virial[1][0], &virial[1][1], &virial[1][2],
254  &virial[2][0], &virial[2][1], &virial[2][2]);
255  if ( iret != 9 ) {
256  for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] = 0;
257  } else {
258  // virial used by NAMD is -'ve of normal convention, so reverse it!
259  for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l ) virial[k][l] *= -1.0;
260  }
261  }
262  fclose(file);
263 
264  // remove force file
265  iret = remove(simParams->extForceFilename);
266  if ( iret ) { NAMD_die(strerror(errno)); }
267 
268  // distribute forces
269 
270  for ( int j=0; j < numSources; ++j ) {
271  ExtCoordMsg *cmsg = coordMsgs[j];
272  coordMsgs[j] = 0;
273  ExtForceMsg *fmsg = new (cmsg->numAtoms, 0) ExtForceMsg;
274  for ( int i=0; i < cmsg->numAtoms; ++i ) {
275  fmsg->force[i] = force[cmsg->coord[i].id];
276  }
277  if ( ! j ) {
278  fmsg->energy = energy;
279  for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l )
280  fmsg->virial[k][l] = virial[k][l];
281  } else {
282  fmsg->energy = 0;
283  for ( int k=0; k<3; ++k ) for ( int l=0; l<3; ++l )
284  fmsg->virial[k][l] = 0;
285  }
286  extProxy[cmsg->sourceNode].recvForce(fmsg);
287  delete cmsg;
288  }
289 
290 }
291 
293  extCompute->saveResults(msg);
294  delete oldmsg;
295  oldmsg = msg;
296 }
297 
299 {
301 
302  ExtForce *results_ptr = msg->force;
303 
304  // add in forces
305  for (ap = ap.begin(); ap != ap.end(); ap++) {
306  Results *r = (*ap).forceBox->open();
307  Force *f = r->f[Results::normal];
308  int numAtoms = (*ap).p->getNumAtoms();
309 
310  int replace = 0;
311  ExtForce *replacementForces = results_ptr;
312  for(int i=0; i<numAtoms; ++i) {
313  if ( results_ptr->replace ) replace = 1;
314  else f[i] += results_ptr->force;
315  ++results_ptr;
316  }
317  if ( replace ) (*ap).p->replaceForces(replacementForces);
318 
319  (*ap).forceBox->close(&r);
320  }
321 
322  reduction->item(REDUCTION_MISC_ENERGY) += msg->energy;
323  reduction->item(REDUCTION_VIRIAL_NORMAL_XX) += msg->virial[0][0];
324  reduction->item(REDUCTION_VIRIAL_NORMAL_XY) += msg->virial[0][1];
325  reduction->item(REDUCTION_VIRIAL_NORMAL_XZ) += msg->virial[0][2];
326  reduction->item(REDUCTION_VIRIAL_NORMAL_YX) += msg->virial[1][0];
327  reduction->item(REDUCTION_VIRIAL_NORMAL_YY) += msg->virial[1][1];
328  reduction->item(REDUCTION_VIRIAL_NORMAL_YZ) += msg->virial[1][2];
329  reduction->item(REDUCTION_VIRIAL_NORMAL_ZX) += msg->virial[2][0];
330  reduction->item(REDUCTION_VIRIAL_NORMAL_ZY) += msg->virial[2][1];
331  reduction->item(REDUCTION_VIRIAL_NORMAL_ZZ) += msg->virial[2][2];
332  reduction->submit();
333 }
334 
335 #include "ComputeExtMgr.def.h"
336 
static Node * Object()
Definition: Node.h:86
void doWork()
Definition: ComputeExt.C:105
Lattice lattice
Definition: ComputeExt.C:43
int ComputeID
Definition: NamdTypes.h:183
char extForcesCommand[256]
void recvCoord(ExtCoordMsg *)
Definition: ComputeExt.C:168
static PatchMap * Object()
Definition: PatchMap.h:27
virtual ~ComputeExt()
Definition: ComputeExt.C:101
void recvForce(ExtForceMsg *)
Definition: ComputeExt.C:292
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
ComputeHomePatchList patchList
BigReal & item(int i)
Definition: ReductionMgr.h:312
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
int replace
Definition: NamdTypes.h:201
BigReal virial[3][3]
Definition: ComputeExt.C:50
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
Charge charge
Definition: NamdTypes.h:54
ComputeExtAtom * coord
Definition: ComputeExt.C:44
ResizeArrayIter< T > end(void) const
void setCompute(ComputeExt *c)
Definition: ComputeExt.C:59
ComputeExt(ComputeID c)
Definition: ComputeExt.C:91
BigReal energy
Definition: ComputeExt.C:49
gridSize z
Force * f[maxNumForces]
Definition: PatchTypes.h:67
BigReal x
Definition: Vector.h:66
void saveResults(ExtForceMsg *)
Definition: ComputeExt.C:298
int numAtoms
Definition: Molecule.h:556
void NAMD_die(const char *err_msg)
Definition: common.C:83
Force force
Definition: NamdTypes.h:202
int sourceNode
Definition: ComputeExt.C:41
char extCoordFilename[128]
ExtForce * force
Definition: ComputeExt.C:51
#define simParams
Definition: Output.C:127
int numAtoms
Definition: ComputeExt.C:42
BigReal y
Definition: Vector.h:66
Vector b() const
Definition: Lattice.h:253
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
gridSize y
void submit(void)
Definition: ReductionMgr.h:323
int b_p() const
Definition: Lattice.h:274
gridSize x
int a_p() const
Definition: Lattice.h:273
Molecule * molecule
Definition: Node.h:176
Vector a() const
Definition: Lattice.h:252
char extForceFilename[128]
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
Position position
Definition: ComputeExt.C:34