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