TclMolInfo.C Source File

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: TclMolInfo.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.98 $       $Date: 2019/10/31 15:31:03 $
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, QMDATA_BUFSTRFMT, 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, 3L*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, 2L*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, 3L*mol->nAtoms);
00719     generic_molinfo_qm_mat("%.6f", "inthessian",  inthessian,  mol->qm_data->get_num_intcoords());
00720     generic_molinfo_qm_mat("%.8f", "normalmodes", normalmodes, 3L*mol->nAtoms);
00721     generic_molinfo_qm_arr("%.2f", "wavenumbers", wavenumbers, 3L*mol->nAtoms);
00722     generic_molinfo_qm_arr("%.6f", "intensities", intensities, 3L*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[3L*i]));
00745         Tcl_ListObjAppendElement(interp, adata,
00746                                  Tcl_NewIntObj(mol->angles[3L*i+1]));
00747         Tcl_ListObjAppendElement(interp, adata,
00748                                  Tcl_NewIntObj(mol->angles[3L*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[4L*i]));
00772         Tcl_ListObjAppendElement(interp, adata,
00773                                  Tcl_NewIntObj(mol->dihedrals[4L*i+1]));
00774         Tcl_ListObjAppendElement(interp, adata,
00775                                  Tcl_NewIntObj(mol->dihedrals[4L*i+2]));
00776         Tcl_ListObjAppendElement(interp, adata,
00777                                  Tcl_NewIntObj(mol->dihedrals[4L*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[4L*i]));
00801         Tcl_ListObjAppendElement(interp, adata,
00802                                  Tcl_NewIntObj(mol->impropers[4L*i+1]));
00803         Tcl_ListObjAppendElement(interp, adata,
00804                                  Tcl_NewIntObj(mol->impropers[4L*i+2]));
00805         Tcl_ListObjAppendElement(interp, adata,
00806                                  Tcl_NewIntObj(mol->impropers[4L*i+3]));
00807         Tcl_ListObjAppendElement(interp, alist, adata);
00808       }
00809       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), alist);
00810     } else if (!strcmp(arg, "crossterms")) {
00811       Tcl_Obj *alist = Tcl_NewListObj(0, NULL);
00812       for (int i=0; i<mol->num_cterms(); i++) {
00813         Tcl_Obj *adata = Tcl_NewListObj(0,NULL);
00814         for (int j=0; j<8; j++) {
00815           Tcl_ListObjAppendElement(interp, adata, Tcl_NewIntObj(mol->cterms[8L*i+j]));
00816         }
00817         Tcl_ListObjAppendElement(interp, alist, adata);
00818       }
00819       Tcl_ListObjAppendElement(interp, Tcl_GetObjResult(interp), alist);
00820     } else {
00821       Tcl_ResetResult(interp);
00822       Tcl_AppendResult(interp, "molinfo: cannot find molinfo attribute '",
00823                        argv[term], "'", NULL);
00824       return TCL_ERROR;
00825     }
00826   }
00827 
00828   return TCL_OK;
00829 }
00830 
00831 
00832 #define generic_molinfo_set_data(name, func1, func2) \
00833 } else if (!strcmp(argv[term], name)) { \
00834 int onoff; \
00835 if (Tcl_GetBoolean(interp, data[term], &onoff) != TCL_OK) return TCL_ERROR; \
00836 if (onoff) { func1 ; } else { func2 ; }
00837 
00838 #define generic_molinfo_simulation_set(name, type)  \
00839 } else if (!strcmp(argv[term], name)) { \
00840 Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00841 if (ts) ts->energy[type] = (float) atof(data[term]);
00842 
00843 #define generic_molinfo_pbc_set(name, type) \
00844 } else if (!strcmp(argv[term], name)) { \
00845 Timestep *ts = mol->get_frame(context.molinfo_frame_number); \
00846 if (ts) { ts->type = (float) atof(data[term]); mol->change_pbc(); }
00847 
00848 // given the list of strings and data, set the right values
00849 static int molinfo_set(ctxt context, int molid, int argc, const char *argv[],
00850                        const char *data[], Tcl_Interp *interp, int frame_num) {
00851   // does the molecule exist?
00852   context.molinfo_get_id = molid;
00853   context.molinfo_get_index = context.mlist->mol_index_from_id(molid);
00854 
00855   if (context.molinfo_get_index == -1) {
00856     char s[10];
00857     sprintf(s, "%d", molid);
00858     Tcl_AppendResult(interp, "molinfo: set: no molecule exists with id ",
00859                      s, NULL);
00860     return TCL_ERROR;
00861   }
00862 
00863   // get the right frame number
00864   switch (frame_num) {
00865     case AtomSel::TS_NOW: 
00866       context.molinfo_frame_number =
00867         context.mlist->molecule(context.molinfo_get_index) -> frame();
00868       break;
00869     case AtomSel::TS_LAST:
00870       context.molinfo_frame_number = 
00871         context.mlist->molecule(context.molinfo_get_index) -> numframes()-1;
00872       break;
00873     default:      
00874       context.molinfo_frame_number = frame_num;
00875   }
00876   Molecule *mol = context.mlist->molecule(context.molinfo_get_index);
00877   
00878   for (int term=0; term<argc; term++) {
00879     context.molinfo_get_string = argv[term];
00880     if (!strcmp(argv[term], "center")) {
00881       float x, y, z;
00882       if (sscanf((const char *)data[term], "%f %f %f", &x, &y, &z) != 3) {
00883         Tcl_AppendResult(interp, 
00884             "molinfo: set center: must have three position elements", NULL);
00885         return TCL_ERROR;
00886       }
00887       mol->change_center(x, y, z);
00888     } else if (!strcmp(argv[term], "center_matrix")) {
00889       float mat[16];
00890       if (!read_matrix(interp, data[term], mat)) return TCL_ERROR;
00891       mol->set_cent_trans(mat[12], mat[13], mat[14]);
00892     } else if (!strcmp(argv[term], "rotate_matrix")) {
00893       Matrix4 mat;
00894       if (!read_matrix(interp, data[term], mat.mat)) return TCL_ERROR;
00895       mol->set_rot(mat);
00896     } else if (!strcmp(argv[term], "scale_matrix")) {
00897       float mat[16];
00898       if (!read_matrix(interp, data[term], mat)) return TCL_ERROR;
00899       mol->set_scale(mat[0]);
00900     } else if (!strcmp(argv[term], "global_matrix")) {
00901       float mat[16];
00902       if (!read_matrix(interp, data[term], mat)) return TCL_ERROR;
00903       mol->set_glob_trans(mat[12], mat[13], mat[14]);
00904     } else if (!strcmp(argv[term], "angles")) {
00905       mol->clear_angles();
00906       mol->angleTypeNames.clear();
00907       mol->set_dataset_flag(BaseMolecule::ANGLES);
00908 
00909       Tcl_Obj *alist = Tcl_NewListObj(0,NULL);
00910       Tcl_SetStringObj(alist, data[term], -1);
00911       Tcl_IncrRefCount(alist); 
00912 
00913       int numangles;
00914       Tcl_Obj **adata;
00915       Tcl_ListObjGetElements(interp, alist, &numangles, &adata);
00916 
00917       for (int i=0; i < numangles; i++) {
00918         int numentries;
00919         Tcl_Obj **a;
00920 
00921         Tcl_ListObjGetElements(interp, adata[i], &numentries, &a);
00922         if (numentries != 4) {
00923           Tcl_AppendResult(interp, "molinfo: incorrect data item for "
00924                            "'set angles' :", Tcl_GetString(adata[i]), NULL);
00925           return TCL_ERROR;
00926         }
00927         int type = mol->angleTypeNames.add_name(Tcl_GetString(a[0]), 0);
00928         int a1, a2, a3;
00929         if (Tcl_GetIntFromObj(interp, a[1], &a1) != TCL_OK) return TCL_ERROR;
00930         if (Tcl_GetIntFromObj(interp, a[2], &a2) != TCL_OK) return TCL_ERROR;
00931         if (Tcl_GetIntFromObj(interp, a[3], &a3) != TCL_OK) return TCL_ERROR;
00932         mol->add_angle(a1, a2, a3, type);
00933       }
00934       // release storage
00935       Tcl_InvalidateStringRep(alist);
00936       Tcl_DecrRefCount(alist); 
00937 
00938     } else if (!strcmp(argv[term], "dihedrals")) { 
00939       mol->clear_dihedrals();
00940       mol->dihedralTypeNames.clear();
00941       mol->set_dataset_flag(BaseMolecule::ANGLES);
00942 
00943       Tcl_Obj *alist = Tcl_NewListObj(0,NULL);
00944       Tcl_SetStringObj(alist, data[term], -1);
00945       Tcl_IncrRefCount(alist); 
00946 
00947       int numdihedrals;
00948       Tcl_Obj **adata;
00949       Tcl_ListObjGetElements(interp, alist, &numdihedrals, &adata);
00950       for (int i=0; i < numdihedrals; i++) {
00951         int numentries;
00952         Tcl_Obj **a;
00953         Tcl_ListObjGetElements(interp, adata[i], &numentries, &a);
00954         if (numentries != 5) {
00955           Tcl_AppendResult(interp, "molinfo: incorrect data item for "
00956                            "'set dihedrals' :", Tcl_GetString(adata[i]), NULL);
00957           return TCL_ERROR;
00958         }
00959         int type = mol->dihedralTypeNames.add_name(Tcl_GetString(a[0]), 0);
00960         int a1, a2, a3, a4;
00961         if (Tcl_GetIntFromObj(interp, a[1], &a1) != TCL_OK) return TCL_ERROR;
00962         if (Tcl_GetIntFromObj(interp, a[2], &a2) != TCL_OK) return TCL_ERROR;
00963         if (Tcl_GetIntFromObj(interp, a[3], &a3) != TCL_OK) return TCL_ERROR;
00964         if (Tcl_GetIntFromObj(interp, a[4], &a4) != TCL_OK) return TCL_ERROR;
00965         mol->add_dihedral(a1, a2, a3, a4, type);
00966       }
00967       // release storage
00968       Tcl_InvalidateStringRep(alist);
00969       Tcl_DecrRefCount(alist); 
00970 
00971     } else if (!strcmp(argv[term], "impropers")) {
00972       mol->clear_impropers();
00973       mol->improperTypeNames.clear();
00974       mol->set_dataset_flag(BaseMolecule::ANGLES);
00975 
00976       Tcl_Obj *alist = Tcl_NewListObj(0,NULL);
00977       Tcl_SetStringObj(alist, data[term], -1);
00978       Tcl_IncrRefCount(alist); 
00979 
00980       int numimpropers;
00981       Tcl_Obj **adata;
00982       Tcl_ListObjGetElements(interp, alist, &numimpropers, &adata);
00983       for (int i=0; i < numimpropers; i++) {
00984         int numentries;
00985         Tcl_Obj **a;
00986         Tcl_ListObjGetElements(interp, adata[i], &numentries, &a);
00987         if (numentries != 5) {
00988           Tcl_AppendResult(interp, "molinfo: incorrect data item for "
00989                            "'set impropers' :", Tcl_GetString(adata[i]), NULL);
00990           return TCL_ERROR;
00991         }
00992         int type = mol->improperTypeNames.add_name(Tcl_GetString(a[0]), 0);
00993         int a1, a2, a3, a4;
00994         if (Tcl_GetIntFromObj(interp, a[1], &a1) != TCL_OK) return TCL_ERROR;
00995         if (Tcl_GetIntFromObj(interp, a[2], &a2) != TCL_OK) return TCL_ERROR;
00996         if (Tcl_GetIntFromObj(interp, a[3], &a3) != TCL_OK) return TCL_ERROR;
00997         if (Tcl_GetIntFromObj(interp, a[4], &a4) != TCL_OK) return TCL_ERROR;
00998         mol->add_improper(a1, a2, a3, a4, type);
00999       }
01000       // release storage
01001       Tcl_InvalidateStringRep(alist);
01002       Tcl_DecrRefCount(alist); 
01003     } else if (!strcmp(argv[term], "crossterms")) {
01004       mol->clear_cterms();
01005       mol->set_dataset_flag(BaseMolecule::CTERMS);
01006 
01007       Tcl_Obj *alist = Tcl_NewListObj(0,NULL);
01008       Tcl_SetStringObj(alist, data[term], -1);
01009       Tcl_IncrRefCount(alist); 
01010 
01011       int numcterms;
01012       Tcl_Obj **adata;
01013       Tcl_ListObjGetElements(interp, alist, &numcterms, &adata);
01014       for (int i=0; i < numcterms; i++) {
01015         int numentries;
01016         Tcl_Obj **a;
01017         Tcl_ListObjGetElements(interp, adata[i], &numentries, &a);
01018         if (numentries != 8) {
01019           Tcl_AppendResult(interp, "molinfo: incorrect data item for "
01020                            "'set crossterms' :", Tcl_GetString(adata[i]), NULL);
01021           return TCL_ERROR;
01022         }
01023         int a1, a2, a3, a4, a5, a6, a7, a8;
01024         if (Tcl_GetIntFromObj(interp, a[0], &a1) != TCL_OK) return TCL_ERROR;
01025         if (Tcl_GetIntFromObj(interp, a[1], &a2) != TCL_OK) return TCL_ERROR;
01026         if (Tcl_GetIntFromObj(interp, a[2], &a3) != TCL_OK) return TCL_ERROR;
01027         if (Tcl_GetIntFromObj(interp, a[3], &a4) != TCL_OK) return TCL_ERROR;
01028         if (Tcl_GetIntFromObj(interp, a[4], &a5) != TCL_OK) return TCL_ERROR;
01029         if (Tcl_GetIntFromObj(interp, a[5], &a6) != TCL_OK) return TCL_ERROR;
01030         if (Tcl_GetIntFromObj(interp, a[6], &a7) != TCL_OK) return TCL_ERROR;
01031         if (Tcl_GetIntFromObj(interp, a[7], &a8) != TCL_OK) return TCL_ERROR;
01032         mol->add_cterm(a1, a2, a3, a4, a5, a6, a7, a8);
01033       }
01034       // release storage
01035       Tcl_InvalidateStringRep(alist);
01036       Tcl_DecrRefCount(alist); 
01037     } else if (!strcmp(argv[term], "frame")) {
01038       // XXX this isn't correctly sending events to the GUI,
01039       // so if you do 'molinfo top set frame 5' the GUI is 
01040       // out of date
01041       //
01042       // This is correct and intended behavior; if you want GUI events,
01043       // the command to use is "animate goto <frame>. -- JRG 12/27/07
01044       mol->override_current_frame(atoi(data[term]));
01045       mol->change_ts();
01046       
01047     generic_molinfo_set_data("active", context.app->molecule_activate(mol->id(), 1), context.app->molecule_activate(mol->id(), 0))
01048     generic_molinfo_set_data("drawn", context.app->molecule_display(mol->id(), 1), context.app->molecule_display(mol->id(), 0))
01049     generic_molinfo_set_data("displayed", context.app->molecule_display(mol->id(), 1), context.app->molecule_display(mol->id(), 0))
01050     generic_molinfo_set_data("fixed", context.app->molecule_fix(mol->id(), 1), context.app->molecule_fix(mol->id(), 0))
01051     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)
01052 
01053     generic_molinfo_simulation_set("bond", TSE_BOND)
01054     generic_molinfo_simulation_set("angle", TSE_ANGLE)
01055     generic_molinfo_simulation_set("dihedral", TSE_DIHE)
01056     generic_molinfo_simulation_set("improper", TSE_IMPR)
01057     generic_molinfo_simulation_set("vdw", TSE_VDW)
01058     generic_molinfo_simulation_set("electrostatic", TSE_COUL)
01059     generic_molinfo_simulation_set("elec", TSE_COUL)
01060     generic_molinfo_simulation_set("hbond", TSE_HBOND)
01061     generic_molinfo_simulation_set("kinetic", TSE_KE)
01062     generic_molinfo_simulation_set("potential", TSE_PE)
01063     generic_molinfo_simulation_set("temperature", TSE_TEMP)
01064     generic_molinfo_simulation_set("temp", TSE_TEMP)
01065     generic_molinfo_simulation_set("energy", TSE_TOTAL)
01066     generic_molinfo_simulation_set("volume", TSE_VOLUME)
01067     generic_molinfo_simulation_set("pressure", TSE_PRESSURE)
01068     generic_molinfo_simulation_set("efield", TSE_EFIELD)
01069     generic_molinfo_simulation_set("urey_bradley", TSE_UREY_BRADLEY)
01070     generic_molinfo_simulation_set("molinfo_restraint", TSE_RESTRAINT)
01071 
01072     } else if (!strcmp(argv[term], "timesteps")) {
01073       Timestep *ts = mol->get_frame(context.molinfo_frame_number);
01074       if (ts) ts->timesteps = atoi(data[term]);
01075 
01076     generic_molinfo_pbc_set("a", a_length)
01077     generic_molinfo_pbc_set("b", b_length)
01078     generic_molinfo_pbc_set("c", c_length)
01079     generic_molinfo_pbc_set("alpha", alpha)
01080     generic_molinfo_pbc_set("beta", beta)
01081     generic_molinfo_pbc_set("gamma", gamma)
01082     generic_molinfo_pbc_set("physical_time", physical_time)
01083 
01084     } else {
01085       Tcl_AppendResult(interp, "molinfo: cannot find molinfo attribute '",
01086                        argv[term], "'", NULL);
01087       return TCL_ERROR;
01088     }
01089   }
01090 
01091   return TCL_OK;
01092 }
01093 
01094 
01095 // Function:  molinfo
01096 // Option  :  molinfo num
01097 //  Returns:   number of molecules
01098 // Option  :  molinfo index <int>
01099 //  Returns:   molecule id of the nth molecule (starting at index = 0)
01100 // Option  :  molinfo list
01101 //  Returns:   list of all molecule ids
01102 // Option  :  molinfo top
01103 //  Returns:   molecule id of the 'top' molecule
01104 // Option  :  molinfo {molecule number} get <data>
01105 //  Returns:   the given data for that molecule
01106 // Option  :  molinfo {molecule number} set <data> <data fields>
01107 //  Does (okay, this isn't a 'info' thing): sets the data field(s)
01108 
01109 int molecule_tcl(ClientData data, Tcl_Interp *interp, int argc, const char *argv[])
01110 {
01111   VMDApp *app = (VMDApp *)data;
01112   // set context variable here
01113   ctxt context;
01114   context.mlist = app->moleculeList;
01115   context.matlist = app->materialList;
01116   context.app = app;
01117 
01118   if (argc == 1) {
01119     Tcl_SetResult(interp,
01120       (char *)
01121       "usage: molinfo <command> [args...]\n\n"
01122       "Commands:"
01123       "\nMolecule IDs:\n"
01124       "  list                  -- lists all existing molecule IDs\n"
01125       "  num                   -- number of loaded molecules\n"
01126       "  top                   -- gets ID of top molecule (or -1 if none)\n"
01127       "  index <n>             -- gets ID of n-th molecule\n"
01128       "\nGetting and Setting Molecular Information:\n" 
01129      // "  keywords              -- returns a list of molinfo keywords\n"     // XXX obsolete???
01130       "  <molid> get <(list of) keywords>\n"
01131       "  <molid> set <(list of) keywords> <(list of) values>\n",
01132       TCL_STATIC);
01133     return TCL_ERROR;
01134   }
01135 
01136   // what does it want?
01137   if (argc == 2) {
01138 
01139 // Option  :  molinfo num
01140 //  Returns:   number of molecules
01141     if (!strcmp(argv[1], "num")) {
01142       char tmpstring[64];
01143       sprintf(tmpstring, "%d", context.mlist->num());
01144       Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
01145       return TCL_OK;
01146     }
01147 
01148 // Option  :  molinfo list
01149 //  Returns:   list of all molecule ids
01150     if (!strcmp(argv[1], "list")) {
01151       if (context.mlist->num() <= 0) {
01152         return TCL_OK;
01153       }
01154       char s[20];
01155       sprintf(s, "%d", context.mlist->molecule(0)->id());
01156       Tcl_AppendResult(interp, s, (char *) NULL);
01157       for (int i=1; i<context.mlist -> num(); i++) {
01158         sprintf(s, "%d", context.mlist->molecule(i)->id());
01159         Tcl_AppendResult(interp, " ", s, (char *) NULL);
01160       }
01161       return TCL_OK;
01162     }
01163 // Option  :  molinfo top
01164 //  Returns:   molecule id of the 'top' molecule
01165     if (!strcmp(argv[1], "top")) {
01166       if (context.mlist->top()) {
01167         char tmpstring[64];
01168         sprintf(tmpstring, "%d", context.mlist->top()->id());
01169         Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
01170       } else {
01171         Tcl_SetResult(interp, (char *) "-1", TCL_STATIC);
01172       }
01173       return TCL_OK;
01174     }
01175 
01176     // otherwise, I don't know
01177     Tcl_AppendResult(interp, "molinfo: couldn't understand '",
01178                      argv[1], "'", NULL);
01179     return TCL_ERROR;
01180   } // end of commands with only one option
01181 
01182   if (argc == 3) { // commands with two options
01183     int val;
01184     if (Tcl_GetInt(interp, argv[2], &val) != TCL_OK) {
01185       return TCL_ERROR;
01186     }
01187 // Option  :  molecule index <int>
01188 //  Returns:   molecule id of the nth molecule (starting at index = 0)
01189     if (!strcmp(argv[1], "index")) {
01190       if (context.mlist->molecule(val)) {
01191         char tmpstring[64];
01192         sprintf(tmpstring, "%d", context.mlist->molecule(val)->id());
01193         Tcl_SetResult(interp, tmpstring, TCL_VOLATILE);
01194       } else {
01195         Tcl_SetResult(interp, (char *) "-1", TCL_STATIC);
01196       }
01197       return TCL_OK;
01198     }
01199     Tcl_AppendResult(interp, "molinfo: couldn't understand '",
01200                      argv[1], "'", NULL);
01201     return TCL_ERROR;
01202   }
01203 // Option  :  molinfo {molecule number} get <data> [frame <number>]
01204 //  Returns:   the given data for that molecule
01205   if ((argc == 4 && !strcmp(argv[2], "get")) ||
01206       (argc == 6 && !strcmp(argv[2], "get") && !strcmp(argv[4], "frame"))) {
01207     int frame_num;
01208     if (argc == 4) {
01209       frame_num = AtomSel::TS_NOW;
01210     } else {
01211       if (AtomSel::get_frame_value(argv[5], &frame_num) != 0) {
01212         Tcl_SetResult(interp, (char *)
01213           "atomselect: bad frame number in input, must be "
01214           "'first', 'last', 'now', or a non-negative number", TCL_STATIC);
01215         return TCL_ERROR;
01216       }
01217     }
01218     int val;
01219     // get the molecule name recursively
01220     if (!strcmp(argv[1], "top")) {
01221       if (Tcl_VarEval(interp, argv[0], " top", NULL) != TCL_OK ||
01222           Tcl_GetInt(interp, Tcl_GetStringResult(interp), &val) != TCL_OK     ) {
01223         return TCL_ERROR;
01224       }
01225     } else {
01226       if (Tcl_GetInt(interp, argv[1], &val) != TCL_OK) {
01227         return TCL_ERROR;
01228       }
01229     }
01230     Tcl_ResetResult(interp);
01231 
01232     // split the data into the various terms
01233     const char **list;
01234     int num_list;
01235     if (Tcl_SplitList(interp, argv[3], &num_list, &list) != TCL_OK) {
01236       return TCL_ERROR;
01237     }
01238     // and return the information
01239     int result = molinfo_get(context, val, num_list, list, interp, frame_num);
01240     ckfree((char *) list); // free of tcl data
01241 
01242     return result;
01243   }
01244 // Option  :  molinfo {molecule number} set <data> <new data> [frame <number>]
01245 //  Does   :   sets the given data for that molecule
01246   if ((argc == 5 && !strcmp(argv[2], "set")) ||
01247       (argc == 7 && !strcmp(argv[2], "set") && !strcmp(argv[5], "frame"))) {
01248     // get the frame number
01249     int frame_num;
01250     if (argc == 5) {
01251       frame_num = AtomSel::TS_NOW;
01252     } else {
01253       if (AtomSel::get_frame_value(argv[6], &frame_num) != 0) {
01254         Tcl_SetResult(interp, (char *)
01255           "atomselect: bad frame number in input, must be "
01256           "'first', 'last', 'now', or a non-negative number", TCL_STATIC);
01257         return TCL_ERROR;
01258       }
01259     }
01260 
01261     int val;
01262     if (!strcmp(argv[1], "top")) {
01263       if (Tcl_VarEval(interp, argv[0], " top", NULL) != TCL_OK ||
01264           Tcl_GetInt(interp, Tcl_GetStringResult(interp), &val) != TCL_OK     ) {
01265         return TCL_ERROR;
01266       }
01267     } else {
01268       if (Tcl_GetInt(interp, argv[1], &val) != TCL_OK) {
01269         return TCL_ERROR;
01270       }
01271     }
01272 
01273     Tcl_ResetResult(interp);
01274 
01275     // make sure the two lists have the same number of terms
01276     const char **list1, **list2;
01277     int num_list1, num_list2;
01278     if (Tcl_SplitList(interp, argv[3], &num_list1, &list1) != TCL_OK) {
01279       return TCL_ERROR;
01280     }
01281     if (Tcl_SplitList(interp, argv[4], &num_list2, &list2) != TCL_OK) {
01282       ckfree((char *)list1); // free of tcl data
01283       return TCL_ERROR;
01284     }
01285     if (num_list1 != num_list2) {
01286       ckfree((char *)list1); // free of tcl data
01287       ckfree((char *)list2); // free of tcl data
01288       Tcl_SetResult(interp, (char *) "molinfo: set: argument and value lists have different sizes", TCL_STATIC);
01289       return TCL_ERROR;
01290     }
01291 
01292     // call the 'set' routine
01293     int result = molinfo_set(context, val, num_list1, list1, list2, interp, frame_num);
01294 
01295     ckfree((char *)list1); // free of tcl data
01296     ckfree((char *)list2); // free of tcl data
01297     return result;
01298   }
01299 
01300   // There's been an error; find out what's wrong
01301   Tcl_SetResult(interp, (char *) "molinfo: called with unknown command", TCL_STATIC);
01302   if (argc >= 3) {
01303     if (!strcmp(argv[2], "get")) {
01304       Tcl_SetResult(interp, (char *) "molinfo: incorrect format for 'get'", 
01305                     TCL_STATIC);
01306     } else if (!strcmp(argv[2], "set")) {
01307       Tcl_SetResult(interp, (char *) "molinfo: incorrect format for 'set'",
01308                     TCL_STATIC);
01309     }
01310   } else if (argc >= 2) {
01311     if (!strcmp(argv[1], "get") || strcmp(argv[1], "set")) {
01312       Tcl_SetResult(interp, (char *) "molinfo: missing molecule number",
01313                     TCL_STATIC);
01314     }
01315   }
01316 
01317   return TCL_ERROR;
01318 }