TclMolInfo.C Source File

00001 /***************************************************************************
00002  *cr                                                                       
00003  *cr            (C) Copyright 1995-2011 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: TclMolInfo.C,v $
00013  *      $Author: akohlmey $     $Locker:  $             $State: Exp $
00014  *      $Revision: 1.93 $       $Date: 2011/12/06 05:04:38 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   This is a helper function for the 'molinfo # get ...' command in
00019  * TclCommands.h .  It is used to access the information known about
00020  * molecules.
00021  *
00022  ***************************************************************************/
00023 
00024 #include <stdlib.h> 
00025 #include "MoleculeList.h"
00026 #include "tcl.h"
00027 #include "Timestep.h"
00028 #include "TclCommands.h" // for my own, external function definitions
00029 #include "Molecule.h"
00030 #include "MaterialList.h"
00031 #include "VMDApp.h"
00032 #include "Inform.h"
00033 #include "QMData.h"
00034 
00037 class ctxt {
00038 public:
00039   int molinfo_get_index;
00040   int molinfo_get_id;
00041   int molinfo_frame_number;
00042   const char *molinfo_get_string;
00043   MoleculeList *mlist;
00044   MaterialList *matlist;
00045   VMDApp *app;
00046  
00047   ctxt() {
00048     molinfo_get_index = -1;
00049     molinfo_get_id = -1;
00050     molinfo_frame_number = -999;
00051     molinfo_get_string = NULL;
00052     mlist = NULL;
00053     matlist = NULL;
00054     app = NULL;
00055   }
00056 };
00057 
00058 static void write_matrix(Tcl_Interp *interp, const float *mat) {
00059   char s[16*25];
00060   sprintf(s, "{%g %g %g %g} {%g %g %g %g} {%g %g %g %g} {%g %g %g %g}", 
00061           mat[0], mat[4],  mat[8], mat[12],                 
00062           mat[1], mat[5],  mat[9], mat[13],                 
00063           mat[2], mat[6], mat[10], mat[14],                
00064           mat[3], mat[7], mat[11], mat[15]);                
00065   Tcl_AppendElement(interp, s);
00066 }
00067 
00068 static int read_matrix(Tcl_Interp *interp, const char *s, float *mat) {
00069   if (sscanf(s, 
00070         " { %f %f %f %f } { %f %f %f %f } { %f %f %f %f } { %f %f %f %f }",
00071         mat+0, mat+4, mat+ 8, mat+12,
00072         mat+1, mat+5, mat+ 9, mat+13,
00073         mat+2, mat+6, mat+10, mat+14,
00074         mat+3, mat+7, mat+11, mat+15) != 16) {
00075     Tcl_AppendResult(interp, "Matrix must contain 16 elements", NULL);
00076     return 0;
00077   }
00078   return 1;
00079 }
00080 
00081 #define generic_molinfo_data(name, func) \
00082 } else if (!strcmp(arg, name)) { \
00083 char buf[20]; sprintf(buf, "%d", func); Tcl_AppendElement(interp, buf);
00084 
00085 #define generic_molinfo_simulation(name, term)  \
00086 } else if (!strcmp(arg, name)) { \
00087 Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00088 if (!ts) Tcl_AppendElement(interp, "0"); \
00089 else { char buf[20]; sprintf(buf, "%f", ts->energy[term]); Tcl_AppendElement(interp, buf); } 
00090 
00091 #define generic_molinfo_pbc(name, term) \
00092 } else if (!strcmp(arg, name)) { \
00093 Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00094 if (!ts) Tcl_AppendElement(interp, "0"); \
00095 else { char buf[20]; sprintf(buf, "%f", ts->term); Tcl_AppendElement(interp, buf); }
00096 
00097 
00098 #define generic_molinfo_wave_int(name, term) \
00099 } else if (!strcmp(arg, name)) { \
00100   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);  \
00101   Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00102   if (ts && ts->qm_timestep) { \
00103     int i; \
00104     char buf[32]; \
00105     for (i=0; i<ts->qm_timestep->get_num_wavef(); i++) { \
00106       sprintf(buf, "%d", ts->qm_timestep->get_##term(i)); \
00107       Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj(buf, -1)); \
00108     } \
00109   } \
00110   Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00111 
00112 #define generic_molinfo_qmts_int(name, term) \
00113 } else if (!strcmp(arg, name)) { \
00114   Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00115   if (!ts || !ts->qm_timestep) \
00116     Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), \
00117                              Tcl_NewIntObj(0)); \
00118   else { \
00119     char buf[20]; sprintf(buf, "%i", ts->qm_timestep->get_##term()); \
00120     Tcl_AppendElement(interp, buf); \
00121   }
00122 
00123 #define generic_molinfo_qmts_arr(type, name, term, n) \
00124 } else if (!strcmp(arg, name)) { \
00125   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);  \
00126   Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00127   if (ts && ts->qm_timestep && ts->qm_timestep->get_##term()) { \
00128     char buf[20]; \
00129     unsigned int i; \
00130     for (i=0; i<(unsigned int)n; i++) { \
00131       sprintf(buf, type, ts->qm_timestep->get_##term()[i]); \
00132       Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj(buf, -1)); \
00133     } \
00134   } \
00135   Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00136 
00137 #define generic_molinfo_qmts_mat(type, name, term, n, m)  \
00138 } else if (!strcmp(arg, name)) { \
00139   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);  \
00140   Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00141   if (ts && ts->qm_timestep && ts->qm_timestep->term) {             \
00142     char buf[20];  \
00143     unsigned int i, j;  \
00144     for (i=0; i<(unsigned int)n; i++) {  \
00145       Tcl_Obj *rowListObj = Tcl_NewListObj(0, NULL);  \
00146       for (j=0; j<(unsigned int)m; j++) {  \
00147         sprintf(buf, type, ts->qm_timestep->term[i*n+j]);  \
00148         Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewStringObj(buf, -1)); \
00149       }  \
00150       Tcl_ListObjAppendElement(interp, tcl_result, rowListObj);  \
00151     }  \
00152   } \
00153   Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00154 
00155 
00156 #define generic_molinfo_qm(type, name, term) \
00157 } else if (!strcmp(arg, name)) {  \
00158   if (!mol->qm_data) \
00159     Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), \
00160                              Tcl_NewIntObj(0)); \
00161   else {  \
00162     char buf[QMDATA_BUFSIZ]; sprintf(buf, type, mol->qm_data->term);  \
00163     Tcl_AppendElement(interp, buf);  \
00164   }
00165 
00166 #define generic_molinfo_qm_string(name, term) \
00167 } else if (!strcmp(arg, name)) {  \
00168   if (mol->qm_data) {  \
00169     char buf[QMDATA_BUFSIZ]; sprintf(buf, "%s", mol->qm_data->term);  \
00170     Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), Tcl_NewStringObj(buf, -1)); \
00171   } else { \
00172     Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), Tcl_NewListObj(0, NULL)); \
00173   }
00174 
00175 #define generic_molinfo_qm_arr(type, name, term, n)  \
00176 } else if (!strcmp(arg, name)) {  \
00177   Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);  \
00178   if (mol->qm_data && mol->qm_data->get_##term()) { \
00179     char buf[20]; \
00180     unsigned int i; \
00181     for (i=0; i<(unsigned int)n; i++) {  \
00182       sprintf(buf, type, mol->qm_data->get_##term()[i]); \
00183       Tcl_ListObjAppendElement(interp, tcl_result, Tcl_NewStringObj(buf, -1)); \
00184     }  \
00185   } \
00186   Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00187 
00188 #define generic_molinfo_qm_mat(type, name, term, n)  \
00189 } else if (!strcmp(arg, name)) {  \
00190   if (mol->qm_data) {  \
00191     if (mol->qm_data->get_##term()) {  \
00192       char buf[20];  \
00193       unsigned int i, j;  \
00194       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);  \
00195       for (i=0; i<(unsigned int)n; i++) {  \
00196         Tcl_Obj *rowListObj = Tcl_NewListObj(0, NULL);  \
00197         for (j=0; j<(unsigned int)n; j++) {  \
00198           sprintf(buf, type, mol->qm_data->get_##term()[i*n+j]);  \
00199           Tcl_ListObjAppendElement(interp, rowListObj, Tcl_NewStringObj(buf, -1)); \
00200         }  \
00201         Tcl_ListObjAppendElement(interp, tcl_result, rowListObj);  \
00202       }  \
00203       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result); \
00204     }  \
00205   }
00206 
00207 
00208 // given the list of strings, return the data
00209 static int molinfo_get(ctxt context, int molid, int argc, const char *argv[],
00210                        Tcl_Interp *interp, int frame_num) {
00211   // does the molecule exist?
00212   context.molinfo_get_id = molid;
00213   context.molinfo_get_index = context.mlist->mol_index_from_id(molid);
00214   if (context.molinfo_get_index == -1) {
00215     char s[16];
00216     sprintf(s, "%d", molid);
00217     Tcl_AppendResult(interp, "molinfo: get: no molecule exists with id ",
00218                      s, NULL);
00219     return TCL_ERROR;
00220   }
00221   Molecule *mol = context.mlist->molecule(context.molinfo_get_index);
00222 
00223   // get the right frame number
00224   switch (frame_num) {
00225     case AtomSel::TS_NOW: 
00226       context.molinfo_frame_number = mol->frame();
00227       break;
00228     case AtomSel::TS_LAST:
00229       context.molinfo_frame_number = mol->numframes()-1;
00230       break;
00231     default:      
00232       context.molinfo_frame_number = frame_num;
00233   }
00234 
00235   for (int term=0; term<argc; term++) {
00236     context.molinfo_get_string = argv[term];
00237 
00238     // skip initial spaces
00239     const char *arg = argv[term];
00240     while (*arg == ' ')
00241       arg++;
00242 
00243     if (!strcmp(arg, "filename")) {
00244       Tcl_Obj *files = Tcl_NewListObj(0, NULL);
00245       for (int i=0; i<mol->num_files(); i++) {
00246         Tcl_ListObjAppendElement(interp, files, 
00247               Tcl_NewStringObj((char *)mol->get_file(i), -1));
00248       }
00249       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), files);
00250     } else if (!strcmp(arg, "index")) {
00251       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00252                                Tcl_NewIntObj(context.molinfo_get_index));
00253     } else if (!strcmp(arg, "id")) {
00254       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00255                                Tcl_NewIntObj(context.molinfo_get_id));
00256     } else if (!strcmp(arg, "numatoms")) {
00257       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00258                                Tcl_NewIntObj(mol->nAtoms));
00259     } else if (!strcmp(arg, "name")) {
00260       Tcl_AppendElement(interp, mol->molname());
00261     } else if (!strcmp(arg, "numreps")) {
00262       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00263                                Tcl_NewIntObj(mol->components()));
00264     } else if (!strcmp(arg, "numframes")) {
00265       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00266                                Tcl_NewIntObj(mol->numframes()));
00267     } else if (!strcmp(arg, "numvolumedata")) {
00268       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00269                                Tcl_NewIntObj(mol->num_volume_data()));
00270     } else if (!strcmp(arg, "last")) {
00271       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00272                                Tcl_NewIntObj(mol->numframes()-1));
00273     } else if (!strcmp(arg, "frame")) {
00274       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), 
00275                                Tcl_NewIntObj(mol->frame()));
00276     } else if (!strcmp(arg, "filespec")) {
00277       Tcl_Obj *specs = Tcl_NewListObj(0, NULL);
00278       for (int i=0; i<mol->num_files(); i++) {
00279         Tcl_ListObjAppendElement(interp, specs, 
00280             Tcl_NewStringObj((char *)mol->get_file_specs(i), -1));
00281       }
00282       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), specs);
00283     } else if (!strcmp(arg, "filetype")) {
00284       Tcl_Obj *types = Tcl_NewListObj(0, NULL);
00285       for (int i=0; i<mol->num_files(); i++) {
00286         Tcl_ListObjAppendElement(interp, types, 
00287             Tcl_NewStringObj((char *)mol->get_type(i), -1));
00288       }
00289       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), types);
00290     } else if (!strcmp(arg, "database")) {
00291       Tcl_Obj *dbs = Tcl_NewListObj(0, NULL);
00292       for (int i=0; i<mol->num_files(); i++) {
00293         Tcl_ListObjAppendElement(interp, dbs, 
00294             Tcl_NewStringObj((char *)mol->get_database(i), -1));
00295       }
00296       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), dbs);
00297     } else if (!strcmp(arg, "accession")) {
00298       Tcl_Obj *acs = Tcl_NewListObj(0, NULL);
00299       for (int i=0; i<mol->num_files(); i++) {
00300         Tcl_ListObjAppendElement(interp, acs, 
00301             Tcl_NewStringObj((char *)mol->get_accession(i), -1));
00302       }
00303       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), acs);
00304     } else if (!strcmp(arg, "remarks")) {
00305       Tcl_Obj *rmk = Tcl_NewListObj(0, NULL);
00306       for (int i=0; i<mol->num_files(); i++) {
00307         Tcl_ListObjAppendElement(interp, rmk, 
00308             Tcl_NewStringObj((char *)mol->get_remarks(i), -1));
00309       }
00310       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), rmk);
00311     } else if (!strcmp(arg, "center")) {
00312       char s[50];
00313       sprintf(s, "%f %f %f", -mol->centt[0], -mol->centt[1], -mol->centt[2]);
00314       Tcl_AppendElement(interp, s);
00315     } else if (!strcmp(arg, "center_matrix")) {
00316       Matrix4 m;
00317       m.translate(mol->centt);
00318       write_matrix(interp, m.mat);
00319     } else if (!strcmp(arg, "rotate_matrix")) {
00320       write_matrix(interp, mol->rotm.mat);
00321     } else if (!strcmp(arg, "scale_matrix")) {
00322       Matrix4 m;
00323       m.scale(mol->scale);
00324       write_matrix(interp, m.mat);
00325     } else if (!strcmp(arg, "global_matrix")) {
00326       Matrix4 m;
00327       m.translate(mol->globt);
00328       write_matrix(interp, m.mat);
00329     } else if (!strcmp(arg, "view_matrix")) {
00330       write_matrix(interp, mol->tm.mat);
00331     } else if (!strncmp(arg, "rep", 3)) {
00332       const char *tmp = arg;
00333       while (*tmp++ != ' ');
00334       while (*tmp++ == ' ');
00335       int repnum = atoi(tmp-1);
00336       if (repnum<0 || repnum >= mol->components()) {
00337         Tcl_AppendResult(interp, arg, " out of range", NULL);
00338         return TCL_ERROR;
00339       }
00340       Tcl_AppendElement(interp, mol->component(repnum)->atomRep->cmdStr);
00341     } else if (!strncmp(arg, "selection", 9)) {
00342       const char *tmp = arg;
00343       while (*tmp++ != ' ');
00344       while (*tmp++ == ' ');
00345       int repnum = atoi(tmp-1);
00346       if (repnum<0 || repnum >= mol->components()) {
00347         Tcl_AppendResult(interp, arg, " out of range", NULL);
00348         return TCL_ERROR;
00349       }
00350       Tcl_AppendElement(interp, mol->component(repnum)->atomSel->cmdStr);
00351     } else if (!strncmp(arg, "color", 5)) {
00352       const char *tmp = arg;
00353       while (*tmp++ != ' ');
00354       while (*tmp++ == ' ');
00355       int repnum = atoi(tmp-1);
00356       if (repnum >= mol->components()) {
00357         Tcl_AppendResult(interp, arg, " out of range", NULL);
00358         return TCL_ERROR;
00359       }
00360       Tcl_AppendElement(interp, mol->component(repnum)->atomColor->cmdStr);
00361     } else if (!strncmp(arg, "material", 8)) {
00362       const char *tmp = arg;
00363       while (*tmp++ != ' ');
00364       while (*tmp++ == ' ');
00365       int repnum = atoi(tmp-1);
00366       if (repnum >= mol->components()) {
00367         Tcl_AppendResult(interp, arg, " out of range", NULL);
00368         return TCL_ERROR;
00369       }
00370       int matind = mol->component(repnum)->curr_material();
00371       Tcl_AppendElement(interp, context.matlist->material_name(matind));
00372 
00373     generic_molinfo_data("active", mol->active)
00374     generic_molinfo_data("drawn", mol->displayed())
00375     generic_molinfo_data("displayed", mol->displayed())
00376     generic_molinfo_data("fixed", mol->fixed())
00377     generic_molinfo_data("top", context.mlist->is_top(context.molinfo_get_index))
00378     generic_molinfo_simulation("bond", TSE_BOND)
00379     generic_molinfo_simulation("angle", TSE_ANGLE)
00380     generic_molinfo_simulation("dihedral", TSE_DIHE)
00381     generic_molinfo_simulation("improper", TSE_IMPR)
00382     generic_molinfo_simulation("vdw", TSE_VDW)
00383     generic_molinfo_simulation("electrostatic", TSE_COUL)
00384     generic_molinfo_simulation("elec", TSE_COUL)
00385     generic_molinfo_simulation("hbond", TSE_HBOND)
00386     generic_molinfo_simulation("kinetic", TSE_KE)
00387     generic_molinfo_simulation("potential", TSE_PE)
00388     generic_molinfo_simulation("temperature", TSE_TEMP)
00389     generic_molinfo_simulation("temp", TSE_TEMP)
00390     generic_molinfo_simulation("energy", TSE_TOTAL)
00391     generic_molinfo_simulation("volume", TSE_VOLUME)
00392     generic_molinfo_simulation("pressure", TSE_PRESSURE)
00393     generic_molinfo_simulation("efield", TSE_EFIELD)
00394     generic_molinfo_simulation("urey_bradley", TSE_UREY_BRADLEY)
00395     generic_molinfo_simulation("molinfo_restraint", TSE_RESTRAINT)
00396 
00397     } else if (!strcmp(arg, "timesteps")) {
00398       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00399       if (!ts) {
00400         Tcl_AppendElement(interp, "0");
00401       } else { 
00402         char buf[20]; 
00403         sprintf(buf, "%d", ts->timesteps);
00404         Tcl_AppendElement(interp, buf);
00405       }
00406 
00407     generic_molinfo_pbc("a", a_length)
00408     generic_molinfo_pbc("b", b_length)
00409     generic_molinfo_pbc("c", c_length)
00410     generic_molinfo_pbc("alpha", alpha)
00411     generic_molinfo_pbc("beta", beta)
00412     generic_molinfo_pbc("gamma", gamma)
00413     generic_molinfo_pbc("physical_time", physical_time)
00414 
00415     // Orbital info
00416     generic_molinfo_qmts_int("numscfiter",  num_scfiter)
00417     generic_molinfo_qmts_int("numwavef", num_wavef);
00418     generic_molinfo_wave_int("numorbitals", num_orbitals);
00419     generic_molinfo_wave_int("multiplicity", multiplicity);
00420     generic_molinfo_qmts_arr("%.12f", "scfenergy", scfenergies, ts->qm_timestep->get_num_scfiter());
00421     generic_molinfo_qmts_arr("%.6f", "gradients", gradients, 3*mol->nAtoms);
00422 
00423     generic_molinfo_qm("%d", "nimag", get_num_imag());
00424     generic_molinfo_qm("%d", "nintcoords", get_num_intcoords());
00425 //generic_molinfo_qm("%d", "ncart", get_num_cartcoords());
00426     generic_molinfo_qm("%d", "numbasis", num_basis);
00427     generic_molinfo_qm("%d", "numshells", get_num_shells());
00428 
00429     generic_molinfo_qm("%d", "num_basis_funcs", num_wave_f);
00430     generic_molinfo_qm("%d", "numelectrons", num_electrons);
00431     generic_molinfo_qm("%d", "totalcharge",  totalcharge);
00432     generic_molinfo_qm("%d", "nproc",   nproc);
00433     generic_molinfo_qm("%d", "memory", memory);
00434     generic_molinfo_qm_string("runtype", get_runtype_string());
00435     generic_molinfo_qm_string("scftype", get_scftype_string());
00436     generic_molinfo_qm_string("basis_string", basis_string);
00437     generic_molinfo_qm_string("runtitle", runtitle);
00438     generic_molinfo_qm_string("geometry", geometry);
00439     generic_molinfo_qm_string("qmversion", version_string);
00440     generic_molinfo_qm_string("qmstatus", get_status_string());
00441 
00442     generic_molinfo_qm_arr("%d", "num_shells_per_atom", num_shells_per_atom, mol->nAtoms);
00443     generic_molinfo_qm_arr("%d", "num_prim_per_shell", num_prim_per_shell, mol->qm_data->get_num_shells());
00444 //    generic_molinfo_qm_arr("%d", "shell_types", shell_types, mol->qm_data->num_basis);
00445 //    generic_molinfo_qm_arr("%d", "angular_momentum", angular_momentum, mol->qm_data->num_wave_f);
00446 //    generic_molinfo_qm_arr("%.6f", "basis_array", basis_array, 2*mol->qm_data->num_basis);
00447 
00448     } else if (!strcmp(arg, "basis")) {
00449       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00450       if (mol->qm_data && mol->qm_data->get_basis()) {             
00451         char buf[20];                       
00452         int i, j, k;
00453         for (i=0; i<mol->nAtoms; i++) {
00454           if (!mol->qm_data->get_basis(i)) break;
00455           Tcl_Obj *atomListObj = Tcl_NewListObj(0, NULL);
00456           for (j=0; j<mol->qm_data->get_basis(i)->numshells; j++) {
00457             const shell_t *shell = mol->qm_data->get_basis(i, j);
00458             Tcl_Obj *shellListObj = Tcl_NewListObj(0, NULL);
00459             Tcl_Obj *exponListObj = Tcl_NewListObj(0, NULL);
00460             Tcl_Obj *coeffListObj = Tcl_NewListObj(0, NULL);
00461             for (k=0; k<shell->numprims; k++) {
00462               sprintf(buf, "%.7f", shell->prim[k].expon);  
00463               Tcl_ListObjAppendElement(interp, exponListObj, Tcl_NewStringObj(buf, -1));
00464               sprintf(buf, "%.12f", shell->prim[k].coeff);  
00465               Tcl_ListObjAppendElement(interp, coeffListObj, Tcl_NewStringObj(buf, -1));
00466             }
00467             sprintf(buf, "%s", mol->qm_data->get_shell_type_str(shell));
00468             Tcl_ListObjAppendElement(interp, shellListObj, Tcl_NewStringObj(buf, -1));
00469             Tcl_ListObjAppendElement(interp, shellListObj, exponListObj);
00470             Tcl_ListObjAppendElement(interp, shellListObj, coeffListObj);
00471             Tcl_ListObjAppendElement(interp, atomListObj, shellListObj);
00472           }
00473           Tcl_ListObjAppendElement(interp, tcl_result, atomListObj);
00474         }
00475       }
00476       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00477 
00478     } else if (!strcmp(arg, "overlap_matrix")) {
00479       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00480       if (ts && ts->qm_timestep &&
00481           mol->qm_data && mol->qm_data->get_basis()) {
00482         Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00483 
00484         float *overlap_matrix = NULL;
00485         float *expandedbasis = NULL;
00486         int *numprims = NULL;
00487         mol->qm_data->expand_basis_array(expandedbasis, numprims);
00488         mol->qm_data->compute_overlap_integrals(ts, expandedbasis,
00489                                                 numprims,
00490                                                 overlap_matrix);
00491         delete [] expandedbasis;
00492         delete [] numprims;
00493         int i, j;
00494         int numwavef = mol->qm_data->num_wave_f;
00495         char buf[20];
00496         for (i=0; i<numwavef; i++) {
00497           for (j=i; j<numwavef; j++) {
00498              sprintf(buf, "%.6f", overlap_matrix[i*numwavef+j]);
00499              Tcl_ListObjAppendElement(interp, tcl_result,
00500                  Tcl_NewStringObj(buf, -1));
00501 
00502            }
00503         }
00504         Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00505         delete [] overlap_matrix;
00506       }
00507 
00508     } else if (!strcmp(arg, "homo")) {
00509       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00510       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00511       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00512           mol->qm_data) {
00513         int iwave;
00514         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00515           Tcl_ListObjAppendElement(interp, tcl_result,
00516              Tcl_NewIntObj(ts->qm_timestep->get_homo(iwave)));
00517         }
00518       }
00519       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00520 
00521     } else if (!strcmp(arg, "numavailorbs")) {
00522       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00523       if (mol->qm_data) {
00524         int i;
00525         for (i=0; i<mol->qm_data->num_wavef_signa; i++) {
00526           Tcl_ListObjAppendElement(interp, tcl_result,
00527              Tcl_NewIntObj(mol->qm_data->get_max_avail_orbitals(i)));
00528         }
00529       }
00530       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00531 
00532     } else if (!strcmp(arg, "wavef_type")) {
00533       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00534       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00535       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00536           mol->qm_data) {
00537         int iwave;
00538         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00539           char *buf;
00540           ts->qm_timestep->get_wavef_typestr(iwave, buf);
00541           Tcl_ListObjAppendElement(interp, tcl_result,
00542                      Tcl_NewStringObj(buf, -1));
00543           delete [] buf;
00544         }
00545       }
00546       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00547 
00548     } else if (!strcmp(arg, "wavef_spin")) {
00549       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00550       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00551           mol->qm_data) {
00552         Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00553         int iwave;
00554         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00555           Tcl_ListObjAppendElement(interp, tcl_result,
00556                  Tcl_NewIntObj(ts->qm_timestep->get_spin(iwave)));
00557         }
00558         Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00559       }
00560 
00561     } else if (!strcmp(arg, "wavef_excitation")) {
00562       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00563       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00564           mol->qm_data) {
00565         Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00566         int iwave;
00567         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00568           Tcl_ListObjAppendElement(interp, tcl_result,
00569                  Tcl_NewIntObj(ts->qm_timestep->get_excitation(iwave)));
00570         }
00571         Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00572       }
00573 
00574     } else if (!strcmp(arg, "wavef_energy")) {
00575       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00576       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00577           mol->qm_data) {
00578         Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00579         int iwave;
00580         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00581           Tcl_ListObjAppendElement(interp, tcl_result,
00582                  Tcl_NewDoubleObj(ts->qm_timestep->get_wave_energy(iwave)));
00583         }
00584         Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00585       }
00586 
00587     } else if (!strcmp(arg, "orbenergies")) {
00588       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00589       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00590       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00591           mol->qm_data) {
00592         char buf[20];                       
00593         int iwave, orb;
00594         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00595           Tcl_Obj *waveListObj = Tcl_NewListObj(0, NULL);
00596           const float *orben = ts->qm_timestep->get_orbitalenergy(iwave);
00597           if (orben) {
00598             int norbitals = ts->qm_timestep->get_num_orbitals(iwave);
00599             for (orb=0; orb<norbitals; orb++) {
00600               sprintf(buf, "%.12f", orben[orb]);
00601               Tcl_ListObjAppendElement(interp, waveListObj, Tcl_NewStringObj(buf, -1));
00602             }
00603           }
00604           Tcl_ListObjAppendElement(interp, tcl_result, waveListObj);
00605         }
00606       }
00607       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00608 
00609     } else if (!strcmp(arg, "orboccupancies")) {
00610       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00611       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00612       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00613           mol->qm_data) {
00614         char buf[20];                       
00615         int iwave, orb;
00616         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00617           Tcl_Obj *waveListObj = Tcl_NewListObj(0, NULL);
00618           const float *occ = ts->qm_timestep->get_occupancies(iwave);
00619           int norbitals = ts->qm_timestep->get_num_orbitals(iwave);
00620           for (orb=0; orb<norbitals; orb++) {
00621             sprintf(buf, "%.12f", occ[orb]);
00622             Tcl_ListObjAppendElement(interp, waveListObj, Tcl_NewStringObj(buf, -1));
00623           }
00624           Tcl_ListObjAppendElement(interp, tcl_result, waveListObj);
00625         }
00626       }
00627       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00628 
00629     } else if (!strcmp(arg, "wavefunction")) {
00630       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00631       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00632       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00633           mol->qm_data) {
00634         char buf[20];                       
00635         int iwave, orb, i;
00636         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00637           Tcl_Obj *waveListObj = Tcl_NewListObj(0, NULL);
00638           int norbitals = ts->qm_timestep->get_num_orbitals(iwave);
00639           for (orb=0; orb<norbitals; orb++) {
00640             Tcl_Obj *orbListObj = Tcl_NewListObj(0, NULL);
00641             const float *wave_f = ts->qm_timestep->get_wavecoeffs(iwave) + 
00642                                   orb*norbitals;
00643             for (i=0; i<ts->qm_timestep->get_num_coeffs(iwave); i++) {
00644               sprintf(buf, "%.12f", wave_f[i]);
00645               Tcl_ListObjAppendElement(interp, orbListObj, Tcl_NewStringObj(buf, -1));
00646             }
00647             Tcl_ListObjAppendElement(interp, waveListObj, orbListObj);
00648           }
00649           Tcl_ListObjAppendElement(interp, tcl_result, waveListObj);
00650         }
00651       }
00652       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00653        
00654     } else if (!strcmp(arg, "wavef_tree")) {
00655       Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00656       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00657       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_wavef() &&
00658           mol->qm_data && mol->qm_data->get_basis()) {
00659         char buf[20];                       
00660         int iwave, orb, i, j, k;
00661         for (iwave=0; iwave<ts->qm_timestep->get_num_wavef(); iwave++) {
00662           Tcl_Obj *waveListObj = Tcl_NewListObj(0, NULL);
00663           int norbitals = ts->qm_timestep->get_num_orbitals(iwave);
00664           for (orb=0; orb<norbitals; orb++) {
00665             Tcl_Obj *orbListObj = Tcl_NewListObj(0, NULL);
00666             for (i=0; i<mol->nAtoms; i++) {
00667               Tcl_Obj *atomListObj = Tcl_NewListObj(0, NULL);
00668               for (j=0; j<mol->qm_data->get_basis(i)->numshells; j++) {
00669                 const shell_t *shell = mol->qm_data->get_basis(i, j);
00670                 Tcl_Obj *shellListObj = Tcl_NewListObj(0, NULL);
00671                 Tcl_Obj *angListObj  = Tcl_NewListObj(0, NULL);
00672                 Tcl_Obj *waveListObj = Tcl_NewListObj(0, NULL);
00673                 const float *wave_f = ts->qm_timestep->get_wavecoeffs(iwave) + 
00674                                       orb*norbitals +
00675                                       mol->qm_data->get_wave_offset(i, j);
00676                 for (k=0; k<shell->num_cart_func; k++) {
00677                   char *s = mol->qm_data->get_angular_momentum_str(i, j, k);
00678                   Tcl_ListObjAppendElement(interp, angListObj, Tcl_NewStringObj(s, -1));
00679                   delete [] s;
00680                   sprintf(buf, "%.12f", wave_f[k]);
00681                   Tcl_ListObjAppendElement(interp, waveListObj, Tcl_NewStringObj(buf, -1));
00682                 }
00683                 sprintf(buf, "%s", mol->qm_data->get_shell_type_str(shell));
00684                 Tcl_ListObjAppendElement(interp, shellListObj, Tcl_NewStringObj(buf, -1));
00685                 Tcl_ListObjAppendElement(interp, shellListObj, angListObj);
00686                 Tcl_ListObjAppendElement(interp, shellListObj, waveListObj);
00687                 Tcl_ListObjAppendElement(interp, atomListObj, shellListObj);
00688               }
00689               Tcl_ListObjAppendElement(interp, orbListObj, atomListObj);
00690             }
00691             Tcl_ListObjAppendElement(interp, waveListObj, orbListObj);
00692           }
00693           Tcl_ListObjAppendElement(interp, tcl_result, waveListObj);
00694         }
00695       }
00696       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00697         
00698     } else if (!strcmp(arg, "qmcharges")) {
00699       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
00700       if (ts && ts->qm_timestep && ts->qm_timestep->get_num_charge_sets()) {
00701         Tcl_Obj *tcl_result = Tcl_NewListObj(0, NULL);
00702         int i, j;
00703         for (i=0; i<ts->qm_timestep->get_num_charge_sets(); i++) {
00704           Tcl_Obj *chargesetObj = Tcl_NewListObj(0, NULL);
00705           Tcl_Obj *atomListObj = Tcl_NewListObj(0, NULL);
00706           const double *chargeset = ts->qm_timestep->get_charge_set(i);
00707           for (j=0; j<mol->nAtoms; j++) {
00708             Tcl_ListObjAppendElement(interp, atomListObj,
00709                Tcl_NewDoubleObj(chargeset[j]));
00710           }
00711           Tcl_ListObjAppendElement(interp, chargesetObj,
00712             Tcl_NewStringObj(ts->qm_timestep->get_charge_type_str(i), -1));
00713           Tcl_ListObjAppendElement(interp, chargesetObj, atomListObj);
00714           Tcl_ListObjAppendElement(interp, tcl_result, chargesetObj);
00715         }
00716         Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), tcl_result);
00717       }
00718     generic_molinfo_qm_mat("%.6f", "carthessian", carthessian, 3*mol->nAtoms);
00719     generic_molinfo_qm_mat("%.6f", "inthessian",  inthessian,  mol->qm_data->get_num_intcoords());
00720     generic_molinfo_qm_mat("%.8f", "normalmodes", normalmodes, 3*mol->nAtoms);
00721     generic_molinfo_qm_arr("%.2f", "wavenumbers", wavenumbers, 3*mol->nAtoms);
00722     generic_molinfo_qm_arr("%.6f", "intensities", intensities, 3*mol->nAtoms);
00723     generic_molinfo_qm_arr("%d",   "imagmodes",   imagmodes,   mol->qm_data->get_num_imag());
00724 
00725     } else if (!strcmp(arg, "angles")) {
00726       Tcl_Obj *alist = Tcl_NewListObj(0, NULL);
00727       for (int i=0; i<mol->num_angles(); i++) {
00728         Tcl_Obj *adata = Tcl_NewListObj(0,NULL);
00729         int type = -1;
00730         const char *atname;
00731                     
00732         if (mol->angleTypes.num() > 0) {
00733            type = mol->get_angletype(i);
00734         }
00735 
00736         if (type < 0)
00737           atname = "unknown";
00738         else
00739           atname = mol->angleTypeNames.name(type);
00740         
00741         Tcl_ListObjAppendElement(interp, adata, 
00742                                  Tcl_NewStringObj((char *)atname, -1));
00743         Tcl_ListObjAppendElement(interp, adata,
00744                                  Tcl_NewIntObj(mol->angles[3*i]));
00745         Tcl_ListObjAppendElement(interp, adata,
00746                                  Tcl_NewIntObj(mol->angles[3*i+1]));
00747         Tcl_ListObjAppendElement(interp, adata,
00748                                  Tcl_NewIntObj(mol->angles[3*i+2]));
00749         Tcl_ListObjAppendElement(interp, alist, adata);
00750       }
00751       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), alist);
00752     } else if (!strcmp(arg, "dihedrals")) {
00753       Tcl_Obj *alist = Tcl_NewListObj(0, NULL);
00754       for (int i=0; i<mol->num_dihedrals(); i++) {
00755         Tcl_Obj *adata = Tcl_NewListObj(0,NULL);
00756         int type = -1;
00757         const char *atname;
00758                     
00759         if (mol->dihedralTypes.num() > 0) {
00760           type = mol->get_dihedraltype(i);
00761         }
00762 
00763         if (type < 0)
00764           atname = "unknown";
00765         else
00766           atname = mol->dihedralTypeNames.name(type);
00767           
00768         Tcl_ListObjAppendElement(interp, adata, 
00769                                  Tcl_NewStringObj((char *)atname, -1));
00770         Tcl_ListObjAppendElement(interp, adata,
00771                                  Tcl_NewIntObj(mol->dihedrals[4*i]));
00772         Tcl_ListObjAppendElement(interp, adata,
00773                                  Tcl_NewIntObj(mol->dihedrals[4*i+1]));
00774         Tcl_ListObjAppendElement(interp, adata,
00775                                  Tcl_NewIntObj(mol->dihedrals[4*i+2]));
00776         Tcl_ListObjAppendElement(interp, adata,
00777                                  Tcl_NewIntObj(mol->dihedrals[4*i+3]));
00778         Tcl_ListObjAppendElement(interp, alist, adata);
00779       }
00780       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), alist);
00781     } else if (!strcmp(arg, "impropers")) {
00782       Tcl_Obj *alist = Tcl_NewListObj(0, NULL);
00783       for (int i=0; i<mol->num_impropers(); i++) {
00784         Tcl_Obj *adata = Tcl_NewListObj(0,NULL);
00785         int type = -1;
00786         const char *atname;
00787                     
00788         if (mol->improperTypes.num() > 0) {
00789           type = mol->get_impropertype(i);
00790         }
00791 
00792         if (type < 0)
00793           atname = "unknown";
00794         else
00795           atname = mol->improperTypeNames.name(type);
00796           
00797         Tcl_ListObjAppendElement(interp, adata, 
00798                                  Tcl_NewStringObj((char *)atname, -1));
00799         Tcl_ListObjAppendElement(interp, adata,
00800                                  Tcl_NewIntObj(mol->impropers[4*i]));
00801         Tcl_ListObjAppendElement(interp, adata,
00802                                  Tcl_NewIntObj(mol->impropers[4*i+1]));
00803         Tcl_ListObjAppendElement(interp, adata,
00804                                  Tcl_NewIntObj(mol->impropers[4*i+2]));
00805         Tcl_ListObjAppendElement(interp, adata,
00806                                  Tcl_NewIntObj(mol->impropers[4*i+3]));
00807         Tcl_ListObjAppendElement(interp, alist, adata);
00808       }
00809       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), alist);
00810     } else {
00811       Tcl_ResetResult(interp);
00812       Tcl_AppendResult(interp, "molinfo: cannot find molinfo attribute '",
00813                        argv[term], "'", NULL);
00814       return TCL_ERROR;
00815     }
00816   }
00817 
00818   return TCL_OK;
00819 }
00820 
00821 
00822 #define generic_molinfo_set_data(name, func1, func2) \
00823 } else if (!strcmp(argv[term], name)) { \
00824 int onoff; \
00825 if (Tcl_GetBoolean(interp, data[term], &onoff) != TCL_OK) return TCL_ERROR; \
00826 if (onoff) { func1 ; } else { func2 ; }
00827 
00828 #define generic_molinfo_simulation_set(name, type)  \
00829 } else if (!strcmp(argv[term], name)) { \
00830 Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00831 if (ts) ts->energy[type] = (float) atof(data[term]);
00832 
00833 #define generic_molinfo_pbc_set(name, type) \
00834 } else if (!strcmp(argv[term], name)) { \
00835 Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00836 if (ts) { ts->type = (float) atof(data[term]); mol->change_pbc(); }
00837 
00838 // given the list of strings and data, set the right values
00839 static int molinfo_set(ctxt context, int molid, int argc, const char *argv[],
00840                        const char *data[], Tcl_Interp *interp, int frame_num) {
00841   // does the molecule exist?
00842   context.molinfo_get_id = molid;
00843   context.molinfo_get_index = context.mlist->mol_index_from_id(molid);
00844 
00845   if (context.molinfo_get_index == -1) {
00846     char s[10];
00847     sprintf(s, "%d", molid);
00848     Tcl_AppendResult(interp, "molinfo: set: no molecule exists with id ",
00849                      s, NULL);
00850     return TCL_ERROR;
00851   }
00852 
00853   // get the right frame number
00854   switch (frame_num) {
00855     case AtomSel::TS_NOW: 
00856       context.molinfo_frame_number =
00857         context.mlist->molecule(context.molinfo_get_index) -> frame();
00858       break;
00859     case AtomSel::TS_LAST:
00860       context.molinfo_frame_number = 
00861         context.mlist->molecule(context.molinfo_get_index) -> numframes()-1;
00862       break;
00863     default:      
00864       context.molinfo_frame_number = frame_num;
00865   }
00866   Molecule *mol = context.mlist->molecule(context.molinfo_get_index);
00867   
00868   for (int term=0; term<argc; term++) {
00869     context.molinfo_get_string = argv[term];
00870     if (!strcmp(argv[term], "center")) {
00871       float x, y, z;
00872       if (sscanf((const char *)data[term], "%f %f %f", &x, &y, &z) != 3) {
00873         Tcl_AppendResult(interp, 
00874             "molinfo: set center: must have three position elements", NULL);
00875         return TCL_ERROR;
00876       }
00877       mol->change_center(x, y, z);
00878     } else if (!strcmp(argv[term], "center_matrix")) {
00879       float mat[16];
00880       if (!read_matrix(interp, data[term], mat)) return TCL_ERROR;
00881       mol->set_cent_trans(mat[12], mat[13], mat[14]);
00882     } else if (!strcmp(argv[term], "rotate_matrix")) {
00883       Matrix4 mat;
00884       if (!read_matrix(interp, data[term], mat.mat)) return TCL_ERROR;
00885       mol->set_rot(mat);
00886     } else if (!strcmp(argv[term], "scale_matrix")) {
00887       float mat[16];
00888       if (!read_matrix(interp, data[term], mat)) return TCL_ERROR;
00889       mol->set_scale(mat[0]);
00890     } else if (!strcmp(argv[term], "global_matrix")) {
00891       float mat[16];
00892       if (!read_matrix(interp, data[term], mat)) return TCL_ERROR;
00893       mol->set_glob_trans(mat[12], mat[13], mat[14]);
00894     } else if (!strcmp(argv[term], "angles")) {
00895       mol->clear_angles();
00896       mol->angleTypeNames.clear();
00897       mol->set_dataset_flag(BaseMolecule::ANGLES);
00898 
00899       Tcl_Obj *alist = Tcl_NewListObj(0,NULL);
00900       Tcl_SetStringObj(alist, data[term], -1);
00901       Tcl_IncrRefCount(alist); 
00902 
00903       int numangles;
00904       Tcl_Obj **adata;
00905       Tcl_ListObjGetElements(interp, alist, &numangles, &adata);
00906 
00907       for (int i=0; i < numangles; i++) {
00908         int numentries;
00909         Tcl_Obj **a;
00910 
00911         Tcl_ListObjGetElements(interp, adata[i], &numentries, &a);
00912         if (numentries != 4) {
00913           Tcl_AppendResult(interp, "molinfo: incorrect data item for "
00914                            "'set angles' :", Tcl_GetString(adata[i]), NULL);
00915           return TCL_ERROR;
00916         }
00917         int type = mol->angleTypeNames.add_name(Tcl_GetString(a[0]), 0);
00918         int a1, a2, a3;
00919         if (Tcl_GetIntFromObj(interp, a[1], &a1) != TCL_OK) return TCL_ERROR;
00920         if (Tcl_GetIntFromObj(interp, a[2], &a2) != TCL_OK) return TCL_ERROR;
00921         if (Tcl_GetIntFromObj(interp, a[3], &a3) != TCL_OK) return TCL_ERROR;
00922         mol->add_angle(a1, a2, a3, type);
00923       }
00924       // release storage
00925       Tcl_InvalidateStringRep(alist);
00926       Tcl_DecrRefCount(alist); 
00927 
00928     } else if (!strcmp(argv[term], "dihedrals")) { 
00929       mol->clear_dihedrals();
00930       mol->dihedralTypeNames.clear();
00931       mol->set_dataset_flag(BaseMolecule::ANGLES);
00932 
00933       Tcl_Obj *alist = Tcl_NewListObj(0,NULL);
00934       Tcl_SetStringObj(alist, data[term], -1);
00935       Tcl_IncrRefCount(alist); 
00936 
00937       int numdihedrals;
00938       Tcl_Obj **adata;
00939       Tcl_ListObjGetElements(interp, alist, &numdihedrals, &adata);
00940       for (int i=0; i < numdihedrals; i++) {
00941         int numentries;
00942         Tcl_Obj **a;
00943         Tcl_ListObjGetElements(interp, adata[i], &numentries, &a);
00944         if (numentries != 5) {
00945           Tcl_AppendResult(interp, "molinfo: incorrect data item for "
00946                            "'set dihedrals' :", Tcl_GetString(adata[i]), NULL);
00947           return TCL_ERROR;
00948         }
00949         int type = mol->dihedralTypeNames.add_name(Tcl_GetString(a[0]), 0);
00950         int a1, a2, a3, a4;
00951         if (Tcl_GetIntFromObj(interp, a[1], &a1) != TCL_OK) return TCL_ERROR;
00952         if (Tcl_GetIntFromObj(interp, a[2], &a2) != TCL_OK) return TCL_ERROR;
00953         if (Tcl_GetIntFromObj(interp, a[3], &a3) != TCL_OK) return TCL_ERROR;
00954         if (Tcl_GetIntFromObj(interp, a[4], &a4) != TCL_OK) return TCL_ERROR;
00955         mol->add_dihedral(a1, a2, a3, a4, type);
00956       }
00957       // release storage
00958       Tcl_InvalidateStringRep(alist);
00959       Tcl_DecrRefCount(alist); 
00960 
00961     } else if (!strcmp(argv[term], "impropers")) {
00962       mol->clear_impropers();
00963       mol->improperTypeNames.clear();
00964       mol->set_dataset_flag(BaseMolecule::ANGLES);
00965 
00966       Tcl_Obj *alist = Tcl_NewListObj(0,NULL);
00967       Tcl_SetStringObj(alist, data[term], -1);
00968       Tcl_IncrRefCount(alist); 
00969 
00970       int numimpropers;
00971       Tcl_Obj **adata;
00972       Tcl_ListObjGetElements(interp, alist, &numimpropers, &adata);
00973       for (int i=0; i < numimpropers; i++) {
00974         int numentries;
00975         Tcl_Obj **a;
00976         Tcl_ListObjGetElements(interp, adata[i], &numentries, &a);
00977         if (numentries != 5) {
00978           Tcl_AppendResult(interp, "molinfo: incorrect data item for "
00979                            "'set impropers' :", Tcl_GetString(adata[i]), NULL);
00980           return TCL_ERROR;
00981         }
00982         int type = mol->improperTypeNames.add_name(Tcl_GetString(a[0]), 0);
00983         int a1, a2, a3, a4;
00984         if (Tcl_GetIntFromObj(interp, a[1], &a1) != TCL_OK) return TCL_ERROR;
00985         if (Tcl_GetIntFromObj(interp, a[2], &a2) != TCL_OK) return TCL_ERROR;
00986         if (Tcl_GetIntFromObj(interp, a[3], &a3) != TCL_OK) return TCL_ERROR;
00987         if (Tcl_GetIntFromObj(interp, a[4], &a4) != TCL_OK) return TCL_ERROR;
00988         mol->add_improper(a1, a2, a3, a4, type);
00989       }
00990       // release storage
00991       Tcl_InvalidateStringRep(alist);
00992       Tcl_DecrRefCount(alist); 
00993 
00994     } else if (!strcmp(argv[term], "frame")) {
00995       // XXX this isn't correctly sending events to the GUI,
00996       // so if you do 'molinfo top set frame 5' the GUI is 
00997       // out of date
00998       //
00999       // This is correct and intended behavior; if you want GUI events,
01000       // the command to use is "animate goto <frame>. -- JRG 12/27/07
01001       mol->override_current_frame(atoi(data[term]));
01002       mol->change_ts();
01003       
01004     generic_molinfo_set_data("active", context.app->molecule_activate(mol->id(), 1), context.app->molecule_activate(mol->id(), 0))
01005     generic_molinfo_set_data("drawn", context.app->molecule_display(mol->id(), 1), context.app->molecule_display(mol->id(), 0))
01006     generic_molinfo_set_data("displayed", context.app->molecule_display(mol->id(), 1), context.app->molecule_display(mol->id(), 0))
01007     generic_molinfo_set_data("fixed", context.app->molecule_fix(mol->id(), 1), context.app->molecule_fix(mol->id(), 0))
01008     generic_molinfo_set_data("top", context.app->molecule_make_top(mol->id()), Tcl_SetResult(interp, (char *) "Cannot set 'top' to false.", TCL_STATIC); return TCL_ERROR)
01009 
01010     generic_molinfo_simulation_set("bond", TSE_BOND)
01011     generic_molinfo_simulation_set("angle", TSE_ANGLE)
01012     generic_molinfo_simulation_set("dihedral", TSE_DIHE)
01013     generic_molinfo_simulation_set("improper", TSE_IMPR)
01014     generic_molinfo_simulation_set("vdw", TSE_VDW)
01015     generic_molinfo_simulation_set("electrostatic", TSE_COUL)
01016     generic_molinfo_simulation_set("elec", TSE_COUL)
01017     generic_molinfo_simulation_set("hbond", TSE_HBOND)
01018     generic_molinfo_simulation_set("kinetic", TSE_KE)
01019     generic_molinfo_simulation_set("potential", TSE_PE)
01020     generic_molinfo_simulation_set("temperature", TSE_TEMP)
01021     generic_molinfo_simulation_set("temp", TSE_TEMP)
01022     generic_molinfo_simulation_set("energy", TSE_TOTAL)
01023     generic_molinfo_simulation_set("volume", TSE_VOLUME)
01024     generic_molinfo_simulation_set("pressure", TSE_PRESSURE)
01025     generic_molinfo_simulation_set("efield", TSE_EFIELD)
01026     generic_molinfo_simulation_set("urey_bradley", TSE_UREY_BRADLEY)
01027     generic_molinfo_simulation_set("molinfo_restraint", TSE_RESTRAINT)
01028 
01029     } else if (!strcmp(argv[term], "timesteps")) {
01030       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
01031       if (ts) ts->timesteps = atoi(data[term]);
01032 
01033     generic_molinfo_pbc_set("a", a_length)
01034     generic_molinfo_pbc_set("b", b_length)
01035     generic_molinfo_pbc_set("c", c_length)
01036     generic_molinfo_pbc_set("alpha", alpha)
01037     generic_molinfo_pbc_set("beta", beta)
01038     generic_molinfo_pbc_set("gamma", gamma)
01039     generic_molinfo_pbc_set("physical_time", physical_time)
01040 
01041     } else {
01042       Tcl_AppendResult(interp, "molinfo: cannot find molinfo attribute '",
01043                        argv[term], "'", NULL);
01044       return TCL_ERROR;
01045     }
01046   }
01047 
01048   return TCL_OK;
01049 }
01050 
01051 
01052 // Function:  molinfo
01053 // Option  :  molinfo num
01054 //  Returns:   number of molecules
01055 // Option  :  molinfo index <int>
01056 //  Returns:   molecule id of the nth molecule (starting at index = 0)
01057 // Option  :  molinfo list
01058 //  Returns:   list of all molecule ids
01059 // Option  :  molinfo top
01060 //  Returns:   molecule id of the 'top' molecule
01061 // Option  :  molinfo {molecule number} get <data>
01062 //  Returns:   the given data for that molecule
01063 // Option  :  molinfo {molecule number} set <data> <data fields>
01064 //  Does (okay, this isn't a 'info' thing): sets the data field(s)
01065 
01066 int molecule_tcl(ClientData data, Tcl_Interp *interp, int argc, const char *argv[])
01067 {
01068   VMDApp *app = (VMDApp *)data;
01069   // set context variable here
01070   ctxt context;
01071   context.mlist = app->moleculeList;
01072   context.matlist = app->materialList;
01073   context.app = app;
01074 
01075   if (argc == 1) {
01076     Tcl_SetResult(interp,
01077       (char *)
01078       "usage: molinfo <command> [args...]\n\n"
01079       "Commands:"
01080       "\nMolecule IDs:\n"
01081       "  list                  -- lists all existing molecule IDs\n"
01082       "  num                   -- number of loaded molecules\n"
01083       "  top                   -- gets ID of top molecule (or -1 if none)\n"
01084       "  index <n>             -- gets ID of n-th molecule\n"
01085       "\nGetting and Setting Molecular Information:\n" 
01086      // "  keywords              -- returns a list of molinfo keywords\n"     // XXX obsolete???
01087       "  <molid> get <(list of) keywords>\n"
01088       "  <molid> set <(list of) keywords> <(list of) values>\n",
01089       TCL_STATIC);
01090     return TCL_ERROR;
01091   }
01092 
01093   // what does it want?
01094   if (argc == 2) {
01095 
01096 // Option  :  molinfo num
01097 //  Returns:   number of molecules
01098     if (!strcmp(argv[1], "num")) {
01099       char tmpstring[64];
01100       sprintf(tmpstring, "%d", context.mlist->num());
01101       Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
01102       return TCL_OK;
01103     }
01104 
01105 // Option  :  molinfo list
01106 //  Returns:   list of all molecule ids
01107     if (!strcmp(argv[1], "list")) {
01108       if (context.mlist->num() <= 0) {
01109         return TCL_OK;
01110       }
01111       char s[20];
01112       sprintf(s, "%d", context.mlist->molecule(0)->id());
01113       Tcl_AppendResult(interp, s, (char *) NULL);
01114       for (int i=1; i<context.mlist -> num(); i++) {
01115         sprintf(s, "%d", context.mlist->molecule(i)->id());
01116         Tcl_AppendResult(interp, " ", s, (char *) NULL);
01117       }
01118       return TCL_OK;
01119     }
01120 // Option  :  molinfo top
01121 //  Returns:   molecule id of the 'top' molecule
01122     if (!strcmp(argv[1], "top")) {
01123       if (context.mlist->top()) {
01124         char tmpstring[64];
01125         sprintf(tmpstring, "%d", context.mlist->top()->id());
01126         Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
01127       } else {
01128         Tcl_SetResult(interp, (char *) "-1", TCL_STATIC);
01129       }
01130       return TCL_OK;
01131     }
01132 
01133     // otherwise, I don't know
01134     Tcl_AppendResult(interp, "molinfo: couldn't understand '",
01135                      argv[1], "'", NULL);
01136     return TCL_ERROR;
01137   } // end of commands with only one option
01138 
01139   if (argc == 3) { // commands with two options
01140     int val;
01141     if (Tcl_GetInt(interp, argv[2], &val) != TCL_OK) {
01142       return TCL_ERROR;
01143     }
01144 // Option  :  molecule index <int>
01145 //  Returns:   molecule id of the nth molecule (starting at index = 0)
01146     if (!strcmp(argv[1], "index")) {
01147       if (context.mlist->molecule(val)) {
01148         char tmpstring[64];
01149         sprintf(tmpstring, "%d", context.mlist->molecule(val)->id());
01150         Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
01151       } else {
01152         Tcl_SetResult(interp, (char *) "-1", TCL_STATIC);
01153       }
01154       return TCL_OK;
01155     }
01156     Tcl_AppendResult(interp, "molinfo: couldn't understand '",
01157                      argv[1], "'", NULL);
01158     return TCL_ERROR;
01159   }
01160 // Option  :  molinfo {molecule number} get <data> [frame <number>]
01161 //  Returns:   the given data for that molecule
01162   if ((argc == 4 && !strcmp(argv[2], "get")) ||
01163       (argc == 6 && !strcmp(argv[2], "get") && !strcmp(argv[4], "frame"))) {
01164     int frame_num;
01165     if (argc == 4) {
01166       frame_num = AtomSel::TS_NOW;
01167     } else {
01168       if (AtomSel::get_frame_value(argv[5], &frame_num) != 0) {
01169         Tcl_SetResult(interp, (char *)
01170           "atomselect: bad frame number in input, must be "
01171           "'first', 'last', 'now', or a non-negative number", TCL_STATIC);
01172         return TCL_ERROR;
01173       }
01174     }
01175     int val;
01176     // get the molecule name recursively
01177     if (!strcmp(argv[1], "top")) {
01178       if (Tcl_VarEval(interp, argv[0], " top", NULL) != TCL_OK ||
01179           Tcl_GetInt(interp, Tcl_GetStringResult(interp), &val) != TCL_OK     ) {
01180         return TCL_ERROR;
01181       }
01182     } else {
01183       if (Tcl_GetInt(interp, argv[1], &val) != TCL_OK) {
01184         return TCL_ERROR;
01185       }
01186     }
01187     Tcl_ResetResult(interp);
01188 
01189     // split the data into the various terms
01190     const char **list;
01191     int num_list;
01192     if (Tcl_SplitList(interp, argv[3], &num_list, &list) != TCL_OK) {
01193       return TCL_ERROR;
01194     }
01195     // and return the information
01196     int result = molinfo_get(context, val, num_list, list, interp, frame_num);
01197     ckfree((char *) list); // free of tcl data
01198 
01199     return result;
01200   }
01201 // Option  :  molinfo {molecule number} set <data> <new data> [frame <number>]
01202 //  Does   :   sets the given data for that molecule
01203   if ((argc == 5 && !strcmp(argv[2], "set")) ||
01204       (argc == 7 && !strcmp(argv[2], "set") && !strcmp(argv[5], "frame"))) {
01205     // get the frame number
01206     int frame_num;
01207     if (argc == 5) {
01208       frame_num = AtomSel::TS_NOW;
01209     } else {
01210       if (AtomSel::get_frame_value(argv[6], &frame_num) != 0) {
01211         Tcl_SetResult(interp, (char *)
01212           "atomselect: bad frame number in input, must be "
01213           "'first', 'last', 'now', or a non-negative number", TCL_STATIC);
01214         return TCL_ERROR;
01215       }
01216     }
01217 
01218     int val;
01219     if (!strcmp(argv[1], "top")) {
01220       if (Tcl_VarEval(interp, argv[0], " top", NULL) != TCL_OK ||
01221           Tcl_GetInt(interp, Tcl_GetStringResult(interp), &val) != TCL_OK     ) {
01222         return TCL_ERROR;
01223       }
01224     } else {
01225       if (Tcl_GetInt(interp, argv[1], &val) != TCL_OK) {
01226         return TCL_ERROR;
01227       }
01228     }
01229 
01230     Tcl_ResetResult(interp);
01231 
01232     // make sure the two lists have the same number of terms
01233     const char **list1, **list2;
01234     int num_list1, num_list2;
01235     if (Tcl_SplitList(interp, argv[3], &num_list1, &list1) != TCL_OK) {
01236       return TCL_ERROR;
01237     }
01238     if (Tcl_SplitList(interp, argv[4], &num_list2, &list2) != TCL_OK) {
01239       ckfree((char *)list1); // free of tcl data
01240       return TCL_ERROR;
01241     }
01242     if (num_list1 != num_list2) {
01243       ckfree((char *)list1); // free of tcl data
01244       ckfree((char *)list2); // free of tcl data
01245       Tcl_SetResult(interp, (char *) "molinfo: set: argument and value lists have different sizes", TCL_STATIC);
01246       return TCL_ERROR;
01247     }
01248 
01249     // call the 'set' routine
01250     int result = molinfo_set(context, val, num_list1, list1, list2, interp, frame_num);
01251 
01252     ckfree((char *)list1); // free of tcl data
01253     ckfree((char *)list2); // free of tcl data
01254     return result;
01255   }
01256 
01257   // There's been an error; find out what's wrong
01258   Tcl_SetResult(interp, (char *) "molinfo: called with unknown command", TCL_STATIC);
01259   if (argc >= 3) {
01260     if (!strcmp(argv[2], "get")) {
01261       Tcl_SetResult(interp, (char *) "molinfo: incorrect format for 'get'", 
01262                     TCL_STATIC);
01263     } else if (!strcmp(argv[2], "set")) {
01264       Tcl_SetResult(interp, (char *) "molinfo: incorrect format for 'set'",
01265                     TCL_STATIC);
01266     }
01267   } else if (argc >= 2) {
01268     if (!strcmp(argv[1], "get") || strcmp(argv[1], "set")) {
01269       Tcl_SetResult(interp, (char *) "molinfo: missing molecule number",
01270                     TCL_STATIC);
01271     }
01272   }
01273 
01274   return TCL_ERROR;
01275 }