ComputeTclBC.C

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

Generated on Tue Sep 26 01:17:12 2017 for NAMD by  doxygen 1.4.7