1 /* MACHINE GENERATED FILE, DO NOT EDIT! */
2 
3 #define VMDPLUGIN molfile_moldenplugin
4 #define STATIC_PLUGIN 1
5 
6 /***************************************************************************
7  *cr
8  *cr            (C) Copyright 1995-2016 The Board of Trustees of the
9  *cr                        University of Illinois
10  *cr                         All Rights Reserved
11  *cr
12  ***************************************************************************/
13 
14 /***************************************************************************
15  * RCS INFORMATION:
16  *
17  *      $RCSfile: moldenplugin.c,v $
18  *      $Author: johns $       $Locker:  $             $State: Exp $
19  *      $Revision: 1.40 $       $Date: 2016/11/28 05:01:54 $
20  *
21  ***************************************************************************/
22 
23 /* This is a plugin that will read input from a MOLDEN
24 ** generated output file
25 ** some more details will go here soon
26 ** NOTE: The current version of the plugin relies
27 ** on the fact that the [Atom] field comes before
28 ** the [Geometries] field */
29 
30 
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <ctype.h>
34 #include <string.h>
35 
36 #include <string.h>
37 
38 #if defined(_AIX)
39 #include <strings.h>
40 #endif
41 
42 #if defined(WIN32) || defined(WIN64)
43 #define strcasecmp stricmp
44 #endif
45 
46 #include "molfile_plugin.h"
47 #include "unit_conversion.h"
48 #include "periodic_table.h"
49 #include "qmplugin.h"
50 
51 
52 #define ALLOCATE(array, type, size) \
53   array = (type *)calloc(size, sizeof(type)); \
54   if (array == NULL) { \
55     fprintf(stderr, "moldenplugin) Memory allocation for %s failed!\n", #array); \
56     return FALSE; \
57   }
58 
59 #define GET_LINE(x,y) if (!fgets(x, sizeof(x), y)) return FALSE
60 
61 
62 /* I could use a flag already present in qmdata_t to
63  * indicate a trajectory, but I'm using moldendata_t
64  * to demonstrate the use of
65  * void *format_specific_data;
66  * in qmdata_t as a means to store data specific to
67  * the plugin.
68  */
69 typedef struct {
70   long filepos_atoms;   /* [ATOMS] section */
71   long filepos_geomxyz; /* [GEOMETRIES] XYZ section */
72   long filepos_gto;     /* [GTO] section */
73   long filepos_mo;      /* [MO] section */
74   char units[16];
75   int coordsonly;
76 } moldendata_t;
77 
78 
79 
80 /* Read the basis set data */
81 static int get_basis (qmdata_t *);
82 static int shelltype_int(char *type);
83 static int fill_basis_arrays(qmdata_t *data);
84 
85 static int read_geom_block(qmdata_t *data);
86 static int read_molecular_orbitals(qmdata_t *data);
87 static int read_wave_coeffs(FILE *file, qm_wavefunction_t *wave);
88 static int count_orbitals(qmdata_t *data);
89 
90 
91 /*********************************************************
92  *
93  * Open file and return number of atoms.
94  *
95  * In order to determine the # atom we have to read into
96  * the [ATOMS]/[GEOMETRIES] sections where the atoms are
97  * defined.
98  * There can be either [ATOMS] or [GEOMETRIES] or both.
99  * While [GEOMETRIES] has the # atoms as its first line,
100  * we actually have to count lines in [ATOMS].
101  *
102  * We assume no particular order for the sections and scan
103  * the entire file for the according keywords. In order to
104  * save time when the section contents are actually read
105  * we store the file pointers to the beginning of each
106  * section in moldendata_t.
107  *
108  *********************************************************/
open_molden_read(const char * filename,const char * filetype,int * natoms)109 static void *open_molden_read(const char *filename,
110                               const char *filetype,
111                               int *natoms) {
112   FILE *fd;
113   qmdata_t *data = NULL;
114   moldendata_t *moldendata;
115   char buffer[1024];
116   char keystring[20];
117 
118   fd = fopen(filename, "rb");
119   if (!fd) return NULL;
120 
121   /* allocate memory for main QM data structure */
122   data = init_qmdata();
123   if (!data) return NULL;
124 
125   data->file = fd;
126 
127   /* allocate GAMESS specific data */
128   moldendata = (moldendata_t *)calloc(1, sizeof(moldendata_t));
129   if (!moldendata) return NULL;
130 
131   data->format_specific_data = moldendata;
132 
133 
134   /* Read first line */
135   if (!fgets(buffer,1024,data->file)) return NULL;
136 
137   /* Check if the file is MOLDEN format */
138   if (!strcmp(strtoupper(trimleft(trimright(buffer))), "[MOLDEN FORMAT]")) {
139     printf("moldenplugin) Detected MOLDEN file format!\n");
140   } else {
141     printf("moldenplugin) The file is not in MOLDEN format!\n");
142     return NULL;
143   }
144 
145   eatwhitelines(data->file);
146 
147 
148   /* Identify the different sections. */
149   while (fgets(buffer,1024,data->file)) {
150 
151     /* Get key string (ignoring empty lines) */
152     if (!sscanf(buffer, "%s", keystring)) continue;
153 
154     /* Quick test avoiding the uppercase transformation */
155     if (keystring[0]!='[') continue;
156 
157     /* Make keystring upper case */
158     strtoupper(keystring);
159 
160     if (!strcmp(keystring, "[5D]") || !strcmp(keystring, "[5D7F]") ||
161 	!strcmp(keystring, "[7F]") || !strcmp(keystring, "[5D10F]") ||
162 	!strcmp(keystring, "[9G]")) {
163       printf("moldenplugin) Spherical harmonic basis found %s. \n", keystring);
164       printf("moldenplugin)   Currently VMD handles only basis sets/wave functions\n");
165       printf("moldenplugin)   with cartesian Gaussian functions.\n");
166       printf("moldenplugin)   Loading coordinates only.\n");
167       moldendata->coordsonly = 1;
168     }
169 
170     if (!strcmp(keystring, "[ATOMS]")) {
171       char *s;
172       long prevline=ftell(fd);
173       printf("moldenplugin) Found [ATOMS] section ...\n");
174       moldendata->filepos_atoms = ftell(data->file);
175 
176       if (!sscanf(buffer, "%*s %s", moldendata->units)) {
177         printf("moldenplugin) Missing units in [ATOMS] section!\n");
178         return NULL;
179       }
180 
181       /* start counting the atoms;
182        * read until I hit the first line that starts with a "["
183        * bracket */
184       (*natoms) = 0;
185       s = fgets(buffer, 1024, fd);
186 
187       /* Here we assume that the [ATOMS] section goes
188        * on until the next empty line or another section
189        * starts, i.e. there is a "[" or we encounter EOF */
190       while (trimleft(buffer)[0]!='[' && s!=NULL && !iswhiteline(buffer)) {
191         (*natoms)++;
192         prevline = ftell(fd);
193         s = fgets(buffer, 1024, fd);
194       }
195       data->numatoms = *natoms;
196       data->num_frames = 1;
197 
198       /* Go back to the previous line, it might contain
199        * the next keyword */
200       fseek(data->file, prevline, SEEK_SET);
201     }
202 
203     else if (!strcmp(keystring, "[GEOMETRIES]")) {
204       if (!strcmp(trimright(buffer), "[GEOMETRIES] XYZ")) {
205         printf("moldenplugin) Found [GEOMETRIES] XYZ section ...\n");
206 
207         moldendata->filepos_geomxyz = ftell(data->file);
208 
209         /* The first line of the XYZ type [GEOMETRIES] input
210          * contains the number of atoms. */
211         if (fscanf(data->file, "%d", natoms) != 1) {
212           printf("moldenplugin) No # atoms found in [GEOMETRIES] section!\n");
213           return NULL;
214         }
215         data->numatoms = *natoms;
216 
217         /* Jump back to the beginning of the section */
218         fseek(data->file, moldendata->filepos_geomxyz, SEEK_SET);
219 
220         /* Count # frames */
221         data->num_frames = 0;
222         do {
223           int natm = 0;
224           fscanf(data->file, "%d", &natm);
225           if (natm!=data->numatoms) break;
226           eatline(data->file, 1);
227 
228           data->filepos_array = (long*)realloc(data->filepos_array,
229                                                (data->num_frames+1)*sizeof(long));
230           data->filepos_array[data->num_frames] = ftell(data->file);
231 
232           /* Skip title line + numatoms lines */
233           eatline(data->file, 1+data->numatoms);
234           if (feof(data->file)) break;
235 
236           data->num_frames++;
237         } while (1);
238 
239         printf("moldenplugin) Found %d frames\n", data->num_frames);
240       } else if (!strcmp(trimright(buffer), "[GEOMETRIES] ZMAT")) {
241         printf("moldenplugin) [GEOMETRIES] ZMAT not supported!\n");
242       }
243     }
244 
245     else if (!strcmp(keystring,"[GTO]")) {
246       printf("moldenplugin) Found [GTO] section ...\n");
247       moldendata->filepos_gto = ftell(data->file);
248     }
249 
250     else if (!strcmp(keystring,"[MO]")) {
251       printf("moldenplugin) Found [MO] section ...\n");
252       moldendata->filepos_mo = ftell(data->file);
253     }
254 
255   };
256 
257   return data;
258 }
259 
260 
261 
262 /*********************************************************
263  *
264  * Read geometry from file
265  *
266  * The [ATOMS] section provides atom name, atomic number
267  * and coordinates while the [GEOMETRIES] XYZ section
268  * provides atom type and coordinates. Trajectories can
269  * only be specified using [GEOMETRIES].
270  * In case we only have [GEOMETRIES] the atomic number
271  * will have to be deduced from the atom name.
272  *
273  *********************************************************/
read_molden_structure(void * mydata,int * optflags,molfile_atom_t * atoms)274 static int read_molden_structure(void *mydata, int *optflags,
275                                  molfile_atom_t *atoms)
276 {
277   int i;
278   char buffer[1024];
279   char atname[1024];
280   int num, atomicnum;
281   molfile_atom_t *atom;
282   qmdata_t *data = (qmdata_t *)mydata;
283   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
284 
285   ALLOCATE(data->atoms, qm_atom_t, data->numatoms);
286 
287   /* atomic number is provided by plugin.
288    * (This is required for QM plugins!) */
289   *optflags = MOLFILE_ATOMICNUMBER;
290 
291   /* [ATOMS] section */
292   if (moldendata->filepos_atoms) {
293     float unitfac = 1.f;
294 
295     /* If the units are given in AU we have to convert them.       */
296     /* Note: Also recognize parenthesized units emitted by Molcas. */
297     if (!strcmp(moldendata->units, "AU") ||
298         !strcmp(moldendata->units, "(AU)")) {
299       unitfac = BOHR_TO_ANGS;
300     }
301 
302     /* Jump to beginning of [ATOMS] section. */
303     fseek(data->file, moldendata->filepos_atoms, SEEK_SET);
304 
305     /* Read the atom types, names, atomic numbers
306      * as well as x,y,z coordinates */
307     for (i=0; i<data->numatoms; i++) {
308       float x,y,z;
309       atom = atoms+i;
310 
311       if (!fgets(buffer,1024,data->file)) return MOLFILE_ERROR;
312 
313       sscanf(buffer,"%s %d %d %f %f %f", atname, &num,
314              &atomicnum, &x, &y, &z);
315 
316       /* populate data structure for VMD */
317       strncpy(atom->name, atname,     sizeof(atom->name));
318       strncpy(atom->type, atom->name, sizeof(atom->type));
319       atom->atomicnumber = atomicnum;
320       atom->resname[0] = '\0';
321       atom->resid = 1;
322       atom->chain[0] = '\0';
323       atom->segid[0] = '\0';
324 
325       /* keep local copy */
326       strncpy(data->atoms[i].type, atname, sizeof(data->atoms[i].type));
327       data->atoms[i].atomicnum = atomicnum;
328       data->atoms[i].x = x*unitfac;
329       data->atoms[i].y = y*unitfac;
330       data->atoms[i].z = z*unitfac;
331     }
332     data->num_frames_read = 1;
333 
334     return MOLFILE_SUCCESS;
335   }
336 
337   /* [GEOMETRIES] XYZ section */
338   if (moldendata->filepos_geomxyz) {
339 
340     /* Jump to beginning of [GEOMETRIES] section. */
341     fseek(data->file, moldendata->filepos_geomxyz, SEEK_SET);
342     eatline(data->file, 2);
343 
344     /* Read block from file */
345     for (i=0; i<data->numatoms; i++) {
346       atom = atoms+i;
347       if (!fgets(buffer,1024,data->file)) return MOLFILE_ERROR;
348       sscanf(buffer,"%s %*f %*f %*f", atname);
349 
350       strncpy(atom->type, atname, sizeof(atom->type));
351       strncpy(atom->name, atname, sizeof(atom->name));
352       atom->atomicnumber = get_pte_idx_from_string(atname);
353       atom->resname[0] = '\0';
354       atom->resid = 1;
355       atom->chain[0] = '\0';
356       atom->segid[0] = '\0';
357       data->atoms[i].atomicnum = atom->atomicnumber;
358     }
359     data->num_frames_read = 0;
360 
361     return MOLFILE_SUCCESS;
362   }
363 
364   printf("Sorry, could not obtain structure information \n");
365   printf("from either the [ATOMS] or [GEOMETRIES] section! \n");
366   printf("Please check your MOLDEN output file! \n");
367   return MOLFILE_ERROR;
368 }
369 
370 
371 /***********************************************************
372  *
373  * Read atoms for one frame from [GEOMETRIES] section.
374  *
375  ***********************************************************/
read_geom_block(qmdata_t * data)376 static int read_geom_block(qmdata_t *data) {
377   int i;
378   char buffer[1024];
379   float x,y,z;
380 
381   /* Skip title line */
382   eatline(data->file, 1);
383 
384   for (i=0; i<data->numatoms; i++) {
385     if (!fgets(buffer,1024,data->file)) return 0;
386     sscanf(buffer,"%*s %f %f %f", &x, &y, &z);
387     data->atoms[i].x = x;
388     data->atoms[i].y = y;
389     data->atoms[i].z = z;
390   }
391 
392   return 1;
393 }
394 
395 
396 /***********************************************************
397  *
398  * Provide non-QM metadata for next timestep.
399  * Required by the plugin interface.
400  *
401  ***********************************************************/
read_timestep_metadata(void * mydata,molfile_timestep_metadata_t * meta)402 static int read_timestep_metadata(void *mydata,
403                                   molfile_timestep_metadata_t *meta) {
404   meta->count = -1;
405   meta->has_velocities = 0;
406 
407   return MOLFILE_SUCCESS;
408 }
409 
410 
411 /***********************************************************
412  *
413  * We are not reading the coefficients themselves,
414  * because that could require a large amount of memory.
415  *
416  ***********************************************************/
read_qm_timestep_metadata(void * mydata,molfile_qm_timestep_metadata_t * meta)417 static int read_qm_timestep_metadata(void *mydata,
418                                     molfile_qm_timestep_metadata_t *meta) {
419   qmdata_t *data = (qmdata_t *)mydata;
420   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
421 
422   if (data->num_frames_sent >= data->num_frames) {
423     /* All frames were sent. */
424     return MOLFILE_ERROR;
425   }
426 
427   /* Can't send metadata if only coordinates were read */
428   if (moldendata->coordsonly) return MOLFILE_ERROR;
429 
430   /* Count the number of cartesian basis functions in
431      the basis set */
432   if (data->num_frames_sent == data->num_frames-1) {
433     int i;
434     qm_timestep_t *cur_ts;
435 
436     if (!count_orbitals(data)) return MOLFILE_ERROR;
437 
438     /* get a pointer to the current qm timestep */
439     cur_ts = data->qm_timestep;
440 
441     for (i=0; (i<MOLFILE_MAXWAVEPERTS && i<cur_ts->numwave); i++) {
442       meta->num_orbitals_per_wavef[i] = cur_ts->wave[i].num_orbitals;
443       meta->has_occup_per_wavef[i]    = cur_ts->wave[i].has_occup;
444       meta->has_orben_per_wavef[i]    = cur_ts->wave[i].has_orben;
445     }
446     meta->wavef_size   = data->wavef_size;
447     meta->num_wavef    = cur_ts->numwave;
448     meta->num_scfiter  = cur_ts->num_scfiter;
449     meta->has_gradient = FALSE;
450     meta->num_charge_sets = 0;
451   }
452 
453   return MOLFILE_SUCCESS;
454 }
455 
456 
457 
458 /***********************************************************
459  *
460  * Provides VMD with the data of the next timestep.
461  *
462  ***********************************************************/
read_timestep(void * mydata,int natoms,molfile_timestep_t * ts,molfile_qm_metadata_t * qm_metadata,molfile_qm_timestep_t * qm_ts)463 static int read_timestep(void *mydata, int natoms,
464        molfile_timestep_t *ts, molfile_qm_metadata_t *qm_metadata,
465                          molfile_qm_timestep_t *qm_ts) {
466   int i;
467   qmdata_t *data = (qmdata_t *)mydata;
468   qm_timestep_t *cur_ts;
469 
470   if (data->num_frames_sent >= data->num_frames) {
471     /* All frames were sent. */
472     return MOLFILE_ERROR;
473   }
474 
475   if (data->num_frames_sent == data->num_frames_read) {
476     /* Read next coordinate block from file */
477     fseek(data->file, data->filepos_array[data->num_frames_read], SEEK_SET);
478     read_geom_block(data);
479 
480     /*printf("moldenplugin) Read frame %d\n", data->num_frames_read); */
481     data->num_frames_read++;
482   }
483 
484 
485   /* Copy the coordinates */
486   for (i=0; i<natoms; i++) {
487     ts->coords[3*i  ] = data->atoms[i].x;
488     ts->coords[3*i+1] = data->atoms[i].y;
489     ts->coords[3*i+2] = data->atoms[i].z;
490   }
491 
492   /*printf("moldenplugin) Sent frame %d\n", data->num_frames_sent); */
493   data->num_frames_sent++;
494 
495   /* In MOLDEN the MOs are listed only for the last frame */
496   if (data->num_frames_sent == data->num_frames) {
497 
498     /* get a convenient pointer to the current qm timestep */
499     cur_ts = data->qm_timestep;
500 
501     read_molecular_orbitals(data);
502 
503     /* store the wave function and orbital energies */
504     if (cur_ts != NULL && cur_ts->wave != NULL) {
505       for (i=0; i<cur_ts->numwave; i++) {
506         qm_wavefunction_t *wave = &cur_ts->wave[i];
507         qm_ts->wave[i].type         = wave->type;
508         qm_ts->wave[i].spin         = wave->spin;
509         qm_ts->wave[i].excitation   = wave->exci;
510         qm_ts->wave[i].multiplicity = wave->mult;
511         qm_ts->wave[i].energy       = wave->energy;
512         strncpy(qm_ts->wave[i].info, wave->info, MOLFILE_BUFSIZ);
513 
514         if (wave->wave_coeffs) {
515           memcpy(qm_ts->wave[i].wave_coeffs, wave->wave_coeffs,
516                  wave->num_orbitals*data->wavef_size*sizeof(float));
517         }
518         if (wave->orb_energies) {
519           memcpy(qm_ts->wave[i].orbital_energies, wave->orb_energies,
520                  wave->num_orbitals*sizeof(float));
521         }
522         if (wave->has_occup) {
523           memcpy(qm_ts->wave[i].occupancies, wave->orb_occupancies,
524                  wave->num_orbitals*sizeof(float));
525         }
526       }
527     }
528 
529   }
530 
531   return MOLFILE_SUCCESS;
532 }
533 
534 
535 /*****************************************************
536  *
537  * Provide VMD with the sizes of the QM related
538  * data structure arrays that need to be made
539  * available.
540  * Since we cannot determine the basis set meta data
541  * without parsing the whole basis set section, we
542  * read all basis set data here. The data is stored
543  * in the qmdata_t structure for later use in
544  * read_molden_rundata().
545  *
546  *****************************************************/
read_molden_metadata(void * mydata,molfile_qm_metadata_t * metadata)547 static int read_molden_metadata(void *mydata,
548     molfile_qm_metadata_t *metadata) {
549 
550   qmdata_t *data;
551   moldendata_t *moldendata;
552   data = (qmdata_t *)mydata;
553   moldendata = (moldendata_t *)data->format_specific_data;
554 
555 
556   metadata->ncart = 0;
557   metadata->nimag = 0;
558   metadata->nintcoords = 0;
559 
560   metadata->have_sysinfo = 0;
561   metadata->have_carthessian = 0;
562   metadata->have_inthessian = 0;
563   metadata->have_normalmodes = 0;
564 
565   metadata->num_basis_funcs = 0;
566   metadata->num_basis_atoms = 0;
567   metadata->num_shells = 0;
568   metadata->wavef_size = 0;
569 
570   if (!moldendata->coordsonly) {
571     /* Read the basis set */
572     if (!get_basis(data)) return MOLFILE_ERROR;
573 
574     /* orbital + basis set data */
575     metadata->num_basis_funcs = data->num_basis_funcs;
576     metadata->num_basis_atoms = data->num_basis_atoms;
577     metadata->num_shells      = data->num_shells;
578     metadata->wavef_size      = data->wavef_size;
579   }
580 
581   return MOLFILE_SUCCESS;
582 }
583 
584 
585 /******************************************************
586  *
587  * Provide VMD with the static (i.e. non-trajectory)
588  * data. That means we are filling the molfile_plugin
589  * data structures.
590  *
591  ******************************************************/
read_molden_rundata(void * mydata,molfile_qm_t * qm_data)592 static int read_molden_rundata(void *mydata,
593                                molfile_qm_t *qm_data) {
594   qmdata_t *data = (qmdata_t *)mydata;
595   int i;
596   molfile_qm_hessian_t *hessian_data;
597   molfile_qm_basis_t   *basis_data;
598   molfile_qm_sysinfo_t *sys_data;
599 
600   if (!qm_data) return MOLFILE_ERROR;
601 
602 
603   hessian_data = &qm_data->hess;
604   basis_data   = &qm_data->basis;
605   sys_data     = &qm_data->run;
606 
607   sys_data->num_electrons = data->num_electrons;
608   sys_data->totalcharge = data->totalcharge;
609 
610   /* Populate basis set data */
611   if (data->num_basis_funcs) {
612     for (i=0; i<data->num_basis_atoms; i++) {
613       basis_data->num_shells_per_atom[i] = data->num_shells_per_atom[i];
614       basis_data->atomic_number[i] = data->atomicnum_per_basisatom[i];
615     }
616 
617     for (i=0; i<data->num_shells; i++) {
618       basis_data->num_prim_per_shell[i] = data->num_prim_per_shell[i];
619       basis_data->shell_types[i] = data->shell_types[i];
620     }
621 
622     for (i=0; i<2*data->num_basis_funcs; i++) {
623       basis_data->basis[i] = data->basis[i];
624     }
625 
626     /* If we have MOs in the file we must provide the
627      * angular momentum exponents */
628     if (data->angular_momentum) {
629       for (i=0; i<3*data->wavef_size; i++) {
630         basis_data->angular_momentum[i] = data->angular_momentum[i];
631       }
632     }
633   }
634 
635   /* fill in molfile_qm_sysinfo_t */
636   /*sys_data->runtype = data->runtype;
637   sys_data->scftype = data->scftype;
638   sys_data->nproc   = data->nproc;
639   sys_data->num_electrons  = data->num_electrons;
640   sys_data->totalcharge    = data->totalcharge;
641   sys_data->num_occupied_A = data->num_occupied_A;
642   sys_data->num_occupied_B = data->num_occupied_B;
643   sys_data->status         = data->opt_status;
644   */
645   return MOLFILE_SUCCESS;
646 }
647 
648 
649 /**********************************************************
650  *
651  * close file and free memory
652  *
653  **********************************************************/
close_molden_read(void * mydata)654 static void close_molden_read(void *mydata) {
655   int i, j;
656   qmdata_t *data = (qmdata_t *)mydata;
657 
658   fclose(data->file);
659 
660   free(data->atoms);
661   free(data->basis);
662   free(data->shell_types);
663   free(data->atomicnum_per_basisatom);
664   free(data->num_shells_per_atom);
665   free(data->num_prim_per_shell);
666   free(data->angular_momentum);
667 
668   if (data->basis_set) {
669     for(i=0; i<data->num_basis_atoms; i++) {
670       for (j=0; j<data->basis_set[i].numshells; j++) {
671         free(data->basis_set[i].shell[j].prim);
672       }
673       free(data->basis_set[i].shell);
674     }
675     free(data->basis_set);
676   }
677 
678   free(data->format_specific_data);
679   free(data->filepos_array);
680 
681   if (data->qm_timestep != NULL) {
682     for (j=0; j<data->qm_timestep[0].numwave; j++) {
683       free(data->qm_timestep[0].wave[j].wave_coeffs);
684       free(data->qm_timestep[0].wave[j].orb_energies);
685       free(data->qm_timestep[0].wave[j].orb_occupancies);
686     }
687     free(data->qm_timestep[0].wave);
688     free(data->qm_timestep);
689   } else {
690     printf("close_molden_read(): NULL qm_timestep!\n");
691   }
692 
693   free(data);
694 }
695 
696 
697 
698 /* ####################################################### */
699 /*             End of API functions                        */
700 /* The following functions actually do the file parsing.   */
701 /* ####################################################### */
702 
703 
704 /******************************************************
705  *
706  * Format specification of the basis-set consisting of
707  * contracted Gaussian Type Orbitals.
708  *
709  * [GTO]
710  * atom_sequence_number1 0
711  * shell_label number_of_primitives 1.00
712  * exponent_prim_1 contraction_coeff_1 (contraction_coeff_sp_1)
713  * ...
714  * <empty line>
715  * atom_sequence_number2 0
716  * shell_label number_of_primitives 1.00
717  * exponent_prim_1 contraction_coeff_1 (contraction_coeff_1)
718  * ...
719  * <empty line>
720  *
721  * recognized shell_labels: s, p, d, f, sp, g
722  *
723  * For 'sp' shells two contraction coefficients must be given,
724  * for both s and p functions.
725  * The 0 on the shell_number line and the 1.00 on the
726  * shell_label line are no longer functional and can be
727  * ignored.
728  *
729  * All workings with the [GTO] keyword are in Atomic Units.
730  *
731  ******************************************************/
732 
733 
734 /*******************************************************
735  *
736  * Read the basis set data
737  *
738  * Format example:
739  * [GTO]
740  *   1 0
741  *  s    2 1.00
742  *   0.2738503300E+02  0.4301284983E+00
743  *   0.4874522100E+01  0.6789135305E+00
744  *  sp   2 1.00
745  *   0.1136748200E+01  0.4947176920E-01  0.5115407076E+00
746  *   0.2883094000E+00  0.9637824081E+00  0.6128198961E+00
747  *
748  *   2 0
749  *  s    2 1.00
750  *   0.1309756400E+01  0.4301284983E+00
751  *   0.2331360000E+00  0.6789135305E+00
752  *  ...
753  *
754  * qmdata_t provides hierarchical data structures for
755  * the basis set which are convenient for parsing.
756  * The molfile_plugin interface, however, requires flat
757  * arrays, so after reading is done we have to populate
758  * the according arrays using fill_basis_arrays().
759  *
760  *******************************************************/
get_basis(qmdata_t * data)761 static int get_basis(qmdata_t *data) {
762   char buffer[1024];
763   char shelltype[1024];
764   int atomid, numprims;
765   int i, j=0;
766   int numshells;
767   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
768 
769   /* XXX already initialized in open_molden_read() */
770   data->num_shells = 0;
771   data->num_basis_funcs = 0;
772   data->num_basis_atoms = 0;
773 
774   /* initialize basis set the character array */
775   memset(data->basis_string, 0, sizeof(data->basis_string));
776 
777 
778   /* Place file pointer on line after the [GTO] keyword. */
779   fseek(data->file, moldendata->filepos_gto, SEEK_SET);
780 
781   /* Allocate memory for the basis of all atoms. */
782   ALLOCATE(data->basis_set, basis_atom_t, data->numatoms);
783 
784   /* Loop over all atoms. */
785   for (i=0; i<data->numatoms; i++) {
786 
787     if (!fgets(buffer,1024,data->file)) return FALSE;
788     sscanf(buffer,"%d %*d", &atomid);
789 
790     numshells = 0;
791     data->basis_set[i].shell = NULL;
792 
793     /* Read an unknown number of shells */
794     while (1) {
795       shell_t *shell, *shell2=NULL;
796 
797       if (!fgets(buffer,1024,data->file)) return FALSE;
798 
799       /* Empty line signifies beginning of next atom */
800       if (!strlen(trimleft(buffer))) break;
801 
802       /* Get shell type (s, p, d, f, g, sp) and # of primitives */
803       sscanf(buffer,"%s %d %*f", shelltype, &numprims);
804 
805 
806       /* Add new shell(s). */
807       if (!strcasecmp(shelltype, "sp")) {
808         /* Two new shells for SP */
809         data->basis_set[i].shell =
810           (shell_t*)realloc(data->basis_set[i].shell,
811                             (numshells+2)*sizeof(shell_t));
812       } else {
813         /* One new shell for non-SP */
814         data->basis_set[i].shell =
815           (shell_t*)realloc(data->basis_set[i].shell,
816                             (numshells+1)*sizeof(shell_t));
817       }
818 
819       shell  = &(data->basis_set[i].shell[numshells]);
820       memset(shell, 0, sizeof(shell_t));
821       shell->numprims = numprims;
822       shell->type     = shelltype_int(shelltype);
823       shell->prim     = (prim_t*)calloc(numprims, sizeof(prim_t));
824 
825       /* If this is an SP-shell we have to add as separate
826        * S-shell and P-shell. */
827       if (!strcasecmp(shelltype, "sp")) {
828         shell->type      = SP_S_SHELL;
829         shell2 = &(data->basis_set[i].shell[numshells+1]);
830         shell2->numprims = numprims;
831         shell2->type     = SP_P_SHELL;
832         shell2->prim     = (prim_t*)calloc(numprims, sizeof(prim_t));
833       }
834 
835       /* Loop over the primitives */
836       for (j=0; j<numprims; j++) {
837         int nr;
838         double expon=0.f, coeff1, coeff2=0.f;
839         if (!fgets(buffer,1024,data->file)) return FALSE;
840 
841 	/* MOLDEN writes the basis set coefficients using Fortran style notation
842 	 * where the exponential character is 'D' instead of 'E'. Other packages
843 	 * adhere to C-style notation. Unfortunately sscanf() won't recognize
844 	 * Fortran-style numbers. Therefore we have to read the line as string
845 	 * first, convert the numbers by replacing the 'D' and then extract the
846 	 * floats using sscanf(). */
847 	fpexpftoc(buffer);
848 	nr = sscanf(buffer,"%lf %lf %lf", &expon, &coeff1, &coeff2);
849         if (nr<2) {
850           printf("moldenplugin) Bad format in [GTO] section\n");
851           return FALSE;
852         }
853         shell->prim[j].exponent = expon;
854         shell->prim[j].contraction_coeff = coeff1;
855 
856         /* P-shell component of SP-shell */
857         if (!strcasecmp(shelltype, "sp")) {
858           if (nr!=3) {
859             printf("moldenplugin) Bad SP-shell format in [GTO] section\n");
860             return FALSE;
861           }
862           shell2->prim[j].exponent = expon;
863           shell2->prim[j].contraction_coeff = coeff2;
864         }
865       }
866 
867       /* Update # uncontracted basis functions */
868       data->num_basis_funcs += numprims;
869 
870       numshells++;
871 
872       /* Account for SP-shells */
873       if (!strcasecmp(shelltype, "sp")) {
874         numshells++;
875         data->num_basis_funcs += numprims;
876       }
877     }
878 
879     /* Store # shells for current atom */
880     data->basis_set[i].numshells = numshells;
881 
882     /* Update total number of basis functions */
883     data->num_shells += numshells;
884   }
885 
886   /* As far as I know in Molden format the basis has
887    * to be specified for each individual atom, even
888    * if the types are the same. */
889   data->num_basis_atoms = data->numatoms;
890 
891   /* allocate and populate flat arrays needed for molfileplugin */
892   fill_basis_arrays(data);
893 
894   /* Count the number of cartesian basis functions in
895    * the basis set */
896   data->wavef_size = 0;
897   for (i=0; i<data->num_shells; i++) {
898     switch (data->shell_types[i]) {
899     case S_SHELL:
900     case SP_S_SHELL:
901       data->wavef_size += 1;
902       break;
903     case P_SHELL:
904     case SP_P_SHELL:
905       data->wavef_size += 3;
906       break;
907     case D_SHELL:
908       data->wavef_size += 6;
909       break;
910     case F_SHELL:
911       data->wavef_size += 10;
912       break;
913     case G_SHELL:
914       data->wavef_size += 15;
915       break;
916     }
917   }
918 
919   /* If we have MOs in the file we must provide the
920    * angular momentum exponents.
921    * The order of P, D, F en G functions is as follows:
922 
923    *  5D: D 0, D+1, D-1, D+2, D-2
924    *  6D: xx, yy, zz, xy, xz, yz
925 
926    *  7F: F 0, F+1, F-1, F+2, F-2, F+3, F-3
927    * 10F: xxx, yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
928 
929    *  9G: G 0, G+1, G-1, G+2, G-2, G+3, G-3, G+4, G-4
930    * 15G: xxxx yyyy zzzz xxxy xxxz yyyx yyyz zzzx zzzy,
931    *      xxyy xxzz yyzz xxyz yyxz zzxy
932    */
933   ALLOCATE(data->angular_momentum, int, 3*data->wavef_size);
934 
935   j=0;
936   for (i=0; i<data->num_shells; i++) {
937     switch (data->shell_types[i]) {
938     case S_SHELL:
939     case SP_S_SHELL:
940       data->angular_momentum[j  ]=0;
941       data->angular_momentum[j+1]=0;
942       data->angular_momentum[j+2]=0;
943       j += 3;
944       break;
945     case P_SHELL:
946     case SP_P_SHELL:
947       angular_momentum_expon(&data->angular_momentum[j  ], "x");
948       angular_momentum_expon(&data->angular_momentum[j+3], "y");
949       angular_momentum_expon(&data->angular_momentum[j+6], "z");
950       j += 9;
951       break;
952     case D_SHELL:
953       angular_momentum_expon(&data->angular_momentum[j   ], "xx");
954       angular_momentum_expon(&data->angular_momentum[j+3 ], "yy");
955       angular_momentum_expon(&data->angular_momentum[j+6 ], "zz");
956       angular_momentum_expon(&data->angular_momentum[j+9 ], "xy");
957       angular_momentum_expon(&data->angular_momentum[j+12], "xz");
958       angular_momentum_expon(&data->angular_momentum[j+15], "yz");
959       j += 18;
960       break;
961     case F_SHELL:
962       angular_momentum_expon(&data->angular_momentum[j   ], "xxx");
963       angular_momentum_expon(&data->angular_momentum[j+3 ], "yyy");
964       angular_momentum_expon(&data->angular_momentum[j+6 ], "zzz");
965       angular_momentum_expon(&data->angular_momentum[j+9 ], "xyy");
966       angular_momentum_expon(&data->angular_momentum[j+12], "xxy");
967       angular_momentum_expon(&data->angular_momentum[j+15], "xxz");
968       angular_momentum_expon(&data->angular_momentum[j+18], "xzz");
969       angular_momentum_expon(&data->angular_momentum[j+21], "yzz");
970       angular_momentum_expon(&data->angular_momentum[j+24], "yyz");
971       angular_momentum_expon(&data->angular_momentum[j+27], "xyz");
972       j += 30;
973       break;
974     case G_SHELL:
975       angular_momentum_expon(&data->angular_momentum[j   ], "xxxx");
976       angular_momentum_expon(&data->angular_momentum[j+3 ], "yyyy");
977       angular_momentum_expon(&data->angular_momentum[j+6 ], "zzzz");
978       angular_momentum_expon(&data->angular_momentum[j+9 ], "xxxy");
979       angular_momentum_expon(&data->angular_momentum[j+12], "xxxz");
980       angular_momentum_expon(&data->angular_momentum[j+15], "yyyx");
981       angular_momentum_expon(&data->angular_momentum[j+18], "yyyz");
982       angular_momentum_expon(&data->angular_momentum[j+21], "zzzx");
983       angular_momentum_expon(&data->angular_momentum[j+24], "zzzy");
984       angular_momentum_expon(&data->angular_momentum[j+27], "xxyy");
985       angular_momentum_expon(&data->angular_momentum[j+30], "xxzz");
986       angular_momentum_expon(&data->angular_momentum[j+33], "yyzz");
987       angular_momentum_expon(&data->angular_momentum[j+36], "xxyz");
988       angular_momentum_expon(&data->angular_momentum[j+39], "yyxz");
989       angular_momentum_expon(&data->angular_momentum[j+42], "zzxy");
990       j += 45;
991       break;
992     }
993   }
994 
995   return TRUE;
996 }
997 
998 
999 /******************************************************
1000  *
1001  * Populate the flat arrays containing the basis
1002  * set data.
1003  *
1004  ******************************************************/
fill_basis_arrays(qmdata_t * data)1005 static int fill_basis_arrays(qmdata_t *data) {
1006   int i, j, k;
1007   int shellcount = 0;
1008   int primcount = 0;
1009 
1010   float *basis;
1011   int *num_shells_per_atom;
1012   int *num_prim_per_shell;
1013   int *shell_types;
1014   int *atomicnum_per_basisatom;
1015 
1016   /* Count the total number of primitives which
1017    * determines the size of the basis array. */
1018   for(i=0; i<data->num_basis_atoms; i++) {
1019     for (j=0; j<data->basis_set[i].numshells; j++) {
1020       primcount += data->basis_set[i].shell[j].numprims;
1021     }
1022   }
1023   data->num_basis_funcs = primcount;
1024 
1025   /* reserve space for pointer to array containing basis
1026    * info, i.e. contraction coeficients and expansion
1027    * coefficients; need 2 entries per basis function, i.e.
1028    * exponent and contraction coefficient; also,
1029    * allocate space for the array holding the orbital symmetry
1030    * information per primitive Gaussian.
1031    * Finally, initialize the arrays holding the number of
1032    * shells per atom and the number of primitives per shell*/
1033   ALLOCATE(basis,                   float, 2*primcount);
1034   ALLOCATE(shell_types,             int,   data->num_shells);
1035   ALLOCATE(num_shells_per_atom,     int,   data->num_basis_atoms);
1036   ALLOCATE(num_prim_per_shell,      int,   data->num_shells);
1037   ALLOCATE(atomicnum_per_basisatom, int,   data->num_basis_atoms);
1038 
1039 
1040 
1041   /* store pointers in struct qmdata_t */
1042   data->basis = basis;
1043   data->shell_types = shell_types;
1044   data->num_shells_per_atom = num_shells_per_atom;
1045   data->num_prim_per_shell  = num_prim_per_shell;
1046   data->atomicnum_per_basisatom = atomicnum_per_basisatom;
1047 
1048   primcount = 0;
1049   for (i=0; i<data->num_basis_atoms; i++) {
1050     /* assign atomic number */
1051     data->basis_set[i].atomicnum = data->atoms[i].atomicnum;
1052     atomicnum_per_basisatom[i]   = data->atoms[i].atomicnum;
1053 
1054     num_shells_per_atom[i] = data->basis_set[i].numshells;
1055 
1056     for (j=0; j<data->basis_set[i].numshells; j++) {
1057       shell_types[shellcount]        = data->basis_set[i].shell[j].type;
1058       num_prim_per_shell[shellcount] = data->basis_set[i].shell[j].numprims;
1059 
1060       for (k=0; k<data->basis_set[i].shell[j].numprims; k++) {
1061         basis[2*primcount  ] = data->basis_set[i].shell[j].prim[k].exponent;
1062         basis[2*primcount+1] = data->basis_set[i].shell[j].prim[k].contraction_coeff;
1063         primcount++;
1064       }
1065       shellcount++;
1066     }
1067   }
1068 
1069   return TRUE;
1070 }
1071 
1072 
1073 
1074 /**************************************************
1075  *
1076  * Convert shell type from char to int.
1077  * Note that SP_P shells are assigned in get_basis()
1078  *
1079  ************************************************ */
shelltype_int(char * type)1080 static int shelltype_int(char *type) {
1081   int shelltype;
1082   if      (!strcasecmp(type, "sp")) shelltype = SP_SHELL;
1083   else if (!strcasecmp(type, "s"))  shelltype = S_SHELL;
1084   else if (!strcasecmp(type, "p"))  shelltype = P_SHELL;
1085   else if (!strcasecmp(type, "d"))  shelltype = D_SHELL;
1086   else if (!strcasecmp(type, "f"))  shelltype = F_SHELL;
1087   else if (!strcasecmp(type, "g"))  shelltype = G_SHELL;
1088   else shelltype = UNK_SHELL;
1089 
1090   return shelltype;
1091 }
1092 
1093 
1094 /***********************************************************
1095  *
1096  * Parse through the [MO] section and count orbitals and
1097  * the number of wavefunction coefficients per orbital.
1098  *
1099  * Format example of [MO] section:
1100  *  [MO]
1101  *  Sym=   1a         <-- this line is optional
1102  * Ene=   -15.5764
1103  *  Spin= Alpha
1104  * Occup=    2.00000
1105  *    1        -0.00000435
1106  *    2         0.00005919
1107  *    ...
1108  *
1109  ***********************************************************/
count_orbitals(qmdata_t * data)1110 static int count_orbitals(qmdata_t *data) {
1111   int nr;
1112   int num_wave_coeff=0;
1113   float orbenergy, occu;
1114   char spin[1024];
1115   qm_wavefunction_t *wave;
1116   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
1117   int dummy1;
1118   float dummy2;
1119 
1120   /* Place file pointer after [MO] keyword in line containing "Spin". */
1121   fseek(data->file, moldendata->filepos_mo, SEEK_SET);
1122   if (!goto_keyline(data->file, "Spin=", NULL)) {
1123     printf("moldenplugin) Couldn't find keyword 'Spin' in [MO] section!\n");
1124     return FALSE;
1125   }
1126 
1127   nr = fscanf(data->file, " Spin= %s\n", spin);
1128   eatline(data->file, 1);
1129 
1130   /* The first wavefunction should have spin alpha */
1131   strtoupper(spin);
1132   if (strcmp(spin, "ALPHA")) return FALSE;
1133 
1134   /* Removed redundant count of num_wave_coeff as this is equivalent to wavef_size */
1135   num_wave_coeff = data->wavef_size;
1136 
1137   /* For pruned AOs, with only non-zero coeffs, this is redundant */
1138   /*
1139   if (data->wavef_size &&
1140       data->wavef_size != num_wave_coeff) {
1141     printf("moldenplugin) Mismatch between # wavefunction coefficients (%d)\n",
1142            num_wave_coeff);
1143     printf("moldenplugin) and # cart. basis functions (%d)in basis set!\n",
1144            data->wavef_size);
1145     return FALSE;
1146   }
1147   */
1148 
1149   /* Allocate memory for the qm_timestep frame */
1150   data->qm_timestep = (qm_timestep_t *)calloc(1, sizeof(qm_timestep_t));
1151 
1152   /* Add wavefunction for spin alpha */
1153   wave = add_wavefunction(data->qm_timestep);
1154 
1155   wave->spin = SPIN_ALPHA;
1156   wave->type = MOLFILE_WAVE_UNKNOWN;
1157   wave->exci = 0;
1158   wave->mult = 1;
1159   wave->num_coeffs = num_wave_coeff;
1160 
1161   /* Place file pointer on line after the [MO] keyword. */
1162   fseek(data->file, moldendata->filepos_mo, SEEK_SET);
1163 
1164   /* Count MOs */
1165   fscanf(data->file, " Ene= %f\n", &orbenergy);
1166   fscanf(data->file, " Spin= %s\n", spin);
1167   fscanf(data->file, " Occup= %f\n", &occu);
1168 
1169   while (1) {
1170     int check_reads = 2;
1171     wave->num_orbitals++;
1172 
1173     /* skip over MO coeffs */
1174     while(check_reads == 2)
1175     {
1176       check_reads = fscanf(data->file, "%d %f", &dummy1, &dummy2);
1177     }
1178 
1179     nr  = fscanf(data->file, " Ene= %f\n", &orbenergy);
1180     nr += fscanf(data->file, " Spin= %s\n", spin);
1181     nr += fscanf(data->file, " Occup= %f\n", &occu);
1182 
1183     if (nr!=3 || toupper(spin[0])!='A')
1184       break;
1185   }
1186   //printf("found %d MOs!\n",wave->num_orbitals);
1187 
1188 
1189   /* Add wavefunction for spin beta */
1190   if (!strcmp(strtoupper(spin), "BETA")) {
1191     wave = add_wavefunction(data->qm_timestep);
1192     wave->spin = SPIN_BETA;
1193     wave->type = MOLFILE_WAVE_UNKNOWN;
1194     wave->exci = 0;
1195     wave->mult = 1;
1196     wave->num_coeffs = num_wave_coeff;
1197     wave->num_orbitals = 1;
1198 
1199     while (1) {
1200       int check_reads = 2;
1201       wave->num_orbitals++;
1202 
1203       /* skip over MO coeffs */
1204       while(check_reads == 2)
1205       {
1206         check_reads = fscanf(data->file, "%d %f", &dummy1, &dummy2);
1207       }
1208 
1209       nr  = fscanf(data->file, " Ene= %f\n", &orbenergy);
1210       nr += fscanf(data->file, " Spin= %s\n", spin);
1211       nr += fscanf(data->file, " Occup= %f\n", &occu);
1212 
1213       if (nr!=3 || toupper(spin[0])!='B' ||
1214           wave->num_orbitals>=num_wave_coeff) break;
1215     }
1216   }
1217 
1218   return TRUE;
1219 }
1220 
1221 
read_molecular_orbitals(qmdata_t * data)1222 static int read_molecular_orbitals(qmdata_t *data) {
1223   moldendata_t *moldendata = (moldendata_t *)data->format_specific_data;
1224   qm_wavefunction_t *wave;
1225 
1226   if (!data->qm_timestep || moldendata->coordsonly) return FALSE;
1227 
1228   /* Place file pointer on line after the [MO] keyword. */
1229   fseek(data->file, moldendata->filepos_mo, SEEK_SET);
1230 
1231   wave = &data->qm_timestep->wave[0];
1232   ALLOCATE(wave->wave_coeffs, float, wave->num_coeffs*wave->num_orbitals);
1233   // DEBUG
1234   /*
1235   printf("moldenplugin) num_coeffs   = %d\n", wave->num_coeffs);
1236   printf("moldenplugin) num_orbitals = %d\n", wave->num_orbitals);
1237   printf("moldenplugin) num_wave     = %d\n", data->qm_timestep->numwave);
1238   */
1239 
1240 
1241   /* Read wavefunction coefficients for spin alpha */
1242   if (!read_wave_coeffs(data->file, wave)) return FALSE;
1243 
1244   if (data->qm_timestep->numwave==1) return TRUE;
1245 
1246   /* Read wavefunction coefficients for spin beta */
1247   wave = &data->qm_timestep->wave[1];
1248   ALLOCATE(wave->wave_coeffs, float, wave->num_coeffs*wave->num_orbitals);
1249 
1250   if (!read_wave_coeffs(data->file, wave)) return FALSE;
1251 
1252   return TRUE;
1253 }
1254 
read_wave_coeffs(FILE * file,qm_wavefunction_t * wave)1255 static int read_wave_coeffs(FILE *file, qm_wavefunction_t *wave) {
1256   int i, j, nr;
1257   char buffer[1024];
1258   int AOid;
1259   float wf_coeff;
1260   char keystring[10];
1261   float *wave_coeffs = wave->wave_coeffs;
1262 
1263   /* This works for pruned and unpruned (all zero coeffs are preserved) MO representaitions.
1264    * Clearly this data redundancy should be avoided on systems with a few hundred atoms. */
1265 
1266   /* set all coeffs to zero */
1267   for (i=0; i<wave->num_orbitals; i++)
1268     for (j=0; j<wave->num_coeffs; j++)
1269         wave_coeffs[i*wave->num_coeffs + j] = 0.0;
1270 
1271   /* each molecular orbital must have at least 1 non-zero coeff in its represenation */
1272   /* eat Ene= Spin= Occup=  lines */
1273   eatline(file, 3);
1274   for (i=0; i<wave->num_orbitals; i++) {
1275     while(1) {
1276       int nr2;
1277       if (!fgets(buffer,1024,file)) return FALSE;
1278       nr = sscanf(buffer,"%d %f", &AOid, &wf_coeff);
1279       wave_coeffs[i*wave->num_coeffs+AOid-1] = wf_coeff;
1280 
1281       // DEBUG
1282       //printf("moldenplugin) %d,%d: %d %f\n", i, AOid-1, AOid, wave_coeffs[i*wave->num_coeffs+AOid-1]);
1283       nr2 = sscanf(buffer, "%s", keystring);
1284       if(!strcmp(keystring,"Ene=")||nr2==-1)
1285         break;
1286 
1287       if (nr==0) {
1288         printf("moldenplugin) Error reading wavefunction coefficients!\n");
1289         return FALSE;
1290       }
1291     }
1292     eatline(file, 2);
1293   }
1294   return TRUE;
1295 }
1296 
1297 /*************************************************************
1298  *
1299  * plugin registration
1300  *
1301  **************************************************************/
1302 static molfile_plugin_t plugin;
1303 
VMDPLUGIN_init()1304 VMDPLUGIN_API int VMDPLUGIN_init() {
1305   memset(&plugin, 0, sizeof(molfile_plugin_t));
1306   plugin.abiversion = vmdplugin_ABIVERSION;
1307   plugin.type = MOLFILE_PLUGIN_TYPE;
1308   plugin.name = "molden";
1309   plugin.prettyname = "Molden";
1310   plugin.author = "Markus Dittrich, Jan Saam, Alexey Titov";
1311   plugin.majorv = 0;
1312   plugin.minorv = 10;
1313   plugin.is_reentrant = VMDPLUGIN_THREADSAFE;
1314   plugin.filename_extension = "molden";
1315   plugin.open_file_read = open_molden_read;
1316   plugin.read_structure = read_molden_structure;
1317 
1318   plugin.read_timestep_metadata    = read_timestep_metadata;
1319   plugin.read_timestep             = read_timestep;
1320   plugin.read_qm_timestep_metadata = read_qm_timestep_metadata;
1321 
1322   plugin.read_qm_metadata = read_molden_metadata;
1323   plugin.read_qm_rundata  = read_molden_rundata;
1324 
1325   plugin.close_file_read = close_molden_read;
1326   return VMDPLUGIN_SUCCESS;
1327 }
1328 
VMDPLUGIN_register(void * v,vmdplugin_register_cb cb)1329 VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) {
1330   (*cb)(v, (vmdplugin_t *)&plugin);
1331   return VMDPLUGIN_SUCCESS;
1332 }
1333 
VMDPLUGIN_fini()1334 VMDPLUGIN_API int VMDPLUGIN_fini() {
1335   return VMDPLUGIN_SUCCESS;
1336 }
1337 
1338