TclMeasure.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: TclMeasure.C,v $
00013 * $Author: dgomes $ $Locker: $ $State: Exp $
00014 * $Revision: 1.186 $ $Date: 2024年05月16日 14:56:56 $
00015 *
00016 ***************************************************************************
00017 * DESCRIPTION:
00018 * These are essentially just Tcl wrappers for the measure commands in
00019 * Measure.C.
00020 *
00021 ***************************************************************************/
00022 
00023 #include <stdlib.h>
00024 #include <tcl.h>
00025 #include <ctype.h>
00026 #include <math.h>
00027 #include "TclCommands.h"
00028 #include "AtomSel.h"
00029 #include "Matrix4.h"
00030 #include "SymbolTable.h"
00031 #include "VMDApp.h"
00032 #include "MoleculeList.h"
00033 #include "utilities.h"
00034 #include "config.h"
00035 #include "Measure.h"
00036 #include "MeasureSymmetry.h"
00037 #include "SpatialSearch.h"
00038 #include "Atom.h"
00039 #include "Molecule.h"
00040 
00041 // needed by VolInterior commands
00042 #include "QuickSurf.h"
00043 #include "MDFF.h"
00044 #include "CUDAMDFF.h"
00045 #include "MeasureVolInterior.h"
00046 
00047 
00048 // Get weights from a Tcl obj. Data must hold sel->selected items, or natoms.
00049 // If there is no obj, or if it's "none", return a list of ones.
00050 // If the obj matches an atom selection keyword/field, and the field returns 
00051 // floating point data, return that data, otherwise return error.
00052 // Otherwise, the obj must be a list of floats to use for the weights,
00053 // with either a size == sel->selected, or size == natoms. 
00054 //
00055 // NOTE: this routine cannot be used on selections that change between frames.
00056 //
00057 int tcl_get_weights(Tcl_Interp *interp, VMDApp *app, AtomSel *sel, 
00058 Tcl_Obj *weight_obj, float *data) {
00059 char *weight_string = NULL;
00060 
00061 if (!sel) return MEASURE_ERR_NOSEL;
00062 if (!app->molecule_valid_id(sel->molid())) return MEASURE_ERR_NOMOLECULE;
00063 if (weight_obj)
00064 weight_string = Tcl_GetStringFromObj(weight_obj, NULL);
00065 
00066 if (!weight_string || !strcmp(weight_string, "none")) {
00067 for (int i=0; i<sel->selected; i++) {
00068 data[i] = 1.0;
00069 }
00070 return MEASURE_NOERR;
00071 }
00072 
00073 // if a selection string was given, check the symbol table
00074 SymbolTable *atomSelParser = app->atomSelParser; 
00075 
00076 // weights must return floating point values, so the symbol must not 
00077 // be a singleword, so macro is NULL.
00078 atomsel_ctxt context(atomSelParser, 
00079 app->moleculeList->mol_from_id(sel->molid()), 
00080 sel->which_frame, NULL);
00081 
00082 int fctn = atomSelParser->find_attribute(weight_string);
00083 if (fctn >= 0) {
00084 // the keyword exists, so get the data if it is type float, otherwise fail
00085 if (atomSelParser->fctns.data(fctn)->returns_a != SymbolTableElement::IS_FLOAT) {
00086 Tcl_AppendResult(interp, 
00087 "weight attribute must have floating point values", NULL);
00088 return MEASURE_ERR_BADWEIGHTPARM; // can't understand weight parameter 
00089 }
00090 
00091 double *tmp_data = new double[sel->num_atoms];
00092 atomSelParser->fctns.data(fctn)->keyword_double(&context, 
00093 sel->num_atoms, tmp_data, sel->on);
00094 
00095 for (int i=sel->firstsel, j=0; i<=sel->lastsel; i++) {
00096 if (sel->on[i])
00097 data[j++] = (float)tmp_data[i];
00098 }
00099 
00100 delete [] tmp_data;
00101 return MEASURE_NOERR;
00102 }
00103 
00104 // Determine if weights are a Tcl list with the right number of atoms
00105 int list_num;
00106 Tcl_Obj **list_data;
00107 if (Tcl_ListObjGetElements(interp, weight_obj, &list_num, &list_data) 
00108 != TCL_OK) {
00109 return MEASURE_ERR_BADWEIGHTPARM;
00110 }
00111 if (list_num != sel->selected && list_num != sel->num_atoms) 
00112 return MEASURE_ERR_BADWEIGHTNUM;
00113 
00114 for (int i=0, j=0; i<list_num; i++) {
00115 double tmp_data;
00116 
00117 if (Tcl_GetDoubleFromObj(interp, list_data[i], &tmp_data) != TCL_OK) 
00118 return MEASURE_ERR_NONNUMBERPARM;
00119 
00120 if (list_num == sel->selected) {
00121 data[i] = (float)tmp_data; // one weight list item per selected atom
00122 } else {
00123 if (sel->on[i]) {
00124 data[j++] = (float)tmp_data; // one weight list item for each atom
00125 }
00126 }
00127 }
00128 
00129 return MEASURE_NOERR;
00130 }
00131 
00132 
00133 int atomsel_default_weights(AtomSel *sel, float *weights) {
00134 if (sel->firstsel > 0)
00135 memset(&weights[0], 0, sel->firstsel * sizeof(float)); 
00136 
00137 for (int i=sel->firstsel; i<=sel->lastsel; i++) {
00138 // Use the standard "on" check instead of typecasting the array elements
00139 weights[i] = sel->on[i] ? 1.0f : 0.0f;
00140 }
00141 
00142 if (sel->lastsel < (sel->num_atoms - 1))
00143 memset(&weights[sel->lastsel+1], 0, ((sel->num_atoms - 1) - sel->lastsel) * sizeof(float)); 
00144 
00145 return 0;
00146 }
00147 
00148 
00149 int get_weights_from_tcl_list(Tcl_Interp *interp, VMDApp *app, AtomSel *sel,
00150 Tcl_Obj *weights_obj, float *weights) {
00151 int list_num = 0;
00152 Tcl_Obj **list_elems = NULL;
00153 if (Tcl_ListObjGetElements(interp, weights_obj, &list_num, &list_elems)
00154 != TCL_OK) {
00155 return MEASURE_ERR_BADWEIGHTPARM;
00156 }
00157 if (list_num != sel->num_atoms) {
00158 return MEASURE_ERR_BADWEIGHTNUM;
00159 }
00160 for (int i = 0; i < sel->num_atoms; i++) {
00161 double tmp_data = 0.0;
00162 if (Tcl_GetDoubleFromObj(interp, list_elems[i], &tmp_data) != TCL_OK) {
00163 return TCL_ERROR;
00164 }
00165 weights[i] = static_cast<float>(tmp_data);
00166 }
00167 return 0;
00168 }
00169 
00170 
00171 // This routine eliminates the need to include SymbolTable.h,
00172 // and allows the caller to determine if a per-atom attribute/field
00173 // exists or not.
00174 int get_attribute_index(VMDApp *app, char const *string) {
00175 SymbolTable *atomSelParser = app->atomSelParser;
00176 return atomSelParser->find_attribute(string);
00177 }
00178 
00179 
00180 int get_weights_from_attribute(VMDApp *app, AtomSel *sel,
00181 char const *weights_string, float *weights) {
00182 SymbolTable *atomSelParser = app->atomSelParser;
00183 int fctn = atomSelParser->find_attribute(weights_string);
00184 // first, check to see that the function returns floats.
00185 // if it doesn't, it makes no sense to use it as a weight
00186 if (atomSelParser->fctns.data(fctn)->returns_a !=
00187 SymbolTableElement::IS_FLOAT) {
00188 return MEASURE_ERR_BADWEIGHTPARM; // can't understand weight parameter
00189 }
00190 atomsel_ctxt context(atomSelParser,
00191 app->moleculeList->mol_from_id(sel->molid()),
00192 sel->which_frame, NULL);
00193 double *tmp_data = new double[sel->num_atoms];
00194 atomSelParser->fctns.data(fctn)->keyword_double(&context, sel->num_atoms,
00195 tmp_data, sel->on);
00196 
00197 if (sel->firstsel > 0)
00198 memset(&weights[0], 0, sel->firstsel * sizeof(float)); 
00199 
00200 for (int i=sel->firstsel; i<=sel->lastsel; i++) {
00201 weights[i] = sel->on[i] ? static_cast<float>(tmp_data[i]) : 0.0f;
00202 }
00203 
00204 if (sel->lastsel < (sel->num_atoms - 1))
00205 memset(&weights[sel->lastsel+1], 0, ((sel->num_atoms - 1) - sel->lastsel) * sizeof(float)); 
00206 
00207 delete [] tmp_data;
00208 return 0;
00209 }
00210 
00211 
00212 // get the atom index re-ordering list for use by measure_fit
00213 int tcl_get_orders(Tcl_Interp *interp, int selnum, 
00214 Tcl_Obj *order_obj, int *data) {
00215 int list_num;
00216 Tcl_Obj **list_data;
00217 
00218 if (Tcl_ListObjGetElements(interp, order_obj, &list_num, &list_data)
00219 != TCL_OK) {
00220 return MEASURE_ERR_NOSEL;
00221 }
00222 
00223 // see if this is a Tcl list with the right number of atom indices
00224 if (list_num != selnum) return MEASURE_ERR_NOSEL;
00225 
00226 for (int i=0; i<list_num; i++) {
00227 if (Tcl_GetIntFromObj(interp, list_data[i], &data[i]) != TCL_OK)
00228 return MEASURE_ERR_NONNUMBERPARM;
00229 
00230 // order indices are 0-based
00231 if (data[i] < 0 || data[i] >= selnum)
00232 return MEASURE_ERR_BADORDERINDEX; // order index is out of range
00233 }
00234 
00235 return MEASURE_NOERR;
00236 }
00237 
00238 
00239 //
00240 // Function: vmd_measure_center
00241 // Parameters: <selection> // computes with weight == 1
00242 // Parameters: <selection> weight [none | atom field | list] 
00243 // computes with the weights based on the following:
00244 // none => weights all 1
00245 // atom field => value from atomSelParser (eg
00246 // mass => use weight based on mass
00247 // index => use weight based on atom index (0 to n-1)
00248 // list => use list to get weights for each atom. 
00249 // The list can have length == number of selected atoms,
00250 // or it can have length == the total number of atoms
00251 // 
00252 // Examples: 
00253 // vmd_measure_center atomselect12
00254 // vmd_measure_center atomselect12 weight mass
00255 // vmd_measure_center atomselect12 {12 13} [atomselect top "index 2 3"]
00256 // If no weight is given, no weight term is used (computes center of number)
00257 //
00258 static int vmd_measure_center(VMDApp *app, int argc, Tcl_Obj *const objv[], Tcl_Interp *interp) {
00259 if (argc != 2 && argc != 4 ) {
00260 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<sel> [weight <weights>]");
00261 return TCL_ERROR;
00262 }
00263 if (argc == 4 && strcmp(Tcl_GetStringFromObj(objv[2],NULL), "weight")) {
00264 Tcl_SetResult(interp, (char *) "measure center: parameter can only be 'weight'", TCL_STATIC);
00265 return TCL_ERROR;
00266 }
00267 
00268 // get the selection
00269 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
00270 if (!sel) {
00271 Tcl_SetResult(interp, (char *) "measure center: no atom selection", TCL_STATIC);
00272 return TCL_ERROR;
00273 }
00274 
00275 // get the weight
00276 float *weight = new float[sel->selected];
00277 int ret_val=0;
00278 if (argc == 2) { // only from atom selection, so weight is 1
00279 ret_val = tcl_get_weights(interp, app, sel, NULL, weight);
00280 } else {
00281 ret_val = tcl_get_weights(interp, app, sel, objv[3], weight);
00282 }
00283 if (ret_val < 0) {
00284 Tcl_AppendResult(interp, "measure center: ", measure_error(ret_val), NULL);
00285 delete [] weight;
00286 return TCL_ERROR;
00287 }
00288 
00289 // compute the center of "mass"
00290 float com[3];
00291 const float *framepos = sel->coordinates(app->moleculeList);
00292 ret_val = measure_center(sel, framepos, weight, com);
00293 delete [] weight;
00294 if (ret_val < 0) {
00295 Tcl_AppendResult(interp, "measure center: ", measure_error(ret_val), NULL);
00296 return TCL_ERROR;
00297 }
00298 
00299 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00300 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(com[0]));
00301 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(com[1]));
00302 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(com[2]));
00303 Tcl_SetObjResult(interp, tcl_result);
00304 
00305 return TCL_OK;
00306 }
00307 
00308 
00309 static int vmd_measure_centerperresidue(VMDApp *app, int argc, Tcl_Obj *const objv[], Tcl_Interp *interp) {
00310 int ret_val = 0;
00311 
00312 // check argument counts for valid combinations
00313 if (argc != 2 && argc != 4 ) {
00314 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<sel> [weight <weights>]");
00315 return TCL_ERROR;
00316 }
00317 
00318 // the only valid optional parameter is "weight"
00319 if (argc == 4 && strcmp(Tcl_GetStringFromObj(objv[2], NULL), "weight")) {
00320 Tcl_SetResult(interp, (char *) "measure centerperresidue: parameter can only be 'weight'", TCL_STATIC);
00321 return TCL_ERROR;
00322 }
00323 
00324 // get the selection
00325 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1], NULL));
00326 if (!sel) {
00327 Tcl_SetResult(interp, (char *) "measure centerperresidue: no atom selection", TCL_STATIC);
00328 return TCL_ERROR;
00329 }
00330 
00331 // get the weight
00332 float *weight = new float[sel->selected];
00333 if (argc == 2) {
00334 // only from atom selection, so weights are all set to 1
00335 ret_val = tcl_get_weights(interp, app, sel, NULL, weight);
00336 } else {
00337 // from a per-atom field or a list
00338 ret_val = tcl_get_weights(interp, app, sel, objv[3], weight);
00339 }
00340 if (ret_val < 0) {
00341 Tcl_AppendResult(interp, "measure centerperresidue: ", 
00342 measure_error(ret_val), NULL);
00343 delete [] weight;
00344 return TCL_ERROR;
00345 }
00346 
00347 // compute the center of "mass"
00348 float *com = new float[3*sel->selected];
00349 const float *framepos = sel->coordinates(app->moleculeList);
00350 ret_val = measure_center_perresidue(app->moleculeList, sel, 
00351 framepos, weight, com);
00352 delete [] weight;
00353 if (ret_val < 0) {
00354 Tcl_AppendResult(interp, "measure centerperresidue: ", 
00355 measure_error(ret_val), NULL);
00356 return TCL_ERROR;
00357 }
00358 
00359 // generate lists of CoMs from results
00360 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00361 for (int i=0; i<ret_val; i++) {
00362 Tcl_Obj *m = Tcl_NewListObj(0, NULL);
00363 for (int j=0; j<3; j++) {
00364 Tcl_ListObjAppendElement(interp, m, Tcl_NewDoubleObj(com[3*i+j]));
00365 }
00366 Tcl_ListObjAppendElement(interp, tcl_result, m);
00367 }
00368 Tcl_SetObjResult(interp, tcl_result);
00369 delete [] com;
00370 
00371 return TCL_OK;
00372 }
00373 
00374 
00375 // measure sum of weights for selected atoms
00376 static int vmd_measure_sumweights(VMDApp *app, int argc, Tcl_Obj *const objv[], Tcl_Interp *interp) {
00377 if (argc != 4) {
00378 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<sel> weight <weights>");
00379 return TCL_ERROR;
00380 }
00381 if (strcmp(Tcl_GetStringFromObj(objv[2],NULL), "weight")) {
00382 Tcl_SetResult(interp, (char *) "measure sumweights: parameter can only be 'weight'", TCL_STATIC);
00383 return TCL_ERROR;
00384 }
00385 
00386 // get the selection
00387 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL)
00388 );
00389 if (!sel) {
00390 Tcl_SetResult(interp, (char *) "measure sumweights: no atom selection", TCL_STATIC);
00391 return TCL_ERROR;
00392 }
00393 
00394 // get the weight
00395 float *weight = new float[sel->selected];
00396 int ret_val = 0;
00397 ret_val = tcl_get_weights(interp, app, sel, objv[3], weight);
00398 if (ret_val < 0) {
00399 Tcl_AppendResult(interp, "measure center: ", measure_error(ret_val), NULL);
00400 delete [] weight;
00401 return TCL_ERROR;
00402 }
00403 
00404 // compute the sum of the weights
00405 float weightsum=0;
00406 ret_val = measure_sumweights(sel, sel->selected, weight, &weightsum);
00407 delete [] weight;
00408 if (ret_val < 0) {
00409 Tcl_AppendResult(interp, "measure center: ", measure_error(ret_val), NULL);
00410 return TCL_ERROR;
00411 }
00412 
00413 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00414 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(weightsum));
00415 Tcl_SetObjResult(interp, tcl_result);
00416 
00417 return TCL_OK;
00418 }
00419 
00420 
00421 // Function: vmd_measure_avpos <selection> first <first> last <last> step <step>
00422 // Returns: the average position of the selected atoms over the selected frames
00423 // Example: measure avpos atomselect76 0 20 1
00424 static int vmd_measure_avpos(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00425 int first = 0; // start with first frame by default
00426 int last = -1; // finish with last frame by default
00427 int step = 1; // use all frames by default
00428 
00429 if (argc < 2 || argc > 8) {
00430 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<sel> [first <first>] [last <last>] [step <step>]");
00431 return TCL_ERROR;
00432 }
00433 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL)
00434 );
00435 if (!sel) {
00436 Tcl_AppendResult(interp, "measure avpos: no atom selection", NULL);
00437 return TCL_ERROR;
00438 }
00439 
00440 int i;
00441 for (i=2; i<argc; i+=2) {
00442 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
00443 if (!strupncmp(argvcur, "first", CMDLEN)) {
00444 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
00445 Tcl_AppendResult(interp, "measure avpos: bad first frame value", NULL);
00446 return TCL_ERROR;
00447 }
00448 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
00449 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
00450 Tcl_AppendResult(interp, "measure avpos: bad last frame value", NULL);
00451 return TCL_ERROR;
00452 }
00453 } else if (!strupncmp(argvcur, "step", CMDLEN)) {
00454 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK) {
00455 Tcl_AppendResult(interp, "measure avpos: bad frame step value", NULL);
00456 return TCL_ERROR;
00457 }
00458 } else {
00459 Tcl_AppendResult(interp, "measure rmsdmat_qcp: invalid syntax, no such keyword: ", argvcur, NULL);
00460 return TCL_ERROR;
00461 }
00462 }
00463 
00464 float *avpos = new float[3L*sel->selected];
00465 int ret_val = measure_avpos(sel, app->moleculeList, first, last, step, avpos);
00466 if (ret_val < 0) {
00467 Tcl_AppendResult(interp, "measure avpos: ", measure_error(ret_val), NULL);
00468 delete [] avpos;
00469 return TCL_ERROR;
00470 }
00471 
00472 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00473 for (i=0; i<sel->selected; i++) {
00474 Tcl_Obj *atom = Tcl_NewListObj(0, NULL);
00475 Tcl_ListObjAppendElement(interp, atom, Tcl_NewDoubleObj(avpos[i*3L ]));
00476 Tcl_ListObjAppendElement(interp, atom, Tcl_NewDoubleObj(avpos[i*3L + 1]));
00477 Tcl_ListObjAppendElement(interp, atom, Tcl_NewDoubleObj(avpos[i*3L + 2]));
00478 Tcl_ListObjAppendElement(interp, tcl_result, atom);
00479 }
00480 
00481 Tcl_SetObjResult(interp, tcl_result);
00482 delete [] avpos;
00483 
00484 return TCL_OK;
00485 }
00486 
00487 
00488 // Function: vmd_measure_dipole <selection>
00489 // Returns: the dipole moment for the selected atoms
00490 static int vmd_measure_dipole(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00491 const char *opt;
00492 int unitsdebye=0; // default units are elementary charges/Angstrom
00493 int usecenter=1; // remove net charge at the center of mass (-1), geometrical center (1), don't (0)
00494 
00495 if ((argc < 2) || (argc > 4)) {
00496 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<sel> [-elementary|-debye] [-geocenter|-masscenter|-origincenter]");
00497 return TCL_ERROR;
00498 }
00499 AtomSel *sel;
00500 sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1], NULL));
00501 if (!sel) {
00502 Tcl_AppendResult(interp, "measure dipole: no atom selection", NULL);
00503 return TCL_ERROR;
00504 }
00505 
00506 int i;
00507 for (i=0; i < (argc-2); ++i) {
00508 opt = Tcl_GetStringFromObj(objv[2+i], NULL);
00509 if (!strcmp(opt, "-debye"))
00510 unitsdebye=1; 
00511 if (!strcmp(opt, "-elementary"))
00512 unitsdebye=0; 
00513 
00514 if (!strcmp(opt, "-geocenter"))
00515 usecenter=1; 
00516 if (!strcmp(opt, "-masscenter"))
00517 usecenter=-1; 
00518 if (!strcmp(opt, "-origincenter"))
00519 usecenter=0; 
00520 }
00521 
00522 float dipole[3];
00523 int ret_val = measure_dipole(sel, app->moleculeList, dipole, unitsdebye, usecenter);
00524 if (ret_val < 0) {
00525 Tcl_AppendResult(interp, "measure dipole: ", measure_error(ret_val), NULL);
00526 return TCL_ERROR;
00527 }
00528 
00529 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00530 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(dipole[0]));
00531 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(dipole[1]));
00532 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(dipole[2]));
00533 Tcl_SetObjResult(interp, tcl_result);
00534 
00535 return TCL_OK;
00536 }
00537 
00538 
00539 
00540 // Function: vmd_measure_dihed {4 atoms as {<atomid> ?<molid>?}} ?molid <default molid>?
00541 // ?frame [f|all|last]? | ?first <first>? ?last <last>?
00542 // Returns: the dihedral angle for the specified atoms
00543 static int vmd_measure_dihed(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00544 int first=-1, last=-1, frame=-1;
00545 if (argc < 2) {
00546 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "{{<atomid1> [<molid1>]} {<atomid2> [<molid2>]} {<atomid3> [<molid3>]} {<atomid4> [<molid4>]}} [molid <default molid>] [frame <frame|all|last> | first <first> last <last>]");
00547 return TCL_ERROR;
00548 }
00549 
00550 int molid[4], atmid[4], defmolid = -1;
00551 bool allframes = false;
00552 
00553 // Get the geometry type dihed/imprp
00554 char *geomname = Tcl_GetStringFromObj(objv[0],NULL);
00555 
00556 // Read the atom list
00557 int numatms;
00558 Tcl_Obj **data;
00559 if (Tcl_ListObjGetElements(interp, objv[1], &numatms, &data) != TCL_OK) {
00560 Tcl_AppendResult(interp, " measure ", geomname, ": bad syntax", NULL);
00561 Tcl_AppendResult(interp, " Usage: measure ", geomname, " {{<atomid1> [<molid1>]} {<atomid2> [<molid2>]} {<atomid3> [<molid3>]} {<atomid4> [<molid4>]}} [molid <default molid>] [frame <frame|all|last> | first <first> last <last>]", NULL);
00562 return TCL_ERROR;
00563 }
00564 
00565 if (numatms != 4) {
00566 Tcl_AppendResult(interp, " measure dihed: must specify exactly four atoms in a list", NULL);
00567 return TCL_ERROR;
00568 }
00569 
00570 if (argc > 3) {
00571 int i;
00572 for (i=2; i<argc; i+=2) {
00573 char *argvcur = Tcl_GetStringFromObj(objv[i], NULL);
00574 if (!strupncmp(argvcur, "molid", CMDLEN)) {
00575 if (Tcl_GetIntFromObj(interp, objv[i+1], &defmolid) != TCL_OK) {
00576 Tcl_AppendResult(interp, " measure ", geomname, ": bad molid", NULL);
00577 return TCL_ERROR;
00578 }
00579 } else if (!strupncmp(argvcur, "first", CMDLEN)) {
00580 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
00581 Tcl_AppendResult(interp, " measure ", geomname, ": bad first frame value", NULL);
00582 return TCL_ERROR;
00583 }
00584 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
00585 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
00586 Tcl_AppendResult(interp, " measure ", geomname, ": bad last frame value", NULL);
00587 return TCL_ERROR;
00588 }
00589 } else if (!strupncmp(argvcur, "frame", CMDLEN)) {
00590 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "all")) {
00591 allframes = true;
00592 } else if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "last")) {
00593 frame=-2;
00594 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &frame) != TCL_OK) {
00595 Tcl_AppendResult(interp, " measure ", geomname, ": bad frame value", NULL);
00596 return TCL_ERROR;
00597 }
00598 } else {
00599 Tcl_AppendResult(interp, "measure ", geomname, ": invalid syntax, no such keyword: ", argvcur, NULL);
00600 return TCL_ERROR;
00601 }
00602 }
00603 }
00604 
00605 if ((allframes || frame>=0) && (first>=0 || last>=0)) {
00606 Tcl_AppendResult(interp, "measure ", geomname, ": Ambiguous syntax: You cannot specify a frame AND a frame range (using first or last).", NULL);
00607 Tcl_AppendResult(interp, "\nUsage:\nmeasure ", geomname, " {<atomid1> [<molid1>]} {<atomid2> [<molid2>]} {<atomid3> [<molid3>]} {<atomid4> [<molid4>]}} [molid <default molid>] [frame <frame|all|last> | first <first> last <last>]", NULL);
00608 return TCL_ERROR; 
00609 }
00610 
00611 if (allframes) first=0;
00612 
00613 // If no molecule was specified use top as default
00614 if (defmolid<0) defmolid = app->molecule_top();
00615 
00616 // Assign atom IDs and molecule IDs
00617 int i,numelem;
00618 Tcl_Obj **atmmol;
00619 for (i=0; i<numatms; i++) {
00620 if (Tcl_ListObjGetElements(interp, data[i], &numelem, &atmmol) != TCL_OK) {
00621 return TCL_ERROR;
00622 }
00623 
00624 if (!numelem) {
00625 Tcl_AppendResult(interp, " measure ", geomname, ": empty atom index", NULL);
00626 return TCL_ERROR;
00627 }
00628 
00629 if (Tcl_GetIntFromObj(interp, atmmol[0], atmid+i) != TCL_OK) {
00630 Tcl_AppendResult(interp, " measure ", geomname, ": bad atom index", NULL);
00631 return TCL_ERROR;
00632 }
00633 
00634 if (numelem==2) {
00635 if (Tcl_GetIntFromObj(interp, atmmol[1], molid+i) != TCL_OK) {
00636 Tcl_AppendResult(interp, " measure ", geomname, ": bad molid", NULL);
00637 return TCL_ERROR;
00638 }
00639 } else molid[i] = defmolid;
00640 }
00641 
00642 // Compute the value
00643 ResizeArray<float> gValues(1024);
00644 int ret_val;
00645 ret_val = measure_geom(app->moleculeList, molid, atmid, &gValues, 
00646 frame, first, last, defmolid, MEASURE_DIHED);
00647 if (ret_val < 0) {
00648 Tcl_AppendResult(interp, measure_error(ret_val), NULL);
00649 return TCL_ERROR;
00650 }
00651 
00652 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00653 int numvalues = gValues.num();
00654 for (int count = 0; count < numvalues; count++) {
00655 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(gValues[count]));
00656 }
00657 Tcl_SetObjResult(interp, tcl_result);
00658 
00659 return TCL_OK;
00660 }
00661 
00662 
00663 // Function: vmd_measure_angle {3 atoms as {<atomid> ?<molid>?}} ?molid <default molid>? 
00664 // ?frame [f|all|last]? | ?first <first>? ?last <last>?
00665 // Returns: the bond angle for the specified atoms
00666 static int vmd_measure_angle(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00667 int first=-1, last=-1, frame=-1;
00668 
00669 if(argc<2) {
00670 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "{{<atomid1> [<molid1>]} {<atomid2> [<molid2>]} {<atomid3> [<molid3>]}} [molid <default molid>] [frame <frame|all|last> | first <first> last <last>]");
00671 return TCL_ERROR;
00672 }
00673 
00674 int molid[3], atmid[3], defmolid = -1;
00675 bool allframes = false;
00676 
00677 // Read the atom list
00678 int numatms;
00679 Tcl_Obj **data;
00680 if (Tcl_ListObjGetElements(interp, objv[1], &numatms, &data) != TCL_OK) {
00681 Tcl_AppendResult(interp, " measure bond: bad syntax", NULL);
00682 Tcl_AppendResult(interp, " Usage: measure angle {{<atomid1> [<molid1>]} {<atomid2> [<molid2>]} {<atomid3> [<molid3>]}} [molid <default molid>] [frame <frame|all|last> | first <first> last <last>]", NULL);
00683 return TCL_ERROR;
00684 }
00685 
00686 if (numatms != 3) {
00687 Tcl_AppendResult(interp, " measure angle: must specify exactly three atoms in a list", NULL);
00688 return TCL_ERROR;
00689 }
00690 
00691 if (argc > 3) {
00692 for (int i=2; i<argc; i+=2) {
00693 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
00694 if (!strupncmp(argvcur, "molid", CMDLEN)) {
00695 if (Tcl_GetIntFromObj(interp, objv[i+1], &defmolid) != TCL_OK) {
00696 Tcl_AppendResult(interp, " measure angle: bad molid", NULL);
00697 return TCL_ERROR;
00698 }
00699 } else if (!strupncmp(argvcur, "first", CMDLEN)) {
00700 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
00701 Tcl_AppendResult(interp, " measure angle: bad first frame value", NULL);
00702 return TCL_ERROR;
00703 }
00704 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
00705 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
00706 Tcl_AppendResult(interp, " measure angle: bad last frame value", NULL);
00707 return TCL_ERROR;
00708 }
00709 } else if (!strupncmp(argvcur, "frame", CMDLEN)) {
00710 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "all")) {
00711 allframes = true;
00712 } else if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "last")) {
00713 frame=-2;
00714 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &frame) != TCL_OK) {
00715 Tcl_AppendResult(interp, " measure angle: bad frame value", NULL);
00716 return TCL_ERROR;
00717 }
00718 } else {
00719 Tcl_AppendResult(interp, " measure angle: invalid syntax, no such keyword: ", argvcur, NULL);
00720 return TCL_ERROR;
00721 }
00722 }
00723 }
00724 
00725 if ((allframes || frame>=0) && (first>=0 || last>=0)) {
00726 Tcl_AppendResult(interp, "measure angle: Ambiguous syntax: You cannot specify a frame AND a frame range (using first or last).", NULL);
00727 Tcl_AppendResult(interp, "\nUsage:\nmeasure angle {<atomid1> [<molid1>]} {<atomid2> [<molid2>]} {<atomid3> [<molid3>]}} [molid <default molid>] [frame <frame|all|last> | first <first> last <last>]", NULL);
00728 return TCL_ERROR; 
00729 }
00730 
00731 if (allframes) first=0;
00732 
00733 // If no molecule was specified use top as default
00734 if (defmolid<0) defmolid = app->molecule_top();
00735 
00736 // Assign atom IDs and molecule IDs
00737 int i,numelem;
00738 Tcl_Obj **atmmol;
00739 for (i=0; i<numatms; i++) {
00740 if (Tcl_ListObjGetElements(interp, data[i], &numelem, &atmmol) != TCL_OK) {
00741 return TCL_ERROR;
00742 }
00743 
00744 if (!numelem) {
00745 Tcl_AppendResult(interp, " measure angle: empty atom index", NULL);
00746 return TCL_ERROR;
00747 }
00748 
00749 if (Tcl_GetIntFromObj(interp, atmmol[0], atmid+i) != TCL_OK) {
00750 Tcl_AppendResult(interp, " measure angle: bad atom index", NULL);
00751 return TCL_ERROR;
00752 }
00753 
00754 if (numelem==2) {
00755 if (Tcl_GetIntFromObj(interp, atmmol[1], molid+i) != TCL_OK) {
00756 Tcl_AppendResult(interp, " measure angle: bad molid", NULL);
00757 return TCL_ERROR;
00758 }
00759 } else molid[i] = defmolid;
00760 }
00761 
00762 // Compute the value
00763 ResizeArray<float> gValues(1024);
00764 int ret_val;
00765 ret_val = measure_geom(app->moleculeList, molid, atmid, &gValues, 
00766 frame, first, last, defmolid, MEASURE_ANGLE);
00767 if (ret_val<0) {
00768 printf("ERROR\n %s\n", measure_error(ret_val));
00769 Tcl_AppendResult(interp, measure_error(ret_val), NULL);
00770 return TCL_ERROR;
00771 }
00772 
00773 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00774 int numvalues = gValues.num();
00775 for (int count = 0; count < numvalues; count++) {
00776 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(gValues[count]));
00777 }
00778 Tcl_SetObjResult(interp, tcl_result);
00779 
00780 return TCL_OK;
00781 }
00782 
00783 
00784 // vmd_measure_bond {2 atoms as {<atomid> ?<molid>?}} ?molid <molid>? ?frame [f|all|last]? | ?first <first>? ?last <last>?
00785 // Returns the bond length for the specified atoms
00786 static int vmd_measure_bond(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00787 int first=-1, last=-1, frame=-1;
00788 
00789 if (argc<2) {
00790 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "{{<atomid1> [<molid1>]} {<atomid2> [<molid2>]}} [molid <default molid>] [frame <frame|all|last> | first <first> last <last>]]");
00791 return TCL_ERROR;
00792 }
00793 
00794 int molid[2], atmid[2], defmolid = -1;
00795 bool allframes = false;
00796 
00797 // Read the atom list
00798 int numatms=0;
00799 Tcl_Obj **data;
00800 if (Tcl_ListObjGetElements(interp, objv[1], &numatms, &data) != TCL_OK) {
00801 Tcl_AppendResult(interp, " measure bond: bad syntax", NULL);
00802 Tcl_AppendResult(interp, " Usage: measure bond {{<atomid1> [<molid1>]} {<atomid2> [<molid2>]}} [molid <default>] [frame <frame|all|last> | first <first> last <last>]]", NULL);
00803 return TCL_ERROR;
00804 }
00805 
00806 if (numatms != 2) {
00807 Tcl_AppendResult(interp, " measure bond: must specify exactly two atoms in a list", NULL);
00808 return TCL_ERROR;
00809 }
00810 
00811 if (argc > 3) {
00812 for (int i=2; i<argc; i+=2) {
00813 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
00814 if (!strupncmp(argvcur, "molid", CMDLEN)) {
00815 if (Tcl_GetIntFromObj(interp, objv[i+1], &defmolid) != TCL_OK) {
00816 Tcl_AppendResult(interp, " measure bond: bad molid", NULL);
00817 return TCL_ERROR;
00818 }
00819 } else if (!strupncmp(argvcur, "first", CMDLEN)) {
00820 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
00821 Tcl_AppendResult(interp, " measure bond: bad first frame value", NULL);
00822 return TCL_ERROR;
00823 }
00824 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
00825 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
00826 Tcl_AppendResult(interp, " measure bond: bad last frame value", NULL);
00827 return TCL_ERROR;
00828 }
00829 } else if (!strupncmp(argvcur, "frame", CMDLEN)) {
00830 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "all")) {
00831 allframes = true;
00832 } else if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "last")) {
00833 frame=-2;
00834 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &frame) != TCL_OK) {
00835 Tcl_AppendResult(interp, " measure bond: bad frame value", NULL);
00836 return TCL_ERROR;
00837 }
00838 } else {
00839 Tcl_AppendResult(interp, " measure bond: invalid syntax, no such keyword: ", argvcur, NULL);
00840 return TCL_ERROR;
00841 }
00842 }
00843 }
00844 
00845 if ((allframes || frame>=0) && (first>=0 || last>=0)) {
00846 Tcl_AppendResult(interp, "measure bond: Ambiguous syntax: You cannot specify a frame AND a frame range (using first or last).", NULL);
00847 Tcl_AppendResult(interp, "\nUsage:\nmeasure bond {{<atomid1> [<molid1>]} {<atomid2> [<molid2>]}} [molid <default>] [frame <frame|all|last> | first <first> last <last>]", NULL);
00848 return TCL_ERROR; 
00849 }
00850 
00851 if (allframes) first=0;
00852 
00853 // If no molecule was specified use top as default
00854 if (defmolid<0) defmolid = app->molecule_top();
00855 
00856 // Assign atom IDs and molecule IDs
00857 int i, numelem;
00858 Tcl_Obj **atmmol;
00859 for (i=0; i<numatms; i++) {
00860 if (Tcl_ListObjGetElements(interp, data[i], &numelem, &atmmol) != TCL_OK) {
00861 return TCL_ERROR;
00862 }
00863 
00864 if (!numelem) {
00865 Tcl_AppendResult(interp, " measure bond: empty atom index", NULL);
00866 return TCL_ERROR;
00867 }
00868 
00869 if (Tcl_GetIntFromObj(interp, atmmol[0], atmid+i) != TCL_OK) {
00870 Tcl_AppendResult(interp, " measure bond: bad atom index", NULL);
00871 return TCL_ERROR;
00872 }
00873 
00874 if (numelem==2) {
00875 if (Tcl_GetIntFromObj(interp, atmmol[1], molid+i) != TCL_OK) {
00876 Tcl_AppendResult(interp, " measure bond: bad molid", NULL);
00877 return TCL_ERROR;
00878 }
00879 } else molid[i] = defmolid;
00880 }
00881 
00882 // Compute the value
00883 ResizeArray<float> gValues(1024);
00884 int ret_val;
00885 ret_val = measure_geom(app->moleculeList, molid, atmid, &gValues, 
00886 frame, first, last, defmolid, MEASURE_BOND);
00887 if (ret_val < 0) {
00888 printf("ERROR\n %s\n", measure_error(ret_val));
00889 Tcl_AppendResult(interp, measure_error(ret_val), NULL);
00890 return TCL_ERROR;
00891 }
00892 
00893 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00894 int numvalues = gValues.num();
00895 for (int count = 0; count < numvalues; count++) {
00896 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(gValues[count]));
00897 }
00898 Tcl_SetObjResult(interp, tcl_result);
00899 
00900 return TCL_OK;
00901 }
00902 
00903 
00904 // vmd_measure_rmsfperresidue <selection> first <first> last <last> step <step>
00905 // Returns: position variance of the selected atoms over the selected frames, 
00906 // returning one value per residue in the selection.
00907 // Example: measure rmsfperresidue atomselect76 0 20 1
00908 static int vmd_measure_rmsfperresidue(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00909 int first = 0; // start with first frame by default
00910 int last = -1; // finish with last frame by default
00911 int step = 1; // use all frames by default
00912 
00913 if (argc < 2 || argc > 8) {
00914 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<sel> [first <first>] [last <last>] [step <step>]");
00915 return TCL_ERROR;
00916 }
00917 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL)
00918 );
00919 if (!sel) {
00920 Tcl_AppendResult(interp, "measure rmsfperresidue: no atom selection", NULL);
00921 return TCL_ERROR;
00922 }
00923 
00924 int i;
00925 for (i=2; i<argc; i+=2) {
00926 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
00927 if (!strupncmp(argvcur, "first", CMDLEN)) {
00928 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
00929 Tcl_AppendResult(interp, "measure rmsfperresidue: bad first frame value", NULL);
00930 return TCL_ERROR;
00931 }
00932 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
00933 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
00934 Tcl_AppendResult(interp, "measure rmsfperresidue: bad last frame value", NULL);
00935 return TCL_ERROR;
00936 }
00937 } else if (!strupncmp(argvcur, "step", CMDLEN)) {
00938 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK) {
00939 Tcl_AppendResult(interp, "measure rmsfperresidue: bad frame step value", NULL);
00940 return TCL_ERROR;
00941 }
00942 } else {
00943 Tcl_AppendResult(interp, "measure rmsfperresidue: invalid syntax, no such keyword: ", argvcur, NULL);
00944 return TCL_ERROR;
00945 }
00946 }
00947 
00948 float *rmsf = new float[sel->selected];
00949 int ret_val = measure_rmsf_perresidue(sel, app->moleculeList, first, last, step, rmsf);
00950 if (ret_val < 0) {
00951 Tcl_AppendResult(interp, "measure rmsfperresidue: ", measure_error(ret_val), NULL);
00952 delete [] rmsf;
00953 return TCL_ERROR;
00954 }
00955 
00956 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00957 for (i=0; i<ret_val; i++) {
00958 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(rmsf[i]));
00959 }
00960 Tcl_SetObjResult(interp, tcl_result);
00961 delete [] rmsf;
00962 
00963 return TCL_OK;
00964 }
00965 
00966 
00967 // Function: vmd_measure_rmsf <selection> first <first> last <last> step <step>
00968 // Returns: position variance of the selected atoms over the selected frames
00969 // Example: measure rmsf atomselect76 0 20 1
00970 static int vmd_measure_rmsf(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00971 int first = 0; // start with first frame by default
00972 int last = -1; // finish with last frame by default
00973 int step = 1; // use all frames by default
00974 
00975 if (argc < 2 || argc > 8) {
00976 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<sel> [first <first>] [last <last>] [step <step>]");
00977 return TCL_ERROR;
00978 }
00979 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL)
00980 );
00981 if (!sel) {
00982 Tcl_AppendResult(interp, "measure rmsf: no atom selection", NULL);
00983 return TCL_ERROR;
00984 }
00985 
00986 int i;
00987 for (i=2; i<argc; i+=2) {
00988 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
00989 if (!strupncmp(argvcur, "first", CMDLEN)) {
00990 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
00991 Tcl_AppendResult(interp, "measure rmsf: bad first frame value", NULL);
00992 return TCL_ERROR;
00993 }
00994 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
00995 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
00996 Tcl_AppendResult(interp, "measure rmsf: bad last frame value", NULL);
00997 return TCL_ERROR;
00998 }
00999 } else if (!strupncmp(argvcur, "step", CMDLEN)) {
01000 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK) {
01001 Tcl_AppendResult(interp, "measure rmsf: bad frame step value", NULL);
01002 return TCL_ERROR;
01003 }
01004 } else {
01005 Tcl_AppendResult(interp, "measure rmsf: invalid syntax, no such keyword: ", argvcur, NULL);
01006 return TCL_ERROR;
01007 }
01008 }
01009 
01010 float *rmsf = new float[sel->selected];
01011 int ret_val = measure_rmsf(sel, app->moleculeList, first, last, step, rmsf);
01012 if (ret_val < 0) {
01013 Tcl_AppendResult(interp, "measure rmsf: ", measure_error(ret_val), NULL);
01014 delete [] rmsf;
01015 return TCL_ERROR;
01016 }
01017 
01018 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01019 for (i=0; i<sel->selected; i++) {
01020 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(rmsf[i]));
01021 }
01022 Tcl_SetObjResult(interp, tcl_result);
01023 
01024 delete [] rmsf;
01025 return TCL_OK;
01026 }
01027 
01028 
01029 // measure radius of gyration for selected atoms
01030 static int vmd_measure_rgyr(VMDApp *app, int argc, Tcl_Obj *const objv[], Tcl_Interp *interp) {
01031 if (argc != 2 && argc != 4) {
01032 Tcl_WrongNumArgs(interp, 2, objv-1,
01033 (char *)"<selection> [weight <weights>]");
01034 return TCL_ERROR;
01035 }
01036 if (argc == 4 && strcmp(Tcl_GetStringFromObj(objv[2],NULL), "weight")) {
01037 Tcl_SetResult(interp, (char *) "measure rgyr: parameter can only be 'weight'", TCL_STATIC);
01038 return TCL_ERROR;
01039 }
01040 
01041 // get the selection
01042 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
01043 if (!sel) {
01044 Tcl_SetResult(interp, (char *) "measure rgyr: no atom selection", TCL_STATIC);
01045 return TCL_ERROR;
01046 }
01047 
01048 // get the weight
01049 float *weight = new float[sel->selected];
01050 int ret_val = 0;
01051 if (argc == 2) { // only from atom selection, so weight is 1
01052 ret_val = tcl_get_weights(interp, app, sel, NULL, weight);
01053 } else {
01054 ret_val = tcl_get_weights(interp, app, sel, objv[3], weight);
01055 }
01056 if (ret_val < 0) {
01057 Tcl_AppendResult(interp, "measure rgyr: ", measure_error(ret_val), NULL);
01058 delete [] weight;
01059 return TCL_ERROR;
01060 }
01061 
01062 float rgyr;
01063 ret_val = measure_rgyr(sel, app->moleculeList, weight, &rgyr);
01064 delete [] weight;
01065 if (ret_val < 0) {
01066 Tcl_AppendResult(interp, "measure rgyr: ", measure_error(ret_val), NULL);
01067 return TCL_ERROR;
01068 }
01069 Tcl_SetObjResult(interp, Tcl_NewDoubleObj(rgyr));
01070 return TCL_OK;
01071 }
01072 
01073 
01074 // Function: vmd_measure_minmax <selection>
01075 // Returns: the cartesian range of a selection (min/max){x,y,z}
01076 // Example: vmd_measure_minmax atomselect76
01077 // {-5 0 0} {15 10 11.2}
01078 static int vmd_measure_minmax(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01079 const float *radii = NULL;
01080 if (argc != 2 && argc != 3) {
01081 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<selection> [-withradii]");
01082 return TCL_ERROR;
01083 }
01084 if (argc == 3 && strcmp(Tcl_GetStringFromObj(objv[2],NULL), "-withradii")) {
01085 Tcl_SetResult(interp, (char *) "measure minmax: parameter can only be '-withradii'", TCL_STATIC);
01086 return TCL_ERROR;
01087 }
01088 
01089 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
01090 if (!sel) {
01091 Tcl_AppendResult(interp, "measure minmax: no atom selection", NULL);
01092 return TCL_ERROR;
01093 }
01094 
01095 float min_coord[3], max_coord[3];
01096 const float *framepos = sel->coordinates(app->moleculeList);
01097 if (!framepos) return TCL_ERROR;
01098 
01099 // get atom radii if requested
01100 if (argc == 3) {
01101 Molecule *mol = app->moleculeList->mol_from_id(sel->molid());
01102 radii = mol->extraflt.data("radius");
01103 } 
01104 
01105 int ret_val = measure_minmax(sel->num_atoms, sel->on, framepos, radii, 
01106 min_coord, max_coord);
01107 if (ret_val < 0) {
01108 Tcl_AppendResult(interp, "measure minmax: ", measure_error(ret_val), NULL);
01109 return TCL_ERROR;
01110 }
01111 
01112 Tcl_Obj *list1 = Tcl_NewListObj(0, NULL);
01113 Tcl_Obj *list2 = Tcl_NewListObj(0, NULL);
01114 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01115 
01116 Tcl_ListObjAppendElement(interp, list1, Tcl_NewDoubleObj(min_coord[0]));
01117 Tcl_ListObjAppendElement(interp, list1, Tcl_NewDoubleObj(min_coord[1]));
01118 Tcl_ListObjAppendElement(interp, list1, Tcl_NewDoubleObj(min_coord[2]));
01119 
01120 Tcl_ListObjAppendElement(interp, list2, Tcl_NewDoubleObj(max_coord[0]));
01121 Tcl_ListObjAppendElement(interp, list2, Tcl_NewDoubleObj(max_coord[1]));
01122 Tcl_ListObjAppendElement(interp, list2, Tcl_NewDoubleObj(max_coord[2]));
01123 
01124 Tcl_ListObjAppendElement(interp, tcl_result, list1);
01125 Tcl_ListObjAppendElement(interp, tcl_result, list2);
01126 Tcl_SetObjResult(interp, tcl_result);
01127 
01128 return TCL_OK;
01129 }
01130 
01131 
01132 static int vmd_measure_rmsdperresidue(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01133 if (argc !=3 && argc != 5) {
01134 Tcl_WrongNumArgs(interp, 2, objv-1, 
01135 (char *)"<sel1> <sel2> [weight <weights>]");
01136 return TCL_ERROR;
01137 }
01138 
01139 // get the selections
01140 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
01141 AtomSel *sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
01142 if (!sel1 || !sel2) {
01143 Tcl_AppendResult(interp, "measure rmsd: no atom selection", NULL);
01144 return TCL_ERROR;
01145 }
01146 
01147 if (sel1->selected != sel2->selected) {
01148 Tcl_AppendResult(interp, "measure rmsd: selections must have the same number of atoms", NULL);
01149 return TCL_ERROR;
01150 }
01151 if (!sel1->selected) {
01152 Tcl_AppendResult(interp, "measure rmsd: no atoms selected", NULL);
01153 return TCL_ERROR;
01154 }
01155 float *weight = new float[sel1->selected];
01156 int ret_val = 0;
01157 if (argc == 3) {
01158 ret_val = tcl_get_weights(interp, app, sel1, NULL, weight);
01159 } else {
01160 ret_val = tcl_get_weights(interp, app, sel1, objv[4], weight);
01161 }
01162 if (ret_val < 0) {
01163 Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val), NULL);
01164 delete [] weight;
01165 return TCL_ERROR;
01166 }
01167 
01168 // compute the rmsd
01169 float *rmsd = new float[sel1->selected];
01170 int rc = measure_rmsd_perresidue(sel1, sel2, app->moleculeList, 
01171 sel1->selected, weight, rmsd);
01172 if (rc < 0) {
01173 Tcl_AppendResult(interp, "measure rmsd: ", measure_error(rc), NULL);
01174 delete [] weight;
01175 delete [] rmsd;
01176 return TCL_ERROR;
01177 }
01178 delete [] weight;
01179 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01180 int i;
01181 for (i=0; i<rc; i++) {
01182 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(rmsd[i]));
01183 }
01184 delete [] rmsd;
01185 Tcl_SetObjResult(interp, tcl_result);
01186 
01187 return TCL_OK;
01188 }
01189 
01190 // Function: vmd_measure_rmsd <selection1> <selection2> 
01191 // {weight [none|atom value|string array]}
01192 //
01193 // Returns the RMSD between the two selection, taking the weight (if
01194 // any) into account. If number of elements in the weight != num_atoms
01195 // in sel1 then (num in weight = num selected in sel1 = num selected in
01196 // sel2) else (num in weight = total num in sel1 = total num in sel2).
01197 // The weights are taken from the FIRST selection, if needed
01198 // 
01199 // Examples:
01200 // set sel1 [atomselect 0 all]
01201 // set sel2 [atomselect 1 all]
01202 // measure rmsd $sel1 $sel2
01203 // measure rmsd $sel1 $sel2 weight mass
01204 // set sel3 [atomselect 0 "index 3 4 5"]
01205 // set sel4 [atomselect 1 "index 8 5 9"] # gets turned to 5 8 9
01206 // measure rmsd $sel3 $sel4 weight occupancy
01207 //
01208 static int vmd_measure_rmsd(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp)
01209 {
01210 if (argc !=3 && argc != 5) {
01211 Tcl_WrongNumArgs(interp, 2, objv-1, 
01212 (char *)"<sel1> <sel2> [weight <weights>]");
01213 return TCL_ERROR;
01214 }
01215 
01216 // get the selections
01217 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
01218 AtomSel *sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
01219 if (!sel1 || !sel2) {
01220 Tcl_AppendResult(interp, "measure rmsd: no atom selection", NULL);
01221 return TCL_ERROR;
01222 }
01223 
01224 if (sel1->selected != sel2->selected) {
01225 Tcl_AppendResult(interp, "measure rmsd: selections must have the same number of atoms", NULL);
01226 return TCL_ERROR;
01227 }
01228 if (!sel1->selected) {
01229 Tcl_AppendResult(interp, "measure rmsd: no atoms selected", NULL);
01230 return TCL_ERROR;
01231 }
01232 
01233 float *weight = new float[sel1->selected];
01234 int ret_val = 0;
01235 if (argc == 3) {
01236 ret_val = tcl_get_weights(interp, app, sel1, NULL, weight);
01237 } else {
01238 ret_val = tcl_get_weights(interp, app, sel1, objv[4], weight);
01239 }
01240 if (ret_val < 0) {
01241 Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val), NULL);
01242 delete [] weight;
01243 return TCL_ERROR;
01244 }
01245 
01246 // compute the rmsd
01247 float rmsd = 0;
01248 const float *x = sel1->coordinates(app->moleculeList);
01249 const float *y = sel2->coordinates(app->moleculeList);
01250 if (!x || !y) {
01251 delete [] weight;
01252 return TCL_ERROR;
01253 }
01254 
01255 ret_val = measure_rmsd(sel1, sel2, sel1->selected, x, y, weight, &rmsd);
01256 delete [] weight;
01257 if (ret_val < 0) {
01258 Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val), NULL);
01259 return TCL_ERROR;
01260 }
01261 Tcl_SetObjResult(interp, Tcl_NewDoubleObj(rmsd));
01262 
01263 return TCL_OK;
01264 }
01265 
01266 
01268 // measure rmsd_qcp $sel1 $sel2 [weight <weights>]
01269 static int vmd_measure_rmsd_qcp(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01270 if (argc !=3 && argc != 5) {
01271 Tcl_WrongNumArgs(interp, 2, objv-1, 
01272 (char *)"<sel1> <sel2> [weight <weights>]");
01273 return TCL_ERROR;
01274 }
01275 // get the selections
01276 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
01277 AtomSel *sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
01278 if (!sel1 || !sel2) {
01279 Tcl_AppendResult(interp, "measure rmsd: no atom selection", NULL);
01280 return TCL_ERROR;
01281 }
01282 
01283 if (sel1->selected != sel2->selected) {
01284 Tcl_AppendResult(interp, "measure rmsd: selections must have the same number of atoms", NULL);
01285 return TCL_ERROR;
01286 }
01287 if (!sel1->selected) {
01288 Tcl_AppendResult(interp, "measure rmsd: no atoms selected", NULL);
01289 return TCL_ERROR;
01290 }
01291 float *weight = new float[sel1->selected];
01292 {
01293 int ret_val;
01294 if (argc == 3) {
01295 ret_val = tcl_get_weights(interp, app, sel1, NULL, weight);
01296 } else {
01297 ret_val = tcl_get_weights(interp, app, sel1, objv[4], weight);
01298 }
01299 if (ret_val < 0) {
01300 Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val), NULL);
01301 delete [] weight;
01302 return TCL_ERROR;
01303 }
01304 }
01305 
01306 // compute the rmsd
01307 {
01308 float rmsd = 0;
01309 const float *x = sel1->coordinates(app->moleculeList);
01310 const float *y = sel2->coordinates(app->moleculeList);
01311 if (!x || !y) {
01312 delete [] weight;
01313 return TCL_ERROR;
01314 }
01315 int ret_val = measure_rmsd_qcp(app, sel1, sel2, sel1->selected, x, y, weight, &rmsd);
01316 delete [] weight;
01317 if (ret_val < 0) {
01318 Tcl_AppendResult(interp, "measure rmsd: ", measure_error(ret_val), NULL);
01319 return TCL_ERROR;
01320 }
01321 Tcl_SetObjResult(interp, Tcl_NewDoubleObj(rmsd));
01322 }
01323 return TCL_OK;
01324 }
01325 
01326 
01328 // measure rmsdmat_qcp $sel1 [weight <weights>] start s end e step s
01329 static int vmd_measure_rmsdmat_qcp(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01330 int first = 0; // start with first frame by default
01331 int last = -1; // finish with last frame by default
01332 int step = 1; // use all frames by default
01333 
01334 if (argc < 2 || argc > 9) {
01335 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<sel> [weight <weights>] [first <first>] [last <last>] [step <step>]");
01336 return TCL_ERROR;
01337 }
01338 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL)
01339 );
01340 if (!sel) {
01341 Tcl_AppendResult(interp, "measure rmsdmat_qcp: no atom selection", NULL);
01342 return TCL_ERROR;
01343 }
01344 
01345 int i;
01346 for (i=2; i<argc; i+=2) {
01347 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
01348 if (!strupncmp(argvcur, "first", CMDLEN)) {
01349 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
01350 Tcl_AppendResult(interp, "measure rmsdmat_qcp: bad first frame value", NULL);
01351 return TCL_ERROR;
01352 }
01353 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
01354 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
01355 Tcl_AppendResult(interp, "measure rmsdmat_qcp: bad last frame value", NULL);
01356 return TCL_ERROR;
01357 }
01358 } else if (!strupncmp(argvcur, "step", CMDLEN)) {
01359 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK) {
01360 Tcl_AppendResult(interp, "measure rmsdmat_qcp: bad frame step value", NULL);
01361 return TCL_ERROR;
01362 }
01363 } else {
01364 Tcl_AppendResult(interp, "measure rmsdmat_qcp: invalid syntax, no such keyword: ", argvcur, NULL);
01365 return TCL_ERROR;
01366 }
01367 }
01368 
01369 float *weight = NULL;
01370 if (0) {
01371 weight = new float[sel->selected];
01372 int ret_val;
01373 if (argc == 2) {
01374 ret_val = tcl_get_weights(interp, app, sel, NULL, weight);
01375 } else {
01376 ret_val = tcl_get_weights(interp, app, sel, objv[3], weight);
01377 }
01378 
01379 if (ret_val < 0) {
01380 Tcl_AppendResult(interp, "measure rmsdmat_qcp: ", measure_error(ret_val),
01381 NULL);
01382 delete [] weight;
01383 return TCL_ERROR;
01384 }
01385 }
01386 
01387 //If no last frame is given, set it to be the last frame so the framecount is correct.
01388 if (last < 0)
01389 last = app->molecule_numframes(sel->molid()) - 1;
01390 int framecount = (last - first + 1) / step;
01391 float *rmsdmat = (float *) calloc(1, framecount * framecount * sizeof(float));
01392 
01393 // compute the rmsd matrix
01394 int ret_val = measure_rmsdmat_qcp(app, sel, app->moleculeList, 
01395 sel->selected, weight, 
01396 first, last, step, rmsdmat);
01397 delete [] weight;
01398 
01399 if (ret_val < 0) {
01400 Tcl_AppendResult(interp, "measure rmsdmat_qcp: ", measure_error(ret_val), NULL);
01401 return TCL_ERROR;
01402 }
01403 
01404 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01405 long j;
01406 for (j=0; j<framecount; j++) {
01407 Tcl_Obj *rmsdlist = Tcl_NewListObj(0, NULL);
01408 for (i=0; i<framecount; i++) {
01409 Tcl_ListObjAppendElement(interp, rmsdlist, 
01410 Tcl_NewDoubleObj(rmsdmat[j*framecount + i]));
01411 }
01412 Tcl_ListObjAppendElement(interp, tcl_result, rmsdlist);
01413 }
01414 Tcl_SetObjResult(interp, tcl_result);
01415 
01416 free(rmsdmat);
01417 
01418 return TCL_OK;
01419 }
01420 
01421 
01423 // Experimental out-of-core variant 
01424 //
01425 // measure rmsdmat_qcp $sel1 [weight <weights>] start s end e step s
01426 static int vmd_measure_rmsdmat_qcp_ooc(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01427 int first = 0; // start with first frame by default
01428 int last = -1; // finish with last frame by default
01429 int step = 1; // use all frames by default
01430 int i;
01431 float *weight = NULL;
01432 int nfiles = 0;
01433 const char **trjfileset = NULL;
01434 
01435 if (argc < 2) {
01436 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<sel> [weight <weights>] [first <first>] [last <last>] [step <step>] files [file1] [fileN]");
01437 return TCL_ERROR;
01438 }
01439 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL)
01440 );
01441 if (!sel) {
01442 Tcl_AppendResult(interp, "measure rmsdmat_qcp: no atom selection", NULL);
01443 return TCL_ERROR;
01444 }
01445 
01446 for (i=2; i<argc; i+=2) {
01447 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
01448 if (!strupncmp(argvcur, "first", CMDLEN)) {
01449 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
01450 Tcl_AppendResult(interp, "measure rmsdmat_qcp_ooc: bad first frame value", NULL);
01451 return TCL_ERROR;
01452 }
01453 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
01454 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
01455 Tcl_AppendResult(interp, "measure rmsdmat_qcp_ooc: bad last frame value", NULL);
01456 return TCL_ERROR;
01457 }
01458 } else if (!strupncmp(argvcur, "step", CMDLEN)) {
01459 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK) {
01460 Tcl_AppendResult(interp, "measure rmsdmat_qcp_ooc: bad frame step value", NULL);
01461 return TCL_ERROR;
01462 }
01463 } else if (!strupncmp(argvcur, "files", CMDLEN)) {
01464 int list_num;
01465 Tcl_Obj **list_data;
01466 if (Tcl_ListObjGetElements(interp, objv[i+1], &list_num, &list_data) != TCL_OK) {
01467 Tcl_AppendResult(interp, "measure rmsdmat_qcp_ooc: bad trajectory file list", NULL);
01468 return TCL_ERROR;
01469 }
01470 
01471 int f;
01472 nfiles = list_num;
01473 trjfileset = (const char **) calloc(1, nfiles * sizeof(const char *));
01474 for (f=0; f<nfiles; f++) {
01475 trjfileset[f] = Tcl_GetStringFromObj(list_data[f], NULL);
01476 printf("File[%d] '%s'\n", f, trjfileset[f]);
01477 }
01478 } else {
01479 Tcl_AppendResult(interp, "measure rmsdmat_qcp_ooc: invalid syntax, no such keyword: ", argvcur, NULL);
01480 return TCL_ERROR;
01481 }
01482 }
01483 
01484 #if 0
01485 if (0) {
01486 weight = new float[sel->selected];
01487 int ret_val;
01488 if (argc == 2) {
01489 ret_val = tcl_get_weights(interp, app, sel, NULL, weight);
01490 } else {
01491 ret_val = tcl_get_weights(interp, app, sel, objv[3], weight);
01492 }
01493 
01494 if (ret_val < 0) {
01495 Tcl_AppendResult(interp, "measure rmsdmat_qcp: ", measure_error(ret_val),
01496 NULL);
01497 delete [] weight;
01498 return TCL_ERROR;
01499 }
01500 }
01501 #endif
01502 
01503 // compute the rmsd matrix
01504 float *rmsdmat = NULL;
01505 int framecount = 0; // returned after work is done
01506 int ret_val = measure_rmsdmat_qcp_ooc(app, sel, app->moleculeList, 
01507 nfiles, trjfileset,
01508 sel->selected, weight, 
01509 first, last, step, 
01510 framecount, rmsdmat);
01511 delete [] weight;
01512 
01513 if (ret_val < 0) {
01514 Tcl_AppendResult(interp, "measure rmsdmat_qcp_ooc: ", measure_error(ret_val), NULL);
01515 return TCL_ERROR;
01516 }
01517 
01518 #if 0
01519 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01520 long j;
01521 for (j=0; j<framecount; j++) {
01522 Tcl_Obj *rmsdlist = Tcl_NewListObj(0, NULL);
01523 for (i=0; i<framecount; i++) {
01524 Tcl_ListObjAppendElement(interp, rmsdlist, 
01525 Tcl_NewDoubleObj(rmsdmat[j*framecount + i]));
01526 }
01527 Tcl_ListObjAppendElement(interp, tcl_result, rmsdlist);
01528 }
01529 Tcl_SetObjResult(interp, tcl_result);
01530 #endif
01531 
01532 if (trjfileset != NULL)
01533 free(trjfileset);
01534 
01535 if (rmsdmat != NULL)
01536 free(rmsdmat);
01537 
01538 return TCL_OK;
01539 }
01540 
01541 
01542 
01544 // measure fit $sel1 $sel2 [weight <weights>][
01545 static int vmd_measure_fit(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp)
01546 {
01547 AtomSel *sel1, *sel2;
01548 int *order = NULL;
01549 float *weight = NULL;
01550 int rc;
01551 
01552 if (argc != 3 && argc != 5 
01553 && argc != 7) {
01554 Tcl_WrongNumArgs(interp, 2, objv-1, 
01555 (char *)"<sel1> <sel2> [weight <weights>] [order <index list>]");
01556 return TCL_ERROR;
01557 } else if (argc == 5
01558 && strcmp("weight", Tcl_GetStringFromObj(objv[3], NULL)) 
01559 && strcmp("order", Tcl_GetStringFromObj(objv[3], NULL))) {
01560 Tcl_WrongNumArgs(interp, 2, objv-1, 
01561 (char *)"<sel1> <sel2> [weight <weights>] [order <index list>]");
01562 return TCL_ERROR;
01563 } else if (argc == 7 && 
01564 (strcmp("weight", Tcl_GetStringFromObj(objv[3], NULL)) || 
01565 strcmp("order", Tcl_GetStringFromObj(objv[5], NULL)))) {
01566 Tcl_WrongNumArgs(interp, 2, objv-1, 
01567 (char *)"<sel1> <sel2> [weight <weights>] [order <index list>]");
01568 return TCL_ERROR;
01569 }
01570 
01571 // get the selections
01572 sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
01573 sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
01574 
01575 if (!sel1 || !sel2) {
01576 Tcl_AppendResult(interp, "measure fit: no atom selection", NULL);
01577 return TCL_ERROR;
01578 }
01579 
01580 int num = sel1->selected;
01581 if (!num) {
01582 Tcl_AppendResult(interp, "measure fit: no atoms selected", NULL);
01583 return TCL_ERROR;
01584 }
01585 if (num != sel2->selected) {
01586 Tcl_AppendResult(interp, "measure fit: selections must have the same number of atoms", NULL);
01587 return TCL_ERROR;
01588 }
01589 
01590 // get the weights
01591 weight = new float[num]; 
01592 if (argc > 3 && !strcmp("weight", Tcl_GetStringFromObj(objv[3], NULL))) {
01593 // get user weight parameter
01594 rc = tcl_get_weights(interp, app, sel1, objv[4], weight);
01595 } else {
01596 // default to weights of 1.0
01597 rc = tcl_get_weights(interp, app, sel1, NULL, weight);
01598 }
01599 if (rc < 0) {
01600 Tcl_AppendResult(interp, "measure fit: ", measure_error(rc), NULL);
01601 delete [] weight;
01602 return TCL_ERROR;
01603 }
01604 
01605 
01606 int orderparam = 0;
01607 if (argc == 5 && !strcmp("order", Tcl_GetStringFromObj(objv[3], NULL))) {
01608 orderparam = 4;
01609 } else if (argc == 7 && !strcmp("order", Tcl_GetStringFromObj(objv[5], NULL))) {
01610 orderparam = 6;
01611 }
01612 
01613 if (orderparam != 0) {
01614 // get the atom order
01615 order = new int[num];
01616 rc = tcl_get_orders(interp, num, objv[orderparam], order);
01617 if (rc < 0) {
01618 Tcl_AppendResult(interp, "measure fit: ", measure_error(rc), NULL);
01619 delete [] order;
01620 return TCL_ERROR;
01621 }
01622 }
01623 
01624 // compute the transformation matrix
01625 Matrix4 T;
01626 const float *x = sel1->coordinates(app->moleculeList);
01627 const float *y = sel2->coordinates(app->moleculeList);
01628 
01629 int ret_val = MEASURE_ERR_NOMOLECULE;
01630 if (x && y) 
01631 ret_val = measure_fit(sel1, sel2, x, y, weight, order, &T);
01632 
01633 delete [] weight;
01634 
01635 if (order != NULL)
01636 delete [] order;
01637 
01638 if (ret_val < 0) {
01639 Tcl_AppendResult(interp, "measure fit: ", measure_error(ret_val), NULL);
01640 return TCL_ERROR;
01641 }
01642 
01643 // and return the matrix
01644 tcl_append_matrix(interp, T.mat);
01645 return TCL_OK;
01646 }
01647 
01649 // measure inverse 4x4matrix
01650 static int vmd_measure_inverse(int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01651 // make there there is exactly one matrix
01652 if (argc != 2) {
01653 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<matrix>");
01654 return TCL_ERROR;
01655 }
01656 // Get the first matrix
01657 Matrix4 inv;
01658 if (tcl_get_matrix("measure inverse: ",interp,objv[1],inv.mat) != TCL_OK) {
01659 return TCL_ERROR;
01660 }
01661 if (inv.inverse()) {
01662 Tcl_AppendResult(interp, "Singular Matrix; inverse not computed", NULL);
01663 return TCL_ERROR;
01664 }
01665 tcl_append_matrix(interp, inv.mat);
01666 return TCL_OK;
01667 }
01668 
01669 // Find all atoms p in sel1 and q in sel2 within the cutoff. 
01670 static int vmd_measure_contacts(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01671 
01672 // Cutoff, and either one or two atom selections
01673 if (argc != 3 && argc != 4) {
01674 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<cutoff> <sel1> [<sel2>]");
01675 return TCL_ERROR;
01676 }
01677 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
01678 if (!sel1) {
01679 Tcl_AppendResult(interp, "measure contacts: no atom selection", NULL);
01680 return TCL_ERROR;
01681 }
01682 AtomSel *sel2 = NULL;
01683 if (argc == 4) {
01684 sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[3],NULL));
01685 if (!sel2) {
01686 Tcl_AppendResult(interp, "measure contacts: no atom selection", NULL);
01687 return TCL_ERROR;
01688 }
01689 }
01690 if (!sel2) sel2 = sel1;
01691 
01692 double cutoff;
01693 if (Tcl_GetDoubleFromObj(interp, objv[1], &cutoff) != TCL_OK)
01694 return TCL_ERROR;
01695 
01696 const float *pos1 = sel1->coordinates(app->moleculeList);
01697 const float *pos2 = sel2->coordinates(app->moleculeList);
01698 if (!pos1 || !pos2) {
01699 Tcl_AppendResult(interp, "measure contacts: error, molecule contains no coordinates", NULL);
01700 return TCL_ERROR;
01701 }
01702 Molecule *mol1 = app->moleculeList->mol_from_id(sel1->molid());
01703 Molecule *mol2 = app->moleculeList->mol_from_id(sel2->molid());
01704 
01705 GridSearchPair *pairlist = vmd_gridsearch3(pos1, sel1->num_atoms, sel1->on, pos2, sel2->num_atoms, sel2->on, (float) cutoff, -1, (sel1->num_atoms + sel2->num_atoms) * 27L);
01706 GridSearchPair *p, *tmp;
01707 Tcl_Obj *list1 = Tcl_NewListObj(0, NULL);
01708 Tcl_Obj *list2 = Tcl_NewListObj(0, NULL);
01709 for (p=pairlist; p != NULL; p=tmp) {
01710 // throw out pairs that are already bonded
01711 MolAtom *a1 = mol1->atom(p->ind1);
01712 if (mol1 != mol2 || !a1->bonded(p->ind2)) {
01713 Tcl_ListObjAppendElement(interp, list1, Tcl_NewIntObj(p->ind1));
01714 Tcl_ListObjAppendElement(interp, list2, Tcl_NewIntObj(p->ind2));
01715 }
01716 tmp = p->next;
01717 free(p);
01718 }
01719 Tcl_Obj *result = Tcl_NewListObj(0, NULL);
01720 Tcl_ListObjAppendElement(interp, result, list1);
01721 Tcl_ListObjAppendElement(interp, result, list2);
01722 Tcl_SetObjResult(interp, result);
01723 return TCL_OK;
01724 }
01725 
01726 
01727 // measure g(r) for two selections, with delta, rmax, usepbc, first/last/step 
01728 // frame parameters the code will compute the normalized histogram.
01729 static int vmd_measure_gofr(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01730 int i;
01731 // initialize optional arguments to default values
01732 double rmax=10.0;
01733 double delta=0.1;
01734 int usepbc=0;
01735 int selupdate=0;
01736 int first=-1, last=-1, step=1;
01737 int rc;
01738 
01739 // argument error message
01740 const char *argerrmsg = "<sel1> <sel2> [delta <value>] [rmax <value>] [usepbc <bool>] [selupdate <bool>] [first <first>] [last <last>] [step <step>]";
01741 
01742 // Two atom selections and optional keyword/value pairs.
01743 if ((argc < 3) || (argc > 17) || (argc % 2 == 0) ) {
01744 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
01745 return TCL_ERROR;
01746 }
01747 
01748 // check atom selections
01749 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1], NULL));
01750 if (!sel1) {
01751 Tcl_AppendResult(interp, "measure gofr: invalid first atom selection", NULL);
01752 return TCL_ERROR;
01753 }
01754 
01755 AtomSel *sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2], NULL));
01756 if (!sel2) {
01757 Tcl_AppendResult(interp, "measure gofr: invalid second atom selection", NULL);
01758 return TCL_ERROR;
01759 }
01760 
01761 // parse optional arguments
01762 for (i=3; i<argc; i+=2) {
01763 const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
01764 if (i==(argc-1)) {
01765 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
01766 return TCL_ERROR;
01767 }
01768 if (!strcmp(opt, "delta")) {
01769 if (Tcl_GetDoubleFromObj(interp, objv[i+1], &delta) != TCL_OK)
01770 return TCL_ERROR;
01771 if (delta <= 0.0) {
01772 Tcl_AppendResult(interp, "measure gofr: invalid 'delta' value", NULL);
01773 return TCL_ERROR;
01774 }
01775 } else if (!strcmp(opt, "rmax")) {
01776 if (Tcl_GetDoubleFromObj(interp, objv[i+1], &rmax) != TCL_OK)
01777 return TCL_ERROR;
01778 if (rmax <= 0.0) {
01779 Tcl_AppendResult(interp, "measure gofr: invalid 'rmax' value", NULL);
01780 return TCL_ERROR;
01781 }
01782 } else if (!strcmp(opt, "usepbc")) {
01783 if (Tcl_GetBooleanFromObj(interp, objv[i+1], &usepbc) != TCL_OK)
01784 return TCL_ERROR;
01785 } else if (!strcmp(opt, "selupdate")) {
01786 if (Tcl_GetBooleanFromObj(interp, objv[i+1], &selupdate) != TCL_OK)
01787 return TCL_ERROR;
01788 } else if (!strcmp(opt, "first")) {
01789 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK)
01790 return TCL_ERROR;
01791 } else if (!strcmp(opt, "last")) {
01792 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK)
01793 return TCL_ERROR;
01794 } else if (!strcmp(opt, "step")) {
01795 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK)
01796 return TCL_ERROR;
01797 } else { // unknown keyword.
01798 Tcl_AppendResult(interp, "unknown keyword '", opt, "'. usage: measure gofr ", argerrmsg, NULL);
01799 return TCL_ERROR;
01800 }
01801 }
01802 
01803 // allocate and initialize histogram arrays
01804 int count_h = (int)(rmax / delta + 1.0);
01805 double *gofr = new double[count_h];
01806 double *numint = new double[count_h];
01807 double *histog = new double[count_h];
01808 int *framecntr = new int[3];
01809 
01810 // do the gofr calculation
01811 rc = measure_gofr(sel1, sel2, app->moleculeList,
01812 count_h, gofr, numint, histog,
01813 (float) delta,
01814 first, last, step, framecntr,
01815 usepbc, selupdate);
01816 
01817 // XXX: this needs a 'case' structure to provide more meaninful error messages.
01818 if (rc != MEASURE_NOERR) { 
01819 Tcl_AppendResult(interp, "measure gofr: error during g(r) calculation.", NULL);
01820 return TCL_ERROR;
01821 }
01822 
01823 // convert the results of the lowlevel call to tcl lists
01824 // and build a list from them as return value.
01825 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01826 Tcl_Obj *tcl_rlist = Tcl_NewListObj(0, NULL);
01827 Tcl_Obj *tcl_gofr = Tcl_NewListObj(0, NULL);
01828 Tcl_Obj *tcl_numint = Tcl_NewListObj(0, NULL);
01829 Tcl_Obj *tcl_histog = Tcl_NewListObj(0, NULL);
01830 Tcl_Obj *tcl_frames = Tcl_NewListObj(0, NULL);
01831 
01832 // build lists with results ready for plotting
01833 for (i=0; i<count_h; i++) { 
01834 Tcl_ListObjAppendElement(interp, tcl_rlist, Tcl_NewDoubleObj(delta * ((double)i + 0.5)));
01835 Tcl_ListObjAppendElement(interp, tcl_gofr, Tcl_NewDoubleObj(gofr[i]));
01836 Tcl_ListObjAppendElement(interp, tcl_numint, Tcl_NewDoubleObj(numint[i]));
01837 Tcl_ListObjAppendElement(interp, tcl_histog, Tcl_NewDoubleObj(histog[i]));
01838 }
01839 
01840 // build list with number of frames: 
01841 // total, skipped and processed (one entry for each algorithm).
01842 Tcl_ListObjAppendElement(interp, tcl_frames, Tcl_NewIntObj(framecntr[0]));
01843 Tcl_ListObjAppendElement(interp, tcl_frames, Tcl_NewIntObj(framecntr[1]));
01844 Tcl_ListObjAppendElement(interp, tcl_frames, Tcl_NewIntObj(framecntr[2]));
01845 
01846 // build final list-of-lists as return value
01847 Tcl_ListObjAppendElement(interp, tcl_result, tcl_rlist);
01848 Tcl_ListObjAppendElement(interp, tcl_result, tcl_gofr);
01849 Tcl_ListObjAppendElement(interp, tcl_result, tcl_numint);
01850 Tcl_ListObjAppendElement(interp, tcl_result, tcl_histog);
01851 Tcl_ListObjAppendElement(interp, tcl_result, tcl_frames);
01852 Tcl_SetObjResult(interp, tcl_result);
01853 
01854 delete [] gofr;
01855 delete [] numint;
01856 delete [] histog;
01857 delete [] framecntr;
01858 return TCL_OK;
01859 }
01860 
01861 
01862 // measure g(r) for two selections, with delta, rmax, usepbc, first/last/step 
01863 // frame parameters the code will compute the normalized histogram.
01864 static int vmd_measure_rdf(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01865 int i;
01866 // initialize optional arguments to default values
01867 double rmax=10.0;
01868 double delta=0.1;
01869 int usepbc=0;
01870 int selupdate=0;
01871 int first=-1, last=-1, step=1;
01872 int rc;
01873 
01874 // argument error message
01875 const char *argerrmsg = "<sel1> <sel2> [delta <value>] [rmax <value>] [usepbc <bool>] [selupdate <bool>] [first <first>] [last <last>] [step <step>]";
01876 
01877 // Two atom selections and optional keyword/value pairs.
01878 if ((argc < 3) || (argc > 17) || (argc % 2 == 0) ) {
01879 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
01880 return TCL_ERROR;
01881 }
01882 
01883 // check atom selections
01884 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1], NULL));
01885 if (!sel1) {
01886 Tcl_AppendResult(interp, "measure rdf: invalid first atom selection", NULL);
01887 return TCL_ERROR;
01888 }
01889 
01890 AtomSel *sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2], NULL));
01891 if (!sel2) {
01892 Tcl_AppendResult(interp, "measure rdf: invalid second atom selection", NULL);
01893 return TCL_ERROR;
01894 }
01895 
01896 // parse optional arguments
01897 for (i=3; i<argc; i+=2) {
01898 const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
01899 if (i==(argc-1)) {
01900 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
01901 return TCL_ERROR;
01902 }
01903 if (!strcmp(opt, "delta")) {
01904 if (Tcl_GetDoubleFromObj(interp, objv[i+1], &delta) != TCL_OK)
01905 return TCL_ERROR;
01906 if (delta <= 0.0) {
01907 Tcl_AppendResult(interp, "measure rdf: invalid 'delta' value", NULL);
01908 return TCL_ERROR;
01909 }
01910 } else if (!strcmp(opt, "rmax")) {
01911 if (Tcl_GetDoubleFromObj(interp, objv[i+1], &rmax) != TCL_OK)
01912 return TCL_ERROR;
01913 if (rmax <= 0.0) {
01914 Tcl_AppendResult(interp, "measure rdf: invalid 'rmax' value", NULL);
01915 return TCL_ERROR;
01916 }
01917 } else if (!strcmp(opt, "usepbc")) {
01918 if (Tcl_GetBooleanFromObj(interp, objv[i+1], &usepbc) != TCL_OK)
01919 return TCL_ERROR;
01920 } else if (!strcmp(opt, "selupdate")) {
01921 if (Tcl_GetBooleanFromObj(interp, objv[i+1], &selupdate) != TCL_OK)
01922 return TCL_ERROR;
01923 } else if (!strcmp(opt, "first")) {
01924 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK)
01925 return TCL_ERROR;
01926 } else if (!strcmp(opt, "last")) {
01927 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK)
01928 return TCL_ERROR;
01929 } else if (!strcmp(opt, "step")) {
01930 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK)
01931 return TCL_ERROR;
01932 } else { // unknown keyword.
01933 Tcl_AppendResult(interp, "unknown keyword '", opt, "'. usage: measure rdf ", argerrmsg, NULL);
01934 return TCL_ERROR;
01935 }
01936 }
01937 
01938 // allocate and initialize histogram arrays
01939 int count_h = (int)(rmax / delta + 1.0);
01940 double *gofr = new double[count_h];
01941 double *numint = new double[count_h];
01942 double *histog = new double[count_h];
01943 int *framecntr = new int[3];
01944 
01945 // do the gofr calculation
01946 rc = measure_rdf(app, sel1, sel2, app->moleculeList,
01947 count_h, gofr, numint, histog,
01948 (float) delta,
01949 first, last, step, framecntr,
01950 usepbc, selupdate);
01951 
01952 // XXX: this needs a 'case' structure to provide more meaninful error messages.
01953 if (rc != MEASURE_NOERR) { 
01954 Tcl_AppendResult(interp, "measure rdf: error during rdf calculation.", NULL);
01955 return TCL_ERROR;
01956 }
01957 
01958 // convert the results of the lowlevel call to tcl lists
01959 // and build a list from them as return value.
01960 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
01961 Tcl_Obj *tcl_rlist = Tcl_NewListObj(0, NULL);
01962 Tcl_Obj *tcl_gofr = Tcl_NewListObj(0, NULL);
01963 Tcl_Obj *tcl_numint = Tcl_NewListObj(0, NULL);
01964 Tcl_Obj *tcl_histog = Tcl_NewListObj(0, NULL);
01965 Tcl_Obj *tcl_frames = Tcl_NewListObj(0, NULL);
01966 
01967 // build lists with results ready for plotting
01968 for (i=0; i<count_h; i++) { 
01969 Tcl_ListObjAppendElement(interp, tcl_rlist, Tcl_NewDoubleObj(delta * ((double)i + 0.5)));
01970 Tcl_ListObjAppendElement(interp, tcl_gofr, Tcl_NewDoubleObj(gofr[i]));
01971 Tcl_ListObjAppendElement(interp, tcl_numint, Tcl_NewDoubleObj(numint[i]));
01972 Tcl_ListObjAppendElement(interp, tcl_histog, Tcl_NewDoubleObj(histog[i]));
01973 }
01974 
01975 // build list with number of frames: 
01976 // total, skipped and processed (one entry for each algorithm).
01977 Tcl_ListObjAppendElement(interp, tcl_frames, Tcl_NewIntObj(framecntr[0]));
01978 Tcl_ListObjAppendElement(interp, tcl_frames, Tcl_NewIntObj(framecntr[1]));
01979 Tcl_ListObjAppendElement(interp, tcl_frames, Tcl_NewIntObj(framecntr[2]));
01980 
01981 // build final list-of-lists as return value
01982 Tcl_ListObjAppendElement(interp, tcl_result, tcl_rlist);
01983 Tcl_ListObjAppendElement(interp, tcl_result, tcl_gofr);
01984 Tcl_ListObjAppendElement(interp, tcl_result, tcl_numint);
01985 Tcl_ListObjAppendElement(interp, tcl_result, tcl_histog);
01986 Tcl_ListObjAppendElement(interp, tcl_result, tcl_frames);
01987 Tcl_SetObjResult(interp, tcl_result);
01988 
01989 delete [] gofr;
01990 delete [] numint;
01991 delete [] histog;
01992 delete [] framecntr;
01993 return TCL_OK;
01994 }
01995 
01996 
01998 static int vmd_measure_cluster(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
01999 int i,j;
02000 // initialize optional arguments to default values
02001 int algorithm=0; // will be MEASURE_CLUSTER_QT when finished
02002 int likeness=MEASURE_DIST_FITRMSD;
02003 int numcluster=5;
02004 double cutoff=1.0;
02005 float *weights=NULL;
02006 int selupdate=0;
02007 int first=0, last=-1, step=1;
02008 int rc;
02009 
02010 // argument error message
02011 const char *argerrmsg = "<sel> [num <#clusters>] [distfunc <flag>] "
02012 "[cutoff <cutoff>] [first <first>] [last <last>] [step <step>] "
02013 "[selupdate <bool>] [weight <weights>]";
02014 
02015 // Two atom selections and optional keyword/value pairs.
02016 if ((argc < 2) || (argc > 19) || ((argc-1) % 2 == 0) ) {
02017 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
02018 return TCL_ERROR;
02019 }
02020 
02021 // check atom selection
02022 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1], NULL));
02023 if (!sel) {
02024 Tcl_AppendResult(interp, "measure cluster: invalid atom selection", NULL);
02025 return TCL_ERROR;
02026 }
02027 if (!app->molecule_valid_id(sel->molid())) return MEASURE_ERR_NOMOLECULE;
02028 
02029 // parse optional arguments
02030 for (i=2; i<argc; i+=2) {
02031 const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
02032 if (i==(argc-1)) {
02033 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
02034 return TCL_ERROR;
02035 }
02036 if (!strcmp(opt, "num")) {
02037 if (Tcl_GetIntFromObj(interp, objv[i+1], &numcluster) != TCL_OK)
02038 return TCL_ERROR;
02039 if (numcluster < 1) {
02040 Tcl_AppendResult(interp, "measure cluster: invalid 'num' value (cannot be smaller than 1)", NULL);
02041 return TCL_ERROR;
02042 }
02043 } else if (!strcmp(opt, "cutoff")) {
02044 if (Tcl_GetDoubleFromObj(interp, objv[i+1], &cutoff) != TCL_OK)
02045 return TCL_ERROR;
02046 if (cutoff <= 0.0) {
02047 Tcl_AppendResult(interp, "measure cluster: invalid 'cutoff' value (should be larger than 0.0)", NULL);
02048 return TCL_ERROR;
02049 }
02050 } else if (!strcmp(opt, "distfunc")) {
02051 char *argstr = Tcl_GetStringFromObj(objv[i+1], NULL);
02052 if (!strcmp(argstr,"rmsd")) {
02053 likeness = MEASURE_DIST_RMSD;
02054 } else if (!strcmp(argstr,"fitrmsd")) {
02055 likeness = MEASURE_DIST_FITRMSD;
02056 } else if (!strcmp(argstr,"rgyrd")) {
02057 likeness = MEASURE_DIST_RGYRD;
02058 } else {
02059 Tcl_AppendResult(interp, "measure cluster: unknown distance function (supported are 'rmsd', 'rgyrd' and 'fitrmsd')", NULL);
02060 return TCL_ERROR;
02061 }
02062 } else if (!strcmp(opt, "selupdate")) {
02063 if (Tcl_GetBooleanFromObj(interp, objv[i+1], &selupdate) != TCL_OK)
02064 return TCL_ERROR;
02065 } else if (!strcmp(opt, "weight")) {
02066 // NOTE: we cannot use tcl_get_weights here, since we have to 
02067 // get a full (all atoms) list of weights as we may be updating
02068 // the selection in the process. Also we don't support explicit
02069 // lists of weights for now.
02070 const char *weight_string = Tcl_GetStringFromObj(objv[i+1], NULL);
02071 weights = new float[sel->num_atoms];
02072 
02073 if (!weight_string || !strcmp(weight_string, "none")) {
02074 for (j=0; j<sel->num_atoms; j++) weights[j]=1.0f;
02075 } else {
02076 // if a selection string was given, check the symbol table
02077 SymbolTable *atomSelParser = app->atomSelParser; 
02078 // weights must return floating point values, so the symbol must not 
02079 // be a singleword, so macro is NULL.
02080 atomsel_ctxt context(atomSelParser,
02081 app->moleculeList->mol_from_id(sel->molid()), 
02082 sel->which_frame, NULL);
02083 
02084 int fctn = atomSelParser->find_attribute(weight_string);
02085 if (fctn >= 0) {
02086 // the keyword exists, so get the data
02087 // first, check to see that the function returns floats.
02088 // if it doesn't, it makes no sense to use it as a weight
02089 if (atomSelParser->fctns.data(fctn)->returns_a != SymbolTableElement::IS_FLOAT) {
02090 Tcl_AppendResult(interp, "weight attribute must have floating point values", NULL);
02091 delete [] weights;
02092 return MEASURE_ERR_BADWEIGHTPARM; // can't understand weight parameter 
02093 }
02094 
02095 double *tmp_data = new double[sel->num_atoms];
02096 int *all_on = new int[sel->num_atoms];
02097 for (j=0; j<sel->num_atoms; j++) all_on[j]=1;
02098 
02099 atomSelParser->fctns.data(fctn)->keyword_double(
02100 &context, sel->num_atoms, tmp_data, all_on);
02101 
02102 for (j=0; j<sel->num_atoms; j++) weights[j] = (float)tmp_data[j];
02103 
02104 // clean up.
02105 delete [] tmp_data;
02106 delete [] all_on;
02107 }
02108 }
02109 } else if (!strcmp(opt, "first")) {
02110 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK)
02111 return TCL_ERROR;
02112 } else if (!strcmp(opt, "last")) {
02113 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK)
02114 return TCL_ERROR;
02115 } else if (!strcmp(opt, "step")) {
02116 if (Tcl_GetIntFromObj(interp, objv[i+1], &step) != TCL_OK)
02117 return TCL_ERROR;
02118 } else { // unknown keyword.
02119 Tcl_AppendResult(interp, "unknown keyword '", opt, "'. usage: measure cluster ", argerrmsg, NULL);
02120 return TCL_ERROR;
02121 }
02122 }
02123 
02124 // set default for weights if not already defined
02125 if (!weights) {
02126 weights = new float[sel->num_atoms];
02127 for (j=0; j<sel->num_atoms; j++) weights[j]=1.0f;
02128 }
02129 
02130 // Allocate temporary result storage. we add one more cluster 
02131 // slot for collecting unclustered frames in an additional "cluster".
02132 // NOTE: the individual cluster lists are going to
02133 // allocated in the ancilliary code.
02134 int *clustersize = new int [numcluster+1];
02135 int **clusterlist = new int *[numcluster+1];
02136 
02137 // do the cluster analysis
02138 rc = measure_cluster(sel, app->moleculeList, numcluster, algorithm, likeness, cutoff, 
02139 clustersize, clusterlist, first, last, step, selupdate, weights);
02140 
02141 if (weights) delete [] weights;
02142 
02143 // XXX: this needs a 'case' structure to provide more meaninful error messages.
02144 if (rc != MEASURE_NOERR) { 
02145 Tcl_AppendResult(interp, "measure cluster: error during cluster analysis calculation.", NULL);
02146 return TCL_ERROR;
02147 }
02148 
02149 // convert the results of the lowlevel call to tcl lists
02150 // and build a list from them as return value.
02151 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
02152 
02153 for (i=0; i <= numcluster; ++i) { 
02154 int j;
02155 Tcl_Obj *tcl_clist = Tcl_NewListObj(0, NULL);
02156 for (j=0; j < clustersize[i]; ++j) {
02157 Tcl_ListObjAppendElement(interp, tcl_clist, Tcl_NewIntObj(clusterlist[i][j]));
02158 }
02159 Tcl_ListObjAppendElement(interp, tcl_result, tcl_clist);
02160 }
02161 Tcl_SetObjResult(interp, tcl_result);
02162 
02163 // free temporary result storage
02164 for (i=0; i <= numcluster; ++i)
02165 delete[] clusterlist[i];
02166 
02167 delete[] clusterlist;
02168 delete[] clustersize;
02169 
02170 return TCL_OK;
02171 }
02172 
02174 static int vmd_measure_clustsize(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02175 int i;
02176 // initialize optional arguments to default values
02177 double cutoff=3.0; // arbitrary. took hbond cutoff.
02178 char *storesize=NULL;
02179 char *storenum=NULL;
02180 int usepbc=0;
02181 int minsize=2;
02182 int numshared=1;
02183 int rc;
02184 
02185 // argument error message
02186 const char *argerrmsg = "<sel> [cutoff <float>] [minsize <num>] [numshared <num>] "
02187 "[usepbc <bool>] [storesize <fieldname>] [storenum <fieldname>]";
02188 
02189 // Two atom selections and optional keyword/value pairs.
02190 if ((argc < 2) || (argc > 13) || ((argc-1) % 2 == 0) ) {
02191 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
02192 return TCL_ERROR;
02193 }
02194 
02195 // check atom selection
02196 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1], NULL));
02197 if (!sel) {
02198 Tcl_AppendResult(interp, "measure clustsize: invalid atom selection", NULL);
02199 return TCL_ERROR;
02200 }
02201 if (!app->molecule_valid_id(sel->molid())) return MEASURE_ERR_NOMOLECULE;
02202 
02203 // parse optional arguments
02204 for (i=2; i<argc; i+=2) {
02205 const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
02206 if (i==(argc-1)) {
02207 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)argerrmsg);
02208 return TCL_ERROR;
02209 } else if (!strcmp(opt, "cutoff")) {
02210 if (Tcl_GetDoubleFromObj(interp, objv[i+1], &cutoff) != TCL_OK)
02211 return TCL_ERROR;
02212 if (cutoff <= 0.0) {
02213 Tcl_AppendResult(interp, "measure clustsize: invalid 'cutoff' value", NULL);
02214 return TCL_ERROR;
02215 }
02216 } else if (!strcmp(opt, "minsize")) {
02217 if (Tcl_GetIntFromObj(interp, objv[i+1], &minsize) != TCL_OK)
02218 return TCL_ERROR;
02219 if (minsize < 2) {
02220 Tcl_AppendResult(interp, "measure clustsize: invalid 'minsize' value", NULL);
02221 return TCL_ERROR;
02222 }
02223 } else if (!strcmp(opt, "numshared")) {
02224 if (Tcl_GetIntFromObj(interp, objv[i+1], &numshared) != TCL_OK)
02225 return TCL_ERROR;
02226 if (numshared < 0) {
02227 Tcl_AppendResult(interp, "measure clustsize: invalid 'numshared' value", NULL);
02228 return TCL_ERROR;
02229 }
02230 } else if (!strcmp(opt, "usepbc")) {
02231 if (Tcl_GetBooleanFromObj(interp, objv[i+1], &usepbc) != TCL_OK)
02232 return TCL_ERROR;
02233 } else if (!strcmp(opt, "storenum")) {
02234 storenum = Tcl_GetStringFromObj(objv[i+1], NULL);
02235 } else if (!strcmp(opt, "storesize")) {
02236 storesize = Tcl_GetStringFromObj(objv[i+1], NULL);
02237 } else { // unknown keyword.
02238 Tcl_AppendResult(interp, "unknown keyword '", opt, "'. usage: measure clustsize ", argerrmsg, NULL);
02239 return TCL_ERROR;
02240 }
02241 }
02242 
02243 if (usepbc) { 
02244 Tcl_AppendResult(interp, "measure clustsize: does not support periodic boundaries yet.", NULL);
02245 return TCL_ERROR;
02246 }
02247 
02248 // allocate temporary result storage
02249 // NOTE: the individual cluster lists are going to
02250 // allocated in the ancilliary code.
02251 int num_selected=sel->selected;
02252 int *clustersize = new int[num_selected];
02253 int *clusternum= new int [num_selected];
02254 int *clusteridx= new int [num_selected];
02255 for (i=0; i < num_selected; i++) {
02256 clustersize[i] = 0;
02257 clusternum[i] = -1;
02258 clusteridx[i] = -1;
02259 }
02260 
02261 // do the cluster analysis
02262 rc = measure_clustsize(sel, app->moleculeList, cutoff, 
02263 clustersize, clusternum, clusteridx,
02264 minsize, numshared, usepbc);
02265 
02266 // XXX: this needs a 'case' structure to provide more meaninful error messages.
02267 if (rc != MEASURE_NOERR) { 
02268 Tcl_AppendResult(interp, "measure clustsize: error during cluster size analysis calculation.", NULL);
02269 return TCL_ERROR;
02270 }
02271 
02272 
02273 if (storenum || storesize) {
02274 // field names were given to store the results. check the keywords and so on.
02275 SymbolTable *atomSelParser = app->atomSelParser; 
02276 atomsel_ctxt context(atomSelParser, 
02277 app->moleculeList->mol_from_id(sel->molid()), 
02278 sel->which_frame, NULL);
02279 
02280 // the keyword exists, set the data
02281 if (storenum) {
02282 int fctn = atomSelParser->find_attribute(storenum);
02283 if (fctn >= 0) {
02284 if (atomSelParser->fctns.data(fctn)->returns_a == SymbolTableElement::IS_FLOAT) {
02285 double *tmp_data = new double[sel->num_atoms];
02286 int j=0;
02287 for (int i=sel->firstsel; i<=sel->lastsel; i++) {
02288 if (sel->on[i])
02289 tmp_data[i] = (double) clusternum[j++];
02290 }
02291 atomSelParser->fctns.data(fctn)->set_keyword_double(&context, 
02292 sel->num_atoms,
02293 tmp_data, sel->on);
02294 delete[] tmp_data;
02295 
02296 } else if (atomSelParser->fctns.data(fctn)->returns_a == SymbolTableElement::IS_INT) {
02297 int *tmp_data = new int[sel->num_atoms];
02298 int j=0;
02299 for (int i=sel->firstsel; i<=sel->lastsel; i++) {
02300 if (sel->on[i])
02301 tmp_data[i] = clusternum[j++];
02302 }
02303 atomSelParser->fctns.data(fctn)->set_keyword_int(&context, 
02304 sel->num_atoms,
02305 tmp_data, sel->on);
02306 delete[] tmp_data;
02307 } else {
02308 Tcl_AppendResult(interp, "measure clustsize: storenum field must accept numbers", NULL);
02309 return TCL_ERROR;
02310 }
02311 } else {
02312 Tcl_AppendResult(interp, "measure clustsize: invalid field name for storenum", NULL);
02313 return TCL_ERROR;
02314 }
02315 }
02316 
02317 // the keyword exists, set the data
02318 if (storesize) {
02319 int fctn = atomSelParser->find_attribute(storesize);
02320 if (fctn >= 0) {
02321 if (atomSelParser->fctns.data(fctn)->returns_a == SymbolTableElement::IS_FLOAT) {
02322 double *tmp_data = new double[sel->num_atoms];
02323 int j=0;
02324 for (int i=sel->firstsel; i<=sel->lastsel; i++) {
02325 if (sel->on[i])
02326 tmp_data[i] = (double) clustersize[j++];
02327 }
02328 atomSelParser->fctns.data(fctn)->set_keyword_double(&context, 
02329 sel->num_atoms,
02330 tmp_data, sel->on);
02331 delete[] tmp_data;
02332 
02333 } else if (atomSelParser->fctns.data(fctn)->returns_a == SymbolTableElement::IS_INT) {
02334 int *tmp_data = new int[sel->num_atoms];
02335 int j=0;
02336 for (int i=sel->firstsel; i<=sel->lastsel; i++) {
02337 if (sel->on[i])
02338 tmp_data[i] = clustersize[j++];
02339 }
02340 atomSelParser->fctns.data(fctn)->set_keyword_int(&context, 
02341 sel->num_atoms,
02342 tmp_data, sel->on);
02343 delete[] tmp_data;
02344 } else {
02345 Tcl_AppendResult(interp, "measure clustsize: storenum field must accept numbers", NULL);
02346 return TCL_ERROR;
02347 }
02348 } else {
02349 Tcl_AppendResult(interp, "measure clustsize: invalid field name for storesize", NULL);
02350 return TCL_ERROR;
02351 }
02352 }
02353 } else {
02354 // convert the results of the lowlevel call to tcl lists
02355 // and build a list from them as return value.
02356 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
02357 
02358 Tcl_Obj *tcl_ilist = Tcl_NewListObj(0, NULL);
02359 Tcl_Obj *tcl_clist = Tcl_NewListObj(0, NULL);
02360 Tcl_Obj *tcl_nlist = Tcl_NewListObj(0, NULL);
02361 for (i=0; i<num_selected; i++) { 
02362 Tcl_ListObjAppendElement(interp, tcl_ilist, Tcl_NewIntObj(clusteridx[i]));
02363 Tcl_ListObjAppendElement(interp, tcl_clist, Tcl_NewIntObj(clustersize[i]));
02364 Tcl_ListObjAppendElement(interp, tcl_nlist, Tcl_NewIntObj(clusternum[i]));
02365 }
02366 Tcl_ListObjAppendElement(interp, tcl_result, tcl_ilist);
02367 Tcl_ListObjAppendElement(interp, tcl_result, tcl_clist);
02368 Tcl_ListObjAppendElement(interp, tcl_result, tcl_nlist);
02369 Tcl_SetObjResult(interp, tcl_result);
02370 }
02371 
02372 delete[] clustersize;
02373 delete[] clusternum;
02374 delete[] clusteridx;
02375 
02376 return TCL_OK;
02377 }
02378 
02379 static int vmd_measure_hbonds(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02380 
02381 // Cutoff, angle, and either one or two atom selections
02382 if (argc != 4 && argc != 5) {
02383 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<cutoff> <angle> <selection1> [<selection2>]");
02384 return TCL_ERROR;
02385 }
02386 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[3],NULL));
02387 if (!sel1) {
02388 Tcl_AppendResult(interp, "measure hbonds: invalid first atom selection", NULL);
02389 return TCL_ERROR;
02390 }
02391 
02392 AtomSel *sel2 = NULL;
02393 if (argc == 5) {
02394 sel2 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[4],NULL));
02395 if (!sel2) {
02396 Tcl_AppendResult(interp, "measure hbonds: invalid second atom selection", NULL);
02397 return TCL_ERROR;
02398 }
02399 }
02400 if (sel2 && sel2->molid() != sel1->molid()) {
02401 Tcl_AppendResult(interp, "measure hbonds: error, atom selections must come from same molecule.", NULL);
02402 return TCL_ERROR;
02403 }
02404 double cutoff;
02405 if (Tcl_GetDoubleFromObj(interp, objv[1], &cutoff) != TCL_OK) 
02406 return TCL_ERROR;
02407 
02408 double maxangle;
02409 if (Tcl_GetDoubleFromObj(interp, objv[2], &maxangle) != TCL_OK) 
02410 return TCL_ERROR;
02411 
02412 const float *pos = sel1->coordinates(app->moleculeList);
02413 if (!pos) {
02414 Tcl_AppendResult(interp, "measure bondsearch: error, molecule contains no coordinates", NULL);
02415 return TCL_ERROR;
02416 }
02417 
02418 // XXX the actual code for measuring hbonds doesn't belong here, it should
02419 // be moved into Measure.[Ch] where it really belongs. This file
02420 // only implements the Tcl interface, and should not be doing the
02421 // hard core math, particularly if we want to expose the same
02422 // feature via other scripting interfaces. Also, having a single
02423 // implementation avoids having different Tcl/Python bugs in the 
02424 // long-term. Too late to do anything about this now, but should be
02425 // addressed for the next major version when time allows.
02426 
02427 // XXX This code is close, but not identical to the HBonds code in 
02428 // DrawMolItem. Is there any good reason they aren't identical?
02429 // This version does a few extra tests that the other does not.
02430 
02431 Molecule *mol = app->moleculeList->mol_from_id(sel1->molid());
02432 
02433 int *donlist, *hydlist, *acclist;
02434 int maxsize = 2 * sel1->num_atoms; //This heuristic is based on ice, where there are < 2 hydrogen bonds per atom if hydrogens are in the selection, and exactly 2 if hydrogens are not considered.
02435 donlist = new int[maxsize];
02436 hydlist = new int[maxsize];
02437 acclist = new int[maxsize];
02438 int rc = measure_hbonds(mol, sel1, sel2, cutoff, maxangle, donlist, hydlist, acclist, maxsize);
02439 if (rc > maxsize) {
02440 delete [] donlist;
02441 delete [] hydlist;
02442 delete [] acclist;
02443 maxsize = rc;
02444 donlist = new int[maxsize];
02445 hydlist = new int[maxsize];
02446 acclist = new int[maxsize];
02447 rc = measure_hbonds(mol, sel1, sel2, cutoff, maxangle, donlist, hydlist, acclist, maxsize);
02448 }
02449 if (rc < 0) {
02450 Tcl_AppendResult(interp, "measure hbonds: internal error to measure_hbonds", NULL);
02451 return TCL_ERROR;
02452 }
02453 
02454 Tcl_Obj *newdonlist = Tcl_NewListObj(0, NULL);
02455 Tcl_Obj *newhydlist = Tcl_NewListObj(0, NULL);
02456 Tcl_Obj *newacclist = Tcl_NewListObj(0, NULL);
02457 for (int k = 0; k < rc; k++) {
02458 Tcl_ListObjAppendElement(interp, newdonlist, Tcl_NewIntObj(donlist[k]));
02459 Tcl_ListObjAppendElement(interp, newhydlist, Tcl_NewIntObj(hydlist[k]));
02460 Tcl_ListObjAppendElement(interp, newacclist, Tcl_NewIntObj(acclist[k]));
02461 }
02462 delete [] donlist;
02463 delete [] hydlist;
02464 delete [] acclist;
02465 Tcl_Obj *result = Tcl_NewListObj(0, NULL);
02466 Tcl_ListObjAppendElement(interp, result, newdonlist);
02467 Tcl_ListObjAppendElement(interp, result, newacclist);
02468 Tcl_ListObjAppendElement(interp, result, newhydlist);
02469 Tcl_SetObjResult(interp, result);
02470 return TCL_OK;
02471 }
02472 
02473 static int vmd_measure_sasaperresidue(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02474 
02475 int i;
02476 // srad and one atom selection, plus additional options
02477 if (argc < 3) {
02478 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<srad> <sel> [-restrict <restrictedsel>] [-samples <numsamples>]");
02479 return TCL_ERROR;
02480 }
02481 // parse options
02482 Tcl_Obj *ptsvar = NULL;
02483 AtomSel *restrictsel = NULL;
02484 int nsamples = -1, rescount;
02485 int *sampleptr = NULL;
02486 for (i=3; i<argc-1; i+=2) {
02487 const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
02488 if (!strcmp(opt, "-points")) {
02489 ptsvar = objv[i+1];
02490 } else if (!strcmp(opt, "-restrict")) {
02491 restrictsel = tcl_commands_get_sel(interp, 
02492 Tcl_GetStringFromObj(objv[i+1], NULL));
02493 if (!restrictsel) {
02494 Tcl_AppendResult(interp, "measure sasa: invalid restrict atom selection", NULL);
02495 return TCL_ERROR;
02496 }
02497 } else if (!strcmp(opt, "-samples")) {
02498 if (Tcl_GetIntFromObj(interp, objv[i+1], &nsamples) != TCL_OK)
02499 return TCL_ERROR;
02500 sampleptr = &nsamples;
02501 } else {
02502 Tcl_AppendResult(interp, "measure sasa per residue: unknown option '", opt, "'", 
02503 NULL);
02504 return TCL_ERROR;
02505 }
02506 }
02507 
02508 double srad;
02509 if (Tcl_GetDoubleFromObj(interp, objv[1], &srad) != TCL_OK) 
02510 return TCL_ERROR;
02511 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
02512 if (!sel1) {
02513 Tcl_AppendResult(interp, "measure sasa per residue: invalid first atom selection", NULL);
02514 return TCL_ERROR;
02515 }
02516 
02517 const float *pos = sel1->coordinates(app->moleculeList);
02518 if (!pos) {
02519 Tcl_AppendResult(interp, "measure sasa per residue: error, molecule contains no coordinates", NULL);
02520 return TCL_ERROR;
02521 }
02522 Molecule *mol = app->moleculeList->mol_from_id(sel1->molid());
02523 const float *radius = mol->extraflt.data("radius");
02524 
02525 ResizeArray<float> sasapts;
02526 float *sasavals;
02527 sasavals = new float[sel1->selected];
02528 int rc = measure_sasa_perresidue(sel1, pos, radius, (float) srad, sasavals, &sasapts, 
02529 restrictsel, sampleptr,&rescount,app->moleculeList);
02530 if (rc < 0) {
02531 Tcl_AppendResult(interp, "measure: sasa per residue : ", measure_error(rc), NULL);
02532 delete [] sasavals;
02533 return TCL_ERROR;
02534 }
02535 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
02536 for (int i=0; i<rescount; i++) {
02537 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(sasavals[i]));
02538 }
02539 Tcl_SetObjResult(interp, tcl_result);
02540 if (ptsvar) {
02541 // construct list from sasapts
02542 Tcl_Obj *listobj = Tcl_NewListObj(0, NULL);
02543 i=0;
02544 while (i<sasapts.num()) {
02545 Tcl_Obj *elem = Tcl_NewListObj(0, NULL);
02546 Tcl_ListObjAppendElement(interp, elem, Tcl_NewDoubleObj(sasapts[i++]));
02547 Tcl_ListObjAppendElement(interp, elem, Tcl_NewDoubleObj(sasapts[i++]));
02548 Tcl_ListObjAppendElement(interp, elem, Tcl_NewDoubleObj(sasapts[i++]));
02549 Tcl_ListObjAppendElement(interp, listobj, elem);
02550 }
02551 Tcl_ObjSetVar2(interp, ptsvar, NULL, listobj, 0);
02552 }
02553 delete [] sasavals;
02554 return TCL_OK;
02555 }
02556 
02557 
02558 static int vmd_measure_sasa(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02559 
02560 int i;
02561 // srad and one atom selection, plus additional options
02562 if (argc < 3) {
02563 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<srad> <sel> [-points <varname>] [-restrict <restrictedsel>] [-samples <numsamples>]");
02564 return TCL_ERROR;
02565 }
02566 // parse options
02567 Tcl_Obj *ptsvar = NULL;
02568 AtomSel *restrictsel = NULL;
02569 int nsamples = -1;
02570 int *sampleptr = NULL;
02571 for (i=3; i<argc-1; i+=2) {
02572 const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
02573 if (!strcmp(opt, "-points")) {
02574 ptsvar = objv[i+1];
02575 } else if (!strcmp(opt, "-restrict")) {
02576 restrictsel = tcl_commands_get_sel(interp, 
02577 Tcl_GetStringFromObj(objv[i+1], NULL));
02578 if (!restrictsel) {
02579 Tcl_AppendResult(interp, "measure sasa: invalid restrict atom selection", NULL);
02580 return TCL_ERROR;
02581 }
02582 } else if (!strcmp(opt, "-samples")) {
02583 if (Tcl_GetIntFromObj(interp, objv[i+1], &nsamples) != TCL_OK)
02584 return TCL_ERROR;
02585 sampleptr = &nsamples;
02586 } else {
02587 Tcl_AppendResult(interp, "measure sasa: unknown option '", opt, "'", 
02588 NULL);
02589 return TCL_ERROR;
02590 }
02591 }
02592 
02593 double srad;
02594 if (Tcl_GetDoubleFromObj(interp, objv[1], &srad) != TCL_OK) 
02595 return TCL_ERROR;
02596 AtomSel *sel1 = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[2],NULL));
02597 if (!sel1) {
02598 Tcl_AppendResult(interp, "measure sasa: invalid first atom selection", NULL);
02599 return TCL_ERROR;
02600 }
02601 
02602 const float *pos = sel1->coordinates(app->moleculeList);
02603 if (!pos) {
02604 Tcl_AppendResult(interp, "measure sasa: error, molecule contains no coordinates", NULL);
02605 return TCL_ERROR;
02606 }
02607 Molecule *mol = app->moleculeList->mol_from_id(sel1->molid());
02608 const float *radius = mol->extraflt.data("radius");
02609 
02610 ResizeArray<float> sasapts;
02611 float sasa = 0;
02612 int rc = measure_sasa(sel1, pos, radius, (float) srad, &sasa, &sasapts, 
02613 restrictsel, sampleptr);
02614 if (rc < 0) {
02615 Tcl_AppendResult(interp, "measure: sasa: ", measure_error(rc), NULL);
02616 return TCL_ERROR;
02617 }
02618 Tcl_SetObjResult(interp, Tcl_NewDoubleObj(sasa));
02619 if (ptsvar) {
02620 // construct list from sasapts
02621 Tcl_Obj *listobj = Tcl_NewListObj(0, NULL);
02622 i=0;
02623 while (i<sasapts.num()) {
02624 Tcl_Obj *elem = Tcl_NewListObj(0, NULL);
02625 Tcl_ListObjAppendElement(interp, elem, Tcl_NewDoubleObj(sasapts[i++]));
02626 Tcl_ListObjAppendElement(interp, elem, Tcl_NewDoubleObj(sasapts[i++]));
02627 Tcl_ListObjAppendElement(interp, elem, Tcl_NewDoubleObj(sasapts[i++]));
02628 Tcl_ListObjAppendElement(interp, listobj, elem);
02629 }
02630 Tcl_ObjSetVar2(interp, ptsvar, NULL, listobj, 0);
02631 }
02632 return TCL_OK;
02633 }
02634 
02635 
02636 #if 1
02637 
02638 static int vmd_measure_sasalist(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02639 
02640 int i;
02641 // srad and one atom selection, plus additional options
02642 if (argc < 3) {
02643 Tcl_WrongNumArgs(interp, 2, objv-1, (char *)"<srad> <sel list> [-samples <numsamples>]");
02644 return TCL_ERROR;
02645 }
02646 
02647 // parse options
02648 int nsamples = -1;
02649 int *sampleptr = NULL;
02650 for (i=3; i<argc-1; i+=2) {
02651 const char *opt = Tcl_GetStringFromObj(objv[i], NULL);
02652 
02653 if (!strcmp(opt, "-samples")) {
02654 if (Tcl_GetIntFromObj(interp, objv[i+1], &nsamples) != TCL_OK)
02655 return TCL_ERROR;
02656 sampleptr = &nsamples;
02657 } else {
02658 Tcl_AppendResult(interp, "measure sasa: unknown option '", opt, "'", NULL);
02659 return TCL_ERROR;
02660 }
02661 }
02662 
02663 double srad;
02664 if (Tcl_GetDoubleFromObj(interp, objv[1], &srad) != TCL_OK) 
02665 return TCL_ERROR;
02666 
02667 int numsels;
02668 Tcl_Obj **sel_list;
02669 if (Tcl_ListObjGetElements(interp, objv[2], &numsels, &sel_list) != TCL_OK) {
02670 Tcl_AppendResult(interp, "measure sasalist: bad selection list", NULL);
02671 return TCL_ERROR;
02672 }
02673 
02674 #if 0
02675 printf("measure sasalist: numsels %d\n", numsels);
02676 #endif
02677 
02678 AtomSel **asels = (AtomSel **) calloc(1, numsels * sizeof(AtomSel *));
02679 float *sasalist = (float *) calloc(1, numsels * sizeof(float));
02680 
02681 int s;
02682 for (s=0; s<numsels; s++) {
02683 asels[s] = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(sel_list[s], NULL));
02684 if (!asels[s]) {
02685 printf("measure sasalist: invalid selection %d\n", s);
02686 Tcl_AppendResult(interp, "measure sasalist: invalid atom selection list element", NULL);
02687 return TCL_ERROR;
02688 }
02689 }
02690 
02691 int rc = measure_sasalist(app->moleculeList, (const AtomSel **) asels, 
02692 numsels, (float) srad, sasalist, sampleptr);
02693 free(asels);
02694 if (rc < 0) {
02695 Tcl_AppendResult(interp, "measure: sasalist: ", measure_error(rc), NULL);
02696 return TCL_ERROR;
02697 }
02698 
02699 // construct list from sasa values
02700 Tcl_Obj *listobj = Tcl_NewListObj(0, NULL);
02701 for (i=0; i<numsels; i++) {
02702 Tcl_ListObjAppendElement(interp, listobj, Tcl_NewDoubleObj(sasalist[i]));
02703 }
02704 Tcl_SetObjResult(interp, listobj);
02705 free(sasalist);
02706 
02707 return TCL_OK;
02708 }
02709 
02710 #endif
02711 
02712 
02713 // Function: vmd_measure_energy 
02714 // Returns: the energy for the specified bond/angle/dihed/imprp/vdw/elec
02715 // FIXME -- usage doesn't match user guide
02716 static int vmd_measure_energy(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02717 
02718 if (argc<3) {
02719 Tcl_WrongNumArgs(interp, 2, objv-1,
02720 (char *) "bond|angle|dihed|imprp|vdw|elec {{<atomid1> ?<molid1>?} {<atomid2> ?<molid2>?}} ?molid <default molid>? [?frame <frame|all>? | ?first <first>? ?last <last>?]");
02721 return TCL_ERROR;
02722 }
02723 
02724 int geomtype, reqatoms;
02725 char *geomname = Tcl_GetStringFromObj(objv[1],NULL);
02726 if (!strncmp(geomname, "bond", 4)) {
02727 reqatoms = 2; geomtype = MEASURE_BOND;
02728 } else if (!strncmp(geomname, "angl", 4)) {
02729 reqatoms = 3; geomtype = MEASURE_ANGLE;
02730 } else if (!strncmp(geomname, "dihe", 4)) {
02731 reqatoms = 4; geomtype = MEASURE_DIHED;
02732 } else if (!strncmp(geomname, "impr", 4)) {
02733 reqatoms = 4; geomtype = MEASURE_IMPRP;
02734 } else if (!strncmp(geomname, "vdw", 3)) {
02735 reqatoms = 2; geomtype = MEASURE_VDW;
02736 } else if (!strncmp(geomname, "elec", 4)) {
02737 reqatoms = 2; geomtype = MEASURE_ELECT;
02738 } else {
02739 Tcl_AppendResult(interp, " measure energy: bad syntax (must specify bond|angle|dihed|imprp|vdw|elec)", NULL);
02740 return TCL_ERROR;
02741 }
02742 
02743 int molid[4];
02744 int atmid[4];
02745 int defmolid = -1;
02746 bool allframes = false;
02747 char errstring[200];
02748 
02749 // Read the atom list
02750 int numatms;
02751 Tcl_Obj **data;
02752 if (Tcl_ListObjGetElements(interp, objv[2], &numatms, &data) != TCL_OK) {
02753 Tcl_AppendResult(interp, " measure energy: bad syntax", NULL);
02754 Tcl_AppendResult(interp, " Usage: measure energy bond|angle|dihed|imprp|vdw|elec {{<atomid1> ?<molid1>?} {<atomid2> ?<molid2>?} ...} ?molid <default molid>? [?frame <frame|all>? | ?first <first>? ?last <last>?]", NULL);
02755 return TCL_ERROR;
02756 }
02757 
02758 if (numatms!=reqatoms) {
02759 sprintf(errstring, " measure energy %s: must specify exactly %i atoms in list", geomname, reqatoms);
02760 Tcl_AppendResult(interp, errstring, NULL);
02761 return TCL_ERROR;
02762 }
02763 
02764 int first=-1, last=-1, frame=-1;
02765 double params[6];
02766 memset(params, 0, 6L*sizeof(double));
02767 
02768 if (argc>4) {
02769 int i;
02770 for (i=3; i<argc-1; i+=2) {
02771 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
02772 if (!strupncmp(argvcur, "molid", CMDLEN)) {
02773 if (Tcl_GetIntFromObj(interp, objv[i+1], &defmolid) != TCL_OK) {
02774 Tcl_AppendResult(interp, " measure energy: bad molid", NULL);
02775 return TCL_ERROR;
02776 }
02777 } else if (!strupncmp(argvcur, "k", CMDLEN)) {
02778 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params) != TCL_OK) {
02779 Tcl_AppendResult(interp, " measure energy: bad force constant value", NULL);
02780 return TCL_ERROR;
02781 }
02782 } else if (!strupncmp(argvcur, "x0", CMDLEN)) {
02783 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+1) != TCL_OK) {
02784 Tcl_AppendResult(interp, " measure energy: bad equilibrium value", NULL);
02785 return TCL_ERROR;
02786 }
02787 } else if (!strupncmp(argvcur, "kub", CMDLEN)) {
02788 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+2) != TCL_OK) {
02789 Tcl_AppendResult(interp, " measure energy: bad Urey-Bradley force constant value", NULL);
02790 return TCL_ERROR;
02791 }
02792 } else if (!strupncmp(argvcur, "s0", CMDLEN)) {
02793 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+3) != TCL_OK) {
02794 Tcl_AppendResult(interp, " measure energy: bad Urey-Bradley equilibrium distance", NULL);
02795 return TCL_ERROR;
02796 }
02797 } else if (!strupncmp(argvcur, "n", CMDLEN)) {
02798 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+1) != TCL_OK) {
02799 Tcl_AppendResult(interp, " measure energy: bad dihedral periodicity", NULL);
02800 return TCL_ERROR;
02801 }
02802 } else if (!strupncmp(argvcur, "delta", CMDLEN)) {
02803 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+2) != TCL_OK) {
02804 Tcl_AppendResult(interp, " measure energy: bad dihedral phase shift", NULL);
02805 return TCL_ERROR;
02806 }
02807 } else if (!strupncmp(argvcur, "eps1", CMDLEN)) {
02808 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params) != TCL_OK) {
02809 Tcl_AppendResult(interp, " measure energy: bad vdw well depth", NULL);
02810 return TCL_ERROR;
02811 }
02812 } else if (!strupncmp(argvcur, "rmin1", CMDLEN)) {
02813 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+1) != TCL_OK) {
02814 Tcl_AppendResult(interp, " measure energy: bad vdw equilibrium distance", NULL);
02815 return TCL_ERROR;
02816 }
02817 } else if (!strupncmp(argvcur, "eps2", CMDLEN)) {
02818 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+2) != TCL_OK) {
02819 Tcl_AppendResult(interp, " measure energy: bad vdw well depth", NULL);
02820 return TCL_ERROR;
02821 }
02822 } else if (!strupncmp(argvcur, "rmin2", CMDLEN)) {
02823 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+3) != TCL_OK) {
02824 Tcl_AppendResult(interp, " measure energy: bad vdw equilibrium distance", NULL);
02825 return TCL_ERROR;
02826 }
02827 } else if (!strupncmp(argvcur, "q1", CMDLEN)) {
02828 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params) != TCL_OK) {
02829 Tcl_AppendResult(interp, " measure energy: bad charge value", NULL);
02830 return TCL_ERROR;
02831 }
02832 params[2]=1.0;
02833 } else if (!strupncmp(argvcur, "q2", CMDLEN)) {
02834 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+1) != TCL_OK) {
02835 Tcl_AppendResult(interp, " measure energy: bad charge value", NULL);
02836 return TCL_ERROR;
02837 }
02838 params[3]=1.0;
02839 } else if (!strupncmp(argvcur, "cutoff", CMDLEN)) {
02840 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+4) != TCL_OK) {
02841 Tcl_AppendResult(interp, " measure energy: bad electrostatic cutoff value", NULL);
02842 return TCL_ERROR;
02843 }
02844 } else if (!strupncmp(argvcur, "switchdist", CMDLEN)) {
02845 if (Tcl_GetDoubleFromObj(interp, objv[i+1], params+5) != TCL_OK) {
02846 Tcl_AppendResult(interp, " measure energy: bad switching distance value", NULL);
02847 return TCL_ERROR;
02848 }
02849 } else if (!strupncmp(argvcur, "first", CMDLEN)) {
02850 if (Tcl_GetIntFromObj(interp, objv[i+1], &first) != TCL_OK) {
02851 Tcl_AppendResult(interp, " measure energy: bad first frame value", NULL);
02852 return TCL_ERROR;
02853 }
02854 } else if (!strupncmp(argvcur, "last", CMDLEN)) {
02855 if (Tcl_GetIntFromObj(interp, objv[i+1], &last) != TCL_OK) {
02856 Tcl_AppendResult(interp, " measure energy: bad last frame value", NULL);
02857 return TCL_ERROR;
02858 }
02859 } else if (!strupncmp(argvcur, "frame", CMDLEN)) {
02860 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "all")) {
02861 allframes = true;
02862 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &frame) != TCL_OK) {
02863 Tcl_AppendResult(interp, " measure energy: bad frame value", NULL);
02864 return TCL_ERROR;
02865 }
02866 } else {
02867 Tcl_AppendResult(interp, " measure energy: invalid syntax, no such keyword: ", argvcur, NULL);
02868 return TCL_ERROR;
02869 }
02870 }
02871 }
02872 
02873 if ((allframes || frame>=0) && (first>=0 || last>=0)) {
02874 Tcl_AppendResult(interp, "measure energy: Ambiguous syntax: You cannot specify a frame AND a frame range (using first or last).", NULL);
02875 Tcl_AppendResult(interp, "\nUsage:\nmeasure bond <molid1>/<atomid1> <molid2>/<atomid2> [?frame <frame|all>? | ?first <first>? ?last <last>?]", NULL);
02876 return TCL_ERROR; 
02877 }
02878 
02879 if (allframes) first=0;
02880 
02881 // If no molecule was specified use top as default
02882 if (defmolid<0) defmolid = app->molecule_top();
02883 
02884 // Assign atom IDs and molecule IDs
02885 int i,numelem;
02886 Tcl_Obj **atmmol;
02887 for (i=0; i<numatms; i++) {
02888 if (Tcl_ListObjGetElements(interp, data[i], &numelem, &atmmol) != TCL_OK) {
02889 return TCL_ERROR;
02890 }
02891 
02892 if (Tcl_GetIntFromObj(interp, atmmol[0], atmid+i) != TCL_OK) {
02893 Tcl_AppendResult(interp, " measure energy: bad atom index", NULL);
02894 return TCL_ERROR;
02895 }
02896 
02897 if (numelem==2) {
02898 if (Tcl_GetIntFromObj(interp, atmmol[1], molid+i) != TCL_OK) {
02899 Tcl_AppendResult(interp, " measure energy: bad molid", NULL);
02900 return TCL_ERROR;
02901 }
02902 } else molid[i] = defmolid;
02903 }
02904 
02905 
02906 // Compute the value
02907 ResizeArray<float> gValues(1024);
02908 int ret_val;
02909 ret_val = measure_energy(app->moleculeList, molid, atmid, reqatoms, &gValues, frame, first, last, defmolid, params, geomtype);
02910 if (ret_val<0) {
02911 printf("ERROR\n %s\n", measure_error(ret_val));
02912 Tcl_AppendResult(interp, measure_error(ret_val), NULL);
02913 return TCL_ERROR;
02914 }
02915 
02916 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
02917 int numvalues = gValues.num();
02918 for (int count = 0; count < numvalues; count++) {
02919 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(gValues[count]));
02920 }
02921 Tcl_SetObjResult(interp, tcl_result);
02922 
02923 return TCL_OK;
02924 }
02925 
02926 
02927 //
02928 // Function: vmd_measure_surface <selection> <gridsize> <radius> <depth>
02929 //
02930 static int vmd_measure_surface(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02931 if (argc!=5) {
02932 Tcl_WrongNumArgs(interp, 2, objv-1,
02933 (char *) "<sel> <gridsize> <radius> <depth>");
02934 return TCL_ERROR;
02935 }
02936 
02937 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
02938 if (!sel) {
02939 Tcl_AppendResult(interp, "measure surface: no atom selection", NULL);
02940 return TCL_ERROR;
02941 }
02942 
02943 const float *framepos = sel->coordinates(app->moleculeList);
02944 if (!framepos) return TCL_ERROR;
02945 
02946 double gridsz;
02947 if (Tcl_GetDoubleFromObj(interp, objv[2], &gridsz) != TCL_OK) 
02948 return TCL_ERROR;
02949 
02950 double radius;
02951 if (Tcl_GetDoubleFromObj(interp, objv[3], &radius) != TCL_OK) 
02952 return TCL_ERROR;
02953 
02954 double depth;
02955 if (Tcl_GetDoubleFromObj(interp, objv[4], &depth) != TCL_OK) 
02956 return TCL_ERROR;
02957 
02958 int *surface;
02959 int n_surf;
02960 
02961 int ret_val = measure_surface(sel, app->moleculeList, framepos,
02962 gridsz, radius, depth, &surface, &n_surf);
02963 if (ret_val < 0) {
02964 Tcl_AppendResult(interp, "measure surface: ", measure_error(ret_val));
02965 return TCL_ERROR;
02966 }
02967 
02968 Tcl_Obj *surf_list = Tcl_NewListObj(0, NULL);
02969 
02970 int i;
02971 for(i=0; i < n_surf; i++) {
02972 Tcl_ListObjAppendElement(interp, surf_list, Tcl_NewIntObj(surface[i]));
02973 }
02974 Tcl_SetObjResult(interp, surf_list);
02975 delete [] surface;
02976 
02977 return TCL_OK;
02978 }
02979 
02980 
02981 // Function: vmd_measure_pbc2onc <center> ?molid <default>? ?frame <frame|last>?
02982 // Returns: the transformation matrix to wrap a nonorthogonal pbc unicell
02983 // into an orthonormal cell.
02984 static int vmd_measure_pbc2onc_transform(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
02985 if (argc < 2) {
02986 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<center> [molid <default>] [frame <frame|last>]");
02987 return TCL_ERROR;
02988 }
02989 
02990 // Read the center
02991 int ndim;
02992 Tcl_Obj **centerObj;
02993 if (Tcl_ListObjGetElements(interp, objv[1], &ndim, &centerObj) != TCL_OK) {
02994 Tcl_AppendResult(interp, " measure pbc2onc: bad syntax", NULL);
02995 Tcl_AppendResult(interp, " Usage: measure pbc2onc <center> [molid <default>] [frame <frame|last>]", NULL);
02996 return TCL_ERROR;
02997 }
02998 
02999 if (ndim!=3) {
03000 Tcl_AppendResult(interp, " measure pbc2onc: need three numbers for a vector", NULL);
03001 return TCL_ERROR;
03002 }
03003 
03004 int i;
03005 double tmp;
03006 float center[3];
03007 for (i=0; i<3; i++) {
03008 if (Tcl_GetDoubleFromObj(interp, centerObj[i], &tmp) != TCL_OK) {
03009 Tcl_AppendResult(interp, " measure pbc2onc: non-numeric in center", NULL);
03010 return TCL_ERROR;
03011 }
03012 center[i] = (float)tmp;
03013 }
03014 
03015 int molid = app->molecule_top();
03016 int frame = -2;
03017 if (argc>3) {
03018 for (i=2; i<argc; i+=2) {
03019 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
03020 if (!strupncmp(argvcur, "molid", CMDLEN)) {
03021 // XXX lame test formulation should be corrected
03022 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "top")) {
03023 // top is already default
03024 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &molid) != TCL_OK) {
03025 Tcl_AppendResult(interp, " measure pbc2onc: bad molid", NULL);
03026 return TCL_ERROR;
03027 }
03028 } else if (!strupncmp(argvcur, "frame", CMDLEN)) {
03029 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "last")) {
03030 frame=-1;
03031 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &frame) != TCL_OK) {
03032 Tcl_AppendResult(interp, " measure pbc2onc: bad frame value", NULL);
03033 return TCL_ERROR;
03034 }
03035 } else {
03036 Tcl_AppendResult(interp, " measure pbc2onc: invalid syntax, no such keyword: ", argvcur, NULL);
03037 return TCL_ERROR;
03038 }
03039 }
03040 }
03041 
03042 
03043 Matrix4 transform;
03044 int ret_val = measure_pbc2onc(app->moleculeList, molid, frame, center, transform);
03045 if (ret_val < 0) {
03046 Tcl_AppendResult(interp, "measure pbc2onc: ", measure_error(ret_val), NULL);
03047 return TCL_ERROR;
03048 }
03049 
03050 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
03051 for (i=0; i<4; i++) {
03052 Tcl_Obj *rowListObj = Tcl_NewListObj(0, NULL);
03053 Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewDoubleObj(transform.mat[i+0]));
03054 Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewDoubleObj(transform.mat[i+4]));
03055 Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewDoubleObj(transform.mat[i+8]));
03056 Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewDoubleObj(transform.mat[i+12]));
03057 Tcl_ListObjAppendElement(interp, tcl_result, rowListObj);
03058 }
03059 Tcl_SetObjResult(interp, tcl_result);
03060 
03061 return TCL_OK;
03062 }
03063 
03064 // Function: vmd_measure_pbc_neighbors <center> <cutoff>
03065 // Returns: All image atoms that are within cutoff Angstrom of the pbc unitcell.
03066 // Two lists are returned. The first list holds the atom coordinates while
03067 // the second one is an indexlist mapping the image atoms to the atoms in
03068 // the unitcell.
03069 
03070 // Input: Since the pbc cell <center> is not stored in DCDs and cannot be set in
03071 // VMD it must be provided by the user as the first argument.
03072 // The second argument <cutoff> is the maximum distance (in Angstrom) from
03073 // the PBC unit cell for atoms to be considered.
03074 
03075 // Options: ?sel <selection>? :
03076 // If an atomselection is provided after the keyword 'sel' then only those
03077 // image atoms are returned that are within cutoff of the selected atoms
03078 // of the main cell. In case cutoff is a vector the largest value will be
03079 // used.
03080 
03081 // ?align <matrix>? :
03082 // In case the molecule was aligned you can supply the alignment matrix
03083 // which is then used to correct for the rotation and shift of the pbc cell.
03084 
03085 // ?boundingbox PBC | {<mincoord> <maxcoord>}?
03086 // With this option the atoms are wrapped into a rectangular bounding box.
03087 // If you provide "PBC" as an argument then the bounding box encloses the
03088 // PBC box but then the cutoff is added to the bounding box. Negative 
03089 // values for the cutoff dimensions are allowed and lead to a smaller box.
03090 // Instead you can also provide a custom bounding box in form of the 
03091 // minmax coordinates (list containing two coordinate vectors such as 
03092 // returned by the measure minmax command). Here, again, the cutoff is
03093 // added to the bounding box.
03094 static int vmd_measure_pbc_neighbors(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
03095 if (argc < 3) {
03096 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<center> <cutoff> ?sel <selection>? ?align <matrix>? ?molid <default>? ?frame <frame|last>? ?boundingbox PBC|{<mincoord> <maxcoord>}?");
03097 return TCL_ERROR;
03098 }
03099 
03100 // Read the center
03101 int ndim;
03102 Tcl_Obj **centerObj;
03103 if (Tcl_ListObjGetElements(interp, objv[1], &ndim, &centerObj) != TCL_OK) {
03104 Tcl_AppendResult(interp, " measure pbcneighbors: bad syntax", NULL);
03105 Tcl_AppendResult(interp, " Usage: measure pbcneighbors <center> <cutoff> [<sel <sel>] [align <matrix>] [molid <default>] [frame <frame|last>] [boundingbox <PBC|{<mincoord> <maxcoord>}>]", NULL);
03106 return TCL_ERROR;
03107 }
03108 
03109 if (ndim!=3) {
03110 Tcl_AppendResult(interp, " measure pbcneighbors: need three numbers for a vector", NULL);
03111 return TCL_ERROR;
03112 }
03113 
03114 int i;
03115 double tmp;
03116 float center[3];
03117 for (i=0; i<3; i++) {
03118 if (Tcl_GetDoubleFromObj(interp, centerObj[i], &tmp) != TCL_OK) {
03119 Tcl_AppendResult(interp, " measure pbcneighbors: non-numeric in center", NULL);
03120 return TCL_ERROR;
03121 }
03122 center[i] = float(tmp);
03123 }
03124 
03125 // Read the cutoff
03126 Tcl_Obj **cutoffObj;
03127 if (Tcl_ListObjGetElements(interp, objv[2], &ndim, &cutoffObj) != TCL_OK) {
03128 Tcl_AppendResult(interp, " measure pbcneighbors: bad syntax", NULL);
03129 Tcl_AppendResult(interp, " Usage: measure pbcneighbors <center> <cutoff> [<sel <sel>] [align <matrix>] [molid <default>] [frame <frame|last>] [boundingbox <PBC|{<mincoord> <maxcoord>}>]", NULL);
03130 return TCL_ERROR;
03131 }
03132 
03133 if (ndim!=3 && ndim!=1) {
03134 Tcl_AppendResult(interp, " measure pbcneighbors: need either one or three numbers for cutoff", NULL);
03135 return TCL_ERROR;
03136 }
03137 
03138 float cutoff[3];
03139 for (i=0; i<ndim; i++) {
03140 if (Tcl_GetDoubleFromObj(interp, cutoffObj[i], &tmp) != TCL_OK) {
03141 Tcl_AppendResult(interp, " measure pbcneighbors: non-numeric in cutoff", NULL);
03142 return TCL_ERROR;
03143 }
03144 cutoff[i] = float(tmp);
03145 }
03146 
03147 if (ndim==1) { cutoff[2] = cutoff[1] = cutoff[0]; }
03148 
03149 bool molidprovided=0;
03150 float *boxminmax=NULL;
03151 int molid = app->molecule_top();
03152 int frame = -2;
03153 AtomSel *sel = NULL;
03154 Matrix4 alignment;
03155 if (argc>4) {
03156 for (i=3; i<argc; i+=2) {
03157 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
03158 if (!strupncmp(argvcur, "sel", CMDLEN)) {
03159 // Read the atom selection
03160 sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[i+1],NULL));
03161 if (!sel) {
03162 Tcl_AppendResult(interp, "measure pbcneighbors: invalid atom selection", NULL);
03163 return TCL_ERROR;
03164 }
03165 if (!app->molecule_valid_id(sel->molid())) {
03166 Tcl_AppendResult(interp, "measure pbcneighbors: ",
03167 measure_error(MEASURE_ERR_NOMOLECULE), NULL);
03168 return TCL_ERROR;
03169 }
03170 if (!sel->selected) {
03171 Tcl_AppendResult(interp, "measure pbcneighbors: selection contains no atoms.", NULL);
03172 return TCL_ERROR;
03173 }
03174 } else if (!strupncmp(argvcur, "molid", CMDLEN)) {
03175 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "top")) {
03176 // top is already default
03177 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &molid) != TCL_OK) {
03178 Tcl_AppendResult(interp, " measure pbcneighbors: bad molid", NULL);
03179 return TCL_ERROR;
03180 }
03181 molidprovided = 1;
03182 } else if (!strupncmp(argvcur, "frame", CMDLEN)) {
03183 if (!strupcmp(Tcl_GetStringFromObj(objv[i+1],NULL), "last")) {
03184 frame=-1;
03185 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &frame) != TCL_OK) {
03186 Tcl_AppendResult(interp, " measure pbcneighbors: bad frame value", NULL);
03187 return TCL_ERROR;
03188 }
03189 } else if (!strupncmp(argvcur, "align", CMDLEN)) {
03190 // Get the alignment matrix (as returned by 'measure fit')
03191 if (tcl_get_matrix("measure pbcneighbors: ", interp, objv[i+1], alignment.mat) != TCL_OK) {
03192 return TCL_ERROR;
03193 }
03194 } else if (!strupncmp(argvcur, "boundingbox", CMDLEN)) {
03195 // Determine if the atoms shall be wrapped to the smallest rectangular
03196 // bounding box that still encloses the unitcell plus the given cutoff.
03197 char *argv2 = Tcl_GetStringFromObj(objv[i+1],NULL);
03198 if (!strupncmp(argv2, "on", CMDLEN) || !strupncmp(argv2, "pbc", CMDLEN)) {
03199 boxminmax = new float[6];
03200 compute_pbcminmax(app->moleculeList, molid, frame, center, &alignment,
03201 boxminmax, boxminmax+3);
03202 } else {
03203 // Read the bounding box
03204 int j, k, ncoor;
03205 Tcl_Obj **boxListObj;
03206 if (Tcl_ListObjGetElements(interp, objv[i+1], &ncoor, &boxListObj) != TCL_OK) {
03207 Tcl_AppendResult(interp, " measure pbcneighbors: invalid bounding box parameter", NULL);
03208 return TCL_ERROR;
03209 }
03210 if (ncoor!=2) {
03211 Tcl_AppendResult(interp, " measure pbcneighbors: need 2 points for bounding box", NULL);
03212 return TCL_ERROR;
03213 }
03214 int ndim = 0;
03215 double tmp;
03216 Tcl_Obj **boxObj;
03217 boxminmax = new float[6];
03218 for (j=0; j<2; j++) {
03219 if (Tcl_ListObjGetElements(interp, boxListObj[j], &ndim, &boxObj) != TCL_OK) {
03220 Tcl_AppendResult(interp, " measure pbcneighbors: bad syntax in boundingbox", NULL);
03221 return TCL_ERROR;
03222 }
03223 
03224 for (k=0; k<3; k++) {
03225 if (Tcl_GetDoubleFromObj(interp, boxObj[k], &tmp) != TCL_OK) {
03226 Tcl_AppendResult(interp, " measure pbcneighbors: non-numeric in boundingbox", NULL);
03227 return TCL_ERROR;
03228 }
03229 boxminmax[3L*j+k] = (float)tmp;
03230 }
03231 }
03232 }
03233 
03234 } else {
03235 Tcl_AppendResult(interp, " measure pbcneighbors: invalid syntax, no such keyword: ", argvcur, NULL);
03236 return TCL_ERROR;
03237 }
03238 }
03239 }
03240 
03241 // If no molid was provided explicitely but a selection was given
03242 // then use that molid.
03243 if (sel && !molidprovided) {
03244 molid = sel->molid();
03245 }
03246 
03247 ResizeArray<float> extcoord_array;
03248 ResizeArray<int> indexmap_array;
03249 
03250 int ret_val = measure_pbc_neighbors(app->moleculeList, sel, molid, frame, &alignment, center, cutoff, boxminmax, &extcoord_array, &indexmap_array);
03251 delete [] boxminmax;
03252 
03253 if (ret_val < 0) {
03254 Tcl_AppendResult(interp, "measure pbcneighbors: ", measure_error(ret_val), NULL);
03255 return TCL_ERROR;
03256 }
03257 printf("measure pbcneighbors: %ld neighbor atoms found\n", 
03258 long(indexmap_array.num()));
03259 
03260 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
03261 Tcl_Obj *coorListObj = Tcl_NewListObj(0, NULL);
03262 Tcl_Obj *indexListObj = Tcl_NewListObj(0, NULL);
03263 for (i=0; i<indexmap_array.num(); i++) {
03264 Tcl_Obj *rowListObj = Tcl_NewListObj(0, NULL);
03265 Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewDoubleObj(extcoord_array[3L*i]));
03266 Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewDoubleObj(extcoord_array[3L*i+1]));
03267 Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewDoubleObj(extcoord_array[3L*i+2]));
03268 Tcl_ListObjAppendElement(interp, coorListObj, rowListObj);
03269 Tcl_ListObjAppendElement(interp, indexListObj, Tcl_NewIntObj(indexmap_array[i]));
03270 }
03271 Tcl_ListObjAppendElement(interp, tcl_result, coorListObj);
03272 Tcl_ListObjAppendElement(interp, tcl_result, indexListObj);
03273 Tcl_SetObjResult(interp, tcl_result);
03274 
03275 return TCL_OK;
03276 }
03277 
03278 
03279 
03280 // Function: vmd_measure_inertia <selection> [moments] [eigenvals]
03281 // Returns: The center of mass and the principles axes of inertia for the 
03282 // selected atoms. 
03283 // Options: -moments -- also return the moments of inertia tensor
03284 // -eigenvals -- also return the corresponding eigenvalues
03285 static int vmd_measure_inertia(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
03286 bool moments = FALSE;
03287 bool eigenvals = FALSE;
03288 if (argc < 2 || argc > 4) {
03289 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<selection> [moments] [eigenvals]");
03290 return TCL_ERROR;
03291 }
03292 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
03293 if (!sel) {
03294 Tcl_AppendResult(interp, "measure inertia: no atom selection", NULL);
03295 return TCL_ERROR;
03296 }
03297 
03298 if (!app->molecule_valid_id(sel->molid())) {
03299 Tcl_AppendResult(interp, "measure inertia: ",
03300 measure_error(MEASURE_ERR_NOMOLECULE), NULL);
03301 return TCL_ERROR;
03302 }
03303 
03304 if (argc>2) {
03305 int i;
03306 for (i=2; i<argc; i++) {
03307 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
03308 // Allow syntax with and without leading dash
03309 if (argvcur[0]=='-') argvcur++;
03310 
03311 if (!strupncmp(argvcur, "moments", CMDLEN)) {
03312 moments = TRUE; 
03313 }
03314 else if (!strupncmp(argvcur, "eigenvals", CMDLEN)) {
03315 eigenvals = TRUE; 
03316 }
03317 else {
03318 Tcl_AppendResult(interp, " measure inertia: unrecognized option\n", NULL);
03319 Tcl_AppendResult(interp, " Usage: measure inertia <selection> [moments] [eigenvals]", NULL);
03320 return TCL_ERROR;
03321 }
03322 }
03323 }
03324 
03325 float priaxes[3][3];
03326 float itensor[4][4];
03327 float evalue[3], rcom[3];
03328 int ret_val = measure_inertia(sel, app->moleculeList, NULL, rcom, priaxes, itensor, evalue);
03329 if (ret_val < 0) {
03330 Tcl_AppendResult(interp, "measure inertia: ", measure_error(ret_val), NULL);
03331 return TCL_ERROR;
03332 }
03333 
03334 // return COM and list of 3 principal axes
03335 int i;
03336 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
03337 
03338 Tcl_Obj *rcomObj = Tcl_NewListObj(0, NULL);
03339 Tcl_ListObjAppendElement(interp, rcomObj, Tcl_NewDoubleObj(rcom[0]));
03340 Tcl_ListObjAppendElement(interp, rcomObj, Tcl_NewDoubleObj(rcom[1]));
03341 Tcl_ListObjAppendElement(interp, rcomObj, Tcl_NewDoubleObj(rcom[2]));
03342 Tcl_ListObjAppendElement(interp, tcl_result, rcomObj);
03343 
03344 Tcl_Obj *axesListObj = Tcl_NewListObj(0, NULL);
03345 for (i=0; i<3; i++) {
03346 Tcl_Obj *axesObj = Tcl_NewListObj(0, NULL);
03347 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(priaxes[i][0]));
03348 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(priaxes[i][1]));
03349 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(priaxes[i][2]));
03350 Tcl_ListObjAppendElement(interp, axesListObj, axesObj);
03351 }
03352 Tcl_ListObjAppendElement(interp, tcl_result, axesListObj);
03353 
03354 if (moments) {
03355 Tcl_Obj *momListObj = Tcl_NewListObj(0, NULL);
03356 for (i=0; i<3; i++) {
03357 Tcl_Obj *momObj = Tcl_NewListObj(0, NULL);
03358 Tcl_ListObjAppendElement(interp, momObj, Tcl_NewDoubleObj(itensor[i][0]));
03359 Tcl_ListObjAppendElement(interp, momObj, Tcl_NewDoubleObj(itensor[i][1]));
03360 Tcl_ListObjAppendElement(interp, momObj, Tcl_NewDoubleObj(itensor[i][2]));
03361 Tcl_ListObjAppendElement(interp, momListObj, momObj);
03362 }
03363 Tcl_ListObjAppendElement(interp, tcl_result, momListObj);
03364 }
03365 
03366 if (eigenvals) {
03367 Tcl_Obj *eigvListObj = Tcl_NewListObj(0, NULL);
03368 Tcl_ListObjAppendElement(interp, eigvListObj, Tcl_NewDoubleObj(evalue[0]));
03369 Tcl_ListObjAppendElement(interp, eigvListObj, Tcl_NewDoubleObj(evalue[1]));
03370 Tcl_ListObjAppendElement(interp, eigvListObj, Tcl_NewDoubleObj(evalue[2]));
03371 Tcl_ListObjAppendElement(interp, tcl_result, eigvListObj);
03372 }
03373 Tcl_SetObjResult(interp, tcl_result);
03374 
03375 return TCL_OK;
03376 }
03377 
03378 
03379 // measure symmetry <sel> [plane|I|C<n>|S<n> [<vector>]] [-tol <value>]
03380 // [-nobonds] [-verbose <level>]
03381 //
03382 // This function evaluates the molecular symmetry of an atom selection.
03383 // The underlying algorithm finds the symmetry elements such as 
03384 // inversion center, mirror planes, rotary axes and rotary reflections.
03385 // Based on the found elements it guesses the underlying point group.
03386 // The guess is fairly robust and can handle molecules whose
03387 // coordinatesthat deviate to a certain extent from the ideal symmetry.
03388 // The closest match with the highest symmetry will be returned.
03389 //
03390 // Options:
03391 // --------
03392 // -tol <value>
03393 // Allows one to control tolerance of the algorithm when
03394 // considering wether something is symmetric or not.
03395 // A smaller value signifies a lower tolerance, the default
03396 // is 0.1.
03397 // -nobonds 
03398 // If this flag is set then the bond order and orientation
03399 // are not considered when comparing structures.
03400 // -verbose <level>
03401 // Controls the amount of console output.
03402 // A level of 0 means no output, 1 gives some statistics at the
03403 // end of the search (default). Level 2 gives additional info
03404 // about each stage, level 3 more output for each iteration
03405 // and 3, 4 yield yet more additional info.
03406 // -I
03407 // Instead of guessing the symmetry pointgroup of the selection
03408 // determine if the selection's center off mass represents an
03409 // inversion center. The returned value is a score between 0
03410 // and 1 where 1 denotes a perfect match.
03411 // -plane <vector>
03412 // Instead of guessing the symmetry pointgroup of the selection
03413 // determine if the plane with the defined by its normal
03414 // <vector> is a mirror plane of the selection. The
03415 // returned value is a score between 0 and 1 where 1 denotes
03416 // a perfect match.
03417 // -Cn | -Sn <vector>
03418 // Instead of guessing the symmetry pointgroup of the selection
03419 // determine if the rotation or rotary reflection axis Cn/Sn
03420 // with order n defined by <vector> exists for the
03421 // selection. E.g., if you want to query wether the Y-axis
03422 // has a C3 rotational symmetry you specify "C3 {0 1 0}".
03423 // The returned value is a score between 0 and 1 where 1
03424 // denotes a perfect match.
03425 //
03426 // Result:
03427 // -------
03428 // The return value is a TCL list of pairs consisting of a label
03429 // string and a value or list. For each label the data following
03430 // it are described below:
03431 //
03432 // * [pointgroup] The guessed pointgroup:
03433 // * [order] Point group order (the n from above)
03434 // * [elements] Summary of found elements.
03435 // For instance {(i) (C3) 3*(C2) (S6) 3*(sigma)} for D3d.
03436 // * [missing] Elements missing with respect to ideal number of
03437 // elements (same format as above). If this is not an empty
03438 // list then something has gone awfully wrong with the symmetry
03439 // finding algorithm.
03440 // * [additional] Additional elements that would not be expected
03441 // for this point group (same format as above). If this is not
03442 // an empty list then something has gone awfully wrong with the
03443 // symmetry finding algorithm.
03444 // * [com] Center of mass of the selection based on the idealized
03445 // coordinates.
03446 // * [inertia] List of the three axes of inertia, the eigenvalues
03447 // of the moments of inertia tensor and a list of three 0/1 flags 
03448 // specifying for each axis wether it is unique or not.
03449 // * [inversion] Flag 0/1 signifying if there is an inversion center.
03450 // * [axes] Normalized vectors defining rotary axes
03451 // * [rotreflect] Normalized vectors defining rotary reflections
03452 // * [planes] Normalized vectors defining mirror planes.
03453 // * [ideal] Idealized symmetric coordinates for all atoms of
03454 // the selection. The coordinates are listed in the order of 
03455 // increasing atom indices (same order asa returned by the
03456 // atomselect command ``get {x y z}''). Thus you can use the list
03457 // to set the atoms of your selection to the ideal coordinates
03458 // * [unique] Index list defining a set of atoms with unique
03459 // coordinates
03460 // * [orient] Matrix that aligns molecule with GAMESS standard
03461 // orientation
03462 //
03463 // If a certain item is not present (e.g. no planes or no axes)
03464 // then the corresponding value is an empty list.
03465 // The pair format allows to use the result as a TCL array for
03466 // convenient access of the different return items.
03467 
03468 static int vmd_measure_symmetry(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
03469 if (argc<2) {
03470 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<sel> [element [<vector>]] [-tol <value>] [-nobonds] [-verbose <level>]");
03471 return TCL_ERROR;
03472 }
03473 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
03474 if (!sel) {
03475 Tcl_AppendResult(interp, "measure symmetry: no atom selection", NULL);
03476 return TCL_ERROR;
03477 }
03478 if (!app->molecule_valid_id(sel->molid())) {
03479 Tcl_AppendResult(interp, "measure symmetry: ",
03480 measure_error(MEASURE_ERR_NOMOLECULE), NULL);
03481 return TCL_ERROR;
03482 }
03483 
03484 double sigma = 0.1;
03485 float axis[3];
03486 int checkelement = 0;
03487 int checkbonds = 1;
03488 int order = 0;
03489 int verbose = 1;
03490 int impose = 0;
03491 int imposeinvers = 0;
03492 int numimposeplan = 0;
03493 int numimposeaxes = 0;
03494 int numimposerref = 0;
03495 float *imposeplan = NULL;
03496 float *imposeaxes = NULL;
03497 float *imposerref = NULL;
03498 int *imposeaxesord = NULL;
03499 int *imposerreford = NULL;
03500 AtomSel *idealsel = NULL;
03501 
03502 if (argc>2) {
03503 int bailout = 0;
03504 int i;
03505 for (i=2; (i<argc && !bailout); i++) {
03506 char *argvcur = Tcl_GetStringFromObj(objv[i],NULL);
03507 // Allow syntax with and without leading dash
03508 if (argvcur[0]=='-') argvcur++;
03509 
03510 if (!strupncmp(argvcur, "tol", CMDLEN)) {
03511 if (Tcl_GetDoubleFromObj(interp, objv[i+1], &sigma) != TCL_OK) {
03512 Tcl_AppendResult(interp, " measure symmetry: bad tolerance value", NULL);
03513 bailout = 1; continue;
03514 }
03515 i++;
03516 }
03517 else if (!strupncmp(argvcur, "verbose", CMDLEN)) {
03518 if (Tcl_GetIntFromObj(interp, objv[i+1], &verbose) != TCL_OK) {
03519 Tcl_AppendResult(interp, " measure symmetry: bad verbosity level value", NULL);
03520 bailout = 1; continue;
03521 }
03522 i++;
03523 }
03524 else if (!strupncmp(argvcur, "nobonds", CMDLEN)) {
03525 checkbonds = 0;
03526 }
03527 else if (!strupcmp(argvcur, "I")) {
03528 checkelement = 1;
03529 }
03530 else if (!strupncmp(argvcur, "plane", CMDLEN)) {
03531 if (tcl_get_vector(Tcl_GetStringFromObj(objv[i+1],NULL), axis, interp)!=TCL_OK) {
03532 bailout = 1; continue;
03533 }
03534 checkelement = 2;
03535 i++;
03536 }
03537 else if (!strupncmp(argvcur, "C", 1) || !strupncmp(argvcur, "S", 1)) {
03538 char *begptr = argvcur+1;
03539 char *endptr;
03540 order = strtol(begptr, &endptr, 10);
03541 if (endptr==begptr || *endptr!='0円') {
03542 Tcl_AppendResult(interp, " measure symmetry: bad symmetry element format (must be I, C*, S*, plane, where * is the order). ", NULL);
03543 bailout = 1; continue;
03544 }
03545 
03546 if (tcl_get_vector(Tcl_GetStringFromObj(objv[i+1],NULL), axis, interp)!=TCL_OK) {
03547 bailout = 1; continue;
03548 }
03549 
03550 if (!strupncmp(argvcur, "C", 1)) checkelement = 3;
03551 if (!strupncmp(argvcur, "S", 1)) checkelement = 4;
03552 i++;
03553 }
03554 else if (!strupncmp(argvcur, "idealsel", CMDLEN)) {
03555 idealsel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[i+1],NULL));
03556 if (!sel) {
03557 Tcl_AppendResult(interp, "measure symmetry: no atom selection for idealized coordinates", NULL);
03558 bailout = 1; continue;
03559 }
03560 if (idealsel->molid()!=sel->molid()) {
03561 Tcl_AppendResult(interp, "measure symmetry: selection and idealsel must be from the same molecule", NULL);
03562 bailout = 1; continue;
03563 }
03564 if (idealsel->selected<sel->selected) {
03565 Tcl_AppendResult(interp, "measure symmetry: selection must be a subset of idealsel", NULL);
03566 bailout = 1; continue;
03567 }
03568 // make sure sel is subset of idealsel
03569 int j;
03570 for (j=0; j<sel->num_atoms; j++) {
03571 if (sel->on[j] && !idealsel->on[j]) {
03572 Tcl_AppendResult(interp, "measure symmetry: selection must be a subset of idealsel", NULL);
03573 bailout = 1; continue;
03574 }
03575 }
03576 
03577 i++;
03578 }
03579 else if (!strupncmp(argvcur, "imposeinversion", CMDLEN)) {
03580 imposeinvers = 1;
03581 impose = 1;
03582 }
03583 else if (!strupncmp(argvcur, "imposeplanes", CMDLEN)) {
03584 // Format: list of normal vectors {{x y z} {x y z} ...}
03585 int nelem;
03586 Tcl_Obj **elemListObj;
03587 if (i+1>=argc ||
03588 Tcl_ListObjGetElements(interp, objv[i+1], &nelem, &elemListObj) != TCL_OK) {
03589 Tcl_AppendResult(interp, " measure symmetry: bad syntax for imposeplanes option", NULL);
03590 bailout = 1; continue;
03591 }
03592 float *elem = new float[3L*nelem];
03593 int k;
03594 for (k=0; k<nelem; k++) {
03595 int nobj;
03596 Tcl_Obj **vecObj;
03597 if (Tcl_ListObjGetElements(interp, elemListObj[k], &nobj, &vecObj) != TCL_OK) {
03598 delete [] elem;
03599 Tcl_AppendResult(interp, " measure symmetry: bad syntax for imposeplanes option", NULL);
03600 bailout = 1; continue;
03601 }
03602 if (nobj!=3) {
03603 delete [] elem;
03604 Tcl_AppendResult(interp, " measure symmetry imposeplanes: vector must have 3 elements", NULL);
03605 bailout = 1; continue;
03606 }
03607 
03608 int m;
03609 for (m=0; m<3; m++) {
03610 double d;
03611 if (Tcl_GetDoubleFromObj(interp, vecObj[m], &d) != TCL_OK) {
03612 delete [] elem;
03613 bailout = 1; continue;
03614 }
03615 elem[3L*k+m] = (float)d;
03616 }
03617 }
03618 if (imposeplan) delete [] imposeplan;
03619 imposeplan = elem;
03620 numimposeplan = nelem;
03621 impose = 1;
03622 i++;
03623 }
03624 else if (!strupncmp(argvcur, "imposeaxes", CMDLEN) ||
03625 !strupncmp(argvcur, "imposerotref", CMDLEN)) {
03626 // Format: list of axes and orders {{x y z} 3 {x y z} 2 ...}
03627 int nelem;
03628 Tcl_Obj **elemListObj;
03629 if (i+1>=argc ||
03630 Tcl_ListObjGetElements(interp, objv[i+1], &nelem, &elemListObj) != TCL_OK ||
03631 nelem%2) {
03632 Tcl_AppendResult(interp, " measure symmetry: bad syntax for imposeaxes/imposerotref option", NULL);
03633 bailout = 1; continue;
03634 }
03635 nelem /= 2;
03636 
03637 if (nelem<=0) {
03638 i++;
03639 continue;
03640 }
03641 float *elem = new float[3L*nelem];
03642 int *axorder = new int[nelem];
03643 int k;
03644 for (k=0; k<nelem; k++) {
03645 int nobj;
03646 Tcl_Obj **vecObj;
03647 if (Tcl_ListObjGetElements(interp, elemListObj[2L*k], &nobj, &vecObj) != TCL_OK) {
03648 delete [] elem;
03649 delete [] axorder;
03650 Tcl_AppendResult(interp, " measure symmetry impose: bad syntax for axis vector", NULL);
03651 bailout = 1; continue;
03652 }
03653 if (Tcl_GetIntFromObj(interp, elemListObj[2L*k+1], &axorder[k]) != TCL_OK) {
03654 delete [] elem;
03655 delete [] axorder;
03656 bailout = 1; continue;
03657 }
03658 if (nobj!=3) {
03659 delete [] elem;
03660 delete [] axorder;
03661 Tcl_AppendResult(interp, " measure symmetry impose: axis vector must have 3 elements", NULL);
03662 bailout = 1; continue;
03663 }
03664 
03665 int m;
03666 for (m=0; m<3; m++) {
03667 double d;
03668 if (Tcl_GetDoubleFromObj(interp, vecObj[m], &d) != TCL_OK) {
03669 delete [] elem;
03670 delete [] axorder;
03671 bailout = 1; continue;
03672 }
03673 elem[3L*k+m] = (float)d;
03674 }
03675 }
03676 if (!strupncmp(argvcur, "imposeaxes", CMDLEN)) {
03677 if (imposeaxes) delete [] imposeaxes;
03678 if (imposeaxesord) delete [] imposeaxesord;
03679 imposeaxes = elem;
03680 imposeaxesord = axorder;
03681 numimposeaxes = nelem;
03682 } else if (!strupncmp(argvcur, "imposerotref", CMDLEN)) {
03683 if (imposerref) delete [] imposerref;
03684 if (imposerreford) delete [] imposerreford;
03685 imposerref = elem;
03686 imposerreford = axorder;
03687 numimposerref = nelem;
03688 }
03689 
03690 impose = 1;
03691 i++;
03692 }
03693 else {
03694 Tcl_AppendResult(interp, " measure symmetry: unrecognized option ", NULL);
03695 Tcl_AppendResult(interp, argvcur);
03696 Tcl_AppendResult(interp, ".\n Usage: measure symmetry <selection> [element [<vector>]] [-tol <value>] [-nobonds] [-verbose <level>]", NULL);
03697 bailout = 1; continue;
03698 }
03699 }
03700 
03701 if (bailout) {
03702 if (imposeplan) delete [] imposeplan;
03703 if (imposeaxes) delete [] imposeaxes;
03704 if (imposerref) delete [] imposerref;
03705 if (imposeaxesord) delete [] imposeaxesord;
03706 if (imposerreford) delete [] imposerreford;
03707 return TCL_ERROR;
03708 }
03709 }
03710 
03711 // Initialize
03712 Symmetry sym = Symmetry(sel, app->moleculeList, verbose);
03713 
03714 // Set tolerance for atomic overlapping
03715 sym.set_overlaptol(float(sigma));
03716 
03717 // Wether to include bonds into the analysis
03718 sym.set_checkbonds(checkbonds);
03719 
03720 if (checkelement) {
03721 // We are interested only in the score for a specific
03722 // symmetry element.
03723 float overlap = 0.0;
03724 if (checkelement==1) {
03725 overlap = sym.score_inversion();
03726 }
03727 else if (checkelement==2) {
03728 overlap = sym.score_plane(axis);
03729 }
03730 else if (checkelement==3) {
03731 overlap = sym.score_axis(axis, order);
03732 }
03733 else if (checkelement==4) {
03734 overlap = sym.score_rotary_reflection(axis, order);
03735 }
03736 
03737 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
03738 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(overlap));
03739 Tcl_SetObjResult(interp, tcl_result);
03740 return TCL_OK;
03741 }
03742 
03743 if (impose) {
03744 sym.impose(imposeinvers, numimposeplan, imposeplan,
03745 numimposeaxes, imposeaxes, imposeaxesord,
03746 numimposerref, imposerref, imposerreford);
03747 if (imposeplan) delete [] imposeplan;
03748 if (imposeaxes) delete [] imposeaxes;
03749 if (imposerref) delete [] imposerref;
03750 if (imposeaxesord) delete [] imposeaxesord;
03751 if (imposerreford) delete [] imposerreford;
03752 }
03753 
03754 int ret_val = sym.guess(float(sigma));
03755 if (ret_val < 0) {
03756 Tcl_AppendResult(interp, "measure symmetry: ", measure_error(ret_val), NULL);
03757 return TCL_ERROR;
03758 }
03759 
03760 int natoms = sel->selected;
03761 Symmetry *s = &sym;
03762 
03763 if (idealsel) {
03764 Symmetry *isym = new Symmetry(idealsel, app->moleculeList, verbose);
03765 isym->set_overlaptol(float(sigma));
03766 isym->set_checkbonds(checkbonds);
03767 int j;
03768 float *plane = new float[3L*sym.numplanes()];
03769 for (j=0; j<sym.numplanes(); j++) {
03770 vec_copy(&plane[3L*j], sym.plane(j));
03771 }
03772 int *axisorder = new int[sym.numaxes()];
03773 float *axis = new float[3L*sym.numaxes()];
03774 for (j=0; j<sym.numaxes(); j++) {
03775 axisorder[j] = sym.get_axisorder(j);
03776 vec_copy(&axis[3L*j], sym.axis(j));
03777 }
03778 int *rrorder = new int[sym.numrotreflect()];
03779 float *rraxis = new float[3L*sym.numrotreflect()];
03780 for (j=0; j<sym.numrotreflect(); j++) {
03781 rrorder[j] = sym.get_rotreflectorder(j);
03782 vec_copy(&rraxis[3L*j], sym.rotreflect(j));
03783 }
03784 // XXX must check if sel is subset of idealsel
03785 int k=0, m=0;
03786 for (j=0; j<sel->num_atoms; j++) {
03787 if (idealsel->on[j]) {
03788 if (sel->on[j]) {
03789 vec_copy(isym->idealpos(k), sym.idealpos(m));
03790 m++;
03791 }
03792 k++;
03793 }
03794 }
03795 isym->impose(sym.has_inversion(),
03796 sym.numplanes(), plane,
03797 sym.numaxes(), axis, axisorder,
03798 sym.numrotreflect(), rraxis, rrorder);
03799 
03800 ret_val = isym->guess(float(sigma));
03801 if (ret_val < 0) {
03802 Tcl_AppendResult(interp, "measure symmetry: ", measure_error(ret_val), NULL);
03803 return TCL_ERROR;
03804 }
03805 
03806 natoms = idealsel->selected;
03807 s = isym;
03808 }
03809 
03810 int pgorder;
03811 char pointgroup[6];
03812 s->get_pointgroup(pointgroup, &pgorder);
03813 
03814 int i;
03815 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
03816 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("pointgroup", -1));
03817 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj(pointgroup, -1));
03818 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("order", -1));
03819 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewIntObj(pgorder));
03820 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("rmsd", -1));
03821 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(s->get_rmsd()));
03822 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("elements", -1));
03823 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj(s->get_element_summary(), -1));
03824 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("missing", -1));
03825 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj(s->get_missing_elements(), -1));
03826 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("additional", -1));
03827 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj(s->get_additional_elements(), -1));
03828 
03829 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("com", -1));
03830 Tcl_Obj *invObj = Tcl_NewListObj(0, NULL);
03831 float *com = s->center_of_mass();
03832 Tcl_ListObjAppendElement(interp, invObj, Tcl_NewDoubleObj(com[0]));
03833 Tcl_ListObjAppendElement(interp, invObj, Tcl_NewDoubleObj(com[1]));
03834 Tcl_ListObjAppendElement(interp, invObj, Tcl_NewDoubleObj(com[2]));
03835 Tcl_ListObjAppendElement(interp, tcl_result, invObj);
03836 
03837 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("inertia", -1));
03838 Tcl_Obj *inertListListObj = Tcl_NewListObj(0, NULL);
03839 float *inertia = s->get_inertia_axes();
03840 float *eigenval = s->get_inertia_eigenvals();
03841 int *unique = s->get_inertia_unique();
03842 for (i=0; i<3; i++) {
03843 Tcl_Obj *inertObj = Tcl_NewListObj(0, NULL);
03844 Tcl_ListObjAppendElement(interp, inertObj, Tcl_NewDoubleObj(inertia[3L*i]));
03845 Tcl_ListObjAppendElement(interp, inertObj, Tcl_NewDoubleObj(inertia[3L*i+1]));
03846 Tcl_ListObjAppendElement(interp, inertObj, Tcl_NewDoubleObj(inertia[3L*i+2]));
03847 Tcl_Obj *inertListObj = Tcl_NewListObj(0, NULL);
03848 Tcl_ListObjAppendElement(interp, inertListObj, inertObj);
03849 Tcl_ListObjAppendElement(interp, inertListObj, Tcl_NewDoubleObj(eigenval[i]));
03850 Tcl_ListObjAppendElement(interp, inertListObj, Tcl_NewIntObj(unique[i]));
03851 Tcl_ListObjAppendElement(interp, inertListListObj, inertListObj);
03852 }
03853 Tcl_ListObjAppendElement(interp, tcl_result, inertListListObj);
03854 
03855 
03856 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("inversion", -1));
03857 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewIntObj(s->has_inversion()));
03858 
03859 
03860 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("axes", -1));
03861 Tcl_Obj *axesListListObj = Tcl_NewListObj(0, NULL);
03862 for (i=0; i<s->numaxes(); i++) {
03863 Tcl_Obj *axesObj = Tcl_NewListObj(0, NULL);
03864 //printf("Tcl %i: %.2f %.2f %.2f\n", i, s->axis(i)[0], s->axis(i)[1], s->axis(i)[2]);
03865 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(s->axis(i)[0]));
03866 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(s->axis(i)[1]));
03867 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(s->axis(i)[2]));
03868 Tcl_Obj *axesListObj = Tcl_NewListObj(0, NULL);
03869 Tcl_ListObjAppendElement(interp, axesListObj, axesObj);
03870 Tcl_ListObjAppendElement(interp, axesListObj, Tcl_NewIntObj(s->get_axisorder(i)));
03871 int axistype = s->get_axistype(i);
03872 if (axistype & PRINCIPAL_AXIS && 
03873 !(!s->is_spherical_top() && (axistype & PERPENDICULAR_AXIS))) {
03874 Tcl_ListObjAppendElement(interp, axesListObj, Tcl_NewStringObj("principal", -1));
03875 } else if (!s->is_spherical_top() && (axistype & PERPENDICULAR_AXIS)) {
03876 Tcl_ListObjAppendElement(interp, axesListObj, Tcl_NewStringObj("perpendicular", -1));
03877 } else {
03878 Tcl_ListObjAppendElement(interp, axesListObj, Tcl_NewListObj(0, NULL));
03879 }
03880 Tcl_ListObjAppendElement(interp, axesListListObj, axesListObj);
03881 }
03882 Tcl_ListObjAppendElement(interp, tcl_result, axesListListObj);
03883 
03884 
03885 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("rotreflect", -1));
03886 Tcl_Obj *rraxesListListObj = Tcl_NewListObj(0, NULL);
03887 for (i=0; i<s->numrotreflect(); i++) {
03888 Tcl_Obj *axesObj = Tcl_NewListObj(0, NULL);
03889 //printf("Tcl %i: %.2f %.2f %.2f\n", i, s->axis(i)[0], s->axis(i)[1], s->axis(i)[2]);
03890 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(s->rotreflect(i)[0]));
03891 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(s->rotreflect(i)[1]));
03892 Tcl_ListObjAppendElement(interp, axesObj, Tcl_NewDoubleObj(s->rotreflect(i)[2]));
03893 Tcl_Obj *rraxesListObj = Tcl_NewListObj(0, NULL);
03894 Tcl_ListObjAppendElement(interp, rraxesListObj, axesObj);
03895 Tcl_ListObjAppendElement(interp, rraxesListObj, Tcl_NewIntObj(s->get_rotreflectorder(i)));
03896 Tcl_ListObjAppendElement(interp, rraxesListObj, Tcl_NewIntObj(s->get_rotrefltype(i)));
03897 Tcl_ListObjAppendElement(interp, rraxesListListObj, rraxesListObj);
03898 }
03899 Tcl_ListObjAppendElement(interp, tcl_result, rraxesListListObj);
03900 
03901 
03902 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("planes", -1));
03903 Tcl_Obj *planesListListObj = Tcl_NewListObj(0, NULL);
03904 for (i=0; i<s->numplanes(); i++) {
03905 Tcl_Obj *planesObj = Tcl_NewListObj(0, NULL);
03906 //printf("Tcl %i: %.2f %.2f %.2f\n", i, s->plane(i)[0], s->plane(i)[1], s->plane(i)[2]);
03907 Tcl_ListObjAppendElement(interp, planesObj, Tcl_NewDoubleObj(s->plane(i)[0]));
03908 Tcl_ListObjAppendElement(interp, planesObj, Tcl_NewDoubleObj(s->plane(i)[1]));
03909 Tcl_ListObjAppendElement(interp, planesObj, Tcl_NewDoubleObj(s->plane(i)[2]));
03910 Tcl_Obj *planesListObj = Tcl_NewListObj(0, NULL);
03911 Tcl_ListObjAppendElement(interp, planesListObj, planesObj);
03912 switch (s->get_planetype(i)) {
03913 case 1:
03914 Tcl_ListObjAppendElement(interp, planesListObj, Tcl_NewStringObj("vertical", -1));
03915 break;
03916 case 3:
03917 Tcl_ListObjAppendElement(interp, planesListObj, Tcl_NewStringObj("dihedral", -1));
03918 break;
03919 case 4:
03920 Tcl_ListObjAppendElement(interp, planesListObj, Tcl_NewStringObj("horizontal", -1));
03921 break;
03922 default:
03923 Tcl_ListObjAppendElement(interp, planesListObj, Tcl_NewListObj(0, NULL));
03924 }
03925 Tcl_ListObjAppendElement(interp, planesListListObj, planesListObj);
03926 }
03927 Tcl_ListObjAppendElement(interp, tcl_result, planesListListObj);
03928 
03929 
03930 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("ideal", -1));
03931 Tcl_Obj *idealcoorListObj = Tcl_NewListObj(0, NULL);
03932 for (i=0; i<natoms; i++) {
03933 Tcl_Obj *idealcoorObj = Tcl_NewListObj(0, NULL);
03934 //printf("Tcl %i: %.2f %.2f %.2f\n", i, s->plane(i)[0], s->plane(i)[1], s->plane(i)[2]);
03935 Tcl_ListObjAppendElement(interp, idealcoorObj, Tcl_NewDoubleObj(s->idealpos(i)[0]));
03936 Tcl_ListObjAppendElement(interp, idealcoorObj, Tcl_NewDoubleObj(s->idealpos(i)[1]));
03937 Tcl_ListObjAppendElement(interp, idealcoorObj, Tcl_NewDoubleObj(s->idealpos(i)[2]));
03938 Tcl_ListObjAppendElement(interp, idealcoorListObj, idealcoorObj);
03939 }
03940 Tcl_ListObjAppendElement(interp, tcl_result, idealcoorListObj);
03941 
03942 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("unique", -1));
03943 Tcl_Obj *uniquecoorListObj = Tcl_NewListObj(0, NULL);
03944 for (i=0; i<natoms; i++) {
03945 if (!s->get_unique_atom(i)) continue; 
03946 Tcl_ListObjAppendElement(interp, uniquecoorListObj, Tcl_NewIntObj(i));
03947 }
03948 Tcl_ListObjAppendElement(interp, tcl_result, uniquecoorListObj);
03949 
03950 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj("orient", -1));
03951 Matrix4 *orient = s->get_orientation();
03952 Tcl_Obj *matrixObj = Tcl_NewListObj(0, NULL);
03953 for (i=0; i<4; i++) {
03954 Tcl_Obj *rowObj = Tcl_NewListObj(0, NULL);
03955 Tcl_ListObjAppendElement(interp, rowObj, Tcl_NewDoubleObj(orient->mat[i+0]));
03956 Tcl_ListObjAppendElement(interp, rowObj, Tcl_NewDoubleObj(orient->mat[i+4]));
03957 Tcl_ListObjAppendElement(interp, rowObj, Tcl_NewDoubleObj(orient->mat[i+8]));
03958 Tcl_ListObjAppendElement(interp, rowObj, Tcl_NewDoubleObj(orient->mat[i+12]));
03959 Tcl_ListObjAppendElement(interp, matrixObj, rowObj);
03960 }
03961 Tcl_ListObjAppendElement(interp, tcl_result, matrixObj);
03962 
03963 Tcl_SetObjResult(interp, tcl_result);
03964 if (idealsel) {
03965 delete s;
03966 }
03967 return TCL_OK;
03968 
03969 }
03970 
03971 // Function: vmd_measure_transoverlap <selection>
03972 // Returns: The structural overlap of a selection with a copy of itself
03973 // that is transformed according to a given transformation matrix.
03974 // The normalized sum over all gaussian function values of the pair
03975 // distances between atoms in the original and the transformed
03976 // selection.
03977 // Input: the atom selection and the transformation matrix.
03978 // Option: with -sigma you can specify the sigma value of the overlap 
03979 // gaussian function. The default is 0.1 Angstrom.
03980 static int vmd_measure_trans_overlap(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
03981 if (argc!=3 && argc!=5) {
03982 Tcl_WrongNumArgs(interp, 2, objv-1, (char *) "<sel> <matrix> [-sigma <value>]");
03983 return TCL_ERROR;
03984 }
03985 AtomSel *sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
03986 if (!sel) {
03987 Tcl_AppendResult(interp, "measure transoverlap: no atom selection", NULL);
03988 return TCL_ERROR;
03989 }
03990 if (!app->molecule_valid_id(sel->molid())) {
03991 Tcl_AppendResult(interp, "measure transoverlap: ",
03992 measure_error(MEASURE_ERR_NOMOLECULE), NULL);
03993 return TCL_ERROR;
03994 }
03995 
03996 // Get the transformation matrix
03997 Matrix4 trans;
03998 if (tcl_get_matrix("measure transoverlap: ",interp,objv[2], trans.mat) != TCL_OK) {
03999 return TCL_ERROR;
04000 }
04001 
04002 double sigma = 0.1;
04003 if (argc==5) {
04004 if (!strupncmp(Tcl_GetStringFromObj(objv[3],NULL), "-sigma", CMDLEN)) {
04005 if (Tcl_GetDoubleFromObj(interp, objv[4], &sigma) != TCL_OK) {
04006 Tcl_AppendResult(interp, " measure transoverlap: bad sigma value", NULL);
04007 return TCL_ERROR;
04008 }
04009 } else {
04010 Tcl_AppendResult(interp, " measure transoverlap: unrecognized option\n", NULL);
04011 Tcl_AppendResult(interp, " Usage: measure transoverlap <sel> <matrix> [-sigma <value>]", NULL);
04012 return TCL_ERROR;
04013 }
04014 }
04015 
04016 int maxnatoms = 100;// XXX
04017 float overlap;
04018 int ret_val = measure_trans_overlap(sel, app->moleculeList, &trans,
04019 float(sigma), NOSKIP_IDENTICAL, 
04020 maxnatoms, overlap);
04021 if (ret_val < 0) {
04022 Tcl_AppendResult(interp, "measure transoverlap: ", measure_error(ret_val), NULL);
04023 return TCL_ERROR;
04024 }
04025 
04026 Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
04027 Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewDoubleObj(overlap));
04028 Tcl_SetObjResult(interp, tcl_result);
04029 
04030 return TCL_OK;
04031 }
04032 
04033 
04034 
04035 //
04036 // Raycasting brute-force approach by Juan R. Perilla <juan@perilla.me>
04037 //
04038 int vmd_measure_volinterior(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
04039 int verbose = 0;
04040 int recognized=0;
04041 if (argc < 3) {
04042 // " options: --allframes (average over all frames)\n"
04043 Tcl_SetResult(interp, (char *) "usage: volinterior "
04044 " <selection1> [options]\n"
04045 "-res (default 10.0)\n"
04046 "-spacing (default res/3)\n"
04047 "-loadmap (load synth map)\n"
04048 "-probmap (compute and load probability map)\n"
04049 " --> -count_pmap (return percentile-based voxel counts)\n"
04050 "-discretize [float cutoff] (creates discrete map with cutoff probability)\n"
04051 "-mol molid [default: top]\n"
04052 "-vol volid [0]\n"
04053 "-nrays number of rays to cast [6]\n"
04054 "-isovalue [float: 1.0]\n"
04055 "-verbose\n"
04056 "-overwrite volid\n",
04057 TCL_STATIC);
04058 return TCL_ERROR;
04059 }
04060 
04061 //atom selection
04062 AtomSel *sel = NULL;
04063 sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
04064 if (!sel) {
04065 Tcl_AppendResult(interp, "volinterior: no atom selection.", NULL);
04066 return TCL_ERROR;
04067 }
04068 if (!sel->selected) {
04069 Tcl_AppendResult(interp, "volinterior: no atoms selected.", NULL);
04070 return TCL_ERROR;
04071 }
04072 if (!app->molecule_valid_id(sel->molid())) {
04073 Tcl_AppendResult(interp, "invalid mol id.", NULL);
04074 return TCL_ERROR;
04075 }
04076 int loadsynthmap = 0;
04077 int loadpmap=0; 
04078 int count_pmap=0;
04079 int print_raydirs=0; 
04080 int molid = -1;
04081 int volid = -1;
04082 int volidOverwrite = -1;
04083 int overwrite=0;
04084 int nrays=6;
04085 long Nout = 0;
04086 float isovalue=1.0;
04087 float proCutoff=0.90f;
04088 int process=0;
04089 float radscale;
04090 double resolution = 10.0;
04091 double gspacing;
04092 MoleculeList *mlist = app->moleculeList;
04093 Molecule *mymol = mlist->mol_from_id(sel->molid());
04094 
04095 for (int i=2; i < argc; i++) {
04096 recognized=0;
04097 char *opt = Tcl_GetStringFromObj(objv[i], NULL);
04098 if (!strcmp(opt, "volinterior")) {
04099 recognized=1;
04100 }
04101 if (!strcmp(opt, "-mol")) {
04102 if (i == argc-1) {
04103 Tcl_AppendResult(interp, "No molid specified",NULL);
04104 return TCL_ERROR;
04105 } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &molid) != TCL_OK) {
04106 Tcl_AppendResult(interp, "\n molid incorrectly specified",NULL);
04107 return TCL_ERROR;
04108 }
04109 recognized=1;
04110 }
04111 if (!strcmp(opt, "-probmap")) {
04112 loadpmap=1;
04113 recognized=1;
04114 }
04115 if (!strcmp(opt, "-checkrays")) {
04116 print_raydirs=1;
04117 recognized=1;
04118 }
04119 if (!strcmp(opt, "-count_pmap")) {
04120 count_pmap=1;
04121 recognized=1;
04122 }
04123 if (!strcmp(opt, "-discretize")) {
04124 double tcutoff;
04125 if (i==argc-1) {
04126 Tcl_AppendResult(interp, "No cutoff probability given for -discretize flag",NULL);
04127 return TCL_ERROR;
04128 } else if (Tcl_GetDoubleFromObj(interp, objv[i+1], &tcutoff)!=TCL_OK) {
04129 Tcl_AppendResult(interp, "\n Cutoff probability for -discretize flag incorrectly specified.",NULL);
04130 return TCL_ERROR;
04131 }
04132 proCutoff=(float)tcutoff;
04133 recognized=1;
04134 process=1;
04135 }
04136 if (!strcmp(opt, "-vol")) {
04137 if (i == argc-1){
04138 Tcl_AppendResult(interp, "No volume id specified",NULL);
04139 return TCL_ERROR;
04140 } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &volid) != TCL_OK) {
04141 Tcl_AppendResult(interp, "\n volume id incorrectly specified",NULL);
04142 return TCL_ERROR;
04143 }
04144 recognized=1;
04145 }
04146 
04147 if (!strcmp(opt, "-overwrite")) {
04148 if (i == argc-1){
04149 Tcl_AppendResult(interp, "No volume id specified",NULL);
04150 return TCL_ERROR;
04151 } else if (Tcl_GetIntFromObj(interp, objv[i+1], &volidOverwrite) != TCL_OK) {
04152 Tcl_AppendResult(interp, "\n volume id incorrectly specified",NULL);
04153 return TCL_ERROR;
04154 }
04155 overwrite=1;
04156 recognized=1;
04157 }
04158 
04159 // Calculate synthmap
04160 if (!strcmp(opt, "-loadmap")) {
04161 loadsynthmap = 1;
04162 recognized=1;
04163 }
04164 if (!strcmp(opt, "-mol")) {
04165 if (i == argc-1) {
04166 Tcl_AppendResult(interp, "No molid specified",NULL);
04167 return TCL_ERROR;
04168 } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &molid) != TCL_OK) {
04169 Tcl_AppendResult(interp, "\n molid incorrectly specified",NULL);
04170 return TCL_ERROR;
04171 }
04172 recognized=1;
04173 }
04174 if (!strcmp(opt, "-res")) {
04175 if (i == argc-1) {
04176 Tcl_AppendResult(interp, "No resolution specified",NULL);
04177 return TCL_ERROR;
04178 } else if (Tcl_GetDoubleFromObj(interp, objv[i+1], &resolution) != TCL_OK) { 
04179 Tcl_AppendResult(interp, "\nResolution incorrectly specified",NULL);
04180 return TCL_ERROR;
04181 }
04182 recognized=1;
04183 }
04184 if (!strcmp(opt, "-spacing")) {
04185 if (i == argc-1) {
04186 Tcl_AppendResult(interp, "No spacing specified",NULL);
04187 return TCL_ERROR;
04188 } else if ( Tcl_GetDoubleFromObj(interp, objv[i+1], &gspacing) != TCL_OK) {
04189 Tcl_AppendResult(interp, "\nspacing incorrectly specified",NULL);
04190 return TCL_ERROR;
04191 }
04192 recognized=1;
04193 }
04194 
04195 if (!strcmp(opt,"-nrays")) {
04196 if (i == argc-1) {
04197 Tcl_AppendResult(interp, "No number of rays specified",NULL);
04198 return TCL_ERROR;
04199 } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &nrays) != TCL_OK) {
04200 Tcl_AppendResult(interp, "\n number of rays incorrectly specified",NULL);
04201 return TCL_ERROR;
04202 }
04203 recognized=1;
04204 }
04205 
04206 if (!strcmp(opt,"-isovalue")) {
04207 double tmp;
04208 if (i == argc-1) {
04209 Tcl_AppendResult(interp,"No isovalue specified",NULL);
04210 return TCL_ERROR;
04211 } else if (Tcl_GetDoubleFromObj(interp,objv[i+1],&tmp) != TCL_OK) {
04212 Tcl_AppendResult(interp,"\n Isovalue incorrectly specified", NULL);
04213 return TCL_ERROR;
04214 }
04215 isovalue=(float) tmp;
04216 recognized=1;
04217 }
04218 
04219 // if (!strcmp(opt, "-o")) {
04220 // if (i == argc-1) {
04221 // // Tcl_AppendResult(interp, "No output file specified",NULL);
04222 // // return TCL_ERROR;
04223 // if (verbose) 
04224 // printf("No output file specified.");
04225 // } else {
04226 // outputmap = Tcl_GetStringFromObj(objv[i+1], NULL);
04227 // }
04228 // }
04229 
04230 if (!strcmp(opt, "-verbose") || (getenv("VMDVOLINVERBOSE") != NULL)) {
04231 verbose = 1;
04232 recognized=1;
04233 }
04234 
04235 if (recognized==0) {
04236 // XXX what case is this supposed to handle? Dangling number parms?
04237 if (isdigit(opt[0]) || isfloat(opt)) { 
04238 continue;
04239 } else {
04240 Tcl_AppendResult(interp,"\n unrecognized option: '",opt,"'",NULL);
04241 return TCL_ERROR;
04242 }
04243 }
04244 }
04245 
04246 if ((process||count_pmap) && !loadpmap) {
04247 Tcl_AppendResult(interp,"\n Cannot give -process or -count_pmap flags given without -probmap!",NULL);
04248 return TCL_ERROR;
04249 }
04250 
04251 Molecule *volmolOverwrite = NULL;
04252 VolumetricData *volmapOverwrite = NULL;
04253 Molecule *pmolOverwrite=NULL;
04254 VolumetricData *pmapOverwrite=NULL;
04255 if (overwrite==1) {
04256 if (loadpmap!=1) {
04257 volmolOverwrite = mlist->mol_from_id(molid);
04258 volmapOverwrite = volmolOverwrite->modify_volume_data(volidOverwrite);
04259 if(volmapOverwrite==NULL) {
04260 Tcl_AppendResult(interp, "\n no overwrite volume correctly specified",NULL);
04261 return TCL_ERROR;
04262 }
04263 } else {
04264 pmolOverwrite=mlist->mol_from_id(molid);
04265 pmapOverwrite=pmolOverwrite->modify_volume_data(volidOverwrite);
04266 if (pmapOverwrite==NULL) {
04267 Tcl_AppendResult(interp,"\n no pmap overwrite volume correctly specified",NULL);
04268 return TCL_ERROR;
04269 }
04270 }
04271 }
04272 
04273 VolumetricData *pmapT = NULL;
04274 VolumetricData *procmap = NULL;
04275 VolumetricData *pmap = NULL; 
04276 VolumetricData *target = NULL;
04277 VolumetricData *volmapA = NULL;
04278 
04279 
04280 const float *framepos = sel->coordinates(app->moleculeList);
04281 const float *radii = mymol->radius();
04282 radscale= 0.2f * float(resolution);
04283 if (gspacing == 0) {
04284 gspacing=resolution*0.33;
04285 }
04286 
04287 int quality=0;
04288 QuickSurf *qs = new QuickSurf();
04289 if (resolution >= 9)
04290 quality = 0;
04291 else
04292 quality = 3;
04293 
04294 if (molid != -1 && volid != -1 ) {
04295 Molecule *volmol = mlist->mol_from_id(molid);
04296 volmapA = (VolumetricData *) volmol->get_volume_data(volid);
04297 if (volmapA==NULL) {
04298 Tcl_AppendResult(interp, "\n no target volume correctly specified", NULL);
04299 return TCL_ERROR;
04300 }
04301 } else {
04302 if (verbose) printf("\n Calculating grid ... \n");
04303 int cuda_err=-1;
04304 #if defined(VMDCUDA)
04305 cuda_err = vmd_cuda_calc_density(sel, app->moleculeList, quality,
04306 radscale, gspacing, &volmapA, 
04307 NULL, NULL, verbose);
04308 #endif
04309 if (cuda_err == -1 ) {
04310 if (verbose) printf("Using CPU version of the code ... \n");
04311 volmapA = qs->calc_density_map(sel, mymol, framepos, radii,
04312 quality, radscale, float(gspacing));
04313 }
04314 if (verbose) printf("Done.\n");
04315 }
04316 
04317 if (volmapA == NULL) {
04318 Tcl_AppendResult(interp, "\n no test volume or molid correctly specified",NULL);
04319 return TCL_ERROR;
04320 }
04321 target=CreateEmptyGrid(volmapA);
04322 
04323 if (loadpmap == 1) {
04324 pmapT = CreateProbGrid(volmapA);
04325 }
04326 
04327 // Normalize dir vector 
04328 // TODO: Let user define direction vectors
04329 
04330 //vec_normalize(rayDir);
04331 
04332 // Create Grid
04333 #if 0
04334 float rayDir[3] = {1,0,0};
04335 Nout += RaycastGrid(volmapA,target,isovalue,rayDir);
04336 rayDir[0] = -1; rayDir[1]=0; rayDir[2]=0;
04337 Nout += RaycastGrid(volmapA,target,isovalue,rayDir);
04338 rayDir[0] = 0; rayDir[1]=1; rayDir[2]=0;
04339 Nout += RaycastGrid(volmapA,target,isovalue,rayDir);
04340 rayDir[0] = 0; rayDir[1]=-1; rayDir[2]=0;
04341 Nout += RaycastGrid(volmapA,target,isovalue,rayDir);
04342 rayDir[0] = 0; rayDir[1]=0; rayDir[2]=1;
04343 Nout += RaycastGrid(volmapA,target,isovalue,rayDir);
04344 rayDir[0] = 0; rayDir[1]=0; rayDir[2]=-1;
04345 Nout += RaycastGrid(volmapA,target,isovalue,rayDir);
04346 #endif
04347 #if 0
04348 if (verbose) printf("Marking down boundary.\n");
04349 long nVoxIso=markIsoGrid(volmapA, target,isovalue);
04350 if (verbose) printf("Casting rays ...\n");
04351 float rayDir[3] = {0.1f,0.1f,1.0f};
04352 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04353 if (verbose) printf("Dir 1 %ld \n",Nout);
04354 rayDir[0] = -0.1f; rayDir[1]=-0.1f; rayDir[2]=-1.0f;
04355 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04356 if (verbose) printf("Dir 2 %ld \n",Nout);
04357 rayDir[0] = 0.1f; rayDir[1]=1.0f; rayDir[2]=0.1f;
04358 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04359 if (verbose) printf("Dir 3 %ld \n",Nout);
04360 rayDir[0] = -0.1f; rayDir[1]=-1.0f; rayDir[2]=-0.1f;
04361 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04362 if (verbose) printf("Dir 4 %ld \n",Nout);
04363 rayDir[0] = 1.0f; rayDir[1]=0.1f; rayDir[2]=0.1f;
04364 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04365 if (verbose) printf("Dir 5 %ld \n",Nout);
04366 rayDir[0] = -1.0f; rayDir[1]=-0.1f; rayDir[2]=-0.1f;
04367 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04368 if (verbose) printf("Dir 6 %ld \n",Nout);
04369 rayDir[0] = -0.5f; rayDir[1]=-0.5f; rayDir[2]=-0.5f;
04370 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04371 if (verbose) printf("Dir 7 %ld \n",Nout);
04372 rayDir[0] = 0.5f; rayDir[1]=0.5f; rayDir[2]=0.5f;
04373 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04374 if (verbose) printf("Dir 9 %ld \n",Nout);
04375 rayDir[0] = -0.5f; rayDir[1]=0.5f; rayDir[2]=0.5f;
04376 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04377 if (verbose) printf("Dir 10 %ld \n",Nout);
04378 rayDir[0] = 0.5f; rayDir[1]=-0.5f; rayDir[2]=0.5f;
04379 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04380 if (verbose) printf("Dir 11 %ld \n",Nout);
04381 rayDir[0] = 0.5f; rayDir[1]=0.5f; rayDir[2]=-0.5f;
04382 Nout += volin_threaded(volmapA,target,isovalue,rayDir);
04383 if (verbose) printf("Dir 12 %ld \n",Nout);
04384 #endif
04385 if (verbose) printf("Marking down boundary.\n");
04386 long nVoxIso=markIsoGrid(volmapA, target,isovalue);
04387 float rayDir[3] = {0.1f,0.1f,1.0f};
04388 static const float RAND_MAX_INV = 1.0f/float(VMD_RAND_MAX);
04389 
04390 float *poisson_set = new float[(nrays+1)*2];
04391 int numcandidates=40;
04392 int converged=poisson_sample_on_sphere(poisson_set,nrays,numcandidates,verbose);
04393 if (converged==0) {
04394 if (verbose) {
04395 printf("Poisson disk sampler failed for measure volinterior!\n");
04396 printf("... continuing with standard ray direction generator.\n");
04397 }
04398 }
04399 
04400 if (print_raydirs && converged)
04401 print_xyz(poisson_set,nrays);
04402 
04403 if (verbose) printf("Casting rays...\n");
04404 // Seed the random number generator before each calculation. 
04405 //vmd_srandom(103467829);
04406 // Cast nrays uniformly over the sphere
04407 
04408 // XXX Use of sincosf() now in place. It would be much more work-efficient
04409 // to build up batches of rays and pass them into the multithreaded
04410 // work routines, as launching the worker threads has a non-trivial
04411 // cost. It would be better if they plow through many/all rays before
04412 // returning, since the control loops here aren't checking for any
04413 // conditions to be met other than that we've processed all of the rays.
04414 
04415 float sine1, cosine1, sine2, cosine2;
04416 for (int ray=0; ray < nrays; ray++) { 
04417 float u1 = (float) vmd_random();
04418 float u2 = (float) vmd_random();
04419 float z = 2.0f*u1*RAND_MAX_INV -1.0f;
04420 float phi = (float) (2.0f*VMD_PI*u2*RAND_MAX_INV);
04421 float R = sqrtf(1.0f-z*z);
04422 if (converged==1) {
04423 sincosf(poisson_set[ray*2+0],&sine1,&cosine1);
04424 sincosf(poisson_set[ray*2+1],&sine2,&cosine2);
04425 rayDir[0]=cosine1*sine2;
04426 rayDir[1]=sine1*sine2;
04427 rayDir[2]=cosine2;
04428 } else {
04429 sincosf(phi,&sine1,&cosine1);
04430 rayDir[0]=R*cosine1;
04431 rayDir[1]=R*sine1;
04432 rayDir[2]=z;
04433 }
04434 
04435 // XXX rather than launching threads per ray, we should be launching
04436 // the worker threads and processing _all_ of the rays...
04437 // That would also lead toward an implementation that could
04438 // eventually be migrated to the GPU for higher performance.
04439 if (loadpmap == 1) {
04440 Nout += volin_threaded_prob(volmapA, target, pmapT, isovalue, rayDir);
04441 } else {
04442 Nout += volin_threaded(volmapA, target, isovalue, rayDir);
04443 }
04444 
04445 if (verbose) printf("Ray(%d) (%4.3f %4.3f %4.3f). Voxels : %ld \n",ray+1,rayDir[0],rayDir[1],rayDir[2],Nout); 
04446 }
04447 if (verbose) printf("Done.\n");
04448 
04449 delete [] poisson_set;
04450 
04451 if (verbose) printf("Counting voxels above isovalue ...\n");
04452 
04453 Nout=countIsoGrids(target,EXTERIORVOXEL);
04454 long nIn=countIsoGrids(target,INTERIORVOXEL);
04455 
04456 if (loadsynthmap == 1 ) {
04457 app->molecule_add_volumetric(molid, volmapA->name, volmapA->origin, 
04458 volmapA->xaxis, volmapA->yaxis, volmapA->zaxis,
04459 volmapA->xsize, volmapA->ysize, volmapA->zsize,
04460 volmapA->data);
04461 volmapA->compute_volume_gradient();
04462 volmapA->data = NULL; // prevent destruction of density array; 
04463 }
04464 
04465 if (overwrite != 1 ) {
04466 if (loadpmap==1) {
04467 if (verbose) printf("Normalizing probability map ...\n");
04468 pmap = normalize_pmap(pmapT, nrays);
04469 app->molecule_add_volumetric(molid, pmap->name, pmap->origin,
04470 pmap->xaxis, pmap->yaxis, pmap->zaxis,
04471 pmap->xsize, pmap->ysize, pmap->zsize, 
04472 pmap->data);
04473 pmap->compute_volume_gradient();
04474 } else {
04475 app->molecule_add_volumetric(molid, target->name, target->origin, 
04476 target->xaxis, target->yaxis, target->zaxis,
04477 target->xsize, target->ysize, target->zsize,
04478 target->data);
04479 target->compute_volume_gradient();
04480 target->data = NULL; // prevent destruction of density array; 
04481 }
04482 } else {
04483 if (loadpmap==1) {
04484 if (verbose) printf("Normalizing probability map ...\n");
04485 pmap = normalize_pmap(pmapT, nrays);
04486 if (verbose) printf("Overwriting volume...\n");
04487 pmapOverwrite->name=pmap->name;
04488 pmapOverwrite->xsize=pmap->xsize;
04489 pmapOverwrite->ysize=pmap->ysize;
04490 pmapOverwrite->zsize=pmap->zsize;
04491 pmapOverwrite->data=pmap->data;
04492 pmap->data=NULL;
04493 pmapOverwrite->compute_volume_gradient();
04494 if (verbose) printf("Done!\n");
04495 } else {
04496 if (verbose) printf("Overwriting volume ... \n");
04497 volmapOverwrite->name = target -> name;
04498 volmapOverwrite->xsize = target->xsize;
04499 volmapOverwrite->ysize = target->ysize;
04500 volmapOverwrite->zsize = target->zsize;
04501 volmapOverwrite->data = target->data;
04502 target->data=NULL;
04503 volmapOverwrite->compute_volume_gradient();
04504 if (verbose) printf("Done.\n ");
04505 } 
04506 }
04507 
04508 if ((loadpmap&&process)&&pmap->data!=NULL) {
04509 //Tcl_AppendResult(interp,"Processing probability map with cutoff: ",proCutoff,"...",NULL);
04510 procmap=process_pmap(pmap, proCutoff);
04511 app->molecule_add_volumetric(molid, procmap->name, procmap->origin,
04512 procmap->xaxis, procmap->yaxis, procmap->zaxis,
04513 procmap->xsize, procmap->ysize, procmap->zsize, procmap->data);
04514 procmap->compute_volume_gradient();
04515 procmap->data=NULL;
04516 if (verbose) printf("Loaded processed map\n");
04517 } else if ((loadpmap&&process)&&pmapOverwrite->data!=NULL) {
04518 procmap=process_pmap(pmapOverwrite, proCutoff);
04519 app->molecule_add_volumetric(molid, procmap->name, procmap->origin,
04520 procmap->xaxis, procmap->yaxis, procmap->zaxis,
04521 procmap->xsize, procmap->ysize, procmap->zsize, procmap->data);
04522 procmap->compute_volume_gradient();
04523 procmap->data=NULL;
04524 if (verbose) printf("Loaded processed map (they don't overwrite).\n");
04525 }
04526 
04527 if (count_pmap==1 && pmap->data!=NULL) {
04528 if (loadpmap!=1) {
04529 printf("-count_pmap flag specified without -probmap!\n");
04530 return TCL_ERROR;
04531 } else {
04532 long Pvol=0;
04533 float currProb=0.0f;
04534 Tcl_Obj *analysis=Tcl_NewListObj(0,NULL);
04535 for (int c=0; c<10; c++) {
04536 currProb+=0.1f;
04537 if (c==9) currProb=((float)nrays-((float)nrays*0.01f))/(float)nrays;
04538 Pvol=vol_probability(pmap, currProb, float(gspacing));
04539 Tcl_ListObjAppendElement(interp, analysis, Tcl_NewLongObj(Pvol));
04540 }
04541 Tcl_SetObjResult(interp,analysis);
04542 if (molid != -1 && volid != -1 ) {
04543 target->data = NULL; // prevent destruction of density array;
04544 pmap->data = NULL;
04545 } else {
04546 delete pmapT;
04547 delete qs;
04548 delete volmapA;
04549 } 
04550 
04551 return TCL_OK;
04552 }
04553 } else if (count_pmap==1 && overwrite==1 && pmapOverwrite->data!=NULL) {
04554 long Pvol=0;
04555 float currProb=0.0f;
04556 Tcl_Obj *analysis=Tcl_NewListObj(0,NULL);
04557 for (int c=0; c<10; c++) {
04558 currProb+=0.1f;
04559 if (c==9) currProb=((float)nrays-((float)nrays*0.01f))/(float)nrays;
04560 Pvol=vol_probability(pmapOverwrite, currProb, float(gspacing));
04561 Tcl_ListObjAppendElement(interp, analysis, Tcl_NewLongObj(Pvol));
04562 }
04563 Tcl_SetObjResult(interp,analysis);
04564 if (molid != -1 && volid != -1 ) {
04565 target->data = NULL; // prevent destruction of density array;
04566 pmap->data = NULL;
04567 } else {
04568 delete pmapT;
04569 delete qs;
04570 delete volmapA;
04571 }
04572 
04573 return TCL_OK;
04574 }
04575 
04576 if (molid != -1 && volid != -1 ) {
04577 target->data = NULL; // prevent destruction of density array;
04578 procmap->data=NULL;
04579 pmap->data = NULL;
04580 } else {
04581 delete pmapT;
04582 delete qs;
04583 delete volmapA;
04584 }
04585 
04586 long Ntotal = target->xsize * target->ysize * target->zsize;
04587 //printf("VolIn: %ld external voxels.\n",Nout);
04588 
04589 Tcl_Obj *tcl_result = Tcl_NewListObj(0,NULL);
04590 Tcl_ListObjAppendElement(interp,tcl_result,Tcl_NewLongObj(Ntotal));
04591 Tcl_ListObjAppendElement(interp,tcl_result,Tcl_NewLongObj(Nout));
04592 Tcl_ListObjAppendElement(interp,tcl_result,Tcl_NewLongObj(nIn));
04593 Tcl_ListObjAppendElement(interp,tcl_result,Tcl_NewLongObj(nVoxIso));
04594 Tcl_SetObjResult(interp,tcl_result);
04595 return TCL_OK;
04596 }
04597 
04598 int obj_measure(ClientData cd, Tcl_Interp *interp, int argc,
04599 Tcl_Obj * const objv[]) {
04600 
04601 if (argc < 2) {
04602 Tcl_SetResult(interp, 
04603 (char *) "usage: measure <command> [args...]\n"
04604 "\nMeasure Commands:\n"
04605 " avpos <sel> [first <first>] [last <last>] [step <step>] -- average position\n"
04606 " center <sel> [weight <weights>] -- geometrical (or weighted) center\n"
04607 " centerperresidue <sel> [weight <weights>] -- geometrical center for every residue in sel\n"
04608 " cluster <sel> [num <#clusters>] [distfunc <flag>] [cutoff <cutoff>]\n"
04609 " [first <first>] [last <last>] [step <step>] [selupdate <bool>]\n"
04610 " [weight <weights>]\n"
04611 " -- perform a cluster analysis (cluster similar timesteps)\n"
04612 " clustsize <sel> [cutoff <float>] [minsize <num>] [numshared <num>]\n"
04613 " [usepbc <bool>] [storesize <fieldname>] [storenum <fieldname>]\n"
04614 " -- perform a cluster size analysis (find clusters of atoms)\n"
04615 " contacts <cutoff> <sel1> [<sel2>] -- list contacts\n" 
04616 " dipole <sel> [-elementary|-debye] [-geocenter|-masscenter|-origincenter]\n"
04617 " -- dipole moment\n"
04618 " fit <sel1> <sel2> [weight <weights>] [order <index list>]\n"
04619 " -- transformation matrix from selection 1 to 2\n"
04620 " gofr <sel1> <sel2> [delta <value>] [rmax <value>] [usepbc <bool>]\n"
04621 " [selupdate <bool>] [first <first>] [last <last>] [step <step>]\n"
04622 " -- atomic pair distribution function g(r)\n"
04623 " hbonds <cutoff> <angle> <sel1> [<sel2>]\n"
04624 " -- list donors, acceptors, hydrogens involved in hydrogen bonds\n"
04625 " inverse <matrix> -- inverse matrix\n"
04626 " inertia <sel> [-moments] [-eigenvals] -- COM and principle axes of inertia\n"
04627 " minmax <sel> [-withradii] -- bounding box\n"
04628 " rgyr <sel> [weight <weights>] -- radius of gyration\n"
04629 " rmsd <sel1> <sel2> [weight <weights>] -- RMS deviation\n"
04630 " rmsdperresidue <sel1> <sel2> [weight <weights>] -- deviation per residue\n"
04631 " rmsf <sel> [first <first>] [last <last>] [step <step>] -- RMS fluctuation\n"
04632 " rmsfperresidue <sel> [first <first>] [last <last>] [step <step>] -- fluctuation per residue\n"
04633 " sasa <srad> <sel> [-points <varname>] [-restrict <restrictedsel>]\n"
04634 " [-samples <numsamples>] -- solvent-accessible surface area\n"
04635 " sasaperresidue <srad> <sel> [-points <varname>] [-restrict <restrictedsel>]\n"
04636 " [-samples <numsamples>] -- solvent-accessible surface area per residue\n"
04637 " sumweights <sel> weight <weights> -- sum of selected weights\n"
04638 " bond {{<atomid1> [<molid1>]} {<atomid2> [<molid2>]}}\n"
04639 " [molid <default mol>] [frame <frame|all|last> | first <first> last <last>]\n"
04640 " -- bond length between atoms 1 and 2\n"
04641 " angle {{<atomid1> [<molid1>]} {<atomid2> [<molid2>]} {<atomid3> [<molid3>]}}\n"
04642 " [molid <default mol>] [frame <frame|all|last> | first <first> last <last>]\n"
04643 " -- angle between atoms 1-3\n"
04644 " dihed {{<atomid1> [<molid1>]} ... {<atomid4> [<molid4>]}}\n"
04645 " [molid <default mol>] [frame <frame|all|last> | first <first> last <last>]\n"
04646 " -- dihedral angle between atoms 1-4\n"
04647 " imprp {{<atomid1> [<molid1>]} ... {<atomid4> [<molid4>]}}\n"
04648 " [molid <default mol>] [frame <frame|all|last> | first <first> last <last>]\n"
04649 " -- improper angle between atoms 1-4\n"
04650 // FIXME: Complete 'measure energy' usage info here?
04651 " energy bond|angle|dihed|impr|vdw|elec -- compute energy\n"
04652 " surface <sel> <gridsize> <radius> <thickness> -- surface of selection\n"
04653 " pbc2onc <center> [molid <default>] [frame <frame|last>]\n"
04654 " -- transformation matrix to wrap a nonorthogonal PBC unit cell\n"
04655 " into an orthonormal cell\n"
04656 " pbcneighbors <center> <cutoff> [sel <sel>] [align <matrix>] [molid <default>]\n"
04657 " [frame <frame|last>] [boundingbox <PBC|{<mincoord> <maxcoord>}>]\n"
04658 " -- all image atoms that are within cutoff Angstrom of the pbc unit cell\n"
04659 " symmetry <sel> [element [<vector>]] [-tol <value>] [-nobonds] [-verbose <level>]\n"
04660 " transoverlap <sel> <matrix> [-sigma <value>]\n"
04661 " -- overlap of a structure with a transformed copy of itself\n"
04662 " volinterior <sel> [options]\n"
04663 " -- type 'measure volinterior' for usage\n",
04664 TCL_STATIC);
04665 return TCL_ERROR;
04666 }
04667 VMDApp *app = (VMDApp *)cd;
04668 char *argv1 = Tcl_GetStringFromObj(objv[1],NULL);
04669 if (!strupncmp(argv1, "avpos", CMDLEN))
04670 return vmd_measure_avpos(app, argc-1, objv+1, interp);
04671 else if (!strupncmp(argv1, "center", CMDLEN))
04672 return vmd_measure_center(app, argc-1, objv+1, interp);
04673 #if 1
04674 // XXX in test
04675 else if (!strupncmp(argv1, "centerperresidue", CMDLEN))
04676 return vmd_measure_centerperresidue(app, argc-1, objv+1, interp);
04677 #endif
04678 else if (!strupncmp(argv1, "cluster", CMDLEN))
04679 return vmd_measure_cluster(app, argc-1, objv+1, interp);
04680 else if (!strupncmp(argv1, "clustsize", CMDLEN))
04681 return vmd_measure_clustsize(app, argc-1, objv+1, interp);
04682 else if (!strupncmp(argv1, "contacts", CMDLEN))
04683 return vmd_measure_contacts(app, argc-1, objv+1, interp);
04684 else if (!strupncmp(argv1, "dipole", CMDLEN))
04685 return vmd_measure_dipole(app, argc-1, objv+1, interp);
04686 else if (!strupncmp(argv1, "fit", CMDLEN)) 
04687 return vmd_measure_fit(app, argc-1, objv+1, interp);
04688 else if (!strupncmp(argv1, "minmax", CMDLEN))
04689 return vmd_measure_minmax(app, argc-1, objv+1, interp);
04690 else if (!strupncmp(argv1, "gofr", CMDLEN))
04691 return vmd_measure_gofr(app, argc-1, objv+1, interp);
04692 else if (!strupncmp(argv1, "rdf", CMDLEN))
04693 return vmd_measure_rdf(app, argc-1, objv+1, interp);
04694 else if (!strupncmp(argv1, "hbonds", CMDLEN))
04695 return vmd_measure_hbonds(app, argc-1, objv+1, interp);
04696 else if (!strupncmp(argv1, "inverse", CMDLEN)) 
04697 return vmd_measure_inverse(argc-1, objv+1, interp);
04698 else if (!strupncmp(argv1, "inertia", CMDLEN)) 
04699 return vmd_measure_inertia(app, argc-1, objv+1, interp);
04700 else if (!strupncmp(argv1, "rgyr", CMDLEN))
04701 return vmd_measure_rgyr(app, argc-1, objv+1, interp);
04702 else if (!strupncmp(argv1, "rmsd", CMDLEN))
04703 return vmd_measure_rmsd(app, argc-1, objv+1, interp);
04704 #if 1
04705 // XXX in test
04706 else if (!strupncmp(argv1, "rmsdperresidue", CMDLEN))
04707 return vmd_measure_rmsdperresidue(app, argc-1, objv+1, interp);
04708 #endif
04709 #if 1
04710 else if (!strupncmp(argv1, "rmsd_qcp", CMDLEN))
04711 return vmd_measure_rmsd_qcp(app, argc-1, objv+1, interp);
04712 else if (!strupncmp(argv1, "rmsdmat_qcp", CMDLEN))
04713 return vmd_measure_rmsdmat_qcp(app, argc-1, objv+1, interp);
04714 else if (!strupncmp(argv1, "rmsdmat_qcp_ooc", CMDLEN))
04715 return vmd_measure_rmsdmat_qcp_ooc(app, argc-1, objv+1, interp);
04716 #endif
04717 else if (!strupncmp(argv1, "rmsf", CMDLEN))
04718 return vmd_measure_rmsf(app, argc-1, objv+1, interp);
04719 #if 1
04720 // XXX in test
04721 else if (!strupncmp(argv1, "rmsfperresidue", CMDLEN))
04722 return vmd_measure_rmsfperresidue(app, argc-1, objv+1, interp);
04723 #endif
04724 else if (!strupncmp(argv1, "sasa", CMDLEN))
04725 return vmd_measure_sasa(app, argc-1, objv+1, interp);
04726 else if (!strupncmp(argv1, "sasaperresidue", CMDLEN))
04727 return vmd_measure_sasaperresidue(app, argc-1, objv+1, interp);
04728 #if 1
04729 else if (!strupncmp(argv1, "sasalist", CMDLEN))
04730 return vmd_measure_sasalist(app, argc-1, objv+1, interp);
04731 #endif
04732 else if (!strupncmp(argv1, "sumweights", CMDLEN))
04733 return vmd_measure_sumweights(app, argc-1, objv+1, interp);
04734 else if (!strupncmp(argv1, "imprp", CMDLEN))
04735 return vmd_measure_dihed(app, argc-1, objv+1, interp);
04736 else if (!strupncmp(argv1, "dihed", CMDLEN))
04737 return vmd_measure_dihed(app, argc-1, objv+1, interp);
04738 else if (!strupncmp(argv1, "angle", CMDLEN))
04739 return vmd_measure_angle(app, argc-1, objv+1, interp);
04740 else if (!strupncmp(argv1, "bond", CMDLEN))
04741 return vmd_measure_bond(app, argc-1, objv+1, interp);
04742 else if (!strupncmp(argv1, "energy", CMDLEN))
04743 return vmd_measure_energy(app, argc-1, objv+1, interp);
04744 else if (!strupncmp(argv1, "pbc2onc", CMDLEN))
04745 return vmd_measure_pbc2onc_transform(app, argc-1, objv+1, interp);
04746 else if (!strupncmp(argv1, "pbcneighbors", CMDLEN))
04747 return vmd_measure_pbc_neighbors(app, argc-1, objv+1, interp);
04748 else if (!strupncmp(argv1, "surface", CMDLEN))
04749 return vmd_measure_surface(app, argc-1, objv+1, interp);
04750 else if (!strupncmp(argv1, "transoverlap", CMDLEN))
04751 return vmd_measure_trans_overlap(app, argc-1, objv+1, interp);
04752 else if (!strupncmp(argv1, "symmetry", CMDLEN))
04753 return vmd_measure_symmetry(app, argc-1, objv+1, interp);
04754 else if (!strupncmp(argv1, "volinterior", CMDLEN))
04755 return vmd_measure_volinterior(app, argc-1, objv+1, interp);
04756 
04757 Tcl_SetResult(interp, (char *) "Type 'measure' for summary of usage\n", TCL_VOLATILE);
04758 return TCL_OK;
04759 }
04760 
04761 
04762