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