NAMD
ComputeTclBC.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "ComputeTclBC.h"
9 #include "Node.h"
10 #include "SimParameters.h"
11 #include "Patch.h"
12 #include "Molecule.h"
13 
14 #ifdef NAMD_TCL
15 #include <tcl.h>
16 #endif
17 #include "TclCommands.h"
18 
19 #define WRAPMODE_PATCH 0
20 #define WRAPMODE_INPUT 1
21 #define WRAPMODE_CELL 2
22 #define WRAPMODE_NEAREST 3
23 
25  : ComputeHomePatches(c), ap(patchList) {
26 
29 
30  wrapmode = WRAPMODE_PATCH;
31 
32  drops.resize(Node::Object()->molecule->numAtoms);
33  cleardrops();
34 
35 #ifdef NAMD_TCL
36  if ( CkMyRank() && ! simParams->tclIsThreaded ) {
37  NAMD_die("Sorry, tclBC requires TCL to be built with --enable-threads to use multiple threads per process.");
38  }
39 
40  interp = Tcl_CreateInterp();
41  tcl_vector_math_init(interp);
42 
43  Tcl_CreateCommand(interp, "print", Tcl_print,
44  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
45  Tcl_CreateCommand(interp, "wrapmode", Tcl_wrapmode,
46  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
47  Tcl_CreateObjCommand(interp, "cleardrops", Tcl_cleardrops,
48  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
49 
50  // run script to define calcforces, etc.
51  if ( simParams->tclBCScript ) {
52  int code = Tcl_Eval(interp,simParams->tclBCScript);
53  const char *result = Tcl_GetStringResult(interp);
54  if (result && *result != 0) CkPrintf("TCL: %s\n",result);
55  if (code != TCL_OK) {
56  const char *errorInfo = Tcl_GetVar(interp,"errorInfo",0);
57  NAMD_die(errorInfo ? errorInfo : "Unknown Tcl error");
58  }
59  } else NAMD_bug("tclBCScript pointer was NULL");
60 
61  // don't want these available until calcforces call
62  Tcl_CreateObjCommand(interp, "dropatom", Tcl_dropatom,
63  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
64  Tcl_CreateObjCommand(interp, "nextatom", Tcl_nextatom,
65  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
66  Tcl_CreateObjCommand(interp, "getcoord", Tcl_getcoord,
67  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
68  Tcl_CreateObjCommand(interp, "getcell", Tcl_getcell,
69  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
70  Tcl_CreateObjCommand(interp, "getmass", Tcl_getmass,
71  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
72  Tcl_CreateObjCommand(interp, "getcharge", Tcl_getcharge,
73  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
74  Tcl_CreateObjCommand(interp, "getid", Tcl_getid,
75  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
76  Tcl_CreateObjCommand(interp, "addforce", Tcl_addforce,
77  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
78  Tcl_CreateObjCommand(interp, "addenergy", Tcl_addenergy,
79  (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
80 
81 #else
82 
83  NAMD_die("Sorry, tclBC is not available; built without TCL.");
84 
85 #endif
86 
87 }
88 
89 
91 #ifdef NAMD_TCL
92  Tcl_DeleteInterp(interp);
93 #endif
94  delete reduction;
95 }
96 
97 
99 
101  lattice = &(patchList[0].p->lattice);
102  const int step = patchList[0].p->flags.step;
103  char cmd[128];
104 
105  energy = 0;
106  n_atom = -1; // set initial flags for iteration by nextatom
107 
108 #ifdef NAMD_TCL
109  sprintf(cmd,"calcforces %d %d %s",step,hasPatchZero,simParams->tclBCArgs);
110  int code = Tcl_Eval(interp,cmd);
111  if (code != TCL_OK) {
112  const char *errorInfo = Tcl_GetVar(interp,"errorInfo",0);
113  NAMD_die(errorInfo ? errorInfo : "Unknown Tcl error");
114  }
115  if (n_atom != -2) {
116  NAMD_die("tclBCScript failed to call nextatom until failure");
117  }
118 #endif
119 
120  reduction->item(REDUCTION_BC_ENERGY) += energy;
121  reduction->submit();
122 
123 }
124 
125 #ifdef NAMD_TCL
126 
127 int ComputeTclBC::Tcl_print(ClientData,
128  Tcl_Interp *, int argc, const char *argv[]) {
129  Tcl_DString msg;
130  Tcl_DStringInit(&msg);
131  for ( int i = 1; i < argc; ++i ) {
132  Tcl_DStringAppend(&msg," ",-1);
133  Tcl_DStringAppend(&msg,argv[i],-1);
134  }
135  CkPrintf("TCL:%s\n",Tcl_DStringValue(&msg));
136  Tcl_DStringFree(&msg);
137  return TCL_OK;
138 }
139 
140 int ComputeTclBC::Tcl_wrapmode(ClientData clientData,
141  Tcl_Interp *interp, int argc, const char *argv[]) {
142  if (argc != 2) {
143  Tcl_SetResult(interp,(char*)"usage: wrapmode patch|input|cell|nearest",
144  TCL_VOLATILE);
145  return TCL_ERROR;
146  }
147  ComputeTclBC *self = (ComputeTclBC *)clientData;
148 
149  if ( ! strcmp(argv[1],"patch") ) self->wrapmode = WRAPMODE_PATCH;
150  else if ( ! strcmp(argv[1],"input") ) self->wrapmode = WRAPMODE_INPUT;
151  else if ( ! strcmp(argv[1],"cell") ) self->wrapmode = WRAPMODE_CELL;
152  else if ( ! strcmp(argv[1],"nearest") ) self->wrapmode = WRAPMODE_NEAREST;
153  else {
154  Tcl_SetResult(interp,(char*)"usage: wrapmode patch|input|cell|nearest",
155  TCL_VOLATILE);
156  return TCL_ERROR;
157  }
158 
159  return TCL_OK;
160 }
161 
162 int ComputeTclBC::Tcl_cleardrops(ClientData clientData,
163  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
164  if (objc != 1) {
165  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
166  return TCL_ERROR;
167  }
168  ComputeTclBC *self = (ComputeTclBC *)clientData;
169 
170  self->cleardrops();
171 
172  return TCL_OK;
173 }
174 
175 int ComputeTclBC::Tcl_dropatom(ClientData clientData,
176  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
177  if (objc != 1) {
178  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
179  return TCL_ERROR;
180  }
181  ComputeTclBC *self = (ComputeTclBC *)clientData;
182  if ( self->n_atom <= 0 ) {
183  Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
184  return TCL_ERROR;
185  }
186 
187  self->drops[self->fullatoms[self->i_atom].id] = 1;
188 
189  return TCL_OK;
190 }
191 
192 int ComputeTclBC::Tcl_nextatom(ClientData clientData,
193  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
194  if (objc != 1) {
195  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
196  return TCL_ERROR;
197  }
198  ComputeTclBC *self = (ComputeTclBC *)clientData;
199 
200  // n_atom = -2 after all atoms processed
201  if (self->n_atom < -1) {
202  Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(0)));
203  return TCL_OK;
204  }
205 
206  // assume n_atom = -1 before first call
207  do while ( self->n_atom < 0 || ++self->i_atom >= self->n_atom ) {
208  if ( self->n_atom < 0 ) { // first call
209  self->ap = self->ap.begin();
210  } else {
211  (*(self->ap)).positionBox->close(&(self->atoms));
212  (*(self->ap)).forceBox->close(&((*(self->ap)).r));
213  self->ap++;
214  }
215  if ( self->ap == self->ap.end() ) {
216  self->n_atom = -2; // set error condition
217  Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(0)));
218  return TCL_OK;
219  }
220  self->i_atom = -1;
221  self->n_atom = (*(self->ap)).p->getNumAtoms();
222  self->fullatoms = (*(self->ap)).p->getAtomList().begin();
223  self->atoms = (*(self->ap)).positionBox->open();
224  (*(self->ap)).r = (*(self->ap)).forceBox->open();
225  self->forces = (*(self->ap)).r->f[Results::normal];
226  } while ( self->drops[self->fullatoms[self->i_atom].id] );
227 
228  Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(1)));
229  return TCL_OK;
230 }
231 
232 int ComputeTclBC::Tcl_getcoord(ClientData clientData,
233  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
234  if (objc != 1) {
235  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
236  return TCL_ERROR;
237  }
238  ComputeTclBC *self = (ComputeTclBC *)clientData;
239  if ( self->n_atom <= 0 ) {
240  Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
241  return TCL_ERROR;
242  }
243 
244  int i = self->i_atom;
245  Position pos = self->atoms[i].position;
246  switch ( self->wrapmode ) {
247  case WRAPMODE_PATCH:
248  break;
249  case WRAPMODE_INPUT:
250  pos = self->lattice->reverse_transform(pos,self->fullatoms[i].transform);
251  break;
252  case WRAPMODE_CELL:
253  pos += self->lattice->wrap_delta(pos);
254  break;
255  case WRAPMODE_NEAREST:
256  pos += self->lattice->wrap_nearest_delta(pos);
257  break;
258  }
259 
260  Tcl_Obj *newlist = Tcl_NewListObj(0, NULL);
261  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(pos.x)));
262  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(pos.y)));
263  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(pos.z)));
264  Tcl_SetObjResult(interp, newlist);
265  return TCL_OK;
266 }
267 
268 int ComputeTclBC::Tcl_getcell(ClientData clientData,
269  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
270  if (objc != 1) {
271  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
272  return TCL_ERROR;
273  }
274  ComputeTclBC *self = (ComputeTclBC *)clientData;
275 
276  Tcl_Obj *newcell = Tcl_NewListObj(0, NULL);
277  Tcl_Obj *newlist = Tcl_NewListObj(0, NULL);
278  Vector o(self->lattice->origin()); // always have origin
279  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(o.x)));
280  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(o.y)));
281  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(o.z)));
282  Tcl_ListObjAppendElement(interp, newcell, newlist);
283  if (self->lattice->a_p()) { // only if periodic
284  newlist = Tcl_NewListObj(0, NULL);
285  Vector a(self->lattice->a());
286  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(a.x)));
287  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(a.y)));
288  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(a.z)));
289  Tcl_ListObjAppendElement(interp, newcell, newlist);
290  }
291  if (self->lattice->b_p()) { // only if periodic
292  newlist = Tcl_NewListObj(0, NULL);
293  Vector b(self->lattice->b());
294  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(b.x)));
295  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(b.y)));
296  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(b.z)));
297  Tcl_ListObjAppendElement(interp, newcell, newlist);
298  }
299  if (self->lattice->c_p()) { // only if periodic
300  newlist = Tcl_NewListObj(0, NULL);
301  Vector c(self->lattice->c());
302  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(c.x)));
303  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(c.y)));
304  Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(c.z)));
305  Tcl_ListObjAppendElement(interp, newcell, newlist);
306  }
307  Tcl_SetObjResult(interp, newcell);
308  return TCL_OK;
309 }
310 
311 int ComputeTclBC::Tcl_getmass(ClientData clientData,
312  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
313  if (objc != 1) {
314  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
315  return TCL_ERROR;
316  }
317  ComputeTclBC *self = (ComputeTclBC *)clientData;
318  if ( self->n_atom <= 0 ) {
319  Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
320  return TCL_ERROR;
321  }
322 
323  int i = self->i_atom;
324  Tcl_SetObjResult(interp, Tcl_NewDoubleObj((double)(self->fullatoms[i].mass)));
325  return TCL_OK;
326 }
327 
328 int ComputeTclBC::Tcl_getcharge(ClientData clientData,
329  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
330  if (objc != 1) {
331  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
332  return TCL_ERROR;
333  }
334  ComputeTclBC *self = (ComputeTclBC *)clientData;
335  if ( self->n_atom <= 0 ) {
336  Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
337  return TCL_ERROR;
338  }
339 
340  int i = self->i_atom;
341  Tcl_SetObjResult(interp, Tcl_NewDoubleObj((double)(self->atoms[i].charge)));
342  return TCL_OK;
343 }
344 
345 int ComputeTclBC::Tcl_getid(ClientData clientData,
346  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
347  if (objc != 1) {
348  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
349  return TCL_ERROR;
350  }
351  ComputeTclBC *self = (ComputeTclBC *)clientData;
352  if ( self->n_atom <= 0 ) {
353  Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
354  return TCL_ERROR;
355  }
356 
357  int i = self->i_atom;
358  Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(self->fullatoms[i].id + 1)));
359  return TCL_OK;
360 }
361 
362 int ComputeTclBC::Tcl_addforce(ClientData clientData,
363  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
364  if (objc != 2) {
365  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
366  return TCL_ERROR;
367  }
368 
369  Tcl_Obj **force; int fnum; double x,y,z;
370  if (Tcl_ListObjGetElements(interp, objv[1], &fnum, &force) != TCL_OK) {
371  return TCL_ERROR;
372  }
373  if ( (fnum != 3) ||
374  (Tcl_GetDoubleFromObj(interp, force[0],&x) != TCL_OK) ||
375  (Tcl_GetDoubleFromObj(interp, force[1],&y) != TCL_OK) ||
376  (Tcl_GetDoubleFromObj(interp, force[2],&z) != TCL_OK) ) {
377  Tcl_SetResult(interp,(char*)"force not a vector",TCL_VOLATILE);
378  return TCL_ERROR;
379  }
380 
381  ComputeTclBC *self = (ComputeTclBC *)clientData;
382  if ( self->n_atom <= 0 ) {
383  Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
384  return TCL_ERROR;
385  }
386  int i = self->i_atom;
387  self->forces[i].x += x;
388  self->forces[i].y += y;
389  self->forces[i].z += z;
390 
391  return TCL_OK;
392 }
393 
394 int ComputeTclBC::Tcl_addenergy(ClientData clientData,
395  Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
396  if (objc != 2) {
397  Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
398  return TCL_ERROR;
399  }
400 
401  double energy;
402  if ( Tcl_GetDoubleFromObj(interp, objv[1], &energy) != TCL_OK ) {
403  Tcl_SetResult(interp,(char*)"energy not a number",TCL_VOLATILE);
404  return TCL_ERROR;
405  }
406 
407  ComputeTclBC *self = (ComputeTclBC *)clientData;
408  self->energy += energy;
409 
410  return TCL_OK;
411 }
412 
413 #endif
414 
static Node * Object()
Definition: Node.h:86
int ComputeID
Definition: NamdTypes.h:183
#define WRAPMODE_NEAREST
Definition: ComputeTclBC.C:22
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
ComputeHomePatchList patchList
BigReal & item(int i)
Definition: ReductionMgr.h:312
char * tclBCScript
BigReal z
Definition: Vector.h:66
#define WRAPMODE_PATCH
Definition: ComputeTclBC.C:19
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
if(ComputeNonbondedUtil::goMethod==2)
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
#define WRAPMODE_INPUT
Definition: ComputeTclBC.C:20
virtual ~ComputeTclBC()
Definition: ComputeTclBC.C:90
char tclBCArgs[NAMD_FILENAME_BUFFER_SIZE]
void NAMD_bug(const char *err_msg)
Definition: common.C:123
gridSize z
BigReal x
Definition: Vector.h:66
void doWork()
Definition: ComputeTclBC.C:98
void NAMD_die(const char *err_msg)
Definition: common.C:83
#define simParams
Definition: Output.C:127
int tcl_vector_math_init(Tcl_Interp *interp)
Definition: TclCommands.C:299
void resize(int i)
Definition: ResizeArray.h:84
BigReal y
Definition: Vector.h:66
gridSize y
void submit(void)
Definition: ReductionMgr.h:323
gridSize x
#define WRAPMODE_CELL
Definition: ComputeTclBC.C:21
ComputeTclBC(ComputeID c)
Definition: ComputeTclBC.C:24