Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

gamessplugin.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2016 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: gamessplugin.c,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.211 $       $Date: 2016/11/28 05:01:54 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* *******************************************************
00019  *
00020  *          G A M E S S     P L U G I N 
00021  *
00022  * Read GAMESS-US log files.
00023  * This file reader is fairly robust and should be able to 
00024  * parse GAMESS logfiles from a wide range of versions
00025  * (2000-2009 tested). 
00026  * From the variety of data that can be found in GAMESS
00027  * logfiles we read coordinates, basis set, wavefunctions, 
00028  * gradients, hessian, charges, frequencies and a number 
00029  * of data describing the type of calculation.
00030  *
00031  * ********************************************************/
00032 
00033 
00034 /**********************************************************
00035 
00036  
00037 FUNCTION CALL CHAIN
00038 ===================
00039 
00040 Below is an overview about the hierarchy of function calls.
00041 Not all functions or possible branches are listed but the
00042 general picture is covered.
00043 
00044 The top level consists of functions defined by the 
00045 molfile_plugin interface: First VMD calls open_gamess_read()
00046 then it requests information about atoms and topology using 
00047 read_gamess_structure(). Note that no coordinates are provided
00048 at this step. Timestep independent data such as info about the
00049 calculation method or the basis set are provided through 
00050 read_gamess_rundata(). Since the allocation of these
00051 arrays populated for this purpose is done by VMD rather than
00052 the plugin, VMD needs to know the sizes beforehand. This is 
00053 achieved by calling read_qm_metadata() before read_qm_rundata().
00054 
00055 Next, in order to obtain the info for all timesteps VMD will call
00056 read_timestep_metadata(), read_qm_timestep_metadata(), and
00057 read_timestep() repeatedly until read_timestep() returns
00058 MOLFILE_ERROR indicating the end of the trajectory.
00059 Here too, the metadata are transferred before the main chunk
00060 of data.
00061 
00062 Finally, VMD calls close_gamess_read() which frees the temporary
00063 memory used by the plugin and closes the file.
00064 
00065 
00066 VMD
00067  |__ open_gamess_read()
00068  |    |__ parse_static_data()
00069  |        |__ get_proc_mem()
00070  |        |__ get_basis_options()
00071  |        |__ get_runtitle()
00072  |        |__ get_contrl()
00073  |        |__ get_input_structure()
00074  |        |    |__ get_coordinates()
00075  |        |
00076  |        |__ get_basis()
00077  |        |    |__ read_shell_primitives()
00078  |        |    |__ fill_basis_arrays()
00079  |        |
00080  |        |__ get_basis_stats()
00081  |        |__ get_properties_input()
00082  |        |__ get_int_coords()
00083  |        |__ get_symmetry()
00084  |        |__ get_guess_options()
00085  |        |__ get_mcscf()
00086  |        |__ analyze_traj()
00087  |        |__ read_first_frame()
00088  |        |    |__get_traj_frame()
00089  |        |        |__get_coordinates()
00090  |        |        |__get_scfdata()
00091  |        |        |__check_add_wavefunctions()
00092  |        |        |   |__add_wavefunction()
00093  |        |        |   |__get_wavefunction()
00094  |        |        |   |   |__read_coeff_block()
00095  |        |        |   |       |__angular_momentum_expon()
00096  |        |        |   |
00097  |        |        |   |__del_wavefunction()
00098  |        |        |
00099  |        |        |__get_population()
00100  |        |        |__get_gradient()
00101  |        |
00102  |        |__ get_final_properties()
00103  |             |__ get_population()
00104  |             |__ get_esp_charges()
00105  |             |__ get_final_gradient()
00106  |             |__ get_int_hessian()
00107  |             |__ get_cart_hessian()
00108  |             |__ get_normal_modes()
00109  |             |__ read_localized_orbitals()
00110  |        
00111  |
00112  |__ read_gamess_structure()
00113  |
00114  |__ read_gamess_metadata()
00115  |
00116  |__ read_gamess_rundata()
00117  |
00118  |
00119  |  DO:                         <----.
00120  |__ read_timestep_metadata()        |
00121  |__ read_qm_timestep_metadata()     |
00122  |__ read_timestep()                 |
00123  |                                   |
00124  |  WHILE(FRAMES)               -----' 
00125  |
00126  |
00127  |__ close_gamess_read()
00128 
00129  
00130 
00131 PARSING STRATEGY
00132 ================
00133 
00134 Because we potentially have to read quite a bit into the
00135 logfile in order to obtain the number of atoms required by 
00136 open_gamess_read(), we just parse the whole file and store 
00137 all timestep independent data. This process is managed by 
00138 parse_static_data(). Some of the static data are at the end
00139 of the file. Function analyze_traj() find the end of the
00140 trajectory and records the file pointer for the beginning
00141 of each frame. Thus, reading the frames when they are requested 
00142 later by read_timestep() is much faster without having to  
00143 keep large amountss of data in memory.
00144 
00145 
00146 **********************************************************/
00147 
00148 #include <stdlib.h>
00149 #include <stdarg.h>
00150 #include <string.h>
00151 #include <ctype.h>
00152 #include <errno.h>
00153 #include <time.h>
00154 
00155 #include "qmplugin.h"
00156 #include "unit_conversion.h"
00157 #include "periodic_table.h"
00158 
00159 #define ANGSTROM 0
00160 #define BOHR     1
00161 
00162 
00163 /* #define DEBUGGING 1 */
00164 
00165 /*
00166  * Error reporting macro for use in DEBUG mode
00167  */
00168 
00169 #ifdef DEBUGGING
00170 #define PRINTERR fprintf(stderr, "\n In file %s, line %d: \n %s \n \n", \
00171                             __FILE__, __LINE__, strerror(errno))
00172 #else
00173 #define PRINTERR (void)(0)
00174 #endif
00175 
00176 /*
00177  * Error reporting macro for the multiple fgets calls in
00178  * the code
00179  */
00180 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE
00181 
00182 
00183 #define NOTFOUND 0
00184 #define FOUND    1
00185 #define STOPPED  2
00186 
00187 
00188 #define GAMESSUNKNOWN      0
00189 #define GAMESSPRE20050627  1
00190 #define GAMESSPOST20050627 2
00191 #define FIREFLY8PRE6695    3
00192 #define FIREFLY8POST6695   4
00193 
00194 /* Data specific to parsing GAMESS files */
00195 typedef struct 
00196 {
00197   int  version; /* here we track the GAMESS versions, since the
00198                  * file format has changed with 
00199                  * version 27 JUN 2005 (R2);
00200                  * version = 1  : pre-27 JUN 2005 (R2)
00201                  * version = 2  : 27 JUN 2005 (R2) or later
00202                  * version = 0  : version unrecognized;
00203                  *                file will not be read */
00204 
00205   int have_pcgamess; /* this flag is set to 1 if the output
00206                       * file is recognized as a PC Gamess output
00207                       * file; we might need to introduce a few
00208                       * switches in the code depending on if
00209                       * the log file is plain Gamess or PC Gamess
00210                       */
00211 
00212   int have_fmo;      /* set to 1 if this is an FMO calculation */
00213 
00214 } gmsdata;
00215 
00216 
00217 /* ######################################################## */
00218 /* declaration/documentation of internal (static) functions */
00219 /* ######################################################## */
00220 
00221 /* this routine is the main gamess log file
00222  * parser responsible for static, i.e. 
00223  * non-trajectory information */
00224 static int parse_static_data(qmdata_t *, int *);
00225 
00226 /* for debugging */
00227 static void print_input_data(qmdata_t *);
00228 
00229 /* this routine checks if the current run is an
00230  * actual GAMESS run; returns true/false */
00231 static int have_gamess(qmdata_t *, gmsdata *);
00232 
00233 /* this function reads the number of processors requested */
00234 static int get_proc_mem(qmdata_t *, gmsdata *);
00235 
00236 /* Parse the $BASIS options*/
00237 static int get_basis_options(qmdata_t *);
00238 
00239 /* Determine the run title line */
00240 static int get_runtitle(qmdata_t *);
00241 
00242 /* Read the input atom definitions and geometry */
00243 static int get_input_structure(qmdata_t *data, gmsdata *gms);
00244 
00245 /* Read basis set and orbital statistics such as
00246  * # of shells, # of A/B orbitals, # of electrons, 
00247  * multiplicity and total charge */
00248 static int get_basis_stats(qmdata_t *);
00249 
00250 /* Read the contrl group for firefly calc and check for
00251  * supported RUNTYPes. Terminate the plugin
00252  * if an unsupported one is encountered. */
00253 static int get_contrl_firefly(qmdata_t *);
00254 
00255 /* Read the contrl group and check for
00256  * supported RUNTYPes. Terminate the plugin
00257  * if an unsupported one is encountered. */
00258 static int get_contrl(qmdata_t *);
00259 
00260 /* Read input parameters regarding calculation of 
00261  * certain molecular properties such as electrostatic
00262  * moments and the MEP. */
00263 static int get_properties_input(qmdata_t *);
00264 
00265 /* Read symmetry point group and highest axis */
00266 static int get_symmetry(qmdata_t *);
00267 
00268 /* Read in the $GUESS options */
00269 static int get_guess_options(qmdata_t *);
00270 
00271 /* Read MCSCF data */
00272 static int get_mcscf(qmdata_t *data);
00273 
00274 /* the function get_initial_info provides the atom number,
00275  * coordinates, and atom types and stores them
00276  * temporarily. */ 
00277 static int get_final_properties (qmdata_t *);
00278 
00279 /* Read atom coordinate block */
00280 static int get_coordinates(FILE *file, qm_atom_t **atoms, int unit,
00281                            int *numatoms);
00282 
00283 /* Read coordinates from $FMOXYZ section in the INPUT CARD
00284  * listing at the beginning of the file. */
00285 static int get_fmoxyz(FILE *file, qm_atom_t **atoms, int unit,
00286                       int *numatoms);
00287 
00288 /* Read the basis set data */
00289 static int get_basis (qmdata_t *);
00290 
00291 /* Read all primitives for the current shell */
00292 static int read_shell_primitives(qmdata_t *, prim_t **prim,
00293                                  char *shelltype, int icoeff, int pcgamess);
00294 
00295 /* convert shell type from char to int */
00296 static int shelltype_int(char type);
00297 
00298 /* Populate the flat arrays containing the basis set data */
00299 static int fill_basis_arrays(qmdata_t *);
00300 
00301 /* Read the first trajectory frame. */
00302 static int read_first_frame(qmdata_t *);
00303 
00304 /* Read next trajectory frame. */
00305 static int get_traj_frame(qmdata_t *, qm_atom_t *, int);
00306 
00307 /* returns 1 if the optimization has converged */
00308 static int analyze_traj(qmdata_t *, gmsdata *);
00309 
00310 /* read the number of scf iterations and the scf energies
00311  * for the current timestep. */
00312 static int get_scfdata(qmdata_t *, qm_timestep_t *);
00313 
00314 /* Reads a set of wavefunctions for the current timestep.
00315  * (typically alpha and beta spin wavefunctions) */
00316 static int check_add_wavefunctions(qmdata_t *data,
00317                                    qm_timestep_t *ts);
00318 
00319 /* Parse the wavefunction. */
00320 static int get_wavefunction(qmdata_t *, qm_timestep_t *, qm_wavefunction_t *);
00321 
00322 /* Read the wavefunction coefficients. */
00323 static int read_coeff_block(FILE *file, int wavef_size,
00324                               float *wave_coeff, int *angular_momentum);
00325 
00326 /* Read localized orbitals (Boys/Ruedenberg/Pipek) */
00327 static int read_localized_orbitals(qmdata_t *data);
00328 
00329 /* Read population analysis (Mulliken and Lowdin charges) */
00330 static int get_population(qmdata_t *, qm_timestep_t *);
00331 
00332 /* Read the energy gradient for each atom. */
00333 static int get_gradient(qmdata_t *, qm_timestep_t *ts);
00334 
00335 /* Read energy gradient from final traj step. */
00336 static int get_final_gradient(qmdata_t *, qm_timestep_t *ts);
00337 
00338 /* Read ESP charges. */
00339 static int get_esp_charges(qmdata_t *data);
00340 
00341 /* For runtyp=HESSIAN, this subroutine scans the file for 
00342  * the hessian matrix in internal coordinates 
00343  * as well as the internal coordinate information */
00344 static int get_int_coords(qmdata_t *);
00345 
00346 /* get Hessian matrix in internal coordinates */
00347 static int get_int_hessian(qmdata_t *);
00348 
00349 /* For runtyp=HESSIAN, this subroutine scans the file for 
00350  * the cartesian hessian matrix */ 
00351 static int get_cart_hessian(qmdata_t *);
00352 
00353 /* For runtyp=HESSIAN, this subroutine reads the frequencies
00354  * and intensities of the normal modes */
00355 static int get_normal_modes(qmdata_t *);
00356 
00357 
00358 
00359 /* ######################################################## */
00360 /* Functions that are needed by the molfile_plugin          */
00361 /* interface to provide VMD with the parsed data            */
00362 /* ######################################################## */
00363 
00364 
00365 /***************************************************************
00366  *
00367  * Called by VMD to open the GAMESS logfile and get the number
00368  * of atoms.
00369  * We are also reading all the static (i.e. non-trajectory)
00370  * data here since we have to parse a bit to get the atom count
00371  * anyway. These data will then be provided to VMD by
00372  * read_gamess_metadata() and read_gamess_rundata().
00373  *
00374  * *************************************************************/
00375 static void *open_gamess_read(const char *filename, 
00376                   const char *filetype, int *natoms) {
00377 
00378   FILE *fd;
00379   qmdata_t *data = NULL;
00380   gmsdata *gms;
00381 
00382   /* open the input file */
00383   fd = fopen(filename, "rb");
00384  
00385   if (!fd) {
00386     PRINTERR;
00387     return NULL;
00388   }
00389 
00390   /* allocate and initialize main data structure */
00391   data = init_qmdata();
00392 
00393   /* make sure memory was allocated properly */
00394   if (data == NULL) {
00395     PRINTERR;
00396     return NULL;
00397   }
00398 
00399   /* allocate GAMESS specific data */
00400   gms = (gmsdata *)calloc(1,sizeof(gmsdata));
00401   data->format_specific_data = gms;
00402 
00403   gms->version = 0;
00404   gms->have_pcgamess = 0;
00405   gms->have_fmo = 0;
00406 
00407 
00408   /* store file pointer and filename in gamess struct */
00409   data->file = fd;
00410 
00411 
00412   /* check if the file is GAMESS format; if yes
00413    * parse it, if no exit */
00414   if (have_gamess(data, gms)==TRUE) {
00415     if (gms->have_pcgamess) {
00416       printf("gamessplugin) Warning: PC GAMESS/FIREFLY is not yet fully supported!\n");
00417     //  return NULL;
00418     }
00419     /* if we're dealing with an unsupported GAMESS
00420      * version, we better quit */
00421     if (gms->version==GAMESSUNKNOWN) {
00422       printf("gamessplugin) GAMESS version %s not supported. \n",
00423              data->version_string);
00424       return NULL;
00425     }
00426 
00427     /* get the non-trajectory information from the log file */    
00428     if (parse_static_data(data, natoms) == FALSE) 
00429       return NULL;
00430   }
00431   else {
00432     printf("gamessplugin) This seems to not be a GAMESS logfile.\n");
00433     return NULL;
00434   }
00435 
00436   return data;
00437 }
00438 
00439 
00440 /************************************************************
00441  * 
00442  * Provide VMD with the structure of the molecule, i.e the
00443  * atoms coordinates names, etc.
00444  *
00445  *************************************************************/
00446 static int read_gamess_structure(void *mydata, int *optflags, 
00447                       molfile_atom_t *atoms) 
00448 {
00449   qmdata_t *data = (qmdata_t *)mydata;
00450   qm_atom_t *cur_atom;
00451   molfile_atom_t *atom;
00452   int i = 0;
00453  
00454   /* optional atomic number provided */
00455   *optflags = MOLFILE_ATOMICNUMBER;
00456   /*  if (data->have_mulliken) 
00457    *optflags |= MOLFILE_CHARGE;*/
00458 
00459   /* all the information I need has already been read in
00460    * via the initial scan and I simply need to copy 
00461    * everything from the temporary arrays into the 
00462    * proper VMD arrays.
00463    * Since there are no atom names in the GAMESS output
00464    * I use the atom type here --- maybe there is a better
00465    * way to do this !!?? */
00466 
00467   /* get initial pointer for atom array */
00468   cur_atom = data->atoms;
00469 
00470   for(i=0; i<data->numatoms; i++) {
00471     atom = atoms+i;
00472     strncpy(atom->name, cur_atom->type, sizeof(atom->name)); 
00473     strncpy(atom->type, cur_atom->type, sizeof(atom->type));
00474     strncpy(atom->resname,"", sizeof(atom->resname)); 
00475     atom->resid = 1;
00476     atom->chain[0] = '\0';
00477     atom->segid[0] = '\0';
00478     atom->atomicnumber = cur_atom->atomicnum;
00479 #ifdef DEBUGGING
00480     printf("gamessplugin) atomicnum[%d] = %d\n", i, atom->atomicnumber);
00481 #endif
00482 
00483     /* if (data->have_mulliken)
00484       atom->charge = data->qm_timestep->mulliken_charges[i];
00485     */
00486     cur_atom++;
00487   }
00488  
00489   return MOLFILE_SUCCESS; 
00490 }
00491 
00492 
00493 /*****************************************************
00494  *
00495  * provide VMD with the sizes of the QM related
00496  * data structure arrays that need to be made
00497  * available
00498  *
00499  *****************************************************/
00500 static int read_gamess_metadata(void *mydata, 
00501     molfile_qm_metadata_t *metadata) {
00502 
00503   qmdata_t *data = (qmdata_t *)mydata;
00504 
00505   if (data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
00506     metadata->ncart = (3*data->numatoms);
00507     metadata->nimag = data->nimag;             
00508    
00509     if (data->have_internals) {
00510       metadata->nintcoords = data->nintcoords; 
00511     } else {
00512       metadata->nintcoords = 0;
00513     }
00514   }
00515   else {
00516     metadata->ncart = 0;
00517     metadata->nimag = 0;
00518     metadata->nintcoords = 0;
00519   }
00520 
00521   /* orbital data */
00522   metadata->num_basis_funcs = data->num_basis_funcs;
00523   metadata->num_basis_atoms = data->num_basis_atoms;
00524   metadata->num_shells      = data->num_shells;
00525   metadata->wavef_size      = data->wavef_size;  
00526 
00527 #if vmdplugin_ABIVERSION > 11
00528   /* system and run info */
00529   metadata->have_sysinfo = 1;
00530 
00531   /* hessian info */
00532   metadata->have_carthessian = data->have_cart_hessian;
00533   metadata->have_inthessian  = data->have_int_hessian;
00534 
00535   /* normal mode info */
00536   metadata->have_normalmodes = data->have_normal_modes;
00537 #endif
00538 
00539   return MOLFILE_SUCCESS;
00540 }
00541 
00542 
00543 /******************************************************
00544  * 
00545  * Provide VMD with the static (i.e. non-trajectory)
00546  * data. That means we are filling the molfile_plugin
00547  * data structures.
00548  *
00549  ******************************************************/
00550 static int read_gamess_rundata(void *mydata, 
00551                                molfile_qm_t *qm_data) {
00552 
00553   qmdata_t *data = (qmdata_t *)mydata;
00554   int i, j;
00555   int ncart;
00556   molfile_qm_hessian_t *hessian_data = &qm_data->hess;
00557   molfile_qm_basis_t   *basis_data   = &qm_data->basis;
00558   molfile_qm_sysinfo_t *sys_data     = &qm_data->run;
00559 
00560   /* fill in molfile_qm_hessian_t */
00561   if (data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
00562     ncart = 3*data->numatoms;
00563 
00564     /* Hessian matrix in cartesian coordinates */
00565     if (data->have_cart_hessian) {
00566       for (i=0; i<ncart; i++) {
00567         for (j=0; j<=i; j++) {
00568           hessian_data->carthessian[ncart*i+j] = data->carthessian[ncart*i+j];
00569           hessian_data->carthessian[ncart*j+i] = data->carthessian[ncart*i+j];
00570         }
00571       }
00572     }
00573 
00574     /* Hessian matrix in internal coordinates */
00575     if (data->have_int_hessian) {
00576       for (i=0; i<(data->nintcoords)*(data->nintcoords); i++) {
00577         hessian_data->inthessian[i] = data->inthessian[i];
00578       }
00579     }
00580 
00581     /* wavenumbers, intensities, normal modes */
00582     if (data->have_normal_modes) {
00583       for (i=0; i<ncart*ncart; i++) {
00584         hessian_data->normalmodes[i] = data->normal_modes[i];
00585       }
00586       for (i=0; i<ncart; i++) {
00587         hessian_data->wavenumbers[i] = data->wavenumbers[i];
00588         hessian_data->intensities[i] = data->intensities[i];
00589       }
00590     }
00591 
00592     /* imaginary modes */
00593     for (i=0; i<data->nimag; i++) {
00594       /*printf("imag_modes[%d]=%d\n", i, data->imag_modes[i]);*/
00595       hessian_data->imag_modes[i] = data->imag_modes[i];
00596     }
00597   }
00598 
00599   /* fill in molfile_qm_sysinfo_t */
00600   sys_data->runtype = data->runtype;
00601   sys_data->scftype = data->scftype;
00602   sys_data->nproc   = data->nproc;
00603   sys_data->num_electrons  = data->num_electrons;
00604   sys_data->totalcharge    = data->totalcharge;
00605   sys_data->num_occupied_A = data->num_occupied_A;
00606   sys_data->num_occupied_B = data->num_occupied_B;
00607   sys_data->status         = data->status;
00608 
00609 
00610   strncpy(sys_data->basis_string, data->basis_string,
00611           sizeof(sys_data->basis_string));
00612 
00613   sys_data->memory = 0; /* XXX fixme */
00614 
00615   strncpy(sys_data->runtitle, data->runtitle, sizeof(sys_data->runtitle));
00616   strncpy(sys_data->geometry, data->geometry, sizeof(sys_data->geometry));
00617   strncpy(sys_data->version_string, data->version_string,
00618           sizeof(sys_data->version_string));
00619 
00620 #if vmdplugin_ABIVERSION > 11
00621   /* fill in molfile_qm_basis_t */
00622   if (data->num_basis_funcs) {
00623     for (i=0; i<data->num_basis_atoms; i++) {
00624       basis_data->num_shells_per_atom[i] = data->num_shells_per_atom[i];
00625       basis_data->atomic_number[i] = data->atomicnum_per_basisatom[i];
00626     }
00627     
00628     for (i=0; i<data->num_shells; i++) {
00629       basis_data->num_prim_per_shell[i] = data->num_prim_per_shell[i];
00630       basis_data->shell_types[i] = data->shell_types[i];
00631     }
00632     
00633     for (i=0; i<2*data->num_basis_funcs; i++) {
00634       basis_data->basis[i] = data->basis[i];
00635     }
00636 
00637     for (i=0; i<3*data->wavef_size; i++) {
00638       basis_data->angular_momentum[i] = data->angular_momentum[i];
00639     }
00640   }
00641 #endif
00642  
00643   return MOLFILE_SUCCESS;
00644 }
00645 
00646 
00647 
00648 
00649 #if vmdplugin_ABIVERSION > 11
00650 
00651 /***********************************************************
00652  * Provide non-QM metadata for next timestep. 
00653  * Required by the plugin interface.
00654  ***********************************************************/
00655 static int read_timestep_metadata(void *mydata,
00656                                   molfile_timestep_metadata_t *meta) {
00657   meta->count = -1;
00658   meta->has_velocities = 0;
00659 
00660   return MOLFILE_SUCCESS;
00661 }
00662 
00663 /***********************************************************
00664  * Provide QM metadata for next timestep. 
00665  * This actually triggers reading the entire next timestep
00666  * since we have to parse the whole timestep anyway in order
00667  * to get the metadata. So we store the read data locally
00668  * and hand them to VMD when requested by read_timestep().
00669  *
00670  ***********************************************************/
00671 static int read_qm_timestep_metadata(void *mydata,
00672                                     molfile_qm_timestep_metadata_t *meta) {
00673   int have = 0;
00674 
00675   qmdata_t *data = (qmdata_t *)mydata;
00676 
00677   meta->count = -1; /* Don't know the number of frames yet */
00678 
00679   if (data->num_frames_read > data->num_frames_sent) {
00680     have = 1;
00681   }
00682   else if (data->num_frames_read < data->num_frames) {
00683     /*printf("gamessplugin) Probing timestep %d\n", data->num_frames_read);*/
00684 
00685     have = get_traj_frame(data, data->atoms, data->numatoms);
00686   }
00687 
00688   if (have) {
00689     int i;
00690     qm_timestep_t *cur_ts;
00691 
00692     /* get a pointer to the current qm timestep */
00693     cur_ts = data->qm_timestep+data->num_frames_sent;
00694     
00695     for (i=0; (i<MOLFILE_MAXWAVEPERTS && i<cur_ts->numwave); i++) {
00696       meta->num_orbitals_per_wavef[i] = cur_ts->wave[i].num_orbitals;
00697       meta->has_occup_per_wavef[i]    = cur_ts->wave[i].has_occup;
00698       meta->has_orben_per_wavef[i]    = cur_ts->wave[i].has_orben;
00699     }
00700     meta->wavef_size      = data->wavef_size;
00701     meta->num_wavef       = cur_ts->numwave;
00702     meta->num_scfiter     = cur_ts->num_scfiter;
00703     meta->num_charge_sets = cur_ts->have_mulliken +
00704       cur_ts->have_lowdin + cur_ts->have_esp;
00705     if (cur_ts->gradient) meta->has_gradient = TRUE;
00706 
00707   } else {
00708     meta->has_gradient = FALSE;
00709     meta->num_scfiter  = 0;
00710     meta->num_orbitals_per_wavef[0] = 0;
00711     meta->has_occup_per_wavef[0] = FALSE;
00712     meta->num_wavef = 0;
00713     meta->wavef_size = 0;
00714     meta->num_charge_sets = 0;
00715     data->trajectory_done = TRUE;
00716   }
00717 
00718   return MOLFILE_SUCCESS;
00719 }
00720 
00721 
00722 /***********************************************************
00723  *
00724  * This function provides the data of the next timestep.
00725  * Here we actually don't read the data from file, that had
00726  * to be done already upon calling read_timestep_metadata().
00727  * Instead we copy the stuff from the local data structure
00728  * into the one's provided by VMD.
00729  *
00730  ***********************************************************/
00731 static int read_timestep(void *mydata, int natoms, 
00732        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
00733                          molfile_qm_timestep_t *qm_ts) 
00734 {
00735   qmdata_t *data = (qmdata_t *)mydata;
00736   qm_timestep_t *cur_ts;
00737   int offset;
00738   int i = 0;
00739   int num_charge_sets = 0;
00740 
00741   if (data->trajectory_done == TRUE) return MOLFILE_ERROR;
00742 
00743 
00744   /* copy the coordinates */
00745   for (i=0; i<natoms; i++) {
00746     ts->coords[3*i  ] = data->atoms[i].x;
00747     ts->coords[3*i+1] = data->atoms[i].y;
00748     ts->coords[3*i+2] = data->atoms[i].z; 
00749   }    
00750   
00751   /* get a convenient pointer to the current qm timestep */
00752   cur_ts = data->qm_timestep+data->num_frames_sent;
00753 
00754   /* store the SCF energies */
00755   for (i=0; i<cur_ts->num_scfiter; i++) {
00756     qm_ts->scfenergies[i] = cur_ts->scfenergies[i];
00757   }
00758 
00759   /* store gradients */
00760   if (cur_ts->gradient) {
00761     for (i=0; i<3*natoms; i++) {
00762       qm_ts->gradient[i] = cur_ts->gradient[i];
00763     }
00764   }
00765 
00766   /* store charge sets*/
00767   if (cur_ts->have_mulliken) {
00768     offset = num_charge_sets*data->numatoms;
00769     for (i=0; i<data->numatoms; i++) {
00770       qm_ts->charges[offset+i] = cur_ts->mulliken_charges[i];
00771     }
00772     qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_MULLIKEN;
00773     num_charge_sets++;
00774   }
00775 
00776   if (cur_ts->have_lowdin) {
00777     offset = num_charge_sets*data->numatoms;
00778     for (i=0; i<data->numatoms; i++) {
00779       qm_ts->charges[offset+i] = cur_ts->lowdin_charges[i];
00780     }
00781     qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_LOWDIN;
00782     num_charge_sets++;
00783   }
00784   if (cur_ts->have_esp) {
00785     offset = num_charge_sets*data->numatoms;
00786     for (i=0; i<data->numatoms; i++) {
00787       qm_ts->charges[offset+i] = cur_ts->esp_charges[i];
00788     }
00789     qm_ts->charge_types[num_charge_sets] = MOLFILE_QMCHARGE_ESP;
00790     num_charge_sets++;
00791   }
00792 
00793 
00794   /* store the wave function and orbital energies */
00795   if (cur_ts->wave) {
00796     for (i=0; i<cur_ts->numwave; i++) {
00797       qm_wavefunction_t *wave = &cur_ts->wave[i];
00798       qm_ts->wave[i].type         = wave->type;
00799       qm_ts->wave[i].spin         = wave->spin;
00800       qm_ts->wave[i].excitation   = wave->exci;
00801       qm_ts->wave[i].multiplicity = wave->mult;
00802       qm_ts->wave[i].energy       = wave->energy;
00803       strncpy(qm_ts->wave[i].info, wave->info, MOLFILE_BUFSIZ);
00804 
00805       if (wave->wave_coeffs) {
00806         memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs,
00807                wave->num_orbitals*data->wavef_size*sizeof(float));
00808       }
00809       if (wave->orb_energies) {
00810         memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
00811                wave->num_orbitals*sizeof(float));
00812       }
00813       if (wave->has_occup) {
00814         memcpy(qm_ts->wave[i].occupancies, wave->orb_occupancies,
00815                wave->num_orbitals*sizeof(float));
00816       }
00817     }
00818   }
00819 
00820   if (data->runtype == MOLFILE_RUNTYPE_ENERGY || 
00821       data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
00822     /* We have only a single point */
00823     data->trajectory_done = TRUE;
00824   }
00825 
00826   data->num_frames_sent++;
00827 
00828   return MOLFILE_SUCCESS;
00829 }
00830 #endif
00831 
00832 
00833 
00834 /**********************************************************
00835  *
00836  * close file and free memory
00837  *
00838  **********************************************************/
00839 static void close_gamess_read(void *mydata) {
00840 
00841   qmdata_t *data = (qmdata_t *)mydata;
00842   int i, j;
00843   fclose(data->file);
00844 
00845   free(data->atoms);
00846   free(data->basis);
00847   free(data->shell_types);
00848   free(data->atomicnum_per_basisatom);
00849   free(data->num_shells_per_atom);
00850   free(data->num_prim_per_shell);
00851   free(data->bonds);
00852   free(data->angles);
00853   free(data->dihedrals);
00854   free(data->impropers);
00855   free(data->internal_coordinates);
00856   free(data->bond_force_const);
00857   free(data->angle_force_const);
00858   free(data->dihedral_force_const);
00859   free(data->improper_force_const);
00860   free(data->inthessian);
00861   free(data->carthessian);
00862   free(data->wavenumbers);
00863   free(data->intensities);
00864   free(data->normal_modes);
00865   free(data->imag_modes);
00866   free(data->angular_momentum);
00867   free(data->filepos_array);
00868 
00869   if (data->basis_set) {
00870     for(i=0; i<data->num_basis_atoms; i++) {
00871       for (j=0; j<data->basis_set[i].numshells; j++) {
00872         free(data->basis_set[i].shell[j].prim);
00873       }
00874       free(data->basis_set[i].shell);
00875     } 
00876     free(data->basis_set);
00877   }
00878 
00879   for (i=0; i<data->num_frames; i++) {
00880     free(data->qm_timestep[i].scfenergies);
00881     free(data->qm_timestep[i].gradient);
00882     free(data->qm_timestep[i].mulliken_charges);
00883     free(data->qm_timestep[i].lowdin_charges);
00884     free(data->qm_timestep[i].esp_charges);
00885     for (j=0; j<data->qm_timestep[i].numwave; j++) {
00886       free(data->qm_timestep[i].wave[j].wave_coeffs);
00887       free(data->qm_timestep[i].wave[j].orb_energies);
00888       free(data->qm_timestep[i].wave[j].orb_occupancies);
00889     }
00890     free(data->qm_timestep[i].wave);
00891   }
00892   free(data->qm_timestep);
00893   free(data->format_specific_data);
00894   free(data);
00895 }
00896 
00897 /* ####################################################### */
00898 /*             End of API functions                        */
00899 /* The following functions actually do the file parsing.   */
00900 /* ####################################################### */
00901 
00902 
00903 
00904 /********************************************************
00905  *
00906  * Main gamess log file parser responsible for static,  
00907  * i.e. non-trajectory information.
00908  *
00909  ********************************************************/
00910 static int parse_static_data(qmdata_t *data, int *natoms) 
00911 {
00912   /* Cast GAMESS specific data */
00913   gmsdata *gms = (gmsdata *)data->format_specific_data;
00914 
00915 
00916   /* Read # of procs and amount of requested memory */
00917   get_proc_mem(data, gms);
00918 
00919   /* Read the $BASIS options */
00920   if (!get_basis_options(data))   return FALSE;
00921 
00922   /* Read the run title */
00923   if (!get_runtitle(data))        return FALSE;
00924 
00925   /* Read the $CONTRL group;
00926    * It actually appears after the basis stats in the
00927    * output file but we jump ahead and read it here
00928    * because we need the units (ANGS/BOHR) before
00929    * reading the input structure. */
00930   if (gms->have_pcgamess){
00931     if (!get_contrl_firefly(data))          return FALSE;
00932   }
00933   else {
00934     if (!get_contrl(data))          return FALSE;
00935   }
00936 
00937   /* Read the input atom definitions and geometry */
00938   if (!get_input_structure(data, gms)) return FALSE;
00939 
00940   /* Read the basis set */
00941   if (!get_basis(data))           return FALSE; 
00942 
00943   /* Read the number of orbitals, electrons, 
00944    * charge, multiplicity, ... */
00945   if (!get_basis_stats(data))     return FALSE;
00946 
00947   /* Read input parameters regarding calculation of 
00948    * certain molecular properties such as electrostatic
00949    * moments and the MEP. */
00950   if (!get_properties_input(data)) return FALSE;
00951 
00952   /* Read internal coordinates */
00953   get_int_coords(data);
00954 
00955   /* Read symmetry point group and highest axis */
00956   if (!get_symmetry(data))         return FALSE;
00957 
00958   /* Read the $GUESS options */
00959   get_guess_options(data);
00960 
00961   if (data->scftype==MOLFILE_SCFTYPE_MCSCF) {
00962     if (!get_mcscf(data))          return FALSE;
00963   }
00964 
00965   /* Find the end of the trajectory and count the
00966    * frames on the way.
00967    * If no regular end is found we won't find any
00968    * properties to read and return. */
00969   if (!analyze_traj(data, gms)) {
00970     printf("gamessplugin) WARNING: Truncated or abnormally terminated file!\n\n");
00971   }
00972 
00973 
00974   /* provide VMD with the proper number of atoms */
00975   *natoms = data->numatoms;
00976 
00977   /* Read the first frame*/
00978   read_first_frame(data);
00979 
00980   /* Read the properties at the end of a calculation */
00981   get_final_properties(data);
00982 
00983 #ifdef DEBUGGING 
00984   printf("gamessplugin) num_frames_read = %d\n", data->num_frames_read);
00985   printf("gamessplugin) num_frames_sent = %d\n", data->num_frames_sent);
00986 
00987   /* Test print the parsed data in same format as logfile */
00988   print_input_data(data);
00989 #endif
00990 
00991   return TRUE;
00992 }
00993 
00994 
00995 #define TORF(x) (x ? "T" : "F")
00996 
00997 static void print_input_data(qmdata_t *data) {
00998   int i, j, k;
00999   int primcount=0;
01000   int shellcount=0;
01001 
01002   printf("\nDATA READ FROM FILE:\n\n");
01003   printf(" %10s WORDS OF MEMORY AVAILABLE\n", data->memory);
01004   printf("\n");
01005   printf("     BASIS OPTIONS\n");
01006   printf("     -------------\n");
01007   printf("%s\n", data->basis_string);
01008   printf("\n\n\n");
01009   printf("     RUN TITLE\n");
01010   printf("     ---------\n");
01011   printf(" %s\n", data->runtitle);
01012   printf("\n");
01013   printf(" THE POINT GROUP OF THE MOLECULE IS %s\n", "XXX");
01014   printf(" THE ORDER OF THE PRINCIPAL AXIS IS %5i\n", 0);
01015   printf("\n");
01016   printf(" YOUR FULLY SUBSTITUTED Z-MATRIX IS\n");
01017   printf("\n");
01018   printf(" THE MOMENTS OF INERTIA ARE (AMU-ANGSTROM**2)\n");
01019   printf(" IXX=%10.3f   IYY=%10.3f   IZZ=%10.3f\n", 0.0, 0.0, 0.0);
01020   printf("\n");
01021   printf(" ATOM      ATOMIC                      COORDINATES (BOHR)\n");
01022   printf("           CHARGE         X                   Y                   Z\n");
01023   for (i=0; i<data->numatoms; i++) {
01024     printf(" %-8s %6d", data->atoms[i].type, data->atoms[i].atomicnum);
01025     
01026     printf("%17.10f",   ANGS_TO_BOHR*data->atoms[i].x);
01027     printf("%20.10f",   ANGS_TO_BOHR*data->atoms[i].y);
01028     printf("%20.10f\n", ANGS_TO_BOHR*data->atoms[i].z);
01029   }
01030   printf("\n");
01031   printf("     ATOMIC BASIS SET\n");
01032   printf("     ----------------\n");
01033   printf(" THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN UNNORMALIZED\n");
01034   printf(" THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED TO UNITY\n");
01035   printf("\n");
01036   printf("  SHELL TYPE  PRIMITIVE        EXPONENT          CONTRACTION COEFFICIENT(S)\n");
01037   printf("\n");
01038 
01039 #if 0
01040   for (i=0; i<data->numatoms; i++) {
01041     printf("%-8s\n\n", data->atoms[i].type);
01042     printf("\n");
01043     printf("nshells=%d\n", data->num_shells_per_atom[i]);
01044 
01045     for (j=0; j<data->num_shells_per_atom[i]; j++) {
01046       printf("nprim=%d\n", data->num_prim_per_shell[shellcount]);
01047 
01048       for (k=0; k<data->num_prim_per_shell[shellcount]; k++) {
01049         printf("%6d   %d %7d %22f%22f\n", j, data->shell_types[shellcount],
01050                primcount+1, data->basis[2*primcount], data->basis[2*primcount+1]);
01051         primcount++;
01052       }
01053 
01054       printf("\n");
01055       shellcount++;
01056     }
01057   }
01058 #endif
01059   printf("gamessplugin) =================================================================\n");
01060   for (i=0; i<data->num_basis_atoms; i++) {
01061     printf("%-8s (%10s)\n\n", data->atoms[i].type, data->basis_set[i].name);
01062     printf("\n");
01063 
01064     for (j=0; j<data->basis_set[i].numshells; j++) {
01065 
01066       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
01067         printf("%6d   %d %7d %22f%22f\n", j,
01068                data->basis_set[i].shell[j].type,
01069                primcount+1,
01070                data->basis_set[i].shell[j].prim[k].exponent,
01071                data->basis_set[i].shell[j].prim[k].contraction_coeff);
01072         primcount++;
01073       }
01074 
01075       printf("\n");
01076       shellcount++;
01077     }
01078   }
01079   printf("\n");
01080   printf(" TOTAL NUMBER OF BASIS SET SHELLS             =%5d\n", data->num_shells);
01081   printf(" NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =%5d\n", data->wavef_size);
01082   printf(" NUMBER OF ELECTRONS                          =%5d\n", data->num_electrons);
01083   printf(" CHARGE OF MOLECULE                           =%5d\n", data->totalcharge);
01084   printf(" SPIN MULTIPLICITY                            =%5d\n", data->multiplicity);
01085   printf(" NUMBER OF OCCUPIED ORBITALS (ALPHA)          =%5d\n", data->num_occupied_A);
01086   printf(" NUMBER OF OCCUPIED ORBITALS (BETA )          =%5d\n", data->num_occupied_B);
01087   printf(" TOTAL NUMBER OF ATOMS                        =%5i\n", data->numatoms);
01088   printf("\n");
01089 }
01090 
01091 
01092 
01093 /**********************************************************
01094  *
01095  * this subroutine checks if the provided files is
01096  * actually a GAMESS file;
01097  *
01098  **********************************************************/
01099 static int have_gamess(qmdata_t *data, gmsdata *gms) 
01100 {
01101   char word[3][BUFSIZ];
01102   char buffer[BUFSIZ];
01103   char versionstr[BUFSIZ];
01104   int day, year;
01105   char month[BUFSIZ], rev[BUFSIZ];
01106   int i = 0;
01107   int program;
01108   int ver,build;
01109   buffer[0] = '\0';
01110   for (i=0; i<3; i++) word[i][0] = '\0';
01111 
01112 
01113   /* check if the file is GAMESS format */
01114   program = goto_keyline(data->file,
01115                           "PC GAMESS version",
01116                           "GAMESS VERSION =", 
01117                           "Firefly version",NULL);
01118   if (program==1) {
01119     gms->have_pcgamess = 1;
01120     gms->version = 1;
01121     strcpy(data->version_string, "PC GAMESS ");
01122   } else if (program==2) {
01123     gms->have_pcgamess = 0;
01124     strcpy(data->version_string, "GAMESS ");
01125   } else if (program==3) {
01126     gms->have_pcgamess = 1;
01127     gms->version = FIREFLY8PRE6695;
01128     strcpy(data->version_string, "Firefly ");
01129   } else {
01130     printf("gamessplugin) This is no GAMESS/PCGAMESS/Firefly logfile!\n");
01131     return FALSE;
01132   }
01133 
01134   GET_LINE(buffer, data->file);
01135 
01136   if (gms->have_pcgamess) {
01137     if (strstr(buffer,"version") != NULL) {
01138       strncpy(versionstr, strstr(buffer,"version")+8, 16);
01139       *strchr(versionstr, ' ') = '\0';
01140       sscanf(buffer, "%*s %*s %*s %*s %*s %*s %d", &build);
01141       sscanf(versionstr, "%1d%*s", &ver);
01142       printf("gamessplugin) Firefly build = %d %d\n", 
01143          ver,build);
01144       if (ver >= 8 && build >= 6695)
01145         gms->version = FIREFLY8POST6695;
01146       else
01147         gms->version = FIREFLY8PRE6695;
01148     }
01149   } else {
01150     /* extract the version number if possible; otherwise
01151      * return empty string */
01152     if (strstr(buffer,"=") != NULL) {
01153       strncpy(versionstr, strstr(buffer,"=")+2, 16);
01154       versionstr[16] = '\0';
01155     }
01156     
01157     /* determine if we're dealing with pre-"27 JUN 2005"
01158      * version */
01159     sscanf(versionstr, "%d %s %d %s", &day, month, &year, rev);
01160     
01161     if ( ( year >= 2006 ) ||
01162          ( year == 2005 && !strcmp(month,"JUN") ) ||
01163          ( year == 2005 && !strcmp(month,"NOV") ) ||
01164          ( year == 2005 && !strcmp(month,"DEC") ) )
01165       {
01166         gms->version = GAMESSPOST20050627;
01167       } else { 
01168         gms->version = GAMESSPRE20050627;
01169       }
01170   }
01171 
01172   strcat(data->version_string, versionstr);
01173 
01174   printf("gamessplugin) Version = %s\n", 
01175          data->version_string);
01176 
01177   return TRUE;
01178 }
01179 
01180 
01181 /**********************************************************
01182  *
01183  * this subroutine reads the number of procs and the amount
01184  * of memory requested
01185  *
01186  **********************************************************/
01187 static int get_proc_mem(qmdata_t *data, gmsdata *gms) {
01188 
01189   char word[4][BUFSIZ];
01190   char buffer[BUFSIZ];
01191   char *temp;
01192   int nproc;
01193   int i;
01194 
01195   buffer[0] = '\0';
01196   for (i=0; i<3; i++) word[i][0] = '\0';
01197 
01198   rewind(data->file);
01199 
01200   
01201   /* scan for the number of processors; here we need
01202    * distinguish between vanilla Gamess and PC Gamess */
01203   if (gms->have_pcgamess == 1) {
01204     /* XXX for now we fake ncpu = 1 until we know exactly
01205      *     how the output format looks like */
01206     nproc = 1;
01207     do {
01208       GET_LINE(buffer, data->file);
01209       sscanf(buffer,"%s %d %s",&word[0][0],&nproc,&word[1][0]);
01210         if (!strcmp(&word[0][0],"PARALLEL") &&
01211                !strcmp(&word[0][1],"RUNNING")) {
01212             sscanf(buffer,"%*s %*s %*s %*s %*s %d %*s %*s",&nproc);
01213             break;
01214         }
01215       } while (strcmp(&word[0][0],"ECHO") || 
01216                strcmp(&word[1][0],"THE") );
01217 
01218   }
01219   else {
01220     do {
01221       GET_LINE(buffer, data->file);
01222       sscanf(buffer,"%s %d %s",&word[0][0],&nproc,&word[1][0]);
01223 
01224       if (!strcmp(&word[0][0],"Initiating") &&
01225           !strcmp(&word[1][0],"compute")) {
01226         break;
01227       }
01228 
01229       else if (!strcmp(&word[0][0],"Initiating") &&
01230                !strcmp(&word[1][0],"processes")) {
01231         break;
01232       }
01233 
01234       /* Some versions */
01235       else if (!strcmp(&word[0][0],"PARALLEL") &&
01236                !strcmp(&word[0][1],"RUNNING")) {
01237         sscanf(buffer,"%*s %*s %*s %*s %d %*s",&nproc);
01238         break;
01239       }
01240 
01241     } while (strcmp(&word[0][0],"ECHO") || 
01242              strcmp(&word[1][0],"THE") );
01243   }
01244   
01245   /* store the number of processors */
01246   data->nproc = nproc;
01247 
01248   
01249   /* scan for the amount of memory requested */
01250   do {
01251     GET_LINE(buffer, data->file);
01252     sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);
01253 
01254   } while( strcmp(&word[0][0],"$SYSTEM") || 
01255            strcmp(&word[1][0],"OPTIONS") );
01256 
01257   eatline(data->file, 1);
01258 
01259 
01260   /* next line contains the amount of memory requested,
01261    * vanilla Gamess and PC Gamess need separate treatment */
01262   if (gms->have_pcgamess == 1) {
01263     GET_LINE(buffer, data->file);
01264 
01265     /* store it */
01266     if ((temp = strstr(buffer,"MEMORY=")+8)==NULL) return FALSE;
01267     strncpy(data->memory,trimright(temp),sizeof(data->memory));
01268   }
01269   else {
01270     GET_LINE(buffer, data->file);
01271     sscanf(buffer,"%s %s %s",&word[0][0],&word[1][0],&word[2][0]);
01272 
01273     /* store it */
01274     strncpy(data->memory,&word[2][0],sizeof(data->memory));
01275   }
01276 
01277   printf("gamessplugin) GAMESS used %d compute processes \n", nproc);
01278   printf("gamessplugin) GAMESS used %s words of memory \n", data->memory);
01279 
01280   return TRUE;
01281 }
01282 
01283 
01284 /**********************************************************
01285  *
01286  * Extract the $BASIS options
01287  *
01288  **********************************************************/
01289 static int get_basis_options(qmdata_t *data) {
01290 
01291   char buffer[BUFSIZ];
01292   char diffuse[BUFSIZ];
01293   char polarization[BUFSIZ];
01294   int ngauss;
01295 
01296   buffer[0] = '\0';
01297   diffuse[0] = '\0';
01298   polarization[0] = '\0';
01299 
01300   /* The $BASIS section is somewhere at the beginning of
01301    * the file. Rewind to be sure not to miss it. */
01302   rewind(data->file);
01303 
01304   /* start scanning */
01305   if (pass_keyline(data->file, "BASIS OPTIONS",
01306                     "RUN TITLE") != FOUND) {
01307     /* No Basis options section found
01308      * (basis was entered explicitly) */
01309     return TRUE;
01310   }
01311 
01312   eatline(data->file, 1);
01313 
01314 
01315   /* the first string in the current line contains the
01316    * GBASIS used; copy it over into the gbasis variable
01317    * of qmdata_t */
01318   GET_LINE(buffer, data->file);
01319   sscanf(buffer," GBASIS=%s IGAUSS= %d", data->gbasis, &ngauss);
01320  
01321 
01322   /* in case we're using a pople style basis set, i.e. 
01323    * GBASIS=N311,N31,N21 or STO we also scan for the number 
01324    * of gaussians, as well as p,d,f and diffuse functions
01325    * and use this info to assemble a "basis set string" */
01326   if ( !strncmp(data->gbasis,"N311",sizeof(data->gbasis)) ||
01327        !strncmp(data->gbasis,"N31",sizeof(data->gbasis)) ||
01328        !strncmp(data->gbasis,"N21",sizeof(data->gbasis)) ||
01329        !strncmp(data->gbasis,"STO",sizeof(data->gbasis)) ) 
01330   {
01331     int npfunc, ndfunc, nffunc;
01332     int diffs=FALSE, diffsp=FALSE;
01333     char torf;
01334 
01335     /* the next line gives us the d,f and diffuse sp
01336      * functions */
01337     GET_LINE(buffer, data->file);
01338     if (sscanf(buffer," NDFUNC= %d NFFUNC= %d DIFFSP= %c",
01339                &ndfunc, &nffunc, &torf) != 3) {
01340       sscanf(buffer," NDFUNC= %d DIFFSP= %c", &ndfunc, &torf);
01341     }
01342 
01343     /* convert GAMESS' .TRUE./.FALSE. for DIFFSP into 1/0 */
01344     if (torf=='T') diffsp = TRUE;
01345 
01346 
01347     /* the next line gives us the p and diffuse s
01348      * functions */
01349     GET_LINE(buffer, data->file);
01350     sscanf(buffer," NPFUNC= %d DIFFS= %c", &npfunc, &torf);
01351 
01352     /* convert GAMESS' .TRUE./.FALSE. for DIFFSP into 1/0 */
01353     if (torf=='T') diffs = TRUE;
01354 
01355 
01356     /* now we need some logic to assemble this info into
01357      * some nice looking string a la "6-31G*" */
01358 
01359     /* create the diffuse function string */
01360     if (diffs && diffsp) {
01361         strncpy(diffuse,"++",sizeof(diffuse));
01362     }
01363     else if (diffsp) {
01364         strncpy(diffuse,"+",sizeof(diffuse));
01365     }
01366     else {
01367         strncpy(diffuse,"",sizeof(diffuse));
01368     }
01369 
01370 
01371     /* create the polarization function string */
01372     if (npfunc>0 && ndfunc>0 && nffunc>0) {
01373       sprintf(polarization, "(%dp,%dd,%df)", npfunc, ndfunc, nffunc);
01374     } else if (npfunc>0 && ndfunc>0) {
01375       sprintf(polarization, "(%dp,%dd)", npfunc, ndfunc);
01376     }
01377     else if (npfunc>0) {
01378       sprintf(polarization, "(%dp)", npfunc);
01379     }
01380     else if (ndfunc>0) {
01381       sprintf(polarization, "(%dd)", ndfunc);
01382     } 
01383     else {
01384       strncpy(polarization, "", sizeof(polarization));
01385     } 
01386 
01387     /* assemble the bits */ 
01388     if (!strcmp(data->gbasis, "STO")) {
01389       sprintf(data->basis_string, "STO-%dG%s%s",
01390               ngauss, diffuse, polarization);
01391     }
01392     else {
01393       sprintf(data->basis_string, "%d-%s%sG%s",
01394               ngauss, (data->gbasis+1), diffuse, 
01395               polarization);
01396     }      
01397   }
01398 
01399   /* cc-pVnZ and cc-pCVnZ */
01400   else if (!strncmp(data->gbasis, "CC",  2)) {
01401     strcpy(data->basis_string, "cc-p");
01402     if (strlen(data->gbasis)==4 && data->gbasis[3]=='C') {
01403       strcat(data->basis_string, "C");
01404     }
01405     strcat(data->basis_string, "V");
01406     strncat(data->basis_string, &data->gbasis[2], 1);
01407     strcat(data->basis_string, "Z");
01408   }
01409 
01410   /* aug-cc-pVnZ and aug-cc-pCVnZ */
01411   else if (!strncmp(data->gbasis, "ACC", 3)) {
01412     strcpy(data->basis_string, "aug-cc-p");
01413     if (strlen(data->gbasis)==5 && data->gbasis[4]=='C') {
01414       strcat(data->basis_string, "C");
01415     }
01416     strcat(data->basis_string, "V");
01417     strncat(data->basis_string, &data->gbasis[3], 1);
01418     strcat(data->basis_string, "Z");
01419   }
01420 
01421   /* for non Pople style basis sets we just use the GBASIS
01422    * for the basis string;
01423    * TODO: make the basis_string more comprehensive for non
01424    *       pople-style basis sets */
01425   else {
01426     strncpy(data->basis_string,data->gbasis,
01427             sizeof(data->basis_string));
01428   }
01429 
01430   return TRUE;
01431 }
01432 
01433 
01434 
01435 /**********************************************************
01436  *
01437  * Extract the run title line
01438  *
01439  **********************************************************/
01440 static int get_runtitle(qmdata_t *data) {
01441 
01442   char buffer[BUFSIZ];
01443 
01444   if (pass_keyline(data->file, "RUN TITLE",
01445                    "THE POINT GROUP") != FOUND) {
01446     /* This is most likely a broken file, but who knows. 
01447      * Since we don't really care about the title string
01448      * we go on here without error. */
01449     data->runtitle[0] = '\0';
01450     return TRUE;
01451   }
01452 
01453   GET_LINE(buffer, data->file);
01454   strncpy(data->runtitle,trimright(buffer),sizeof(data->runtitle));
01455 
01456   return TRUE;
01457 } 
01458 
01459 
01460 /**********************************************************
01461  *
01462  * Read the input atom definitions and geometry
01463  *
01464  **********************************************************/
01465 static int get_input_structure(qmdata_t *data, gmsdata *gms) {
01466   char buffer[BUFSIZ];
01467   char units[BUFSIZ];
01468   int numatoms = -1;
01469   int bohr;
01470   long filepos;
01471   filepos = ftell(data->file);
01472 
01473   /* See if we find the "ATOM      ATOMIC ..." line before 
01474    * any of the three stopstrings mrking the beginning of
01475    * possible following sections. */
01476   if (goto_keyline(data->file,
01477          "ATOM      ATOMIC                      COORDINATES (",
01478          "INTERNUCLEAR DISTANCES",
01479          "ATOMIC BASIS SET",
01480          "$CONTRL OPTIONS", NULL) == FOUND) {
01481 
01482     GET_LINE(buffer, data->file);
01483     sscanf(buffer, " ATOM      ATOMIC  %*s  %s", units);
01484     eatline(data->file, 1);
01485 
01486   } else {
01487     /* This is probably an FMO calc.; if so, set flag. */
01488     fseek(data->file, filepos, SEEK_SET);
01489     if (pass_keyline(data->file,
01490                      "The Fragment Molecular Orbital (FMO) method.", 
01491                      NULL)) {
01492       gms->have_fmo = 1;
01493       printf("gamessplugin) Fragment Molecular Orbital (FMO) method.\n");
01494     }
01495 
01496     /* We didn't find the normal input section.
01497      * Let's see i we can find coordinates for the first
01498      * frame of a trajectory. */
01499     fseek(data->file, filepos, SEEK_SET);
01500     if (pass_keyline(data->file,
01501                       "BEGINNING GEOMETRY SEARCH POINT NSERCH=   0",
01502                       NULL) &&
01503         goto_keyline(data->file, "COORDINATES OF ALL ATOMS", NULL)) {
01504       GET_LINE(buffer, data->file);
01505       sscanf(buffer, " COORDINATES OF ALL ATOMS ARE %s", units);
01506       eatline(data->file, 2);
01507 
01508     } else {
01509       /* As last resort look for FMO coordinates in the
01510        * INPUT CARD section: */
01511 
01512       /* But first we have to get the units from the $CONTRL section */
01513       rewind(data->file);
01514       if (!pass_keyline(data->file, "$CONTRL OPTIONS", NULL)) {
01515         printf("gamessplugin) Missing $CONTRL OPTIONS section!\n");
01516         return FALSE;
01517       }
01518       goto_keyline(data->file, "UNITS =", NULL);
01519       GET_LINE(buffer, data->file);
01520       sscanf(strstr(buffer, "UNITS ="), "%s", units);
01521       bohr = !strcmp(units, "BOHR");
01522       
01523       /* Find beginning of $FMOXYZ input card */
01524       rewind(data->file);
01525       if (!pass_keyline(data->file, "INPUT CARD> $fmoxyz", 
01526                         "INPUT CARD> $FMOXYZ")) {
01527         printf("gamessplugin) No atom coordinates found!\n");
01528         return FALSE;
01529       }
01530             
01531       /* Read the $FMOXYZ coordinates */     
01532       if (!get_fmoxyz(data->file, &data->atoms, bohr, &numatoms)) {
01533         printf("gamessplugin) Could not read coordinates from $FMOXYZ input!\n");
01534         return FALSE;
01535       } else {
01536         printf("gamessplugin) Fragment Molecular Orbital (FMO) method.\n");
01537         gms->have_fmo = 1;
01538         data->numatoms = numatoms;
01539         return TRUE;
01540       }
01541     }
01542   }
01543 
01544   /* If we reached this point we have found either a regular
01545    * input coordinate block or the coordinate block of the
01546    * first trajectory frame. */
01547 
01548   /* test if coordinate units are Bohr */
01549   bohr = !strcmp(units, "(BOHR)");
01550 
01551   /* Read the coordinate block */
01552   if (get_coordinates(data->file, &data->atoms, bohr, &numatoms))
01553     data->num_frames_read = 0;
01554   else {
01555     printf("gamessplugin) Bad atom coordinate block!\n");
01556     return FALSE;
01557   }
01558 
01559   fseek(data->file, filepos, SEEK_SET);
01560 
01561   /* store number of atoms in data structure */
01562   data->numatoms = numatoms;
01563 
01564   return TRUE; 
01565 }
01566 
01567 
01568 /**********************************************************
01569  *
01570  * Read an atom coordinate block.
01571  *
01572  * Example:
01573  *  F         9.0    3.04259     -0.07605       0.00000
01574  *  N         7.0    0.03017      0.38347       0.00000
01575  * ...
01576  *
01577  **********************************************************/
01578 static int get_coordinates(FILE *file, qm_atom_t **atoms, int unit,
01579                            int *numatoms) {
01580   int i = 0;
01581   int growarray = 0;
01582 
01583   if (*numatoms<0) {
01584     *atoms = (qm_atom_t*)calloc(1, sizeof(qm_atom_t));
01585     growarray = 1;
01586   }
01587 
01588   /* Read in the coordinates until an empty line is reached.
01589    * We expect 5 entries per line */
01590   while (1) {
01591     char buffer[BUFSIZ];
01592     char atname[BUFSIZ];
01593     float atomicnum;
01594     float x,y,z, dum;
01595     int n;
01596     qm_atom_t *atm;
01597 
01598     GET_LINE(buffer, file);
01599 
01600     /* For FMO there is an additional atom index in the
01601      * second column. Try both variants: */
01602     n = sscanf(buffer,"%s %f %f %f %f %f",atname,&dum,&atomicnum,&x,&y,&z);
01603     if (n!=6) {
01604       n = sscanf(buffer,"%s %f %f %f %f",atname,&atomicnum,&x,&y,&z);
01605     }
01606     if (n!=5 && n!=6) break;
01607 
01608     if (growarray && i>0) {
01609       *atoms = (qm_atom_t*)realloc(*atoms, (i+1)*sizeof(qm_atom_t));
01610     }
01611     atm = (*atoms)+i;
01612 
01613     strncpy(atm->type, atname, sizeof(atm->type));
01614     atm->atomicnum = floor(atomicnum+0.5); /* nuclear charge */
01615     /*printf("coor: %s %d %f %f %f\n", atm->type, atm->atomicnum, x, y, z);*/
01616    
01617     /* if coordinates are in Bohr convert them to Angstrom */
01618     if (unit==BOHR) {
01619       x *= BOHR_TO_ANGS;
01620       y *= BOHR_TO_ANGS;
01621       z *= BOHR_TO_ANGS;
01622     }
01623     
01624     atm->x = x;
01625     atm->y = y;
01626     atm->z = z; 
01627     i++;
01628   }
01629 
01630   /* If file is broken off in the middle of the coordinate block 
01631    * we cannot use this frame. */
01632   if (*numatoms>=0 && *numatoms!=i) {
01633     (*numatoms) = i;
01634     return FALSE;
01635   }
01636 
01637   (*numatoms) = i;
01638   return TRUE;
01639 }
01640 
01641 
01642 /* Read coordinates from $FMOXYZ section in the INPUT CARD
01643  * listing at the beginning of the file. This is a method
01644  * of last resort used only for FMO calculations where no
01645  * coordinates are printed. */
01646 static int get_fmoxyz(FILE *file, qm_atom_t **atoms, int unit,
01647                            int *numatoms) {
01648   int i = 0;
01649   int growarray = 0;
01650 
01651   if (*numatoms<0) {
01652     *atoms = (qm_atom_t*)calloc(1, sizeof(qm_atom_t));
01653     growarray = 1;
01654   }
01655 
01656   /* Read in the coordinates until an empty line is reached.
01657    * We expect 5 entries per line */
01658   while (1) {
01659     char buffer[BUFSIZ];
01660     char atname[BUFSIZ], element[BUFSIZ];
01661     float x,y,z;
01662     int n;
01663     qm_atom_t *atm;
01664 
01665     GET_LINE(buffer, file);
01666 
01667     /* skip " INPUT CARD>" at the beginning of the line */
01668     n = sscanf(buffer+12,"%s %s %f %f %f",atname,element,&x,&y,&z);
01669 
01670     if (n!=5) break;
01671 
01672     if (growarray && i>0) {
01673       *atoms = (qm_atom_t*)realloc(*atoms, (i+1)*sizeof(qm_atom_t));
01674     }
01675     atm = (*atoms)+i;
01676 
01677     strncpy(atm->type, atname, sizeof(atm->type));
01678     if (isalpha(element[0]))
01679       atm->atomicnum = get_pte_idx_from_string(element);
01680     else if (isdigit(element[0])) 
01681       atm->atomicnum = floor(element[0]+0.5); /* nuclear charge */
01682     else break;
01683 
01684     /* If coordinates are in Bohr convert them to Angstrom */
01685     if (unit==BOHR) {
01686       x *= BOHR_TO_ANGS;
01687       y *= BOHR_TO_ANGS;
01688       z *= BOHR_TO_ANGS;
01689     }
01690     
01691     atm->x = x;
01692     atm->y = y;
01693     atm->z = z; 
01694     i++;
01695   }
01696 
01697   /* If file is broken off in the middle of the coordinate block 
01698    * we cannot use this frame. */
01699   if (*numatoms>=0 && *numatoms!=i) return FALSE;
01700 
01701   (*numatoms) = i;
01702   return TRUE;
01703 }
01704 
01705 
01706 /**********************************************************
01707  *
01708  * Read data from the $CONTRL card from FIREFLY
01709  *
01710  **********************************************************/
01711 static int get_contrl_firefly(qmdata_t *data) {
01712 
01713   char word[3][BUFSIZ];
01714   char buffer[BUFSIZ];
01715   char *temp;
01716   long filepos;
01717   filepos = ftell(data->file);
01718 
01719   word[0][0] = '\0';
01720   word[1][0] = '\0';
01721   word[2][0] = '\0';
01722   buffer[0] = '\0';
01723 
01724   //printf("gamessplugin) Getting CONTRL group data for Firefly \n");
01725 
01726   /* start scanning; currently we support
01727    * RUNTYP = ENERGY, OPTIMIZE, SADPOINT, HESSIAN, SURFACE */
01728   if (!pass_keyline(data->file, "$CONTRL OPTIONS", NULL)) {
01729     fseek(data->file, filepos, SEEK_SET);
01730     return FALSE;
01731   }
01732 
01733   eatline(data->file, 1);
01734 
01735   /* current line contains SCFTYP, RUNTYP, EXETYP info; scan it */
01736   GET_LINE(buffer, data->file);
01737   sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);
01738 
01739   /* check for supported RUNTYPs */
01740   if      (!strcmp(&word[1][0],"RUNTYP=ENERGY")) {
01741     data->runtype = MOLFILE_RUNTYPE_ENERGY;
01742   }
01743   else if (!strcmp(&word[1][0],"RUNTYP=OPTIMIZE")) {
01744     data->runtype = MOLFILE_RUNTYPE_OPTIMIZE;
01745   }
01746   else if (!strcmp(&word[1][0],"RUNTYP=SADPOINT")) {
01747     data->runtype = MOLFILE_RUNTYPE_SADPOINT;
01748   }
01749   else if (!strcmp(&word[1][0],"RUNTYP=HESSIAN")) {
01750     data->runtype = MOLFILE_RUNTYPE_HESSIAN;
01751   }
01752   else if (!strcmp(&word[1][0],"RUNTYP=SURFACE")) {
01753     data->runtype = MOLFILE_RUNTYPE_SURFACE;
01754   }
01755   else if (!strcmp(&word[1][0],"RUNTYP=GRADIENT")) {
01756     data->runtype = MOLFILE_RUNTYPE_GRADIENT;
01757   }
01758   else if (!strcmp(&word[1][0],"RUNTYP=MEX")) {
01759     data->runtype = MOLFILE_RUNTYPE_MEX;
01760   }
01761   else {
01762 #ifdef DEBUGGING
01763     printf("gamessplugin) The %s is currently not supported \n",
01764            &word[1][0]);
01765 #endif
01766     data->runtype = MOLFILE_RUNTYPE_UNKNOWN;
01767   }
01768   printf("gamessplugin) File generated via %s \n",&word[1][0]);
01769 
01770 
01771   /* determine SCFTYP */
01772   if (!strcmp(&word[0][0],"SCFTYP=RHF")) {
01773     data->scftype = MOLFILE_SCFTYPE_RHF;
01774   }
01775   else if (!strcmp(&word[0][0],"SCFTYP=UHF")) {
01776     data->scftype = MOLFILE_SCFTYPE_UHF;
01777   }
01778   else if (!strcmp(&word[0][0],"SCFTYP=ROHF")) {
01779     data->scftype = MOLFILE_SCFTYPE_ROHF;
01780   }
01781   else if (!strcmp(&word[0][0],"SCFTYP=GVB")) {
01782     data->scftype = MOLFILE_SCFTYPE_GVB;
01783   }
01784   else if (!strcmp(&word[0][0],"SCFTYP=MCSCF")) {
01785     data->scftype = MOLFILE_SCFTYPE_MCSCF;
01786   }
01787   else if (!strcmp(&word[0][0],"SCFTYP=NONE")) {
01788     data->scftype = MOLFILE_SCFTYPE_NONE;
01789   }
01790   else {
01791     /* if we don't find a supported SCFTYP we bomb out; this
01792      * might be a little drastic */
01793     printf("gamessplugin) %s is currently not supported \n",
01794            &word[0][0]);
01795     return FALSE;
01796   }
01797   printf("gamessplugin) Type of wavefunction used %s \n",
01798          &word[0][0]);
01799 
01800   /* scan for MPLEVL, LOCAL and UNITS; */
01801   GET_LINE(buffer, data->file);
01802   sscanf(buffer,"%s %s %*s %s",&word[0][0],&word[1][0],&word[2][0]);
01803   data->mplevel = atoi(&word[1][0]);
01804 
01805   /* scan for MULT, ICHARG and MAXIT; */
01806   GET_LINE(buffer, data->file);
01807   //sscanf(buffer,"%*s %*s %*s %*s %s %s",&word[0][0],&word[1][0]);
01808 
01809   /* find the coordinate type in next line */
01810   while ( (temp=strstr(buffer,"COORD =")) == NULL ) {
01811     GET_LINE(buffer, data->file);;
01812   }
01813   strncpy(data->geometry, trimright(temp+7), sizeof(data->geometry)); 
01814   printf("gamessplugin) Coordinate type used is %s \n", data->geometry);
01815 
01816   while ( (temp=strstr(buffer,"CITYP =")) == NULL ) {
01817     GET_LINE(buffer, data->file);;
01818   }
01819   strncpy(buffer, trimright(temp+7), 8); 
01820 
01821   /* determine CITYP */
01822   if      (!strcmp(buffer,"NONE"))  data->citype = CI_NONE;
01823   else if (!strcmp(buffer,"CIS"))   data->citype = CI_CIS;
01824   else if (!strcmp(buffer,"ALDET")) data->citype = CI_ALDET;
01825   else if (!strcmp(buffer,"ORMAS")) data->citype = CI_ORMAS;
01826   else if (!strcmp(buffer,"GUGA"))  data->citype = CI_GUGA;
01827   else if (!strcmp(buffer,"FSOCI")) data->citype = CI_FSOCI;
01828   else if (!strcmp(buffer,"GENCI")) data->citype = CI_GENCI;
01829   else                                    data->citype = CI_UNKNOWN;
01830   printf("gamessplugin) CI method %s \n",buffer);
01831 
01832   GET_LINE(buffer, data->file);
01833   sscanf(buffer,"%s %*s",&word[0][0]);
01834 
01835   /* scan for DFTTYP, TDDFT info; */
01836   if (!strncmp(&word[0][0],"DFTTYP=", 7)) {
01837     printf("gamessplugin) Density functional used is %s \n",&word[0][7]);
01838     GET_LINE(buffer, data->file);
01839   }
01840 
01841 
01842 
01843   fseek(data->file, filepos, SEEK_SET);
01844   return TRUE;
01845 }
01846 
01847 
01848 /**********************************************************
01849  *
01850  * Read data from the $CONTRL card from GAMESS
01851  *
01852  **********************************************************/
01853 static int get_contrl(qmdata_t *data) {
01854 
01855   char word[3][BUFSIZ];
01856   char buffer[BUFSIZ];
01857   char *temp;
01858   long filepos;
01859   filepos = ftell(data->file);
01860 
01861   word[0][0] = '\0';
01862   word[1][0] = '\0';
01863   word[2][0] = '\0';
01864   buffer[0] = '\0';
01865 
01866 
01867   /* start scanning; currently we support
01868    * RUNTYP = ENERGY, OPTIMIZE, SADPOINT, HESSIAN, SURFACE */
01869   if (!pass_keyline(data->file, "$CONTRL OPTIONS", NULL)) {
01870     fseek(data->file, filepos, SEEK_SET);
01871     return FALSE;
01872   }
01873 
01874   eatline(data->file, 1);
01875 
01876   /* current line contains SCFTYP, RUNTYP, EXETYP info; scan it */
01877   GET_LINE(buffer, data->file);
01878   sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);
01879 
01880   /* check for supported RUNTYPs */
01881   if      (!strcmp(&word[1][0],"RUNTYP=ENERGY")) {
01882     data->runtype = MOLFILE_RUNTYPE_ENERGY;
01883   }
01884   else if (!strcmp(&word[1][0],"RUNTYP=OPTIMIZE")) {
01885     data->runtype = MOLFILE_RUNTYPE_OPTIMIZE;
01886   }
01887   else if (!strcmp(&word[1][0],"RUNTYP=SADPOINT")) {
01888     data->runtype = MOLFILE_RUNTYPE_SADPOINT;
01889   }
01890   else if (!strcmp(&word[1][0],"RUNTYP=HESSIAN")) {
01891     data->runtype = MOLFILE_RUNTYPE_HESSIAN;
01892   }
01893   else if (!strcmp(&word[1][0],"RUNTYP=SURFACE")) {
01894     data->runtype = MOLFILE_RUNTYPE_SURFACE;
01895   }
01896   else if (!strcmp(&word[1][0],"RUNTYP=GRADIENT")) {
01897     data->runtype = MOLFILE_RUNTYPE_GRADIENT;
01898   }
01899   else if (!strcmp(&word[1][0],"RUNTYP=MEX")) {
01900     data->runtype = MOLFILE_RUNTYPE_MEX;
01901   }
01902   else {
01903 #ifdef DEBUGGING
01904     printf("gamessplugin) The %s is currently not supported \n",
01905            &word[1][0]);
01906 #endif
01907     data->runtype = MOLFILE_RUNTYPE_UNKNOWN;
01908   }
01909   printf("gamessplugin) File generated via %s \n",&word[1][0]);
01910 
01911 
01912   /* determine SCFTYP */
01913   if (!strcmp(&word[0][0],"SCFTYP=RHF")) {
01914     data->scftype = MOLFILE_SCFTYPE_RHF;
01915   }
01916   else if (!strcmp(&word[0][0],"SCFTYP=UHF")) {
01917     data->scftype = MOLFILE_SCFTYPE_UHF;
01918   }
01919   else if (!strcmp(&word[0][0],"SCFTYP=ROHF")) {
01920     data->scftype = MOLFILE_SCFTYPE_ROHF;
01921   }
01922   else if (!strcmp(&word[0][0],"SCFTYP=GVB")) {
01923     data->scftype = MOLFILE_SCFTYPE_GVB;
01924   }
01925   else if (!strcmp(&word[0][0],"SCFTYP=MCSCF")) {
01926     data->scftype = MOLFILE_SCFTYPE_MCSCF;
01927   }
01928   else if (!strcmp(&word[0][0],"SCFTYP=NONE")) {
01929     data->scftype = MOLFILE_SCFTYPE_NONE;
01930   }
01931   else {
01932     /* if we don't find a supported SCFTYP we bomb out; this
01933      * might be a little drastic */
01934     printf("gamessplugin) %s is currently not supported \n",
01935            &word[0][0]);
01936     return FALSE;
01937   }
01938   printf("gamessplugin) Type of wavefunction used %s \n",
01939          &word[0][0]);
01940 
01941   /* scan for MPLEVL, CITYP, CCTYP, VBTYP info; */
01942   GET_LINE(buffer, data->file);
01943   sscanf(buffer,"%s %s %*s %s",&word[0][0],&word[1][0],&word[2][0]);
01944 
01945   if (!strcmp(&word[0][0],"MPLEVL=")) {
01946     /* Moller-Plesset perturbation level */
01947     printf("gamessplugin) MP perturbation level %s \n",&word[1][0]);
01948     data->mplevel = atoi(&word[1][0]);
01949 
01950     /* determine CITYP */
01951     if      (!strcmp(&word[2][0],"=NONE"))  data->citype = CI_NONE;
01952     else if (!strcmp(&word[2][0],"=CIS"))   data->citype = CI_CIS;
01953     else if (!strcmp(&word[2][0],"=ALDET")) data->citype = CI_ALDET;
01954     else if (!strcmp(&word[2][0],"=ORMAS")) data->citype = CI_ORMAS;
01955     else if (!strcmp(&word[2][0],"=GUGA"))  data->citype = CI_GUGA;
01956     else if (!strcmp(&word[2][0],"=FSOCI")) data->citype = CI_FSOCI;
01957     else if (!strcmp(&word[2][0],"=GENCI")) data->citype = CI_GENCI;
01958     else                                    data->citype = CI_UNKNOWN;
01959     printf("gamessplugin) CI method %s \n",&word[2][1]);
01960 
01961     GET_LINE(buffer, data->file);
01962     sscanf(buffer,"%s %s",&word[0][0],&word[1][0]);
01963   }
01964 
01965   /* scan for DFTTYP, TDDFT info; */
01966   if (!strncmp(&word[0][0],"DFTTYP=", 7)) {
01967     printf("gamessplugin) Density functional used is %s \n",&word[0][7]);
01968     GET_LINE(buffer, data->file);
01969   }
01970 
01971 
01972   /* find the coordinate type in next line */
01973   while ( (temp=strstr(buffer,"COORD =")) == NULL ) {
01974     GET_LINE(buffer, data->file);;
01975   }
01976   strncpy(data->geometry, trimright(temp+7), sizeof(data->geometry)); 
01977   printf("gamessplugin) Coordinate type used is %s \n", data->geometry);
01978 
01979   fseek(data->file, filepos, SEEK_SET);
01980   return TRUE;
01981 }
01982 
01983 
01984 /* Read input parameters regarding calculation of 
01985  * certain molecular properties such as electrostatic
01986  * moments and the MEP. */
01987 static int get_properties_input(qmdata_t *data) {
01988   /* TODO!! */
01989   return TRUE;
01990 }
01991 
01992 
01993 /* Read symmetry point group and highest axis.
01994  * Currently these values are not used yet, but if this
01995  * section was not found the file is corrupt. */
01996 static int get_symmetry(qmdata_t *data) {
01997   char buffer[BUFSIZ];
01998   char *sep, tmp[BUFSIZ];
01999 
02000   //will need to go back somewhat in future to
02001   //check for XMCQDPT calculation
02002   long filepos = ftell(data->file);
02003 
02004   if (goto_keyline(data->file, "THE POINT GROUP IS",
02005                     "1 ELECTRON INTEGRALS", NULL) != FOUND) {
02006     printf("gamessplugin) No symmetry info found!\n");
02007     return FALSE;
02008   }
02009 
02010   GET_LINE(buffer, data->file);
02011   sscanf(buffer," THE POINT GROUP IS %s", data->pointgroup);
02012 
02013   sep = strchr(data->pointgroup, ',');
02014   if (sep) *sep = '\0';
02015   trimright(data->pointgroup);
02016 
02017   sep = strstr(buffer, "NAXIS=") + 6;
02018   strncpy(tmp, sep, 2); tmp[2] = '\0';
02019   data->naxis = atoi(tmp);
02020 
02021   sep = strstr(buffer, "ORDER=") + 6;
02022   sscanf(sep, "%d", &data->order);
02023 
02024   printf("gamessplugin) Point group = %s, naxis = %d, order = %d\n",
02025          data->pointgroup, data->naxis, data->order);
02026 
02027   fseek(data->file, filepos, SEEK_SET);
02028 
02029   return TRUE;
02030 }
02031 
02032 
02033 /* Read MCSCF input data */
02034 static int get_mcscf(qmdata_t *data) {
02035   gmsdata *gms = (gmsdata *)data->format_specific_data;
02036   char buffer[BUFSIZ];
02037   char *temp;
02038   long filepos;
02039   int tmp;
02040 
02041   filepos = ftell(data->file);
02042 
02043   if (gms->have_pcgamess){
02044       if (pass_keyline(data->file,"XMCQDPT INPUT PARAMETERS",
02045                         "DONE SETTING UP THE RUN") != FOUND) {
02046          //STANDARD MCSCF in GAMESS
02047          if(pass_keyline(data->file, "MCSCF CALCULATION",
02048                        "ITER     TOTAL ENERGY") != FOUND)
02049                 return FALSE;
02050 
02051           if (goto_keyline(data->file, "-CORE-    -INTERNAL-  -EXTERNAL-",
02052                            "ITER     TOTAL ENERGY", NULL) != FOUND)
02053             return FALSE;
02054 
02055           while ( (temp=strstr(buffer,"NFZC=")) == NULL ) {
02056             GET_LINE(buffer, data->file);
02057           }
02058           strncpy(buffer, trimright(temp+6), 5); 
02059           sscanf(buffer, "%d", &data->mcscf_num_core);
02060 
02061           while ( (temp=strstr(buffer,"NMCC=")) == NULL ) {
02062             GET_LINE(buffer, data->file);
02063           }
02064           strncpy(buffer, trimright(temp+6), 5); 
02065           sscanf(buffer, "%d", &tmp);
02066           data-> mcscf_num_core += tmp;
02067           printf("gamessplugin) Number of MCSCF core orbitals = %d\n",
02068              data->mcscf_num_core);
02069       }
02070       else{
02071           //XMCQDPT
02072           while ( (temp=strstr(buffer,"# OF FROZEN CORE ORBITALS")) == NULL ) {
02073             GET_LINE(buffer, data->file);
02074           }
02075           sscanf(buffer, "%*s %*s %*s %*s %*s %*s %d",&data->mcscf_num_core);
02076 
02077           GET_LINE(buffer,data->file);
02078           sscanf(buffer, "%*s %*s %*s %*s %*s %*s %d",&tmp);
02079           data->mcscf_num_core += tmp;
02080           printf("gamessplugin) Number of MCSCF core orbitals = %d\n",
02081              data->mcscf_num_core);
02082           printf("gamessplugin) XMCQDPT2 not supported.\n");
02083           //set scftype to none since XMCQDPT2 not supported
02084           data->scftype = MOLFILE_SCFTYPE_NONE;
02085 
02086       } 
02087   } 
02088   else {
02089       if (pass_keyline(data->file, "MCSCF CALCULATION",
02090                        "ITER     TOTAL ENERGY") != FOUND)
02091         return FALSE;
02092 
02093       if (goto_keyline(data->file, "NUMBER OF CORE ORBITALS",
02094                        "ITER     TOTAL ENERGY", NULL) != FOUND)
02095         return FALSE;
02096 
02097       GET_LINE(buffer, data->file);
02098       sscanf(buffer," NUMBER OF CORE ORBITALS          = %d",
02099              &data->mcscf_num_core);
02100 
02101       printf("gamessplugin) Number of MCSCF core orbitals = %d\n",
02102          data->mcscf_num_core);
02103   }
02104 
02105   fseek(data->file, filepos, SEEK_SET);
02106   return TRUE;
02107 }
02108 
02109 
02110 /* Read the first trajectory frame. */
02111 static int read_first_frame(qmdata_t *data) {
02112   /* The angular momentum is populated in get_wavefunction 
02113    * which is called by get_traj_frame(). We have obtained
02114    * the array size wavef_size already from the basis set
02115    * statistics */
02116   data->angular_momentum = (int*)calloc(3*data->wavef_size, sizeof(int));
02117 
02118   /* Try reading the first frame. 
02119    * If there is only one frame then also read the
02120    * final wavefunction. */
02121   if (!get_traj_frame(data, data->atoms, data->numatoms)) {
02122     return FALSE;
02123   }
02124 
02125   return TRUE;
02126 }
02127 
02128 /******************************************************
02129  *
02130  * Reads the info printed after the geometry search
02131  * has finished or whatever analysis was done in a 
02132  * single point run, e.g. ESP charges, Hessian, etc.
02133  * Rewinds to the beginning of the search when done,
02134  * because we read this part at in the initial phase
02135  * and might have to look for additional timesteps
02136  * later.
02137  *
02138  ******************************************************/
02139 static int get_final_properties(qmdata_t *data) {
02140   qm_timestep_t *ts;
02141   long filepos;
02142   filepos = ftell(data->file);
02143   ts = data->qm_timestep + data->num_frames-1;
02144 
02145   /* Go to end of trajectory */
02146   fseek(data->file, data->end_of_traj, SEEK_SET);
02147 
02148   printf("gamessplugin) Reading final properties section (timestep %d):\n",
02149          data->num_frames-1);
02150   printf("gamessplugin) ===============================================\n");
02151 
02152   /* Read population analysis (Mulliken and Lowdin charges),
02153    * but only if wasn't read already (while parsing the first
02154    * timestep). */
02155   if (!ts->have_mulliken && get_population(data, ts)) {
02156     printf("gamessplugin) Mulliken charges found\n");
02157   }
02158 
02159   if (get_esp_charges(data)) {
02160     printf("gamessplugin) ESP charges found\n");
02161   } 
02162 
02163   if (data->runtype == MOLFILE_RUNTYPE_GRADIENT ||
02164       data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02165     if (get_final_gradient(data, ts)) {
02166       printf("gamessplugin) Final gradient found\n");
02167     }
02168   }
02169 
02170 
02171   if (data->runtype == MOLFILE_RUNTYPE_HESSIAN) {
02172     /* try reading the hessian matrix in internal and
02173      * cartesian coordinates as well as the internal
02174      * coordinates together with their associated
02175      * force constants */
02176     
02177     if (!get_int_hessian(data)) {
02178       printf("gamessplugin) No internal Hessian matrix found.\n");
02179     }
02180     if (!get_cart_hessian(data)) {
02181       printf("gamessplugin) \n");
02182       printf("gamessplugin) Could not determine the cartesian \n");
02183       printf("gamessplugin) Hessian matrix!! \n");
02184       printf("gamessplugin) \n");
02185     }
02186 
02187     /* read the wavenumbers, intensities of the normal modes 
02188      * as well as the modes themselves */
02189     if (!get_normal_modes(data)) {
02190       printf("gamessplugin) \n");
02191       printf("gamessplugin) No normal modes found.\n");
02192       printf("gamessplugin) \n");
02193     }
02194   }
02195 
02196   /* Read localized orbitals if there are any */
02197   read_localized_orbitals(data);
02198 
02199 
02200   fseek(data->file, filepos, SEEK_SET);
02201   return TRUE; 
02202 }
02203 
02204 
02205 /* Read localized orbitals (Boys/Ruedenberg/Pipek) */
02206 static int read_localized_orbitals(qmdata_t *data) {
02207   int i;
02208   qm_timestep_t *ts;
02209   qm_wavefunction_t *wavef;
02210 
02211   /* Move past the listing of the canonical MOs */
02212   pass_keyline(data->file, "ENERGY COMPONENTS", NULL);
02213 
02214   ts = data->qm_timestep + data->num_frames-1;
02215 
02216   for (i=0; i<2; i++) {
02217     wavef = add_wavefunction(ts);
02218 
02219     if (get_wavefunction(data, ts, wavef) == FALSE ||
02220         (wavef->type!=MOLFILE_WAVE_BOYS &&
02221          wavef->type!=MOLFILE_WAVE_PIPEK &&
02222          wavef->type!=MOLFILE_WAVE_RUEDEN)) {
02223       del_wavefunction(ts);
02224       return FALSE;
02225     }
02226     else {
02227       char typestr[64];
02228       if (wavef->spin==SPIN_ALPHA) {
02229         strcpy(typestr, "alpha");
02230       }
02231       else if (wavef->spin==SPIN_BETA) {
02232         strcpy(typestr, "beta");
02233       }
02234       wavef->mult = data->multiplicity;
02235       wavef->energy = ts->scfenergies[ts->num_scfiter-1];
02236 
02237       printf("gamessplugin) Localized orbitals (%s) found for timestep %d\n",
02238              typestr, data->num_frames-1);
02239     }
02240   }
02241 
02242   return TRUE;
02243 }
02244 
02245 
02246 
02247 /********************************************************
02248  *
02249  * Read basis set and orbital statistics such as
02250  * # of shells, # of A/B orbitals, # of electrons, 
02251  * multiplicity and total charge
02252  *
02253  ********************************************************/
02254 static int get_basis_stats(qmdata_t *data) {
02255 
02256   gmsdata *gms = (gmsdata *)data->format_specific_data;
02257 
02258   char buffer[BUFSIZ];
02259   char word[7][BUFSIZ];
02260   int i;
02261 
02262   buffer[0] = '\0';
02263   for (i=0; i<7; i++) word[i][0] = '\0';
02264 
02265   /* look for the orbital/charge/... info section */
02266   if(gms->have_pcgamess){
02267     if (!pass_keyline(data->file, "TOTAL NUMBER OF SHELLS", NULL)){
02268         printf("ERROR!\n");
02269         return FALSE;
02270      }
02271       /* # cartesian gaussian function = wavefunction size */
02272       GET_LINE(buffer, data->file);
02273       sscanf(buffer,"%*s %*s %*s %*s %*s %*s %d",
02274              &(data->wavef_size));
02275 
02276       /* number of electrons */
02277       GET_LINE(buffer, data->file);
02278       sscanf(buffer,"%*s %*s %*s %*s %d",
02279              &(data->num_electrons));
02280 
02281       /* charge of the molecule */
02282       GET_LINE(buffer, data->file);
02283       sscanf(buffer,"%*s %*s %*s %*s %d",
02284              &(data->totalcharge));
02285 
02286       /* Multiplicity of the molecule.
02287        * Multiplicity is actually defined per wavefunction
02288        * but in some cases where there's only one wavefunction
02289        * in the output or they all have the same multiplicity
02290        * it will not be printed with the wavefunction.
02291        * Thus we use this one as default value. */
02292       GET_LINE(buffer, data->file);
02293       sscanf(buffer,"%*s %*s %*s %d",
02294              &(data->multiplicity));
02295 
02296       /* number of A orbitals */
02297       /* Note the different number of items per line for A/B orbitals
02298        * due to "(ALPHA)" and "(BETA )" !! */
02299       GET_LINE(buffer, data->file);
02300       sscanf(buffer,"%*s %*s %*s %*s %*s %*s %d",
02301              &(data->num_occupied_A)); 
02302         
02303       /* number of B orbitals */
02304       GET_LINE(buffer, data->file);
02305       sscanf(buffer,"%*s %*s %*s %*s %*s %*s %*s %d",
02306              &(data->num_occupied_B)); 
02307 
02308   }
02309   else {
02310     if (!pass_keyline(data->file, "TOTAL NUMBER OF BASIS", NULL))
02311         return FALSE;
02312 
02313       /* # cartesian gaussian function = wavefunction size */
02314       GET_LINE(buffer, data->file);
02315       sscanf(buffer,"%*s %*s %*s %*s %*s %*s %*s %d",
02316              &(data->wavef_size));
02317 
02318       /* number of electrons */
02319       GET_LINE(buffer, data->file);
02320       sscanf(buffer,"%*s %*s %*s %*s %d",
02321              &(data->num_electrons));
02322 
02323       /* charge of the molecule */
02324       GET_LINE(buffer, data->file);
02325       sscanf(buffer,"%*s %*s %*s %*s %d",
02326              &(data->totalcharge));
02327 
02328       /* Multiplicity of the molecule.
02329        * Multiplicity is actually defined per wavefunction
02330        * but in some cases where there's only one wavefunction
02331        * in the output or they all have the same multiplicity
02332        * it will not be printed with the wavefunction.
02333        * Thus we use this one as default value. */
02334       GET_LINE(buffer, data->file);
02335       sscanf(buffer,"%*s %*s %*s %d",
02336              &(data->multiplicity));
02337 
02338       /* number of A orbitals */
02339       /* Note the different number of items per line for A/B orbitals
02340        * due to "(ALPHA)" and "(BETA )" !! */
02341       GET_LINE(buffer, data->file);
02342       sscanf(buffer,"%*s %*s %*s %*s %*s %*s %d",
02343              &(data->num_occupied_A)); 
02344         
02345       /* number of B orbitals */
02346       GET_LINE(buffer, data->file);
02347       sscanf(buffer,"%*s %*s %*s %*s %*s %*s %*s %d",
02348              &(data->num_occupied_B)); 
02349   }
02350 
02351 
02352   printf("gamessplugin) Number of Electrons: %d \n",
02353       data->num_electrons);
02354 
02355   printf("gamessplugin) Charge of Molecule : %d \n",
02356       data->totalcharge);
02357 
02358   printf("gamessplugin) Multiplicity of Wavefunction: %d \n",
02359       data->multiplicity);
02360 
02361   printf("gamessplugin) Number of occupied A / B orbitals: %d / %d \n",\
02362       data->num_occupied_A, data->num_occupied_B);
02363 
02364   printf("gamessplugin) Number of gaussian basis functions: %d \n",\
02365       data->wavef_size);
02366 
02367  
02368   return TRUE;
02369 }
02370 
02371 
02372 
02373 /*******************************************************
02374  *
02375  * Reads in the $GUESS options.
02376  *
02377  * ******************************************************/
02378 static int get_guess_options(qmdata_t *data)
02379 {
02380   char word[BUFSIZ];
02381   char buffer[BUFSIZ];
02382   long filepos;
02383   filepos = ftell(data->file);
02384 
02385   /* initialize buffers */
02386   buffer[0] = '\0';
02387   word[0]   = '\0';
02388 
02389   /* parse for GUESS field */
02390   if (pass_keyline(data->file, "GUESS OPTIONS",
02391                    "2 ELECTRON INTEGRALS") != FOUND) {
02392     printf("gamessplugin) No GUESS OPTIONS found.\n");
02393     fseek(data->file, filepos, SEEK_SET);
02394 
02395     /* This section id not mandatory, there are a few
02396        calculation types the don't print it, so we
02397        always return TRUE*/
02398     return TRUE;
02399   }
02400 
02401   /* next line contains all we need */
02402   eatline(data->file, 1);
02403   GET_LINE(buffer, data->file);
02404   sscanf(buffer," GUESS %s NORB",&word[0]);
02405 
02406   /* the first character is '=', we skip it */
02407   strncpy(data->guess,&word[1], sizeof(data->guess));
02408 
02409   printf("gamessplugin) Run was performed with GUESS = %s \n",
02410           data->guess);
02411 
02412  /* Since this block occurs in the middle of first frame
02413   * we need to rewind. */
02414   fseek(data->file, filepos, SEEK_SET);
02415 
02416   return TRUE;
02417 }
02418 
02419 
02420 
02421 /*******************************************************
02422  *
02423  * Read the basis set data into hierarchical structures.
02424  * L-shells are expanded into an S-shell and a P-shell.
02425  *
02426  * ******************************************************/
02427 int get_basis(qmdata_t *data) {
02428 
02429   gmsdata *gms = (gmsdata *)data->format_specific_data;
02430 
02431   char buffer[BUFSIZ];
02432   char word[4][BUFSIZ];
02433   int i = 0; 
02434   int success = 0;
02435   int numread, numshells;
02436   shell_t *shell;
02437   long filepos;
02438 
02439   if (!strcmp(data->gbasis, "MNDO") ||
02440       !strcmp(data->gbasis, "AM1")  ||
02441       !strcmp(data->gbasis, "PM3")) {
02442     /* Semiempirical methods are based on STOs.
02443      * The only parameter we need for orbital rendering
02444      * are the exponents zeta for S, P, D,... shells for
02445      * each atom. Since GAMESS doesn't print these values
02446      * we skip reading the basis set and hardcode the
02447      * parameters in tables in VMD. */
02448     return TRUE;
02449   }
02450 
02451   /* Search for "ATOMIC BASIS SET" line */
02452   if (pass_keyline(data->file, "ATOMIC BASIS SET", 
02453                     "$CONTRL OPTIONS") !=FOUND ) {
02454     printf("gamessplugin) No basis set found!\n");
02455 
02456     return FALSE;
02457   }
02458 
02459   /* initialize buffers */
02460   buffer[0] = '\0';
02461   for (i=0; i<3; i++) word[i][0] = '\0';
02462   
02463 
02464   /* skip the next 5 lines */
02465   eatline(data->file, 5);
02466 
02467   /* Allocate space for the basis for all atoms */
02468   /* When the molecule is symmetric the actual number atoms with
02469    * a basis set could be smaller */
02470   data->basis_set = (basis_atom_t*)calloc(data->numatoms, sizeof(basis_atom_t));
02471 
02472 
02473   i = 0; /* basis atom counter */
02474 
02475   do {
02476     prim_t *prim = NULL;
02477     char shelltype;
02478     int numprim = 0;
02479     int icoeff = 0;
02480     filepos = ftell(data->file);
02481     GET_LINE(buffer, data->file);
02482       
02483     /* Count the number of relevant words in the line. */
02484     numread = sscanf(buffer,"%s %s %s %s",&word[0][0], &word[1][0],
02485            &word[2][0], &word[3][0]);
02486 
02487     switch (numread) {
02488       case 1:
02489         /* Next atom */
02490         strcpy(data->basis_set[i].name, &word[0][0]);
02491 
02492         /* skip initial blank line */
02493         eatline(data->file, 1);
02494 
02495         /* read the basis set for the current atom */
02496         shell = (shell_t*)calloc(1, sizeof(shell_t)); 
02497         numshells = 0;
02498 
02499         do {
02500           filepos = ftell(data->file);
02501           numprim = read_shell_primitives(data, &prim, &shelltype, icoeff, gms->have_pcgamess);
02502 
02503           if (numprim>0) {
02504             /* make sure we have eiter S, L, P, D, F or G shells */
02505             if ( (shelltype!='S' && shelltype!='L' && shelltype!='P' && 
02506                   shelltype!='D' && shelltype!='F' && shelltype!='G') ) {
02507               printf("gamessplugin) WARNING ... %c shells are not supported \n", shelltype);
02508             }
02509             
02510             /* create new shell */
02511             if (numshells) {
02512               shell = (shell_t*)realloc(shell, (numshells+1)*sizeof(shell_t));
02513             }
02514             shell[numshells].numprims = numprim;
02515             /* assign a numeric shell type */
02516             shell[numshells].type = shelltype_int(shelltype);
02517             shell[numshells].prim = prim;
02518             data->num_basis_funcs += numprim;
02519 
02520             /* We split L-shells into one S and one P-shell.
02521              * I.e. for L-shells we have to go back read the shell again
02522              * this time using the second contraction coefficients. */
02523             if (shelltype=='L' && !icoeff) {
02524               fseek(data->file, filepos, SEEK_SET);
02525               icoeff++;
02526             } else if (shelltype=='L' && icoeff) {
02527               shell[numshells].type = SP_P_SHELL;
02528               icoeff = 0;  /* reset the counter */
02529             }
02530 
02531             numshells++;
02532           }
02533         } while (numprim);
02534 
02535         /* store shells in atom */
02536         data->basis_set[i].numshells = numshells;
02537         data->basis_set[i].shell = shell;
02538 
02539         /* Update total number of basis functions */
02540         data->num_shells += numshells;
02541         i++;
02542 
02543         /* go back one line so that we can read the name of the
02544          * next atom */
02545         fseek(data->file, filepos, SEEK_SET);
02546 
02547         break;
02548 
02549       case 4:
02550         /* this is the very end of the basis set */
02551         if(gms->have_pcgamess){
02552             if (!strcmp(&word[0][0],"TOTAL")  &&
02553                 !strcmp(&word[1][0],"NUMBER") && 
02554                 !strcmp(&word[2][0],"OF")     &&
02555                 !strcmp(&word[3][0],"SHELLS")) {
02556               success = 1;
02557               /* go back one line so that get_basis_stats()
02558                  can use this line as a keystring. */
02559               fseek(data->file, filepos, SEEK_SET);
02560             }
02561         }
02562         else {
02563             if (!strcmp(&word[0][0],"TOTAL")  &&
02564                 !strcmp(&word[1][0],"NUMBER") && 
02565                 !strcmp(&word[2][0],"OF")     &&
02566                 !strcmp(&word[3][0],"BASIS")) {
02567               success = 1;
02568               /* go back one line so that get_basis_stats()
02569                  can use this line as a keystring. */
02570               fseek(data->file, filepos, SEEK_SET);
02571             }
02572         }
02573         break;
02574     }
02575 
02576   } while (!success);
02577 
02578 
02579   printf("gamessplugin) Parsed %d uncontracted basis functions for %d atoms.\n",
02580          data->num_basis_funcs, i);
02581 
02582   data->num_basis_atoms = i;
02583 
02584 
02585   /* allocate and populate flat arrays needed for molfileplugin */
02586   return fill_basis_arrays(data);
02587 }
02588 
02589 
02590 /**************************************************
02591  *
02592  * Convert shell type from char to int.
02593  *
02594  ************************************************ */
02595 static int shelltype_int(char type) {
02596   int shelltype;
02597 
02598   switch (type) {
02599     case 'L':
02600       /* SP_P shells are assigned in get_basis() */
02601       shelltype = SP_S_SHELL;
02602       break;
02603     case 'S':
02604       shelltype = S_SHELL;
02605       break;
02606     case 'P':
02607       shelltype = P_SHELL;
02608       break;
02609     case 'D':
02610       shelltype = D_SHELL;
02611       break;
02612     case 'F':
02613       shelltype = F_SHELL;
02614       break;
02615     case 'G':
02616       shelltype = G_SHELL;
02617       break;
02618     default:
02619       shelltype = UNK_SHELL;
02620       break;
02621   }
02622 
02623   return shelltype;
02624 }
02625 
02626 
02627 
02628 /******************************************************
02629  *
02630  * Populate the flat arrays containing the basis
02631  * set data.
02632  *
02633  ******************************************************/
02634 static int fill_basis_arrays(qmdata_t *data) {
02635   gmsdata *gms = (gmsdata *)data->format_specific_data;
02636   int i, j, k;
02637   int shellcount = 0;
02638   int primcount = 0;
02639 
02640   float *basis;
02641   int *num_shells_per_atom;
02642   int *num_prim_per_shell;
02643   int *shell_types;
02644   int *atomicnum_per_basisatom;
02645 
02646   /* Count the total number of primitives which
02647    * determines the size of the basis array. */
02648   for(i=0; i<data->num_basis_atoms; i++) {
02649     for (j=0; j<data->basis_set[i].numshells; j++) {
02650       primcount += data->basis_set[i].shell[j].numprims;
02651     }
02652   }
02653 
02654   /* reserve space for pointer to array containing basis
02655    * info, i.e. contraction coeficients and expansion 
02656    * coefficients; need 2 entries per basis function, i.e.
02657    * exponent and contraction coefficient; also,
02658    * allocate space for the array holding the orbital symmetry
02659    * information per primitive Gaussian.
02660    * Finally, initialize the arrays holding the number of 
02661    * shells per atom and the number of primitives per shell*/
02662   basis = (float *)calloc(2*primcount,sizeof(float));
02663 
02664   /* make sure memory was allocated properly */
02665   if (basis == NULL) {
02666     PRINTERR;
02667     return FALSE;
02668   }
02669 
02670   shell_types = (int *)calloc(data->num_shells, sizeof(int));
02671   
02672   /* make sure memory was allocated properly */
02673   if (shell_types == NULL) {
02674     PRINTERR; 
02675     return FALSE;
02676   }
02677 
02678   num_shells_per_atom = (int *)calloc(data->num_basis_atoms, sizeof(int));
02679 
02680   /* make sure memory was allocated properly */
02681   if (num_shells_per_atom == NULL) {
02682     PRINTERR; 
02683     return FALSE;
02684   }
02685 
02686   num_prim_per_shell = (int *)calloc(data->num_shells, sizeof(int));
02687 
02688   /* make sure memory was allocated properly */
02689   if (num_prim_per_shell == NULL) {
02690     PRINTERR;
02691     return FALSE;
02692   }
02693 
02694   atomicnum_per_basisatom = (int *)calloc(data->num_basis_atoms, sizeof(int));
02695 
02696   /* make sure memory was allocated properly */
02697   if (atomicnum_per_basisatom == NULL) {
02698     PRINTERR;
02699     return FALSE;
02700   }
02701 
02702 
02703   /* store pointers in struct qmdata_t */
02704   data->basis = basis;
02705   data->shell_types = shell_types;
02706   data->num_shells_per_atom = num_shells_per_atom;
02707   data->num_prim_per_shell = num_prim_per_shell;
02708   data->atomicnum_per_basisatom = atomicnum_per_basisatom;
02709 
02710   /* Go through all basis set atoms and try to assign the
02711    * atomic numbers. The basis set atoms are specified by 
02712    * name strings (the same as in the coordinate section,
02713    * except for FMO calcs.) and we try to match the names
02714    * from the two lists. The basis set atom list is symmetry
02715    * unique while the coordinate atom list is complete.*/
02716   primcount = 0;
02717   for (i=0; i<data->num_basis_atoms; i++) {
02718     int found = 0;
02719 
02720     /* For this basis atom find a matching atom from the
02721      * coordinate atom list. */
02722     for(j=0; j<data->numatoms; j++) {
02723       char basisname[BUFSIZ];
02724       strcpy(basisname, data->basis_set[i].name);
02725 
02726       /* for FMO calculations we have to strip the "-n" tail
02727        * of the basis atom name. */
02728       if (gms->have_fmo) {
02729         *strchr(basisname, '-') = '\0';
02730       }
02731 
02732       if (!strcmp(data->atoms[j].type, basisname)) {
02733         found = 1;
02734         break;
02735       }
02736     }
02737     if (!found) {
02738       printf("gamessplugin) WARNING: Couldn't find atomic number for basis set atom %s\n",
02739              data->basis_set[i].name);
02740       data->basis_set[i].atomicnum = 0;
02741       atomicnum_per_basisatom[i] = 0;
02742     } else {
02743       /* assign atomic number */
02744       data->basis_set[i].atomicnum = data->atoms[j].atomicnum;
02745       atomicnum_per_basisatom[i]   = data->atoms[j].atomicnum;
02746     }
02747     num_shells_per_atom[i] = data->basis_set[i].numshells;
02748 
02749     for (j=0; j<data->basis_set[i].numshells; j++) {
02750       shell_types[shellcount] = data->basis_set[i].shell[j].type;
02751       num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
02752 
02753       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
02754         basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
02755         basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
02756         primcount++;
02757       }
02758       shellcount++;
02759     }
02760   } 
02761 
02762   return TRUE;
02763 }
02764 
02765 
02766 /******************************************************
02767  *
02768  * read all primitives for the current shell
02769  *
02770  ******************************************************/
02771 static int read_shell_primitives(qmdata_t *data, prim_t **prim, char *shelltype,
02772                                  int icoeff, int pcgamess) {
02773   char buffer[BUFSIZ];
02774   float exponent = 0.0; 
02775   float contract[2] = {0.0, 0.0};
02776   int shell, success;
02777   int primcounter = 0;
02778   (*prim) = (prim_t*)calloc(1, sizeof(prim_t));
02779 
02780   do {
02781     GET_LINE(buffer, data->file);
02782       if (pcgamess)
02783         success = sscanf(buffer,"%d %c %*s %f %f %*s %*s %f", &shell,
02784                        shelltype,
02785                        &exponent, &contract[0], &contract[1]); 
02786 
02787       else
02788         success = sscanf(buffer,"%d %c %*s %f %f %f", &shell,
02789                        shelltype,
02790                        &exponent, &contract[0], &contract[1]); 
02791 
02792     /* store in basis array and increase the counter */ 
02793     switch (success) {
02794       case 4:
02795         if (primcounter) {
02796           *prim = (prim_t*)realloc(*prim, (primcounter+1)*sizeof(prim_t));
02797         }
02798 
02799         /* store exponent */
02800         (*prim)[primcounter].exponent = exponent;
02801           
02802         /* store coefficient */
02803         (*prim)[primcounter].contraction_coeff = contract[0];
02804         
02805         primcounter++;
02806         break;
02807 
02808       case 5:
02809         if (primcounter) {
02810           *prim = (prim_t*)realloc(*prim, (primcounter+1)*sizeof(prim_t));
02811         }
02812 
02813         /* store exponent */
02814         (*prim)[primcounter].exponent = exponent;
02815           
02816         /* store coefficient */
02817         (*prim)[primcounter].contraction_coeff = contract[icoeff];
02818         
02819         primcounter++;
02820         break;
02821 
02822       case -1:
02823         /* otherwise it's an empty line which represents the end of the shell */
02824         break;
02825 
02826       case 1:
02827         /* the user had given the next atom a numeric name */
02828         success = -1;
02829         break;
02830     }
02831 
02832   } while(success>0);
02833 
02834   if (!primcounter) free(*prim);
02835 
02836   return primcounter;
02837 }
02838 
02839 
02840 
02841 /******************************************************
02842  *
02843  * this function extracts the trajectory information
02844  * from the output file
02845  *
02846  * *****************************************************/
02847 static int get_traj_frame(qmdata_t *data, qm_atom_t *atoms,
02848                           int natoms) {
02849   gmsdata *gms = (gmsdata *)data->format_specific_data;
02850   qm_timestep_t *cur_ts;
02851   char buffer[BUFSIZ];
02852   char word[BUFSIZ];
02853   int units;
02854   buffer[0] = '\0';
02855   word[0]   = '\0';
02856 
02857   printf("gamessplugin) Timestep %d:\n", data->num_frames_read);
02858   printf("gamessplugin) ============\n");
02859 
02860   fseek(data->file, data->filepos_array[data->num_frames_read], SEEK_SET);
02861 
02862   /* Read the coordinate block */
02863   if (data->runtype==MOLFILE_RUNTYPE_OPTIMIZE ||
02864       data->runtype==MOLFILE_RUNTYPE_SADPOINT) {
02865     goto_keyline(data->file, "COORDINATES OF ALL ATOMS", NULL);
02866     /* get the units */
02867     GET_LINE(buffer, data->file);
02868     sscanf(buffer, " COORDINATES OF ALL ATOMS ARE %s", word);
02869     units = !strcmp(word, "(BOHR)");
02870     eatline(data->file, 2);
02871 
02872     if (!get_coordinates(data->file, &data->atoms, units, &natoms)) {
02873       printf("gamessplugin) Couldn't find coordinates for timestep %d\n", data->num_frames_read);
02874     }
02875   }
02876   else if (data->runtype==MOLFILE_RUNTYPE_SURFACE) {
02877     if (pass_keyline(data->file, "HAS ENERGY VALUE",
02878                      "...... END OF ONE-ELECTRON INTEGRALS ......")
02879         == FOUND) {
02880       /* Read the coordinate block following 
02881        * ---- SURFACE MAPPING GEOMETRY ---- */
02882       int i, n;
02883       for (i=0; i<natoms; i++) {
02884         char atname[BUFSIZ];
02885         float x,y,z;
02886         GET_LINE(buffer, data->file);
02887         n = sscanf(buffer,"%s %f %f %f", atname, &x,&y,&z);
02888         if (n!=4 || strcmp(atname, data->atoms[i].type)) break;
02889         data->atoms[i].x = x;
02890         data->atoms[i].y = y;
02891         data->atoms[i].z = z;
02892       }
02893       if (i!=natoms) {
02894         printf("gamessplugin) Couldn't read surface mapping geometry for timestep %d\n", data->num_frames_read);
02895       }
02896     }
02897     else {
02898       /* Read the coordinate block following 
02899        * ATOM      ATOMIC                      COORDINATES (BOHR) */
02900       goto_keyline(data->file, "ATOM      ATOMIC", NULL);
02901       /* get the units */
02902       GET_LINE(buffer, data->file);
02903       sscanf(buffer, " ATOM      ATOMIC                      COORDINATES %s", word);
02904       units = !strcmp(word, "(BOHR)");
02905       eatline(data->file, 1);
02906       
02907       if (!get_coordinates(data->file, &data->atoms, units, &natoms)) {
02908         printf("gamessplugin) Couldn't find coordinates for timestep %d\n", data->num_frames_read);
02909       }
02910     }
02911   }
02912   /* XXX could merge this with OPTIMIZE/SADPOINT */
02913   else if (data->runtype==MOLFILE_RUNTYPE_MEX) {
02914     int numuniqueatoms = natoms;
02915     goto_keyline(data->file, "COORDINATES OF SYMMETRY UNIQUE ATOMS", NULL);
02916     /* get the units */
02917     GET_LINE(buffer, data->file);
02918     sscanf(buffer, " COORDINATES OF SYMMETRY UNIQUE ATOMS ARE %s", word);
02919     units = !strcmp(word, "(BOHR)");
02920     eatline(data->file, 2);
02921     if (!get_coordinates(data->file, &data->atoms, units, &numuniqueatoms)) {
02922       printf("gamessplugin) Expanding symmetry unique coordinates for timestep %d\n", data->num_frames_read);
02923 
02924       /* Create images of symmetry unique atoms so that we have
02925        * the full coordinate set. */
02926       symmetry_expand(&data->atoms, numuniqueatoms, natoms,
02927                       data->pointgroup, data->naxis);
02928     }
02929   }
02930 
02931   /* For FMO calculations we read the coordinates only
02932    * because the wavefunctions are printed per fragment
02933    * and VMD requires that there's a wavefunction present
02934    * for each atom.
02935    * A possible workaround would be to pad the wavefunctions
02936    * accordingly and add a wavefunction for each fragment. */
02937   if (gms->have_fmo) {
02938     data->num_frames_read++;
02939     return TRUE;
02940   }
02941 
02942   /* get a convenient pointer to the current qm timestep */
02943   cur_ts = data->qm_timestep + data->num_frames_read;
02944 
02945   /* read the SCF energies */
02946   if (get_scfdata(data, cur_ts) == FALSE) {
02947     printf("gamessplugin) Couldn't find SCF iterations for timestep %d\n",
02948            data->num_frames_read);
02949   }
02950 
02951   /* Try reading canonical alpha/beta wavefunction */
02952   check_add_wavefunctions(data, cur_ts);
02953 
02954 
02955   /* Read population analysis (Mulliken and Lowdin charges)
02956    * only if wasn't read already while parsing the final
02957    * property section. Otherwise we would potentially 
02958    * overwrite the data with empty fields. */
02959   if (!cur_ts->have_mulliken &&
02960       get_population(data, cur_ts)) {
02961     printf("gamessplugin) Mulliken/Loewdin charges found\n");
02962   }
02963 
02964   if (data->citype==CI_GUGA) {
02965     if (pass_keyline(data->file, "CI DENSITY MATRIX AND NATURAL ORBITALS",
02966                        "GRADIENT (HARTREE/BOHR)")) {
02967       int i, numstates=0, state;
02968       qm_wavefunction_t *wave_ci;
02969       goto_keyline(data->file, "NUMBER OF STATES", NULL);
02970       GET_LINE(buffer, data->file);
02971       trimleft(buffer);
02972       sscanf(buffer, " NUMBER OF STATES = %d", &numstates);
02973       printf("gamessplugin) Number of CI states = %d\n", numstates);
02974 
02975       for (i=0; i<numstates; i++) {
02976         float cienergy = 0.f;
02977         goto_keyline(data->file, "CI EIGENSTATE", NULL);
02978         GET_LINE(buffer, data->file);
02979         sscanf(buffer," CI EIGENSTATE %d %*s %*s %*s %f", &state, &cienergy);
02980         printf("gamessplugin) CI energy[%d] = %f\n", state-1, cienergy);
02981 
02982         wave_ci = add_wavefunction(cur_ts);
02983 
02984         if (get_wavefunction(data, cur_ts, wave_ci) == FALSE) {
02985           del_wavefunction(cur_ts);
02986           break;
02987         }
02988         else {
02989           /* canon =-1;*/
02990           wave_ci->exci = state-1;
02991           wave_ci->energy = cienergy;
02992           wave_ci->mult = data->multiplicity;
02993 
02994           printf("gamessplugin) Found %d CI natural orbitals for excited state %d, mult=%d\n",
02995                  wave_ci->num_orbitals, state-1, wave_ci->mult);
02996         }
02997       }
02998     }
02999   }
03000   else if (data->citype==CI_CIS) {
03001     if (pass_keyline(data->file,
03002                      "USING DAVIDSON ALGORITHM TO FIND CIS EIGENVALUES",
03003                      NULL)) {
03004       int i, numstates=0, state;
03005       qm_wavefunction_t *wave_ci;
03006       float *state_energies, *state_spinquant;
03007       goto_keyline(data->file, "NUMBER OF STATES REQUESTED", NULL);
03008       GET_LINE(buffer, data->file);
03009       trimleft(buffer);
03010       sscanf(buffer, " NUMBER OF STATES REQUESTED = %d", &numstates);
03011       printf("gamessplugin) Number of CIS states = %d\n", numstates);
03012 
03013       /* For CIS only the wavefunction for the excited state 
03014        * (specified by IROOT in $CIS group) will be printed.
03015        * Here we read in the energies for all states, and store 
03016        * the energy for the selected state in its wavefunction */
03017       state_energies  = calloc(numstates+1, sizeof(float));
03018       state_spinquant = calloc(numstates+1, sizeof(float));
03019       goto_keyline(data->file, "RHF REFERENCE ENERGY  =", NULL);
03020       GET_LINE(buffer, data->file);
03021       trimleft(buffer);
03022       sscanf(buffer, " RHF REFERENCE ENERGY  = %f", &state_energies[0]);
03023       state_spinquant[0] = 1.f;
03024 
03025       for (i=1; i<=numstates; i++) {
03026         goto_keyline(data->file, "EXCITED STATE", NULL);
03027         GET_LINE(buffer, data->file);
03028         trimleft(buffer);
03029         sscanf(buffer, " EXCITED STATE %*d  ENERGY= %f  S = %f",
03030                &state_energies[i], &state_spinquant[i]);
03031       }
03032 
03033       goto_keyline(data->file,
03034                    "CIS NATURAL ORBITAL OCCUPATION NUMBERS FOR EXCITED STATE",
03035                    NULL);
03036       GET_LINE(buffer, data->file);
03037       trimleft(buffer);
03038       sscanf(buffer,
03039              " CIS NATURAL ORBITAL OCCUPATION NUMBERS FOR EXCITED STATE %d",
03040              &state);
03041       
03042       wave_ci = add_wavefunction(cur_ts);
03043 
03044       if (get_wavefunction(data, cur_ts, wave_ci) == FALSE) {
03045         del_wavefunction(cur_ts);
03046       }
03047       else {
03048         wave_ci->exci = state;
03049         wave_ci->energy = state_energies[state];
03050         wave_ci->mult = 2*(int)state_spinquant[state]+1;
03051         printf("gamessplugin) Found %d CIS natural orbitals for excited state %d\n",
03052                wave_ci->num_orbitals, state);
03053       }
03054 
03055       free(state_energies);
03056       free(state_spinquant);
03057     }
03058   }
03059 
03060 
03061   /* Read the energy gradients (=forces on atoms) */
03062   if (get_gradient(data, cur_ts)) {
03063     printf("gamessplugin) Energy gradient found.\n");
03064   }
03065 
03066 
03067   /* If this is the last frame of the trajectory and the file
03068    * wasn't truncated and the program didn't terminate
03069    * abnormally then read the final wavefunction. */
03070   if ((data->runtype == MOLFILE_RUNTYPE_OPTIMIZE ||
03071        data->runtype == MOLFILE_RUNTYPE_SADPOINT) &&
03072       (data->num_frames_read+1 == data->num_frames &&
03073        (data->status == MOLFILE_QMSTATUS_UNKNOWN || 
03074         data->status == MOLFILE_QMSTATUS_OPT_CONV ||
03075         data->status == MOLFILE_QMSTATUS_OPT_NOT_CONV))) {
03076 
03077     /* We need to jump over the end of the trajectory because 
03078      * this is also the keystring for get_wavefunction() to
03079      * bail out. */
03080     if (data->status == MOLFILE_QMSTATUS_OPT_CONV || 
03081         data->status == MOLFILE_QMSTATUS_OPT_NOT_CONV) {
03082       fseek(data->file, data->end_of_traj, SEEK_SET);
03083     }
03084 
03085     /* Try to read final wavefunction and orbital energies
03086      * A preexisting canonical wavefunction for this timestep
03087      * with the same characteristics (spin, exci, info) will
03088      * be overwritten by the final wavefuntion if it has more
03089      * orbitals. */
03090     check_add_wavefunctions(data, cur_ts);
03091   }
03092 
03093 
03094   /* For MCSCF optimized orbitals no occupancies are given
03095    * but since their occupancies are identical to the ones
03096    * from natural orbitals we can use those. The natural 
03097    * orbitals are always listed right before the optimized
03098    * ones so we simply copy the data over. */
03099   if (cur_ts->numwave>=2 &&
03100       cur_ts->wave[cur_ts->numwave-1].type==MOLFILE_WAVE_MCSCFOPT &&
03101       cur_ts->wave[cur_ts->numwave-2].type==MOLFILE_WAVE_MCSCFNAT) {
03102     int i;
03103     qm_wavefunction_t *waveopt = &cur_ts->wave[cur_ts->numwave-1];
03104     qm_wavefunction_t *wavenat = &cur_ts->wave[cur_ts->numwave-2];
03105     waveopt->orb_occupancies = (float *)calloc(waveopt->num_orbitals,
03106                                                sizeof(float));
03107     /* Only the core and active natural orbitals are listed. 
03108      * We copy the occupancies for these orbitals and pad the
03109      * rest with zeros. */
03110     for (i=0; i<wavenat->num_orbitals; i++) {
03111       waveopt->orb_occupancies[i] = wavenat->orb_occupancies[i];
03112     }
03113     for (i=wavenat->num_orbitals; i<waveopt->num_orbitals; i++) {
03114       waveopt->orb_occupancies[i] = 0.f;
03115     }
03116     waveopt->has_occup = TRUE;
03117   }
03118 
03119 
03120   data->num_frames_read++;
03121 
03122   return TRUE;
03123 }
03124 
03125 
03126 /* Analyze the trajectory.
03127  * Read the parameters controlling geometry search and
03128  * find the end of the trajectory, couinting the frames
03129  * on the way. Store the filepointer for the beginning of
03130  * each frame in *filepos_array. */
03131 static int analyze_traj(qmdata_t *data, gmsdata *gms) {
03132   char buffer[BUFSIZ], nserch[BUFSIZ];
03133   char *line;
03134   long filepos;
03135   filepos = ftell(data->file);
03136 
03137   data->filepos_array = (long* )calloc(1, sizeof(long ));
03138 
03139   if (data->runtype==MOLFILE_RUNTYPE_OPTIMIZE ||
03140       data->runtype==MOLFILE_RUNTYPE_SADPOINT) {
03141     pass_keyline(data->file,
03142                    "PARAMETERS CONTROLLING GEOMETRY SEARCH", NULL);
03143     eatline(data->file, 2);
03144 
03145     GET_LINE(buffer, data->file);
03146     sscanf(buffer, "NSTEP  = %d", &data->max_opt_steps);
03147     eatline(data->file, 3);
03148     GET_LINE(buffer, data->file);
03149     sscanf(buffer, "OPTTOL = %f", &data->opt_tol);
03150 
03151     /* The $STATP options are followed by the coordinates 
03152      * but we can skip them here because we rewind after
03153      * get_guess_options() and try to read them in
03154      * get_traj_frame(). */
03155   }
03156   else if (data->runtype==MOLFILE_RUNTYPE_SURFACE) {
03157     if (pass_keyline(data->file,
03158                      "POTENTIAL SURFACE MAP INPUT", NULL)) {
03159       
03160       int coord1[2];
03161       int mplevel1=-1, mplevel2=-1, nstep1;
03162       float origin1, disp1;
03163       char runtype1[BUFSIZ], runtype2[BUFSIZ];
03164       char scftype1[BUFSIZ], scftype2[BUFSIZ];
03165       char dfttype1[BUFSIZ], dfttype2[BUFSIZ];
03166       char citype1[BUFSIZ],  citype2[BUFSIZ];
03167       char cctype1[BUFSIZ],  cctype2[BUFSIZ];
03168       char *tmp;
03169       int n;
03170         
03171       eatline(data->file, 1);
03172 
03173       GET_LINE(buffer, data->file);
03174       n=sscanf(buffer, " JOB 1 IS RUNTYP=%s SCFTYP=%s CITYP=%s",
03175                runtype1, scftype1, citype1);
03176       if (n==3) {
03177         GET_LINE(buffer, data->file);
03178         sscanf(buffer, " MPLEVL= %d CCTYP=%s, DFTTYP=%s\n",
03179                &mplevel1, dfttype1, cctype1);
03180         GET_LINE(buffer, data->file);
03181       }
03182       n=sscanf(buffer, " JOB 2 IS RUNTYP=%s SCFTYP=%s CITYP=%s",
03183                runtype2, scftype2, citype2);
03184       if (n==3) {
03185         GET_LINE(buffer, data->file);
03186         sscanf(buffer, " MPLEVL= %d CCTYP=%s, DFTTYP=%s\n",
03187                &mplevel2, dfttype2, cctype2);
03188         GET_LINE(buffer, data->file);
03189       }
03190 
03191       sscanf(buffer, " COORD 1 LYING ALONG ATOM PAIR %d %d",
03192              coord1, coord1+1);
03193       GET_LINE(buffer, data->file);
03194       tmp = strstr(buffer, "ORIGIN=") + 7;
03195       sscanf(tmp, "%f", &origin1);
03196       tmp = strstr(buffer, "DISPLACEMENT=") + 13;
03197       sscanf(tmp, "%f", &disp1);
03198       tmp = strstr(buffer, "AND") + 3;
03199       sscanf(tmp, "%d STEPS.", &nstep1);
03200       printf("gamessplugin) origin=%f, displacement=%f nstep=%d\n", origin1, disp1, nstep1);
03201     }
03202   }
03203   else if (data->runtype==MOLFILE_RUNTYPE_MEX) {
03204     char scftype1[BUFSIZ];
03205     char scftype2[BUFSIZ];
03206     rewind(data->file);
03207     if (!pass_keyline(data->file, "$MEX OPTIONS", NULL)) {
03208       printf("gamessplugin) No $MEX OPTIONS found!\n");
03209       return FALSE;
03210     }
03211     eatline(data->file, 2);
03212     GET_LINE(buffer, data->file);
03213     sscanf(strstr(buffer, "SCF1    =")+7, "%s", scftype1);
03214     sscanf(strstr(buffer, "SCF2   =")+7, "%s", scftype2);
03215     printf("gamessplugin) MEX SCF1=%s SCF2=%s\n", scftype1, scftype2);
03216 
03217   }
03218   else {
03219     /* We have just one frame */
03220     data->num_frames = 1;
03221     pass_keyline(data->file, "1 ELECTRON INTEGRALS",
03222                  "ENERGY COMPONENTS");
03223     data->filepos_array[0] = ftell(data->file);
03224 
03225     /* Check wether SCF has converged */
03226     if (pass_keyline(data->file,
03227                      "SCF IS UNCONVERGED, TOO MANY ITERATIONS",
03228                      "ENERGY COMPONENTS")==FOUND) {
03229       printf("gamessplugin) SCF IS UNCONVERGED, TOO MANY ITERATIONS\n");
03230       data->status = MOLFILE_QMSTATUS_SCF_NOT_CONV;
03231     } else {
03232       data->status = MOLFILE_QMSTATUS_OPT_CONV;
03233       fseek(data->file, data->filepos_array[0], SEEK_SET);
03234     }
03235 
03236     pass_keyline(data->file, "ENERGY COMPONENTS", NULL);
03237     data->end_of_traj = ftell(data->file);
03238 
03239     /* Allocate memory for the frame */
03240     data->qm_timestep = (qm_timestep_t *)calloc(1, sizeof(qm_timestep_t));
03241     memset(data->qm_timestep, 0, sizeof(qm_timestep_t));
03242     
03243     return TRUE;
03244   }
03245 
03246   printf("gamessplugin) Analyzing trajectory...\n");
03247   data->status = MOLFILE_QMSTATUS_UNKNOWN;
03248 
03249   while (1) {
03250     if (!fgets(buffer, sizeof(buffer), data->file)) break;
03251     line = trimleft(buffer);
03252 
03253     /* at this point we have to distinguish between
03254      * pre="27 JUN 2005 (R2)" and "27 JUN 2005 (R2)"
03255      * versions since the output format for geometry
03256      * optimizations has changed */
03257     if (gms->version==FIREFLY8POST6695){
03258       strcpy(nserch, "NSERCH=");
03259     }
03260     else if (gms->version==FIREFLY8PRE6695) {
03261       strcpy(nserch, "1NSERCH=");
03262     }
03263     else if (gms->version==GAMESSPRE20050627) {
03264       strcpy(nserch, "1NSERCH=");
03265     }
03266     else if (gms->version==GAMESSPOST20050627) {
03267       strcpy(nserch, "BEGINNING GEOMETRY SEARCH POINT NSERCH=");
03268     }
03269 
03270     if (strstr(line, nserch) ||
03271         strstr(line, "---- SURFACE MAPPING GEOMETRY") ||
03272         strstr(line, "MINIMUM ENERGY CROSSING POINT SEARCH") ||
03273         (data->runtype==MOLFILE_RUNTYPE_MEX && strstr(line, "NSERCH=")==line)) {
03274       printf("gamessplugin) %s", line);
03275 
03276       if (data->num_frames > 0) {
03277         data->filepos_array = (long*)realloc(data->filepos_array,
03278                                 (data->num_frames+1)*sizeof(long));
03279       }
03280       data->filepos_array[data->num_frames] = ftell(data->file);
03281       if (data->runtype==MOLFILE_RUNTYPE_SURFACE) {
03282         int ret = goto_keyline(data->file,
03283                                "ATOM      ATOMIC", "HAS ENERGY VALUE",
03284                                "---- SURFACE MAPPING GEOMETRY ----", NULL);
03285         if (ret>0 && ret<3 &&
03286             (have_keyline(data->file, "...... END OF ONE-ELECTRON INTEGRALS ......",
03287                           "---- SURFACE MAPPING GEOMETRY ----") ||
03288              have_keyline(data->file, "... DONE WITH POTENTIAL SURFACE SCAN",
03289                           "---- SURFACE MAPPING GEOMETRY ----"))) {
03290           data->num_frames++;          
03291         }
03292       }
03293       else if (pass_keyline(data->file, "COORDINATES OF",
03294                             "BEGINNING GEOMETRY SEARCH POINT NSERCH=")==FOUND)
03295       {
03296         /* Make sure that we have at least a complete coordinate
03297            block in order to consider this a new frame. */
03298         if (have_keyline(data->file, "INTERNUCLEAR DISTANCES",
03299                          "1 ELECTRON INTEGRALS") ||
03300             have_keyline(data->file, "1 ELECTRON INTEGRALS",
03301                          "BEGINNING GEOMETRY SEARCH POINT NSERCH=")) {
03302           data->num_frames++;
03303         }
03304       }
03305     }
03306     else if (strstr(line, "***** EQUILIBRIUM GEOMETRY LOCATED") ||
03307              strstr(line, "... DONE WITH POTENTIAL SURFACE SCAN")) {
03308       printf("gamessplugin) ==== End of trajectory (%d frames) ====\n",
03309              data->num_frames);
03310       data->status = MOLFILE_QMSTATUS_OPT_CONV;
03311       break;
03312     }
03313     else if (strstr(line, "***** FAILURE TO LOCATE STATIONARY POINT,")) {
03314       printf("gamessplugin) %s\n", line);
03315       if (strstr(strchr(line, ','), "SCF HAS NOT CONVERGED")) {
03316         data->status = MOLFILE_QMSTATUS_SCF_NOT_CONV;
03317         break;
03318       }
03319       else if (strstr(strchr(line, ','), "TOO MANY STEPS TAKEN")) {
03320         data->status = MOLFILE_QMSTATUS_OPT_NOT_CONV;
03321         break;
03322       }
03323     }
03324   }
03325   
03326   data->end_of_traj = ftell(data->file);
03327   fseek(data->file, filepos, SEEK_SET);
03328 
03329   if (data->status == MOLFILE_QMSTATUS_UNKNOWN) {
03330     /* We didn't find any of the regular key strings,
03331      * the run was most likely broken off and we have an
03332      * incomplete file. */
03333     data->status = MOLFILE_QMSTATUS_FILE_TRUNCATED;
03334   }
03335 
03336 
03337   /* Allocate memory for all frames */
03338   data->qm_timestep = (qm_timestep_t *)calloc(data->num_frames,
03339                                               sizeof(qm_timestep_t));
03340   memset(data->qm_timestep, 0, data->num_frames*sizeof(qm_timestep_t));
03341 
03342 
03343   if (data->status == MOLFILE_QMSTATUS_SCF_NOT_CONV ||
03344       data->status == MOLFILE_QMSTATUS_FILE_TRUNCATED) {
03345     return FALSE;  
03346   }
03347 
03348   return TRUE;
03349 }
03350 
03351 
03352 /***************************************************************
03353  *
03354  * Read the number of scf iterations and the scf energies
03355  * for the current timestep. 
03356  * Assumes that the file pointer is somewhere before this:
03357  * ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE    DIIS ERROR
03358  * 1  0  0      -39.7266993475   -39.7266993475   0.000000118   0.000000000
03359  * 2  1  0      -39.7266991566     0.0000001909   0.000000032   0.000000000
03360  * ...
03361  * then it reads the block up to the next blank line.
03362  * The second argument is a pointer to the qm timestep you want to
03363  * store the data in. Memory for the scfenergies will be
03364  * allocated.
03365  *
03366  ***************************************************************/
03367 static int get_scfdata(qmdata_t *data, qm_timestep_t *ts) {
03368   char buffer[BUFSIZ];
03369   char word[3][BUFSIZ];
03370   long filepos;
03371   int i, epos = -1;
03372   int numread, numiter=0, dum, dum2;
03373   char *line;
03374   float dumf;
03375   filepos = ftell(data->file);
03376 
03377   for (i=0; i<3; i++) word[i][0] = '\0';
03378 
03379   if (!goto_keyline(data->file, "ITER EX", "ITER     TOTAL",
03380                      "ITER    TOTAL", NULL)) {
03381     fseek(data->file, filepos, SEEK_SET);
03382     ts->num_scfiter = 0;
03383     return FALSE;
03384   }
03385 
03386   /* determine in which column the energy is stored */
03387   GET_LINE(buffer, data->file);
03388   numread = sscanf(buffer, "%*s %s %s %s",
03389                    &word[0][0], &word[1][0], &word[2][0]);
03390   for (i=0; i<numread; i++) {
03391     if (!strcmp(&word[i][0], "TOTAL")) epos = i+1;
03392   }
03393    
03394   if (epos<0) {
03395     fseek(data->file, filepos, SEEK_SET);
03396     ts->num_scfiter = 0;
03397     return FALSE;
03398   }
03399 
03400   /* store current file position since we first have to count
03401    * the iterations */
03402   filepos = ftell(data->file);
03403 
03404   /* read until the next blank line and count the iterations */
03405   do {
03406     GET_LINE(buffer, data->file);
03407     line = trimleft(buffer);
03408     numread = sscanf(line,"%d %d %*d %*f", &dum, &dum2);
03409     if (numread==2) numiter++;
03410   } while (strlen(line)>2);
03411 
03412   printf("gamessplugin) %d SCF iterations\n", numiter);
03413 
03414   /* go back and read energies */
03415   fseek(data->file, filepos, SEEK_SET);
03416   
03417 
03418   /* allocate memory for scfenergy array */
03419   ts->scfenergies = (double *)calloc(numiter,sizeof(double));
03420   
03421   i=0;
03422   do {
03423     GET_LINE(buffer, data->file);
03424     line = trimleft(buffer);
03425     numread = sscanf(line,"%d %f %*i %*f", &dum, &dumf);
03426     if (numread==2) {
03427       switch (epos) {
03428       case 1:
03429         sscanf(buffer,"%*d %lf", ts->scfenergies+i);
03430         break;
03431       case 2:
03432         sscanf(buffer,"%*d %*d %lf", ts->scfenergies+i);
03433         break;
03434       case 3:
03435         sscanf(buffer,"%*d %*d %*d %lf", ts->scfenergies+i);
03436         break;
03437       }
03438       i++;
03439     }
03440   } while (strlen(line)>2);
03441 
03442 #if 0
03443   for (i=0; i<numiter; i++) {
03444     printf("scfenergies[%d] = %f\n", i, ts->scfenergies[i]);
03445   }
03446 #endif
03447 
03448   ts->num_scfiter = numiter;
03449   
03450   return TRUE;
03451 }
03452 
03453 
03454 /*********************************************************
03455  *
03456  * Reads a set of wavefunctions for the current timestep.
03457  * These are typically the alpha and beta spin wavefunctions
03458  * or the MCSCF natural and optimized orbitals or the GVB
03459  * canonical orbitals and geminal pairs.
03460  *
03461  **********************************************************/
03462 static int check_add_wavefunctions(qmdata_t *data,
03463                                    qm_timestep_t *ts) {
03464   qm_wavefunction_t *wavef;
03465   int i, n=1;
03466 
03467   if (data->scftype==MOLFILE_SCFTYPE_UHF || 
03468       data->scftype==MOLFILE_SCFTYPE_GVB ||
03469       data->scftype==MOLFILE_SCFTYPE_MCSCF) {
03470     /* Try to read second wavefunction
03471      * (spin beta or GI orbitals or MCSCF optimized orbs) */
03472     n = 2;
03473   }
03474 
03475   for (i=0; i<n; i++) {
03476     /* Allocate memory for new wavefunction */
03477     wavef = add_wavefunction(ts);
03478 
03479     /* Try to read wavefunction and orbital energies */
03480     if (get_wavefunction(data, ts, wavef) == FALSE) {
03481       /* Free the last wavefunction again. */
03482       del_wavefunction(ts);
03483 #ifdef DEBUGGING
03484       printf("gamessplugin) No canonical wavefunction present for timestep %d\n", data->num_frames_read);
03485 #endif
03486       break;
03487 
03488     } else {
03489       char action[32];
03490       char spinstr[32];
03491       strcpy(spinstr, "");
03492       if (data->scftype==MOLFILE_SCFTYPE_UHF) {
03493         if (wavef->spin==SPIN_BETA) {
03494           strcat(spinstr, "spin  beta, ");
03495         } else {
03496           strcat(spinstr, "spin alpha, ");
03497         }
03498       }
03499       
03500       /* The last SCF energy is the energy of this electronic state */
03501       if (ts->scfenergies) {
03502         wavef->energy = ts->scfenergies[ts->num_scfiter-1];
03503       } else {
03504         wavef->energy = 0.f;
03505       }
03506       
03507       /* Multiplicity */
03508       wavef->mult = data->multiplicity;
03509       
03510 
03511       /* String telling wether wavefunction was added, updated
03512        * or ignored. */
03513       strcpy(action, "added");
03514 
03515       /* If there exists a canonical wavefunction of the same spin
03516        * we'll replace it */
03517       if (ts->numwave>1 && wavef->type==MOLFILE_WAVE_CANON) {
03518         int i, found =-1;
03519         for (i=0; i<ts->numwave-1; i++) {
03520           if (ts->wave[i].type==wavef->type &&
03521               ts->wave[i].spin==wavef->spin &&
03522               ts->wave[i].exci==wavef->exci &&
03523               !strncmp(ts->wave[i].info, wavef->info, MOLFILE_BUFSIZ)) {
03524             found = i;
03525             break;
03526           }
03527         }
03528         if (found>=0) {
03529           /* If the new wavefunction has more orbitals we 
03530            * replace the old one for this step. */
03531           if (wavef->num_orbitals > 
03532               ts->wave[found].num_orbitals) {
03533             /* Replace existing wavefunction for this step */
03534             replace_wavefunction(ts, found);
03535             sprintf(action, "%d updated", found);
03536           } else {
03537             /* Delete last wavefunction again */
03538             del_wavefunction(ts);
03539             sprintf(action, "matching %d ignored", found);
03540           }
03541           wavef = &ts->wave[ts->numwave-1];
03542         }
03543       }
03544 
03545       printf("gamessplugin) Wavefunction %s (%s):\n", action, wavef->info);
03546       printf("gamessplugin)   %d orbitals, %sexcitation %d, multiplicity %d\n",
03547              wavef->num_orbitals, spinstr, wavef->exci, wavef->mult);
03548     }
03549   }
03550 
03551   return i;
03552 }
03553 
03554 
03555 /*********************************************************
03556  *
03557  * Finds the next wavefunction, determines its type by
03558  * analyzing the keystring and reads in the wavefunction
03559  * coefficients.
03560  *
03561  **********************************************************/
03562 static int get_wavefunction(qmdata_t *data, qm_timestep_t *ts,
03563                             qm_wavefunction_t *wf)
03564 {
03565   float *orb_enocc;
03566   float *wave_coeff;
03567   char buffer[BUFSIZ];
03568   char word[6][BUFSIZ];
03569   int num_orbitals = 0;
03570   int i = 0, num_values = 0;
03571   long filepos;
03572   char *line;
03573   int have_orbenocc = 0;
03574   int n[5];
03575 
03576   buffer[0] = '\0';
03577   for (i=0; i<6; i++) word[i][0] = '\0';
03578 
03579   if (wf == NULL) {
03580     PRINTERR;       
03581     return FALSE;
03582   }
03583 
03584   wf->has_occup = FALSE;
03585   wf->has_orben = FALSE;
03586   wf->type = MOLFILE_WAVE_UNKNOWN;
03587   wf->spin = SPIN_ALPHA;
03588   wf->exci = 0;
03589   strncpy(wf->info, "unknown", MOLFILE_BUFSIZ);
03590 
03591   /*
03592    * Scan for something like this:
03593 
03594           ------------------
03595           MOLECULAR ORBITALS     <<--- sometimes EIGENVECTORS
03596           ------------------
03597 
03598                       1          2          3          4          5
03599                   -11.0297    -0.9121    -0.5205    -0.5205    -0.5205  <<-- orbital energies (or occupancies)                     A          A          A          A          A   
03600     1  C  1  S    0.991925   0.221431   0.000006  -0.000001   0.000002
03601     2  C  1  S    0.038356  -0.627585  -0.000021   0.000003  -0.000006
03602     3  C  1  X    0.000000  -0.000004   0.338169  -0.030481  -0.460283
03603      ...
03604 
03605                      6          7          8          9
03606                     0.7192     0.7192     0.7193     0.7611
03607                      A          A          A          A   
03608     1  C  1  S    0.000028   0.000012   0.000092   0.252320
03609     2  C  1  S   -0.000183  -0.000077  -0.000594  -1.632834
03610     3  C  1  X   -0.890147   0.062618   0.654017  -0.000154
03611       ...
03612 
03613 
03614      ----------------------------------------------------------------
03615      PROPERTY VALUES FOR THE RHF   SELF-CONSISTENT FIELD WAVEFUNCTION
03616      ----------------------------------------------------------------
03617   * 
03618   */
03619 
03620   /* Remember position in order to go back if no wave function was found */
03621   filepos = ftell(data->file);
03622 
03623   do {
03624     GET_LINE(buffer, data->file);
03625 
03626     line = trimleft(trimright(buffer));
03627 
03628     if      (!strcmp(line, "----- ALPHA SET -----")) {
03629       wf->type = MOLFILE_WAVE_CANON;
03630       strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
03631       pass_keyline(data->file, "EIGENVECTORS", NULL);
03632     }
03633     else if (!strcmp(line, "----- BETA SET -----")) {
03634       wf->type = MOLFILE_WAVE_CANON;
03635       wf->spin = SPIN_BETA;
03636       strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
03637       pass_keyline(data->file, "EIGENVECTORS", NULL);
03638     }
03639     else if (!strcmp(line, "****** BETA ORBITAL LOCALIZATION *****")) {
03640       wf->spin = SPIN_BETA;
03641     }
03642     else if (!strcmp(line, "EIGENVECTORS")) {
03643       wf->type = MOLFILE_WAVE_CANON;
03644       strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
03645     }
03646     else if (!strcmp(line, "MOLECULAR ORBITALS")) {
03647       wf->type = MOLFILE_WAVE_CANON;
03648       strncpy(wf->info, "canonical", MOLFILE_BUFSIZ);
03649     }
03650     else if (!strcmp(line, "THE BOYS LOCALIZED ORBITALS ARE")) {
03651       wf->type = MOLFILE_WAVE_BOYS;
03652       strncpy(wf->info, "Boys localized", MOLFILE_BUFSIZ);
03653     }
03654     else if (!strcmp(line, "THE PIPEK-MEZEY POPULATION LOCALIZED ORBITALS ARE")) {
03655       wf->type = MOLFILE_WAVE_PIPEK;
03656       strncpy(wf->info, "Pipek-Mezey localized", MOLFILE_BUFSIZ);
03657     }
03658     else if (!strcmp(line, "EDMISTON-RUEDENBERG ENERGY LOCALIZED ORBITALS")) {
03659       wf->type = MOLFILE_WAVE_RUEDEN;
03660       strncpy(wf->info, "Ruedenberg localized", MOLFILE_BUFSIZ);
03661     }
03662     else if (!strcmp(line, "GI ORBITALS")) {
03663       wf->type = MOLFILE_WAVE_GEMINAL;
03664       strncpy(wf->info, "GVB geminal pairs", MOLFILE_BUFSIZ);
03665     }
03666     else if (!strcmp(line, "MCSCF NATURAL ORBITALS")) {
03667       wf->type = MOLFILE_WAVE_MCSCFNAT;
03668       strncpy(wf->info, "MCSCF natural orbitals", MOLFILE_BUFSIZ);
03669     }
03670     else if (!strcmp(line, "MCSCF OPTIMIZED ORBITALS")) {
03671       wf->type = MOLFILE_WAVE_MCSCFOPT;
03672       strncpy(wf->info, "MCSCF optimized orbitals", MOLFILE_BUFSIZ);
03673     }
03674     else if (!strcmp(line, "NATURAL ORBITALS IN ATOMIC ORBITAL BASIS")) {
03675       wf->type = MOLFILE_WAVE_CINATUR;
03676       strncpy(wf->info, "CI natural orbitals", MOLFILE_BUFSIZ);
03677     }
03678     else if (!strcmp(line, "CIS NATURAL ORBITALS")) {
03679       wf->type = MOLFILE_WAVE_CINATUR;
03680       strncpy(wf->info, "CIS natural orbitals", MOLFILE_BUFSIZ);
03681     }
03682     // FOR PCGAMESS/FIREFLY
03683     else if (!strcmp(line, "-MCHF- NATURAL ORBITALS")) {
03684       wf->type = MOLFILE_WAVE_MCSCFNAT;
03685       strncpy(wf->info, "MCSCF natural orbitals", MOLFILE_BUFSIZ);
03686     }
03687     else if (!strcmp(line, "-MCHF- OPTIMIZED ORBITALS")) {
03688       wf->type = MOLFILE_WAVE_MCSCFOPT;
03689       strncpy(wf->info, "MCSCF optimized orbitals", MOLFILE_BUFSIZ);
03690     }
03691     else if (!strcmp(line, "ZERO-ORDER QDPT NATURAL ORBITALS")){
03692       //Not yet implemented
03693     }
03694 
03695   } while(wf->type==MOLFILE_WAVE_UNKNOWN &&
03696           strcmp(line, "ENERGY COMPONENTS") &&
03697           strcmp(line, "***** EQUILIBRIUM GEOMETRY LOCATED *****") &&
03698           strcmp(line, "**** THE GEOMETRY SEARCH IS NOT CONVERGED! ****"));
03699 
03700   /* If we reach the last line of the rhf section without finding 
03701    * one of the keywords marking the beginning of a wavefunction
03702    * table then we return.*/
03703   if (wf->type==MOLFILE_WAVE_UNKNOWN) {
03704 #ifdef DEBUGGING
03705     printf("gamessplugin) get_wavefunction(): No wavefunction found!\n");
03706 #endif
03707     fseek(data->file, filepos, SEEK_SET);
03708     return FALSE;
03709   }
03710 
03711   /* Reserve space for arrays storing wavefunction and orbital
03712    * energies. we reserve the max. space (num_orbitals==wavef_size)
03713    * for now and realloc later if, we have less orbitals. */
03714   wave_coeff = (float *)calloc(data->wavef_size*data->wavef_size,
03715                                   sizeof(float)); 
03716   //printf("gamessplugin) wave_coeff() allocated\n");
03717 
03718   if (wave_coeff == NULL) {
03719     PRINTERR;       
03720     return FALSE;
03721   }
03722 
03723   /* orbital energies/occupancies */  
03724   orb_enocc = (float *)calloc(data->wavef_size, sizeof(float));
03725   //printf("gamessplugin) orb_enocc() allocated\n");
03726 
03727   if (orb_enocc == NULL) {
03728     free(orb_enocc);
03729     PRINTERR; 
03730     return FALSE;
03731   }
03732 
03733 
03734   /* store the coeficient pointer */
03735   wf->wave_coeffs  = wave_coeff;
03736 
03737   /* depending on the wavefunction type the line after the
03738      orbital index stores the orbital occupancies or energies */
03739   if (wf->type == MOLFILE_WAVE_CINATUR  ||
03740       wf->type == MOLFILE_WAVE_MCSCFNAT) {
03741     wf->orb_occupancies = orb_enocc;
03742     wf->has_occup = TRUE;
03743   } else {
03744     wf->orb_energies    = orb_enocc;
03745     wf->has_orben = TRUE;
03746   }
03747 
03748   /* skip the next line which here is typically "-------" */
03749   eatline(data->file, 1);
03750 
03751 
03752   while (1) {
03753     int nr, over=0;
03754     float coeff[5], enocc[5];
03755 
03756     if (wf->type == MOLFILE_WAVE_GEMINAL) {
03757       /* Skip over "PAIR x" header line */
03758       pass_keyline(data->file, "PAIR ", NULL);
03759     }
03760 
03761     eatwhitelines(data->file);
03762     filepos = ftell(data->file);
03763 
03764     /* Parse the orbital indexes */
03765     GET_LINE(buffer, data->file);
03766     num_values = sscanf(buffer, "%d %d %d %d %d",
03767                           &n[0], &n[1], &n[2], &n[3], &n[4]);
03768 
03769     /* If there are no orbital indexes then this must be the
03770      * end of the wavefunction coefficient table. */
03771     if (!num_values) {
03772       fseek(data->file, filepos, SEEK_SET);
03773       break;
03774     }
03775 
03776     eatwhitelines(data->file);
03777 
03778     /* Read first line of orbital energies/occupancies */
03779     filepos = ftell(data->file);
03780     GET_LINE(buffer, data->file);
03781     have_orbenocc = sscanf(buffer,"%f %f %f %f %f", &enocc[0],
03782                         &enocc[1], &enocc[2], &enocc[3], &enocc[4]);
03783 
03784     /* Make sure this is not the first line containing coeffs */
03785     nr = sscanf(buffer, " 1 %*s 1 %*s %f %f %f %f %f",
03786                &coeff[0], &coeff[1], &coeff[2], &coeff[3], &coeff[4]);
03787     if (nr==num_values) have_orbenocc = 0;
03788 
03789     if (have_orbenocc) {
03790       /* store the orbital energies in the appropriate arrays 
03791        * read them until we encounter an empty string */
03792       for(i=0; i<num_values; i++) {
03793         orb_enocc[i] = enocc[i];
03794       }
03795       
03796 
03797       /* If we are in the first block we have to distinguish 
03798          between energies and occupancies */
03799       if (wf->type  == MOLFILE_WAVE_MCSCFNAT &&
03800           orb_enocc == wf->orb_occupancies   &&
03801           enocc[0] <= 0.f) {
03802         wf->orb_occupancies = NULL;
03803         wf->has_occup = FALSE;
03804         wf->orb_energies    = orb_enocc;
03805         wf->has_orben = TRUE;
03806       }
03807 
03808       /* increase orbital energy pointer */
03809       orb_enocc = orb_enocc+5;
03810     }      
03811     else {
03812       /* No orbital energies present, go back one line */
03813       fseek(data->file, filepos, SEEK_SET);
03814     }
03815 
03816     num_orbitals += num_values;
03817 
03818     /* Find first line containing coefficients */
03819     filepos = ftell(data->file);
03820     while (fgets(buffer, sizeof(buffer), data->file)) {
03821       trimleft(buffer);
03822       if (strstr(line, "ENERGY COMPONENTS") ||
03823           strstr(line, "---") ||
03824           strstr(line, "...")) {
03825         over = 1; break;
03826       }
03827 
03828       nr = sscanf(buffer, " 1 %*s 1 %*s %f %f %f %f %f",
03829              &coeff[0], &coeff[1], &coeff[2], &coeff[3], &coeff[4]);
03830       if (nr==num_values) break;
03831       filepos = ftell(data->file);
03832     }
03833     fseek(data->file, filepos, SEEK_SET);
03834 
03835     if (over) break;
03836 
03837 
03838     /* Read the wave function coefficient block for up to 5
03839      * orbitals per line. */
03840     if (!read_coeff_block(data->file, data->wavef_size,
03841                           wave_coeff, data->angular_momentum)) {
03842       printf("gamessplugin) Wavefunction coefficient block truncated or ill formatted!\n");
03843       data->status = MOLFILE_QMSTATUS_FILE_TRUNCATED;
03844       return FALSE;
03845     }
03846 
03847 
03848     /* move wavefunction pointer to start of next five orbitals */
03849     if (wf->type == MOLFILE_WAVE_GEMINAL) {
03850       wave_coeff = wave_coeff + 2*data->wavef_size;
03851     } else {
03852       wave_coeff = wave_coeff + 5*data->wavef_size;
03853     }
03854   }
03855 
03856 
03857   if (!num_orbitals) {
03858     printf("gamessplugin) No orbitals in wavefunction!\n");
03859     return FALSE;
03860   }
03861 
03862   /* resize the array to the actual number of read orbitals */
03863   if (data->wavef_size!=num_orbitals) {
03864 
03865     if (wf->has_occup) {
03866       wf->orb_occupancies = (float *)realloc(wf->orb_occupancies,
03867                num_orbitals*sizeof(float));
03868     }
03869     if (wf->has_orben) {
03870       wf->orb_energies = (float *)realloc(wf->orb_energies,
03871                num_orbitals*sizeof(float));
03872     }
03873 
03874     wf->wave_coeffs  = (float *)realloc(wf->wave_coeffs, data->wavef_size*
03875                                         num_orbitals*sizeof(float)); 
03876   }
03877 
03878   /* In case MCSCF natural orbitals are present, then GAMESS
03879      prints the orbital energy for the core orbitals and the
03880      occupancy for the other orbitals. We zero out the non-core
03881      energies to prevent confusion. The orbital occupancies 
03882      are read separately elsewhere. */
03883   if (wf->type == MOLFILE_WAVE_MCSCFNAT &&
03884       wf->has_orben == TRUE ) {
03885     wf->orb_occupancies = (float *)calloc(num_orbitals, sizeof(float));
03886 
03887     for (i=0; i<data->mcscf_num_core; i++) {
03888       wf->orb_occupancies[i] = 2.f;
03889     }
03890 
03891     for (i=data->mcscf_num_core; i<num_orbitals; i++) {
03892       wf->orb_occupancies[i] = wf->orb_energies[i];
03893       wf->orb_energies[i] = 0.f;
03894     }
03895 
03896     wf->has_occup = TRUE;
03897   }
03898 
03899 
03900   /* store the number of orbitals read in */
03901   wf->num_orbitals  = num_orbitals;
03902 
03903   return TRUE;
03904 }
03905 
03906 
03907 /* Read the wave function coefficient block for up to 5
03908  * orbitals per line:
03909  *  1  C  1  S    0.989835   0.155361   0.000000  -0.214258   0.000000
03910  *  2  C  1  S    0.046228  -0.548915   0.000000   0.645267   0.000000
03911  *  3  C  1  X    0.000000   0.000000   1.030974   0.000000   0.000000
03912  */
03913 static int read_coeff_block(FILE *file, int wavef_size,
03914                             float *wave_coeff, int *angular_momentum) {
03915   int i, j;
03916   int truncated = 0;
03917   char buffer[BUFSIZ];
03918 
03919   /* Read a line with coefficients for up to 5 orbitals
03920    * for each cartesian basis function. */
03921   for (i=0; i<wavef_size; i++) {
03922     char type[BUFSIZ];
03923     float coeff[5];
03924     int num_values = 0;
03925     
03926     GET_LINE(buffer, file);
03927     
03928     /* read in the wavefunction coefficients for 5
03929      * orbitals at a time line by line */
03930     num_values = sscanf(buffer,"%*5i%*4s%*2i%4s %f %f %f %f %f", 
03931                         type, &coeff[0], &coeff[1], &coeff[2],
03932                         &coeff[3], &coeff[4]);
03933     
03934     if (num_values==0) {
03935       /* The file must have been truncated! */
03936       truncated = 1;
03937       break;
03938     }
03939 
03940     angular_momentum_expon(&angular_momentum[3*i], type);
03941    
03942     /* Each orbital has data->wavef_size entries, 
03943      * hence we have to use this number as offset when storing 
03944      * them in groups of five. */
03945     for (j=0 ; j<num_values-1; j++) {
03946       wave_coeff[j*wavef_size+i] = coeff[j];
03947     }
03948   }
03949   
03950   if (truncated) return 0;
03951   
03952   return 1;
03953 }
03954 
03955 /* Read the population analysis section.
03956  * Currently we parse only the Mulliken and Lowdin charges
03957  * but we might want to add support for population analysis. */
03958 static int get_population(qmdata_t *data, qm_timestep_t *ts) {
03959   int i;
03960   char buffer[BUFSIZ];
03961   long filepos;
03962   ts->have_mulliken = FALSE;
03963   ts->have_lowdin   = FALSE;
03964   filepos = ftell(data->file);
03965 
03966   if (pass_keyline(data->file,
03967                      "TOTAL MULLIKEN AND LOWDIN ATOMIC POPULATIONS",
03968                      "NSERCH=") != FOUND) {
03969     fseek(data->file, filepos, SEEK_SET);
03970     return FALSE;
03971   }
03972 
03973   /* Read Mulliken charges if present */
03974   ts->mulliken_charges = 
03975     (double *)calloc(data->numatoms, sizeof(double));
03976 
03977   if (!ts->mulliken_charges) {
03978     PRINTERR; 
03979     return FALSE;
03980   }
03981 
03982   ts->lowdin_charges = 
03983     (double *)calloc(data->numatoms, sizeof(double));
03984 
03985   if (!ts->lowdin_charges) {
03986     free(ts->mulliken_charges);
03987     ts->mulliken_charges = NULL;
03988     PRINTERR; 
03989     return FALSE;
03990   }
03991   
03992   eatline(data->file, 1);
03993   
03994   for (i=0; i<data->numatoms; i++) {
03995     int n;
03996     float mullpop, mullcharge, lowpop, lowcharge;
03997     GET_LINE(buffer, data->file);
03998     n = sscanf(buffer,"%*i %*s %f %f %f %f",
03999                    &mullpop, &mullcharge, &lowpop, &lowcharge);
04000     if (n!=4) {
04001       free(ts->mulliken_charges);
04002       free(ts->lowdin_charges);
04003       ts->mulliken_charges = NULL;
04004       ts->lowdin_charges   = NULL;
04005       return FALSE;
04006     }
04007     ts->mulliken_charges[i] = mullcharge;
04008     ts->lowdin_charges[i]   = lowcharge;
04009   }
04010 
04011   if (i!=data->numatoms) {
04012     free(ts->mulliken_charges);
04013     free(ts->lowdin_charges);
04014     ts->mulliken_charges = NULL;
04015     ts->lowdin_charges   = NULL;
04016     return FALSE;
04017   }
04018 
04019   ts->have_mulliken = TRUE;
04020   ts->have_lowdin   = TRUE;
04021   return TRUE;
04022 }
04023 
04024 
04025 /* Read ESP charges.
04026  * XXX Right now we don't distinguish between different type of
04027  * ESP-style charges (CHELPG, CONNOLLY, GEODESIC). 
04028  * This could be solved by reading in the PTSEL keyword in
04029  * the $PDC group. */
04030 static int get_esp_charges(qmdata_t *data) {
04031   int i;
04032   char buffer[BUFSIZ];
04033   long filepos;
04034   /* Store charges in last timestep */
04035   qm_timestep_t *ts = &data->qm_timestep[data->num_frames-1];
04036 
04037   ts->have_esp = FALSE;
04038   filepos = ftell(data->file);
04039 
04040   if (pass_keyline(data->file,
04041            "ATOM                CHARGE    E.S.D.",
04042            "...... END OF PROPERTY EVALUATION ") != FOUND) {
04043     fseek(data->file, filepos, SEEK_SET);
04044     return FALSE;
04045   }
04046 
04047   /* Read ESP charges if present */
04048   ts->esp_charges = 
04049     (double *)calloc(data->numatoms, sizeof(double));
04050 
04051   if (ts->esp_charges == NULL) {
04052     PRINTERR; 
04053     return FALSE;
04054   }
04055 
04056   eatline(data->file, 1);
04057 
04058   for (i=0; i<data->numatoms; i++) {
04059     int n;
04060     double charge;
04061     GET_LINE(buffer, data->file);
04062     n = sscanf(buffer,"%*s %lf ", &charge);
04063     if (n!=1) return FALSE;
04064     ts->esp_charges[i] = charge;
04065   }
04066 
04067   if (i!=data->numatoms) {
04068     
04069     return FALSE;
04070   }
04071 
04072   ts->have_esp = TRUE;
04073   return TRUE;
04074 }
04075 
04076 
04077 /* Read the energy gradient (=force) for each atom */
04078 static int get_gradient(qmdata_t *data, qm_timestep_t *ts) {
04079   int numgrad=0;
04080   int numread;
04081   char buffer[BUFSIZ];
04082   long filepos;
04083 
04084   buffer[0] = '\0';
04085 
04086   /* remember position in order to go back if no forces were found */
04087   filepos = ftell(data->file);
04088 
04089   /* look for GRADIENT section */
04090   if (goto_keyline(data->file, "GRADIENT (HARTREE",
04091                 "***** EQUILIBRIUM GEOMETRY LOCATED", 
04092                 " BEGINNING GEOMETRY SEARCH", NULL) != FOUND) {
04093     fseek(data->file, filepos, SEEK_SET);
04094     return FALSE;
04095   }
04096 
04097   eatline(data->file, 4);
04098 
04099   ts->gradient = (float *)calloc(3*data->numatoms, sizeof(float));
04100 
04101   if (ts->gradient == NULL) {
04102     PRINTERR;       
04103     fseek(data->file, filepos, SEEK_SET);
04104     return FALSE;
04105   }
04106 
04107   /* read the gradient table */
04108   do {
04109     int i;
04110     float dx, dy, dz;
04111     GET_LINE(buffer, data->file);
04112     numread = sscanf(buffer, "%d %*s %*f %f %f %f", &i, &dx, &dy, &dz);
04113     if (numread==4) {
04114       ts->gradient[3*(i-1)  ] = dx;
04115       ts->gradient[3*(i-1)+1] = dy;
04116       ts->gradient[3*(i-1)+2] = dz;
04117       numgrad++;
04118     }
04119   } while(numread==4);
04120 
04121   if (numgrad!=data->numatoms) {
04122     printf("gamessplugin) Found %d gradients for %d atoms!\n",
04123            numgrad, data->numatoms);
04124     fseek(data->file, filepos, SEEK_SET);
04125     return FALSE;
04126   }
04127 
04128   return TRUE;
04129 }
04130 
04131 
04132 /* Read energy gradients from final step
04133  * (different format in final step) */
04134 static int get_final_gradient(qmdata_t *data, qm_timestep_t *ts) {
04135   int numgrad=0;
04136   int numread;
04137   char buffer[BUFSIZ];
04138   long filepos;
04139 
04140   /* remember position in order to go back at the end */
04141   filepos = ftell(data->file);
04142 
04143   if (pass_keyline(data->file,
04144                    "ATOM                 E'X", NULL) != FOUND) {
04145     fseek(data->file, filepos, SEEK_SET);
04146     return FALSE;
04147   }
04148 
04149   ts->gradient = (float *)calloc(3*data->numatoms, sizeof(float));
04150 
04151   if (ts->gradient == NULL) {
04152     PRINTERR;       
04153     fseek(data->file, filepos, SEEK_SET);
04154     return FALSE;
04155   }
04156 
04157   /* read the gradient table */
04158   do {
04159     int i;
04160     float dx, dy, dz;
04161     GET_LINE(buffer, data->file);
04162     numread = sscanf(buffer, "%d %*s %f %f %f", &i, &dx, &dy, &dz);
04163     if (numread==4) {
04164       ts->gradient[3*(i-1)  ] = dx;
04165       ts->gradient[3*(i-1)+1] = dy;
04166       ts->gradient[3*(i-1)+2] = dz;
04167       numgrad++;
04168     }
04169   } while(numread==4);
04170 
04171   /* go back to where search started */
04172   fseek(data->file, filepos, SEEK_SET);
04173 
04174   if (numgrad!=data->numatoms) {
04175     printf("gamessplugin) Number of gradients != number of atoms!\n");
04176     return FALSE;
04177   }
04178 
04179   return TRUE;
04180 }
04181 
04182 
04183 /***********************************************************
04184  *
04185  * Read in wavenumbers and intensities of the normal modes
04186  *
04187  **********************************************************/
04188 static int get_normal_modes(qmdata_t *data) {
04189   char buffer[BUFSIZ];
04190   int i = 0, k = 0, j = 0;
04191   double entry[6]; 
04192   char *token;
04193 
04194   if (!pass_keyline(data->file, "NORMAL COORDINATE ANALYSIS", NULL)) {
04195     return FALSE;
04196   }
04197 
04198   /* initialize array */
04199   memset(entry, 0, sizeof(entry));
04200 
04201     
04202   /* allocate memory for arrays */
04203   data->wavenumbers = 
04204     (float *)calloc(data->numatoms*3,sizeof(float));
04205   if (data->wavenumbers==NULL) {
04206     PRINTERR;
04207     return FALSE;
04208   }
04209 
04210   data->intensities = 
04211     (float *)calloc(data->numatoms*3,sizeof(float));
04212   if (data->intensities==NULL) {
04213     PRINTERR;
04214     return FALSE;
04215   }
04216 
04217   data->imag_modes = 
04218     (int *)calloc(data->numatoms*3,sizeof(int));
04219   if (data->imag_modes==NULL) {
04220     PRINTERR;
04221     return FALSE;
04222   }
04223 
04224   data->normal_modes = 
04225     (float *)calloc((data->numatoms*3)*(data->numatoms*3),
04226                      sizeof(float));
04227   if (data->normal_modes==NULL) {
04228     PRINTERR;
04229     return FALSE;
04230   }
04231 
04232 
04233   /* Example:
04234    *                       1           2           3           4           5
04235    *    FREQUENCY:      4805.54 I   1793.00 I   1317.43 I      2.13        1.99
04236    * REDUCED MASS:      1.09875     1.00796     1.18876     1.04129     1.03921
04237    * IR INTENSITY:    105.04109     2.80788     3.07217     0.01689     0.01587
04238    *
04239    * 1   CARBON       X  0.00000000  0.00000000 -0.11767262 -0.05407389  0.00000000
04240    *                  Y  0.00000000 -0.00345200  0.00000000  0.00000000 -0.05241437
04241    *                  Z -0.08676597  0.00000000  0.00000000  0.00000000  0.00000000
04242    */
04243   for (i=0; i<ceil(data->numatoms*3/5.f); i++) {
04244     int numread = 0;
04245 
04246     if (!goto_keyline(data->file, "FREQUENCY:", NULL)) {
04247       break;
04248     }
04249     GET_LINE(buffer, data->file);
04250 
04251     /* Scan the frequencies; 
04252      * If there are imaginary modes present then the
04253      * frequency is followed by the 'I' which represents
04254      * an additional char token in the line. */
04255 
04256     /* Skip first token "FREQUENCY:" */
04257     token = strtok(buffer, " \t\r\n");
04258     
04259     /* Walk through the remaining tokens */
04260     while ((token = strtok(NULL, " \t\r\n")) != NULL) {
04261       /* Check if token is 'I'.
04262        * If yes, mark previous mode as imaginary. */
04263       if (*token=='I') {
04264         data->imag_modes[data->nimag] = numread-1;
04265         data->nimag++;
04266       } else {
04267         /* save only the first 5 modes - there NEVER should
04268          * be more in any case, but just to make sure
04269          * we don't overrun the array */
04270         if (numread<5) {
04271           data->wavenumbers[i*5+numread] = atof(token);
04272           numread++;
04273         }
04274       }
04275     }
04276 
04277     eatline(data->file, 1);
04278 
04279     /* Read the IR INTENSITIES */
04280     GET_LINE(buffer, data->file);
04281     numread = sscanf(buffer,"%*s %*s %lf %lf %lf %lf %lf", &entry[0],
04282                      &entry[1], &entry[2], &entry[3], &entry[4]);
04283  
04284     for (k=0; k<numread; k++) {
04285       data->intensities[i*5+k] = entry[k]; 
04286     }
04287 
04288     eatline(data->file, 1);
04289 
04290     /* Read the normal mode vectors */
04291     for (k=0; k<data->numatoms; k++) {
04292       /* x */
04293       GET_LINE(buffer, data->file);
04294       numread = sscanf(buffer,"%*s %*s %*s %lf %lf %lf %lf %lf",
04295              &entry[0], &entry[1], &entry[2], &entry[3], &entry[4]);
04296 
04297       for (j=0; j<numread; j++) {
04298         data->normal_modes[3*k + (i*5+j)*3*data->numatoms] = 
04299           entry[j];
04300       }
04301 
04302       /* y */
04303       GET_LINE(buffer, data->file);
04304       numread = sscanf(buffer,"%*s %lf %lf %lf %lf %lf", &entry[0],
04305                        &entry[1],&entry[2], &entry[3],&entry[4]);
04306 
04307       for (j=0; j<numread; j++) {
04308         data->normal_modes[(3*k+1) + (i*5+j)*3*data->numatoms] =
04309           entry[j];
04310       }
04311 
04312       /* z */
04313       GET_LINE(buffer, data->file);
04314       numread = sscanf(buffer,"%*s %lf %lf %lf %lf %lf", &entry[0],
04315                        &entry[1], &entry[2], &entry[3],&entry[4]);
04316 
04317       for (j=0; j<numread; j++) {
04318         data->normal_modes[(3*k+2) + (i*5+j)*3*data->numatoms] = 
04319           entry[j];
04320       }
04321     }
04322   }
04323 
04324 
04325   /* Chop unused part of imag_modes array */
04326   data->imag_modes = 
04327     (int *)realloc(data->imag_modes, data->nimag*sizeof(int));
04328 
04329 /*   free(token); */
04330 
04331   data->have_normal_modes = TRUE;
04332   printf("gamessplugin) Successfully scanned normal modes (%d imag.)\n", data->nimag);
04333 
04334   return TRUE;
04335 }
04336 
04337 
04338 
04339 /***********************************************************
04340  *
04341  * Read the cartesian hessian matrix 
04342  * XXX Does not read blocks with less than 6 entries correctly!
04343  *
04344  * *********************************************************/
04345 static int get_cart_hessian(qmdata_t *data)
04346 {
04347   char buffer[BUFSIZ];
04348   int i,j,k;
04349   float entry[6]; 
04350 
04351   buffer[0] = '\0';
04352   memset(entry, 0, sizeof(entry));
04353 
04354   /* at this point we need to rewind the file, since
04355    * in case that there is no internal Hessian stuff the
04356    * previous call to get_int_coords scanned the file
04357    * until EOF */
04358   rewind(data->file);
04359 
04360   if (pass_keyline(data->file,
04361                    "CARTESIAN FORCE CONSTANT MATRIX",
04362                    NULL) != FOUND) {
04363     return FALSE;
04364   }
04365 
04366   /* skip next 5 lines */
04367   eatline(data->file, 5);
04368 
04369 
04370   /* reserve memory for array; 
04371    * NOTE: this is a lower triangular matrix, but for now
04372    * we save it in an square matrix of dim(3Nx3N) to 
04373    * facilitate element access */
04374   data->carthessian = 
04375     (double *)calloc((data->numatoms*3)*(data->numatoms*3),
04376                      sizeof(double));
04377 
04378   
04379   /* make sure memory was allocated properly */
04380   if (data->carthessian == NULL) {
04381     PRINTERR;
04382     return FALSE;
04383   }
04384 
04385 
04386   /* start scanning; the cartesian hessian matrix is a lower
04387    * triangular matrix, organized in rows of 6 */
04388 
04389   /* read blocks with complete rows of 6 */
04390   for (i=0; i<(int)ceil(data->numatoms/2.f); i++) {
04391     for (j=0; j<(data->numatoms*3)-(i*6); j++) {
04392       GET_LINE(buffer, data->file);
04393  
04394       if (j%3==0) {
04395         sscanf(buffer,"%*s %*s %*c %f %f %f %f %f %f",
04396                &entry[0],&entry[1],
04397                &entry[2],&entry[3],&entry[4],&entry[5]);
04398       }
04399       else {
04400         sscanf(buffer,"%*1s %f %f %f %f %f %f",
04401                &entry[0],&entry[1],&entry[2],&entry[3],&entry[4],
04402                &entry[5]);
04403       }
04404 
04405 
04406       /* save entries (lower triangular matrix) in a 
04407        * square matrix */
04408       for (k=0; k<=(j<5 ? j : 5); k++) {
04409         data->carthessian[(j+i*6)*3*data->numatoms + (k+i*6)] =
04410           entry[k];
04411       }
04412     }
04413 
04414     /* skip the four lines separating the data blocks */
04415     eatline(data->file, 4);
04416   }
04417 
04418   printf("gamessplugin) Scanned Hessian in CARTESIAN coordinates\n");
04419 
04420   data->have_cart_hessian = TRUE;
04421 
04422   return TRUE;
04423 }
04424   
04425   
04426   
04427 /***********************************************************
04428  *
04429  * Read the internal coordinates and rewind to the file
04430  * position where we started the search.
04431  *
04432  **********************************************************/
04433 static int get_int_coords(qmdata_t *data) {
04434 
04435   char word[BUFSIZ];
04436   char buffer[BUFSIZ];
04437   long filepos, beginning;
04438   int first, second, third, fourth;
04439   double value;
04440   int n, i = 0, j = 0, k = 0, l = 0;
04441 
04442   /* remember current filepos so we can jump back */
04443   beginning = ftell(data->file);
04444 
04445   if (pass_keyline(data->file, "INTERNAL COORDINATES",
04446                    "1 ELECTRON INTEGRALS") != FOUND) {
04447     printf("gamessplugin) No internal coordinates found.\n");
04448     fseek(data->file, beginning, SEEK_SET);
04449     return FALSE;
04450   }
04451 
04452   /* skip next 5 lines */
04453   eatline(data->file, 5);
04454 
04455   /* remember current filepos so we can jump back */
04456   filepos = ftell(data->file);
04457 
04458   /* scan the next line */
04459   GET_LINE(buffer, data->file);
04460   n = sscanf(buffer,"%*s %s", word); 
04461 
04462   /* read line by line */
04463   while (n!=-1) {
04464     /* start counting the number of internal coordinates */
04465     data->nintcoords++;
04466 
04467     /* count the number of bonds, angles, dihedrals */
04468     if (!strcmp(word,"STRETCH")) {
04469       data->nbonds++;
04470     }
04471     else if (!strcmp(word,"BEND")) {
04472       data->nangles++;
04473     }
04474     else if (!strcmp(word,"TORSION")) {
04475       data->ndiheds++;
04476     }
04477     else if (!strcmp(word,"PLA.BEND")) {
04478       data->nimprops++;
04479     }
04480 
04481     /* scan next line */
04482     GET_LINE(buffer, data->file);
04483     n = sscanf(buffer,"%*s %s", word); 
04484   }
04485 
04486   /* now that we know the number of bonds, angles, etc.
04487    * we can read and store the internal coordinates */
04488   fseek(data->file, filepos, SEEK_SET);
04489 
04490 
04491   /* reserve memory for the arrays storing the internal
04492    * coordinates and their values */
04493   data->bonds = (int *)calloc(2*data->nbonds,sizeof(int));
04494   data->angles = (int *)calloc(3*data->nangles,sizeof(int));
04495   data->dihedrals = (int *)calloc(4*data->ndiheds,sizeof(int));
04496   data->impropers = (int *)calloc(4*data->nimprops,sizeof(int));
04497   data->internal_coordinates = (double *)calloc(data->nintcoords,
04498         sizeof(double));
04499 
04500 
04501   /* check if we have sufficient memory available */
04502   if ( (data->bonds == NULL) || 
04503        (data->angles == NULL) ||
04504        (data->dihedrals == NULL) || 
04505        (data->internal_coordinates == NULL)) 
04506   {
04507     PRINTERR; 
04508     return FALSE;
04509   }
04510 
04511 
04512   /* now start going through the internal coordinates
04513    * and save them in the appropriate arrays; here
04514    * I drop all safety check since we went through
04515    * this part of the file already once and should
04516    * be good */
04517  
04518   /* scan the STRETCHES */
04519   for (i=0; i<data->nbonds; i++) {
04520     GET_LINE(buffer, data->file);
04521     sscanf(buffer,"%*s %*s %d %d %lf", &first, &second, &value);
04522 
04523     *(data->bonds+2*i)   = first;
04524     *(data->bonds+2*i+1) = second;
04525     *(data->internal_coordinates+i) = value;
04526   }
04527 
04528   /* scan the BENDS */
04529   for (j=0; j<data->nangles; j++) {
04530     GET_LINE(buffer, data->file);
04531     sscanf(buffer,"%*s %*s %d %d %d %lf",
04532            &first, &second, &third, &value);
04533 
04534     *(data->angles+3*j)   = first;
04535     *(data->angles+3*j+1) = second;
04536     *(data->angles+3*j+2) = third;
04537     *(data->internal_coordinates+i+j) = value;
04538   }
04539 
04540   /* scan the TORSIONS */
04541   for (k=0; k<data->ndiheds; k++) {
04542     GET_LINE(buffer, data->file);
04543     sscanf(buffer,"%*s %*s %d %d %d %d %lf",
04544            &first, &second, &third, &fourth, &value);
04545 
04546     *(data->dihedrals+4*k)   = first;
04547     *(data->dihedrals+4*k+1) = second;
04548     *(data->dihedrals+4*k+2) = third;
04549     *(data->dihedrals+4*k+3) = fourth;
04550     *(data->internal_coordinates+i+j+k) = value;
04551   }
04552 
04553   /* scan the IMPROPERS */
04554   for (l=0; l<data->nimprops; l++) {
04555     GET_LINE(buffer, data->file);
04556     sscanf(buffer,"%*s %*s %d %d %d %d %lf",
04557            &first, &second, &third, &fourth, &value);
04558 
04559     *(data->impropers+4*l)   = first;
04560     *(data->impropers+4*l+1) = second;
04561     *(data->impropers+4*l+2) = third;
04562     *(data->impropers+4*l+3) = fourth;
04563     *(data->internal_coordinates+i+j+k+l) = value;
04564   }
04565 
04566   /* Since the internal coordinate section can appear
04567    * before or after the symmetry section we have to
04568    * jump back to the beginning of the search. */
04569   fseek(data->file, beginning, SEEK_SET);
04570 
04571   printf("gamessplugin) Scanned %d INTERNAL coordinates \n",
04572          data->nintcoords);
04573   printf("gamessplugin)    %d BONDS \n",data->nbonds);
04574   printf("gamessplugin)    %d ANGLES \n",data->nangles);
04575   printf("gamessplugin)    %d DIHEDRALS \n",data->ndiheds);
04576   printf("gamessplugin)    %d IMPROPERS \n",data->nimprops);
04577 
04578   data->have_internals = TRUE;
04579   return TRUE;
04580 }
04581 
04582 
04583 
04584 /***********************************************************
04585  *
04586  * Read the the Hessian in internal coordinates
04587  *
04588  **********************************************************/
04589 static int get_int_hessian(qmdata_t *data) {
04590   char buffer[BUFSIZ];
04591   double hess[5];
04592   int i = 0, j = 0, k = 0, l = 0;
04593 
04594   memset(hess, 0, sizeof(hess));
04595 
04596   if (pass_keyline(data->file,
04597                    "HESSIAN MATRIX IN INTERNAL COORDINATES",
04598                    "ENERGY GRADIENT") != FOUND) {
04599     return FALSE;
04600   }
04601   if (pass_keyline(data->file,
04602                      "UNITS ARE HARTREE/",
04603                      "ENERGY GRADIENT") != FOUND) {
04604     return FALSE;
04605   }
04606 
04607   eatline(data->file, 3);
04608   
04609   /* reserve memory for inthessian array */
04610   data->inthessian = 
04611     (double *)calloc((data->nintcoords)*(data->nintcoords),
04612                      sizeof(double));
04613 
04614 
04615   /* make sure memory was allocated properly */
04616   if (data->inthessian == NULL) {
04617     PRINTERR;
04618     return FALSE;
04619   }
04620 
04621 
04622   /* start scanning; GAMESS organized the output of the
04623    * internal HESSIAN in rows of 5 */
04624 
04625   /* read blocks with complete rows of 5 */
04626   for (i=0; i<(int)ceil(data->nintcoords/5.f); i++) {
04627     for (j=0; j<data->nintcoords; j++) {
04628       int numread = 0;
04629 
04630       GET_LINE(buffer, data->file);
04631       numread = sscanf(buffer,"%*d %lf %lf %lf %lf %lf", &hess[0], 
04632              &hess[1], &hess[2], &hess[3], &hess[4]);
04633 
04634       /* save entries */
04635       for (k=0; k<numread; k++) { 
04636         data->inthessian[j*data->nintcoords + i*5+k] = hess[k];
04637       }
04638     }
04639 
04640     /* skip the two lines separating the matrix entries 
04641      * and scan next line */
04642     eatline(data->file, 2);
04643 
04644     GET_LINE(buffer, data->file);
04645   }
04646 
04647 #if 0
04648   /* read the remaining block with less then 5 rows
04649    * if present */
04650   remaining_blocks = data->nintcoords%5;
04651   
04652   if (remaining_blocks!=0) {
04653     for (j=0; j<data->nintcoords; j++) {
04654       GET_LINE(buffer, data->file);
04655       sscanf(buffer,"%*d %lf %lf %lf %lf %lf", &hess[0], 
04656              &hess[1], &hess[2], &hess[3], &hess[4]);
04657 
04658       for (k=0; k<remaining_blocks; k++) { 
04659         *(data->inthessian+(j*data->nintcoords)+(i*5)+k) = hess[k];
04660       }
04661     }
04662   }
04663 #endif
04664 
04665   printf("gamessplugin) Scanned Hessian in INTERNAL coordinates\n");
04666 
04667   /* finally, dump the diagonal elements of the hessian into the
04668    * force constant arrays, after converting the units 
04669    * appropriately;
04670    * BONDS are in HARTREE/BOHR**2
04671    * ANGLES,DIHEDRALS,IMPROPERS are in HARTREE/RADIAN**2 */
04672   
04673   /* allocate dynamic arrays */
04674   data->bond_force_const = 
04675     (double *)calloc(data->nbonds, sizeof(double));
04676 
04677   if (data->bond_force_const==NULL) {
04678     PRINTERR;
04679     return FALSE;
04680   }
04681 
04682 
04683   data->angle_force_const =
04684     (double *)calloc(data->nangles, sizeof(double));
04685 
04686   if (data->angle_force_const==NULL) {
04687     PRINTERR;
04688     return FALSE;
04689   }
04690 
04691 
04692   data->dihedral_force_const =
04693     (double *)calloc(data->ndiheds, sizeof(double));
04694 
04695   if (data->dihedral_force_const==NULL) {
04696     PRINTERR;
04697     return FALSE;
04698   }
04699 
04700 
04701   data->improper_force_const =
04702     (double *)calloc(data->nimprops, sizeof(double));
04703 
04704   if (data->improper_force_const==NULL) {
04705     PRINTERR;
04706     return FALSE;
04707   }
04708 
04709   /* scan the bonds */
04710   for (i=0; i<data->nbonds; i++) {
04711     data->bond_force_const[i] = 
04712       data->inthessian[(i*data->nintcoords)+i] * 
04713       HARTREE_TO_KCAL / BOHR_TO_ANGS / BOHR_TO_ANGS;
04714 
04715     printf("%3d (BOND) %2d - %2d : %f\n", i, 
04716            data->bonds[2*i], data->bonds[2*i+1],
04717            data->bond_force_const[i]);
04718   }
04719   
04720   /* scan the angles */
04721   for (j=i; j<i+(data->nangles); j++) {
04722     data->angle_force_const[j-i] = 
04723       data->inthessian[j*data->nintcoords + j] * HARTREE_TO_KCAL;
04724     
04725     printf("%3d (ANGLE) %2d - %2d - %2d : %f\n", j,
04726            data->angles[3*(j-i)], data->angles[3*(j-i)+1], 
04727            data->angles[3*(j-i)+2], 
04728            data->angle_force_const[j-i]);
04729   }
04730 
04731   /* scan the dihedrals */
04732   for (k=j; k<j+(data->ndiheds); k++) {
04733     data->dihedral_force_const[k-j] = 
04734       data->inthessian[k*data->nintcoords + k] * HARTREE_TO_KCAL;
04735     
04736     printf("%3d (DIHEDRAL) %2d - %2d - %2d - %2d : %f \n", k,
04737            data->dihedrals[4*(k-j)  ], data->dihedrals[4*(k-j)+1],
04738            data->dihedrals[4*(k-j)+2], data->dihedrals[4*(k-j)+3],
04739            data->dihedral_force_const[k-j]);
04740   }
04741 
04742   /* scan the impropers */
04743   for (l=k; l<k+(data->nimprops); l++) {
04744     data->improper_force_const[l-k] = 
04745       data->inthessian[l*data->nintcoords + l] * HARTREE_TO_KCAL;
04746     
04747     printf("%3d (IMPROPERS) %2d - %2d - %2d - %2d : %f \n", l,
04748            data->impropers[4*(l-k)  ], data->impropers[4*(l-k)+1],
04749            data->impropers[4*(l-k)+2], data->impropers[4*(l-k)+3],
04750            data->improper_force_const[l-k]);
04751   }
04752 
04753   data->have_int_hessian = TRUE;
04754   return TRUE;
04755 }
04756 
04757 
04758 #if 0
04759 
04760 /************************************************************
04761  *
04762  * this function animates a given normal mode by means of
04763  * generating mod_num_frames frames away from the equilibrium
04764  * structure in a direction given by the hessiane 
04765  *
04766  ************************************************************/
04767 static int animate_normal_mode(qmdata_t *data, int mode)
04768 {
04769   mode_data *animated_mode = data->animated_mode;
04770   float *normal_modes = data->normal_modes;
04771   float scale = animated_mode->mode_scaling;
04772   int i = 0, k = 0; 
04773   int l = 0, m = 0;
04774   int natoms = data->numatoms;
04775   int num_frames = animated_mode->mode_num_frames;
04776 
04777   /* first sweep to max of interval */
04778   for ( k = 0; k < num_frames+1; ++k)
04779   {
04780     for ( i = 0; i < natoms; ++i)
04781     {
04782       *(animated_mode->mode_frames+(k*natoms*3)+(3*i)) = 
04783           (data->atoms+i)->x * (1+( k*scale * 
04784           (*(normal_modes+(mode*natoms*3)+(3*i)))));
04785 
04786       *(animated_mode->mode_frames+(k*natoms*3)+(3*i+1)) = 
04787           (data->atoms+i)->y * (1+( k*scale * 
04788           (*(normal_modes+(mode*natoms*3)+(3*i+1)))));
04789 
04790       *(animated_mode->mode_frames+(k*natoms*3)+(3*i+2)) = 
04791           (data->atoms+i)->z * (1+( k*scale * 
04792           (*(normal_modes+(mode*natoms*3)+(3*i+2)))));
04793     }
04794   }
04795 
04796 
04797   /* second sweep all the way back to min of interval */
04798   for ( l = 0; l < 2*num_frames+1; ++l)
04799   {
04800     for ( i = 0; i < natoms; ++i)
04801     {
04802       *(animated_mode->mode_frames+((l+k)*natoms*3)+(3*i)) = 
04803           (data->atoms+i)->x * (1+((int)(num_frames-l)*scale * 
04804           (*(normal_modes+(mode*natoms*3)+(3*i)))));
04805 
04806       *(animated_mode->mode_frames+((l+k)*natoms*3)+(3*i+1)) = 
04807           (data->atoms+i)->y * (1+((int)(num_frames-l)*scale * 
04808           (*(normal_modes+(mode*natoms*3)+(3*i+1)))));
04809 
04810       *(animated_mode->mode_frames+((l+k)*natoms*3)+(3*i+2)) = 
04811           (data->atoms+i)->z * (1+((int)(num_frames-l)*scale * 
04812           (*(normal_modes+(mode*natoms*3)+(3*i+2)))));
04813     }
04814   }
04815 
04816 
04817   /* third sweep back to the native starting structure */
04818   for ( m = 0; m < num_frames+1; ++m)
04819   {
04820     for ( i = 0; i < natoms; ++i)
04821     {
04822       *(animated_mode->mode_frames+((l+k+m)*natoms*3)+(3*i)) = 
04823           (data->atoms+i)->x * (1+((int)(m-num_frames)*scale * 
04824           (*(normal_modes+(mode*natoms*3)+(3*i)))));
04825 
04826       *(animated_mode->mode_frames+((l+k+m)*natoms*3)+(3*i+1)) = 
04827           (data->atoms+i)->y * (1+((int)(m-num_frames)*scale * 
04828           (*(normal_modes+(mode*natoms*3)+(3*i+1)))));
04829 
04830       *(animated_mode->mode_frames+((l+k+m)*natoms*3)+(3*i+2)) = 
04831           (data->atoms+i)->z * (1+((int)(m-num_frames)*scale * 
04832           (*(normal_modes+(mode*natoms*3)+(3*i+2)))));
04833     }
04834   }
04835 
04836   printf("gamessplugin) Successfully animated mode %d \n", mode);
04837 
04838   return TRUE;
04839 }
04840 #endif
04841 
04842 
04843 
04844 /*************************************************************
04845  *
04846  * plugin registration 
04847  *
04848  **************************************************************/
04849 static molfile_plugin_t plugin;
04850 
04851 VMDPLUGIN_API int VMDPLUGIN_init(void) {
04852   memset(&plugin, 0, sizeof(molfile_plugin_t));
04853   plugin.abiversion = vmdplugin_ABIVERSION;
04854   plugin.type = MOLFILE_PLUGIN_TYPE;
04855   plugin.name = "gamess";
04856   plugin.prettyname = "GAMESS";
04857   plugin.author = "Jan Saam, Markus Dittrich, Johan Strumpfer";
04858   plugin.majorv = 1;
04859   plugin.minorv = 2;
04860   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
04861   plugin.filename_extension = "log";
04862   plugin.open_file_read = open_gamess_read;
04863   plugin.read_structure = read_gamess_structure;
04864   plugin.close_file_read = close_gamess_read;
04865 
04866   plugin.read_qm_metadata = read_gamess_metadata;
04867   plugin.read_qm_rundata  = read_gamess_rundata;
04868 
04869 #if vmdplugin_ABIVERSION > 11
04870   plugin.read_timestep_metadata    = read_timestep_metadata;
04871   plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
04872   plugin.read_timestep = read_timestep;
04873 #endif
04874 
04875   return VMDPLUGIN_SUCCESS;
04876 }
04877 
04878 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
04879   (*cb)(v, (vmdplugin_t *)&plugin);
04880   return VMDPLUGIN_SUCCESS;
04881 }
04882 
04883 VMDPLUGIN_API int VMDPLUGIN_fini(void) {
04884   return VMDPLUGIN_SUCCESS;
04885 }

Generated on Sat Dec 14 03:10:20 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002