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

moldenplugin.c

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2016 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  * RCS INFORMATION:
00011  *
00012  *      $RCSfile: moldenplugin.c,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.41 $       $Date: 2016/12/15 15:32:25 $
00015  *
00016  ***************************************************************************/
00017 
00018 /* This is a plugin that will read input from a MOLDEN
00019 ** generated output file 
00020 ** some more details will go here soon 
00021 ** NOTE: The current version of the plugin relies
00022 ** on the fact that the [Atom] field comes before
00023 ** the [Geometries] field */
00024 
00025 
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <ctype.h>
00029 #include <string.h>
00030 
00031 #include <string.h>
00032 
00033 #if defined(_AIX)
00034 #include <strings.h>
00035 #endif
00036 
00037 #if defined(WIN32) || defined(WIN64)
00038 #define strcasecmp stricmp
00039 #endif
00040 
00041 #include "molfile_plugin.h"
00042 #include "unit_conversion.h"
00043 #include "periodic_table.h"
00044 #include "qmplugin.h"
00045 
00046 
00047 #define ALLOCATE(array, type, size) \
00048   array = (type *)calloc(size, sizeof(type)); \
00049   if (array == NULL) { \
00050     fprintf(stderr, "moldenplugin) Memory allocation for %s failed!\n", #array); \
00051     return FALSE; \
00052   }
00053 
00054 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE
00055 
00056 
00057 /* I could use a flag already present in qmdata_t to 
00058  * indicate a trajectory, but I'm using moldendata_t 
00059  * to demonstrate the use of 
00060  * void *format_specific_data;
00061  * in qmdata_t as a means to store data specific to 
00062  * the plugin. 
00063  */
00064 typedef struct {
00065   long filepos_atoms;   /* [ATOMS] section */
00066   long filepos_geomxyz; /* [GEOMETRIES] XYZ section */
00067   long filepos_gto;     /* [GTO] section */
00068   long filepos_mo;      /* [MO] section */
00069   char units[16];
00070   int coordsonly;
00071 } moldendata_t;
00072 
00073 
00074 
00075 /* Read the basis set data */
00076 static int get_basis (qmdata_t *);
00077 static int shelltype_int(char *type);
00078 static int fill_basis_arrays(qmdata_t *data);
00079 
00080 static int read_geom_block(qmdata_t *data);
00081 static int read_molecular_orbitals(qmdata_t *data);
00082 static int read_wave_coeffs(FILE *file, qm_wavefunction_t *wave);
00083 static int count_orbitals(qmdata_t *data);
00084 
00085 
00086 /*********************************************************
00087  *
00088  * Open file and return number of atoms.
00089  *
00090  * In order to determine the # atom we have to read into
00091  * the [ATOMS]/[GEOMETRIES] sections where the atoms are
00092  * defined.
00093  * There can be either [ATOMS] or [GEOMETRIES] or both.
00094  * While [GEOMETRIES] has the # atoms as its first line,
00095  * we actually have to count lines in [ATOMS].
00096  *
00097  * We assume no particular order for the sections and scan
00098  * the entire file for the according keywords. In order to
00099  * save time when the section contents are actually read
00100  * we store the file pointers to the beginning of each
00101  * section in moldendata_t.
00102  *
00103  *********************************************************/
00104 static void *open_molden_read(const char *filename,
00105                               const char *filetype,
00106                               int *natoms) {
00107   FILE *fd;
00108   qmdata_t *data = NULL;
00109   moldendata_t *moldendata;
00110   char buffer[1024];
00111   char keystring[20];
00112 
00113   fd = fopen(filename, "rb");
00114   if (!fd) return NULL;
00115   
00116   /* allocate memory for main QM data structure */
00117   data = init_qmdata();
00118   if (!data) return NULL;
00119 
00120   data->file = fd;
00121 
00122   /* allocate GAMESS specific data */
00123   moldendata = (moldendata_t *)calloc(1, sizeof(moldendata_t));
00124   if (!moldendata) return NULL;
00125 
00126   data->format_specific_data = moldendata;
00127 
00128 
00129   /* Read first line */
00130   if (!fgets(buffer,1024,data->file)) return NULL;
00131 
00132   /* Check if the file is MOLDEN format */
00133   if (!strcmp(strtoupper(trimleft(trimright(buffer))), "[MOLDEN FORMAT]")) {
00134     printf("moldenplugin) Detected MOLDEN file format!\n");
00135   } else {
00136     printf("moldenplugin) The file is not in MOLDEN format!\n");
00137     return NULL;
00138   }
00139 
00140   eatwhitelines(data->file);
00141 
00142 
00143   /* Identify the different sections. */
00144   while (fgets(buffer,1024,data->file)) {
00145 
00146     /* Get key string (ignoring empty lines) */
00147     if (!sscanf(buffer, "%s", keystring)) continue;
00148 
00149     /* Quick test avoiding the uppercase transformation */
00150     if (keystring[0]!='[') continue;
00151 
00152     /* Make keystring upper case */
00153     strtoupper(keystring);
00154 
00155     if (!strcmp(keystring, "[5D]") || !strcmp(keystring, "[5D7F]") ||
00156         !strcmp(keystring, "[7F]") || !strcmp(keystring, "[5D10F]") ||
00157         !strcmp(keystring, "[9G]")) {
00158       printf("moldenplugin) Spherical harmonic basis found %s. \n", keystring);
00159       printf("moldenplugin)   Currently VMD handles only basis sets/wave functions\n");
00160       printf("moldenplugin)   with cartesian Gaussian functions.\n");
00161       printf("moldenplugin)   Loading coordinates only.\n");
00162       moldendata->coordsonly = 1;
00163     }
00164 
00165     if (!strcmp(keystring, "[ATOMS]")) {
00166       char *s;
00167       long prevline=ftell(fd);
00168       printf("moldenplugin) Found [ATOMS] section ...\n");
00169       moldendata->filepos_atoms = ftell(data->file);
00170 
00171       if (!sscanf(buffer, "%*s %s", moldendata->units)) {
00172         printf("moldenplugin) Missing units in [ATOMS] section!\n");
00173         return NULL;
00174       }
00175 
00176       /* start counting the atoms; 
00177        * read until I hit the first line that starts with a "["
00178        * bracket */      
00179       (*natoms) = 0; 
00180       s = fgets(buffer, 1024, fd);
00181 
00182       /* Here we assume that the [ATOMS] section goes
00183        * on until the next empty line or another section
00184        * starts, i.e. there is a "[" or we encounter EOF */
00185       while (trimleft(buffer)[0]!='[' && s!=NULL && !iswhiteline(buffer)) {
00186         (*natoms)++;
00187         prevline = ftell(fd);
00188         s = fgets(buffer, 1024, fd);
00189       }
00190       data->numatoms = *natoms;
00191       data->num_frames = 1;
00192 
00193       /* Go back to the previous line, it might contain
00194        * the next keyword */
00195       fseek(data->file, prevline, SEEK_SET);
00196     }
00197 
00198     else if (!strcmp(keystring, "[GEOMETRIES]")) {
00199       if (!strcmp(trimright(buffer), "[GEOMETRIES] XYZ")) {
00200         printf("moldenplugin) Found [GEOMETRIES] XYZ section ...\n");
00201 
00202         moldendata->filepos_geomxyz = ftell(data->file);
00203 
00204         /* The first line of the XYZ type [GEOMETRIES] input
00205          * contains the number of atoms. */
00206         if (fscanf(data->file, "%d", natoms) != 1) {
00207           printf("moldenplugin) No # atoms found in [GEOMETRIES] section!\n");
00208           return NULL;
00209         }
00210         data->numatoms = *natoms;
00211 
00212         /* Jump back to the beginning of the section */
00213         fseek(data->file, moldendata->filepos_geomxyz, SEEK_SET);
00214 
00215         /* Count # frames */
00216         data->num_frames = 0;
00217         do {
00218           int natm = 0;
00219           fscanf(data->file, "%d", &natm);
00220           if (natm!=data->numatoms) break;
00221           eatline(data->file, 1);
00222 
00223           data->filepos_array = (long*)realloc(data->filepos_array,
00224                                                (data->num_frames+1)*sizeof(long));
00225           data->filepos_array[data->num_frames] = ftell(data->file);
00226 
00227           /* Skip title line + numatoms lines */
00228           eatline(data->file, 1+data->numatoms);
00229           if (feof(data->file)) break;
00230 
00231           data->num_frames++;
00232         } while (1);
00233 
00234         printf("moldenplugin) Found %d frames\n", data->num_frames);
00235       } else if (!strcmp(trimright(buffer), "[GEOMETRIES] ZMAT")) {
00236         printf("moldenplugin) [GEOMETRIES] ZMAT not supported!\n");
00237       }
00238     }
00239 
00240     else if (!strcmp(keystring,"[GTO]")) {
00241       printf("moldenplugin) Found [GTO] section ...\n");
00242       moldendata->filepos_gto = ftell(data->file);
00243     }
00244 
00245     else if (!strcmp(keystring,"[MO]")) {
00246       printf("moldenplugin) Found [MO] section ...\n");
00247       moldendata->filepos_mo = ftell(data->file);
00248     }
00249 
00250   };
00251   
00252   return data;
00253 }
00254 
00255 
00256 
00257 /*********************************************************
00258  *
00259  * Read geometry from file
00260  *
00261  * The [ATOMS] section provides atom name, atomic number
00262  * and coordinates while the [GEOMETRIES] XYZ section
00263  * provides atom type and coordinates. Trajectories can
00264  * only be specified using [GEOMETRIES].
00265  * In case we only have [GEOMETRIES] the atomic number
00266  * will have to be deduced from the atom name.
00267  *
00268  *********************************************************/
00269 static int read_molden_structure(void *mydata, int *optflags, 
00270                                  molfile_atom_t *atoms) 
00271 {
00272   int i;
00273   char buffer[1024];
00274   char atname[1024];
00275   int num, atomicnum;
00276   molfile_atom_t *atom;
00277   qmdata_t *data = (qmdata_t *)mydata;
00278   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
00279 
00280   ALLOCATE(data->atoms, qm_atom_t, data->numatoms);
00281 
00282   /* atomic number is provided by plugin.
00283    * (This is required for QM plugins!) */
00284   *optflags = MOLFILE_ATOMICNUMBER;
00285 
00286   /* [ATOMS] section */
00287   if (moldendata->filepos_atoms) { 
00288     float unitfac = 1.f;
00289 
00290     /* If the units are given in AU we have to convert them.       */
00291     /* Note: Also recognize parenthesized units emitted by Molcas. */
00292     if (!strcmp(moldendata->units, "AU") ||
00293         !strcmp(moldendata->units, "(AU)")) {
00294       unitfac = BOHR_TO_ANGS;
00295     }
00296     
00297     /* Jump to beginning of [ATOMS] section. */
00298     fseek(data->file, moldendata->filepos_atoms, SEEK_SET);
00299 
00300     /* Read the atom types, names, atomic numbers
00301      * as well as x,y,z coordinates */
00302     for (i=0; i<data->numatoms; i++) {
00303       float x,y,z;
00304       atom = atoms+i;
00305       
00306       if (!fgets(buffer,1024,data->file)) return MOLFILE_ERROR;
00307       
00308       sscanf(buffer,"%s %d %d %f %f %f", atname, &num,
00309              &atomicnum, &x, &y, &z);
00310 
00311       /* populate data structure for VMD */
00312       strncpy(atom->name, atname,     sizeof(atom->name)); 
00313       strncpy(atom->type, atom->name, sizeof(atom->type));
00314       atom->atomicnumber = atomicnum;
00315       atom->resname[0] = '\0';
00316       atom->resid = 1;
00317       atom->chain[0] = '\0';
00318       atom->segid[0] = '\0';
00319 
00320       /* keep local copy */
00321       strncpy(data->atoms[i].type, atname, sizeof(data->atoms[i].type));
00322       data->atoms[i].atomicnum = atomicnum;
00323       data->atoms[i].x = x*unitfac;
00324       data->atoms[i].y = y*unitfac;
00325       data->atoms[i].z = z*unitfac;
00326     }
00327     data->num_frames_read = 1;
00328 
00329     return MOLFILE_SUCCESS;
00330   }
00331 
00332   /* [GEOMETRIES] XYZ section */
00333   if (moldendata->filepos_geomxyz) {
00334 
00335     /* Jump to beginning of [GEOMETRIES] section. */
00336     fseek(data->file, moldendata->filepos_geomxyz, SEEK_SET);
00337     eatline(data->file, 2);
00338 
00339     /* Read block from file */
00340     for (i=0; i<data->numatoms; i++) {
00341       atom = atoms+i;
00342       if (!fgets(buffer,1024,data->file)) return MOLFILE_ERROR;
00343       sscanf(buffer,"%s %*f %*f %*f", atname);
00344 
00345       strncpy(atom->type, atname, sizeof(atom->type));
00346       strncpy(atom->name, atname, sizeof(atom->name)); 
00347       atom->atomicnumber = get_pte_idx_from_string(atname);
00348       atom->resname[0] = '\0';
00349       atom->resid = 1;
00350       atom->chain[0] = '\0';
00351       atom->segid[0] = '\0';
00352       data->atoms[i].atomicnum = atom->atomicnumber;
00353     }
00354     data->num_frames_read = 0;
00355 
00356     return MOLFILE_SUCCESS;
00357   }
00358 
00359   printf("Sorry, could not obtain structure information \n");
00360   printf("from either the [ATOMS] or [GEOMETRIES] section! \n");
00361   printf("Please check your MOLDEN output file! \n"); 
00362   return MOLFILE_ERROR; 
00363 }
00364 
00365 
00366 /***********************************************************
00367  *
00368  * Read atoms for one frame from [GEOMETRIES] section.
00369  *
00370  ***********************************************************/
00371 static int read_geom_block(qmdata_t *data) {
00372   int i;
00373   char buffer[1024];
00374   float x,y,z;
00375 
00376   /* Skip title line */
00377   eatline(data->file, 1);
00378 
00379   for (i=0; i<data->numatoms; i++) {
00380     if (!fgets(buffer,1024,data->file)) return 0;
00381     sscanf(buffer,"%*s %f %f %f", &x, &y, &z);
00382     data->atoms[i].x = x;
00383     data->atoms[i].y = y;
00384     data->atoms[i].z = z;
00385   }
00386 
00387   return 1;
00388 }
00389 
00390 
00391 /***********************************************************
00392  *
00393  * Provide non-QM metadata for next timestep. 
00394  * Required by the plugin interface.
00395  *
00396  ***********************************************************/
00397 static int read_timestep_metadata(void *mydata,
00398                                   molfile_timestep_metadata_t *meta) {
00399   meta->count = -1;
00400   meta->has_velocities = 0;
00401 
00402   return MOLFILE_SUCCESS;
00403 }
00404 
00405 
00406 /***********************************************************
00407  *
00408  * We are not reading the coefficients themselves,
00409  * because that could require a large amount of memory.
00410  *
00411  ***********************************************************/
00412 static int read_qm_timestep_metadata(void *mydata,
00413                                     molfile_qm_timestep_metadata_t *meta) {
00414   qmdata_t *data = (qmdata_t *)mydata;
00415   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
00416 
00417   if (data->num_frames_sent >= data->num_frames) {
00418     /* All frames were sent. */
00419     return MOLFILE_ERROR;
00420   }
00421 
00422   /* Can't send metadata if only coordinates were read */
00423   if (moldendata->coordsonly) return MOLFILE_ERROR;
00424 
00425   /* Count the number of cartesian basis functions in 
00426      the basis set */
00427   if (data->num_frames_sent == data->num_frames-1) {
00428     int i;
00429     qm_timestep_t *cur_ts;
00430 
00431     if (!count_orbitals(data)) return MOLFILE_ERROR;
00432 
00433     /* get a pointer to the current qm timestep */
00434     cur_ts = data->qm_timestep;
00435     
00436     for (i=0; (i<MOLFILE_MAXWAVEPERTS && i<cur_ts->numwave); i++) {
00437       meta->num_orbitals_per_wavef[i] = cur_ts->wave[i].num_orbitals;
00438       meta->has_occup_per_wavef[i]    = cur_ts->wave[i].has_occup;
00439       meta->has_orben_per_wavef[i]    = cur_ts->wave[i].has_orben;
00440     }
00441     meta->wavef_size   = data->wavef_size;
00442     meta->num_wavef    = cur_ts->numwave;
00443     meta->num_scfiter  = cur_ts->num_scfiter;
00444     meta->has_gradient = FALSE;
00445     meta->num_charge_sets = 0;
00446   }
00447 
00448   return MOLFILE_SUCCESS;
00449 }
00450 
00451 
00452 
00453 /***********************************************************
00454  *
00455  * Provides VMD with the data of the next timestep.
00456  *
00457  ***********************************************************/
00458 static int read_timestep(void *mydata, int natoms, 
00459        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
00460                          molfile_qm_timestep_t *qm_ts) {
00461   int i;
00462   qmdata_t *data = (qmdata_t *)mydata;
00463   qm_timestep_t *cur_ts;
00464 
00465   if (data->num_frames_sent >= data->num_frames) {
00466     /* All frames were sent. */
00467     return MOLFILE_ERROR;
00468   }
00469 
00470   if (data->num_frames_sent == data->num_frames_read) {
00471     /* Read next coordinate block from file */
00472     fseek(data->file, data->filepos_array[data->num_frames_read], SEEK_SET);
00473     read_geom_block(data);
00474 
00475     /*printf("moldenplugin) Read frame %d\n", data->num_frames_read); */
00476     data->num_frames_read++;
00477   }
00478 
00479 
00480   /* Copy the coordinates */
00481   for (i=0; i<natoms; i++) {
00482     ts->coords[3*i  ] = data->atoms[i].x;
00483     ts->coords[3*i+1] = data->atoms[i].y;
00484     ts->coords[3*i+2] = data->atoms[i].z; 
00485   }
00486 
00487   /*printf("moldenplugin) Sent frame %d\n", data->num_frames_sent); */
00488   data->num_frames_sent++;
00489 
00490   /* In MOLDEN the MOs are listed only for the last frame */
00491   if (data->num_frames_sent == data->num_frames) {
00492 
00493     /* get a convenient pointer to the current qm timestep */
00494     cur_ts = data->qm_timestep;
00495 
00496     read_molecular_orbitals(data);
00497 
00498     /* store the wave function and orbital energies */
00499     if (cur_ts != NULL && cur_ts->wave != NULL) {
00500       for (i=0; i<cur_ts->numwave; i++) {
00501         qm_wavefunction_t *wave = &cur_ts->wave[i];
00502         qm_ts->wave[i].type         = wave->type;
00503         qm_ts->wave[i].spin         = wave->spin;
00504         qm_ts->wave[i].excitation   = wave->exci;
00505         qm_ts->wave[i].multiplicity = wave->mult;
00506         qm_ts->wave[i].energy       = wave->energy;
00507         strncpy(qm_ts->wave[i].info, wave->info, MOLFILE_BUFSIZ);
00508         
00509         if (wave->wave_coeffs) {
00510           memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs,
00511                  wave->num_orbitals*data->wavef_size*sizeof(float));
00512         }
00513         if (wave->orb_energies) {
00514           memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
00515                  wave->num_orbitals*sizeof(float));
00516         }
00517         if (wave->has_occup) {
00518           memcpy(qm_ts->wave[i].occupancies, wave->orb_occupancies,
00519                  wave->num_orbitals*sizeof(float));
00520         }
00521       }
00522     }
00523 
00524   }
00525 
00526   return MOLFILE_SUCCESS;
00527 }
00528   
00529 
00530 /*****************************************************
00531  *
00532  * Provide VMD with the sizes of the QM related
00533  * data structure arrays that need to be made
00534  * available.
00535  * Since we cannot determine the basis set meta data
00536  * without parsing the whole basis set section, we
00537  * read all basis set data here. The data is stored
00538  * in the qmdata_t structure for later use in
00539  * read_molden_rundata().
00540  *
00541  *****************************************************/
00542 static int read_molden_metadata(void *mydata, 
00543     molfile_qm_metadata_t *metadata) {
00544 
00545   qmdata_t *data;
00546   moldendata_t *moldendata;
00547   data = (qmdata_t *)mydata;
00548   moldendata = (moldendata_t *)data->format_specific_data;
00549 
00550 
00551   metadata->ncart = 0;
00552   metadata->nimag = 0;
00553   metadata->nintcoords = 0;
00554 
00555   metadata->have_sysinfo = 0;
00556   metadata->have_carthessian = 0;
00557   metadata->have_inthessian = 0;
00558   metadata->have_normalmodes = 0;
00559 
00560   metadata->num_basis_funcs = 0;
00561   metadata->num_basis_atoms = 0;
00562   metadata->num_shells = 0;
00563   metadata->wavef_size = 0;
00564 
00565   if (!moldendata->coordsonly) {
00566     /* Read the basis set */
00567     if (!get_basis(data)) return MOLFILE_ERROR; 
00568 
00569     /* orbital + basis set data */
00570     metadata->num_basis_funcs = data->num_basis_funcs;
00571     metadata->num_basis_atoms = data->num_basis_atoms;
00572     metadata->num_shells      = data->num_shells;
00573     metadata->wavef_size      = data->wavef_size;  
00574   }
00575 
00576   return MOLFILE_SUCCESS;
00577 }
00578 
00579 
00580 /******************************************************
00581  * 
00582  * Provide VMD with the static (i.e. non-trajectory)
00583  * data. That means we are filling the molfile_plugin
00584  * data structures.
00585  *
00586  ******************************************************/
00587 static int read_molden_rundata(void *mydata, 
00588                                molfile_qm_t *qm_data) {
00589   qmdata_t *data = (qmdata_t *)mydata;
00590   int i;
00591   molfile_qm_hessian_t *hessian_data;
00592   molfile_qm_basis_t   *basis_data;
00593   molfile_qm_sysinfo_t *sys_data;
00594 
00595   if (!qm_data) return MOLFILE_ERROR;
00596 
00597 
00598   hessian_data = &qm_data->hess;
00599   basis_data   = &qm_data->basis;
00600   sys_data     = &qm_data->run;
00601 
00602   sys_data->num_electrons = data->num_electrons;
00603   sys_data->totalcharge = data->totalcharge;
00604 
00605   /* Populate basis set data */
00606   if (data->num_basis_funcs) {
00607     for (i=0; i<data->num_basis_atoms; i++) {
00608       basis_data->num_shells_per_atom[i] = data->num_shells_per_atom[i];
00609       basis_data->atomic_number[i] = data->atomicnum_per_basisatom[i];
00610     }
00611     
00612     for (i=0; i<data->num_shells; i++) {
00613       basis_data->num_prim_per_shell[i] = data->num_prim_per_shell[i];
00614       basis_data->shell_types[i] = data->shell_types[i];
00615     }
00616     
00617     for (i=0; i<2*data->num_basis_funcs; i++) {
00618       basis_data->basis[i] = data->basis[i];
00619     }
00620 
00621     /* If we have MOs in the file we must provide the 
00622      * angular momentum exponents */
00623     if (data->angular_momentum) {
00624       for (i=0; i<3*data->wavef_size; i++) {
00625         basis_data->angular_momentum[i] = data->angular_momentum[i];
00626       }
00627     }
00628   }
00629 
00630   /* fill in molfile_qm_sysinfo_t */
00631   /*sys_data->runtype = data->runtype;
00632   sys_data->scftype = data->scftype;
00633   sys_data->nproc   = data->nproc;
00634   sys_data->num_electrons  = data->num_electrons;
00635   sys_data->totalcharge    = data->totalcharge;
00636   sys_data->num_occupied_A = data->num_occupied_A;
00637   sys_data->num_occupied_B = data->num_occupied_B;
00638   sys_data->status         = data->opt_status;
00639   */
00640   return MOLFILE_SUCCESS;
00641 }
00642 
00643 
00644 /**********************************************************
00645  *
00646  * close file and free memory
00647  *
00648  **********************************************************/
00649 static void close_molden_read(void *mydata) {
00650   int i, j;
00651   qmdata_t *data = (qmdata_t *)mydata;
00652 
00653   fclose(data->file);
00654 
00655   free(data->atoms);
00656   free(data->basis);
00657   free(data->shell_types);
00658   free(data->atomicnum_per_basisatom);
00659   free(data->num_shells_per_atom);
00660   free(data->num_prim_per_shell);
00661   free(data->angular_momentum);
00662 
00663   if (data->basis_set) {
00664     for(i=0; i<data->num_basis_atoms; i++) {
00665       for (j=0; j<data->basis_set[i].numshells; j++) {
00666         free(data->basis_set[i].shell[j].prim);
00667       }
00668       free(data->basis_set[i].shell);
00669     } 
00670     free(data->basis_set);
00671   }
00672 
00673   free(data->format_specific_data);
00674   free(data->filepos_array);
00675 
00676   if (data->qm_timestep != NULL) {
00677     for (j=0; j<data->qm_timestep[0].numwave; j++) {
00678       free(data->qm_timestep[0].wave[j].wave_coeffs);
00679       free(data->qm_timestep[0].wave[j].orb_energies);
00680       free(data->qm_timestep[0].wave[j].orb_occupancies);
00681     }
00682     free(data->qm_timestep[0].wave);
00683     free(data->qm_timestep);
00684   } else {
00685     printf("close_molden_read(): NULL qm_timestep!\n");
00686   }
00687 
00688   free(data);
00689 }
00690 
00691 
00692 
00693 /* ####################################################### */
00694 /*             End of API functions                        */
00695 /* The following functions actually do the file parsing.   */
00696 /* ####################################################### */
00697 
00698 
00699 /******************************************************
00700  * 
00701  * Format specification of the basis-set consisting of 
00702  * contracted Gaussian Type Orbitals.
00703  *
00704  * [GTO]
00705  * atom_sequence_number1 0
00706  * shell_label number_of_primitives 1.00
00707  * exponent_prim_1 contraction_coeff_1 (contraction_coeff_sp_1)
00708  * ...
00709  * <empty line>
00710  * atom_sequence_number2 0
00711  * shell_label number_of_primitives 1.00
00712  * exponent_prim_1 contraction_coeff_1 (contraction_coeff_1)
00713  * ...
00714  * <empty line>
00715  *
00716  * recognized shell_labels: s, p, d, f, sp, g
00717  *
00718  * For 'sp' shells two contraction coefficients must be given,
00719  * for both s and p functions. 
00720  * The 0 on the shell_number line and the 1.00 on the
00721  * shell_label line are no longer functional and can be
00722  * ignored.
00723  *
00724  * All workings with the [GTO] keyword are in Atomic Units.
00725  *
00726  ******************************************************/
00727 
00728 
00729 /*******************************************************
00730  *
00731  * Read the basis set data
00732  *
00733  * Format example:
00734  * [GTO]
00735  *   1 0
00736  *  s    2 1.00
00737  *   0.2738503300E+02  0.4301284983E+00
00738  *   0.4874522100E+01  0.6789135305E+00
00739  *  sp   2 1.00
00740  *   0.1136748200E+01  0.4947176920E-01  0.5115407076E+00
00741  *   0.2883094000E+00  0.9637824081E+00  0.6128198961E+00
00742  *
00743  *   2 0
00744  *  s    2 1.00
00745  *   0.1309756400E+01  0.4301284983E+00
00746  *   0.2331360000E+00  0.6789135305E+00
00747  *  ...
00748  *
00749  * qmdata_t provides hierarchical data structures for 
00750  * the basis set which are convenient for parsing. 
00751  * The molfile_plugin interface, however, requires flat
00752  * arrays, so after reading is done we have to populate
00753  * the according arrays using fill_basis_arrays().
00754  * 
00755  *******************************************************/
00756 static int get_basis(qmdata_t *data) {
00757   char buffer[1024];
00758   char shelltype[1024];
00759   int atomid, numprims;
00760   int i, j=0;
00761   int numshells;
00762   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
00763 
00764   /* XXX already initialized in open_molden_read() */
00765   data->num_shells = 0;
00766   data->num_basis_funcs = 0;
00767   data->num_basis_atoms = 0;
00768 
00769   /* initialize basis set the character array */
00770   memset(data->basis_string, 0, sizeof(data->basis_string));
00771 
00772 
00773   /* Place file pointer on line after the [GTO] keyword. */
00774   fseek(data->file, moldendata->filepos_gto, SEEK_SET);
00775  
00776   /* Allocate memory for the basis of all atoms. */
00777   ALLOCATE(data->basis_set, basis_atom_t, data->numatoms);
00778 
00779   /* Loop over all atoms. */
00780   for (i=0; i<data->numatoms; i++) {
00781 
00782     if (!fgets(buffer,1024,data->file)) return FALSE;
00783     sscanf(buffer,"%d %*d", &atomid);
00784 
00785     numshells = 0;
00786     data->basis_set[i].shell = NULL;
00787 
00788     /* Read an unknown number of shells */
00789     while (1) {
00790       shell_t *shell, *shell2=NULL;
00791 
00792       if (!fgets(buffer,1024,data->file)) return FALSE;
00793       
00794       /* Empty line signifies beginning of next atom */
00795       if (!strlen(trimleft(buffer))) break;
00796       
00797       /* Get shell type (s, p, d, f, g, sp) and # of primitives */
00798       sscanf(buffer,"%s %d %*f", shelltype, &numprims);
00799 
00800 
00801       /* Add new shell(s). */
00802       if (!strcasecmp(shelltype, "sp")) {
00803         /* Two new shells for SP */
00804         data->basis_set[i].shell =
00805           (shell_t*)realloc(data->basis_set[i].shell,
00806                             (numshells+2)*sizeof(shell_t));
00807       } else {
00808         /* One new shell for non-SP */
00809         data->basis_set[i].shell =
00810           (shell_t*)realloc(data->basis_set[i].shell,
00811                             (numshells+1)*sizeof(shell_t));
00812       }
00813 
00814       shell  = &(data->basis_set[i].shell[numshells]);
00815       memset(shell, 0, sizeof(shell_t));
00816       shell->numprims = numprims;
00817       shell->type     = shelltype_int(shelltype);
00818       shell->prim     = (prim_t*)calloc(numprims, sizeof(prim_t));
00819 
00820       /* If this is an SP-shell we have to add as separate 
00821        * S-shell and P-shell. */
00822       if (!strcasecmp(shelltype, "sp")) {
00823         shell->type      = SP_S_SHELL;
00824         shell2 = &(data->basis_set[i].shell[numshells+1]);
00825         shell2->numprims = numprims;
00826         shell2->type     = SP_P_SHELL;
00827         shell2->prim     = (prim_t*)calloc(numprims, sizeof(prim_t));
00828       }
00829 
00830       /* Loop over the primitives */
00831       for (j=0; j<numprims; j++) {
00832         int nr;
00833         double expon=0.f, coeff1, coeff2=0.f;
00834         if (!fgets(buffer,1024,data->file)) return FALSE;
00835 
00836         /* MOLDEN writes the basis set coefficients using Fortran style notation 
00837          * where the exponential character is 'D' instead of 'E'. Other packages 
00838          * adhere to C-style notation. Unfortunately sscanf() won't recognize 
00839          * Fortran-style numbers. Therefore we have to read the line as string 
00840          * first, convert the numbers by replacing the 'D' and then extract the 
00841          * floats using sscanf(). */
00842         fpexpftoc(buffer);
00843         nr = sscanf(buffer,"%lf %lf %lf", &expon, &coeff1, &coeff2);
00844         if (nr<2) {
00845           printf("moldenplugin) Bad format in [GTO] section\n");
00846           return FALSE;
00847         }
00848         shell->prim[j].exponent = expon;
00849         shell->prim[j].contraction_coeff = coeff1;
00850 
00851         /* P-shell component of SP-shell */
00852         if (!strcasecmp(shelltype, "sp")) {
00853           if (nr!=3) {
00854             printf("moldenplugin) Bad SP-shell format in [GTO] section\n");
00855             return FALSE;
00856           }
00857           shell2->prim[j].exponent = expon;
00858           shell2->prim[j].contraction_coeff = coeff2;
00859         }
00860       }
00861 
00862       /* Update # uncontracted basis functions */
00863       data->num_basis_funcs += numprims;
00864 
00865       numshells++;
00866 
00867       /* Account for SP-shells */
00868       if (!strcasecmp(shelltype, "sp")) {
00869         numshells++;
00870         data->num_basis_funcs += numprims;
00871       }
00872     }
00873 
00874     /* Store # shells for current atom */
00875     data->basis_set[i].numshells = numshells;
00876 
00877     /* Update total number of basis functions */
00878     data->num_shells += numshells;
00879   }
00880 
00881   /* As far as I know in Molden format the basis has
00882    * to be specified for each individual atom, even
00883    * if the types are the same. */
00884   data->num_basis_atoms = data->numatoms;
00885 
00886   /* allocate and populate flat arrays needed for molfileplugin */
00887   fill_basis_arrays(data);
00888 
00889   /* Count the number of cartesian basis functions in 
00890    * the basis set */
00891   data->wavef_size = 0;
00892   for (i=0; i<data->num_shells; i++) {
00893     switch (data->shell_types[i]) {
00894     case S_SHELL:
00895     case SP_S_SHELL:
00896       data->wavef_size += 1;
00897       break;
00898     case P_SHELL:
00899     case SP_P_SHELL:
00900       data->wavef_size += 3;
00901       break;
00902     case D_SHELL:
00903       data->wavef_size += 6;
00904       break;
00905     case F_SHELL:
00906       data->wavef_size += 10;
00907       break;
00908     case G_SHELL:
00909       data->wavef_size += 15;
00910       break;
00911     }
00912   }
00913 
00914   /* If we have MOs in the file we must provide the 
00915    * angular momentum exponents.
00916    * The order of P, D, F en G functions is as follows:
00917 
00918    *  5D: D 0, D+1, D-1, D+2, D-2
00919    *  6D: xx, yy, zz, xy, xz, yz
00920 
00921    *  7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
00922    * 10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
00923 
00924    *  9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
00925    * 15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
00926    *      xxyy xxzz yyzz xxyz yyxz zzxy
00927    */
00928   ALLOCATE(data->angular_momentum, int, 3*data->wavef_size);
00929 
00930   j=0;
00931   for (i=0; i<data->num_shells; i++) {
00932     switch (data->shell_types[i]) {
00933     case S_SHELL:
00934     case SP_S_SHELL:
00935       data->angular_momentum[j  ]=0;
00936       data->angular_momentum[j+1]=0;
00937       data->angular_momentum[j+2]=0;
00938       j += 3;
00939       break;
00940     case P_SHELL:
00941     case SP_P_SHELL:
00942       angular_momentum_expon(&data->angular_momentum[j  ], "x");
00943       angular_momentum_expon(&data->angular_momentum[j+3], "y");
00944       angular_momentum_expon(&data->angular_momentum[j+6], "z");
00945       j += 9;
00946       break;
00947     case D_SHELL:
00948       angular_momentum_expon(&data->angular_momentum[j   ], "xx");
00949       angular_momentum_expon(&data->angular_momentum[j+3 ], "yy");
00950       angular_momentum_expon(&data->angular_momentum[j+6 ], "zz");
00951       angular_momentum_expon(&data->angular_momentum[j+9 ], "xy");
00952       angular_momentum_expon(&data->angular_momentum[j+12], "xz");
00953       angular_momentum_expon(&data->angular_momentum[j+15], "yz");
00954       j += 18;
00955       break;
00956     case F_SHELL:
00957       angular_momentum_expon(&data->angular_momentum[j   ], "xxx");
00958       angular_momentum_expon(&data->angular_momentum[j+3 ], "yyy");
00959       angular_momentum_expon(&data->angular_momentum[j+6 ], "zzz");
00960       angular_momentum_expon(&data->angular_momentum[j+9 ], "xyy");
00961       angular_momentum_expon(&data->angular_momentum[j+12], "xxy");
00962       angular_momentum_expon(&data->angular_momentum[j+15], "xxz");
00963       angular_momentum_expon(&data->angular_momentum[j+18], "xzz");
00964       angular_momentum_expon(&data->angular_momentum[j+21], "yzz");
00965       angular_momentum_expon(&data->angular_momentum[j+24], "yyz");
00966       angular_momentum_expon(&data->angular_momentum[j+27], "xyz");
00967       j += 30;
00968       break;
00969     case G_SHELL:
00970       angular_momentum_expon(&data->angular_momentum[j   ], "xxxx");
00971       angular_momentum_expon(&data->angular_momentum[j+3 ], "yyyy");
00972       angular_momentum_expon(&data->angular_momentum[j+6 ], "zzzz");
00973       angular_momentum_expon(&data->angular_momentum[j+9 ], "xxxy");
00974       angular_momentum_expon(&data->angular_momentum[j+12], "xxxz");
00975       angular_momentum_expon(&data->angular_momentum[j+15], "yyyx");
00976       angular_momentum_expon(&data->angular_momentum[j+18], "yyyz");
00977       angular_momentum_expon(&data->angular_momentum[j+21], "zzzx");
00978       angular_momentum_expon(&data->angular_momentum[j+24], "zzzy");
00979       angular_momentum_expon(&data->angular_momentum[j+27], "xxyy");
00980       angular_momentum_expon(&data->angular_momentum[j+30], "xxzz");
00981       angular_momentum_expon(&data->angular_momentum[j+33], "yyzz");
00982       angular_momentum_expon(&data->angular_momentum[j+36], "xxyz");
00983       angular_momentum_expon(&data->angular_momentum[j+39], "yyxz");
00984       angular_momentum_expon(&data->angular_momentum[j+42], "zzxy");
00985       j += 45;
00986       break;
00987     }
00988   }
00989 
00990   return TRUE;
00991 }
00992 
00993 
00994 /******************************************************
00995  *
00996  * Populate the flat arrays containing the basis
00997  * set data.
00998  *
00999  ******************************************************/
01000 static int fill_basis_arrays(qmdata_t *data) {
01001   int i, j, k;
01002   int shellcount = 0;
01003   int primcount = 0;
01004 
01005   float *basis;
01006   int *num_shells_per_atom;
01007   int *num_prim_per_shell;
01008   int *shell_types;
01009   int *atomicnum_per_basisatom;
01010 
01011   /* Count the total number of primitives which
01012    * determines the size of the basis array. */
01013   for(i=0; i<data->num_basis_atoms; i++) {
01014     for (j=0; j<data->basis_set[i].numshells; j++) {
01015       primcount += data->basis_set[i].shell[j].numprims;
01016     }
01017   }
01018   data->num_basis_funcs = primcount;
01019 
01020   /* reserve space for pointer to array containing basis
01021    * info, i.e. contraction coeficients and expansion 
01022    * coefficients; need 2 entries per basis function, i.e.
01023    * exponent and contraction coefficient; also,
01024    * allocate space for the array holding the orbital symmetry
01025    * information per primitive Gaussian.
01026    * Finally, initialize the arrays holding the number of 
01027    * shells per atom and the number of primitives per shell*/
01028   ALLOCATE(basis,                   float, 2*primcount);
01029   ALLOCATE(shell_types,             int,   data->num_shells);
01030   ALLOCATE(num_shells_per_atom,     int,   data->num_basis_atoms);
01031   ALLOCATE(num_prim_per_shell,      int,   data->num_shells);
01032   ALLOCATE(atomicnum_per_basisatom, int,   data->num_basis_atoms);
01033 
01034 
01035 
01036   /* store pointers in struct qmdata_t */
01037   data->basis = basis;
01038   data->shell_types = shell_types;
01039   data->num_shells_per_atom = num_shells_per_atom;
01040   data->num_prim_per_shell  = num_prim_per_shell;
01041   data->atomicnum_per_basisatom = atomicnum_per_basisatom;
01042 
01043   primcount = 0;
01044   for (i=0; i<data->num_basis_atoms; i++) {
01045     /* assign atomic number */
01046     data->basis_set[i].atomicnum = data->atoms[i].atomicnum;
01047     atomicnum_per_basisatom[i]   = data->atoms[i].atomicnum;
01048 
01049     num_shells_per_atom[i] = data->basis_set[i].numshells;
01050 
01051     for (j=0; j<data->basis_set[i].numshells; j++) {
01052       shell_types[shellcount]        = data->basis_set[i].shell[j].type;
01053       num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
01054 
01055       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
01056         basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
01057         basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
01058         primcount++;
01059       }
01060       shellcount++;
01061     }
01062   } 
01063 
01064   return TRUE;
01065 }
01066 
01067 
01068 
01069 /**************************************************
01070  *
01071  * Convert shell type from char to int.
01072  * Note that SP_P shells are assigned in get_basis()
01073  *
01074  ************************************************ */
01075 static int shelltype_int(char *type) {
01076   int shelltype;
01077   if      (!strcasecmp(type, "sp")) shelltype = SP_SHELL;
01078   else if (!strcasecmp(type, "s"))  shelltype = S_SHELL;
01079   else if (!strcasecmp(type, "p"))  shelltype = P_SHELL;
01080   else if (!strcasecmp(type, "d"))  shelltype = D_SHELL;
01081   else if (!strcasecmp(type, "f"))  shelltype = F_SHELL;
01082   else if (!strcasecmp(type, "g"))  shelltype = G_SHELL;
01083   else shelltype = UNK_SHELL;
01084   
01085   return shelltype;
01086 }
01087 
01088 
01089 /***********************************************************
01090  *
01091  * Parse through the [MO] section and count orbitals and
01092  * the number of wavefunction coefficients per orbital.
01093  *
01094  * Format example of [MO] section:
01095  *  [MO]
01096  *  Sym=   1a         <-- this line is optional
01097  * Ene=   -15.5764
01098  *  Spin= Alpha
01099  * Occup=    2.00000
01100  *    1        -0.00000435
01101  *    2         0.00005919
01102  *    ...
01103  *
01104  ***********************************************************/
01105 static int count_orbitals(qmdata_t *data) {
01106   int nr;
01107   int num_wave_coeff=0;
01108   float orbenergy, occu;
01109   char spin[1024];
01110   qm_wavefunction_t *wave;
01111   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
01112   int dummy1;
01113   float dummy2;
01114 
01115   /* Place file pointer after [MO] keyword in line containing "Spin". */
01116   fseek(data->file, moldendata->filepos_mo, SEEK_SET);
01117   if (!goto_keyline(data->file, "Spin=", NULL)) {
01118     printf("moldenplugin) Couldn't find keyword 'Spin' in [MO] section!\n");
01119     return FALSE;
01120   }
01121  
01122   nr = fscanf(data->file, " Spin= %s\n", spin);
01123   eatline(data->file, 1);
01124 
01125   /* The first wavefunction should have spin alpha */
01126   strtoupper(spin);
01127   if (strcmp(spin, "ALPHA")) return FALSE;
01128 
01129   /* Removed redundant count of num_wave_coeff as this is equivalent to wavef_size */
01130   num_wave_coeff = data->wavef_size;
01131 
01132   /* For pruned AOs, with only non-zero coeffs, this is redundant */
01133   /*
01134   if (data->wavef_size && 
01135       data->wavef_size != num_wave_coeff) {
01136     printf("moldenplugin) Mismatch between # wavefunction coefficients (%d)\n",
01137            num_wave_coeff);
01138     printf("moldenplugin) and # cart. basis functions (%d)in basis set!\n", 
01139            data->wavef_size);
01140     return FALSE;
01141   }
01142   */
01143 
01144   /* Allocate memory for the qm_timestep frame */
01145   data->qm_timestep = (qm_timestep_t *)calloc(1, sizeof(qm_timestep_t));
01146 
01147   /* Add wavefunction for spin alpha */
01148   wave = add_wavefunction(data->qm_timestep);
01149 
01150   wave->spin = SPIN_ALPHA;
01151   wave->type = MOLFILE_WAVE_UNKNOWN;
01152   wave->exci = 0;
01153   wave->mult = 1;
01154   wave->num_coeffs = num_wave_coeff;
01155 
01156   /* Place file pointer on line after the [MO] keyword. */
01157   fseek(data->file, moldendata->filepos_mo, SEEK_SET);
01158 
01159   /* Count MOs */
01160   fscanf(data->file, " Ene= %f\n", &orbenergy);
01161   fscanf(data->file, " Spin= %s\n", spin);
01162   fscanf(data->file, " Occup= %f\n", &occu);
01163 
01164   while (1) {
01165     int check_reads = 2;
01166     wave->num_orbitals++;
01167 
01168     /* skip over MO coeffs */
01169     while(check_reads == 2)
01170     {
01171       check_reads = fscanf(data->file, "%d %f", &dummy1, &dummy2);
01172     }
01173 
01174     nr  = fscanf(data->file, " Ene= %f\n", &orbenergy);
01175     nr += fscanf(data->file, " Spin= %s\n", spin);
01176     nr += fscanf(data->file, " Occup= %f\n", &occu);
01177 
01178     if (nr!=3 || toupper(spin[0])!='A') 
01179       break;
01180   }
01181   //printf("found %d MOs!\n",wave->num_orbitals);
01182 
01183 
01184   /* Add wavefunction for spin beta */
01185   if (!strcmp(strtoupper(spin), "BETA")) {
01186     wave = add_wavefunction(data->qm_timestep);
01187     wave->spin = SPIN_BETA;
01188     wave->type = MOLFILE_WAVE_UNKNOWN;
01189     wave->exci = 0;
01190     wave->mult = 1;
01191     wave->num_coeffs = num_wave_coeff;
01192     wave->num_orbitals = 1;
01193 
01194     while (1) {
01195       int check_reads = 2;
01196       wave->num_orbitals++;
01197 
01198       /* skip over MO coeffs */
01199       while(check_reads == 2)
01200       {
01201         check_reads = fscanf(data->file, "%d %f", &dummy1, &dummy2);
01202       }
01203 
01204       nr  = fscanf(data->file, " Ene= %f\n", &orbenergy);
01205       nr += fscanf(data->file, " Spin= %s\n", spin);
01206       nr += fscanf(data->file, " Occup= %f\n", &occu);
01207 
01208       if (nr!=3 || toupper(spin[0])!='B' ||
01209           wave->num_orbitals>=num_wave_coeff) break;
01210     }
01211   }
01212 
01213   return TRUE;
01214 }
01215 
01216 
01217 static int read_molecular_orbitals(qmdata_t *data) {
01218   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
01219   qm_wavefunction_t *wave;
01220 
01221   if (!data->qm_timestep || moldendata->coordsonly) return FALSE;
01222 
01223   /* Place file pointer on line after the [MO] keyword. */
01224   fseek(data->file, moldendata->filepos_mo, SEEK_SET);
01225 
01226   wave = &data->qm_timestep->wave[0];
01227   ALLOCATE(wave->wave_coeffs, float, wave->num_coeffs*wave->num_orbitals);
01228   // DEBUG
01229   /*
01230   printf("moldenplugin) num_coeffs   = %d\n", wave->num_coeffs);
01231   printf("moldenplugin) num_orbitals = %d\n", wave->num_orbitals);
01232   printf("moldenplugin) num_wave     = %d\n", data->qm_timestep->numwave);
01233   */
01234   
01235 
01236   /* Read wavefunction coefficients for spin alpha */
01237   if (!read_wave_coeffs(data->file, wave)) return FALSE;
01238 
01239   if (data->qm_timestep->numwave==1) return TRUE;
01240 
01241   /* Read wavefunction coefficients for spin beta */
01242   wave = &data->qm_timestep->wave[1];
01243   ALLOCATE(wave->wave_coeffs, float, wave->num_coeffs*wave->num_orbitals);
01244 
01245   if (!read_wave_coeffs(data->file, wave)) return FALSE;
01246 
01247   return TRUE;
01248 }
01249 
01250 static int read_wave_coeffs(FILE *file, qm_wavefunction_t *wave) {
01251   int i, j, nr;
01252   char buffer[1024];
01253   int AOid;
01254   float wf_coeff;
01255   char keystring[10];
01256   float *wave_coeffs = wave->wave_coeffs;
01257 
01258   /* This works for pruned and unpruned (all zero coeffs are preserved) MO representaitions.
01259    * Clearly this data redundancy should be avoided on systems with a few hundred atoms. */
01260 
01261   /* set all coeffs to zero */
01262   for (i=0; i<wave->num_orbitals; i++)
01263     for (j=0; j<wave->num_coeffs; j++)
01264         wave_coeffs[i*wave->num_coeffs + j] = 0.0;
01265 
01266   /* each molecular orbital must have at least 1 non-zero coeff in its represenation */
01267   /* eat Ene= Spin= Occup=  lines */
01268   eatline(file, 3);
01269   for (i=0; i<wave->num_orbitals; i++) {
01270     while(1) {
01271       int nr2;
01272       if (!fgets(buffer,1024,file)) return FALSE;
01273       nr = sscanf(buffer,"%d %f", &AOid, &wf_coeff);
01274       wave_coeffs[i*wave->num_coeffs+AOid-1] = wf_coeff;
01275       
01276       // DEBUG
01277       //printf("moldenplugin) %d,%d: %d %f\n", i, AOid-1, AOid, wave_coeffs[i*wave->num_coeffs+AOid-1]);
01278       nr2 = sscanf(buffer, "%s", keystring);
01279       if(!strcmp(keystring,"Ene=")||nr2==-1) 
01280         break;
01281 
01282       if (nr==0) {
01283         printf("moldenplugin) Error reading wavefunction coefficients!\n");
01284         return FALSE;
01285       }
01286     }
01287     eatline(file, 2);
01288   }
01289   return TRUE;
01290 }
01291 
01292 /*************************************************************
01293  *
01294  * plugin registration 
01295  *
01296  **************************************************************/
01297 static molfile_plugin_t plugin;
01298 
01299 VMDPLUGIN_API int VMDPLUGIN_init() {
01300   memset(&plugin, 0, sizeof(molfile_plugin_t));
01301   plugin.abiversion = vmdplugin_ABIVERSION;
01302   plugin.type = MOLFILE_PLUGIN_TYPE;
01303   plugin.name = "molden";
01304   plugin.prettyname = "Molden";
01305   plugin.author = "Markus Dittrich, Jan Saam, Alexey Titov";
01306   plugin.majorv = 0;
01307   plugin.minorv = 10;
01308   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
01309   plugin.filename_extension = "molden";
01310   plugin.open_file_read = open_molden_read;
01311   plugin.read_structure = read_molden_structure;
01312 
01313   plugin.read_timestep_metadata    = read_timestep_metadata;
01314   plugin.read_timestep             = read_timestep;
01315   plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
01316 
01317   plugin.read_qm_metadata = read_molden_metadata;
01318   plugin.read_qm_rundata  = read_molden_rundata;
01319 
01320   plugin.close_file_read = close_molden_read;
01321   return VMDPLUGIN_SUCCESS;
01322 }
01323 
01324 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
01325   (*cb)(v, (vmdplugin_t *)&plugin);
01326   return VMDPLUGIN_SUCCESS;
01327 }
01328 
01329 VMDPLUGIN_API int VMDPLUGIN_fini() {
01330   return VMDPLUGIN_SUCCESS;
01331 }
01332 

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