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

cpmdlogplugin.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  *          CPMD Output File Reader Plugin
00013  *
00014  * This plugin allows VMD to read some CPMD output files.
00015  * So far only PROPERTIES and specially modified 
00016  * MOLECULAR DYNAMICS BO runs.
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 CPMDLOG_DEBUG
00037 #define CPMDLOG_DEBUG 0
00038 #endif
00039 #define CPMDLOG_BASIS_DEBUG 1
00040 
00041 #if CPMDLOG_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 CPMDLOG_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_cpmd(gaussiandata *);
00112 
00113 /* this function replaces the basis set data from the log file with
00114  * with the equivalent data read from an internal basis set data base.
00115  * for simplicity we use the same format as gaussian. the basis set
00116  * data is expected to be in $VMDDIR/basis/<basis-set-name>.gbs. */
00117 static int get_internal_basis(gaussiandata *);
00118 
00119 /* convert shell symmetry type from char to int */
00120 static int shellsymm_int(char *symm);
00121 
00122 /* Populate the flat arrays containing the basis set data */
00123 static int fill_basis_arrays(gaussiandata *);
00124 
00125 static int read_first_frame(gaussiandata *);
00126 
00127 /* this subroutine scans the output file for
00128  * the trajectory information */
00129 static int get_traj_frame(gaussiandata *);
00130 
00131 /* count the number of readable QM timesteps 
00132  * and collect other information about the
00133  * total trajectory. */
00134 static int find_traj_end(gaussiandata *);
00135 
00136 /* this function parses the input file for the final
00137  * wavefunction and stores it in the appropriate arrays; */
00138 static int get_wavefunction(gaussiandata *, qm_timestep_t *, qm_wavefunction_t *);
00139 
00141 static int get_population(gaussiandata *, qm_timestep_t *);
00142 
00143 /* turn fortran double precision 'D' exponents into c parsable 'E's */
00144 static void fix_fortran_exp(char *string) {
00145   while (*string) {
00146     if ( (*string == 'D') || (*string == 'd')) *string='e';
00147     ++string;
00148   }
00149 }
00150 
00151 /* ######################################################## */
00152 /* Functions that are needed by the molfile_plugin          */
00153 /* interface to provide VMD with the parsed data            */
00154 /* ######################################################## */
00155 
00156 
00157 /***************************************************************
00158  *
00159  * Called by VMD to open the CPMD logfile and get the number
00160  * of atoms.
00161  * We are also reading all the static (i.e. non-trajectory)
00162  * data here since we have to parse a bit to get the atom count
00163  * anyway. These data will then be provided to VMD by
00164  * read_cpmdlog_metadata() and read_cpmdlog_rundata().
00165  *
00166  * *************************************************************/
00167 static void *open_cpmdlog_read(const char *filename, 
00168                   const char *filetype, int *natoms) {
00169 
00170   FILE *fd;
00171   gaussiandata *data;
00172   
00173   /* open the input file */
00174   fd = fopen(filename, "rb");
00175   if (fd == NULL) {
00176     PRINTERR;
00177     return NULL;
00178   }
00179 
00180   /* set up main data structure */
00181   data = (gaussiandata *) calloc(1,sizeof(gaussiandata));
00182   if (data == NULL) return NULL;
00183   
00184   data->runtyp = RUNTYP_UNKNOWN;
00185   data->scftyp = SCFTYP_UNKNOWN;
00186   data->file = fd;
00187   data->file_name = strdup(filename);
00188 
00189   /* check if the file is CPMD format; 
00190    * if yes parse it, if not exit */
00191   if (have_cpmd(data)==TRUE) {
00192     /* if we're dealing with an unsupported CPMD
00193      * version, we better quit. so far we can test 3.9.x-3.13.x */
00194     if ((data->version < 30900) || (data->version > 40000)) {
00195       vmdcon_printf(VMDCON_ERROR,
00196                     "cpmdlogplugin) CPMD version %s is not "
00197                     "(yet) supported. Bailing out.\n",
00198                     data->version_string);
00199       free(data);
00200       return NULL;
00201     }
00202 
00203     /* get the non-trajectory information from the log file */    
00204     if (parse_static_data(data, natoms) == FALSE) {
00205       free(data);
00206       return NULL;
00207     }
00208   }
00209   else {
00210     free(data);
00211     return NULL;
00212   }
00213 
00214   return data;
00215 }
00216 
00217 
00218 /************************************************************
00219  * 
00220  * Provide VMD with the structure of the molecule, i.e the
00221  * atoms coordinates names, etc.
00222  *
00223  *************************************************************/
00224 static int read_cpmdlog_structure(void *mydata, int *optflags, 
00225                       molfile_atom_t *atoms) 
00226 {
00227   gaussiandata *data = (gaussiandata *)mydata;
00228   qm_atom_t *cur_atom;
00229   molfile_atom_t *atom;
00230   int i = 0;
00231  
00232   /* optional data from PTE */
00233   *optflags = MOLFILE_ATOMICNUMBER | MOLFILE_MASS | MOLFILE_RADIUS;
00234 
00235   if (data->have_mulliken) 
00236     *optflags |= MOLFILE_CHARGE;
00237 
00238   /* all the information I need has already been read in
00239    * via the initial scan and I simply need to copy 
00240    * everything from the temporary arrays into the 
00241    * proper VMD arrays. */
00242 
00243   /* get initial pointer for atom array */
00244   cur_atom = data->initatoms;
00245 
00246   for(i=0; i<data->numatoms; i++) {
00247     atom = atoms+i;
00248     strncpy(atom->name, cur_atom->type, sizeof(atom->name)); 
00249     strncpy(atom->type, get_pte_label(cur_atom->atomicnum), sizeof(atom->type));
00250     strncpy(atom->resname,"QM", sizeof(atom->resname)); 
00251     atom->resid = 1;
00252     atom->chain[0] = '\0';
00253     atom->segid[0] = '\0';
00254     atom->atomicnumber = cur_atom->atomicnum;
00255     atom->radius = get_pte_vdw_radius(cur_atom->atomicnum);
00256     /* XXX; check on isotopes. should be possible to read. */
00257     atom->mass   = get_pte_mass(cur_atom->atomicnum);  
00258     if (data->have_mulliken)
00259       atom->charge = data->qm_timestep->mulliken_charges[i];
00260 
00261     cur_atom++;
00262   }
00263  
00264   return MOLFILE_SUCCESS; 
00265 }
00266 
00267 
00268 /*****************************************************
00269  *
00270  * provide VMD with the sizes of the QM related
00271  * data structure arrays that need to be made
00272  * available
00273  *
00274  *****************************************************/
00275 static int read_cpmdlog_metadata(void *mydata, 
00276     molfile_qm_metadata_t *gaussian_metadata) {
00277 
00278   gaussiandata *data = (gaussiandata *)mydata;
00279 
00280   gaussian_metadata->ncart = 0;
00281   gaussian_metadata->nimag = 0;
00282   gaussian_metadata->nintcoords = 0;
00283 
00284   /* orbital data */
00285   gaussian_metadata->num_basis_funcs = data->num_basis_funcs;
00286   gaussian_metadata->num_shells      = data->num_shells;
00287   gaussian_metadata->wavef_size      = data->wavef_size;  
00288 
00289 #if vmdplugin_ABIVERSION > 11
00290   /* charges */
00291   gaussian_metadata->have_esp = 0;
00292   gaussian_metadata->have_carthessian = 0;
00293   gaussian_metadata->have_internals   = 0;
00294   gaussian_metadata->have_normalmodes = FALSE;
00295 #endif
00296 
00297   return MOLFILE_SUCCESS;
00298 }
00299 
00300 
00301 /******************************************************
00302  * 
00303  * Provide VMD with the static (i.e. non-trajectory)
00304  * data. That means we are filling the molfile_plugin
00305  * data structures.
00306  *
00307  ******************************************************/
00308 static int read_cpmdlog_rundata(void *mydata, 
00309                                molfile_qm_t *qm_data) {
00310 
00311   gaussiandata *data = (gaussiandata *)mydata;
00312   int i;
00313 
00314   molfile_qm_basis_t   *basis_data   = &qm_data->basis;
00315   molfile_qm_sysinfo_t *sys_data     = &qm_data->run;
00316 
00317   /* fill in molfile_qm_sysinfo_t */
00318   sys_data->nproc = data->nproc;
00319   sys_data->memory = data->memory; 
00320   sys_data->runtype = data->runtyp;
00321   sys_data->scftype = data->scftyp;
00322   sys_data->totalcharge = data->totalcharge;
00323 /*  sys_data->multiplicity = data->multiplicity; */
00324   sys_data->num_electrons = data->num_electrons;
00325   sys_data->num_occupied_A = data->occ_orbitals_A;
00326   sys_data->num_occupied_B = data->occ_orbitals_B;
00327 
00328   strncpy(sys_data->basis_string, data->basis_string,
00329           sizeof(sys_data->basis_string));
00330   
00331   strncpy(sys_data->runtitle, data->runtitle, sizeof(sys_data->runtitle));
00332   strncpy(sys_data->geometry, data->geometry, sizeof(sys_data->geometry));
00333   strncpy(sys_data->version_string, data->version_string,
00334           sizeof(sys_data->version_string));
00335 
00336 #if vmdplugin_ABIVERSION > 11
00337   /* fill in molfile_qm_basis_t */
00338   if (data->num_basis_funcs) {
00339     for (i=0; i<data->numatoms; i++) {
00340       basis_data->num_shells_per_atom[i] = data->num_shells_per_atom[i];
00341     }
00342     
00343     for (i=0; i<data->num_shells; i++) {
00344       basis_data->num_prim_per_shell[i] = data->num_prim_per_shell[i];
00345       basis_data->shell_symmetry[i] = data->shell_symmetry[i];
00346     }
00347     
00348     for (i=0; i<2*data->num_basis_funcs; i++) {
00349       basis_data->basis[i] = data->basis[i];
00350     }
00351 
00352     for (i=0; i<3*data->wavef_size; i++) {
00353       basis_data->angular_momentum[i] = data->angular_momentum[i];
00354     }
00355   }
00356 #endif
00357  
00358   return MOLFILE_SUCCESS;
00359 }
00360 
00361 
00362 #if vmdplugin_ABIVERSION > 11
00363 
00364 /***********************************************************
00365  * Provide non-QM metadata for next timestep. 
00366  * Required by the plugin interface.
00367  ***********************************************************/
00368 static int read_timestep_metadata(void *mydata,
00369                                   molfile_timestep_metadata_t *meta) {
00370   meta->count = -1;
00371   meta->has_velocities = 0;
00372 
00373   return MOLFILE_SUCCESS;
00374 }
00375 
00376 /***********************************************************
00377  * Provide QM metadata for next timestep. 
00378  * This actually triggers reading the entire next timestep
00379  * since we have to parse the whole timestep anyway in order
00380  * to get the metadata. So we store the read data locally
00381  * and hand them to VMD when requested by read_timestep().
00382  *
00383  ***********************************************************/
00384 static int read_qm_timestep_metadata(void *mydata,
00385                                     molfile_qm_timestep_metadata_t *meta) {
00386   int i, have = 0;
00387   gaussiandata *data = (gaussiandata *)mydata;
00388 
00389 #if CPMDLOG_DEBUG
00390   vmdcon_printf(VMDCON_INFO, 
00391                 "cpmdlogplugin) read_qm_timestep_metadata(): %d/%d/%d\n",
00392                 data->num_frames, 
00393                 data->num_frames_read,
00394                 data->num_frames_sent);
00395 #endif
00396 
00397   meta->count = -1; /* Don't know the number of frames yet */
00398   meta->has_gradient = 0;
00399 
00400   if (data->num_frames_read > data->num_frames_sent) {
00401     have = 1;
00402   } else if (data->num_frames_read < data->num_frames) {
00403 #if CPMDLOG_DEBUG
00404     vmdcon_printf(VMDCON_INFO,
00405                   "cpmdlogplugin) Probing timestep %d\n", 
00406                   data->num_frames_read);
00407 #endif
00408     have = get_traj_frame(data);
00409   }
00410 
00411   if (have) {
00412     /* get a pointer to the current qm timestep */
00413     qm_timestep_t *cur_qm_ts = data->qm_timestep+data->num_frames_sent;
00414 #if CPMDLOG_DEBUG
00415     vmdcon_printf(VMDCON_INFO,
00416                   "cpmdlogplugin) Approved timestep %d\n", 
00417                   data->num_frames_sent);
00418 #endif
00419     meta->num_scfiter  = 0;
00420 
00421     for (i=0; (i<MAX_NUM_WAVE && i<cur_qm_ts->numwave); i++) { 
00422 #if CPMDLOG_DEBUG
00423       vmdcon_printf(VMDCON_INFO,
00424                     "cpmdlogplugin) num_orbitals_per_wavef[%d/%d]=%d\n",
00425                     i+1, cur_qm_ts->numwave, cur_qm_ts->wave[i].num_orbitals);
00426 #endif
00427       meta->num_orbitals_per_wavef[i] = cur_qm_ts->wave[i].num_orbitals;
00428     }
00429     meta->wavef_size = data->wavef_size;
00430 
00431   } else {
00432     meta->num_scfiter = 0;
00433     meta->num_orbitals_per_wavef[0] = 0;
00434     meta->num_wavef = 0;
00435     meta->wavef_size = 0;
00436 
00437     data->end_of_trajectory = TRUE;
00438   }
00439 
00440   return MOLFILE_SUCCESS;
00441 }
00442 
00443 
00444 /***********************************************************
00445  *
00446  * This function provides the data of the next timestep.
00447  * Here we actually don't read the data from file, that had
00448  * to be done already upon calling read_timestep_metadata().
00449  * Instead we copy the stuff from the local data structure
00450  * into the one's provided by VMD.
00451  *
00452  ***********************************************************/
00453 static int read_timestep(void *mydata, int natoms, 
00454        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
00455                          molfile_qm_timestep_t *qm_ts) 
00456 {
00457   gaussiandata *data = (gaussiandata *)mydata;
00458   qm_atom_t *cur_atom;
00459   int i = 0;
00460   qm_timestep_t *cur_qm_ts;
00461 
00462   if (data->end_of_trajectory == TRUE) return MOLFILE_ERROR;
00463 
00464 #if CPMDLOG_DEBUG
00465   vmdcon_printf(VMDCON_INFO,
00466                 "cpmdlogplugin) Sending timestep %d\n", 
00467                 data->num_frames_sent);
00468 #endif
00469 
00470   /* initialize pointer for temporary arrays */
00471   cur_atom = data->initatoms; 
00472   
00473   /* copy the coordinates */
00474   for(i=0; i<natoms; i++) {
00475     ts->coords[3*i  ] = cur_atom->x;
00476     ts->coords[3*i+1] = cur_atom->y;
00477     ts->coords[3*i+2] = cur_atom->z; 
00478     cur_atom++;
00479   }    
00480   
00481   /* get a convenient pointer to the current qm timestep */
00482   cur_qm_ts = data->qm_timestep+data->num_frames_sent;
00483 
00484   /* store the SCF energies */
00485   for (i=0; i<cur_qm_ts->num_scfiter; i++) {
00486     qm_ts->scfenergies[i] = cur_qm_ts->scfenergies[i];
00487   }
00488 
00489   /* store the wave function and orbital energies */
00490   if (cur_qm_ts->wave) {
00491     for (i=0; i<cur_qm_ts->numwave; i++) {
00492       qm_wavefunction_t *wave = &cur_qm_ts->wave[i];
00493       if (wave->wave_coeffs && wave->orb_energies) {
00494         memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs,
00495                wave->num_orbitals*data->wavef_size*sizeof(float));
00496         memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
00497                wave->num_orbitals*sizeof(float));
00498       }
00499     }
00500   }
00501 
00502   if (data->runtyp == RUNTYP_ENERGY || data->runtyp == RUNTYP_HESSIAN) {
00503     /* We have only a single point */
00504     data->end_of_trajectory = TRUE;
00505   }
00506 
00507   data->num_frames_sent++;
00508 
00509   return MOLFILE_SUCCESS;
00510 }
00511 #endif
00512 
00513 
00514 
00519 static void close_cpmdlog_read(void *mydata) {
00520 
00521   gaussiandata *data = (gaussiandata *)mydata;
00522   int i, j;
00523   fclose(data->file);
00524 
00525   free(data->file_name);
00526   free(data->initatoms);
00527   free(data->basis);
00528   free(data->shell_symmetry);
00529   free(data->num_shells_per_atom);
00530   free(data->num_prim_per_shell);
00531   free(data->mulliken_charges);
00532   free(data->internal_coordinates);
00533   free(data->wavenumbers);
00534   free(data->intensities);
00535   free(data->normal_modes);
00536   free(data->angular_momentum);
00537 
00538   if (data->basis_set) {
00539     for(i=0; i<data->numatoms; i++) {
00540       for (j=0; j<data->basis_set[i].numshells; j++) {
00541         free(data->basis_set[i].shell[j].prim);
00542       }
00543       free(data->basis_set[i].shell);
00544     } 
00545     free(data->basis_set);
00546   }
00547 
00548   for (i=0; i<data->num_frames_read; i++) {
00549     free(data->qm_timestep[i].scfenergies);
00550     free(data->qm_timestep[i].gradient);
00551     free(data->qm_timestep[i].mulliken_charges);
00552     for (j=0; j<data->qm_timestep[i].numwave; j++) {
00553       free(data->qm_timestep[i].wave[j].wave_coeffs);
00554       free(data->qm_timestep[i].wave[j].orb_energies);
00555 /*       free(data->qm_timestep[i].wave[j].occupancies); */
00556     }
00557     free(data->qm_timestep[i].wave);
00558   }
00559   free(data->qm_timestep);
00560   
00561   free(data);
00562 }
00563 
00564 /* ####################################################### */
00565 /*             End of API functions                        */
00566 /* The following functions actually do the file parsing.   */
00567 /* ####################################################### */
00568 
00570 static int find_traj_end(gaussiandata *data) {
00571   char buffer[BUFSIZ];
00572   long filepos;
00573   filepos = ftell(data->file);
00574 
00575   while (1) {
00576     if (!fgets(buffer, sizeof(buffer), data->file)) break;
00577 
00578     if (strstr(buffer, "PROJECTION COORDINATES")) {
00579       data->num_frames++;
00580     } 
00581   }
00582   data->opt_status = STATUS_UNKNOWN;
00583 
00584   fseek(data->file, filepos, SEEK_SET);
00585   return FALSE;  
00586 }
00587 
00588 
00589 static int get_final_info(gaussiandata *data) {
00590   long filepos;
00591   filepos = ftell(data->file);
00592 
00593   if (data->runtyp == RUNTYP_OPTIMIZE || 
00594       data->runtyp == RUNTYP_DYNAMICS) {
00595     /* Try to advance to the end of the geometry
00596      * optimization or MD. If no regular end is found we
00597      * won't find any propertiies to read and return. */
00598     if (!find_traj_end(data)) return FALSE;
00599   }
00600 
00601 #if 0
00602   if (get_esp_charges(data)) {
00603     vmdcon_printf(VMDCON_INFO, "gaussianplugin) ESP charges found!\n");
00604   }
00605 #endif
00606 
00607   fseek(data->file, filepos, SEEK_SET);
00608   return TRUE; 
00609 }
00610 
00611 
00612 
00613 /********************************************************
00614  *
00615  * Main gaussian log file parser responsible for static,  
00616  * i.e. non-trajectory information.
00617  *
00618  ********************************************************/
00619 static int parse_static_data(gaussiandata *data, int *natoms) 
00620 {
00621   char buffer[BUFSIZ];
00622   char word[4][MOLFILE_BUFSIZ];
00623   char *vmdbasis;
00624   int  i,n;
00625   int  numatoms, numstates, numelectrons, totalcharge;
00626 
00627   buffer[0] = '\0';
00628   
00629   /* set some defaults */
00630   data->scftyp = SCFTYP_RHF;
00631   data->runtyp = RUNTYP_UNKNOWN;
00632   numelectrons = 0;
00633   data->totalcharge = 0;
00634   data->multiplicity = 1;
00635   data->have_basis=FALSE;       
00636   /* CPMD never outputs basis set info, so we use 
00637    * VSTO-6G unless overridden by environment */
00638   vmdbasis = getenv("VMDDEFBASISSET");
00639   if (vmdbasis == NULL) 
00640     vmdbasis = "VSTO-6G";
00641 
00642   strncpy(data->gbasis, vmdbasis, sizeof(data->gbasis));
00643   strncpy(data->basis_string, "Internal ", sizeof(data->basis_string));
00644   strncat(data->basis_string, vmdbasis, sizeof(data->basis_string) - 10);
00645 
00646 
00647   /* try to find job type parameters within the next 100 lines.*/
00648   for (i=0; i<100; ++i) {
00649     GET_LINE(buffer, data->file);
00650     sscanf(buffer,"%s%s%s",word[0],word[1],word[2]);
00651     if ( (strcmp(word[0],"CALCULATE" ) == 0 &&
00652           strcmp(word[1],"SOME"      ) == 0 && 
00653           strcmp(word[2],"PROPERTIES") == 0 ) ) {
00654       data->runtyp = RUNTYP_PROPERTIES;
00655       break;
00656     } else if ( (strcmp(word[0],"GEOMETRY"    ) == 0 &&
00657                  strcmp(word[1],"OPTIMIZATION") == 0 ) ) {
00658       data->runtyp = RUNTYP_OPTIMIZE;
00659       break;
00660     } else if ( (strcmp(word[1],"MOLECULAR"    ) == 0 &&
00661                  strcmp(word[2],"DYNAMICS") == 0 ) ) {
00662       data->runtyp = RUNTYP_DYNAMICS;
00663       break;
00664       
00665     }
00666   }
00667 
00668   /* XXX: add support for other types later? */
00669   if ((data->runtyp != RUNTYP_PROPERTIES) &&
00670       (data->runtyp != RUNTYP_DYNAMICS) ) {
00671     vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) Run-type: %s is currently "
00672                   "not supported by this plugin.\n", runtypes[data->runtyp]);
00673     return FALSE;
00674   }
00675   
00676   if ((data->runtyp != RUNTYP_PROPERTIES) && (data->version < 31303)) {
00677     vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) Run-type: %s is currently "
00678                   "not supported for outputs of this CPMD version.\n", 
00679                   runtypes[data->runtyp]);
00680     return FALSE;
00681   }
00682   
00683   /* scavange for more setup information */
00684   do {
00685 
00686     GET_LINE(buffer, data->file);
00687     n = sscanf(buffer,"%s%s%s%s",word[0],word[1],word[2],word[3]);
00688 
00689     /* empty line */
00690     if (n < 0) continue;
00691 
00692     /* atom types and initial coordinates */
00693     if ( (strstr(word[0],"************") != 0 &&
00694           strcmp(word[1],"ATOMS"       ) == 0 && 
00695           strstr(word[2],"************") != 0 ) ) {
00696       /* NR TYPE X(bohr) Y(bohr) Z(bohr) MBL */
00697       GET_LINE(buffer, data->file);
00698 
00699       numatoms=0;
00700       data->initatoms=NULL;
00701       while (1) {
00702         qm_atom_t *atm;
00703       
00704         GET_LINE(buffer, data->file);
00705         /* end of ATOMS block */
00706         if (strstr(buffer, "*************************") != NULL)
00707           break;
00708         
00709         data->initatoms=realloc(data->initatoms,(numatoms+1)*sizeof(qm_atom_t));
00710         atm = data->initatoms + numatoms;
00711         
00712         n=sscanf(buffer,"%*d%s%g%g%g", atm->type, &atm->x, &atm->y, &atm->z);
00713         if (n != 4) {
00714           free(data->initatoms);
00715           data->initatoms=NULL;
00716           vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) Failed to parse initial" 
00717                         " coordinates. Stopping.\n");
00718           return FALSE;
00719         }
00720         atm->atomicnum=get_pte_idx(atm->type);
00721         atm->x *= BOHR_TO_ANGS;
00722         atm->y *= BOHR_TO_ANGS;
00723         atm->z *= BOHR_TO_ANGS;
00724         ++numatoms;
00725       };
00726       data->numatoms = numatoms;
00727       *natoms = numatoms;
00728       
00729     } else if ( (strcmp(word[0],"NUMBER" ) == 0 &&
00730                  strcmp(word[2],"STATES:") == 0 ) ) {
00731       numstates=atoi(word[3]);
00732 
00733     } else if ( (strcmp(word[0],"NUMBER" ) == 0 &&
00734                  strcmp(word[2],"ELECTRONS:") == 0 ) ) {
00735       numelectrons=(int) (atof(word[3]) + 0.5); /* XXX */
00736       data->num_electrons = numelectrons;
00737       
00738     } else if (strcmp(word[0],"CHARGE:" ) == 0 ) {
00739       totalcharge =(int) (atof(word[1]) + 0.5); /* XXX */
00740       data->totalcharge = totalcharge;
00741       
00742     } else if (strcmp(word[0],"OCCUPATION" ) == 0 ) {
00743       ; /* XXX. */
00744       
00745     } else if ( (strcmp(word[0],"CELL" ) == 0 &&
00746                  strcmp(word[1],"DIMENSION:") == 0 ) ) {
00747       n=sscanf(buffer,"%*s%*s%g%g%g%g%g%g",data->initcell,data->initcell+1, 
00748                data->initcell+2,data->initcell+3,data->initcell+4,data->initcell+5);
00749       
00750     /*  
00751         PROJECT WAVEFUNCTION ON ATOMIC ORBITALS EVERY           10 STEPS
00752     */
00753     } else if ( (strcmp(word[0],"PROJECT" ) == 0 &&
00754                  strcmp(word[1],"WAVEFUNCTION") == 0 ) ) {
00755         data->have_wavefunction=1;
00756     /*
00757       DIPOLE MOMENT CALCULATION 
00758       STORE DIPOLE MOMENTS EVERY                          10 STEPS
00759       WANNIER FUNCTION DYNAMICS 
00760     */
00761     } else if ( (strcmp(word[0],"WANNIER" ) == 0 &&
00762                  strcmp(word[1],"FUNCTION") == 0 ) ) {
00763         if (data->have_wavefunction) {
00764             data->have_wavefunction=2;
00765         }
00766     }
00767   } while( strcmp(word[0],"INITIALIZATION") || 
00768            strcmp(word[1],"TIME:") );
00769 
00770   if (numstates >= numelectrons) 
00771     data->scftyp = SCFTYP_UHF;
00772 
00773   switch (data->scftyp) {
00774     case SCFTYP_RHF:
00775       data->occ_orbitals_A = numstates;
00776       data->occ_orbitals_B = 0;
00777       break;
00778     case SCFTYP_UHF:            /* XXX: this is most likely wrong. check! */
00779       data->occ_orbitals_A = numstates/2;
00780       data->occ_orbitals_B = numstates/2;
00781       break;
00782     default:
00783       break;
00784   }
00785   
00786   vmdcon_printf(VMDCON_INFO, 
00787                 "cpmdlogplugin) Atoms: %d   Charge: %d   Multiplicity: %d\n", 
00788                 data->numatoms, data->totalcharge, data->multiplicity);
00789 
00790   vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) Run-type: %s, SCF-type: %s\n",
00791                 runtypes[data->runtyp], scftypes[data->scftyp]);
00792   vmdcon_printf(VMDCON_INFO, 
00793                 "cpmdlogplugin) using %s basis set.\n", data->basis_string);
00794 
00795   read_first_frame(data);
00796 
00797   get_final_info(data);
00798   
00799   vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) found %d QM data frames.\n", data->num_frames);
00800 #if CPMDLOG_DEBUG
00801   vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) num_frames_read = %d\n", 
00802                 data->num_frames_read);
00803   vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) num_frames_sent = %d\n", 
00804                 data->num_frames_sent);
00805 #endif
00806   return TRUE;
00807 }
00808 
00809 /**********************************************************
00810  *
00811  * this subroutine checks if the provided files is
00812  * actually a CPMD file and gathers its version code.
00813  *
00814  **********************************************************/
00815 static int have_cpmd(gaussiandata *data) 
00816 {
00817   char word[4][MOLFILE_BUFSIZ];
00818   char buffer[BUFSIZ];
00819   char *ptr;
00820   int i = 0;
00821  
00822   buffer[0] = '\0';
00823   for (i=0; i<3; i++) word[i][0] = '\0';
00824 
00825   /* check if the file is CPMD format 
00826    * CPMD output typically begins with:
00827    *  'PROGRAM CPMD STARTED'
00828    */
00829   i=0; /* check only the first 100 lines */
00830   do {
00831     GET_LINE(buffer, data->file);
00832     sscanf(buffer,"%s%s%s",word[0],word[1],word[2]);
00833     ++i;
00834   } while( (strcmp(word[0],"PROGRAM") || 
00835             strcmp(word[1],"CPMD") || 
00836             strcmp(word[2],"STARTED")) && (i<100) );
00837   if (i>=100) return FALSE;
00838   vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) Analyzing CPMD log file: %s\n",data->file_name);
00839   
00840   /* now read on until we find the block of text with encoded version
00841    * number and compile date. */
00842   i=0; /* check only the next 100 lines */
00843   do {
00844     GET_LINE(buffer, data->file);
00845     sscanf(buffer,"%s%s",word[0],word[1]);
00846   } while ( (i<100) &&  (strcmp(word[0],"VERSION")) );
00847   if (i>=100) return FALSE;
00848 
00849   strcpy(data->version_string,word[1]);
00850 
00851   /* now split version number strings */
00852   ptr=strtok(word[1],"._");
00853   data->version = 10000*atoi(ptr);
00854   ptr=strtok(NULL,"._");
00855   data->version += 100*atoi(ptr);
00856   ptr=strtok(NULL,"._");
00857   data->version += atoi(ptr);
00858   
00859   vmdcon_printf(VMDCON_INFO, 
00860                 "cpmdlogplugin) CPMD version = %s  (Version code: %d)\n",
00861                 data->version_string, data->version);
00862   return TRUE;
00863 }
00864 
00866 static int read_first_frame(gaussiandata *data) {
00867   
00868   data->qm_timestep = NULL;
00869 
00870   /* Read the basis set. */
00871   get_internal_basis(data);
00872 
00873   /* the angular momentum is populated in get_wavefunction 
00874    * which is called by get_traj_frame(). We have obtained
00875    * the array size wavef_size already from the basis set
00876    * statistics */
00877   vmdcon_printf(VMDCON_INFO, 
00878                 "cpmdlogplugin) Reserving storage for %d cartesian atomic basis functions\n",
00879                 data->wavef_size);
00880   SAFE_CALLOC(data->angular_momentum,int,3*data->wavef_size);
00881 
00882   if (data->version < 31303) {
00883 
00884     vmdcon_printf(VMDCON_WARN, 
00885                   "cpmdlogplugin) ##################################"
00886                   "####################################\n");
00887     vmdcon_printf(VMDCON_WARN, 
00888                   "cpmdlogplugin) This version of CPMD does not print "
00889                   "the actual coordinates on properties runs.\n");
00890     vmdcon_printf(VMDCON_WARN, 
00891                   "cpmdlogplugin) Using initial coordinates from the "
00892                   "input file instead.\n");
00893     vmdcon_printf(VMDCON_WARN,
00894                   "cpmdlogplugin) These coodinates are most likely "
00895                   "inconsistent with the projection.\n Try the following:\n\n"
00896                   "set cur [atomselect top all]\n"
00897                   "set new [atomselect [mol new GEOMETRY.xyz] all]\n"
00898                   "$cur set {x y z} [$new get {x y z}]\n"
00899                   "$cur delete; mol delete [$new molid]; $new delete\n\n");
00900     vmdcon_printf(VMDCON_WARN, 
00901                   "cpmdlogplugin) ##################################"
00902                   "####################################\n");
00903     /* Try to read wavefunction and orbital energies */
00904     SAFE_CALLOC(data->qm_timestep,qm_timestep_t,1);
00905     if (get_wavefunction(data, data->qm_timestep, NULL) == FALSE) {
00906       vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) No wavefunction present for timestep %d\n", data->num_frames_read);
00907       free(data->qm_timestep);
00908       data->qm_timestep=NULL;
00909     } else {
00910       vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) Wavefunction found.\n");
00911       data->num_frames_read=1;
00912       data->num_frames=1;
00913     }
00914   } else {
00915     /* don't copy coordinates yet. we have the newer version of CPMD
00916      * that prints them right before each projection. */
00917     data->num_frames = 0;
00918   }
00919   
00920   return TRUE;
00921 }
00922 
00923 
00924 /*******************************************************
00925  *
00926  * this function reads in the basis set data from 
00927  * <basis>.gbs or $VMDDIR/basis/<basis>.gbs
00928  *
00929  * XXX: this is the same function as in gaussianplugin
00930  * ******************************************************/
00931 int get_internal_basis(gaussiandata *data) {
00932 
00933   char *vmddir=NULL;
00934   FILE *fp;
00935   char buffer[BUFSIZ];
00936   char word[3][MOLFILE_BUFSIZ];
00937   char filepath[256];
00938   int  i,n, wavef_size; 
00939 
00940   /* no point in adding a basis set if we already have this information */
00941   if (data->have_basis) return TRUE;
00942 
00943   /* try to open basis set database file. a file in the current
00944    * directory takes priority over what is shipped with VMD. */
00945   sprintf(filepath,"%s.gbs",data->gbasis);
00946   fp=fopen(filepath,"rb");
00947   if (fp == NULL) {
00948       vmddir=getenv("VMDDIR");
00949       if (vmddir == NULL) {
00950           vmddir="/usr/local/lib/vmd";
00951       }
00952       sprintf(filepath,"%s/basis/%s.gbs",vmddir,data->gbasis);
00953       fp=fopen(filepath,"rb");
00954   }
00955 
00956   if (fp == NULL) {
00957     vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) failed to read basis set "
00958                   "from data base file %s\n", filepath);
00959     data->num_shells_per_atom=NULL;
00960     data->have_basis=FALSE;
00961     return FALSE;
00962   } else {
00963     vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) reading basis set "
00964                   "from data base file %s\n", filepath);
00965   }
00966   data->wavef_size=0;
00967   
00968   /* Allocate space for the basis for all atoms */
00969   /* When the molecule is symmetric the actual number atoms with
00970    * a basis set could be smaller */
00971   SAFE_CALLOC(data->basis_set,basis_atom_t,data->numatoms);
00972 
00973   for (i=0; i < data->numatoms; ++i) {
00974     int numshells, numprim;
00975     int numread, ishell;
00976     float scalef;
00977     shell_t *shell;
00978     
00979     /* search for the characteristic first line starting with '****'. */
00980     rewind(fp);
00981     do {
00982       fgets(buffer, sizeof(buffer), fp);
00983       sscanf(buffer,"%s%s",word[0],word[1]);
00984     } while(strcmp(word[0],"****"));
00985   
00986     /* search for an entry for the current atom in the format '<name> 0'. */
00987     do {
00988       fgets(buffer, sizeof(buffer), fp);
00989       if (feof(fp)) {
00990         free(data->basis_set);
00991         data->basis_set=NULL;
00992         data->num_shells_per_atom=NULL;
00993         data->have_basis=FALSE;
00994         vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) EOF in data base "
00995                       "file %s while looking for element %s.\n", filepath, 
00996                       get_pte_label(data->initatoms[i].atomicnum), buffer);
00997         fclose(fp);
00998         return FALSE;
00999       }
01000       n=sscanf(buffer,"%s%s",word[0],word[1]);
01001     } while ( (n != 2) || strcmp(word[1],"0") ||
01002               (strcmp(word[0],get_pte_label(data->initatoms[i].atomicnum))) );
01003     
01004     strncpy(data->basis_set[i].name,data->gbasis,sizeof(data->gbasis));
01005     
01006     numshells=0;
01007     shell=NULL;
01008 
01009     /* read basis set until end of element */
01010     do {
01011 
01012       fgets(buffer, sizeof(buffer), fp);
01013       if (strstr(buffer,"****")) break;
01014       if (ferror(fp)) {
01015         vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) read error in data "
01016                       "base file %s while reading basis of element %s.\n", 
01017                       filepath, get_pte_label(data->initatoms[i].atomicnum));
01018         free(data->basis_set);
01019         data->basis_set=NULL;
01020         return FALSE;
01021       }
01022 
01023       numread=sscanf(buffer,"%s%d%f",word[0],&numprim,&scalef);
01024       if (numread == 3) {
01025 #if CPMDLOG_DEBUG && CPMDLOG_BASIS_DEBUG
01026         vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) atom: %d, element: %s, shell: %d "
01027                       "%s-type shell, %d primitives, scalefactor %f\n", i,
01028                       get_pte_label(data->initatoms[i].atomicnum), numshells+1, 
01029                       word[0], numprim, scalef);
01030 #endif
01031         ;
01032       } else {
01033         vmdcon_printf(VMDCON_ERROR, 
01034                       "cpmdlogplugin) basis set parse error: %s",buffer);
01035         free(data->basis_set);
01036         data->basis_set=NULL;
01037         return FALSE;
01038       }
01039 
01040       ishell=numshells;
01041       ++numshells;
01042       shell=realloc(shell,numshells*sizeof(shell_t));
01043       shell[ishell].numprims=numprim;
01044       shell[ishell].symmetry=shellsymm_int(word[0]);
01045       shell[ishell].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01046       data->num_basis_funcs += numprim;
01047 
01048       switch(shell[ishell].symmetry) {
01049         case S_SHELL:
01050           data->wavef_size += 1;
01051           break;
01052         case P_SHELL:
01053           data->wavef_size += 3;
01054           break;
01055         case SP_S_SHELL:
01056           data->wavef_size += 4;
01057           break;
01058         case D_SHELL:
01059           data->wavef_size += 6;  /* cartesian representation! */
01060           break;
01061         case SPD_S_SHELL:
01062           data->wavef_size += 10;  /* cartesian representation! */
01063           break;
01064         case F_SHELL:
01065           data->wavef_size += 10;  /* cartesian representation! */
01066           break;
01067         default:
01068           break;
01069       }
01070       
01071       if (shell[ishell].symmetry == SP_S_SHELL) {
01072         ++numshells;
01073         shell=realloc(shell,numshells*sizeof(shell_t));
01074         shell[ishell+1].numprims=numprim;
01075         shell[ishell+1].symmetry=SP_P_SHELL;
01076         shell[ishell+1].prim = (prim_t *)calloc(numprim,sizeof(prim_t));
01077         data->num_basis_funcs += numprim;
01078       }
01079 
01080       for (n=0; n<numprim; ++n) {
01081         fgets(buffer, sizeof(buffer), fp);
01082         if (ferror(fp)) {
01083           vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) read error in data "
01084                         "base file %s while reading basis of element %s.\n", 
01085                         filepath, get_pte_label(data->initatoms[i].atomicnum));
01086           free(data->basis_set);
01087           data->basis_set=NULL;
01088           return FALSE;
01089         }
01090         sscanf(buffer,"%s%s%s", word[0],word[1],word[2]);
01091         fix_fortran_exp(word[0]);
01092         shell[ishell].prim[n].exponent=atof(word[0])*scalef*scalef;
01093         fix_fortran_exp(word[1]);
01094         shell[ishell].prim[n].contraction_coeff=atof(word[1]);
01095         if (shell[ishell].symmetry == SP_S_SHELL) {
01096           shell[ishell+1].prim[n].exponent=shell[ishell].prim[n].exponent;
01097           fix_fortran_exp(word[2]);
01098           shell[ishell+1].prim[n].contraction_coeff=atof(word[2]);
01099         }
01100       }
01101     } while(1);
01102 
01103   
01104     /* store shells in atom */
01105     data->basis_set[i].numshells = numshells;
01106     data->basis_set[i].shell = shell;
01107     
01108     /* store the total number of basis functions */
01109     data->num_shells += numshells;
01110   }
01111   
01112   vmdcon_printf(VMDCON_INFO, 
01113                 "cpmdlogplugin) Parsed %d uncontracted basis functions. \n", 
01114                 data->num_basis_funcs);
01115 
01116   /* allocate and populate flat arrays needed for molfileplugin */
01117   data->have_basis = TRUE;
01118   return fill_basis_arrays(data);
01119 }
01120 
01121 
01122 /**************************************************
01123  * Convert shell symmetry type from char to int.
01124  * XXX: same function exists in gaussianplugin.
01125  ************************************************ */
01126 static int shellsymm_int(char *symm) {
01127   int shell_symmetry;
01128 
01129   switch (toupper(symm[0])) {
01130     case 'S':
01131       if (symm[1] == '\0') {
01132         shell_symmetry = S_SHELL;
01133       } else if (toupper(symm[1]) == 'P') {
01134         if (symm[2] == '\0') {
01135           shell_symmetry = SP_S_SHELL;
01136         } else if (toupper(symm[1]) == 'D') {
01137           shell_symmetry = SPD_S_SHELL;
01138         } else {
01139           shell_symmetry = UNK_SHELL;
01140         } 
01141       } else {
01142         shell_symmetry = UNK_SHELL;
01143       }
01144       break;
01145     case 'L':
01146       shell_symmetry = SP_S_SHELL;
01147       break;
01148     case 'M': 
01149       shell_symmetry = SP_P_SHELL;
01150       break;
01151     case 'P':
01152       shell_symmetry = P_SHELL;
01153       break;
01154     case 'D':
01155       shell_symmetry = D_SHELL;
01156       break;
01157     case 'F':
01158       shell_symmetry = F_SHELL;
01159       break;
01160     case 'G':
01161       shell_symmetry = G_SHELL;
01162       break;
01163     default:
01164       shell_symmetry = UNK_SHELL;
01165       break;
01166   }
01167 
01168   return shell_symmetry;
01169 }
01170 
01171 
01173 static int fill_basis_arrays(gaussiandata *data) {
01174   int i, j, k;
01175   int shellcount = 0;
01176   int primcount = 0 ;
01177   float *basis;
01178   int *num_shells_per_atom;
01179   int *num_prim_per_shell;
01180   int *shell_symmetry;
01181 
01182   /* reserve space for pointer to array containing basis
01183    * info, i.e. contraction coeficients and expansion 
01184    * coefficients; need 2 entries per primitive gaussian, i.e.
01185    * exponent and contraction coefficient; also,
01186    * allocate space for the array holding the orbital symmetry
01187    * information per primitive gaussian.
01188    * Finally, initialize the arrays holding the number of 
01189    * shells per atom and the number of primitives per shell*/
01190   SAFE_CALLOC(basis,float,2*data->num_basis_funcs);
01191   SAFE_CALLOC(shell_symmetry,int,data->num_shells);
01192   SAFE_CALLOC(num_shells_per_atom,int,data->numatoms);
01193   SAFE_CALLOC(num_prim_per_shell,int,data->num_shells);
01194   
01195   /* place pointers into struct gaussiandata */
01196   data->basis = basis;
01197   data->shell_symmetry = shell_symmetry;
01198   data->num_shells_per_atom = num_shells_per_atom;
01199   data->num_prim_per_shell = num_prim_per_shell;
01200 
01201   for(i=0; i<data->numatoms; i++) {
01202     num_shells_per_atom[i] = data->basis_set[i].numshells;
01203 
01204     for (j=0; j<data->basis_set[i].numshells; j++) {
01205       shell_symmetry[shellcount] = data->basis_set[i].shell[j].symmetry;
01206       num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
01207 
01208       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
01209         basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
01210         basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
01211         primcount++;
01212       }
01213       shellcount++;
01214     }
01215   } 
01216   return TRUE;
01217 }
01218 
01221 static int get_traj_frame(gaussiandata *data) {
01222   qm_timestep_t *cur_qm_ts;
01223   char buffer[BUFSIZ];
01224   char word[4][MOLFILE_BUFSIZ];
01225   int n,i;
01226   
01227 
01228   buffer[0] = '\0';
01229   word[0][0] = '\0';
01230   word[1][0] = '\0';
01231   word[2][0] = '\0';
01232   word[3][0] = '\0';
01233 
01234 #if CPMDLOG_DEBUG
01235   vmdcon_printf(VMDCON_INFO, 
01236                 "cpmdlogplugin) Timestep %d: ====\n", 
01237                 data->num_frames_read);
01238 #endif
01239 
01240   /* allocate more memory for the timestep array */
01241   data->qm_timestep = 
01242     (qm_timestep_t *)realloc(data->qm_timestep, 
01243                              (data->num_frames_read+1)*sizeof(qm_timestep_t));
01244 
01245   /* get a convenient pointer to the current qm timestep */
01246   cur_qm_ts = data->qm_timestep+data->num_frames_read;
01247   memset(cur_qm_ts, 0, sizeof(qm_timestep_t));
01248 
01249   /* search for data */
01250   while (1) {
01251 
01252     GET_LINE(buffer,data->file);
01253     n = sscanf(buffer,"%s%s%s%s",word[0],word[1],word[2],word[3]);
01254     
01255     /* empty or uninteresting line */
01256     if (n < 3) continue;
01257 
01258     /* coordinates relevant for projection.
01259      * this needs a CPMD version > 3.13.2
01260  ********************* PROJECTION COORDINATES ********************
01261    NR   TYPE        X(bohr)        Y(bohr)        Z(bohr)    
01262     1      B       8.536664      10.525423      10.202766
01263     2      B      10.195353       8.831898      12.528872
01264     3      H      10.467571       8.964546       9.816927
01265     4      H       8.327684      10.371750      12.680045
01266     5      H      12.222474       9.664531      13.239192
01267     6      H       9.579094       6.722748      12.777072
01268     7      H       6.516216       9.875605       9.562778
01269     8      H       9.404536      12.649641       9.952504
01270  ****************************************************************
01271  */
01272     if ( (strstr(word[0],"************") != 0 &&
01273           strcmp(word[1],"PROJECTION"  ) == 0 && 
01274           strcmp(word[2],"COORDINATES" ) == 0 && 
01275           strstr(word[3],"************") != 0 ) ) {
01276 
01277       /* NR TYPE X(bohr) Y(bohr) Z(bohr) */
01278       GET_LINE(buffer, data->file);
01279 
01280       if (data->initatoms==NULL) {
01281         vmdcon_printf(VMDCON_ERROR,"why is initatoms NULL?\n");
01282         return FALSE;
01283       }
01284 
01285       for (i=0; i < data->numatoms; ++i) {
01286         qm_atom_t *atm;
01287       
01288         GET_LINE(buffer, data->file);
01289         atm = data->initatoms + i;
01290         
01291         n=sscanf(buffer,"%*d%*s%g%g%g", &atm->x, &atm->y, &atm->z);
01292         
01293         if (n != 3) {
01294           vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) Failed to parse "
01295                         "projection coordinates. Stopping.\n");
01296           return FALSE;
01297         }
01298         atm->atomicnum=get_pte_idx(atm->type);
01299         atm->x *= BOHR_TO_ANGS;
01300         atm->y *= BOHR_TO_ANGS;
01301         atm->z *= BOHR_TO_ANGS;
01302       }
01303       buffer[0] = '\0';
01304 #if CPMDLOG_DEBUG
01305       vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) Coordinates found for timestep %d\n", data->num_frames_read);
01306 #endif
01307     } else if ( (strcmp(word[0],"WAVEFUNCTIONS") == 0) &&
01308          (strcmp(word[2],"ATOMIC"       ) == 0) &&
01309          (strcmp(word[3],"ORBITAL"      ) == 0) ) {
01310 
01311       /* Try to read wavefunction and orbital energies */
01312       if (get_wavefunction(data, cur_qm_ts, NULL) == FALSE) {
01313           vmdcon_printf(VMDCON_WARN, "cpmdlogplugin) No wavefunction present for timestep %d\n", data->num_frames_read);
01314           /* XXX: add flag to ignore wfn */
01315       }
01316       if (data->have_wavefunction == 2) {
01317 #if 0          
01318           /* Try to read localized wavefunctions  */
01319           /* XXX: used offset or something. */
01320           if (get_wavefunction(data, cur_qm_ts+XXX) == FALSE) {
01321               vmdcon_printf(VMDCON_WARN, "cpmdlogplugin) No localized wavefunction present for timestep %d\n", data->num_frames_read);
01322           }
01323 #endif
01324       }
01325       buffer[0] = '\0';
01326       break; /* XXX */
01327     } else if ( (strcmp(word[0],"POPULATION") == 0) &&
01328                 (strcmp(word[1],"ANALYSIS"  ) == 0) &&
01329                 (strcmp(word[3],"PROJECTED" ) == 0) ) {
01330 #if 0
01331       if (get_population(data, cur_qm_ts)) {
01332         vmdcon_printf(VMDCON_INFO, "cpmdlogplugin) Mulliken charges found\n");
01333       }
01334 #endif
01335       buffer[0] = '\0';
01336       
01337     } else if ( strcmp(word[0],"*"     ) == 0 &&
01338                 strcmp(word[1],"TIMING") == 0 && 
01339                 strcmp(word[2],"*"     ) == 0 ) {
01340       data->end_of_trajectory=FALSE;
01341       break;
01342     } else if (feof(data->file)) {
01343       data->end_of_trajectory=TRUE;
01344       break;
01345     }
01346   }
01347 
01348   /* next timestep or end of file */
01349   data->num_frames_read++;
01350 
01351   return TRUE;
01352 }
01353 
01354 
01355 
01356 /*********************************************************
01357  *
01358  * this function reads the actual wavefunction, which is
01359  * punched at the end of the log file
01360  *
01361  **********************************************************/
01362 static int get_wavefunction(gaussiandata *data, qm_timestep_t *ts, 
01363                             qm_wavefunction_t *wf)
01364 {
01365   float *orbital_energies;
01366   float *wave_function;
01367   char buffer[BUFSIZ];
01368 #define ORBSPERBLOCK 8
01369   char word[ORBSPERBLOCK+1][MOLFILE_BUFSIZ];
01370   int orbital_counter = 0;
01371   int i = 0, j = 0, num_values = 0;
01372   int total_orbs = data->occ_orbitals_A + data->occ_orbitals_B;
01373   int num_orbs;
01374 
01375   if (wf == NULL) {
01376     PRINTERR;
01377     return FALSE;
01378   }
01379 
01380   wf->type = MOLFILE_WAVE_CANON;
01381   wf->spin = SPIN_ALPHA;
01382   
01383   buffer[0] = '\0';
01384   for (i=0; i <= ORBSPERBLOCK ; i++) word[i][0] = '\0';
01385   /*
01386    * Scan for something like this:
01387       ORBITAL      1       2       3       4       5       6       7       8
01388   COMPLETNESS    0.972   0.951   0.972   0.972   0.951   0.972   0.951   0.972
01389   OCCUPATION     2.000   2.000   2.000   2.000   2.000   2.000   2.000   2.000
01390   1   B  S      -0.064  -0.163   0.064   0.365   0.163  -0.365   0.163  -0.365
01391          Px      0.000  -0.170   0.000   0.000  -0.170   0.000  -0.170   0.000
01392          Pz      0.086   0.227  -0.086   0.181  -0.227  -0.181  -0.227  -0.181
01393          Py      0.043   0.000   0.043   0.411   0.000   0.411   0.000   0.411
01394   2   B  S       0.365  -0.163  -0.365  -0.064   0.163   0.064   0.163   0.064
01395          Px      0.000  -0.170   0.000   0.000  -0.170   0.000  -0.170   0.000
01396          Pz     -0.181  -0.227   0.181  -0.086   0.227   0.086   0.227   0.086
01397          Py     -0.411   0.000  -0.411  -0.043   0.000  -0.043   0.000  -0.043
01398   3   H  S      -0.049   0.200   0.049  -0.049   0.540   0.049   0.540   0.049
01399   4   H  S      -0.049  -0.540   0.049  -0.049  -0.200   0.049  -0.200   0.049
01400   5   H  S       0.520   0.078   0.041   0.056  -0.078   0.031  -0.078   0.031
01401     ...
01402    */
01403 
01404 
01405   /* Reserve space for arrays storing wavefunction and orbital
01406    * energies. For the wavefunction we reserve the number
01407    * of alpha orbitals squared, for unrestricted twice as much.
01408    * Accordingly, for the energies I use the number of A orbitals 
01409    * In gaussian the number of alpha and beta orbitals is always 
01410    * the same. */
01411 
01412 #if 0
01413   SAFE_CALLOC(wave_function,float,data->wavef_size*num_orbs);
01414   SAFE_CALLOC(orbital_energies,float,num_orbs);
01415   ts->wave_function    = wave_function;
01416   ts->orbital_energies = orbital_energies;
01417 
01418 
01419   while (orbital_counter < num_orbs) {
01420     /* read up to line of orbital indices */
01421     do {
01422       GET_LINE(buffer, data->file);
01423       sscanf(buffer,"%s",word[0]);
01424     } while (strcmp(word[0],"ORBITAL"));
01425 
01426     GET_LINE(buffer, data->file);           /* completeness */
01427     num_orbs = sscanf(buffer,"%*s%s%s%s%s%s%s%s%s",word[0],word[1],
01428                       word[2],word[3],word[4],word[5],word[6],word[7]);
01429 
01430     /* we don't have orbital energies here, but we
01431      * can use it for the completenes parameter. */
01432     for(i=0; i<num_orbs; i++) 
01433       *(orbital_energies+i) = atof(word[i]);
01434 
01435     GET_LINE(buffer, data->file);           /* occupation */
01436              
01437     /* step orbital energy pointer */
01438     orbital_energies = orbital_energies+num_orbs;
01439     orbital_counter += num_orbs;
01440 
01441     /* read in the wavefunction */
01442     for (i=0; i<data->wavef_size; i++) {
01443       int xexp=0, yexp=0, zexp=0;
01444 
01445       /* read in the wavefunction coefficients for up 
01446        * to 8 orbitals at a time line by line */
01447       GET_LINE(buffer, data->file);
01448       num_values = sscanf(buffer+7,"%4s%s%s%s%s%s%s%s%s", word[0], 
01449                           word[1], word[2],word[3], word[4], 
01450                           word[5], word[6], word[7], word[8]);
01451 
01452       /* handle magenetic quantum number. in cartesian basis the
01453        * labels are: S, PX, PY, PZ, DXX, DXY, DXZ, DYY, DYZ, DZZ, ...*/
01454       for (j=1; j<strlen(word[0]); j++) {
01455         switch (toupper(word[0][j])) {
01456           case 'X':
01457             xexp++;
01458             break;
01459           case 'Y':
01460             yexp++;
01461             break;
01462           case 'Z':
01463             zexp++;
01464             break;
01465             /* if we have pure d/f-functions the nomenclature changes to 
01466              * 'D 0', 'D-1', 'D+1', 'D-2', 'D+2' */
01467           case '+': /* fallthrough */
01468           case '-': /* fallthrough */
01469           case '0': /* fallthrough */
01470           case '1': /* fallthrough */
01471           case '2': /* fallthrough */
01472             vmdcon_printf(VMDCON_ERROR, "cpmdlogplugin) pure basis function "
01473                           "detected: '%s'. bailing out...\n", word[0]);
01474             free(ts->wave_function);
01475             wave_function=NULL;
01476             free(ts->orbital_energies);
01477             ts->orbital_energies=NULL;
01478             return FALSE;
01479             break;
01480           default:
01481             /* do nothing */
01482             break;
01483         }
01484       }
01485       data->angular_momentum[3*i  ] = xexp;
01486       data->angular_momentum[3*i+1] = yexp;
01487       data->angular_momentum[3*i+2] = zexp;
01488 #if CPMDLOG_DEBUG && CPMDLOG_BASIS_DEBUG
01489       fprintf(stderr,"%s:%d orbital %d/%d  function %d/%s: %d %d %d\n", 
01490               __FILE__, __LINE__, num_orbs, orbital_counter, 
01491               i, word[0], xexp, yexp, zexp);
01492 #endif
01493 
01494       /* each orbital has data->wavef_size entries, 
01495        * hence we have to use this number as offset when storing 
01496        * them in groups of five */
01497       for (j=0 ; j<num_values-1; j++) {
01498         wave_function[j*data->wavef_size+i] = atof(word[j+1]);
01499       }
01500     }
01501     /* move wavefunction pointer to start of next block of orbitals */
01502     wave_function = wave_function + num_orbs*data->wavef_size;
01503     /* XXX: FIXME test with open shell run */
01504     if ((data->scftyp == SCFTYP_UHF) && (orbital_counter==data->wavef_size)) {
01505       GET_LINE(buffer, data->file);
01506     }
01507   }
01508 
01509   /* store the number of orbitals read in */
01510   ts->wavef_size = orbital_counter;
01511 #endif
01512 
01513 #if CPMDLOG_DEBUG
01514   vmdcon_printf(VMDCON_INFO, 
01515                 "cpmdlogplugin) Number of orbitals scanned: %d \n",
01516                 orbital_counter);
01517 #endif
01518   return TRUE;
01519 }
01520 
01521 
01522 /* Read the population analysis section.
01523  * Currently we parse only the Mulliken charges
01524  * but we might want to add support for populations
01525  * and for Lowdin analysis. */
01526 static int get_population(gaussiandata *data, qm_timestep_t *ts) {
01527   int i;
01528   char buffer[BUFSIZ];
01529   data->have_mulliken = FALSE;
01530 
01531 
01532   /* Read Mulliken charges if present */
01533   ts->mulliken_charges = 
01534     (double *)calloc(data->numatoms, sizeof(double));
01535 
01536   if (!ts->mulliken_charges) {
01537     PRINTERR; 
01538     return FALSE;
01539   }
01540   
01541   for (i=0; i<data->numatoms; i++) {
01542     int n;
01543     float mullpop, mullcharge, lowpop, lowcharge;
01544     GET_LINE(buffer, data->file);
01545     n = sscanf(buffer,"%*i%*s%f%f%f%f",
01546                &mullpop, &mullcharge, &lowpop, &lowcharge);
01547     if (n!=4) return FALSE;
01548     ts->mulliken_charges[i] = mullcharge;
01549   }
01550 
01551   if (i!=data->numatoms) return FALSE;
01552 
01553   data->have_mulliken = TRUE;
01554   return TRUE;
01555 }
01556 
01557 /*************************************************************
01558  *
01559  * plugin registration 
01560  *
01561  **************************************************************/
01562 
01563 VMDPLUGIN_API int VMDPLUGIN_init(void) {
01564   memset(&plugin, 0, sizeof(molfile_plugin_t));
01565   plugin.abiversion = vmdplugin_ABIVERSION;
01566   plugin.type = MOLFILE_PLUGIN_TYPE;
01567   plugin.name = "cpmdlog";
01568   plugin.prettyname = "CPMD 3.x output ";
01569   plugin.author = "Axel Kohlmeyer";
01570   plugin.majorv = 0;
01571   plugin.minorv = 0;
01572   plugin.is_reentrant = VMDPLUGIN_THREADUNSAFE;
01573   plugin.filename_extension = "out";
01574   plugin.open_file_read = open_cpmdlog_read;
01575   plugin.read_structure = read_cpmdlog_structure;
01576   plugin.close_file_read = close_cpmdlog_read;
01577 
01578   plugin.read_qm_metadata = read_cpmdlog_metadata;
01579   plugin.read_qm_rundata  = read_cpmdlog_rundata;
01580 
01581 #if vmdplugin_ABIVERSION > 10
01582   plugin.read_timestep_metadata    = read_timestep_metadata;
01583 #endif
01584 #if vmdplugin_ABIVERSION > 11
01585   plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
01586   plugin.read_timestep = read_timestep;
01587 #endif
01588 
01589   return VMDPLUGIN_SUCCESS;
01590 }
01591 
01592 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01593   (*cb)(v, (vmdplugin_t *)&plugin);
01594   return VMDPLUGIN_SUCCESS;
01595 }
01596 
01597 VMDPLUGIN_API int VMDPLUGIN_fini(void) {
01598   return VMDPLUGIN_SUCCESS;
01599 }
01600 
01601 
01602 #ifdef TEST_PLUGIN
01603 
01604 int main(int argc, char *argv[]) {
01605   int numatoms, i, j, optflags;
01606   molfile_atom_t *atoms;
01607   molfile_timestep_t timestep;
01608   molfile_metadata_t metadata;
01609   molfile_timestep_metadata_t ts_metadata;
01610   molfile_qm_timestep_metadata_t qm_ts_metadata;
01611   molfile_qm_metadata_t qm_metadata;
01612   molfile_qm_timestep_t qm_ts;
01613  
01614   void *v;
01615 
01616   while (--argc) {
01617     ++argv;
01618     v = open_cpmdlog_read(*argv, "log", &numatoms);
01619     if (!v) {
01620       fprintf(stderr, "open_cpmdlog_read failed for file %s\n", *argv);
01621       return 1;
01622     }
01623     fprintf(stderr, "open_cpmdlog_read succeeded for file %s\n", *argv);
01624     fprintf(stderr, "number of atoms: %d\n", numatoms);
01625     atoms = (molfile_atom_t *)malloc(sizeof(molfile_atom_t)*numatoms);
01626     read_cpmdlog_structure(v,&optflags, atoms);
01627 
01628     i = 0;
01629     timestep.coords = (float *)malloc(3*sizeof(float)*numatoms);
01630 
01631     
01632     while (1) {
01633       memset(&ts_metadata, 0, sizeof(molfile_timestep_metadata_t));
01634       read_timestep_metadata(v, &ts_metadata);
01635       memset(&qm_metadata, 0, sizeof(molfile_qm_metadata_t));
01636       read_qm_timestep_metadata(v, &qm_ts_metadata);
01637       qm_ts.scfenergies = (double *)malloc(qm_ts_metadata.num_scfiter*sizeof(double));
01638       qm_ts.wave        = (molfile_qm_wavefunction_t *)malloc(qm_ts_metadata.num_wavef
01639                                                               *sizeof(molfile_qm_wavefunction_t));
01640       memset(qm_ts.wave, 0, qm_ts_metadata.num_wavef*sizeof(molfile_qm_wavefunction_t));
01641       for (j=0; (j<10 && j<qm_ts_metadata.num_wavef); j++) {
01642         qm_ts.wave[j].wave_coeffs = (float *) malloc(qm_ts_metadata.num_orbitals_per_wavef[j]
01643                                                      * qm_ts_metadata.wavef_size * sizeof(float));
01644         qm_ts.wave[j].orbital_energies = (float *) malloc(qm_ts_metadata.num_orbitals_per_wavef[j]*sizeof(float));
01645       }
01646       qm_ts.gradient = (float *) malloc(3*numatoms*sizeof(float));
01647       if (read_timestep(v, numatoms, &timestep, &qm_metadata, &qm_ts)) break;
01648       /* do something with data */
01649       /* XXX */
01650 
01651       free(qm_ts.gradient);
01652       for (j=0; (j<10 && j<qm_ts_metadata.num_wavef); j++) {
01653         free(qm_ts.wave[j].wave_coeffs);
01654         free(qm_ts.wave[j].orbital_energies);
01655       }
01656       free(qm_ts.wave);
01657       free(qm_ts.scfenergies);
01658       i++;
01659     }
01660     fprintf(stderr, "ended read_timestep on frame %d\n", i);
01661     close_cpmdlog_read(v);
01662   }
01663   return 0;
01664 }
01665 
01666 #endif

Generated on Fri Mar 29 03:10:04 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002