TclVolMap.C

00001 /***************************************************************************
00002 *cr 
00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 
00004 *cr University of Illinois 
00005 *cr All Rights Reserved 
00006 *cr 
00007 ***************************************************************************/
00008 
00009 /***************************************************************************
00010 * RCS INFORMATION:
00011 *
00012 * $RCSfile: TclVolMap.C,v $
00013 * $Author: johns $ $Locker: $ $State: Exp $
00014 * $Revision: 1.123 $ $Date: 2019年11月06日 23:36:09 $
00015 *
00016 ***************************************************************************
00017 * DESCRIPTION:
00018 * These are essentially just Tcl wrappers for the volmap commands in
00019 * VolMapCreate.C.
00020 *
00021 ***************************************************************************/
00022 
00023 #include <math.h>
00024 #include <stdlib.h>
00025 #include <tcl.h>
00026 #include "TclCommands.h"
00027 #include "AtomSel.h"
00028 #include "VMDApp.h"
00029 #include "MoleculeList.h"
00030 #include "config.h"
00031 #include "Measure.h"
00032 #include "VolumetricData.h"
00033 #include "VolMapCreate.h"
00034 #include "Inform.h"
00035 #include "MeasureSymmetry.h"
00036 
00037 // XXX Todo List: 
00038 // - Make this into an independent Tcl function 
00039 // (e.g.: "volmap new density <blah...>"
00040 // - Make friends with as-of-yet non-existing VMD volumetric map I/O
00041 // - parse and allow user to set useful params
00042 
00043 // Function: volmap <maptype> <selection> [-weight <none|atom value|string array>] 
00044 // Options for all maptypes:
00045 // -allframes
00046 // -res <res>
00047 // Options for "density":
00048 // -weight <none|atom value|string array>
00049 // -scale <radius scale>
00050 
00051 
00052 static int parse_minmax_args(Tcl_Interp *interp, int arg_minmax,
00053 Tcl_Obj * const objv[], double *minmax) {
00054 int num, i;
00055 Tcl_Obj **data, **vecmin, **vecmax;
00056 
00057 // get the list containing minmax
00058 if (Tcl_ListObjGetElements(interp, objv[arg_minmax], &num, &data) != TCL_OK) {
00059 Tcl_SetResult(interp, (char *)"volmap: could not read parameter (-minmax)", TCL_STATIC);
00060 return TCL_ERROR;
00061 }
00062 if (num != 2) {
00063 Tcl_SetResult(interp, (char *)"volmap: minmax requires a list with two vectors (-minmax)", TCL_STATIC);
00064 return TCL_ERROR;
00065 }
00066 // get the list containing min
00067 if (Tcl_ListObjGetElements(interp, data[0], &num, &vecmin) != TCL_OK) {
00068 return TCL_ERROR;
00069 }
00070 if (num != 3) {
00071 Tcl_SetResult(interp, (char *)"volmap: the first vector does not contain 3 elements (-minmax)", TCL_STATIC);
00072 return TCL_ERROR;
00073 }
00074 // get the list containing max
00075 if (Tcl_ListObjGetElements(interp, data[1], &num, &vecmax) != TCL_OK) {
00076 return TCL_ERROR;
00077 }
00078 if (num != 3) {
00079 Tcl_SetResult(interp, (char *)"volmap: the second vector does not contain 3 elements (-minmax)", TCL_STATIC);
00080 return TCL_ERROR;
00081 }
00082 
00083 // read min
00084 for (i=0; i<3; i++)
00085 if (Tcl_GetDoubleFromObj(interp, vecmin[i], minmax+i) != TCL_OK)
00086 return TCL_ERROR;
00087 
00088 // read max
00089 for (i=0; i<3; i++)
00090 if (Tcl_GetDoubleFromObj(interp, vecmax[i], minmax+i+3) != TCL_OK)
00091 return TCL_ERROR;
00092 
00093 // Check that range is valid...
00094 if (minmax[0] >= minmax[3] || minmax[1] >= minmax[4] || minmax[2] >= minmax[5]) {
00095 Tcl_SetResult(interp, (char *)"volmap: invalid minmax range (-minmax)", TCL_STATIC);
00096 return TCL_ERROR;
00097 }
00098 
00099 return TCL_OK;
00100 }
00101 
00102 
00103 static int vmd_volmap_ils(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00104 bool bad_usage = false;
00105 bool bailout = false;
00106 double cutoff = 10.;
00107 double resolution = 1.;
00108 double minmax[6];
00109 bool pbc = false;
00110 bool pbcbox = false;
00111 float *pbccenter = NULL;
00112 
00113 int maskonly = 0;
00114 bool export_to_file = false;
00115 int molid = -1;
00116 char *filebase = NULL;
00117 char *filename = NULL;
00118 
00119 int nprobecoor = 0;
00120 int nprobevdw = 0;
00121 int nprobecharge = 0;
00122 float *probe_vdwrmin = NULL;
00123 float *probe_vdweps = NULL;
00124 float *probe_coords = NULL;
00125 float *probe_charge = NULL;
00126 double temperature = 300.0;
00127 double maxenergy = 150.0;
00128 int subres = 3;
00129 int num_conf = 0; // number of ligand rotamers
00130 AtomSel *probesel = NULL;
00131 AtomSel *alignsel = NULL;
00132 Matrix4 *transform = NULL;
00133 
00134 if (argc<3) {
00135 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<molid> <minmax> [options...]");
00136 return TCL_ERROR;
00137 }
00138 
00139 // Get the molecule ID
00140 if (!strcmp(Tcl_GetStringFromObj(objv[1], NULL), "top"))
00141 molid = app->molecule_top(); 
00142 else 
00143 Tcl_GetIntFromObj(interp, objv[1], &molid);
00144 
00145 if (!app->molecule_valid_id(molid)) {
00146 Tcl_AppendResult(interp, "volmap: molecule specified for ouput is invalid. (-mol)", NULL);
00147 return TCL_ERROR;
00148 }
00149 
00150 // Get the minmax box or the pbcbox keyword
00151 if (!strcmp(Tcl_GetStringFromObj(objv[2], NULL), "pbcbox")) {
00152 pbcbox = true;
00153 }
00154 else if (parse_minmax_args(interp, 2, objv, minmax)!=TCL_OK) {
00155 Tcl_AppendResult(interp, "volmap: no atoms selected.", NULL);
00156 }
00157 
00158 int first = 0;
00159 int nframes = app->molecule_numframes(molid);
00160 int last = nframes-1;
00161 
00162 // Parse optional arguments
00163 int arg;
00164 for (arg=3; arg<argc && !bad_usage && !bailout; arg++) {
00165 // Arguments common to all volmap types
00166 if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-res")) {
00167 if (arg+1>=argc) { bad_usage=true; break; }
00168 Tcl_GetDoubleFromObj(interp, objv[arg+1], &resolution);
00169 if (resolution <= 0.f) {
00170 Tcl_AppendResult(interp, "volmap ils: resolution must be positive. (-res)", NULL);
00171 bailout = true; break;
00172 }
00173 arg++;
00174 }
00175 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-probesel")) {
00176 if (arg+1>=argc) { bad_usage=true; break; }
00177 // Get atom selection for probe
00178 probesel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[arg+1],NULL));
00179 if (!probesel) {
00180 Tcl_AppendResult(interp, "volmap ils -probesel: no atom selection.", NULL);
00181 bailout = true; break;
00182 }
00183 if (!probesel->selected) {
00184 Tcl_AppendResult(interp, "volmap ils -probesel: no atoms selected.", NULL);
00185 bailout = true; break;
00186 }
00187 if (!app->molecule_valid_id(probesel->molid())) {
00188 Tcl_AppendResult(interp, "volmap ils -probesel: ",
00189 measure_error(MEASURE_ERR_NOMOLECULE), NULL);
00190 bailout = true; break;
00191 }
00192 arg++;
00193 }
00194 // ILS specific arguments
00195 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-first")) {
00196 if (arg+1>=argc) { bad_usage=true; break; }
00197 Tcl_GetIntFromObj(interp, objv[arg+1], &first);
00198 if (first < 0) {
00199 Tcl_AppendResult(interp, "volmap ils: Frame specified with -first must be positive.", NULL);
00200 bailout = true; break;
00201 }
00202 if (first >= nframes) {
00203 Tcl_AppendResult(interp, "volmap ils: Frame specified with -first exceeds number of existing frames.", NULL);
00204 bailout = true; break;
00205 }
00206 arg++;
00207 }
00208 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-last")) {
00209 if (arg+1>=argc) { bad_usage=true; break; }
00210 Tcl_GetIntFromObj(interp, objv[arg+1], &last);
00211 if (last < 0) {
00212 Tcl_AppendResult(interp, "volmap ils: Frame specified with -last must be positive.", NULL);
00213 bailout = true; break;
00214 }
00215 if (last >= nframes) {
00216 Tcl_AppendResult(interp, "volmap ils: Frame specified with -last exceeds number of existing frames.", NULL);
00217 bailout = true; break;
00218 }
00219 arg++;
00220 }
00221 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-cutoff")) {
00222 if (arg+1 >= argc) { bad_usage=true; break; }
00223 Tcl_GetDoubleFromObj(interp, objv[arg+1], &cutoff);
00224 if (cutoff <= 0.) {
00225 Tcl_AppendResult(interp, "volmap ils: cutoff must be positive. (-cutoff)", NULL);
00226 bailout = true; break;
00227 }
00228 arg++;
00229 }
00230 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-dx") ||
00231 !strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-o")) {
00232 if (arg+1>=argc) { bad_usage=true; break; }
00233 const char* outputfilename = Tcl_GetString(objv[arg+1]);
00234 if (outputfilename) {
00235 filebase = new char[strlen(outputfilename)+1];
00236 strcpy(filebase, outputfilename);
00237 export_to_file = true;
00238 }
00239 arg++;
00240 }
00241 
00242 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-probecoor")) {
00243 if (arg+1>=argc) { bad_usage=true; break; }
00244 int i;
00245 double tmp;
00246 Tcl_Obj **listObj;
00247 if (Tcl_ListObjGetElements(interp, objv[arg+1], &nprobecoor, &listObj) != TCL_OK) {
00248 Tcl_AppendResult(interp, " volmap ils: bad syntax in probecoor!", NULL);
00249 bailout = true; break;
00250 }
00251 if (!nprobecoor) {
00252 Tcl_AppendResult(interp, " volmap ils: Empty probecoor list!", NULL);
00253 bailout = true; break;
00254 }
00255 
00256 probe_coords = new float[3L*nprobecoor];
00257 
00258 for (i=0; i<nprobecoor; i++) {
00259 Tcl_Obj **coorObj;
00260 int j, ndim = 0;
00261 if (Tcl_ListObjGetElements(interp, listObj[i], &ndim, &coorObj) != TCL_OK) {
00262 Tcl_AppendResult(interp, " volmap ils: bad syntax in probecoor!", NULL);
00263 bailout = true; break;
00264 }
00265 
00266 if (ndim!=3) {
00267 Tcl_AppendResult(interp, " volmap ils: need three values for each probecoor vector", NULL);
00268 bailout = true; break;
00269 }
00270 
00271 for (j=0; j<3; j++) {
00272 if (Tcl_GetDoubleFromObj(interp, coorObj[j], &tmp) != TCL_OK) {
00273 Tcl_AppendResult(interp, " volmap ils: non-numeric in probecoor", NULL);
00274 bailout = true; break;
00275 }
00276 probe_coords[3L*i+j] = (float)tmp;
00277 }
00278 }
00279 arg++;
00280 }
00281 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-probevdw")) {
00282 if (arg+1>=argc) { bad_usage=true; break; }
00283 double tmp;
00284 Tcl_Obj **listObj;
00285 if (Tcl_ListObjGetElements(interp, objv[arg+1], &nprobevdw, &listObj) != TCL_OK) {
00286 Tcl_AppendResult(interp, " volmap ils: bad syntax in probevdw", NULL);
00287 bailout = true; break;
00288 }
00289 
00290 probe_vdweps = new float[nprobevdw];
00291 probe_vdwrmin = new float[nprobevdw];
00292 
00293 int i;
00294 for (i=0; i<nprobevdw; i++) {
00295 Tcl_Obj **vdwObj;
00296 int ndim = 0;
00297 if (Tcl_ListObjGetElements(interp, listObj[i], &ndim, &vdwObj) != TCL_OK) {
00298 Tcl_AppendResult(interp, " volmap ils: bad syntax in probevdw", NULL);
00299 bailout = true; break;
00300 }
00301 
00302 if (ndim!=2) {
00303 Tcl_AppendResult(interp, " volmap ils: Need two probevdw values (eps, rmin) for each atom", NULL);
00304 bailout = true; break;
00305 }
00306 
00307 if (Tcl_GetDoubleFromObj(interp, vdwObj[0], &tmp) != TCL_OK) {
00308 Tcl_AppendResult(interp, " volmap ils: Non-numeric in probevdw (eps)", NULL);
00309 bailout = true; break;
00310 }
00311 probe_vdweps[i] = (float)tmp;
00312 
00313 if (Tcl_GetDoubleFromObj(interp, vdwObj[1], &tmp) != TCL_OK) {
00314 Tcl_AppendResult(interp, " volmap ils: Non-numeric in probevdw (rmin)", NULL);
00315 bailout = true; break;
00316 }
00317 probe_vdwrmin[i] = (float)tmp; 
00318 }
00319 arg++;
00320 }
00321 
00322 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-probecharge")) {
00323 if (arg+1>=argc) { bad_usage=true; break; }
00324 int i;
00325 double tmp;
00326 Tcl_Obj **listObj;
00327 if (Tcl_ListObjGetElements(interp, objv[arg+1], &nprobecharge, &listObj) != TCL_OK) {
00328 Tcl_AppendResult(interp, " volmap ils: bad syntax in probecharge", NULL);
00329 bailout = true; break;
00330 }
00331 
00332 probe_charge = new float[nprobecharge];
00333 
00334 for (i=0; i<nprobecharge; i++) {
00335 if (Tcl_GetDoubleFromObj(interp, listObj[0], &tmp) != TCL_OK) {
00336 Tcl_AppendResult(interp, " volmap ils: non-numeric in probecharge", NULL);
00337 bailout = true; break;
00338 }
00339 probe_charge[i] = (float)tmp;
00340 }
00341 arg++;
00342 }
00343 
00344 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-orient") ||
00345 !strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-conf")) {
00346 if (arg+1>=argc) { bad_usage=true; break; }
00347 Tcl_GetIntFromObj(interp, objv[arg+1], &num_conf);
00348 if (num_conf < 0) {
00349 Tcl_AppendResult(interp, "volmap ils: invalid -orient parameter", NULL);
00350 bailout = true; break;
00351 }
00352 arg++;
00353 }
00354 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-subres")) {
00355 if (arg+1>=argc) { bad_usage=true; break; }
00356 Tcl_GetIntFromObj(interp, objv[arg+1], &subres);
00357 if (subres < 1) {
00358 Tcl_AppendResult(interp, "volmap ils: invalid -subres parameter", NULL);
00359 bailout = true; break;
00360 }
00361 arg++;
00362 }
00363 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-T")) {
00364 if (arg+1>=argc) { bad_usage=true; break; }
00365 Tcl_GetDoubleFromObj(interp, objv[arg+1], &temperature);
00366 arg++;
00367 } 
00368 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-maxenergy")) {
00369 if (arg+1>=argc) { bad_usage=true; break; }
00370 if (Tcl_GetDoubleFromObj(interp, objv[arg+1], &maxenergy) != TCL_OK) {
00371 Tcl_AppendResult(interp, "volmap ils: invalid -maxenergy parameter", NULL);
00372 bailout = true; break;
00373 }
00374 arg++;
00375 } 
00376 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-alignsel")) {
00377 if (arg+1>=argc) { bad_usage=true; break; }
00378 alignsel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[arg+1],NULL));
00379 if (!alignsel) {
00380 Tcl_AppendResult(interp, "volmap ils: no atom selection for alignment.", NULL);
00381 bailout = true; break;
00382 }
00383 if (!alignsel->selected) {
00384 Tcl_AppendResult(interp, "volmap ils: no atoms selected for alignment.", NULL);
00385 bailout = true; break;
00386 }
00387 arg++;
00388 } 
00389 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-transform")) {
00390 if (arg+1>=argc) { bad_usage=true; break; }
00391 transform = new Matrix4;
00392 tcl_get_matrix("volmap ils: ", interp, objv[arg+1], transform->mat);
00393 arg++;
00394 }
00395 
00396 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-pbc")) {
00397 pbc = true;
00398 }
00399 
00400 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-pbccenter")) {
00401 if (arg+1>=argc) { bad_usage=true; break; }
00402 int i, ndim = 0;
00403 double tmp;
00404 Tcl_Obj **centerObj;
00405 if (Tcl_ListObjGetElements(interp, objv[arg+1], &ndim, &centerObj) != TCL_OK) {
00406 Tcl_AppendResult(interp, " volmap ils: bad syntax in pbccenter", NULL);
00407 bailout = true; break;
00408 }
00409 
00410 if (ndim!=3) {
00411 Tcl_AppendResult(interp, " volmap ils: need three values for vector pbccenter", NULL);
00412 bailout = true; break;
00413 }
00414 
00415 pbccenter = new float[3];
00416 for (i=0; i<3; i++) {
00417 if (Tcl_GetDoubleFromObj(interp, centerObj[i], &tmp) != TCL_OK) {
00418 Tcl_AppendResult(interp, " volmap ils: non-numeric in pbccenter", NULL);
00419 bailout = true; break;
00420 }
00421 pbccenter[i] = (float)tmp;
00422 }
00423 arg++;
00424 pbc = true;
00425 }
00426 
00427 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-maskonly")) {
00428 maskonly = 1;
00429 }
00430 
00431 else {
00432 // unknown arg
00433 Tcl_AppendResult(interp, " volmap ils: unknown argument ",
00434 Tcl_GetStringFromObj(objv[arg], NULL), NULL);
00435 bailout = true;
00436 }
00437 
00438 }
00439 
00440 if (bad_usage) {
00441 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<molid> <minmax> [options...]");
00442 bailout = true;
00443 }
00444 
00445 int num_probe_atoms = 0;
00446 if (!bailout) {
00447 if (probesel) {
00448 if (nprobecoor) {
00449 Tcl_AppendResult(interp, "volmap ils: Must specify either -probesel or -probecoor not both.", NULL);
00450 bailout = true;
00451 }
00452 
00453 // The probe was specified in form of an atomselection
00454 Molecule *probemol = app->moleculeList->mol_from_id(probesel->molid());
00455 const float *radius = probemol->extraflt.data("radius");
00456 if (!radius) {
00457 Tcl_AppendResult(interp, "volmap ils: probe selection contains no VDW radii", NULL);
00458 bailout = true;
00459 }
00460 
00461 const float *occupancy = probemol->extraflt.data("occupancy");
00462 if (!occupancy) {
00463 Tcl_AppendResult(interp, "volmap ils: probe selection contains no VDW epsilon values (in occupancy field)", NULL);
00464 bailout = true;
00465 }
00466 
00467 if (!bailout) {
00468 const float *charge = probemol->extraflt.data("charge");
00469 
00470 num_probe_atoms = probesel->selected;
00471 if (!nprobevdw) {
00472 probe_vdwrmin = new float[num_probe_atoms];
00473 probe_vdweps = new float[num_probe_atoms];
00474 }
00475 probe_charge = new float[num_probe_atoms];
00476 probe_coords = new float[3L*num_probe_atoms];
00477 const float *coords = probesel->coordinates(app->moleculeList);
00478 int i, j=0;
00479 for (i=0; i<probesel->num_atoms; i++) {
00480 if (probesel->on[i]) {
00481 vec_copy(&probe_coords[3L*j], &coords[3L*i]);
00482 
00483 if (!nprobevdw) {
00484 probe_vdweps[j] = -fabsf(occupancy[i]);
00485 probe_vdwrmin[j] = radius[i];
00486 }
00487 if (!nprobecharge) {
00488 if (charge) probe_charge[j] = charge[i];
00489 else probe_charge[j] = 0.f;
00490 }
00491 j++;
00492 }
00493 } 
00494 }
00495 
00496 } else {
00497 // The probe was specified through a coordinate list
00498 if (nprobecoor==0 && nprobevdw==1) {
00499 // No need to specify coodinates of monoatomic probe
00500 num_probe_atoms = 1;
00501 probe_coords = new float[3];
00502 probe_coords[0] = probe_coords[1] = probe_coords[2] = 0.f;
00503 } else {
00504 num_probe_atoms = nprobecoor;
00505 }
00506 
00507 if (!nprobevdw) {
00508 Tcl_AppendResult(interp, "volmap ils: No probe VDW parameters specified.", NULL);
00509 bailout = true; 
00510 }
00511 
00512 if (nprobecharge && nprobecharge!=num_probe_atoms && !bailout) {
00513 Tcl_AppendResult(interp, "volmap ils: # probe charges doesn't match # probe atoms", NULL);
00514 bailout = true;
00515 }
00516 }
00517 
00518 if (num_probe_atoms==0 && !bailout) {
00519 Tcl_AppendResult(interp, "volmap: No probe coordinates specified.", NULL);
00520 bailout = true;
00521 }
00522 
00523 if (nprobevdw && nprobevdw!=num_probe_atoms && !bailout) {
00524 Tcl_AppendResult(interp, "volmap ils: # probe VDW params doesn't match # probe atoms", NULL);
00525 bailout = true;
00526 }
00527 
00528 if (pbc && !alignsel) {
00529 Tcl_AppendResult(interp, "volmap ils: You cannot use -pbc without also "
00530 " providing -alignsel.", NULL);
00531 bailout = true;
00532 }
00533 
00534 //if (num_probe_atoms==1 && num_conf>1 && !bailout) {
00535 //Tcl_AppendResult(interp, "volmap: Specifying -orient for monoatomic probes makes no sense.", NULL);
00536 //bailout = true;
00537 //}
00538 }
00539 
00540 if (bailout) {
00541 if (transform) delete transform;
00542 if (pbccenter) delete [] pbccenter;
00543 if (filebase) delete [] filebase;
00544 if (probe_vdwrmin) delete [] probe_vdwrmin;
00545 if (probe_vdweps) delete [] probe_vdweps;
00546 if (probe_charge) delete [] probe_charge;
00547 if (probe_coords) delete [] probe_coords;
00548 return TCL_ERROR;
00549 }
00550 
00551 // If the probe was provided in form of an atom selection
00552 // determine the symmetry of the probe. We need to know if
00553 // the probe has a tetrahedral point group, the highest
00554 // rotary axis and any C2 axis orthogonal to it.
00555 // If the probe was specified through parameter list instead
00556 // of a selection then the ILS code will try a simple
00557 // symmetry check itself
00558 int tetrahedral_symm = 0;
00559 int order1=0, order2=0;
00560 float symmaxis1[3];
00561 float symmaxis2[3];
00562 if (probesel) {
00563 msgInfo << "Determining probe symmetry:" << sendmsg;
00564 
00565 // Create Symmetry object, verbosity level = 1
00566 Symmetry sym = Symmetry(probesel, app->moleculeList, 1);
00567 
00568 // Set tolerance for atomic overlapping
00569 sym.set_overlaptol(0.05f);
00570 
00571 // Take bonds into account
00572 sym.set_checkbonds(1);
00573 
00574 // Guess the symmetry
00575 int ret_val = sym.guess(0.05f);
00576 if (ret_val < 0) {
00577 Tcl_AppendResult(interp, "measure symmetry: ", measure_error(ret_val), NULL);
00578 return TCL_ERROR;
00579 }
00580 int pgorder;
00581 char pointgroup[6];
00582 sym.get_pointgroup(pointgroup, &pgorder);
00583 
00584 if (pointgroup[0]=='T') tetrahedral_symm = 1;
00585 
00586 if (sym.numaxes()) {
00587 // First symmetry axis
00588 vec_copy(symmaxis1, sym.axis(0));
00589 order1 = sym.get_axisorder(0);
00590 
00591 int i;
00592 for (i=1; i<sym.numaxes(); i++) {
00593 vec_copy(symmaxis2, sym.axis(i));
00594 if (fabs(dot_prod(symmaxis1, symmaxis2)) < DEGTORAD(1.f)) {
00595 // Orthogonal second axis found
00596 order2 = sym.get_axisorder(i);
00597 break;
00598 }
00599 }
00600 }
00601 }
00602 
00603 // 6. Create the volmap
00604 VolMapCreateILS vol(app, molid, first, last, (float)temperature,
00605 (float)resolution, subres, (float)cutoff,
00606 maskonly);
00607 
00608 vol.set_probe(num_probe_atoms, num_conf, probe_coords,
00609 probe_vdwrmin, probe_vdweps, probe_charge);
00610 vol.set_maxenergy(float(maxenergy));
00611 
00612 if (pbc) vol.set_pbc(pbccenter, pbcbox);
00613 if (transform) vol.set_transform(transform);
00614 if (alignsel) vol.set_alignsel(alignsel);
00615 if (!pbcbox) {
00616 vol.set_minmax(float(minmax[0]), float(minmax[1]), float(minmax[2]),
00617 float(minmax[3]), float(minmax[4]), float(minmax[5]));
00618 }
00619 
00620 if (tetrahedral_symm || order1 || order2) {
00621 // Provide info about probe symmetry
00622 vol.set_probe_symmetry(order1, symmaxis1, order2, symmaxis2, tetrahedral_symm);
00623 }
00624 
00625 // Create map...
00626 int ret_val = vol.compute();
00627 
00628 if (ret_val < 0) {
00629 Tcl_AppendResult(interp, "\nvolmap: ERROR ", measure_error(ret_val), NULL);
00630 
00631 if (transform) delete transform;
00632 if (pbccenter) delete [] pbccenter;
00633 if (filebase) delete [] filebase;
00634 if (probe_vdwrmin) delete [] probe_vdwrmin;
00635 if (probe_vdweps) delete [] probe_vdweps;
00636 if (probe_charge) delete [] probe_charge;
00637 if (probe_coords) delete [] probe_coords;
00638 return TCL_ERROR;
00639 }
00640 
00641 int numconf, numorient, numrot;
00642 vol.get_statistics(numconf, numorient, numrot);
00643 
00644 // If the probe was specified in form of a selection and a separate
00645 // molecule then we append a frame for each conformer and set the
00646 // probe coordinates accordingly.
00647 if (probesel && probesel->molid()!=molid) {
00648 float *conformers = NULL;
00649 int numconf = vol.get_conformers(conformers);
00650 Molecule *pmol = app->moleculeList->mol_from_id(probesel->molid());
00651 int i, j;
00652 for (i=0; i<numconf; i++) {
00653 if (i>0) {
00654 pmol->duplicate_frame(pmol->current());
00655 }
00656 float *coor = pmol->get_frame(i)->pos;
00657 int k=0;
00658 for (j=0; j<probesel->num_atoms; j++) { 
00659 if (!probesel->on[j]) continue; //atom is not selected
00660 
00661 vec_copy(&coor[3L*j], &conformers[i*3L*num_probe_atoms+3L*k]);
00662 k++;
00663 }
00664 }
00665 }
00666 
00667 // Export volmap to a file or just add it to the molecule:
00668 if (export_to_file) {
00669 // Add .dx suffix to filebase if it is missing
00670 filename = new char[strlen(filebase)+16];
00671 strcpy(filename, filebase);
00672 char *suffix = strrchr(filename, '.'); // beginning of .dx
00673 if (strcmp(suffix, ".dx")) strcat(filename, ".dx");
00674 
00675 // Write tha map into a dx file
00676 if (!vol.write_map(filename)) {
00677 Tcl_AppendResult(interp, "\nvolmap: ERROR Could not write ils map to file", NULL);
00678 }
00679 
00680 delete[] filename;
00681 
00682 } else {
00683 if (!vol.add_map_to_molecule()) {
00684 Tcl_AppendResult(interp, "\nvolmap: ERROR Could not add ils map to molecule", NULL);
00685 }
00686 }
00687 
00688 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00689 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("numconf", -1));
00690 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(numconf));
00691 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("numorient", -1));
00692 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(numorient));
00693 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("numrot", -1));
00694 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(numrot));
00695 Tcl_SetObjResult(interp, tcl_result);
00696 
00697 if (transform) delete transform;
00698 if (pbccenter) delete [] pbccenter;
00699 if (filebase) delete [] filebase;
00700 if (probe_vdwrmin) delete [] probe_vdwrmin;
00701 if (probe_vdweps) delete [] probe_vdweps;
00702 if (probe_charge) delete [] probe_charge;
00703 if (probe_coords) delete [] probe_coords;
00704 
00705 return TCL_OK;
00706 }
00707 
00708 
00709 static int vmd_volmap_new_fromtype(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00710 
00711 // 1. Figure out which map type we are dealing with:
00712 enum {UNDEF_MAP, DENS_MAP, INTERP_MAP, DIST_MAP, OCCUP_MAP, MASK_MAP,
00713 CPOTENTIAL_MAP, CPOTENTIALMSM_MAP } maptype=UNDEF_MAP;
00714 
00715 int ret_val = 0;
00716 
00717 char *maptype_string=Tcl_GetString(objv[0]); 
00718 
00719 if (!strcmp(maptype_string, "density")) maptype=DENS_MAP;
00720 else if (!strcmp(maptype_string, "interp")) maptype=INTERP_MAP;
00721 else if (!strcmp(maptype_string, "distance")) maptype=DIST_MAP;
00722 else if (!strcmp(maptype_string, "occupancy")) maptype=OCCUP_MAP; 
00723 else if (!strcmp(maptype_string, "mask")) maptype=MASK_MAP; 
00724 else if (!strcmp(maptype_string, "coulomb")) maptype=CPOTENTIAL_MAP;
00725 else if (!strcmp(maptype_string, "coulombpotential")) maptype=CPOTENTIAL_MAP;
00726 else if (!strcmp(maptype_string, "coulombmsm")) maptype=CPOTENTIALMSM_MAP;
00727 
00728 // 2. Get atom selection
00729 if (argc < 2) { 
00730 Tcl_AppendResult(interp, "volmap: no atom selection.", NULL);
00731 return TCL_ERROR;
00732 }
00733 
00734 AtomSel *sel = NULL;
00735 sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
00736 if (!sel) {
00737 Tcl_AppendResult(interp, "volmap: no atom selection.", NULL);
00738 return TCL_ERROR;
00739 }
00740 if (!sel->selected) {
00741 Tcl_AppendResult(interp, "volmap: no atoms selected.", NULL);
00742 return TCL_ERROR;
00743 }
00744 if (!app->molecule_valid_id(sel->molid())) {
00745 Tcl_AppendResult(interp, "volmap: ",
00746 measure_error(MEASURE_ERR_NOMOLECULE), NULL);
00747 return TCL_ERROR;
00748 }
00749 
00750 
00751 // 3. Define and initialize the optional arguments
00752 bool accept_weight = false; // allow a user-specified weight for each atom
00753 bool accept_cutoff = false; // parse a user-specified cutoff distance
00754 bool accept_radius = false; // allow radius multiplicator for density
00755 bool accept_usepoints = false; // allows use of point particles
00756 
00757 bool use_point_particles = false; // for MASK map
00758 bool export_to_file = false;
00759 bool use_all_frames = false;
00760 bool bad_usage = false;
00761 
00762 int export_molecule = -1;
00763 double cutoff = 5.;
00764 double radius_factor = 1.;
00765 double resolution = 1.;
00766 double minmax[6];
00767 
00768 char *filebase = NULL;
00769 char *filename = NULL;
00770 
00771 
00772 // File export options
00773 int checkpoint_freq = 500;
00774 
00775 
00776 // Specify required/accepted options for each maptype as well as default values.
00777 switch(maptype) {
00778 case DENS_MAP:
00779 accept_weight = true;
00780 accept_radius = true;
00781 break;
00782 case INTERP_MAP:
00783 accept_weight = true; 
00784 break;
00785 case DIST_MAP:
00786 accept_cutoff = true;
00787 cutoff = 3.;
00788 break;
00789 case OCCUP_MAP:
00790 accept_usepoints = true;
00791 break;
00792 case MASK_MAP:
00793 accept_cutoff = true;
00794 cutoff = 4.;
00795 break;
00796 case CPOTENTIAL_MAP:
00797 case CPOTENTIALMSM_MAP:
00798 break; 
00799 case UNDEF_MAP:
00800 bad_usage = true;
00801 break; 
00802 }
00803 
00804 
00805 // 4. Parse the command-line
00806 int arg_weight=0, arg_combine=0, arg_minmax=0;
00807 
00808 // Parse the arguments
00809 for (int arg=2; arg<argc && !bad_usage; arg++) {
00810 if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-res")) {
00811 if (arg+1>=argc) bad_usage=true;
00812 Tcl_GetDoubleFromObj(interp, objv[arg+1], &resolution);
00813 if (resolution <= 0.f) {
00814 Tcl_AppendResult(interp, "volmap: resolution must be positive. (-res)", NULL);
00815 return TCL_ERROR;
00816 }
00817 arg++;
00818 }
00819 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-mol")) { // add volmap to mol
00820 if (arg+1 >= argc) bad_usage=true;
00821 if (!strcmp(Tcl_GetStringFromObj(objv[arg+1], NULL), "top"))
00822 export_molecule = app->molecule_top(); 
00823 else 
00824 Tcl_GetIntFromObj(interp, objv[arg+1], &export_molecule);
00825 
00826 if (!app->molecule_valid_id(export_molecule)) {
00827 Tcl_AppendResult(interp, "volmap: molecule specified for ouput is invalid. (-mol)", NULL);
00828 return TCL_ERROR;
00829 }
00830 arg++;
00831 }
00832 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-minmax")) {
00833 arg_minmax=arg+1;
00834 arg++;
00835 if (arg_minmax>=argc) bad_usage=true;
00836 if (parse_minmax_args(interp, arg_minmax, objv, minmax)!=TCL_OK) {
00837 return TCL_ERROR;
00838 }
00839 }
00840 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-checkpoint")) {
00841 if (arg+1 >= argc) bad_usage=true;
00842 Tcl_GetIntFromObj(interp, objv[arg+1], &checkpoint_freq);
00843 if (checkpoint_freq < 0) {
00844 Tcl_AppendResult(interp, "volmap: invalid -checkpoint parameter", NULL);
00845 return TCL_ERROR;
00846 }
00847 arg++;
00848 }
00849 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-allframes")) {
00850 use_all_frames=true;
00851 }
00852 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-dx") ||
00853 !strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-o")) {
00854 if (arg+1>=argc) {bad_usage=true; break;}
00855 const char* outputfilename = Tcl_GetString(objv[arg+1]);
00856 if (outputfilename) {
00857 filebase = new char[strlen(outputfilename)+1];
00858 strcpy(filebase, outputfilename);
00859 export_to_file = true;
00860 }
00861 arg++;
00862 }
00863 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-combine")) {
00864 arg_combine=arg+1;
00865 arg++;
00866 if (arg_combine>=argc) bad_usage=true;
00867 }
00868 else if (accept_usepoints && !strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-points")) {
00869 use_point_particles=true;
00870 }
00871 else if (accept_cutoff && !strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-cutoff")) {
00872 if (arg+1 >= argc) bad_usage=true;
00873 Tcl_GetDoubleFromObj(interp, objv[arg+1], &cutoff);
00874 if (cutoff <= 0.) {
00875 Tcl_AppendResult(interp, "volmap: cutoff must be positive. (-cutoff)", NULL);
00876 return TCL_ERROR;
00877 }
00878 arg++;
00879 }
00880 else if (accept_radius && !strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-radscale")) {
00881 if (arg+1 >= argc) bad_usage=true;
00882 Tcl_GetDoubleFromObj(interp, objv[arg+1], &radius_factor);
00883 if (radius_factor < 0.f) {
00884 Tcl_AppendResult(interp, "volmap: radscale must be positive. (-radscale)", NULL);
00885 return TCL_ERROR;
00886 }
00887 arg++;
00888 }
00889 else if (accept_weight && !strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-weight")) {
00890 if (arg+1>=argc) bad_usage=true;
00891 arg_weight = arg+1; 
00892 arg++;
00893 }
00894 
00895 else //unknown arg
00896 bad_usage=true; 
00897 }
00898 
00899 
00900 if (bad_usage) {
00901 if (maptype == UNDEF_MAP)
00902 Tcl_AppendResult(interp, "volmap: unknown map type.", NULL);
00903 else
00904 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<selection> [options...]");
00905 
00906 return TCL_ERROR;
00907 }
00908 
00909 
00910 // 5. Assign some other optional parameters
00911 
00912 // parse map combination type
00913 VolMapCreate::CombineType combine_type=VolMapCreate::COMBINE_AVG;
00914 if (arg_combine) {
00915 char *combine_str=Tcl_GetString(objv[arg_combine]);
00916 if (!strcmp(combine_str, "avg") || !strcmp(combine_str, "average")) 
00917 combine_type=VolMapCreate::COMBINE_AVG;
00918 else if (!strcmp(combine_str, "max") || !strcmp(combine_str, "maximum")) 
00919 combine_type=VolMapCreate::COMBINE_MAX;
00920 else if (!strcmp(combine_str, "min") || !strcmp(combine_str, "minimum")) 
00921 combine_type=VolMapCreate::COMBINE_MIN;
00922 else if (!strcmp(combine_str, "stdev")) 
00923 combine_type=VolMapCreate::COMBINE_STDEV;
00924 else if (!strcmp(combine_str, "pmf")) 
00925 combine_type=VolMapCreate::COMBINE_PMF;
00926 else {
00927 Tcl_AppendResult(interp, "volmap: -combine argument must be: avg, min, \
00928 max, stdev, pmf", NULL);
00929 return TCL_ERROR;
00930 }
00931 }
00932 
00933 // parse weights
00934 float *weights = NULL;
00935 char *weight_string = NULL;
00936 int weight_mutable = 1;
00937 if (accept_weight) {
00938 weights = new float[sel->num_atoms];
00939 if (arg_weight) {
00940 weight_string = Tcl_GetStringFromObj(objv[arg_weight], NULL);
00941 if (!weight_string || !strcmp(weight_string, "none")) {
00942 ret_val = atomsel_default_weights(sel, weights);
00943 weight_string = NULL;
00944 } else {
00945 // Check if the argument is a VMD attribute
00946 int fctn = get_attribute_index(app, weight_string);
00947 if (fctn >= 0) {
00948 ret_val = get_weights_from_attribute(app, sel, weight_string,
00949 weights);
00950 if (ret_val == MEASURE_ERR_BADWEIGHTPARM) {
00951 Tcl_AppendResult(interp,
00952 "weight attribute must have floating point values",
00953 NULL);
00954 }
00955 } else {
00956 ret_val = get_weights_from_tcl_list(interp, app, sel,
00957 objv[arg_weight], weights);
00958 weight_string = NULL;
00959 weight_mutable = 0; // The Tcl list will not change
00960 }
00961 }
00962 } else {
00963 ret_val = atomsel_default_weights(sel, weights);
00964 }
00965 
00966 if (ret_val < 0) {
00967 Tcl_AppendResult(interp, "volmap: ", measure_error(ret_val), NULL);
00968 delete [] weights;
00969 return TCL_ERROR;
00970 }
00971 }
00972 
00973 
00974 
00975 // 6. Create the volmap
00976 VolMapCreate *volcreate = NULL; 
00977 
00978 // init the map creator objects and set default filenames
00979 switch(maptype) {
00980 case DENS_MAP: 
00981 volcreate = new VolMapCreateDensity(app, sel, (float)resolution,
00982 weights, weight_string,
00983 weight_mutable, (float)radius_factor);
00984 if (!filebase) {
00985 filebase = new char[strlen("density_out.dx")+1];
00986 strcpy(filebase, "density_out.dx");
00987 }
00988 break;
00989 
00990 case INTERP_MAP: 
00991 volcreate = new VolMapCreateInterp(app, sel, (float)resolution,
00992 weights, weight_string,
00993 weight_mutable);
00994 if (!filebase) {
00995 filebase = new char[strlen("interp_out.dx")+1];
00996 strcpy(filebase, "interp_out.dx");
00997 }
00998 break;
00999 
01000 case DIST_MAP:
01001 volcreate = new VolMapCreateDistance(app, sel, (float)resolution, (float)cutoff);
01002 if (!filebase) {
01003 filebase = new char[strlen("distance_out.dx")+1];
01004 strcpy(filebase, "distance_out.dx");
01005 }
01006 break;
01007 
01008 case OCCUP_MAP:
01009 volcreate = new VolMapCreateOccupancy(app, sel, (float)resolution, use_point_particles);
01010 if (!filebase) {
01011 filebase = new char[strlen("occupancy_out.dx")+1];
01012 strcpy(filebase, "occupancy_out.dx");
01013 }
01014 break;
01015 
01016 case MASK_MAP:
01017 volcreate = new VolMapCreateMask(app, sel, (float)resolution, (float)cutoff);
01018 if (!filebase) {
01019 filebase = new char[strlen("mask_out.dx")+1];
01020 strcpy(filebase, "mask_out.dx");
01021 }
01022 break;
01023 
01024 case CPOTENTIAL_MAP:
01025 volcreate = new VolMapCreateCoulombPotential(app, sel, (float)resolution);
01026 if (!filebase) {
01027 filebase = new char[strlen("coulomb_out.dx")+1];
01028 strcpy(filebase, "coulomb_out.dx");
01029 }
01030 break;
01031 
01032 case CPOTENTIALMSM_MAP:
01033 volcreate = new VolMapCreateCoulombPotentialMSM(app, sel, (float)resolution);
01034 if (!filebase) {
01035 filebase = new char[strlen("coulombmsm_out.dx")+1];
01036 strcpy(filebase, "coulombmsm_out.dx");
01037 }
01038 break;
01039 
01040 
01041 default: // silence compiler warnings
01042 break; 
01043 }
01044 
01045 // generate and write out volmap
01046 if (volcreate) {
01047 // Pass parameters common to all volmap types
01048 if (arg_minmax)
01049 volcreate->set_minmax(float(minmax[0]), float(minmax[1]),
01050 float(minmax[2]), float(minmax[3]), 
01051 float(minmax[4]), float(minmax[5]));
01052 
01053 // Setup checkpointing
01054 if (checkpoint_freq) {
01055 char *checkpointname = new char[32+strlen(filebase)+1];
01056 #if defined(_MSC_VER)
01057 char slash = '\\';
01058 #else
01059 char slash = '/';
01060 #endif
01061 char *tailname = strrchr(filebase, slash);
01062 if (!tailname) tailname = filebase;
01063 else tailname = tailname+1;
01064 char *dirname = new char[strlen(filebase)+1];
01065 strcpy(dirname, filebase);
01066 char *sep = strrchr(dirname, slash);
01067 
01068 if (sep) {
01069 *sep = '0円';
01070 sprintf(checkpointname, "%s%ccheckpoint:%s", dirname, slash, tailname);
01071 }
01072 else {
01073 *dirname = '0円';
01074 sprintf(checkpointname, "checkpoint:%s", tailname);
01075 }
01076 delete[] dirname;
01077 
01078 Tcl_AppendResult(interp, "CHECKPOINTNAME = ", checkpointname, NULL);
01079 volcreate->set_checkpoint(checkpoint_freq, checkpointname);
01080 delete[] checkpointname;
01081 }
01082 
01083 // Create map...
01084 ret_val = volcreate->compute_all(use_all_frames, combine_type, NULL);
01085 if (ret_val < 0) {
01086 delete volcreate;
01087 if (weights) delete [] weights;
01088 if (filebase) delete [] filebase;
01089 
01090 Tcl_AppendResult(interp, "\nvolmap: ERROR ", measure_error(ret_val), NULL);
01091 return TCL_ERROR;
01092 }
01093 
01094 // Export volmap to a DX file:
01095 if (export_to_file || export_molecule < 0) {
01096 // add .dx suffix to filebase if it is missing
01097 filename = new char[strlen(filebase)+16];
01098 strcpy(filename, filebase);
01099 char *suffix = strrchr(filename, '.');
01100 if ((suffix != NULL) && !strcmp(suffix,".dx")) 
01101 *suffix = '0円';
01102 strcat(filename, ".dx");
01103 volcreate->write_map(filename);
01104 delete[] filename;
01105 }
01106 
01107 // Export volmap to a molecule:
01108 if (export_molecule >= 0) {
01109 VolumetricData *volmap = volcreate->volmap;
01110 float origin[3], xaxis[3], yaxis[3], zaxis[3];
01111 int i;
01112 for (i=0; i<3; i++) {
01113 origin[i] = (float) volmap->origin[i];
01114 xaxis[i] = (float) volmap->xaxis[i];
01115 yaxis[i] = (float) volmap->yaxis[i];
01116 zaxis[i] = (float) volmap->zaxis[i];
01117 }
01118 int err = app->molecule_add_volumetric(export_molecule, 
01119 (volmap->name) ? volmap->name : "(no name)",
01120 origin, xaxis, yaxis, zaxis,
01121 volmap->xsize, volmap->ysize, volmap->zsize, volmap->data);
01122 if (err != 1) {
01123 Tcl_AppendResult(interp, "ERROR: export of volmap into molecule was unsuccessful!", NULL);
01124 }
01125 else volmap->data=NULL; // avoid data being deleted by volmap's destructor (it is now owned by the molecule)
01126 }
01127 
01128 delete volcreate;
01129 }
01130 
01131 if (weights) delete [] weights;
01132 if (filebase) delete [] filebase;
01133 
01134 return TCL_OK;
01135 }
01136 
01137 
01138 // vec_sub() from utilities.h works with float* only
01139 // here I needed doubles.
01140 #define DOUBLE_VSUB(D, X, Y) \
01141 D[0] = float(X[0]-Y[0]); \
01142 D[1] = float(X[1]-Y[1]); \
01143 D[2] = float(X[2]-Y[2]); 
01144 
01145 // Compare two volumetric maps:
01146 // volmap compare <molid1> <volid1> <molid2> <volid2>
01147 // The two maps must be specified by their molID and volID.
01148 // Prints the min/max vales, the largest difference, the RMSD,
01149 // the RMSD computed inly for the elements that differ and the
01150 // RMSD weighted by the magnitude of the elements compared so
01151 // that smaller values receive a larger weight.
01152 // (For ILS we are interested mainly in the smaller energies)
01153 static int vmd_volmap_compare(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp)
01154 {
01155 int molid1 = -1;
01156 int molid2 = -1;
01157 int mapid1 = -1;
01158 int mapid2 = -1;
01159 
01160 // Get the molecule IDs
01161 if (!strcmp(Tcl_GetStringFromObj(objv[1], NULL), "top"))
01162 molid1 = app->molecule_top(); 
01163 else 
01164 Tcl_GetIntFromObj(interp, objv[1], &molid1);
01165 
01166 if (!strcmp(Tcl_GetStringFromObj(objv[3], NULL), "top"))
01167 molid2 = app->molecule_top(); 
01168 else 
01169 Tcl_GetIntFromObj(interp, objv[3], &molid2);
01170 
01171 if (!app->molecule_valid_id(molid1) || !app->molecule_valid_id(molid2)) {
01172 Tcl_AppendResult(interp, "volmap compare: molecule specified for ouput is invalid. (-mol)", NULL);
01173 return TCL_ERROR;
01174 }
01175 
01176 Molecule *mol1 = app->moleculeList->mol_from_id(molid1);
01177 Molecule *mol2 = app->moleculeList->mol_from_id(molid2);
01178 
01179 // Get volmap IDs
01180 Tcl_GetIntFromObj(interp, objv[2], &mapid1);
01181 Tcl_GetIntFromObj(interp, objv[4], &mapid2);
01182 
01183 if (mapid1<0 || mapid2<0) {
01184 Tcl_AppendResult(interp, "volmap compare: Volmap ID must be positive.", NULL);
01185 return TCL_ERROR;
01186 }
01187 if (mapid1 >= mol1->num_volume_data() ||
01188 mapid2 >= mol2->num_volume_data()) {
01189 Tcl_AppendResult(interp, "volmap compare: Volmap ID does not exist.", NULL);
01190 return TCL_ERROR;
01191 }
01192 
01193 // Parse optional arguments
01194 bool bad_usage = false;
01195 double histinterval = 2.5;
01196 int numbins = 9;
01197 int arg;
01198 for (arg=5; arg<argc && !bad_usage; arg++) {
01199 if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-interval")) {
01200 if (arg+1>=argc) { bad_usage=true; break; }
01201 Tcl_GetDoubleFromObj(interp, objv[arg+1], &histinterval);
01202 if (histinterval <= 0.f) {
01203 Tcl_AppendResult(interp, "volmap compare: histogram interval must be positive. (-interval)", NULL);
01204 return TCL_ERROR;
01205 }
01206 arg++;
01207 }
01208 else if (!strcmp(Tcl_GetStringFromObj(objv[arg], NULL), "-numbins")) {
01209 if (arg+1>=argc) { bad_usage=true; break; }
01210 Tcl_GetIntFromObj(interp, objv[arg+1], &numbins);
01211 if (numbins <= 0) {
01212 Tcl_AppendResult(interp, "volmap compare: histogram bin size must be positive. (-interval)", NULL);
01213 return TCL_ERROR;
01214 }
01215 arg++;
01216 }
01217 else {
01218 // unknown arg
01219 Tcl_AppendResult(interp, " volmap compare: unknown argument ",
01220 Tcl_GetStringFromObj(objv[arg], NULL), NULL);
01221 return TCL_ERROR;
01222 }
01223 
01224 }
01225 if (bad_usage) {
01226 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<molid> <minmax> [options...]");
01227 return TCL_ERROR;
01228 }
01229 
01230 const VolumetricData *vol1 = mol1->get_volume_data(mapid1);
01231 const VolumetricData *vol2 = mol2->get_volume_data(mapid2);
01232 
01233 if (vol1->xsize != vol2->xsize ||
01234 vol1->ysize != vol2->ysize ||
01235 vol1->zsize != vol2->zsize) {
01236 Tcl_AppendResult(interp, "volmap compare: maps have different dimensions.", NULL);
01237 return TCL_ERROR;
01238 }
01239 
01240 float vdiff[3];
01241 DOUBLE_VSUB(vdiff, vol1->origin, vol2->origin);
01242 if (norm(vdiff)>1e-6) {
01243 Tcl_AppendResult(interp, "volmap compare: maps have different origin.", NULL);
01244 }
01245 
01246 DOUBLE_VSUB(vdiff, vol1->xaxis, vol2->xaxis);
01247 if (norm(vdiff)>1e-6) {
01248 Tcl_AppendResult(interp, "volmap compare: maps have different x-axis.", NULL);
01249 }
01250 DOUBLE_VSUB(vdiff, vol1->yaxis, vol2->yaxis);
01251 if (norm(vdiff)>1e-6) {
01252 Tcl_AppendResult(interp, "volmap compare: maps have different y-axis.", NULL);
01253 }
01254 DOUBLE_VSUB(vdiff, vol1->zaxis, vol2->zaxis);
01255 if (norm(vdiff)>1e-6) {
01256 Tcl_AppendResult(interp, "volmap compare: maps have different z-axis.", NULL);
01257 }
01258 
01259 int i;
01260 int numdiff = 0;
01261 float sqsum = 0.f;
01262 float sqsumd = 0.f;
01263 float min1 = 0.f, min2 = 0.f;
01264 float max1 = 0.f, max2 = 0.f;
01265 float maxdiff = 0.f;
01266 int indexmaxdiff = 0;
01267 
01268 for (i=0; i<vol1->gridsize(); i++) {
01269 float v1 = vol1->data[i];
01270 float v2 = vol2->data[i];
01271 float diff = v1-v2;
01272 sqsum += diff*diff;
01273 if (v1<min1) min1 = v1;
01274 if (v2<min2) min2 = v2;
01275 if (v1>max1) max1 = v1;
01276 if (v2>max2) max2 = v2;
01277 if (fabsf(diff)>maxdiff) {
01278 maxdiff = fabsf(diff);
01279 indexmaxdiff = i;
01280 }
01281 if (diff) {
01282 //printf("%g - %g = %g\n", v1, v2, diff);
01283 sqsumd += diff*diff;
01284 numdiff++;
01285 }
01286 }
01287 float rmsd = 0.f;
01288 float diffrmsd = 0.f;
01289 if (sqsum) rmsd = sqrtf(sqsum/vol1->gridsize());
01290 if (sqsumd) diffrmsd = sqrtf(sqsumd/numdiff);
01291 
01292 char tmpstr[2048];
01293 msgInfo << "volmap compare" << sendmsg;
01294 msgInfo << "--------------" << sendmsg;
01295 msgInfo << "Comparing mol "<<molid1<<" -> map "<<mapid1<<"/"<<mol1->num_volume_data()<<sendmsg;
01296 msgInfo << " to mol "<<molid2<<" -> map "<<mapid2<<"/"<<mol2->num_volume_data()<<sendmsg;
01297 msgInfo << "Statistics:" << sendmsg;
01298 sprintf(tmpstr, " %12s | %12s", "MAP 1", "MAP 2");
01299 msgInfo << tmpstr << sendmsg;
01300 sprintf(tmpstr, "min value = %12g | %12g", min1, min2);
01301 msgInfo << tmpstr << sendmsg;
01302 sprintf(tmpstr, "max value = %12g | %12g", max1, max2);
01303 msgInfo << tmpstr << sendmsg;
01304 msgInfo << sendmsg;
01305 sprintf(tmpstr, "# differing elements = %d/%ld", numdiff, vol1->gridsize());
01306 msgInfo << tmpstr << sendmsg;
01307 msgInfo << "max difference:" << sendmsg;
01308 sprintf(tmpstr, " map1[%d] = %g map2[%d] = %g diff = %g",
01309 indexmaxdiff, vol1->data[indexmaxdiff],
01310 indexmaxdiff, vol2->data[indexmaxdiff], maxdiff);
01311 msgInfo << tmpstr << sendmsg;
01312 
01313 // Statistics for the differing elements only:
01314 msgInfo << sendmsg;
01315 sprintf(tmpstr, " RMSD = %12.6f", rmsd);
01316 msgInfo << tmpstr << sendmsg;
01317 sprintf(tmpstr, " diffRMSD = %12.6f (for the set of differing elements)", diffrmsd);
01318 msgInfo << tmpstr << sendmsg;
01319 
01320 // Get weighted RMSD where the differences of smaller values
01321 // are weighted more because the low enegy values are what is 
01322 // of interest in ILS free energy maps.
01323 float wsum = 0.f;
01324 float range = max2-min2;
01325 float wdiffrmsd = 0.f;
01326 if (range) {
01327 sqsum = 0.f;
01328 for (i=0; i<vol1->gridsize(); i++) {
01329 float diff = vol1->data[i]-vol2->data[i];
01330 if (diff) {
01331 float weight = 1.f-(vol2->data[i]-min2)/range;
01332 wsum += weight;
01333 sqsum += diff*diff*weight;
01334 }
01335 }
01336 if (sqsum) wdiffrmsd = sqrtf(sqsum/wsum);
01337 }
01338 
01339 
01340 sprintf(tmpstr, "weighted RMSD = %12.6f (for the set of differing elements)", wdiffrmsd);
01341 msgInfo << tmpstr << sendmsg;
01342 msgInfo << " weight factor w = 1-(E_i-E_min)/(E_max-E_min)" << sendmsg;
01343 msgInfo << sendmsg;
01344 
01345 // Compare error of map 1 relative to map 2, create histogram of error:
01346 // | E_approx - E_exact | / ( E_exact - E_exact_min + 1 ),
01347 // where map 1 is considered the approximation and map 2 is considered exact.
01348 //
01349 // Since energy is arbitrary, we shift from the minimum value 
01350 // (usually around -11 to -6) so that the dominator is always
01351 // greater than or equal to 1. This weights the lower values
01352 // more heavily than the upper values, which is intentional.
01353 //
01354 // The histogram is summed for the intervals
01355 // (-\inf,10) [10,20) [20,30) [30,40) [40,+\inf)
01356 
01357 
01358 // total accumulated error for each bin
01359 float *histo = new float[numbins*sizeof(float)];
01360 // count number of samples per bin
01361 int *num = new int[numbins*sizeof(float)];
01362 // max error for each bin
01363 float *maxEntry = new float[numbins*sizeof(float)];
01364 float *binrmsd = new float[numbins*sizeof(float)];
01365 memset(histo, 0, numbins*sizeof(float));
01366 memset(num, 0, numbins*sizeof(int));
01367 memset(maxEntry, 0, numbins*sizeof(float));
01368 memset(binrmsd, 0, numbins*sizeof(float));
01369 
01370 for (i = 0; i < vol1->gridsize(); i++) {
01371 float e1 = vol1->data[i];
01372 float e2 = vol2->data[i];
01373 float err = fabsf(e1 - e2) / (e2 - min2 + 1);
01374 int index = (int) floorf((e2 - min2) / float(histinterval));
01375 if (index < 0) index = 0;
01376 else if (index >= numbins) index = numbins - 1;
01377 
01378 // check to see if we need to update the max for this bin
01379 if (err > maxEntry[index]) { 
01380 maxEntry[index] = err; 
01381 }
01382 histo[index] += err;
01383 num[index]++;
01384 binrmsd[index] += (e2-e1)*(e2-e1);
01385 }
01386 for (i=0; i<numbins; i++) {
01387 if (binrmsd[i]) binrmsd[i] = sqrtf(binrmsd[i]/num[i]);
01388 }
01389 
01390 // lower boundary of the first bin
01391 float firstbin = floorf(min2/float(histinterval))*float(histinterval);
01392 Tcl_Obj *caption = Tcl_NewListObj(0, NULL);
01393 Tcl_Obj *numEntries = Tcl_NewListObj(0, NULL);
01394 Tcl_Obj *objHisto = Tcl_NewListObj(0, NULL);
01395 Tcl_Obj *objAverage = Tcl_NewListObj(0, NULL);
01396 Tcl_Obj *objMax = Tcl_NewListObj(0, NULL);
01397 Tcl_Obj *objBinRMSD = Tcl_NewListObj(0, NULL);
01398 char label[64];
01399 msgInfo << "Histogram of error in map 1 relative to map 2 " << sendmsg;
01400 msgInfo << sendmsg;
01401 msgInfo << " interval # samples total error avg error"
01402 << " max error rmsd" << sendmsg;
01403 msgInfo << "---------------------------------------------------------"
01404 << "----------------------------" << sendmsg;
01405 
01406 sprintf(label, "[%g,%g)", (double) firstbin, (double) (firstbin+histinterval));
01407 sprintf(tmpstr, "%14s %10d %14.6f %14.6f %14.6f %14.6f", label, num[0],
01408 (double) histo[0], (double) (!num[0]?0:histo[0]/num[0]),
01409 (double) maxEntry[0], (double) binrmsd[0]);
01410 msgInfo << tmpstr << sendmsg;
01411 Tcl_ListObjAppendElement(interp, caption, Tcl_NewStringObj(label, -1));
01412 Tcl_ListObjAppendElement(interp, numEntries, Tcl_NewIntObj(num[0]));
01413 Tcl_ListObjAppendElement(interp, objHisto, Tcl_NewDoubleObj(histo[0]));
01414 Tcl_ListObjAppendElement(interp, objAverage, Tcl_NewDoubleObj(!num[0]?0:histo[0]/num[0]));
01415 Tcl_ListObjAppendElement(interp, objMax, Tcl_NewDoubleObj(maxEntry[0]));
01416 Tcl_ListObjAppendElement(interp, objBinRMSD, Tcl_NewDoubleObj(binrmsd[0]));
01417 for (i = 1; i < numbins-1; i++) {
01418 sprintf(label, "[%g,%g)", (double)(firstbin+i*histinterval), (double)(firstbin+(i+1)*histinterval));
01419 sprintf(tmpstr, "%14s %10d %14.6f %14.6f %14.6f %14.6f", label, num[i],
01420 (double) histo[i], (double) (!num[i]?0:histo[i]/num[i]),
01421 (double) maxEntry[i], (double) binrmsd[i]);
01422 Tcl_ListObjAppendElement(interp, caption, Tcl_NewStringObj(label, -1));
01423 Tcl_ListObjAppendElement(interp, numEntries, Tcl_NewIntObj(num[i]));
01424 Tcl_ListObjAppendElement(interp, objHisto, Tcl_NewDoubleObj(histo[i]));
01425 Tcl_ListObjAppendElement(interp, objAverage, Tcl_NewDoubleObj(!num[i]?0:histo[i]/num[i]));
01426 Tcl_ListObjAppendElement(interp, objMax, Tcl_NewDoubleObj(maxEntry[i]));
01427 Tcl_ListObjAppendElement(interp, objBinRMSD, Tcl_NewDoubleObj(binrmsd[i]));
01428 msgInfo << tmpstr << sendmsg;
01429 }
01430 sprintf(label, "[%g,+infty)", (double)(firstbin+i*histinterval));
01431 sprintf(tmpstr, "%14s %10d %14.6f %14.6f %14.6f %14.6f", label, num[i],
01432 (double) histo[i], (double) (!num[i]?0:histo[i]/num[i]),
01433 (double) maxEntry[i], (double) binrmsd[i]);
01434 Tcl_ListObjAppendElement(interp, caption, Tcl_NewStringObj(label, -1));
01435 Tcl_ListObjAppendElement(interp, numEntries, Tcl_NewIntObj(num[i]));
01436 Tcl_ListObjAppendElement(interp, objHisto, Tcl_NewDoubleObj(histo[i]));
01437 Tcl_ListObjAppendElement(interp, objAverage, Tcl_NewDoubleObj(!num[i]?0:histo[i]/num[i]));
01438 Tcl_ListObjAppendElement(interp, objMax, Tcl_NewDoubleObj(maxEntry[i]));
01439 Tcl_ListObjAppendElement(interp, objBinRMSD, Tcl_NewDoubleObj(binrmsd[i]));
01440 msgInfo << tmpstr << sendmsg;
01441 msgInfo << sendmsg;
01442 
01443 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01444 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("rmsd", -1));
01445 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(rmsd));
01446 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("diffrmsd", -1));
01447 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(diffrmsd));
01448 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("weightedrmsd", -1));
01449 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(wdiffrmsd));
01450 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("maxdiff", -1));
01451 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(maxdiff));
01452 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("numdiff", -1));
01453 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewIntObj(numdiff));
01454 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("caption", -1));
01455 Tcl_ListObjAppendElement(interp, tcl_result, caption);
01456 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("numEntries", -1));
01457 Tcl_ListObjAppendElement(interp, tcl_result, numEntries);
01458 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("histo", -1));
01459 Tcl_ListObjAppendElement(interp, tcl_result, objHisto);
01460 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("avgErr", -1));
01461 Tcl_ListObjAppendElement(interp, tcl_result, objAverage);
01462 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("maxError", -1));
01463 Tcl_ListObjAppendElement(interp, tcl_result, objMax);
01464 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("binRMSD", -1));
01465 Tcl_ListObjAppendElement(interp, tcl_result, objBinRMSD);
01466 Tcl_SetObjResult(interp, tcl_result);
01467 
01468 delete [] histo;
01469 delete [] num;
01470 delete [] maxEntry;
01471 delete [] binrmsd;
01472 return TCL_OK;
01473 }
01474 
01475 int obj_volmap(ClientData cd, Tcl_Interp *interp, int argc, Tcl_Obj * const objv[]) {
01476 
01477 if (argc < 2) {
01478 Tcl_SetResult(interp, (char *)
01479 "usage: volmap <command> <args...>\n"
01480 "\nVolmap Creation:\n"
01481 " volmap <maptype> <selection> [opts...] -- create a new volmap file\n"
01482 " maptypes:\n"
01483 " density -- arbitrary-weight density map [atoms/A^3]\n"
01484 " interp -- arbitrary-weight interpolation map [atoms/A^3]\n"
01485 " distance -- distance nearest atom surface [A]\n"
01486 " occupancy -- percent atomic occupancy of gridpoints [%]\n"
01487 " mask -- binary mask by painting spheres around atoms\n"
01488 " coulomb -- Coulomb electrostatic potential [kT/e] (slow)\n"
01489 " coulombmsm -- Coulomb electrostatic potential [kT/e] (fast)\n"
01490 " ils -- free energy map [kT] computed by implicit ligand sampling\n"
01491 " options common to all maptypes:\n"
01492 " -o <filename> -- output DX format file name (use .dx extension)\n"
01493 " -mol <molid> -- export volmap into the specified mol\n"
01494 " -res <float> -- resolution in A of smallest cube\n"
01495 " -allframes -- compute for all frames of the trajectory\n"
01496 " -combine <arg> -- rule for combining the different frames\n"
01497 " <arg> = avg, min, max, stdev or pmf\n"
01498 " -minmax <list of 2 vectors> -- specify boundary of output grid\n"
01499 " options specific to certain maptypes:\n"
01500 " -points -- use point particles for occupancy\n"
01501 " -cutoff <float> -- distance cutoff for calculations [A]\n"
01502 " -radscale <float> -- premultiply all atomic radii by a factor\n"
01503 " -weight <str/list> -- per atom weights for calculation\n"
01504 " options for ils:\n"
01505 " see documentation\n", NULL);
01506 return TCL_ERROR;
01507 }
01508 
01509 VMDApp *app = (VMDApp *)cd;
01510 char *arg1 = Tcl_GetStringFromObj(objv[1],NULL);
01511 
01512 // If maptype is recognized, proceed with the map creation (vs. yet-unimplemented map operations)...
01513 if (argc > 1 && !strupncmp(arg1, "occupancy", CMDLEN))
01514 return vmd_volmap_new_fromtype(app, argc-1, objv+1, interp);
01515 if (argc > 1 && !strupncmp(arg1, "density", CMDLEN))
01516 return vmd_volmap_new_fromtype(app, argc-1, objv+1, interp);
01517 if (argc > 1 && !strupncmp(arg1, "interp", CMDLEN))
01518 return vmd_volmap_new_fromtype(app, argc-1, objv+1, interp);
01519 if (argc > 1 && !strupncmp(arg1, "distance", CMDLEN))
01520 return vmd_volmap_new_fromtype(app, argc-1, objv+1, interp);
01521 if (argc > 1 && !strupncmp(arg1, "mask", CMDLEN))
01522 return vmd_volmap_new_fromtype(app, argc-1, objv+1, interp);
01523 if (argc > 1 && (!strupncmp(arg1, "coulombpotential", CMDLEN) || 
01524 !strupncmp(arg1, "coulomb", CMDLEN) ||
01525 !strupncmp(arg1, "coulombmsm", CMDLEN)))
01526 return vmd_volmap_new_fromtype(app, argc-1, objv+1, interp); 
01527 
01528 if (argc > 1 && !strupncmp(arg1, "compare", CMDLEN))
01529 return vmd_volmap_compare(app, argc-1, objv+1, interp);
01530 
01531 if (argc > 1 && !strupncmp(arg1, "ils", CMDLEN))
01532 return vmd_volmap_ils(app, argc-1, objv+1, interp);
01533 if (argc > 1 && !strupncmp(arg1, "ligand", CMDLEN))
01534 return vmd_volmap_ils(app, argc-1, objv+1, interp);
01535 
01536 Tcl_SetResult(interp, (char *)"Type 'volmap' for summary of usage\n",NULL);
01537 return TCL_ERROR;
01538 }