1 /***************************************************************************
2 *cr
3 *cr (C) Copyright 1995-2016 The Board of Trustees of the
4 *cr University of Illinois
5 *cr All Rights Reserved
6 *cr
7 ***************************************************************************/
8
9 /***************************************************************************
10 * RCS INFORMATION:
11 *
12 * $RCSfile: qmplugin.h,v $
13 * $Author: johns $ $Locker: $ $State: Exp $
14 * $Revision: 1.23 $ $Date: 2016/11/28 05:01:54 $
15 *
16 ***************************************************************************/
17 /*******************************************************************
18 *
19 * Data structures and utility functions for plugins
20 * reading logfiles from QM packages.
21 *
22 ******************************************************************/
23
24 #ifndef QMPLUGIN_H
25 #define QMPLUGIN_H
26
27 #include <stdio.h>
28 #include <stdarg.h>
29 #include <math.h>
30 #include "molfile_plugin.h"
31
32
33 #define FALSE 0
34 #define TRUE 1
35
36 #define NONE 0
37
38 /* macros for shell types */
39 #define UNK_SHELL -666
40 #define SPD_SHELL -11
41 #define SP_SHELL -10
42 #define SPD_D_SHELL -5
43 #define SPD_P_SHELL -4
44 #define SPD_S_SHELL -3
45 #define SP_S_SHELL -2
46 #define SP_P_SHELL -1
47 #define S_SHELL 0
48 #define P_SHELL 1
49 #define D_SHELL 2
50 #define F_SHELL 3
51 #define G_SHELL 4
52 #define H_SHELL 5
53 #define I_SHELL 6
54
55 #define SPIN_ALPHA 0
56 #define SPIN_BETA 1
57
58 /* XXX the following macros should better be in molfileplugin.h
59 * since the same macros must be defined in VMD too. */
60
61
62 /* macros defining type of CI method (CITYP in GAMESS) */
63 #define CI_UNKNOWN -1
64 #define CI_NONE 0
65 #define CI_CIS 1
66 #define CI_ALDET 2
67 #define CI_ORMAS 3
68 #define CI_GUGA 4
69 #define CI_FSOCI 5
70 #define CI_GENCI 6
71
72 /* Basis set definition for a primitive */
73 typedef struct {
74 float exponent;
75 float contraction_coeff;
76 } prim_t;
77
78 /* Basis set definition for a shell */
79 typedef struct {
80 int numprims; /* number of primitives in this shell */
81 int type; /* S, P, D, F, ...
82 * just for convenience when retrieving info */
83 int wave_offset; /* index into wave_function array */
84 prim_t *prim; /* array of primitives */
85 } shell_t;
86
87 /* Basis set definition for an atom */
88 typedef struct {
89 char name[11]; /* atom name or type */
90 int atomicnum; /* atomic number (nuclear charge) */
91 int numshells; /* number of shells for this atom */
92 shell_t *shell; /* array of shells */
93 } basis_atom_t;
94
95
96 /* Atoms */
97 typedef struct
98 {
99 char type [11]; /* atom name or type */
100
101 int atomicnum; /* atomic number (nuclear charge) */
102
103 float x,y,z; /* coordinates of atom */
104 } qm_atom_t;
105
106
107 /* Wave function */
108 typedef struct {
109 int type; /**< CANONICAL, LOCALIZED, OTHER */
110 int spin; /**< 0 for alpha, 1 for beta */
111 int exci; /**< 0 for ground state, 1,2,3,... for excited states */
112 int mult; /**< spin multiplicity of the electronic state */
113 char info[MOLFILE_BUFSIZ]; /**< string for additional type info */
114
115 int num_orbitals; /**< number of orbitals that was really
116 * present in the output for this step */
117 int num_coeffs; /**< number of coefficients per orbital */
118 int has_orben; /**< flag for orbital energies */
119 int has_occup; /**< flag for orbital occupancies */
120 double energy; /**< total energy for this state */
121 float *wave_coeffs; /**< expansion coefficients for wavefunction in the
122 * form {orbital1(c1),orbital1(c2),.....,orbitalM(cN)} */
123 float *orb_energies; /**< array of orbital energies*/
124 float *orb_occupancies; /**< array of orbital occupancies */
125 } qm_wavefunction_t;
126
127
128 /* Timestep specific data.
129 * (Note that atoms are treated separately) */
130 typedef struct {
131 qm_wavefunction_t *wave;
132 int numwave; /* number of wavefunctions for this ts */
133 float *gradient; /* energy gradient for each atom */
134 int num_scfiter; /* number of SCF iterations */
135
136 double *scfenergies; /* scfenergies per trajectory point */
137
138 /* arrays with atom charges */
139 double *mulliken_charges; /* per-atom Mulliken charges */
140 double *lowdin_charges; /* per-atom Lowdin charges */
141 double *esp_charges; /* per-atom ESP charges */
142 int have_mulliken;
143 int have_lowdin;
144 int have_esp;
145 } qm_timestep_t;
146
147
148 /* Main QM plugin data structure */
149 typedef struct
150 {
151 /* File format specific data.
152 * This pointer must be cast to the according type by the plugin.
153 * Typically, this will be a struct containing various data that
154 * are needed by different functions during the file parsing process
155 * and cannot be sent through the molfile_plugin QM interface since
156 * the underlying data types are specific to the file format being
157 * read. */
158 void *format_specific_data;
159
160 FILE *file; /* the file we are reading */
161
162
163 /******************************************************
164 * calculation metadata (input data)
165 *****************************************************/
166
167 int numatoms; /* number of atoms in structure */
168 int runtype; /* type of calculation
169 * (ENERGY, OPTIMIZE, GRADIENT, ...) */
170 int scftype; /* UHF, RHF, ROHF, ... */
171 int dfttype; /* NONE, B3LYP, ..., */
172 int citype; /* NONE, GUGA, ... */
173
174 int mplevel; /* Moller-Plesset perturbation level */
175
176 char gbasis[10]; /* GBASIS of GAMESS run */
177
178 char basis_string[BUFSIZ]; /* basis name as "nice" string */
179
180 char runtitle[BUFSIZ]; /* title of gamess run */
181
182 char geometry[BUFSIZ]; /* either UNIQUE, CART or ZMP/ZMTMPC */
183 char guess[BUFSIZ]; /* type of guess method used */
184
185 char version_string[BUFSIZ]; /* GAMESS version used for run */
186
187
188 int nproc; /* Number processors used */
189 char memory[256]; /* Amount of memory used, e.g. 1Gb */
190
191 int totalcharge; /* Total charge of the system */
192 int multiplicity; /* Multiplicity of the system */
193 int num_electrons; /* Number of electrons */
194
195 char pointgroup[BUFSIZ]; /* Symmetry point group */
196 int naxis;
197 int order; /* Order of highest axis */
198
199 int mcscf_num_core; /* Number of MCSCF core orbitals
200 * (determines # valid orb energies
201 * for MCSCF natural orbitals) */
202
203 int max_opt_steps; /* Max. number of geom. opt. steps */
204 float opt_tol; /* gradient convergence tolerance,
205 * in Hartree/Bohr. */
206
207
208 /*********************************************************
209 * Basis set data
210 *********************************************************/
211
212 /* this array of floats stores the contraction coefficients
213 * and exponents for the basis functions:
214 * { exp(1), c-coeff(1), exp(2), c-coeff(2), .... }
215 * This holds also for double-zeta basis functions with
216 * exp(i) = exp(j) and c-coeff(i) != c-coeff(j). */
217 float *basis;
218
219 /* hierarchical basis set structures for each atom */
220 basis_atom_t *basis_set;
221
222 /* number of uncontracted basis functions in basis array */
223 int num_basis_funcs;
224
225 /* number of atoms listed in basis set */
226 int num_basis_atoms;
227
228 /* atomic number per atom in basis set */
229 int *atomicnum_per_basisatom;
230
231 /* number of shells per atom in basis set */
232 int *num_shells_per_atom;
233
234 /* the total number of atomic shells */
235 int num_shells;
236
237 /* number of primitives in shell i */
238 int *num_prim_per_shell;
239
240 /* type of each shell */
241 int *shell_types;
242
243 /* number of occupied spin A and B orbitals */
244 int num_occupied_A;
245 int num_occupied_B;
246
247 /* Max. size of the wave_function array per orbital.
248 * I.e. this is also the number of contracted
249 * cartesian gaussian basis functions or the size
250 * of the secular equation.
251 * While the actual # of MOs present can be different
252 * for each frame, this is the maximum number of
253 * possible occupied and virtual orbitals. */
254 int wavef_size;
255
256 /* Array of length 3*num_wave_f containing the exponents
257 * describing the cartesian components of the angular momentum.
258 * E.g. S={0 0 0}, Px={1 0 0}, Dxy={1 1 0}, or Fyyz={0 2 1}. */
259 int *angular_momentum;
260
261 /* Highest shell occuring in basis set */
262 int max_shell;
263
264
265 /******************************************************
266 * normal modes
267 *****************************************************/
268
269 int have_normal_modes; /* TRUE/FALSE flag indicating if we
270 * could properly read normal modes,
271 * wavenumbers and intensities. */
272
273 int nimag; /* Number of imaginary frequencies */
274 int *imag_modes; /* List of imaginary modes */
275
276 float *wavenumbers; /* rotational and translational DoF
277 * are included, but can be removed due
278 * to their zero frequencies */
279 float *intensities; /* Intensities of spectral lines */
280
281 float *normal_modes; /* the normal modes themselves */
282
283
284 /******************************************************
285 * internal coordinate stuff
286 *****************************************************/
287
288 int have_internals; /* TRUE/FALSE flag indicating if we
289 * have internal coordinates */
290
291 int have_cart_hessian; /* TRUE/FALSE flag indicating if we
292 * have a cartesian Hessian matrix */
293
294 int have_int_hessian; /* TRUE/FALSE flag indicating if we
295 * have a Hessian matrix in internal
296 * coordinates */
297
298 int nintcoords; /* Number of internal coordinates */
299 int nbonds; /* Number of bonds */
300 int nangles; /* Number of angles */
301 int ndiheds; /* Number of dihedrals */
302 int nimprops; /* Number of impropers */
303
304 int *bonds; /* bond list (atom tuples) */
305 int *angles; /* angle list (atom triples) */
306 int *dihedrals; /* dihedral list (atom quadrupels) */
307 int *impropers; /* improper list (atom quadrupels) */
308
309 double *internal_coordinates; /* value of internal coordinates */
310
311 /* XXX GAMESS */
312 /* the order of force constants has to match the internal
313 * coordinates in *bonds, *angles, *dihedrals */
314
315 double *bond_force_const; /* force constant for bonds */
316 double *angle_force_const; /* force constant for angles */
317 double *dihedral_force_const; /* force constant for dihedrals */
318 double *improper_force_const; /* force constant for impropers */
319
320 /*******************************************************
321 * Hessian matrices
322 *******************************************************/
323
324 double *carthessian; /* Hessian matrix in cartesian coordinates,
325 * dimension (3*numatoms)*(3*numatoms),
326 * single array of floats
327 * (row(1),row(2),...,row(numatoms))
328 */
329
330 double *inthessian; /* Hessian matrix in internal coordinates,
331 * dimension nintcoords*nintcoords,
332 * single array of floats
333 * (row(1),row(2),...,row(nintcoords))
334 */
335
336 /*******************************************************
337 * Trajectory related data
338 *******************************************************/
339
340 /* per timestep data like wavefunctions and scf iterations */
341 qm_timestep_t *qm_timestep;
342
343 /* array of atoms for the current timestep */
344 qm_atom_t *atoms;
345
346 /* flag to tell if SCF cycle and the geometry search converged.
347 * XXX should distinguish between SCF and geometry convergence
348 * in separate flags */
349 int status;
350
351 /* number of trajectory frames: */
352 int num_frames; /* total # frames */
353 int num_frames_read; /* # frames read in so far */
354 int num_frames_sent; /* # frames read sent to VMD so far */
355
356 /* flag to indicate wether we are done with reading frames */
357 int trajectory_done;
358
359 /* file positions of the beginning of each trajectory frame */
360 long *filepos_array;
361
362 /* file position indicator for the beginning of final section
363 * printed after the last trajectory frame */
364 long end_of_traj;
365
366 } qmdata_t;
367
368
369 /* #######################################################
370 *
371 * Function declarations
372 *
373 * ####################################################### */
374
375 /* Expand a set of symmetry unique atoms by creating images of
376 * the atoms so that we have the full coordinate set. */
377 static int symmetry_expand(qm_atom_t **atoms, int numunique, int natoms,
378 char *pg, int naxis);
379
380 /* Skip n lines at a time */
381 static void eatline(FILE * fd, int n);
382
383 /* Advance to the next non-white line */
384 static void eatwhitelines(FILE *fd);
385
386 /* Trim leading whitespaces from string */
387 static char* trimleft(char *);
388
389 /* Trim trailing whitespaces from string */
390 static char* trimright(char *);
391
392 /* Return 1 if the string consists of whitespace only */
393 static int iswhiteline(char *s);
394
395 /* Convert a string to upper case */
396 static char *strtoupper(char *s);
397
398
399 /* Place file pointer AFTER the line containing one of
400 * the keystrings. */
401 static int pass_keyline(FILE *file, const char *keystring,
402 const char *keystring2);
403
404 /* Place file pointer AT THE BEGINNING of the line containing
405 * a keystring. The keystrings are specified as a list of
406 * const char* function arguments. The last argument must be
407 * NULL in order to terminate the list. */
408 static int goto_keyline(FILE *file, ...);
409
410 /* Check wether keystring1 occurs before keystring2 and
411 * jumps back to beginning of search */
412 static int have_keyline(FILE *file, const char *keystring1,
413 const char *keystring2);
414
415
416 /* Print the current line but don't advance the file pointer. */
417 static void thisline(FILE *file);
418
419 /* Print next nonempty, nonwhite line but do not advance
420 * the file pointer. */
421 static void whereami(FILE *file);
422
423
424
425 /* #######################################################
426 *
427 * Function definitions
428 *
429 * ####################################################### */
430
431
432 /*********************************************************
433 *
434 * Allocates and initiates qmdata_t structure.
435 *
436 *********************************************************/
init_qmdata()437 static qmdata_t* init_qmdata() {
438 /* allocate memory for main data structure */
439 qmdata_t *data = (qmdata_t *) calloc(1, sizeof(qmdata_t));
440 if (data == NULL)
441 return NULL;
442
443 data->runtype = NONE;
444 data->scftype = NONE;
445 data->dfttype = NONE;
446 data->citype = NONE;
447 data->status = MOLFILE_QMSTATUS_UNKNOWN;
448 data->trajectory_done = FALSE;
449 data->have_internals = FALSE;
450 data->have_int_hessian = FALSE;
451 data->have_cart_hessian = FALSE;
452 data->have_normal_modes = FALSE;
453
454 /* initialize some of the character arrays */
455 memset(data->basis_string,0,sizeof(data->basis_string));
456 memset(data->version_string,0,sizeof(data->version_string));
457 memset(data->memory,0,sizeof(data->memory));
458
459 return data;
460 }
461
462
463 /*********************************************************
464 *
465 * functions to manipulate the wavefunction array
466 * in qm_timestep_t.
467 *
468 *********************************************************/
469
470 /* Increase wavefunction array in ts by one. */
add_wavefunction(qm_timestep_t * ts)471 static qm_wavefunction_t* add_wavefunction(qm_timestep_t *ts) {
472 if (ts->numwave) {
473 /* Add a new wavefunction */
474 ts->wave = (qm_wavefunction_t *)realloc(ts->wave,
475 (ts->numwave+1)*sizeof(qm_wavefunction_t));
476 memset(&ts->wave[ts->numwave], 0, sizeof(qm_wavefunction_t));
477 ts->numwave++;
478 } else {
479 /* We have no wavefunction for this timestep so create one */
480 ts->wave = (qm_wavefunction_t *)calloc(1, sizeof(qm_wavefunction_t));
481 ts->numwave = 1;
482 }
483
484 return &ts->wave[ts->numwave-1];
485 }
486
487
488 /* Replace the n-th wavefunction in ts with the last
489 * one and decrease the array length by one. */
replace_wavefunction(qm_timestep_t * ts,int n)490 static void replace_wavefunction(qm_timestep_t *ts, int n) {
491 if (ts->numwave>=2 && n>=0 && n<ts->numwave-1) {
492 qm_wavefunction_t *w1, *w2;
493 w2 = &ts->wave[n];
494 w1 = &ts->wave[ts->numwave-1];
495 free(w2->wave_coeffs);
496 free(w2->orb_energies);
497 free(w2->orb_occupancies);
498 memcpy(w2, w1, sizeof(qm_wavefunction_t));
499 ts->wave = (qm_wavefunction_t *) realloc(ts->wave,
500 (ts->numwave-1)*sizeof(qm_wavefunction_t));
501 ts->numwave--;
502 }
503 }
504
505
506 /* Delete the last wavefunction in ts */
del_wavefunction(qm_timestep_t * ts)507 static void del_wavefunction(qm_timestep_t *ts) {
508 if (ts->numwave) {
509 qm_wavefunction_t *w;
510 w = &ts->wave[ts->numwave-1];
511 free(w->wave_coeffs);
512 free(w->orb_energies);
513 free(w->orb_occupancies);
514 ts->numwave--;
515 ts->wave = (qm_wavefunction_t *)realloc(ts->wave,
516 ts->numwave*sizeof(qm_wavefunction_t));
517 }
518 }
519
520 /* Translate angular momentum string representation into
521 * a triplet angular momentum exponents for X, Y, Z.
522 * Angular momentum strings are a more human readable way
523 * to specify the order of coefficients for wavefunctions.
524 * So if the order of coefficients for D-shells implied
525 * in the QM file is xx, yy, zz, xy, xz, yz then you can
526 * use subsequent calls to this function in order to
527 * translate the strings into the more obscure but machine
528 * friendly angular momentum exponents required by the
529 * molfile_plugin interface.
530 *
531 * Example translations:
532 * S --> {0 0 0}
533 * X --> {1 0 0}
534 * Y --> {0 1 0}
535 * Z --> {0 0 1}
536 * XX --> {2 0 0}
537 * XY --> {1 1 0}
538 * YYZ --> {0 2 1}
539 */
angular_momentum_expon(int * ang_mom_expon,char * ang_mom_str)540 static void angular_momentum_expon(int *ang_mom_expon,
541 char *ang_mom_str) {
542 int i;
543 int xexp=0, yexp=0, zexp=0;
544
545 for (i=0; i<strlen(ang_mom_str); i++) {
546 switch (toupper(ang_mom_str[i])) {
547 case 'X':
548 xexp++;
549 break;
550 case 'Y':
551 yexp++;
552 break;
553 case 'Z':
554 zexp++;
555 break;
556 }
557 }
558 ang_mom_expon[0] = xexp;
559 ang_mom_expon[1] = yexp;
560 ang_mom_expon[2] = zexp;
561 }
562
563 /******************************************************
564 *
565 * matrix and vector functions
566 *
567 ******************************************************/
568
569 /* Degree-to-Radians and Radians-to-Degrees Conversion macros */
570 #define VMD_PI 3.14159265358979323846
571 #define DEGTORAD(a) (a*VMD_PI/180.0)
572 #define RADTODEG(a) (a*180.0/VMD_PI)
573
574 /* clears the matrix (resets it to identity) */
identity(float mat[16])575 static void identity(float mat[16]) {
576 memset(mat, 0, 16*sizeof(float));
577 mat[0]=1.0f;
578 mat[5]=1.0f;
579 mat[10]=1.0f;
580 mat[15]=1.0f;
581 }
582
583
584 /* Print a matrix for debugging purpose */
print_matrix4(const float mat[16])585 static void print_matrix4(const float mat[16]) {
586 int i, j;
587 for (i=0; i<4; i++) {
588 for (j=0; j<4; j++) {
589 printf("%f ", mat[4*j+i]);
590 }
591 printf("\n");
592 }
593 printf("\n");
594 }
595
596
597 /* premultiply the matrix by the given matrix */
multmatrix(const float * m1,float m2[16])598 static void multmatrix(const float *m1, float m2[16]) {
599 int i, j;
600 float tmp[4];
601
602 for (j=0; j<4; j++) {
603 tmp[0] = m2[j];
604 tmp[1] = m2[4+j];
605 tmp[2] = m2[8+j];
606 tmp[3] = m2[12+j];
607 for (i=0; i<4; i++) {
608 m2[4*i+j] = m1[4*i]*tmp[0] + m1[4*i+1]*tmp[1] +
609 m1[4*i+2]*tmp[2] + m1[4*i+3]*tmp[3];
610 }
611 }
612 }
613
614
615 /* performs a rotation around an axis (char == 'x', 'y', or 'z')
616 * angle is in degrees */
rot(float a,char axis,float mat[16])617 static void rot(float a, char axis, float mat[16]) {
618 double angle;
619 float m[16];
620 identity(m); /* create identity matrix */
621
622 angle = (double)DEGTORAD(a);
623
624 if (axis == 'x') {
625 m[0] = 1.f;
626 m[5] = (float)cos(angle);
627 m[10] = m[5];
628 m[6] = (float)sin(angle);
629 m[9] = -m[6];
630 } else if (axis == 'y') {
631 m[0] = (float)cos(angle);
632 m[5] = 1.f;
633 m[10] = m[0];
634 m[2] = (float) -sin(angle);
635 m[8] = -m[2];
636 } else if (axis == 'z') {
637 m[0] = (float)cos(angle);
638 m[5] = m[0];
639 m[10] = 1.f;
640 m[1] = (float)sin(angle);
641 m[4] = -m[1];
642 }
643
644 multmatrix(m, mat);
645 }
646
647
648 /* scale a matrix */
scale(float s,float m[16])649 static void scale(float s, float m[16]) {
650 float t[16];
651 identity(t); /* create identity matrix */
652 t[0] = s;
653 t[5] = s;
654 t[10] = s;
655 multmatrix(t, m);
656 }
657
658
659 /* reflect through mirror plane */
mirror(char axis,float m[16])660 static void mirror(char axis, float m[16]) {
661 scale(-1.f, m);
662 rot(180, axis, m);
663 }
664
665
666 /* multiplies a 3D point (first arg) by the Matrix, returns in second arg */
multpoint3d(const float * mat,const float opoint[3],float npoint[3])667 static void multpoint3d(const float *mat, const float opoint[3], float npoint[3]) {
668 float tmp[3];
669 float itmp3 = 1.0f / (opoint[0]*mat[3] + opoint[1]*mat[7] +
670 opoint[2]*mat[11] + mat[15]);
671 tmp[0] = itmp3*opoint[0];
672 tmp[1] = itmp3*opoint[1];
673 tmp[2] = itmp3*opoint[2];
674 npoint[0]=tmp[0]*mat[0] + tmp[1]*mat[4] + tmp[2]*mat[ 8] + itmp3*mat[12];
675 npoint[1]=tmp[0]*mat[1] + tmp[1]*mat[5] + tmp[2]*mat[ 9] + itmp3*mat[13];
676 npoint[2]=tmp[0]*mat[2] + tmp[1]*mat[6] + tmp[2]*mat[10] + itmp3*mat[14];
677 }
678
679
680 /* subtract 3rd vector from 2nd and put into 1st
681 * in other words, a = b - c */
qm_vec_sub(float * a,const float * b,const float * c)682 static void qm_vec_sub(float *a, const float *b, const float *c) {
683 a[0]=b[0]-c[0];
684 a[1]=b[1]-c[1];
685 a[2]=b[2]-c[2];
686 }
687
688
689 /* length of vector */
norm(const float * vect)690 static float norm(const float *vect) {
691 #if defined(_MSC_VER) || defined(__MINGW32__)
692 return sqrt(vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2]);
693 #else
694 return sqrtf(vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2]);
695 #endif
696 }
697
698
699 /******************************************************
700 *
701 * Symmetry
702 *
703 ******************************************************/
704
705 /* Expand a set of symmetry unique atoms by creating images of
706 * the atoms so that we have the full coordinate set.
707 * The atom array, the number of unique atoms, the total
708 * number of atoms, the pointgroup string and the order of
709 * the highest axis mut be provided. The atoms array will be resized
710 * and extended with the image atoms accordingly. */
symmetry_expand(qm_atom_t ** atoms,int numunique,int natoms,char * pg,int naxis)711 static int symmetry_expand(qm_atom_t **atoms, int numunique, int natoms,
712 char *pg, int naxis) {
713 int i, j;
714 float *unique, *image, *expanded;
715 int *indexmap;
716 int numexp;
717 float m[16];
718
719 printf("gamessplugin) Expanding %d coordinates for pointgroup %s, naxis=%d\n",
720 numunique, pg, naxis);
721
722 unique = calloc(3*numunique, sizeof(float));
723 image = calloc(3*numunique, sizeof(float));
724 indexmap = calloc(numunique, sizeof(int));
725 for (i=0; i<numunique; i++) {
726 unique[3*i ] = (*atoms)[i].x;
727 unique[3*i+1] = (*atoms)[i].y;
728 unique[3*i+2] = (*atoms)[i].z;
729 indexmap[i] = i;
730 /* printf("unique[%d]={%f %f %f}\n", i, unique[3*i], */
731 /* unique[3*i+1],unique[3*i+2]); */
732 }
733
734 /* Define the generating symmetry operations. */
735 /* XXX this is just a start,
736 * many molecules with higher symmetries will need
737 * two generating operations, e.g. a rotation and a
738 * reflection. */
739 identity(m);
740 if (!strcmp(pg, "CI")) {
741 scale(-1.f, m);
742 }
743 else if (!strcmp(pg, "CS")) {
744 mirror('y', m);
745 }
746 else if (!strcmp(pg, "CN")) {
747 rot(360.f/naxis, 'z', m);
748 }
749 else if (!strcmp(pg, "CNV")) {
750 rot(360.f/naxis, 'z', m);
751 }
752 else if (!strcmp(pg, "CNH")) {
753 rot(360.f/naxis, 'z', m);
754 }
755 else if (!strcmp(pg, "DNH")) {
756 rot(360.f/naxis, 'z', m);
757 }
758
759 for (i=0; i<numunique; i++) {
760 multpoint3d(m, &unique[3*i], &image[3*i]);
761 /* printf("image[%d]={%f %f %f}\n", i, image[3*i ], */
762 /* image[3*i+1], image[3*i+2]); */
763 }
764
765 expanded = calloc(3*numunique, sizeof(float));
766 memcpy(expanded, unique, 3*numunique* sizeof(float));
767 numexp=numunique;
768 for (i=0; i<numunique; i++) {
769 int found = 0;
770 for (j=0; j<numunique; j++) {
771 float d[3];
772 qm_vec_sub(d, &image[3*i], &unique[3*j]);
773 /* printf("%d,%d norm(d)=%f\n", i, j, norm(d)); */
774 if (norm(d)<0.001) {
775 found = 1;
776 break;
777 }
778 }
779
780 if (!found) {
781 expanded = realloc((float*)expanded, 3*(numexp+1)*sizeof(float));
782 indexmap = realloc((int*)indexmap, (numexp+1)*sizeof(int));
783 expanded[3*numexp ] = image[3*i ];
784 expanded[3*numexp+1] = image[3*i+1];
785 expanded[3*numexp+2] = image[3*i+2];
786 indexmap[numexp] = i;
787 numexp++;
788 }
789 }
790
791 /* for (i=0; i<numexp; i++) { */
792 /* printf("expanded[%d]={%f %f %f}\n", i, expanded[3*i], */
793 /* expanded[3*i+1], expanded[3*i+2]); */
794 /* } */
795
796 free(unique);
797 free(image);
798
799 if (natoms && numexp!=natoms) {
800 printf("gamessplugin) Couldn't expand symmetry unique atoms.\n");
801 free(expanded);
802 free(indexmap);
803 return FALSE;
804 }
805
806 /* XXX handling the natoms==0 case is futile unless we return
807 * the new number of atoms, so that the caller knows how
808 * many atoms there are. */
809 if (!natoms)
810 *atoms = calloc(numexp, sizeof(qm_atom_t));
811 else
812 *atoms = realloc((qm_atom_t*)(*atoms), numexp*sizeof(qm_atom_t));
813
814 for (i=numunique; i<numexp; i++) {
815 (*atoms)[i].x = expanded[3*i ];
816 (*atoms)[i].y = expanded[3*i+1];
817 (*atoms)[i].z = expanded[3*i+2];
818 (*atoms)[i].atomicnum = (*atoms)[indexmap[i]].atomicnum;
819 strncpy((*atoms)[i].type, (*atoms)[indexmap[i]].type, 10);
820 }
821
822 free(expanded);
823 free(indexmap);
824 return TRUE;
825 }
826
827
828 /******************************************************
829 *
830 * string functions and other parsing helpers
831 *
832 ******************************************************/
833
834 /* skip n lines at a time */
eatline(FILE * fd,int n)835 static void eatline(FILE * fd, int n)
836 {
837 int i;
838 for (i=0; i<n; i++) {
839 char readbuf[1025];
840 fgets(readbuf, 1024, fd);
841 }
842 }
843
844
845 /* Skip all following consecutive white lines. */
eatwhitelines(FILE * file)846 static void eatwhitelines(FILE *file) {
847 char buffer[BUFSIZ];
848 long filepos;
849 filepos = ftell(file);
850 while (fgets(buffer, sizeof(buffer), file)) {
851 if (strlen(trimright(buffer))) {
852 fseek(file, filepos, SEEK_SET);
853 break;
854 }
855 filepos = ftell(file);
856 }
857 }
858
859
860 /* Returns a pointer to the first non-whitespace
861 * character in a string.
862 * The input string s must be null-terminated.
863 */
trimleft(char * the_string)864 static char* trimleft(char* the_string)
865 {
866 char *new_string = the_string;
867 while ( (*new_string=='\n' || *new_string==' ' || *new_string=='\t') &&
868 (*new_string != '\0'))
869 {
870 new_string++;
871 }
872
873 return new_string;
874 }
875
876
877 /* Places a NULL character after the last non-whitespace
878 * character in a string and returns the pointer to the
879 * beginning of the string.
880 * The input string s must be null-terminated.
881 */
trimright(char * s)882 static char* trimright(char* s)
883 {
884 int i;
885 for (i=strlen(s)-1; i>=0; i--) {
886 if (!isspace(s[i])) break;
887 }
888 s[i+1] = '\0';
889
890 return s;
891 }
892
893
894 /* Return 1 if the string contains only whitespace. */
iswhiteline(char * s)895 static int iswhiteline(char *s) {
896 return (!strlen(trimleft(s)));
897 }
898
899
900 /* Convert a string to upper case */
strtoupper(char * s)901 static char *strtoupper(char *s) {
902 int i;
903 int sz = strlen(s);
904
905 if (s != NULL) {
906 for(i=0; i<sz; i++)
907 s[i] = toupper(s[i]);
908 }
909
910 return s;
911 }
912
913
914
915 /* Places file pointer AFTER the line containing one of the
916 * keystrings. If keystring2 is NULL then it will be ignored. */
pass_keyline(FILE * file,const char * keystring,const char * keystring2)917 static int pass_keyline(FILE *file, const char *keystring,
918 const char *keystring2) {
919 char buffer[BUFSIZ];
920 char *line;
921 int found = 0;
922 long filepos;
923 filepos = ftell(file);
924
925 do {
926 if (!fgets(buffer, sizeof(buffer), file)) {
927 fseek(file, filepos, SEEK_SET);
928 return 0;
929 }
930 line = trimleft(buffer);
931 if (strstr(line, keystring)) {
932 found = 1;
933 break;
934 }
935 else if (keystring2 && strstr(line, keystring2)) {
936 found = 2;
937 break;
938 }
939 } while (1);
940
941 if (!found) {
942 fseek(file, filepos, SEEK_SET);
943 }
944
945 return found;
946 }
947
948
949 /* Advances the file pointer until the first appearance
950 * of a line containing one of the given keystrings.
951 * The keystrings are specified as a list of const char*
952 * function arguments. The last argument must be NULL
953 * signifying the end of the keystring argument list.
954 * Returns the 1-based index of the keystring from the
955 * argument list for which a match was found.
956 * The file pointer will be placed AT THE BEGINNING of
957 * the line containing the matched keystring.
958 * If no keystring was found before EOF then the file
959 * is rewound to the position where the search started
960 * and 0 is returned.
961 */
goto_keyline(FILE * file,...)962 static int goto_keyline(FILE *file, ...) {
963 char buffer[BUFSIZ];
964 const char *keystring;
965 int found=0, loop=1;
966 long filepos, curline;
967 va_list argptr;
968 filepos = ftell(file);
969
970 /* loop over lines */
971 while (loop) {
972 int narg = 0;
973 curline = ftell(file);
974 if (!fgets(buffer, sizeof(buffer), file)) {
975 fseek(file, filepos, SEEK_SET);
976
977 return 0;
978 }
979
980 /* loop over the list of search strings */
981 va_start(argptr, file);
982 while ((keystring = va_arg(argptr, const char*)) != NULL) {
983 /* search for keystring in line buffer */
984 if (strstr(buffer, keystring)) {
985 found = narg+1;
986 /* rewind to beginning of current line */
987 fseek(file, curline, SEEK_SET);
988 loop = 0;
989 break;
990 }
991 narg++;
992 }
993 va_end (argptr);
994 };
995
996 if (!found) {
997 /* no match, rewind to beginning of search */
998 fseek(file, filepos, SEEK_SET);
999 }
1000
1001 return found;
1002 }
1003
1004 /* Check wether keystring1 occurs before keystring2 and
1005 * jumps back to beginning of search */
have_keyline(FILE * file,const char * keystring1,const char * keystring2)1006 static int have_keyline(FILE *file, const char *keystring1,
1007 const char *keystring2) {
1008 int found;
1009 long filepos;
1010 filepos = ftell(file);
1011
1012 found = pass_keyline(file, keystring1, keystring2);
1013
1014 fseek(file, filepos, SEEK_SET);
1015 return found;
1016 }
1017
1018
1019 /* Some MOLDEN files use Fortran-style notation where 'D' is used instead of
1020 * 'E' as exponential character.
1021 * If you pass this function a C string that contains floating point numbers
1022 * in Fortran-style scientific notation, it will find the "D" characters
1023 * that correspond to just the floating point numbers and fix them by
1024 * changing them to "E". Once converted, one should be able to
1025 * call scanf() on the string to parse it normally. The code modifies
1026 * the string in-place, which should work well for the molden plugin. */
1027
fpexpftoc(char * ftocstr)1028 static int fpexpftoc(char *ftocstr) {
1029 int convcnt = 0;
1030 int len = strlen(ftocstr);
1031
1032 // Compute string length, minus three chars, since any floating point
1033 // number in scientific notation is going to have at least a "-" or "+",
1034 // and one or more integer digits following the "E" or "D" character for
1035 // exponent notation.
1036 int lenm2 = len - 2;
1037
1038 // Replace "D" exponential characters with "E", so we can parse with
1039 // ANSI stdio routines.
1040 // Our loop starts at the second character in the string since the shortest
1041 // possible FP number in scientific notation would be "1D+01" or something
1042 // similar.
1043 int i;
1044 for (i=1; i<lenm2; i++) {
1045 // Check for a 'D' character. If we find one, then we look to see that the
1046 // preceding character is numeric, that the following character is +/-,
1047 // and that the character following +/- is also numeric. If all of the
1048 // necessary conditions are true then we replace the 'D' with an 'E'.
1049 // The strict checking should allow us to safely convert entire strings
1050 // of characters where there are also integers, and other kinds of strings
1051 // on the same line.
1052 if (ftocstr[i] == 'D' &&
1053 ((ftocstr[i+1] == '-') || (ftocstr[i+1] == '+')) &&
1054 ((ftocstr[i-1] >= '0') && (ftocstr[i-1] <= '9')) &&
1055 ((ftocstr[i+2] >= '0') && (ftocstr[i+2] <= '9'))) {
1056 ftocstr[i] = 'E';
1057 convcnt++;
1058 }
1059 }
1060
1061 return convcnt;
1062 }
1063
1064
1065
1066
1067 /*****************************************************
1068 *
1069 * Debugging functions
1070 *
1071 *****************************************************/
1072
1073 /* Print the current line in the format
1074 * "HERE) <contents of line>".
1075 * Doesn't advance the file pointer.
1076 */
thisline(FILE * file)1077 static void thisline(FILE *file) {
1078 char buffer[BUFSIZ];
1079 long filepos;
1080 filepos = ftell(file);
1081 if (!fgets(buffer, sizeof(buffer), file)) {
1082 if (feof(file)) printf("HERE) EOF\n");
1083 else printf("HERE) ????\n");
1084 return;
1085 }
1086 printf("HERE) %s\n", buffer);
1087 fseek(file, filepos, SEEK_SET);
1088 }
1089
1090
1091 /* Print all lines up to and including the next line that contains
1092 * non-whitespace characters. The file pointer will be restored to
1093 * its previous position.
1094 * Output format: "HERE) <contents of line>".
1095 * If the end of file has been reached print "HERE) EOF".
1096 */
whereami(FILE * file)1097 static void whereami(FILE *file) {
1098 char buffer[BUFSIZ];
1099 char *line;
1100 long filepos;
1101 filepos = ftell(file);
1102 do {
1103 if (!fgets(buffer, sizeof(buffer), file)) {
1104 if (feof(file)) printf("HERE) EOF\n");
1105 else printf("HERE) ????\n");
1106 return;
1107 }
1108 line = trimleft(buffer);
1109 printf("HERE) %s\n", buffer);
1110 } while (!strlen(line));
1111
1112 fseek(file, filepos, SEEK_SET);
1113 }
1114
1115
1116 #endif
1117