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

gaussianplugin.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 /* *******************************************************
00011  *
00012  *          Gaussian Logfile Reader Plugin
00013  *
00014  * This plugin allows VMD to read Gaussian log files.
00015  * The main purpose is to import Wavefunction data.
00016  * It is modeled after the corresponding GAMESS plugin.
00017  *
00018  * ********************************************************/
00019 #include <stdlib.h>
00020 #include <string.h>
00021 #include <ctype.h>
00022 #include <errno.h>
00023 #include <time.h>
00024 #include <math.h>
00025 
00026 #include "gaussianplugin.h"
00027 #include "periodic_table.h"
00028 #include "unit_conversion.h"
00029 
00030 #define THISPLUGIN plugin
00031 #include "vmdconio.h"
00032  
00033 /*
00034  * Error reporting macro for use in DEBUG mode
00035  */
00036 #ifndef GAUSSIAN_DEBUG
00037 #define GAUSSIAN_DEBUG 0
00038 #endif
00039 #define GAUSSIAN_BASIS_DEBUG 1
00040 
00041 #if GAUSSIAN_DEBUG
00042 #define PRINTERR vmdcon_printf(VMDCON_ERROR,                            \
00043                                "\n In file %s, line %d: \n %s \n \n",   \
00044                                __FILE__, __LINE__, strerror(errno))
00045 #else
00046 #define PRINTERR (void)(0)
00047 #endif
00048 
00049 /*
00050  * Error reporting macro for the multiple fgets calls in
00051  * the code
00052  */
00053 #if GAUSSIAN_DEBUG
00054 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE ;   \
00055     else fprintf(stderr,"%s:%d %s",__FILE__, __LINE__, x)
00056 #else
00057 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE
00058 #endif
00059 
00060 /* make sure pointers are NULLed after free(3)ing them. */
00061 #define SAFE_FREE(ptr) free(ptr); ptr=NULL
00062 /* calloc with test of success */
00063 #define SAFE_CALLOC(ptr,type,count)                                 \
00064   ptr = (type *)calloc(count,sizeof(type));                         \
00065   if (ptr == NULL) {                                                \
00066     PRINTERR;                                                       \
00067     return MOLFILE_ERROR;                                           \
00068   }
00069 
00070 #define UNK_SHELL -666
00071 #define SPD_D_SHELL -5
00072 #define SPD_P_SHELL -4
00073 #define SPD_S_SHELL -3
00074 #define SP_S_SHELL  -2
00075 #define SP_P_SHELL  -1
00076 #define S_SHELL 0
00077 #define P_SHELL 1
00078 #define D_SHELL 2
00079 #define F_SHELL 3
00080 #define G_SHELL 4
00081 #define H_SHELL 5
00082 #define I_SHELL 6
00083 
00084 #define SPIN_ALPHA 0
00085 #define SPIN_BETA  1
00086 
00087 #define STATUS_UNKNOWN       -1
00088 #define STATUS_CONVERGED      0
00089 #define STATUS_SCF_NOT_CONV   1
00090 #define STATUS_TOO_MANY_STEPS 2
00091 #define STATUS_BROKEN_OFF     3
00092 
00093 static const char *runtypes[] = { 
00094   "(unknown)", "ENERGY", "OPTIMIZE", "SADPOINT", "HESSIAN", 
00095   "SURFACE", "DYNAMICS", "PROPERTIES", NULL};
00096 
00097 static const char *scftypes[] = { 
00098   "(unknown)", "RHF", "UHF", "ROHF", "GVB", "MCSCF", "FF", NULL };
00099 
00100 
00101 /* ######################################################## */
00102 /* declaration/documentation of internal (static) functions */
00103 /* ######################################################## */
00104 
00107 static int parse_static_data(gaussiandata *, int *);
00108 
00111 static int have_gaussian(gaussiandata *);
00112 
00114 static int get_proc_mem(gaussiandata *);
00115 
00117 static int get_basis_options(gaussiandata *);
00118 
00120 static int get_runtitle(gaussiandata *);
00121 
00123 static int get_input_structure(gaussiandata *);
00125 static int get_coordinates(FILE *file, qm_atom_t *atoms, int numatoms);
00127 static int get_int_coordinates(FILE *file, qm_atom_t *atoms, int numatoms);
00128 
00130 static int get_contrl(gaussiandata *);
00131 
00134 static int get_final_info(gaussiandata *);
00135 
00136 /* in the function get_basis we parse the basis function section to
00137  * determine the number of basis functions and contraction
00138  * coefficients. For Pople/Huzinaga style basis sets these numbers are in
00139  * principle fixed, and could hence be provided by the the plugin
00140  * itself; however, the user might define his own basis/contraction
00141  * coeffients and hence reading them from the input file seem to be
00142  * more general. We can still override it with an "internal" basis
00143  * set */
00144 static int get_basis(gaussiandata *);
00145 
00146 /* this function replaces the basis set data from the log file with
00147  * with the equivalent data read from an internal basis set data base.
00148  * for simplicity we use the same format as gaussian. the basis set
00149  * data is expected to be in $VMDDIR/basis/<basis-set-name>.gbs. */
00150 static int get_internal_basis(gaussiandata *);
00151 
00152 /* convert shell symmetry type from char to int */
00153 static int shellsymm_int(char *symm);
00154 
00155 /* Populate the flat arrays containing the basis set data */
00156 static int fill_basis_arrays(gaussiandata *);
00157 
00158 static int read_first_frame(gaussiandata *);
00159 
00160 /* this subroutine scans the output file for
00161  * the trajectory information */
00162 static int get_traj_frame(gaussiandata *);
00163 
00164 /* returns 1 if the optimization has converged */
00165 static int find_traj_end(gaussiandata *);
00166 
00167 /* this function parses the input file for the final
00168  * wavefunction and stores it in the appropriate arrays; */
00169 static int get_wavefunction(gaussiandata *, qm_timestep_t *);
00170 
00172 static int get_population(gaussiandata *, qm_timestep_t *);
00173 
00174 /* turn fortran double precision 'D' exponents into c parsable 'E's */
00175 static void fix_fortran_exp(char *string) {
00176   while (*string) {
00177     if ( (*string == 'D') || (*string == 'd')) *string='e';
00178     ++string;
00179   }
00180 }
00181 
00182 /* ######################################################## */
00183 /* Functions that are needed by the molfile_plugin          */
00184 /* interface to provide VMD with the parsed data            */
00185 /* ######################################################## */
00186 
00187 
00188 /***************************************************************
00189  *
00190  * Called by VMD to open the Gaussian logfile and get the number
00191  * of atoms.
00192  * We are also reading all the static (i.e. non-trajectory)
00193  * data here since we have to parse a bit to get the atom count
00194  * anyway. These data will then be provided to VMD by
00195  * read_gaussian_metadata() and read_gaussian_rundata().
00196  *
00197  * *************************************************************/
00198 static void *open_gaussian_read(const char *filename, 
00199                   const char *filetype, int *natoms) {
00200 
00201   FILE *fd;
00202   gaussiandata *data;
00203   
00204   /* open the input file */
00205   fd = fopen(filename, "rb");
00206   if (fd == NULL) {
00207     PRINTERR;
00208     return NULL;
00209   }
00210 
00211   /* set up main data structure */
00212   data = (gaussiandata *) calloc(1,sizeof(gaussiandata));
00213   if (data == NULL) return NULL;
00214   
00215   data->runtyp = RUNTYP_UNKNOWN;
00216   data->scftyp = SCFTYP_UNKNOWN;
00217   data->file = fd;
00218   data->file_name = strdup(filename);
00219 
00220   /* check if the file is Gaussian format; 
00221    * if yes parse it, if not exit */
00222   if (have_gaussian(data)==TRUE) {
00223     /* if we're dealing with an unsupported Gaussian
00224      * version, we better quit. so far we tested g98(g94?),g03 */
00225     if ((data->version < 19940000) || (data->version > 20040000)) {
00226       vmdcon_printf(VMDCON_ERROR,
00227                     "gaussianplugin) Gaussian version %s is not "
00228                     "(yet) supported. Bailing out.\n",
00229                     data->version_string);
00230       free(data);
00231       return NULL;
00232     }
00233 
00234     /* get the non-trajectory information from the log file */    
00235     if (parse_static_data(data, natoms) == FALSE) {
00236       free(data);
00237       return NULL;
00238     }
00239   }
00240   else {
00241     free(data);
00242     return NULL;
00243   }
00244 
00245   return data;
00246 }
00247 
00248 
00249 /************************************************************
00250  * 
00251  * Provide VMD with the structure of the molecule, i.e the
00252  * atoms coordinates names, etc.
00253  *
00254  *************************************************************/
00255 static int read_gaussian_structure(void *mydata, int *optflags, 
00256                       molfile_atom_t *atoms) 
00257 {
00258   gaussiandata *data = (gaussiandata *)mydata;
00259   qm_atom_t *cur_atom;
00260   molfile_atom_t *atom;
00261   int i = 0;
00262  
00263   /* optional data from PTE */
00264   *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS;
00265 
00266   if (data->have_mulliken) 
00267     *optflags |= MOLFILE_CHARGE;
00268 
00269   /* all the information I need has already been read in
00270    * via the initial scan and I simply need to copy 
00271    * everything from the temporary arrays into the 
00272    * proper VMD arrays. */
00273 
00274   /* get initial pointer for atom array */
00275   cur_atom = data->initatoms;
00276 
00277   for(i=0; i<data->numatoms; i++) {
00278     atom = atoms+i;
00279     strncpy(atom->name, cur_atom->type, sizeof(atom->name)); 
00280     strncpy(atom->type, get_pte_label(cur_atom->atomicnum), sizeof(atom->type));
00281     strncpy(atom->resname,"QM", sizeof(atom->resname)); 
00282     atom->resid = 1;
00283     atom->chain[0] = '\0';
00284     atom->segid[0] = '\0';
00285     atom->atomicnumber = cur_atom->atomicnum;
00286     atom->radius = get_pte_vdw_radius(cur_atom->atomicnum);
00287     /* XXX; check for isotopes. should be possible to read. */
00288     atom->mass   = get_pte_mass(cur_atom->atomicnum);  
00289     if (data->have_mulliken)
00290       atom->charge = data->qm_timestep->mulliken_charges[i];
00291 
00292     cur_atom++;
00293   }
00294  
00295   return MOLFILE_SUCCESS; 
00296 }
00297 
00298 
00299 /*****************************************************
00300  *
00301  * provide VMD with the sizes of the QM related
00302  * data structure arrays that need to be made
00303  * available
00304  *
00305  *****************************************************/
00306 static int read_gaussian_metadata(void *mydata, 
00307     molfile_qm_metadata_t *gaussian_metadata) {
00308 
00309   gaussiandata *data = (gaussiandata *)mydata;
00310 
00311   gaussian_metadata->ncart = 0;
00312   gaussian_metadata->nimag = 0;
00313   gaussian_metadata->nintcoords = 0;
00314 
00315   /* orbital data */
00316   gaussian_metadata->num_basis_funcs = data->num_basis_funcs;
00317   gaussian_metadata->num_basis_atoms = data->num_basis_atoms;
00318   gaussian_metadata->num_shells      = data->num_shells;
00319   gaussian_metadata->wavef_size      = data->wavef_size;  
00320 
00321 #if vmdplugin_ABIVERSION > 11
00322   gaussian_metadata->have_sysinfo = 1;
00323   
00324   /* charges */
00325   gaussian_metadata->have_esp = 0;
00326   gaussian_metadata->have_carthessian = 0;
00327   gaussian_metadata->have_internals   = 0;
00328   gaussian_metadata->have_normalmodes = FALSE;
00329 #endif
00330 
00331   return MOLFILE_SUCCESS;
00332 }
00333 
00334 
00335 /******************************************************
00336  * 
00337  * Provide VMD with the static (i.e. non-trajectory)
00338  * data. That means we are filling the molfile_plugin
00339  * data structures.
00340  *
00341  ******************************************************/
00342 static int read_gaussian_rundata(void *mydata, molfile_qm_t *qm_data) {
00343 
00344   gaussiandata *data = (gaussiandata *)mydata;
00345 
00346   molfile_qm_basis_t   *basis_data   = &qm_data->basis;
00347   molfile_qm_sysinfo_t *sys_data     = &qm_data->run;
00348 
00349   /* fill in molfile_qm_sysinfo_t */
00350   sys_data->nproc = data->nproc;
00351   sys_data->memory = data->memory; 
00352   sys_data->runtype = data->runtyp;
00353   sys_data->scftype = data->scftyp;
00354   sys_data->totalcharge = data->totalcharge;
00355 /*  sys_data->multiplicity = data->multiplicity; */
00356   sys_data->num_electrons = data->num_electrons;
00357   sys_data->num_occupied_A = data->occ_orbitals_A;
00358   sys_data->num_occupied_B = data->occ_orbitals_B;
00359 
00360   strncpy(sys_data->basis_string, data->basis_string,
00361           sizeof(sys_data->basis_string));
00362   
00363   strncpy(sys_data->runtitle, data->runtitle, sizeof(sys_data->runtitle));
00364   strncpy(sys_data->geometry, data->geometry, sizeof(sys_data->geometry));
00365   strncpy(sys_data->version_string, data->version_string,
00366           sizeof(sys_data->version_string));
00367 
00368 #if vmdplugin_ABIVERSION > 11
00369   /* fill in molfile_qm_basis_t */
00370   if (data->num_basis_funcs) {
00371     memcpy(basis_data->num_shells_per_atom, data->num_shells_per_atom, 
00372            data->num_basis_atoms*sizeof(int));
00373     memcpy(basis_data->num_prim_per_shell, data->num_prim_per_shell, 
00374            data->num_shells*sizeof(int));
00375     memcpy(basis_data->shell_symmetry, data->shell_symmetry, 
00376            data->num_shells*sizeof(int));
00377     memcpy(basis_data->basis, data->basis, 2*data->num_basis_funcs*sizeof(float));
00378     memcpy(basis_data->angular_momentum, data->angular_momentum, 
00379            3*data->wavef_size*sizeof(int));
00380   }
00381 #endif
00382  
00383   return MOLFILE_SUCCESS;
00384 }
00385 
00386 
00387 #if vmdplugin_ABIVERSION > 11
00388 
00389 /***********************************************************
00390  * Provide non-QM metadata for next timestep. 
00391  * Required by the plugin interface.
00392  ***********************************************************/
00393 static int read_timestep_metadata(void *mydata,
00394                                   molfile_timestep_metadata_t *meta) {
00395   meta->count = -1;
00396   meta->has_velocities = 0;
00397 
00398   return MOLFILE_SUCCESS;
00399 }
00400 
00401 /***********************************************************
00402  * Provide QM metadata for next timestep. 
00403  * This actually triggers reading the entire next timestep
00404  * since we have to parse the whole timestep anyway in order
00405  * to get the metadata. So we store the read data locally
00406  * and hand them to VMD when requested by read_timestep().
00407  *
00408  ***********************************************************/
00409 static int read_qm_timestep_metadata(void *mydata,
00410                                     molfile_qm_timestep_metadata_t *meta) {
00411   int i, have = 0;
00412   gaussiandata *data = (gaussiandata *)mydata;
00413 #if GAUSSIAN_DEBUG
00414   vmdcon_printf(VMDCON_INFO, 
00415                 "gaussianplugin) read_qm_timestep_metadata(): %d/%d/%d\n",
00416                 data->num_frames, 
00417                 data->num_frames_read,
00418                 data->num_frames_sent);
00419 #endif
00420 
00421   meta->count = -1; /* Don't know the number of frames yet */
00422   meta->has_gradient = 0;
00423 
00424   if (data->num_frames_read > data->num_frames_sent) {
00425     have = 1;
00426   } else if (data->num_frames_read < data->num_frames) {
00427 #if GAUSSIAN_DEBUG
00428     vmdcon_printf(VMDCON_INFO,
00429                   "gaussianplugin) Probing timestep %d\n",
00430                   data->num_frames_read);
00431 #endif
00432     have = get_traj_frame(data);
00433   }
00434 
00435   if (have) {
00436     /* get a pointer to the current qm timestep */
00437     qm_timestep_t *cur_qm_ts = data->qm_timestep+data->num_frames_sent;
00438 #if GAUSSIAN_DEBUG
00439     vmdcon_printf(VMDCON_INFO,
00440                   "gaussianplugin) Approved timestep %d\n", 
00441                   data->num_frames_sent);
00442 #endif
00443     meta->num_scfiter  = cur_qm_ts->num_scfiter;
00444     for (i=0; (i<MAX_NUM_WAVE && i<cur_qm_ts->numwave); i++) { 
00445 #if GAUSSIAN_DEBUG
00446       vmdcon_printf(VMDCON_INFO,
00447                     "gaussianplugin) num_orbitals_per_wavef[%d/%d]=%d\n",
00448                     i+1, cur_qm_ts->numwave, cur_qm_ts->wave[i].num_orbitals);
00449 #endif
00450       meta->num_orbitals_per_wavef[i] = cur_qm_ts->wave[i].num_orbitals;
00451     }
00452     meta->num_wavef  = cur_qm_ts->numwave;
00453     meta->wavef_size = data->wavef_size;
00454   } else {
00455     meta->num_scfiter  = 0;
00456     meta->num_orbitals_per_wavef[0] = 0;
00457     meta->num_wavef = 0;
00458     meta->wavef_size = 0;
00459 
00460     data->end_of_trajectory = TRUE;
00461   }
00462 
00463   return MOLFILE_SUCCESS;
00464 }
00465 
00466 
00467 /***********************************************************
00468  *
00469  * This function provides the data of the next timestep.
00470  * Here we actually don't read the data from file, that had
00471  * to be done already upon calling read_timestep_metadata().
00472  * Instead we copy the stuff from the local data structure
00473  * into the one's provided by VMD.
00474  *
00475  ***********************************************************/
00476 static int read_timestep(void *mydata, int natoms, 
00477        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
00478                          molfile_qm_timestep_t *qm_ts) 
00479 {
00480   gaussiandata *data = (gaussiandata *)mydata;
00481   qm_atom_t *cur_atom;
00482   int i = 0;
00483   qm_timestep_t *cur_qm_ts;
00484 
00485   if (data->end_of_trajectory == TRUE) return MOLFILE_ERROR;
00486 
00487 #if GAUSSIAN_DEBUG
00488   vmdcon_printf(VMDCON_INFO,
00489                 "gaussianplugin) Sending timestep %d\n", 
00490                 data->num_frames_sent);
00491 #endif
00492 
00493   /* initialize pointer for temporary arrays */
00494   cur_atom = data->initatoms; 
00495   
00496   /* copy the coordinates */
00497   for(i=0; i<natoms; i++) {
00498     ts->coords[3*i  ] = cur_atom->x;
00499     ts->coords[3*i+1] = cur_atom->y;
00500     ts->coords[3*i+2] = cur_atom->z; 
00501     cur_atom++;
00502   }    
00503   
00504   /* get a convenient pointer to the current qm timestep */
00505   cur_qm_ts = data->qm_timestep+data->num_frames_sent;
00506 
00507   /* store the SCF energies */
00508   for (i=0; i<cur_qm_ts->num_scfiter; i++) {
00509     qm_ts->scfenergies[i] = cur_qm_ts->scfenergies[i];
00510   }
00511 
00512   /* store the wave function and orbital energies */
00513   if (cur_qm_ts->wave) {
00514     for (i=0; i<cur_qm_ts->numwave; i++) {
00515       qm_wavefunction_t *wave = &cur_qm_ts->wave[i];
00516       if (wave->wave_coeffs && wave->orb_energies) {
00517         memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs,
00518                wave->num_orbitals*data->wavef_size*sizeof(float));
00519         memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
00520                wave->num_orbitals*sizeof(float));
00521       }
00522     }
00523   }
00524 
00525   if (data->runtyp == RUNTYP_ENERGY || data->runtyp == RUNTYP_HESSIAN) {
00526     /* We have only a single point */
00527     data->end_of_trajectory = TRUE;
00528   }
00529 
00530   data->num_frames_sent++;
00531 
00532   return MOLFILE_SUCCESS;
00533 }
00534 #endif
00535 
00536 
00537 
00542 static void close_gaussian_read(void *mydata) {
00543 
00544   gaussiandata *data = (gaussiandata *)mydata;
00545   int i, j;
00546   fclose(data->file);
00547 
00548   free(data->file_name);
00549   free(data->initatoms);
00550   free(data->basis);
00551   free(data->shell_symmetry);
00552   free(data->num_shells_per_atom);
00553   free(data->num_prim_per_shell);
00554   free(data->mulliken_charges);
00555   free(data->internal_coordinates);
00556   free(data->wavenumbers);
00557   free(data->intensities);
00558   free(data->normal_modes);
00559   free(data->angular_momentum);
00560 
00561   if (data->basis_set) {
00562     for(i=0; i<data->numatoms; i++) {
00563       for (j=0; j<data->basis_set[i].numshells; j++) {
00564         free(data->basis_set[i].shell[j].prim);
00565       }
00566       free(data->basis_set[i].shell);
00567     } 
00568     free(data->basis_set);
00569   }
00570 
00571   for (i=0; i<data->num_frames_read; i++) {
00572     free(data->qm_timestep[i].scfenergies);
00573     free(data->qm_timestep[i].gradient);
00574     free(data->qm_timestep[i].mulliken_charges);
00575     for (j=0; j<data->qm_timestep[i].numwave; j++) {
00576       free(data->qm_timestep[i].wave[j].wave_coeffs);
00577       free(data->qm_timestep[i].wave[j].orb_energies);
00578 /*       free(data->qm_timestep[i].wave[j].occupancies); */
00579     }
00580     free(data->qm_timestep[i].wave);
00581   }
00582   free(data->qm_timestep);
00583   
00584   free(data);
00585 }
00586 
00587 /* ####################################################### */
00588 /*             End of API functions                        */
00589 /* The following functions actually do the file parsing.   */
00590 /* ####################################################### */
00591 
00592 
00593 
00594 /********************************************************
00595  *
00596  * Main gaussian log file parser responsible for static,  
00597  * i.e. non-trajectory information.
00598  *
00599  ********************************************************/
00600 static int parse_static_data(gaussiandata *data, int *natoms) 
00601 {
00602   /* Read # of procs and amount of requested memory */
00603   if (!get_proc_mem(data))        return FALSE;
00604 
00605   /* Read the basis options */
00606   if (!get_basis_options(data))   return FALSE;
00607 
00608   /* Read the route section and try to determine
00609    * the job type. exit if unsupported. */
00610   if (!get_contrl(data))          return FALSE;
00611 
00612   /* Read the run title */
00613   if (!get_runtitle(data))        return FALSE;
00614 
00615   /* Read the input atom definitions and geometry */
00616   if (!get_input_structure(data)) return FALSE;
00617   /* provide VMD with the proper number of atoms */
00618   *natoms = data->numatoms;
00619   /* read first set of coordinates and basis set data */
00620   read_first_frame(data);
00621   
00622   /* */
00623   get_final_info(data);
00624 
00625   vmdcon_printf(VMDCON_INFO, "gaussianplugin) found %d QM data frames.\n", data->num_frames);
00626 #if GAUSSIAN_DEBUG
00627   vmdcon_printf(VMDCON_INFO, "gaussianplugin) num_frames_read = %d\n", data->num_frames_read);
00628   vmdcon_printf(VMDCON_INFO, "gaussianplugin) num_frames_sent = %d\n", data->num_frames_sent);
00629 #endif
00630   return TRUE;
00631 }
00632 
00633 /**********************************************************
00634  *
00635  * this subroutine checks if the provided files is
00636  * actually a Gaussian file;
00637  *
00638  **********************************************************/
00639 static int have_gaussian(gaussiandata *data) 
00640 {
00641   char word[4][MOLFILE_BUFSIZ];
00642   char buffer[BUFSIZ];
00643   char *ptr;
00644   int i = 0;
00645  
00646   buffer[0] = '\0';
00647   for (i=0; i<3; i++) word[i][0] = '\0';
00648 
00649   /* check if the file is Gaussian format 
00650    * Gaussian output typically begins with:
00651    * 'Entering Gaussian System' */
00652   i=0; /* check only the first 100 lines */
00653   do {
00654     GET_LINE(buffer, data->file);
00655     sscanf(buffer,"%s%s%s",word[0],word[1],word[2]);
00656     ++i;
00657   } while( (strcmp(word[0],"Entering") || 
00658             strcmp(word[1],"Gaussian") || 
00659             strcmp(word[2],"System,")) && (i<100) );
00660   if (i>=100) return FALSE;
00661   vmdcon_printf(VMDCON_INFO, "gaussianplugin) Analyzing Gaussian log file: %s\n",data->file_name);
00662   
00663   /* now read on until we find the block of text with encoded version
00664    * number and compile date. */
00665   i=0; /* check only the next 100 lines */
00666   do {
00667     GET_LINE(buffer, data->file);
00668     buffer[20] = '\0';          /* length of the version block varies. */
00669   } while ( (i<100) && strcmp(buffer,
00670       " *******************"));
00671   if (i>=100) return FALSE;
00672 
00673   /* one more line... */
00674   GET_LINE(buffer, data->file);
00675   sscanf(buffer,"%s%s%s%s",word[0],word[1],
00676          word[2],word[3]);
00677 
00678   data->version = atoi(word[1]) * 10000;
00679   if (data->version > 700000) {
00680       data->version += 19000000;
00681   } else {
00682       data->version += 20000000;
00683   }
00684   strcpy(data->version_string,word[2]);
00685 
00686   ptr=strrchr(word[2],'-');
00687   if (ptr != NULL) { 
00688       /* extract revision and patchlevel from G##Rev%.## word.*/
00689       ptr += 7;
00690       i =  (*ptr) - 'A' + 1;
00691       data->version += i*100;
00692       ++ptr; ++ptr;
00693       data->version += atoi(ptr);
00694   }
00695   
00696   vmdcon_printf(VMDCON_INFO, 
00697                 "gaussianplugin) Gaussian version = %s  (Version code: %d)\n",
00698                 data->version_string, data->version);
00699   vmdcon_printf(VMDCON_INFO,
00700                 "gaussianplugin) Compiled on      = %s \n", word[3]);
00701 
00702 
00703   return TRUE;
00704 }
00705 
00706 
00707 /**********************************************************
00708  *
00709  * this subroutine reads the number of procs and the amount
00710  * of memory requested
00711  *
00712  **********************************************************/
00713 static int get_proc_mem(gaussiandata *data) {
00714 
00715   char word[5][MOLFILE_BUFSIZ];
00716   char buffer[BUFSIZ];
00717   int nproc,maxmem;
00718   int i, link;
00719 
00720   buffer[0] = '\0';
00721   for (i=0; i<3; i++) word[i][0] = '\0';
00722 
00723   rewind(data->file);
00724 
00725   /* set some defaults */
00726   nproc = 1;
00727   maxmem = -1;
00728   
00729   do {
00730       GET_LINE(buffer, data->file);
00731       sscanf(buffer,"%s%s%s%*s%s%*s%*s%*s%*s%*s%s",
00732              word[0],word[1],word[2],word[3],word[4]);
00733 
00734       /* max dynamical memory. */ 
00735       if ((strcmp(word[0],"Leave") == 0) &&
00736           (strcmp(word[1],"Link")  == 0)) {
00737         link = atoi(word[2]);
00738         /* gaussian uses real*8 words internally. convert to MByte. */
00739         if (link > 1) maxmem=atoi(word[4])/128/1024;
00740       }
00741 
00742       /* number of SMP cpus */
00743       if ( (strcmp(word[0],"Will") == 0) &&
00744            (strcmp(word[1],"use")  == 0) &&
00745            (strcmp(word[2],"up")  == 0) ) {
00746         nproc = atoi(word[3]);
00747       }
00748 
00749       /* detect if have read too far */ 
00750       if (((strcmp(word[0],"Standard") == 0) ||
00751            (strcmp(word[0],"Z-Matrix") == 0) ||
00752            (strcmp(word[0],"Input") == 0) ) &&
00753           (strcmp(word[1],"orientation:")  == 0)) {
00754         /* can not detect memory used */
00755         maxmem=0;
00756       }
00757   } while (maxmem < 0);
00758 
00759   /* store findings */
00760   data->nproc = nproc;
00761   data->memory = maxmem;
00762   if (maxmem) 
00763     vmdcon_printf(VMDCON_INFO, 
00764                   "gaussianplugin) Gaussian used %2d SMP process(es), "
00765                   "% 6d Mbytes of memory \n", nproc, maxmem);
00766 
00767   return TRUE;
00768 }
00769 
00770 
00771 /**********************************************************
00772  *
00773  * Extract basis set options
00774  *
00775  **********************************************************/
00776 static int get_basis_options(gaussiandata *data) {
00777 
00778   char word[5][MOLFILE_BUFSIZ];
00779   char buffer[BUFSIZ], *ptr;
00780   int i = 0, nfunc, nprim, nume_a, nume_b;
00781 
00782   buffer[0] = '\0';
00783   for (i=0; i<3; i++) word[i][0] = '\0';
00784 
00785   /* to be safe let's rewind the file */
00786   rewind(data->file);
00787 
00788   /* scanning for basis set string */
00789   nume_a=-1;
00790   nfunc=-1;
00791   nprim=-1;
00792   do {
00793     GET_LINE(buffer, data->file);
00794     i=sscanf(buffer,"%s%s%s",word[0],word[1],word[2]);
00795     if (i==3) {
00796 
00797       if ( (strcmp(word[0],"Standard") == 0) &&
00798            (strcmp(word[1],"basis:") == 0) ) {
00799 
00800         ptr = &buffer[0] + strlen(buffer) - 1;
00801         while(*ptr==' ') --ptr;
00802         *ptr='\0';
00803         strncpy(data->gbasis, word[2], 10);
00804         strncpy(data->basis_string, buffer+17, MOLFILE_BUFSIZ);
00805 
00806         /* make sure gbasis is uppercase. */
00807         ptr=data->gbasis;
00808         for (;*ptr;++ptr) *ptr=toupper(*ptr);
00809 
00810       } else if ( (strcmp(word[0],"General") == 0) &&
00811                   (strcmp(word[1],"basis") == 0) ) {
00812 
00813         /* General basis read from cards */
00814 
00815         ptr = &buffer[0] + strlen(buffer) - 1;
00816         while(*ptr==' ') --ptr;
00817         *ptr='\0';
00818 
00819         strncpy(data->gbasis, "GEN", 4);
00820         strncpy(data->basis_string, buffer, MOLFILE_BUFSIZ);
00821       
00822       } else if ( (strcmp(word[0],"AO") == 0) 
00823                   && (strcmp(word[1],"basis") == 0) ) {
00824         /* inline basis definition. count number of atomic basis 
00825          * functions and primitive gaussians */
00826         int numcenter, numgauss, numbasis, numshell;
00827         numcenter=numgauss=numbasis=numshell=0;
00828         GET_LINE(buffer, data->file);
00829         do {
00830           int numprim=0;
00831           ++numcenter;
00832           do {
00833             /* angular momentum, number of primitive gaussians
00834              * and first entry */
00835             GET_LINE(buffer, data->file);
00836             i=sscanf(buffer,"%s%d",word[0], &numprim);
00837             if (i==2) {
00838               ++numshell;
00839               if (strcmp(word[0],"S") == 0) {
00840                 numbasis += 1;    /* simple s-shell*/
00841                 numgauss += numprim;
00842               } else if (strcmp(word[0],"P") == 0) {
00843                 numbasis += 3;    /* simple p-shell */
00844                 numgauss += numprim;
00845               } else if (strcmp(word[0],"SP") == 0) {
00846                 numbasis += 1+3;  /* combined sp-shell, stored individually */
00847                 numgauss += 2*numprim;
00848                 numshell += 1;
00849               } else if (strcmp(word[0],"D") == 0) {
00850                 numbasis += 6;    /* cartesian d-shell. pure will be converted */
00851                 numgauss += numprim;
00852               } else if (strcmp(word[0],"SPD") == 0) {
00853                 numbasis += 1+3+6; /* combined s,p,d shell */
00854                 numgauss += 3*numprim;
00855                 numshell += 2;
00856               } else if (strcmp(word[0],"F") == 0) {
00857                 numbasis += 10;   /* cartesian f-shell. pure will be converted */
00858                 numgauss += numprim;
00859               } else if (strcmp(word[0],"G") == 0) {
00860                 numbasis += 15;   /* cartesian g-shell. pure will be converted */
00861                 numgauss += numprim;
00862               } else {
00863                 vmdcon_printf(VMDCON_ERROR, "gaussianplugin) support for %s-"
00864                               "shells is not yet programmed.\n", word[0]);
00865                 return MOLFILE_ERROR;
00866               }
00867             } else if (i==1) {
00868               if (strcmp(word[0],"****") == 0) {
00869                 break;
00870               }
00871             } else {
00872               return MOLFILE_ERROR;
00873             }
00874             /* skip over lines with primitives */
00875             for (i=0; i < numprim; ++i) 
00876               GET_LINE(buffer, data->file);
00877 
00878 #if GAUSSIAN_DEBUG
00879             vmdcon_printf(VMDCON_INFO, "numcenter:% 4d  numbasis:% 4d  numgauss:% 4d  "
00880                     "numshell:% 4d\n", numcenter, numbasis, numgauss, numshell);
00881 #endif        
00882           } while (strcmp(word[0],"****"));
00883           GET_LINE(buffer, data->file);
00884           i=sscanf(buffer,"%s%s",word[0], word[1]);
00885         } while (i == 2);
00886         data->wavef_size = numbasis;
00887         data->num_shells = numshell;
00888         data->num_basis_funcs = numgauss;
00889         data->num_basis_atoms = numcenter;
00890       } else {
00891         i=sscanf(buffer,"%d%s%s%d%s",&nfunc,word[0],word[1],
00892                  &nprim,word[2]);
00893         if (i==5) {
00894           if ((strncmp(word[0],"basis",5) == 0)     &&
00895               (strncmp(word[1],"functions",9) == 0) &&
00896               (strncmp(word[2],"primitive",9) == 0) ) {
00897             GET_LINE(buffer, data->file);
00898             sscanf(buffer,"%d%s%s%d%s",&nume_a,word[0],word[1],
00899                    &nume_b,word[2]);
00900           }
00901         }
00902       }
00903     }
00904   } while (nume_a < 0);
00905 
00906   vmdcon_printf(VMDCON_INFO, 
00907                 "gaussianplugin) %d shells with %d/%d basis functions, "
00908                 "%d/%d primitive gaussians\n", data->num_shells,
00909                 nfunc, data->wavef_size, nprim, data->num_basis_funcs);
00910   vmdcon_printf(VMDCON_INFO, 
00911                 "gaussianplugin) %d QM atoms with %d alpha electrons, "
00912                 "%d beta electron\n", data->num_basis_atoms, nume_a, nume_b);
00913 
00914   data->num_orbitals   = nfunc;
00915   data->occ_orbitals_A = nume_a;
00916   data->occ_orbitals_B = nume_b;
00917   data->num_electrons  = nume_a + nume_b;
00918   data->multiplicity   = (nume_a+nume_b)/2-(nume_a>nume_b?nume_b:nume_a)+1;
00919   
00920   return TRUE;
00921 }
00922 
00923 
00924 /**********************************************************
00925  *
00926  * Extract the run title line
00927  *
00928  **********************************************************/
00929 static int get_runtitle(gaussiandata *data) {
00930 
00931   char buffer[BUFSIZ];
00932   char *temp;
00933   size_t offs;
00934   
00935   /* look for RUN TITLE section */
00936   do {
00937     GET_LINE(buffer, data->file);
00938   } while (strncmp(buffer," -", 2) );
00939   
00940   GET_LINE(buffer, data->file);
00941   do {
00942     char s; 
00943 
00944     offs=strlen(buffer);
00945     temp=&buffer[offs];
00946 
00947     while (*temp == '\n' || *temp == '\r' || *temp == '\0') --temp;
00948     s = *temp;
00949     fgets(temp,BUFSIZ-offs,data->file);
00950     *temp = s;
00951   } while ( strncmp(temp+1,"--",2) );
00952 
00953   offs=strlen(buffer);
00954   temp=&buffer[offs-1];
00955   while (*temp == '-' || *temp == '\n' || *temp == '\r') --temp;
00956   ++temp;  *temp='\0';
00957   strncpy(data->runtitle,buffer,sizeof(data->runtitle));
00958   
00959   return TRUE;
00960 } 
00961 
00962 
00963 /* Read the input atom definitions and geometry */
00964 static int get_input_structure(gaussiandata *data) {
00965   char buffer[BUFSIZ];
00966   char word[4][MOLFILE_BUFSIZ];
00967   int i, numatoms;
00968 
00969   buffer[0] = '\0';
00970   for (i=0; i<4; i++) word[i][0] = '\0';
00971   
00972   GET_LINE(buffer, data->file);
00973   sscanf(buffer,"%s%s",word[0],word[1]);
00974   
00975   if ( ( (strcmp(word[0],"Symbolic") == 0) &&
00976          (strcmp(word[1],"Z-matrix:") == 0) ) ||
00977        ( (strcmp(word[0],"Redundant") == 0) &&
00978          (strcmp(word[1],"internal") == 0)  ) || 
00979        ( (strcmp(word[0],"Z-Matrix") == 0) &&
00980          (strcmp(word[1],"taken") == 0) ) ) {
00981 
00982     /* skip over line with checkpoint file name */
00983     if ( ( (strcmp(word[0],"Redundant") == 0) &&
00984            (strcmp(word[1],"internal") == 0)  ) || 
00985          ( (strcmp(word[0],"Z-Matrix") == 0) &&
00986            (strcmp(word[1],"taken") == 0) ) ) {
00987       GET_LINE(buffer, data->file);
00988     }
00989       
00990     /* charge and multiplicity */
00991     GET_LINE(buffer, data->file);
00992     sscanf(buffer,"%*s%*s%d%*s%*s%d", &(data->totalcharge),
00993            &(data->multiplicity));
00994 
00995     /* parse coordinates from input */
00996     numatoms=0;
00997     data->initatoms=NULL;
00998     i=1;
00999     do {
01000       char *ptr;
01001       qm_atom_t *atm;
01002       
01003       i=1;
01004       GET_LINE(buffer, data->file);
01005       ptr = strtok(buffer," ,\t\n");
01006       
01007       if (ptr == NULL) break;
01008       if (strcmp(ptr, "Variables:") == 0) break;
01009       if (strcmp(ptr, "Recover") == 0) break;
01010       if (strcmp(ptr, "The") == 0) break;
01011       if (strcmp(ptr, "Leave") == 0) break;
01012       
01013       /* for now we only read the atom label */
01014       data->initatoms=realloc(data->initatoms,(numatoms+1)*sizeof(qm_atom_t));
01015       atm = data->initatoms + numatoms;
01016       strncpy(atm->type, ptr, sizeof(atm->type));
01017       atm->atomicnum=get_pte_idx(ptr);
01018       ++numatoms;
01019     } while ( i >= 0 );
01020     /* TODO */
01021   } else {
01022     vmdcon_printf(VMDCON_ERROR,
01023                   "gaussianplugin) ERROR, cannot parse input coordinates.\n");
01024     return FALSE;
01025   }
01026   
01027   /* store number of atoms in data structure */
01028   data->numatoms = numatoms;
01029   vmdcon_printf(VMDCON_INFO, 
01030                 "gaussianplugin) Atoms: %d   Charge: %d   Multiplicity: %d\n", 
01031                 numatoms, data->totalcharge, data->multiplicity);
01032   return TRUE; 
01033 }
01034 
01035 
01036 /**********************************************************
01037  *
01038  * Read data from the Route section.
01039  *
01040  * XXX: there currently do not seem to be any provisions
01041  * XXX: for multi-step jobs, i.e. jobs that either have
01042  * XXX: multiple run types in one calculations or perform 
01043  * XXX: several calculation steps subsequently...
01044  *
01045  **********************************************************/
01046 static int get_contrl(gaussiandata *data) {
01047 
01048   char buffer[BUFSIZ];
01049   const char *vmdbasis;
01050   char *temp;
01051   size_t offs;
01052 
01053   buffer[0] = '\0';
01054   rewind(data->file);
01055   
01056   /* try to find route section. */
01057   do {
01058     GET_LINE(buffer, data->file);
01059   } while( strncmp(buffer," #", 2) );
01060 
01061   /* append to string in buffer until we reach the end of 
01062    * the route section. this way we'll have the whole info
01063    * in one long string.
01064    * Gaussian writes this out with format (x,A70), so 
01065    * joining lines need some special tricks.  */
01066   do {
01067     char s; 
01068 
01069     offs=strlen(buffer);
01070     temp=&buffer[offs];
01071 
01072     while (*temp == '\n' || *temp == '\r' || *temp == '\0') --temp;
01073     s = *temp;
01074     fgets(temp,BUFSIZ-offs,data->file);
01075     *temp = s;
01076   } while ( strncmp(temp+1,"--",2) );
01077 
01078   offs=strlen(buffer);
01079   temp=&buffer[offs-1];
01080   while (*temp == '-' || *temp == '\n' || *temp == '\r') --temp;
01081 
01082   /* some more magic is required.
01083    * make sure we are zero terminated and have a trailing blank.
01084    * the latter allows to search for keywords with strstr without
01085    * getting false positives on keyword strings that are contained
01086    * in other keywords */
01087   ++temp;  *temp=' '; ++temp; *temp='\0';
01088 
01089   /* convert to upper case */
01090   temp=&buffer[0];
01091   while (*temp++) *temp = toupper(*temp);
01092 
01093   /* now we are ready to look for useful information */
01094 
01095   /* will we have wavefunction output? */
01096   if (strstr(buffer," IOP(6/7=3) ")) {
01097     data->have_wavefunction=TRUE;
01098   } else {
01099     data->have_wavefunction=FALSE;
01100   }
01101 
01102   /* is there a basis set in "input format" */
01103   if (strstr(buffer," GFINPUT ")) {
01104     data->have_basis=TRUE;
01105   } else {
01106     data->have_basis=FALSE;
01107   }
01108 
01109   /* cartesian d-basis functions */
01110   if (strstr(buffer," 6D ")) {
01111     data->have_cart_basis |= 1;
01112   }
01113   /* cartesian f-basis functions */
01114   if (strstr(buffer," 10F ")) {
01115     data->have_cart_basis |= 2;
01116   }
01117   /* pure d-basis functions */
01118   if (strstr(buffer," 5D ")) {
01119     data->have_cart_basis &= ~1;
01120   }
01121   /* pure f-basis functions */
01122   if (strstr(buffer," 7F ")) {
01123     data->have_cart_basis &= ~2;
01124   }
01125 
01126   /* find scf calculation type. XXX: might be safer to detect from output. */
01127   if ((strstr(buffer," ROHF/")) ||
01128       (strstr(buffer," ROHF ")) ||
01129       (strstr(buffer," ROMP"))) {
01130     data->scftyp = SCFTYP_ROHF;
01131   } else if (data->multiplicity != 1) {
01132     data->scftyp = SCFTYP_UHF;
01133   } else {
01134     data->scftyp = SCFTYP_RHF;
01135   }
01136         
01137   /* for semi-empirical, we set the basis to valence-STO-3G. */
01138   if ((strstr(buffer," AM1/")) ||
01139       (strstr(buffer," AM1 ")) ||
01140       (strstr(buffer," PM3/")) ||
01141       (strstr(buffer," PM3 ")) ||
01142       (strstr(buffer," MNDO/")) ||
01143       (strstr(buffer," MNDO "))) {
01144 
01145     vmdbasis = getenv("VMDDEFBASISSET");
01146     if (vmdbasis == NULL) 
01147       vmdbasis = "VSTO-3G";
01148     
01149     if(strlen(data->gbasis) == 0)
01150       strncpy(data->gbasis, vmdbasis, sizeof(data->gbasis));
01151 
01152     if(strlen(data->basis_string) == 0) {
01153       strncpy(data->basis_string, "Internal ", sizeof(data->basis_string));
01154       strncat(data->basis_string, vmdbasis, sizeof(data->basis_string) - 10);
01155 
01156       if(data->have_cart_basis & 1)
01157         strcat(data->basis_string," 6D");
01158       else
01159         strcat(data->basis_string," 5D");
01160       if(data->have_cart_basis & 2)
01161         strcat(data->basis_string," 10F");
01162       else
01163         strcat(data->basis_string," 7F");
01164     }
01165   }
01166 
01167   /* find run type. XXX: might be safer to detect from output. */
01168   data->runtyp = RUNTYP_ENERGY;  /* default */
01169   if ((strstr(buffer," FOPT ")) || 
01170       (strstr(buffer," FOPT=")) ||
01171       (strstr(buffer," FOPT(")) ||
01172       (strstr(buffer," OPT=")) ||
01173       (strstr(buffer," OPT(")) ||
01174       (strstr(buffer," OPT "))) {
01175     data->runtyp = RUNTYP_OPTIMIZE;
01176   } 
01177   if (strstr(buffer," FREQ ")) {
01178     data->runtyp = RUNTYP_HESSIAN;
01179   }
01180   if (strstr(buffer," SCAN ")) {
01181     data->runtyp = RUNTYP_SURFACE;
01182   }
01183         
01184   /* make sure we have some basis set string
01185    * and allow to set it from within VMD. */
01186   vmdbasis = getenv("VMDDEFBASISSET");
01187   if(strlen(data->gbasis) == 0) {
01188     if (vmdbasis == NULL) {
01189       strncpy(data->gbasis, "(unknown)", sizeof(data->gbasis));
01190       strncpy(data->basis_string, "(unknown)", sizeof(data->basis_string));
01191     } else {
01192       strncpy(data->gbasis, vmdbasis, sizeof(data->gbasis));
01193       strncpy(data->basis_string, "Internal ", sizeof(data->basis_string));
01194       strncat(data->basis_string, vmdbasis, sizeof(data->basis_string) - 10);
01195 
01196       if(data->have_cart_basis & 1)
01197         strcat(data->basis_string," 6D");
01198       else
01199         strcat(data->basis_string," 5D");
01200       if(data->have_cart_basis & 2)
01201         strcat(data->basis_string," 10F");
01202       else
01203         strcat(data->basis_string," 7F");
01204     }
01205   }
01206 
01207   /* TODO: add more. e.g. Opt=(ModRedundant), IRC */
01208 
01209   /* print (some of) our findings */
01210   vmdcon_printf(VMDCON_INFO, "gaussianplugin) Run-type: %s, SCF-type: %s\n",
01211                 runtypes[data->runtyp], scftypes[data->scftyp]);
01212   vmdcon_printf(VMDCON_INFO, 
01213                 "gaussianplugin) using %s basis set.\n", data->basis_string);
01214   
01215   return TRUE;
01216 }
01217 
01218 static int read_first_frame(gaussiandata *data) {
01219   data->qm_timestep = NULL;
01220 
01221   /* the angular momentum is populated in get_wavefunction 
01222    * which is called by get_traj_frame(). We have obtained
01223    * the array size wavef_size already from the basis set
01224    * statistics */
01225   vmdcon_printf(VMDCON_INFO, 
01226                 "gaussianplugin) preparing for %d atomic basis functions "
01227                 "per wavefunction.\n", data->wavef_size);
01228   SAFE_CALLOC(data->angular_momentum,int,3*data->wavef_size);
01229 
01230   if (!get_traj_frame(data)) {
01231     return FALSE;
01232   }
01233 
01234   data->num_frames = 1;
01235   return TRUE;
01236 }
01237 
01238 /******************************************************
01239  *
01240  * Reads the info printed after the geometry search
01241  * has finished or whatever analysis was done in a 
01242  * single point run, e.g. ESP charges, Hessian, etc.
01243  * Rewinds to the beginning of the search when done,
01244  * because we read this part at in the initial phase
01245  * and might have to look for additional timesteps
01246  * later.
01247  *
01248  ******************************************************/
01249 static int get_final_info(gaussiandata *data) {
01250   long filepos;
01251   filepos = ftell(data->file);
01252 
01253   if (data->runtyp == RUNTYP_OPTIMIZE || 
01254       data->runtyp == RUNTYP_SADPOINT ||
01255       data->runtyp == RUNTYP_SURFACE) {
01256     /* Try to advance to the end of the geometry
01257      * optimization. If no regular end is found we
01258      * won't find any propertiies to read and return. */
01259     if (!find_traj_end(data)) return FALSE;
01260   }
01261 
01262   if (data->runtyp == RUNTYP_HESSIAN || data->runtyp == RUNTYP_SURFACE) {
01263     vmdcon_printf(VMDCON_WARN, "gaussianplugin) this run type is not fully supported\n");
01264   }
01265     
01266   fseek(data->file, filepos, SEEK_SET);
01267   return TRUE; 
01268 }
01269 
01270 
01271 static int get_coordinates(FILE *file, qm_atom_t *atoms, int numatoms) {
01272   char buffer[BUFSIZ];
01273   int atomicnum;
01274   float x,y,z;
01275   int i,n;
01276 
01277   /* we look for:
01278                               Input orientation:                          
01279  ---------------------------------------------------------------------
01280  Center     Atomic     Atomic              Coordinates (Angstroms)
01281  Number     Number      Type              X           Y           Z
01282  ---------------------------------------------------------------------
01283    */
01284   do {
01285     GET_LINE(buffer, file);
01286   } while ((strstr(buffer,"Input orientation:") == NULL) &&
01287            (strstr(buffer,"Z-Matrix orientation:") == NULL));
01288   GET_LINE(buffer, file);
01289   GET_LINE(buffer, file);
01290   GET_LINE(buffer, file);
01291   GET_LINE(buffer, file);
01292   
01293   for (i=0; i < numatoms; ++i) {
01294     GET_LINE(buffer, file);
01295     n = sscanf(buffer,"%*d%d%*d%f%f%f",&atomicnum,&x,&y,&z);
01296     if (n!=4) return FALSE;
01297     atoms[i].x=x;
01298     atoms[i].y=y;
01299     atoms[i].z=z;
01300   }
01301   return TRUE;
01302 }
01303 
01306 static int get_int_coordinates(FILE *file, qm_atom_t *atoms, int numatoms) {
01307   char buffer[BUFSIZ];
01308   int atomicnum;
01309   float x,y,z;
01310   int i,n;
01311 
01312   /* we look for:
01313                               Standard orientation:                          
01314  ---------------------------------------------------------------------
01315  Center     Atomic     Atomic              Coordinates (Angstroms)
01316  Number     Number      Type              X           Y           Z
01317  ---------------------------------------------------------------------
01318    */
01319   do {
01320     GET_LINE(buffer, file);
01321   } while ((strstr(buffer,"Standard orientation:") == NULL) &&
01322            (strstr(buffer,"Z-Matrix orientation:") == NULL));
01323   GET_LINE(buffer, file);
01324   GET_LINE(buffer, file);
01325   GET_LINE(buffer, file);
01326   GET_LINE(buffer, file);
01327   
01328   for (i=0; i < numatoms; ++i) {
01329     GET_LINE(buffer, file);
01330     n = sscanf(buffer,"%*d%d%*d%f%f%f",&atomicnum,&x,&y,&z);
01331     if (n!=4) return FALSE;
01332     atoms[i].x=x;
01333     atoms[i].y=y;
01334     atoms[i].z=z;
01335   }
01336   return TRUE;
01337 }
01338 
01339 
01340 
01341 /*******************************************************
01342  *
01343  * this function reads in the basis set data 
01344  *
01345  * ******************************************************/
01346 /* typical data looks like this: 
01347  ****
01348    <atomnum> 0
01349    <shellsymm> <#prims> <scalefactor> 0.0
01350    <#prims> lines of <coeff(i)> <exp(i>
01351    <shellsymm> <#prims> <scalefactor> 0.0
01352    <#prims> lines of <coeff(i)> <exp(i>
01353     ...
01354  ****    
01355 
01356   1 0
01357  S   6 1.00       0.000000000000
01358       0.1941330000D+05  0.1851598923D-02
01359       0.2909420000D+04  0.1420619174D-01
01360       0.6613640000D+03  0.6999945928D-01
01361       0.1857590000D+03  0.2400788603D+00
01362       0.5919430000D+02  0.4847617180D+00
01363       0.2003100000D+02  0.3351998050D+00
01364  SP   6 1.00       0.000000000000
01365       0.3394780000D+03 -0.2782170105D-02  0.4564616191D-02
01366       0.8101010000D+02 -0.3604990135D-01  0.3369357188D-01
01367       0.2587800000D+02 -0.1166310044D+00  0.1397548834D+00
01368       0.9452210000D+01  0.9683280364D-01  0.3393617168D+00
01369       0.3665660000D+01  0.6144180231D+00  0.4509206237D+00
01370       0.1467460000D+01  0.4037980152D+00  0.2385858009D+00
01371  */
01372 
01373 int get_basis(gaussiandata *data) {
01374 
01375   char buffer[BUFSIZ];
01376   char word[3][MOLFILE_BUFSIZ];
01377   int i; 
01378 
01379   /* no point in searching through the log file,
01380    * if we cannot have this information */
01381   if (!data->have_basis) return FALSE;
01382 
01383   /* search for characteristic line. but only in the next 1000 lines */
01384   i=0;
01385   if(data->version < 20030000) {            /* g98 */
01386     do {
01387       GET_LINE(buffer, data->file);
01388       sscanf(buffer,"%s%s",word[0],word[1]);
01389       ++i;
01390     } while( (strcmp(word[0],"Basis") || 
01391               strcmp(word[1],"set"))  && (i<1000) );
01392   } else {                                  /* g03 */
01393     do {
01394       GET_LINE(buffer, data->file);
01395       sscanf(buffer,"%s%s",word[0],word[1]);
01396       ++i;
01397     } while( (strcmp(word[0],"AO") || 
01398               strcmp(word[1],"basis"))  && (i<1000) );
01399   }
01400   if (i>=1000) {
01401     /* flag that we have no basis set data */
01402     data->num_shells_per_atom=NULL;
01403     data->have_basis=FALSE;
01404     return MOLFILE_ERROR;
01405   }
01406   
01407   /* Allocate space for the basis for all atoms */
01408   /* When the molecule is symmetric the actual number atoms with
01409    * a basis set could be smaller */
01410   SAFE_CALLOC(data->basis_set,basis_atom_t,data->num_basis_atoms);
01411 
01412   for (i=0; i< data->num_basis_atoms; ++i) {
01413     int numshells, numprim;
01414     int numread, ishell;
01415     float scalef;
01416     shell_t *shell;
01417     
01418     /* this line should be '<atomindex> 0'. gaussian allows much more
01419      * flexible input (it is a zero terminated list of indices or labels), 
01420      * but upon writing to the log files there is exactly one basis set 
01421      * specification per atom. */
01422     GET_LINE(buffer, data->file);
01423     if (atoi(buffer) != (i+1) ) {
01424       vmdcon_printf(VMDCON_WARN, 
01425                     "gaussianplugin) basis set atom counter mismatch: %d vs %s\n",
01426                     atoi(buffer), i+1);
01427       SAFE_FREE(data->basis_set);
01428       data->num_shells_per_atom=NULL;
01429       return MOLFILE_ERROR;
01430     }
01431     strncpy(data->basis_set[i].name,data->gbasis,10);
01432     
01433     numshells=0;
01434     shell=NULL;
01435     GET_LINE(buffer, data->file);
01436     while( strncmp(buffer," ****",5) ) {
01437       int n;
01438 
01439       numread=sscanf(buffer,"%s%d%f",word[0],&numprim,&scalef);
01440       if (numread == 3) {
01441 #if GAUSSIAN_DEBUG && GAUSSIAN_BASIS_DEBUG
01442         vmdcon_printf(VMDCON_INFO, "gaussianplugin) atom: %d, element: %s, shell: %d "
01443                       "%s-type shell, %d primitives, scalefactor %f\n", i,
01444                       get_pte_label(data->initatoms[i].atomicnum), numshells+1, 
01445                       word[0], numprim, scalef);
01446 #endif
01447         ;
01448       } else {
01449         vmdcon_printf(VMDCON_WARN, 
01450                       "gaussianplugin) basis set parse error: %s",buffer);
01451         free(data->basis_set);
01452         data->have_basis=FALSE;
01453         data->basis_set=NULL;
01454         data->num_shells_per_atom=NULL;
01455         return FALSE;
01456       }
01457       ishell=numshells;
01458       ++numshells;
01459       shell=realloc(shell,numshells*sizeof(shell_t));
01460       shell[ishell].numprims=numprim;
01461       shell[ishell].symmetry=shellsymm_int(word[0]);
01462       shell[ishell].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01463       if (shell[ishell].symmetry == SP_S_SHELL) {
01464         ++numshells;
01465         shell=realloc(shell,numshells*sizeof(shell_t));
01466         shell[ishell+1].numprims=numprim;
01467         shell[ishell+1].symmetry=SP_P_SHELL;
01468         shell[ishell+1].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01469       }
01470 
01471       for (n=0; n<numprim; ++n) {
01472         GET_LINE(buffer, data->file);
01473         sscanf(buffer,"%s%s%s", word[0],word[1],word[2]);
01474         fix_fortran_exp(word[0]);
01475         shell[ishell].prim[n].exponent=atof(word[0])*scalef*scalef;
01476         fix_fortran_exp(word[1]);
01477         shell[ishell].prim[n].contraction_coeff=atof(word[1]);
01478         if (shell[ishell].symmetry == SP_S_SHELL) {
01479           shell[ishell+1].prim[n].exponent=shell[ishell].prim[n].exponent;
01480           fix_fortran_exp(word[2]);
01481           shell[ishell+1].prim[n].contraction_coeff=atof(word[2]);
01482         }
01483       }
01484       GET_LINE(buffer, data->file);
01485     }
01486 
01487     /* store shells in atom */
01488     data->basis_set[i].numshells = numshells;
01489     data->basis_set[i].shell = shell;
01490   }
01491 
01492   vmdcon_printf(VMDCON_INFO, 
01493                 "gaussianplugin) Parsed %d uncontracted basis functions with %d shells.\n",
01494                 data->num_basis_funcs, data->num_shells);
01495 
01496 #if GAUSSIAN_DEBUG
01497   for (i=0; i<data->num_basis_atoms; i++) {
01498     int j,k,primcount,shellcount;
01499     primcount=0;
01500     shellcount=0;
01501     
01502     printf("%-8s (%10s)\n\n", data->initatoms[i].type, data->basis_set[i].name);
01503     printf("\n");
01504 
01505     for (j=0; j<data->basis_set[i].numshells; j++) {
01506 
01507       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
01508         printf("%6d   %d %7d %22f%22f\n", j,
01509                data->basis_set[i].shell[j].symmetry,
01510                primcount+1,
01511                data->basis_set[i].shell[j].prim[k].exponent,
01512                data->basis_set[i].shell[j].prim[k].contraction_coeff);
01513         primcount++;
01514       }
01515 
01516       printf("\n");
01517       shellcount++;
01518     }
01519   }
01520 #endif
01521 
01522   /* allocate and populate flat arrays needed for molfileplugin */
01523   return fill_basis_arrays(data);
01524 }
01525 
01526 
01527 /*******************************************************
01528  *
01529  * this function reads in the basis set data from 
01530  * <basis>.gbs or $VMDDIR/basis/<basis>.gbs
01531  *
01532  * ******************************************************/
01533 int get_internal_basis(gaussiandata *data) {
01534 
01535   char *vmddir = NULL;
01536   FILE *fp;
01537   char buffer[BUFSIZ];
01538   char word[3][MOLFILE_BUFSIZ];
01539   char filepath[256];
01540   int  i,n; 
01541 
01542   /* no point in adding a basis set if we already have this information */
01543   if (data->have_basis) return TRUE;
01544 
01545   /* try to open basis set database file. */
01546   sprintf(filepath,"%s.gbs",data->gbasis);
01547   fp=fopen(filepath,"rb");
01548   if (fp == NULL) {
01549     vmddir=getenv("VMDDIR");
01550     if (vmddir == NULL) {
01551       vmddir="/usr/local/lib/vmd";
01552     }
01553     sprintf(filepath,"%s/basis/%s.gbs",vmddir,data->gbasis);
01554     fp=fopen(filepath,"rb");
01555   }
01556   
01557   if (fp == NULL) {
01558     vmdcon_printf(VMDCON_ERROR, "gaussianplugin) failed to read basis set "
01559                   "from data base file %s\n", filepath);
01560     data->num_shells_per_atom=NULL;
01561     data->have_basis=FALSE;
01562     return FALSE;
01563   } else {
01564     vmdcon_printf(VMDCON_INFO, "gaussianplugin) reading basis set "
01565                   "from data base file %s\n", filepath);
01566   }
01567   
01568   /* Allocate space for the basis for all atoms */
01569   /* When the molecule is symmetric the actual number atoms with
01570    * a basis set could be smaller */
01571   SAFE_CALLOC(data->basis_set,basis_atom_t,data->num_basis_atoms);
01572 
01573   for (i=0; i < data->num_basis_atoms; ++i) {
01574     int numshells, numprim;
01575     int numread, ishell;
01576     float scalef;
01577     shell_t *shell;
01578     
01579     /* search for the characteristic first line starting with '****'. */
01580     rewind(fp);
01581     do {
01582       fgets(buffer, sizeof(buffer), fp);
01583       sscanf(buffer,"%s%s",word[0],word[1]);
01584     } while(strcmp(word[0],"****"));
01585   
01586     /* search for an entry for the current atom in the format '<name> 0'. */
01587     do {
01588       fgets(buffer, sizeof(buffer), fp);
01589       if (feof(fp)) {
01590         free(data->basis_set);
01591         data->basis_set=NULL;
01592         data->num_shells_per_atom=NULL;
01593         data->have_basis=FALSE;
01594         vmdcon_printf(VMDCON_ERROR, "gaussianplugin) EOF in data base "
01595                       "file %s while looking for element %s.\n", filepath, 
01596                       get_pte_label(data->initatoms[i].atomicnum), buffer);
01597         fclose(fp);
01598         return FALSE;
01599       }
01600       n=sscanf(buffer,"%s%s",word[0],word[1]);
01601     } while ( (n != 2) || strcmp(word[1],"0") ||
01602               (strcmp(word[0],get_pte_label(data->initatoms[i].atomicnum))) );
01603     
01604     strncpy(data->basis_set[i].name,data->gbasis,sizeof(data->gbasis));
01605     
01606     numshells=0;
01607     shell=NULL;
01608 
01609     /* read basis set until end of element */
01610     do {
01611 
01612       fgets(buffer, sizeof(buffer), fp);
01613       if (strstr(buffer,"****")) break;
01614       if (ferror(fp)) {
01615         vmdcon_printf(VMDCON_ERROR, "gaussianplugin) read error in data "
01616                       "base file %s while reading basis of element %s.\n", 
01617                       filepath, get_pte_label(data->initatoms[i].atomicnum));
01618         free(data->basis_set);
01619         data->basis_set=NULL;
01620         return FALSE;
01621       }
01622 
01623       numread=sscanf(buffer,"%s%d%f",word[0],&numprim,&scalef);
01624       if (numread == 3) {
01625 #if GAUSSIAN_DEBUG && GAUSSIAN_BASIS_DEBUG
01626         vmdcon_printf(VMDCON_INFO, "gaussianplugin) atom: %d, element: %s, shell: %d "
01627                       "%s-type shell, %d primitives, scalefactor %f\n", i,
01628                       get_pte_label(data->initatoms[i].atomicnum), numshells+1, 
01629                       word[0], numprim, scalef);
01630 #endif
01631         ;
01632       } else {
01633         vmdcon_printf(VMDCON_WARN, 
01634                       "gaussianplugin) basis set parse error: %s",buffer);
01635         free(data->basis_set);
01636         data->basis_set=NULL;
01637         return FALSE;
01638       }
01639 
01640       ishell=numshells;
01641       ++numshells;
01642       shell=realloc(shell,numshells*sizeof(shell_t));
01643       shell[ishell].numprims=numprim;
01644       shell[ishell].symmetry=shellsymm_int(word[0]);
01645       shell[ishell].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01646       if (shell[ishell].symmetry == SP_S_SHELL) {
01647         ++numshells;
01648         shell=realloc(shell,numshells*sizeof(shell_t));
01649         shell[ishell+1].numprims=numprim;
01650         shell[ishell+1].symmetry=SP_P_SHELL;
01651         shell[ishell+1].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01652       }
01653 
01654       for (n=0; n<numprim; ++n) {
01655         fgets(buffer, sizeof(buffer), fp);
01656         if (ferror(fp)) {
01657           vmdcon_printf(VMDCON_ERROR, "gaussianplugin) read error in data "
01658                         "base file %s while reading basis of element %s.\n", 
01659                         filepath, get_pte_label(data->initatoms[i].atomicnum));
01660           free(data->basis_set);
01661           data->basis_set=NULL;
01662           return FALSE;
01663         }
01664         sscanf(buffer,"%s%s%s", word[0],word[1],word[2]);
01665         fix_fortran_exp(word[0]);
01666         shell[ishell].prim[n].exponent=atof(word[0])*scalef*scalef;
01667         fix_fortran_exp(word[1]);
01668         shell[ishell].prim[n].contraction_coeff=atof(word[1]);
01669         if (shell[ishell].symmetry == SP_S_SHELL) {
01670           shell[ishell+1].prim[n].exponent=shell[ishell].prim[n].exponent;
01671           fix_fortran_exp(word[2]);
01672           shell[ishell+1].prim[n].contraction_coeff=atof(word[2]);
01673         }
01674       }
01675     } while(1);
01676 
01677   
01678     /* store shells in atom */
01679     data->basis_set[i].numshells = numshells;
01680     data->basis_set[i].shell = shell;
01681     
01682     /* store the total number of basis functions */
01683     data->num_shells += numshells;
01684   }
01685 
01686   vmdcon_printf(VMDCON_INFO, 
01687                 "gaussianplugin) Parsed %d uncontracted basis functions. \n",
01688                 data->num_basis_funcs);
01689 
01690   /* allocate and populate flat arrays needed for molfileplugin */
01691   data->have_basis = TRUE;
01692   return fill_basis_arrays(data);
01693 }
01694 
01695 
01696 /**************************************************
01697  *
01698  * Convert shell symmetry type from char to int.
01699  *
01700  ************************************************ */
01701 static int shellsymm_int(char *symm) {
01702   int shell_symmetry;
01703 
01704   switch (toupper(symm[0])) {
01705     case 'S':
01706       if (symm[1] == '\0') {
01707         shell_symmetry = S_SHELL;
01708       } else if (toupper(symm[1]) == 'P') {
01709         if (symm[2] == '\0') {
01710           shell_symmetry = SP_S_SHELL;
01711         } else if (toupper(symm[1]) == 'D') {
01712           shell_symmetry = SPD_S_SHELL;
01713         } else {
01714           shell_symmetry = UNK_SHELL;
01715         } 
01716       } else {
01717         shell_symmetry = UNK_SHELL;
01718       }
01719       break;
01720     case 'L':
01721       shell_symmetry = SP_S_SHELL;
01722       break;
01723     case 'M': 
01724       shell_symmetry = SP_P_SHELL;
01725       break;
01726     case 'P':
01727       shell_symmetry = P_SHELL;
01728       break;
01729     case 'D':
01730       shell_symmetry = D_SHELL;
01731       break;
01732     case 'F':
01733       shell_symmetry = F_SHELL;
01734       break;
01735     case 'G':
01736       shell_symmetry = G_SHELL;
01737       break;
01738     default:
01739       shell_symmetry = UNK_SHELL;
01740       break;
01741   }
01742 
01743   return shell_symmetry;
01744 }
01745 
01746 
01748 static int fill_basis_arrays(gaussiandata *data) {
01749   int i, j, k;
01750   int shellcount, primcount;
01751   float *basis;
01752   int *num_shells_per_atom;
01753   int *num_prim_per_shell;
01754   int *shell_symmetry;
01755 
01756   /* reserve space for pointer to array containing basis
01757    * info, i.e. contraction coeficients and expansion 
01758    * coefficients; need 2 entries per primitive gaussian, i.e.
01759    * exponent and contraction coefficient; also,
01760    * allocate space for the array holding the orbital symmetry
01761    * information per primitive Gaussian.
01762    * Finally, initialize the arrays holding the number of 
01763    * shells per atom and the number of primitives per shell*/
01764   SAFE_CALLOC(basis,float,2*data->num_basis_funcs);
01765   SAFE_CALLOC(num_shells_per_atom,int,data->num_basis_atoms);
01766   SAFE_CALLOC(shell_symmetry,int,data->num_shells);
01767   SAFE_CALLOC(num_prim_per_shell,int,data->num_shells);
01768 
01769   /* place pointers into struct gaussiandata */
01770   data->basis = basis;
01771   data->shell_symmetry = shell_symmetry;
01772   data->num_shells_per_atom = num_shells_per_atom;
01773   data->num_prim_per_shell  = num_prim_per_shell;
01774 
01775   shellcount=primcount=0;
01776   
01777   for(i=0; i<data->num_basis_atoms; i++) {
01778     num_shells_per_atom[i] = data->basis_set[i].numshells;
01779 
01780     for (j=0; j<data->basis_set[i].numshells; j++) {
01781       shell_symmetry[shellcount] = data->basis_set[i].shell[j].symmetry;
01782       num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
01783 
01784       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
01785         basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
01786         basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
01787         primcount++;
01788       }
01789       shellcount++;
01790     }
01791   } 
01792   return TRUE;
01793 }
01794 
01795 
01796 
01799 static int get_traj_frame(gaussiandata *data) {
01800   qm_timestep_t *cur_qm_ts;
01801   int i;
01802   long fpos;
01803   char buffer[BUFSIZ];
01804   char word[MOLFILE_BUFSIZ];  
01805   buffer[0] = '\0';
01806   word[0] = '\0';
01807 
01808 #if GAUSSIAN_DEBUG
01809   vmdcon_printf(VMDCON_INFO, 
01810                 "gaussianplugin) Timestep %d: =======================\n", 
01811                 data->num_frames_read);
01812 #endif
01813 
01814   fpos=ftell(data->file);       /* XXX: */
01815   if (!get_coordinates(data->file, data->initatoms, data->numatoms)) {
01816     vmdcon_printf(VMDCON_WARN, 
01817                   "gaussianplugin) Couldn't find input orientation coordinates"
01818                   " for timestep %d\n", data->num_frames_read);
01819     vmdcon_printf(VMDCON_WARN, 
01820                   "gaussianplugin) Trying internal coordinates instead.\n");
01821 
01822     fseek(data->file, fpos, SEEK_SET); /* XXX: */
01823     if(!get_int_coordinates(data->file, data->initatoms, data->numatoms)) {
01824       vmdcon_printf(VMDCON_ERROR, 
01825                     "gaussianplugin) Couldn't find any coordinates.\n", 
01826                     data->num_frames_read);
01827       return FALSE;
01828     } 
01829   }
01830   
01831   /* allocate more memory for the timestep array */
01832   data->qm_timestep = 
01833     (qm_timestep_t *)realloc(data->qm_timestep, 
01834                              (data->num_frames_read+1)*sizeof(qm_timestep_t));
01835 
01836   /* get a convenient pointer to the current qm timestep and clear it. */
01837   cur_qm_ts = data->qm_timestep+data->num_frames_read;
01838   memset(cur_qm_ts, 0, sizeof(qm_timestep_t));
01839 
01840   /* XXX: this is a guess only. we need to know this for certain */
01841   cur_qm_ts->numwave = 1;
01842   if (data->scftyp == SCFTYP_UHF)
01843     cur_qm_ts->numwave = 2;
01844 
01845   /* Read the basis set. if available */
01846   if (data->num_frames_read == 0) {
01847     get_basis(data);
01848     get_internal_basis(data);
01849   }
01850   
01851   SAFE_CALLOC(cur_qm_ts->wave,qm_wavefunction_t,cur_qm_ts->numwave);
01852   for (i=0; i < cur_qm_ts->numwave; ++i) {
01853     cur_qm_ts->wave[i].idtag = i;
01854     SAFE_CALLOC(cur_qm_ts->wave[i].orb_indices,int,data->wavef_size);
01855     SAFE_CALLOC(cur_qm_ts->wave[i].orb_energies,float,data->wavef_size);
01856     SAFE_CALLOC(cur_qm_ts->wave[i].occupancies,float,data->wavef_size);
01857     SAFE_CALLOC(cur_qm_ts->wave[i].wave_coeffs,float,
01858                 data->wavef_size*data->wavef_size);
01859   }
01860   
01861   /* Try to read wavefunction and orbital energies */
01862   if (get_wavefunction(data, cur_qm_ts) == FALSE) {
01863     vmdcon_printf(VMDCON_WARN, "gaussianplugin) No wavefunction present for timestep %d\n", data->num_frames_read);
01864     /* free storage */
01865     for (i=0; i < cur_qm_ts->numwave; ++i) {
01866       free(cur_qm_ts->wave[i].wave_coeffs);
01867       free(cur_qm_ts->wave[i].orb_energies);
01868       free(cur_qm_ts->wave[i].occupancies);
01869     }
01870     free(cur_qm_ts->wave);
01871     cur_qm_ts->wave=NULL;
01872     cur_qm_ts->numwave=0;
01873   } else {
01874     vmdcon_printf(VMDCON_INFO, "gaussianplugin) Wavefunction found for timestep %d\n", data->num_frames_read);
01875   }
01876 
01877 #if 0
01878   if (get_population(data, cur_qm_ts)) {
01879     vmdcon_printf(VMDCON_INFO, "gaussianplugin) Mulliken charges found\n");
01880   }
01881 #endif
01882 
01883   data->num_frames_read++;
01884 
01885   return TRUE;
01886 }
01887 
01888 
01889 /* Look for the "EQUILIBRIUM GEOMETRY LOCATED" line thereby
01890  * advancing the file pointer so that the final info block
01891  * can be parsed.
01892  * If we don't find this line before the next geometry
01893  * the file pointer will be set back to where the search
01894  * started. */
01895 static int find_traj_end(gaussiandata *data) {
01896   char buffer[BUFSIZ];
01897   long filepos;
01898   filepos = ftell(data->file);
01899 
01900   while (1) {
01901     if (!fgets(buffer, sizeof(buffer), data->file)) break;
01902 
01903     if (strstr(buffer, "Berny optimization.")) {
01904       data->num_frames++;
01905     } else if (strstr(buffer, "Optimization completed.")) {
01906 #if GAUSSIAN_DEBUG
01907       vmdcon_printf(VMDCON_INFO, "gaussianplugin) ==== End of trajectory. ====\n");
01908 #endif
01909       data->opt_status = STATUS_CONVERGED;
01910       return TRUE;
01911     } else if (strstr(buffer, "Optimization stopped.")) {
01912 #if GAUSSIAN_DEBUG
01913       vmdcon_printf(VMDCON_INFO, "gaussianplugin) ==== End of trajectory. ====\n");
01914 #endif
01915       data->opt_status = STATUS_TOO_MANY_STEPS;
01916       return TRUE;
01917     } else if (strstr(buffer, "Convergence failure -- run terminated.")) {
01918       data->opt_status = STATUS_SCF_NOT_CONV;
01919       return FALSE;
01920     } 
01921   }
01922 
01923   /* We didn't find any of the regular key strings,
01924    * the run was most likely broken off and we have an
01925    * incomplete file. */
01926   data->opt_status = STATUS_BROKEN_OFF;
01927 
01928   fseek(data->file, filepos, SEEK_SET);
01929   return FALSE;  
01930 }
01931 
01932 /*********************************************************
01933  *
01934  * this function reads the actual wavefunction, which is
01935  * punched at the end of the log file
01936  *
01937  **********************************************************/
01938 static int get_wavefunction(gaussiandata *data, qm_timestep_t *ts)
01939 {
01940   long filepos;
01941   char buffer[BUFSIZ];
01942   char word[6][MOLFILE_BUFSIZ];
01943   float *orb_erg, *orb_occ;
01944   int orbital_counter;
01945   int i, j, num_values, num_orbs, *orb_idx;
01946   qm_wavefunction_t *wf;
01947   
01948   i=j=num_values=num_orbs=orbital_counter=0;
01949 
01950   buffer[0] = '\0';
01951   for (i=0; i<6; i++) word[i][0] = '\0';
01952 
01953   wf = ts->wave;
01954   if (wf == NULL || ts->numwave == 0) {
01955     PRINTERR;
01956     return FALSE;
01957   }
01958 
01959   /* no point in searching for wavefunction info, if there cannot be any. */
01960   if (!data->have_wavefunction) return FALSE;
01961   
01962   /* default values. to be changed if needed */
01963   wf->type = MOLFILE_WAVE_CANON;
01964   wf->spin = SPIN_ALPHA;
01965   wf->cartesian = 1;
01966   orb_erg = wf->orb_energies;
01967   orb_occ = wf->occupancies;
01968   orb_idx = wf->orb_indices;
01969   
01970   /*
01971    * the following output requires  IOP(6/7=3) in the route section.
01972    *
01973    * Scan for something like this (g98 closed shell):
01974      Molecular Orbital Coefficients
01975                            1         2         3         4         5
01976                            O         O         O         O         O
01977      EIGENVALUES --    -7.61925  -7.61868  -0.88486  -0.63877  -0.55750
01978    1 1   B  1S          0.71623  -0.69260  -0.13195  -0.12199   0.00000
01979    2        2S          0.01884  -0.01744   0.20755   0.19054   0.00000
01980    3        2PX         0.00000   0.00000   0.00000   0.00000   0.24930
01981 
01982    * or this (g03 closed shell):
01983 
01984      Molecular Orbital Coefficients
01985                            1         2         3         4         5
01986                        (SGG)--O  (SGU)--O  (SGG)--O  (PIU)--O  (PIU)--O
01987      EIGENVALUES --    -1.30558  -1.16543  -0.63335  -0.56287  -0.51096
01988    1 1   O  1S          0.68764   0.70168   0.16476   0.00000   0.00000
01989    2        1PX         0.00000   0.00000   0.00000   0.70711   0.00000
01990    3        1PY         0.00000   0.00000   0.00000   0.00000   0.70711
01991 
01992 
01993    * or this (g03 open shell):
01994 
01995    Alpha Molecular Orbital Coefficients
01996                            1         2         3         4         5
01997                        (SGG)--O  (SGU)--O  (SGG)--O  (SGU)--O  (PIU)--O
01998      EIGENVALUES --   -20.82303 -20.82274  -1.55752  -1.30463  -0.78343
01999    1 1   O  1S          0.70309   0.70312  -0.15839  -0.17089   0.00000
02000    2        2S          0.01545   0.01540   0.38998   0.42214   0.00000
02001 
02002 
02003    */
02004 
02005   /* remember position in order to go back if no wave function was found */
02006   filepos = ftell(data->file);
02007 
02008   do {
02009     GET_LINE(buffer, data->file);
02010     /* check if we searched too far. XXX need more tests here. */
02011     if (strstr(buffer,"Input orientation") ) {
02012       fseek(data->file, filepos, SEEK_SET);
02013       return FALSE;
02014     }
02015   } while(!strstr(buffer,"Molecular Orbital Coefficients"));
02016 
02017   while (orbital_counter < data->num_orbitals) {
02018     /* read up to line of orbital energies */
02019     GET_LINE(buffer, data->file);           /* orbital index */
02020     num_orbs = sscanf(buffer,"%s%s%s%s%s",word[0],word[1],word[2],word[3],word[4]);
02021     /* XXX: we need to keep track of these numbers, since the wavefunction may
02022             have only a subset of orbitals, and we want to know where the 
02023             frontier orbitals are located. with the similarity reordering 
02024             in VMD this information will be essential, if somebody compares 
02025             the .log file and the orbital rep. */
02026 
02027     GET_LINE(buffer, data->file); /* occupied or virtual orbital (+ orbital symm) */
02028     sscanf(buffer,"%s%s%s%s%s", word[0],word[1],word[2],word[3],word[4]);
02029     for (i=0; i<num_orbs; i++) {
02030       j=strlen(word[i]);
02031       orb_occ[i] = (word[i][j-1] == 'O') ? 1.0f : 0.0f;
02032     }
02033     
02034     GET_LINE(buffer, data->file); /* eigenvalues */
02035     sscanf(buffer,"%*s%*s%s%s%s%s%s",word[0],word[1],word[2],word[3],word[4]);
02036     for (i=0; i<num_orbs; i++) 
02037       orb_erg[i] = atof(word[i]);
02038       
02039     /* step counters and pointers */
02040     orb_erg += num_orbs;
02041     orb_occ += num_orbs;
02042 
02043     /* now read in the wavefunction */
02044     for (i=0; i<data->num_orbitals; i++) {
02045       int xexp=0, yexp=0, zexp=0;
02046 
02047       /* read in the wavefunction coefficients for up 
02048        * to 5 orbitals at a time line by line */
02049       GET_LINE(buffer, data->file);
02050       num_values = sscanf(buffer+12,"%4s%s%s%s%s%s", 
02051                           word[0], word[1], word[2],
02052                           word[3], word[4], word[5]);
02053 
02054       /* handle magenetic quantum number. in cartesian basis the
02055        * labels are: S, PX, PY, PZ, DXX, DXY, DXZ, DYY, DYZ, DZZ, ...*/
02056       for (j=1; j<strlen(word[0]); j++) {
02057         switch (word[0][j]) {
02058           case 'X':
02059             xexp++;
02060             break;
02061           case 'Y':
02062             yexp++;
02063             break;
02064           case 'Z':
02065             zexp++;
02066             break;
02067             /* if we have pure d/f-functions the nomenclature changes to 
02068              * 'D 0', 'D-1', 'D+1', 'D-2', 'D+2' */
02069           case '+': /* fallthrough */
02070           case '-': /* fallthrough */
02071           case '0': /* fallthrough */
02072           case '1': /* fallthrough */
02073           case '2': /* fallthrough */
02074             if (wf->cartesian) {
02075               wf->cartesian=0;    /* flag it, so we can convert it later. */
02076               vmdcon_printf(VMDCON_ERROR, "gaussianplugin) pure basis function "
02077                             "detected: '%s'. those are not supported yet. bailing out...\n", word[0]);
02078               return FALSE;
02079             }
02080             break;
02081           default:
02082             /* do nothing */
02083             break;
02084         }
02085       }
02086       data->angular_momentum[3*i  ] = xexp;
02087       data->angular_momentum[3*i+1] = yexp;
02088       data->angular_momentum[3*i+2] = zexp;
02089 #if GAUSSIAN_DEBUG && GAUSSIAN_BASIS_DEBUG && 0
02090       vmdcon_printf(VMDCON_INFO,"%s:%d orbitals %d/%d/%d/%d  shell %d/%d/%s: %d %d %d\n", 
02091                     __FILE__, __LINE__, i, data->num_orbitals, num_orbs, orbital_counter, 
02092                     i, data->wavef_size, word[0], xexp, yexp, zexp);
02093 #endif
02094 
02095       /* each orbital has data->wavef_size entries when converted to 
02096        * cartesian spherical harmonics. to ease conversion we use this 
02097        * number as offset when storing them in groups. */
02098       for (j=0 ; j<num_orbs; j++) {
02099         wf->wave_coeffs[(orbital_counter+j)*data->wavef_size+i] = atof(&word[j+1][0]);
02100       }
02101     }
02102     orbital_counter += num_orbs;
02103   }
02104 
02105   /* store the number of orbitals read in */
02106   wf->num_orbitals = orbital_counter;
02107   vmdcon_printf(VMDCON_INFO, 
02108                 "gaussianplugin) Number of orbitals scanned: %d \n",
02109                 orbital_counter);
02110   return TRUE;
02111 }
02112 
02113 /* Read the population analysis section.
02114  * Currently we parse only the Mulliken charges
02115  * but we might want to add support for populations
02116  * and for Lowdin analysis. */
02117 static int get_population(gaussiandata *data, qm_timestep_t *ts) {
02118 #if 0
02119   int i;
02120   char buffer[BUFSIZ];
02121   data->have_mulliken = FALSE;
02122 
02123 
02124   /* Read Mulliken charges if present */
02125   ts->mulliken_charges = 
02126     (double *)calloc(data->num_basis_atoms, sizeof(double));
02127 
02128   if (!ts->mulliken_charges) {
02129     PRINTERR; 
02130     return FALSE;
02131   }
02132   
02133   for (i=0; i<data->num_basis_atoms; i++) {
02134     int n;
02135     float mullpop, mullcharge, lowpop, lowcharge;
02136     GET_LINE(buffer, data->file);
02137     n = sscanf(buffer,"%*i%*s%f%f%f%f",
02138                &mullpop, &mullcharge, &lowpop, &lowcharge);
02139     if (n!=4) return FALSE;
02140     ts->mulliken_charges[i] = mullcharge;
02141   }
02142 
02143   if (i!=data->numatoms) return FALSE;
02144 
02145   data->have_mulliken = TRUE;
02146 
02147 #if GAUSSIAN_DEBUG
02148   vmdcon_printf(VMDCON_INFO, 
02149                 "gaussianplugin) Number of orbitals scanned: %d \n",
02150                 orbital_counter);
02151 #endif
02152 
02153 #endif
02154   return TRUE;
02155 }
02156 
02157 
02158 /*************************************************************
02159  *
02160  * plugin registration 
02161  *
02162  **************************************************************/
02163 
02164 VMDPLUGIN_API int VMDPLUGIN_init(void) {
02165   memset(&plugin, 0, sizeof(molfile_plugin_t));
02166   plugin.abiversion = vmdplugin_ABIVERSION;
02167   plugin.type = MOLFILE_PLUGIN_TYPE;
02168   plugin.name = "gaussian";
02169   plugin.prettyname = "Gaussian Logfile (g94,g98,g03)";
02170   plugin.author = "Axel Kohlmeyer, Markus Dittrich, Jan Saam";
02171   plugin.majorv = 0;
02172   plugin.minorv = 2;
02173   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
02174   plugin.filename_extension = "log";
02175   plugin.open_file_read = open_gaussian_read;
02176   plugin.read_structure = read_gaussian_structure;
02177   plugin.close_file_read = close_gaussian_read;
02178 
02179   plugin.read_qm_metadata = read_gaussian_metadata;
02180   plugin.read_qm_rundata  = read_gaussian_rundata;
02181 
02182 #if vmdplugin_ABIVERSION > 11
02183   plugin.read_timestep_metadata    = read_timestep_metadata;
02184   plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
02185   plugin.read_timestep = read_timestep;
02186 #endif
02187   return VMDPLUGIN_SUCCESS;
02188 }
02189 
02190 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
02191   (*cb)(v, (vmdplugin_t *)&plugin);
02192   return VMDPLUGIN_SUCCESS;
02193 }
02194 
02195 VMDPLUGIN_API int VMDPLUGIN_fini(void) {
02196   return VMDPLUGIN_SUCCESS;
02197 }
02198 
02199 
02200 #ifdef TEST_PLUGIN
02201 
02202 int main(int argc, char *argv[]) {
02203   int numatoms, i, j, optflags;
02204   molfile_atom_t *atoms;
02205   molfile_timestep_t timestep;
02206   molfile_metadata_t metadata;
02207   molfile_timestep_metadata_t ts_metadata;
02208   molfile_qm_timestep_metadata_t qm_ts_metadata;
02209   molfile_qm_metadata_t qm_metadata;
02210   molfile_qm_timestep_t qm_ts;
02211     
02212   void *v;
02213 
02214   while (--argc) {
02215     ++argv;
02216     v = open_gaussian_read(*argv, "log", &numatoms);
02217     if (!v) {
02218       fprintf(stderr, "open_gaussian_read failed for file %s\n", *argv);
02219       return 1;
02220     }
02221     fprintf(stderr, "open_gaussian_read succeeded for file %s\n", *argv);
02222     fprintf(stderr, "number of atoms: %d\n", numatoms);
02223     atoms = (molfile_atom_t *)malloc(sizeof(molfile_atom_t)*numatoms);
02224     read_gaussian_structure(v,&optflags, atoms);
02225 
02226     i = 0;
02227     timestep.coords = (float *)malloc(3*sizeof(float)*numatoms);
02228     
02229     while (1) {
02230       memset(&ts_metadata, 0, sizeof(molfile_timestep_metadata_t));
02231       read_timestep_metadata(v, &ts_metadata);
02232       memset(&qm_metadata, 0, sizeof(molfile_qm_metadata_t));
02233       read_qm_timestep_metadata(v, &qm_ts_metadata);
02234       qm_ts.scfenergies = (double *)malloc(qm_ts_metadata.num_scfiter*sizeof(double));
02235       qm_ts.wave        = (molfile_qm_wavefunction_t *)malloc(qm_ts_metadata.num_wavef
02236                                                               *sizeof(molfile_qm_wavefunction_t));
02237       memset(qm_ts.wave, 0, qm_ts_metadata.num_wavef*sizeof(molfile_qm_wavefunction_t));
02238       for (j=0; (j<10 && j<qm_ts_metadata.num_wavef); j++) {
02239         qm_ts.wave[j].wave_coeffs = (float *) malloc(qm_ts_metadata.num_orbitals_per_wavef[j]
02240                                                      * qm_ts_metadata.wavef_size * sizeof(float));
02241         qm_ts.wave[j].orbital_energies = (float *) malloc(qm_ts_metadata.num_orbitals_per_wavef[j]*sizeof(float));
02242       }
02243       qm_ts.gradient = (float *) malloc(3*numatoms*sizeof(float));
02244       if (read_timestep(v, numatoms, &timestep, &qm_metadata, &qm_ts)) break;
02245       /* do something with data */
02246       /* XXX */
02247 
02248       free(qm_ts.gradient);
02249       for (j=0; (j<10 && j<qm_ts_metadata.num_wavef); j++) {
02250         free(qm_ts.wave[j].wave_coeffs);
02251         free(qm_ts.wave[j].orbital_energies);
02252       }
02253       free(qm_ts.wave);
02254       free(qm_ts.scfenergies);
02255       i++;
02256     }
02257     fprintf(stderr, "ended read_timestep on frame %d\n", i);
02258     close_gaussian_read(v);
02259   }
02260   return 0;
02261 }
02262 #endif

Generated on Wed Nov 11 03:06:24 2020 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002