NAMD: ComputeTclBC.C Source File

00001 
00007 #include "InfoStream.h"
00008 #include "ComputeTclBC.h"
00009 #include "Node.h"
00010 #include "SimParameters.h"
00011 #include "Patch.h"
00012 #include "Molecule.h"
00013 
00014 #ifdef NAMD_TCL
00015 #include <tcl.h>
00016 #endif
00017 #include "TclCommands.h"
00018 
00019 #define WRAPMODE_PATCH 0
00020 #define WRAPMODE_INPUT 1
00021 #define WRAPMODE_CELL 2
00022 #define WRAPMODE_NEAREST 3
00023 
00024 ComputeTclBC::ComputeTclBC(ComputeID c)
00025   : ComputeHomePatches(c), ap(patchList) {
00026 
00027   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00028   SimParameters *simParams = Node::Object()->simParameters;
00029 
00030   wrapmode = WRAPMODE_PATCH;
00031 
00032   drops.resize(Node::Object()->molecule->numAtoms);
00033   cleardrops();
00034 
00035 #ifdef NAMD_TCL
00036   if ( CkMyRank() && ! simParams->tclIsThreaded ) {
00037     NAMD_die("Sorry, tclBC requires TCL to be built with --enable-threads to use multiple threads per process.");
00038   }
00039 
00040   interp = Tcl_CreateInterp();
00041   tcl_vector_math_init(interp);
00042 
00043   Tcl_CreateCommand(interp, "print", Tcl_print,
00044     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00045   Tcl_CreateCommand(interp, "wrapmode", Tcl_wrapmode,
00046     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00047   Tcl_CreateObjCommand(interp, "cleardrops", Tcl_cleardrops,
00048     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00049 
00050   // run script to define calcforces, etc.
00051   if ( simParams->tclBCScript ) {
00052     int code = Tcl_Eval(interp,simParams->tclBCScript);
00053     const char *result = Tcl_GetStringResult(interp);
00054     if (result && *result != 0) CkPrintf("TCL: %s\n",result);
00055     if (code != TCL_OK) {
00056       const char *errorInfo = Tcl_GetVar(interp,"errorInfo",0);
00057       NAMD_die(errorInfo ? errorInfo : "Unknown Tcl error");
00058     }
00059   } else NAMD_bug("tclBCScript pointer was NULL");
00060 
00061   // don't want these available until calcforces call
00062   Tcl_CreateObjCommand(interp, "dropatom", Tcl_dropatom,
00063     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00064   Tcl_CreateObjCommand(interp, "nextatom", Tcl_nextatom,
00065     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00066   Tcl_CreateObjCommand(interp, "getcoord", Tcl_getcoord,
00067     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00068   Tcl_CreateObjCommand(interp, "getcell", Tcl_getcell,
00069     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00070   Tcl_CreateObjCommand(interp, "getmass", Tcl_getmass,
00071     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00072   Tcl_CreateObjCommand(interp, "getcharge", Tcl_getcharge,
00073     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00074   Tcl_CreateObjCommand(interp, "getid", Tcl_getid,
00075     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00076   Tcl_CreateObjCommand(interp, "addforce", Tcl_addforce,
00077     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00078   Tcl_CreateObjCommand(interp, "addenergy", Tcl_addenergy,
00079     (ClientData) this, (Tcl_CmdDeleteProc *) NULL);
00080 
00081 #else
00082 
00083   NAMD_die("Sorry, tclBC is not available; built without TCL.");
00084 
00085 #endif
00086 
00087 }
00088 
00089 
00090 ComputeTclBC::~ComputeTclBC() {
00091 #ifdef NAMD_TCL
00092   Tcl_DeleteInterp(interp);
00093 #endif
00094   delete reduction;
00095 }
00096 
00097 
00098 void ComputeTclBC::doWork() {
00099 
00100   SimParameters *simParams = Node::Object()->simParameters;
00101   lattice = &(patchList[0].p->lattice);
00102   const int step = patchList[0].p->flags.step;
00103   char cmd[128];
00104 
00105   energy = 0;
00106   n_atom = -1;  // set initial flags for iteration by nextatom
00107 
00108 #ifdef NAMD_TCL
00109   sprintf(cmd,"calcforces %d %d %s",step,hasPatchZero,simParams->tclBCArgs);
00110   int code = Tcl_Eval(interp,cmd);
00111   if (code != TCL_OK) {
00112     const char *errorInfo = Tcl_GetVar(interp,"errorInfo",0);
00113     NAMD_die(errorInfo ? errorInfo : "Unknown Tcl error");
00114   }
00115   if (n_atom != -2) {
00116     NAMD_die("tclBCScript failed to call nextatom until failure");
00117   }
00118 #endif
00119 
00120   reduction->item(REDUCTION_BC_ENERGY) += energy;
00121   reduction->submit();
00122 
00123 }
00124 
00125 #ifdef NAMD_TCL
00126 
00127 int ComputeTclBC::Tcl_print(ClientData,
00128         Tcl_Interp *, int argc, const char *argv[]) {
00129   Tcl_DString msg;
00130   Tcl_DStringInit(&msg);
00131   for ( int i = 1; i < argc; ++i ) {
00132     Tcl_DStringAppend(&msg," ",-1);
00133     Tcl_DStringAppend(&msg,argv[i],-1);
00134   }
00135   CkPrintf("TCL:%s\n",Tcl_DStringValue(&msg));
00136   Tcl_DStringFree(&msg);
00137   return TCL_OK;
00138 }
00139 
00140 int ComputeTclBC::Tcl_wrapmode(ClientData clientData,
00141         Tcl_Interp *interp, int argc, const char *argv[]) {
00142   if (argc != 2) {
00143     Tcl_SetResult(interp,(char*)"usage: wrapmode patch|input|cell|nearest",
00144                                                                 TCL_VOLATILE);
00145     return TCL_ERROR;
00146   }
00147   ComputeTclBC *self = (ComputeTclBC *)clientData;
00148 
00149   if ( ! strcmp(argv[1],"patch") ) self->wrapmode = WRAPMODE_PATCH;
00150   else if ( ! strcmp(argv[1],"input") ) self->wrapmode = WRAPMODE_INPUT;
00151   else if ( ! strcmp(argv[1],"cell") ) self->wrapmode = WRAPMODE_CELL;
00152   else if ( ! strcmp(argv[1],"nearest") ) self->wrapmode = WRAPMODE_NEAREST;
00153   else {
00154     Tcl_SetResult(interp,(char*)"usage: wrapmode patch|input|cell|nearest",
00155                                                                 TCL_VOLATILE);
00156     return TCL_ERROR;
00157   }
00158 
00159   return TCL_OK;
00160 }
00161 
00162 int ComputeTclBC::Tcl_cleardrops(ClientData clientData,
00163         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00164   if (objc != 1) {
00165     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00166     return TCL_ERROR;
00167   }
00168   ComputeTclBC *self = (ComputeTclBC *)clientData;
00169 
00170   self->cleardrops();
00171 
00172   return TCL_OK;
00173 }
00174 
00175 int ComputeTclBC::Tcl_dropatom(ClientData clientData,
00176         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00177   if (objc != 1) {
00178     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00179     return TCL_ERROR;
00180   }
00181   ComputeTclBC *self = (ComputeTclBC *)clientData;
00182   if ( self->n_atom <= 0 ) {
00183     Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
00184     return TCL_ERROR;
00185   }
00186 
00187   self->drops[self->fullatoms[self->i_atom].id] = 1;
00188 
00189   return TCL_OK;
00190 }
00191 
00192 int ComputeTclBC::Tcl_nextatom(ClientData clientData,
00193         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00194   if (objc != 1) {
00195     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00196     return TCL_ERROR;
00197   }
00198   ComputeTclBC *self = (ComputeTclBC *)clientData;
00199 
00200   // n_atom = -2 after all atoms processed
00201   if (self->n_atom < -1) {
00202     Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(0)));
00203     return TCL_OK;
00204   }
00205 
00206   // assume n_atom = -1 before first call
00207   do while ( self->n_atom < 0 || ++self->i_atom >= self->n_atom ) {
00208     if ( self->n_atom < 0 ) {  // first call
00209       self->ap = self->ap.begin();
00210     } else {
00211       (*(self->ap)).positionBox->close(&(self->atoms));
00212       (*(self->ap)).forceBox->close(&((*(self->ap)).r));
00213       self->ap++;
00214     }
00215     if ( self->ap == self->ap.end() ) {
00216       self->n_atom = -2;  // set error condition
00217       Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(0)));
00218       return TCL_OK;
00219     }
00220     self->i_atom = -1;
00221     self->n_atom = (*(self->ap)).p->getNumAtoms();
00222     self->fullatoms = (*(self->ap)).p->getAtomList().begin();
00223     self->atoms = (*(self->ap)).positionBox->open();
00224     (*(self->ap)).r = (*(self->ap)).forceBox->open();
00225     self->forces = (*(self->ap)).r->f[Results::normal];
00226   } while ( self->drops[self->fullatoms[self->i_atom].id] );
00227 
00228   Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(1)));
00229   return TCL_OK;
00230 }
00231 
00232 int ComputeTclBC::Tcl_getcoord(ClientData clientData,
00233         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00234   if (objc != 1) {
00235     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00236     return TCL_ERROR;
00237   }
00238   ComputeTclBC *self = (ComputeTclBC *)clientData;
00239   if ( self->n_atom <= 0 ) {
00240     Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
00241     return TCL_ERROR;
00242   }
00243 
00244   int i = self->i_atom;
00245   Position pos = self->atoms[i].position;
00246   switch ( self->wrapmode ) {
00247   case WRAPMODE_PATCH:
00248     break;
00249   case WRAPMODE_INPUT:
00250     pos = self->lattice->reverse_transform(pos,self->fullatoms[i].transform);
00251     break;
00252   case WRAPMODE_CELL:
00253     pos += self->lattice->wrap_delta(pos);
00254     break;
00255   case WRAPMODE_NEAREST:
00256     pos += self->lattice->wrap_nearest_delta(pos);
00257     break;
00258   }
00259 
00260   Tcl_Obj *newlist = Tcl_NewListObj(0, NULL);
00261   Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(pos.x)));
00262   Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(pos.y)));
00263   Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(pos.z)));
00264   Tcl_SetObjResult(interp, newlist);
00265   return TCL_OK;
00266 }
00267 
00268 int ComputeTclBC::Tcl_getcell(ClientData clientData,
00269         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00270   if (objc != 1) {
00271     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00272     return TCL_ERROR;
00273   }
00274   ComputeTclBC *self = (ComputeTclBC *)clientData;
00275 
00276   Tcl_Obj *newcell = Tcl_NewListObj(0, NULL);
00277   Tcl_Obj *newlist = Tcl_NewListObj(0, NULL);
00278   Vector o(self->lattice->origin());  // always have origin
00279   Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(o.x)));
00280   Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(o.y)));
00281   Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(o.z)));
00282   Tcl_ListObjAppendElement(interp, newcell, newlist);
00283   if (self->lattice->a_p()) {  // only if periodic
00284     newlist = Tcl_NewListObj(0, NULL);
00285     Vector a(self->lattice->a());
00286     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(a.x)));
00287     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(a.y)));
00288     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(a.z)));
00289     Tcl_ListObjAppendElement(interp, newcell, newlist);
00290   }
00291   if (self->lattice->b_p()) {  // only if periodic
00292     newlist = Tcl_NewListObj(0, NULL);
00293     Vector b(self->lattice->b());
00294     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(b.x)));
00295     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(b.y)));
00296     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(b.z)));
00297     Tcl_ListObjAppendElement(interp, newcell, newlist);
00298   }
00299   if (self->lattice->c_p()) {  // only if periodic
00300     newlist = Tcl_NewListObj(0, NULL);
00301     Vector c(self->lattice->c());
00302     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(c.x)));
00303     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(c.y)));
00304     Tcl_ListObjAppendElement(interp, newlist, Tcl_NewDoubleObj((double)(c.z)));
00305     Tcl_ListObjAppendElement(interp, newcell, newlist);
00306   }
00307   Tcl_SetObjResult(interp, newcell);
00308   return TCL_OK;
00309 }
00310 
00311 int ComputeTclBC::Tcl_getmass(ClientData clientData,
00312         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00313   if (objc != 1) {
00314     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00315     return TCL_ERROR;
00316   }
00317   ComputeTclBC *self = (ComputeTclBC *)clientData;
00318   if ( self->n_atom <= 0 ) {
00319     Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
00320     return TCL_ERROR;
00321   }
00322 
00323   int i = self->i_atom;
00324   Tcl_SetObjResult(interp, Tcl_NewDoubleObj((double)(self->fullatoms[i].mass)));
00325   return TCL_OK;
00326 }
00327 
00328 int ComputeTclBC::Tcl_getcharge(ClientData clientData,
00329         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00330   if (objc != 1) {
00331     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00332     return TCL_ERROR;
00333   }
00334   ComputeTclBC *self = (ComputeTclBC *)clientData;
00335   if ( self->n_atom <= 0 ) {
00336     Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
00337     return TCL_ERROR;
00338   }
00339 
00340   int i = self->i_atom;
00341   Tcl_SetObjResult(interp, Tcl_NewDoubleObj((double)(self->atoms[i].charge)));
00342   return TCL_OK;
00343 }
00344 
00345 int ComputeTclBC::Tcl_getid(ClientData clientData,
00346         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00347   if (objc != 1) {
00348     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00349     return TCL_ERROR;
00350   }
00351   ComputeTclBC *self = (ComputeTclBC *)clientData;
00352   if ( self->n_atom <= 0 ) {
00353     Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
00354     return TCL_ERROR;
00355   }
00356 
00357   int i = self->i_atom;
00358   Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(self->fullatoms[i].id + 1)));
00359   return TCL_OK;
00360 }
00361 
00362 int ComputeTclBC::Tcl_addforce(ClientData clientData,
00363         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00364   if (objc != 2) {
00365     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00366     return TCL_ERROR;
00367   }
00368 
00369   Tcl_Obj **force;  int fnum;  double x,y,z;
00370   if (Tcl_ListObjGetElements(interp, objv[1], &fnum, &force) != TCL_OK) {
00371     return TCL_ERROR;
00372   }
00373   if ( (fnum != 3) ||
00374        (Tcl_GetDoubleFromObj(interp, force[0],&x) != TCL_OK) ||
00375        (Tcl_GetDoubleFromObj(interp, force[1],&y) != TCL_OK) ||
00376        (Tcl_GetDoubleFromObj(interp, force[2],&z) != TCL_OK) ) {
00377     Tcl_SetResult(interp,(char*)"force not a vector",TCL_VOLATILE);
00378     return TCL_ERROR;
00379   }
00380 
00381   ComputeTclBC *self = (ComputeTclBC *)clientData;
00382   if ( self->n_atom <= 0 ) {
00383     Tcl_SetResult(interp,(char*)"no atom available",TCL_VOLATILE);
00384     return TCL_ERROR;
00385   }
00386   int i = self->i_atom;
00387   self->forces[i].x += x;
00388   self->forces[i].y += y;
00389   self->forces[i].z += z;
00390 
00391   return TCL_OK;
00392 }
00393 
00394 int ComputeTclBC::Tcl_addenergy(ClientData clientData,
00395         Tcl_Interp *interp, int objc, Tcl_Obj * const objv[]) {
00396   if (objc != 2) {
00397     Tcl_SetResult(interp,(char*)"wrong # args",TCL_VOLATILE);
00398     return TCL_ERROR;
00399   }
00400 
00401   double energy;
00402   if ( Tcl_GetDoubleFromObj(interp, objv[1], &energy) != TCL_OK ) {
00403     Tcl_SetResult(interp,(char*)"energy not a number",TCL_VOLATILE);
00404     return TCL_ERROR;
00405   }
00406 
00407   ComputeTclBC *self = (ComputeTclBC *)clientData;
00408   self->energy += energy;
00409 
00410   return TCL_OK;
00411 }
00412 
00413 #endif
00414