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
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
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;
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
00201 if (self->n_atom < -1) {
00202 Tcl_SetObjResult(interp, Tcl_NewIntObj((long)(0)));
00203 return TCL_OK;
00204 }
00205
00206
00207 do while ( self->n_atom < 0 || ++self->i_atom >= self->n_atom ) {
00208 if ( self->n_atom < 0 ) {
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;
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());
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()) {
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()) {
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()) {
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