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

qmplugin.h

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: qmplugin.h,v $
00013  *      $Author: johns $       $Locker:  $             $State: Exp $
00014  *      $Revision: 1.25 $       $Date: 2018/12/12 03:41:21 $
00015  *
00016  ***************************************************************************/
00017 /*******************************************************************
00018  * 
00019  *  Data structures and utility functions for plugins
00020  *  reading logfiles from QM packages.
00021  *  
00022  ******************************************************************/
00023 
00024 #ifndef QMPLUGIN_H
00025 #define QMPLUGIN_H
00026 
00027 #include <stdio.h>
00028 #include <stdarg.h>
00029 #include <math.h>
00030 #include "molfile_plugin.h"
00031 
00032 
00033 #define FALSE 0
00034 #define TRUE  1
00035 
00036 #define NONE  0
00037 
00038 /* macros for shell types */
00039 #define UNK_SHELL -666
00040 #define SPD_SHELL   -11
00041 #define SP_SHELL    -10
00042 #define SPD_D_SHELL -5
00043 #define SPD_P_SHELL -4
00044 #define SPD_S_SHELL -3
00045 #define SP_S_SHELL -2
00046 #define SP_P_SHELL -1
00047 #define S_SHELL 0
00048 #define P_SHELL 1
00049 #define D_SHELL 2
00050 #define F_SHELL 3
00051 #define G_SHELL 4
00052 #define H_SHELL 5
00053 #define I_SHELL 6
00054 
00055 #define SPIN_ALPHA  0
00056 #define SPIN_BETA   1
00057 
00058 /* XXX the following macros should better be in molfileplugin.h
00059  * since the same macros must be defined in VMD too. */
00060 
00061 
00062 /* macros defining type of CI method (CITYP in GAMESS) */
00063 #define CI_UNKNOWN -1
00064 #define CI_NONE     0
00065 #define CI_CIS      1
00066 #define CI_ALDET    2
00067 #define CI_ORMAS    3
00068 #define CI_GUGA     4
00069 #define CI_FSOCI    5
00070 #define CI_GENCI    6
00071 
00072 /* Basis set definition for a primitive */
00073 typedef struct {
00074   float exponent;
00075   float contraction_coeff;
00076 } prim_t;
00077 
00078 /* Basis set definition for a shell */
00079 typedef struct {
00080   int numprims;      /* number of primitives in this shell */
00081   int type;          /* S, P, D, F, ...
00082                       * just for convenience when retrieving info */
00083   int wave_offset;   /* index into wave_function array */
00084   prim_t *prim;      /* array of primitives */
00085 } shell_t;
00086 
00087 /* Basis set definition for an atom */
00088 typedef struct {
00089   char name[11];  /* atom name or type */
00090   int atomicnum;  /* atomic number (nuclear charge) */
00091   int numshells;  /* number of shells for this atom */
00092   shell_t *shell; /* array of shells */
00093 } basis_atom_t;
00094 
00095 
00096 /* Atoms */
00097 typedef struct 
00098 {
00099   char type [11]; /* atom name or type */
00100 
00101   int atomicnum;  /* atomic number (nuclear charge) */
00102 
00103   float x,y,z;    /* coordinates of atom */
00104 } qm_atom_t;
00105 
00106 
00107 /* Wave function */
00108 typedef struct {
00109   int   type;           
00110   int   spin;           
00111   int   exci;           
00112   int   mult;           
00113   char info[MOLFILE_BUFSIZ]; 
00115   int   num_orbitals;   
00117   int   num_coeffs;     
00118   int   has_orben;      
00119   int   has_occup;      
00120   double energy;        
00121   float *wave_coeffs;   
00123   float *orb_energies;    
00124   float *orb_occupancies; 
00125 } qm_wavefunction_t;
00126 
00127 
00128 /* Timestep specific data.
00129  * (Note that atoms are treated separately) */
00130 typedef struct {
00131   qm_wavefunction_t *wave;
00132   int     numwave;          /* number of wavefunctions for this ts */
00133   float  *gradient;         /* energy gradient for each atom */
00134   int     num_scfiter;      /* number of SCF iterations */
00135 
00136   double *scfenergies;      /* scfenergies per trajectory point */
00137 
00138   /* arrays with atom charges */
00139   double *mulliken_charges; /* per-atom Mulliken charges */
00140   double *lowdin_charges;   /* per-atom Lowdin charges */
00141   double *esp_charges;      /* per-atom ESP charges */
00142   int   have_mulliken; 
00143   int   have_lowdin; 
00144   int   have_esp; 
00145 } qm_timestep_t;
00146 
00147 
00148 /* Main QM plugin data structure */
00149 typedef struct 
00150 {
00151   /* File format specific data.
00152    * This pointer must be cast to the according type by the plugin.
00153    * Typically, this will be a struct containing various data that
00154    * are needed by different functions during the file parsing process
00155    * and cannot be sent through the molfile_plugin QM interface since 
00156    * the underlying data types are specific to the file format being
00157    * read. */
00158   void *format_specific_data; 
00159 
00160   FILE *file;       /* the file we are reading */
00161 
00162 
00163   /******************************************************
00164    * calculation metadata (input data)
00165    *****************************************************/
00166 
00167   int numatoms;     /* number of atoms in structure */
00168   int runtype;      /* type of calculation 
00169                      * (ENERGY, OPTIMIZE, GRADIENT, ...) */
00170   int scftype;      /* UHF, RHF, ROHF, ... */
00171   int dfttype;      /* NONE, B3LYP, ...,   */
00172   int citype;       /* NONE, GUGA, ...     */
00173 
00174   int mplevel;      /* Moller-Plesset perturbation level */
00175 
00176   char gbasis[10];  /* GBASIS of GAMESS run */
00177 
00178   char basis_string[BUFSIZ]; /* basis name as "nice" string */
00179 
00180   char runtitle[BUFSIZ];  /* title of gamess run */
00181 
00182   char geometry[BUFSIZ];  /* either UNIQUE, CART or ZMP/ZMTMPC */
00183   char guess[BUFSIZ];     /* type of guess method used */
00184 
00185   char version_string[BUFSIZ]; /* GAMESS version used for run */
00186 
00187 
00188   int  nproc;          /* Number processors used */
00189   char memory[256];    /* Amount of memory used, e.g. 1Gb */
00190 
00191   int totalcharge;     /* Total charge of the system */
00192   int multiplicity;    /* Multiplicity of the system */
00193   int num_electrons;   /* Number of electrons */
00194 
00195   char pointgroup[BUFSIZ]; /* Symmetry point group */
00196   int naxis;
00197   int order;               /* Order of highest axis */
00198 
00199   int mcscf_num_core;  /* Number of MCSCF core orbitals 
00200                         * (determines # valid orb energies
00201                         * for MCSCF natural orbitals) */
00202 
00203   int max_opt_steps;   /* Max. number of geom. opt. steps */
00204   float opt_tol;       /* gradient convergence tolerance,
00205                         * in Hartree/Bohr. */
00206 
00207 
00208   /*********************************************************
00209    * Basis set data
00210    *********************************************************/
00211 
00212   /* this array of floats stores the contraction coefficients
00213    * and exponents for the basis functions:
00214    * { exp(1), c-coeff(1), exp(2), c-coeff(2), .... }
00215    * This holds also for double-zeta basis functions with
00216    * exp(i) = exp(j) and c-coeff(i) != c-coeff(j). */
00217   float *basis;
00218 
00219   /* hierarchical basis set structures for each atom */
00220   basis_atom_t *basis_set;
00221 
00222   /* number of uncontracted basis functions in basis array */
00223   int num_basis_funcs;
00224 
00225   /* number of atoms listed in basis set */
00226   int num_basis_atoms;
00227 
00228   /* atomic number per atom in basis set */
00229   int *atomicnum_per_basisatom;
00230 
00231   /* number of shells per atom in basis set */
00232   int *num_shells_per_atom;
00233 
00234   /* the total number of atomic shells */
00235   int num_shells;
00236 
00237   /* number of primitives in shell i */
00238   int *num_prim_per_shell;
00239 
00240   /* type of each shell */
00241   int *shell_types; 
00242 
00243   /* number of occupied spin A and B orbitals */
00244   int num_occupied_A;
00245   int num_occupied_B;
00246 
00247   /* Max. size of the wave_function array per orbital.
00248    * I.e. this is also the number of contracted
00249    * cartesian gaussian basis functions or the size
00250    * of the secular equation.
00251    * While the actual # of MOs present can be different
00252    * for each frame, this is the maximum number of 
00253    * possible occupied and virtual orbitals. */
00254   int wavef_size;
00255 
00256   /* Array of length 3*num_wave_f containing the exponents 
00257    * describing the cartesian components of the angular momentum. 
00258    * E.g. S={0 0 0}, Px={1 0 0}, Dxy={1 1 0}, or Fyyz={0 2 1}. */
00259   int *angular_momentum;
00260 
00261   /* Highest shell occuring in basis set */
00262   int max_shell;
00263 
00264 
00265   /******************************************************
00266    * normal modes
00267    *****************************************************/
00268 
00269   int have_normal_modes; /* TRUE/FALSE flag indicating if we
00270                           * could properly read normal modes,
00271                           * wavenumbers and intensities. */
00272 
00273   int  nimag;          /* Number of imaginary frequencies */
00274   int *imag_modes;     /* List of imaginary modes */
00275 
00276   float *wavenumbers;  /* rotational and translational DoF 
00277                         * are included, but can be removed due
00278                         * to their zero frequencies */
00279   float *intensities;  /* Intensities of spectral lines */
00280 
00281   float *normal_modes; /* the normal modes themselves */
00282 
00283 
00284   /******************************************************
00285    * internal coordinate stuff
00286    *****************************************************/
00287 
00288   int have_internals;  /* TRUE/FALSE flag indicating if we
00289                         * have internal coordinates */
00290 
00291   int have_cart_hessian; /* TRUE/FALSE flag indicating if we
00292                           * have a cartesian Hessian matrix */
00293 
00294   int have_int_hessian;  /* TRUE/FALSE flag indicating if we
00295                           * have a Hessian matrix in internal
00296                           * coordinates */
00297 
00298   int nintcoords;    /* Number of internal coordinates */
00299   int nbonds;        /* Number of bonds */
00300   int nangles;       /* Number of angles */
00301   int ndiheds;       /* Number of dihedrals */
00302   int nimprops;      /* Number of impropers */
00303 
00304   int *bonds;        /* bond list (atom tuples) */
00305   int *angles;       /* angle list (atom triples) */
00306   int *dihedrals;    /* dihedral list (atom quadrupels) */
00307   int *impropers;    /* improper list (atom quadrupels) */
00308 
00309   double *internal_coordinates; /* value of internal coordinates */ 
00310   
00311   /* XXX GAMESS */
00312   /* the order of force constants has to match the internal
00313    * coordinates in *bonds, *angles, *dihedrals */
00314 
00315   double *bond_force_const;     /* force constant for bonds */
00316   double *angle_force_const;    /* force constant for angles */
00317   double *dihedral_force_const; /* force constant for dihedrals */
00318   double *improper_force_const; /* force constant for impropers */
00319 
00320   /*******************************************************
00321    * Hessian matrices
00322    *******************************************************/
00323 
00324   double *carthessian;  /* Hessian matrix in cartesian coordinates,
00325                          * dimension (3*numatoms)*(3*numatoms),
00326                          * single array of floats 
00327                          * (row(1),row(2),...,row(numatoms))
00328                          */
00329 
00330   double *inthessian;  /* Hessian matrix in internal coordinates,
00331                         * dimension nintcoords*nintcoords,
00332                         * single array of floats 
00333                         * (row(1),row(2),...,row(nintcoords))
00334                         */
00335 
00336   /*******************************************************
00337    * Trajectory related data
00338    *******************************************************/
00339 
00340   /* per timestep data like wavefunctions and scf iterations */
00341   qm_timestep_t *qm_timestep;
00342 
00343   /* array of atoms for the current timestep */
00344   qm_atom_t *atoms;
00345 
00346   /* flag to tell if SCF cycle and the geometry search converged.
00347    * XXX should distinguish between SCF and geometry convergence
00348    * in separate flags */
00349   int status;
00350 
00351   /* number of trajectory frames: */
00352   int num_frames;       /* total # frames */
00353   int num_frames_read;  /* # frames read in so far */
00354   int num_frames_sent;  /* # frames read sent to VMD so far */
00355 
00356   /* flag to indicate wether we are done with reading frames */
00357   int trajectory_done;
00358 
00359   /* file positions of the beginning of each trajectory frame */
00360   long *filepos_array;
00361 
00362   /* file position indicator for the beginning of final section 
00363    * printed after the last trajectory frame */
00364   long end_of_traj;
00365 
00366 } qmdata_t;
00367 
00368 
00369 /* #######################################################
00370  *
00371  * Function declarations
00372  *
00373  * ####################################################### */
00374 
00375 /* Expand a set of symmetry unique atoms by creating images of
00376  * the atoms so that we have the full coordinate set. */
00377 static int symmetry_expand(qm_atom_t **atoms, int numunique, int natoms, 
00378                             char *pg, int naxis);
00379 
00380 /* Skip n lines at a time */
00381 static void eatline(FILE * fd, int n);
00382 
00383 /* Advance to the next non-white line */
00384 static void eatwhitelines(FILE *fd);
00385 
00386 /* Trim leading whitespaces from string */
00387 static char* trimleft(char *);
00388 
00389 /* Trim trailing whitespaces from string */
00390 static char* trimright(char *);
00391 
00392 /* Return 1 if the string consists of whitespace only */
00393 static int iswhiteline(char *s);
00394 
00395 /* Convert a string to upper case */
00396 static char *strtoupper(char *s);
00397 
00398 
00399 /* Place file pointer AFTER the line containing one of
00400  * the keystrings. */
00401 static int pass_keyline(FILE *file, const char *keystring,
00402                         const char *keystring2);
00403 
00404 /* Place file pointer AT THE BEGINNING of the line containing 
00405  * a keystring. The keystrings are specified as a list of
00406  * const char* function arguments. The last argument must be 
00407  * NULL in order to terminate the list. */
00408 static int goto_keyline(FILE *file, ...);
00409 
00410 /* Check wether keystring1 occurs before keystring2 and
00411  * jumps back to beginning of search */
00412 static int have_keyline(FILE *file, const char *keystring1,
00413                         const char *keystring2);
00414 
00415 
00416 /* Print the current line but don't advance the file pointer. */
00417 static void thisline(FILE *file);
00418 
00419 /* Print next nonempty, nonwhite line but do not advance
00420  * the file pointer. */
00421 static void whereami(FILE *file);
00422 
00423 
00424 
00425 /* #######################################################
00426  *
00427  * Function definitions
00428  *
00429  * ####################################################### */
00430 
00431 
00432 /*********************************************************
00433  *
00434  * Allocates and initiates qmdata_t structure.
00435  *
00436  *********************************************************/
00437 static qmdata_t* init_qmdata() {
00438   /* allocate memory for main data structure */
00439   qmdata_t *data = (qmdata_t *) calloc(1, sizeof(qmdata_t));
00440   if (data == NULL)
00441     return NULL;
00442 
00443   data->runtype = NONE;
00444   data->scftype = NONE;
00445   data->dfttype = NONE;
00446   data->citype  = NONE;
00447   data->status = MOLFILE_QMSTATUS_UNKNOWN;
00448   data->trajectory_done   = FALSE;
00449   data->have_internals    = FALSE;
00450   data->have_int_hessian  = FALSE;
00451   data->have_cart_hessian = FALSE;
00452   data->have_normal_modes = FALSE;
00453 
00454   /* initialize some of the character arrays */
00455   memset(data->basis_string,0,sizeof(data->basis_string));
00456   memset(data->version_string,0,sizeof(data->version_string));
00457   memset(data->memory,0,sizeof(data->memory));
00458 
00459   return data;
00460 }
00461 
00462 
00463 /*********************************************************
00464  *
00465  * functions to manipulate the wavefunction array 
00466  * in qm_timestep_t.
00467  *
00468  *********************************************************/
00469 
00470 /* Increase wavefunction array in ts by one. */
00471 static qm_wavefunction_t* add_wavefunction(qm_timestep_t *ts) {
00472   if (ts->numwave) {
00473     /* Add a new wavefunction */
00474     ts->wave = (qm_wavefunction_t *) realloc(ts->wave,
00475                         (ts->numwave+1)*sizeof(qm_wavefunction_t));
00476     memset(&ts->wave[ts->numwave], 0, sizeof(qm_wavefunction_t));
00477     ts->numwave++;
00478   } else {
00479     /* We have no wavefunction for this timestep so create one */
00480     ts->wave = (qm_wavefunction_t *) calloc(1, sizeof(qm_wavefunction_t));
00481     ts->numwave = 1;
00482   }
00483 
00484   return &ts->wave[ts->numwave-1];
00485 }
00486 
00487 
00488 /* Replace the n-th wavefunction in ts with the last
00489  * one and decrease the array length by one. */
00490 static void replace_wavefunction(qm_timestep_t *ts, int n) {
00491   if (ts->numwave>=2 && n>=0 && n<ts->numwave-1) {
00492     qm_wavefunction_t *w1, *w2;
00493     w2 = &ts->wave[n];
00494     w1 = &ts->wave[ts->numwave-1];
00495     free(w2->wave_coeffs);
00496     free(w2->orb_energies);
00497     free(w2->orb_occupancies);
00498     memcpy(w2, w1, sizeof(qm_wavefunction_t));
00499     ts->wave = (qm_wavefunction_t *) realloc(ts->wave,
00500                    (ts->numwave-1)*sizeof(qm_wavefunction_t));
00501     ts->numwave--;
00502   }
00503 }
00504 
00505 
00506 /* Delete the last wavefunction in ts */
00507 static void del_wavefunction(qm_timestep_t *ts) {
00508   if (ts->numwave) {
00509     qm_wavefunction_t *w;
00510     w = &ts->wave[ts->numwave-1];
00511     free(w->wave_coeffs);
00512     free(w->orb_energies);
00513     free(w->orb_occupancies);
00514     ts->numwave--;
00515     ts->wave = (qm_wavefunction_t *) realloc(ts->wave,
00516                         ts->numwave*sizeof(qm_wavefunction_t));    
00517   }
00518 }
00519 
00520 /* Translate angular momentum string representation into
00521  * a triplet angular momentum exponents for X, Y, Z.
00522  * Angular momentum strings are a more human readable way
00523  * to specify the order of coefficients for wavefunctions.
00524  * So if the order of coefficients for D-shells implied
00525  * in the QM file is xx, yy, zz, xy, xz, yz then you can
00526  * use subsequent calls to this function in order to 
00527  * translate the strings into the more obscure but machine
00528  * friendly angular momentum exponents required by the 
00529  * molfile_plugin interface.
00530  *
00531  * Example translations:
00532  * S   --> {0 0 0}
00533  * X   --> {1 0 0}
00534  * Y   --> {0 1 0}
00535  * Z   --> {0 0 1}
00536  * XX  --> {2 0 0}
00537  * XY  --> {1 1 0}
00538  * YYZ --> {0 2 1}
00539  */
00540 static void angular_momentum_expon(int  *ang_mom_expon,
00541                                    char *ang_mom_str) {
00542   int i;
00543   int xexp=0, yexp=0, zexp=0;
00544 
00545   for (i=0; i<strlen(ang_mom_str); i++) {
00546     switch (toupper(ang_mom_str[i])) {
00547     case 'X':
00548       xexp++;
00549       break;
00550     case 'Y':
00551       yexp++;
00552       break;
00553     case 'Z':
00554       zexp++;
00555       break;
00556     }
00557   }
00558   ang_mom_expon[0] = xexp;
00559   ang_mom_expon[1] = yexp;
00560   ang_mom_expon[2] = zexp;
00561 }
00562 
00563 /******************************************************
00564  *
00565  * matrix and vector functions
00566  *
00567  ******************************************************/
00568 
00569 /* Degree-to-Radians and Radians-to-Degrees Conversion macros */
00570 #define VMD_PI      3.14159265358979323846
00571 #define DEGTORAD(a)     (a*VMD_PI/180.0)
00572 #define RADTODEG(a)     (a*180.0/VMD_PI)
00573 
00574 /* clears the matrix (resets it to identity) */
00575 static void identity(float mat[16]) {
00576   memset(mat, 0, 16*sizeof(float));
00577   mat[0]=1.0f;
00578   mat[5]=1.0f;
00579   mat[10]=1.0f;
00580   mat[15]=1.0f;
00581 }
00582 
00583 
00584 /* Print a matrix for debugging purpose */
00585 static void print_matrix4(const float mat[16]) {
00586   int i, j;
00587   for (i=0; i<4; i++) {
00588     for (j=0; j<4; j++) {
00589       printf("%f ", mat[4*j+i]);
00590     }
00591     printf("\n");
00592   }
00593   printf("\n");
00594 }
00595 
00596 
00597 /* premultiply the matrix by the given matrix */
00598 static void multmatrix(const float *m1, float m2[16]) {
00599   int i, j;
00600   float tmp[4];
00601 
00602   for (j=0; j<4; j++) {
00603     tmp[0] = m2[j];
00604     tmp[1] = m2[4+j];
00605     tmp[2] = m2[8+j]; 
00606     tmp[3] = m2[12+j];
00607     for (i=0; i<4; i++) {
00608       m2[4*i+j] = m1[4*i]*tmp[0] + m1[4*i+1]*tmp[1] +
00609         m1[4*i+2]*tmp[2] + m1[4*i+3]*tmp[3]; 
00610     }
00611   } 
00612 }
00613 
00614 
00615 /* performs a rotation around an axis (char == 'x', 'y', or 'z')
00616  * angle is in degrees */
00617 static void rot(float a, char axis, float mat[16]) {
00618   double angle;
00619   float m[16];
00620   identity(m);                  /* create identity matrix */
00621 
00622   angle = (double)DEGTORAD(a);
00623 
00624   if (axis == 'x') {
00625     m[0]  = 1.f;
00626     m[5]  = (float)cos(angle);
00627     m[10] = m[5];
00628     m[6]  = (float)sin(angle);
00629     m[9]  = -m[6];
00630   } else if (axis == 'y') {
00631     m[0]  = (float)cos(angle);
00632     m[5]  = 1.f;
00633     m[10] = m[0];
00634     m[2]  = (float) -sin(angle);
00635     m[8]  = -m[2];
00636   } else if (axis == 'z') {
00637     m[0]  = (float)cos(angle);
00638     m[5]  = m[0];
00639     m[10] = 1.f;
00640     m[1]  = (float)sin(angle);
00641     m[4]  = -m[1];
00642   }
00643 
00644   multmatrix(m, mat);
00645 }
00646 
00647 
00648 /* scale a matrix */
00649 static void scale(float s, float m[16]) {
00650   float t[16];
00651   identity(t);          /* create identity matrix */
00652   t[0]  = s;
00653   t[5]  = s;
00654   t[10] = s;
00655   multmatrix(t, m);
00656 }
00657 
00658 
00659 /* reflect through mirror plane */
00660 static void mirror(char axis, float m[16]) {
00661   scale(-1.f, m);
00662   rot(180, axis, m);
00663 }
00664 
00665 
00666 /* multiplies a 3D point (first arg) by the Matrix, returns in second arg */
00667 static void multpoint3d(const float *mat, const float opoint[3], float npoint[3]) {
00668   float tmp[3];
00669   float itmp3 = 1.0f / (opoint[0]*mat[3] + opoint[1]*mat[7] +
00670                         opoint[2]*mat[11] + mat[15]);
00671   tmp[0] = itmp3*opoint[0];
00672   tmp[1] = itmp3*opoint[1];
00673   tmp[2] = itmp3*opoint[2];
00674   npoint[0]=tmp[0]*mat[0] + tmp[1]*mat[4] + tmp[2]*mat[ 8] + itmp3*mat[12];
00675   npoint[1]=tmp[0]*mat[1] + tmp[1]*mat[5] + tmp[2]*mat[ 9] + itmp3*mat[13];
00676   npoint[2]=tmp[0]*mat[2] + tmp[1]*mat[6] + tmp[2]*mat[10] + itmp3*mat[14];
00677 }
00678 
00679 
00680 /* subtract 3rd vector from 2nd and put into 1st
00681  * in other words, a = b - c  */
00682 static void qm_vec_sub(float *a, const float *b, const float *c) {
00683   a[0]=b[0]-c[0];
00684   a[1]=b[1]-c[1];
00685   a[2]=b[2]-c[2];
00686 }
00687 
00688 
00689 /* length of vector */
00690 static float norm(const float *vect) {
00691 #if defined(_MSC_VER) || defined(__MINGW32__)
00692   return sqrt(vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2]);
00693 #else
00694   return sqrtf(vect[0]*vect[0] + vect[1]*vect[1] + vect[2]*vect[2]);
00695 #endif
00696 }
00697 
00698 
00699 /******************************************************
00700  *
00701  * Symmetry
00702  *
00703  ******************************************************/
00704 
00705 /* Expand a set of symmetry unique atoms by creating images of
00706  * the atoms so that we have the full coordinate set.
00707  * The atom array, the number of unique atoms, the total 
00708  * number of atoms, the pointgroup string and the order of
00709  * the highest axis mut be provided. The atoms array will be resized
00710  * and extended with the image atoms accordingly. */
00711 static int symmetry_expand(qm_atom_t **atoms, int numunique, int natoms, 
00712                             char *pg, int naxis) {
00713   int i, j;
00714   float *unique, *image, *expanded;
00715   int *indexmap;
00716   int numexp;
00717   float m[16];
00718 
00719   printf("gamessplugin) Expanding %d coordinates for pointgroup %s, naxis=%d\n",
00720          numunique, pg, naxis);
00721 
00722   unique = (float *) calloc(3*numunique, sizeof(float));
00723   image  = (float *) calloc(3*numunique, sizeof(float));
00724   indexmap = (int *) calloc(numunique, sizeof(int));
00725   for (i=0; i<numunique; i++) {
00726     unique[3*i  ] = (*atoms)[i].x;
00727     unique[3*i+1] = (*atoms)[i].y;
00728     unique[3*i+2] = (*atoms)[i].z;
00729     indexmap[i] = i;
00730 /*     printf("unique[%d]={%f %f %f}\n", i, unique[3*i], */
00731 /*            unique[3*i+1],unique[3*i+2]); */
00732   }
00733 
00734   /* Define the generating symmetry operations. */
00735   /* XXX this is just a start, 
00736    *     many molecules with higher symmetries will need 
00737    *     two generating operations, e.g. a rotation and a
00738    *     reflection. */
00739   identity(m);
00740   if (!strcmp(pg, "CI")) {
00741     scale(-1.f, m);
00742   }
00743   else if (!strcmp(pg, "CS")) {
00744     mirror('y', m);
00745   }
00746   else if (!strcmp(pg, "CN")) {
00747     rot(360.f/naxis, 'z', m);
00748   }
00749   else if (!strcmp(pg, "CNV")) {
00750     rot(360.f/naxis, 'z', m);
00751   }
00752   else if (!strcmp(pg, "CNH")) {
00753     rot(360.f/naxis, 'z', m);
00754   }
00755   else if (!strcmp(pg, "DNH")) {
00756     rot(360.f/naxis, 'z', m);
00757   }
00758 
00759   for (i=0; i<numunique; i++) {
00760     multpoint3d(m, &unique[3*i], &image[3*i]);
00761 /*     printf("image[%d]={%f %f %f}\n", i, image[3*i  ], */
00762 /*            image[3*i+1], image[3*i+2]); */
00763   }
00764 
00765   expanded = (float *) calloc(3*numunique, sizeof(float));
00766   memcpy(expanded, unique, 3*numunique* sizeof(float));
00767   numexp=numunique;
00768   for (i=0; i<numunique; i++) {
00769     int found = 0;
00770     for (j=0; j<numunique; j++) {
00771       float d[3];
00772       qm_vec_sub(d, &image[3*i], &unique[3*j]);
00773       /* printf("%d,%d norm(d)=%f\n", i, j, norm(d)); */
00774       if (norm(d)<0.001) {
00775         found = 1;
00776         break;
00777       }
00778     }
00779 
00780     if (!found) {
00781       expanded = (float*) realloc((float*)expanded, 3*(numexp+1)*sizeof(float));
00782       indexmap = (int*) realloc((int*)indexmap, (numexp+1)*sizeof(int));
00783       expanded[3*numexp  ] = image[3*i  ];
00784       expanded[3*numexp+1] = image[3*i+1];
00785       expanded[3*numexp+2] = image[3*i+2];
00786       indexmap[numexp] = i;
00787       numexp++;
00788     }
00789   }
00790 
00791 /*   for (i=0; i<numexp; i++) { */
00792 /*     printf("expanded[%d]={%f %f %f}\n", i, expanded[3*i], */
00793 /*            expanded[3*i+1], expanded[3*i+2]); */
00794 /*   } */
00795 
00796   free(unique);
00797   free(image);
00798 
00799   if (natoms && numexp!=natoms) {
00800     printf("gamessplugin) Couldn't expand symmetry unique atoms.\n");
00801     free(expanded);
00802     free(indexmap);
00803     return FALSE;
00804   }
00805 
00806   /* XXX handling the natoms==0 case is futile unless we return 
00807    *     the new number of atoms, so that the caller knows how
00808    *     many atoms there are. */
00809   if (!natoms)
00810     *atoms = (qm_atom_t *) calloc(numexp, sizeof(qm_atom_t));
00811   else 
00812     *atoms = (qm_atom_t *) realloc((qm_atom_t*)(*atoms), numexp*sizeof(qm_atom_t));
00813 
00814   for (i=numunique; i<numexp; i++) {
00815     (*atoms)[i].x = expanded[3*i  ];
00816     (*atoms)[i].y = expanded[3*i+1];
00817     (*atoms)[i].z = expanded[3*i+2];
00818     (*atoms)[i].atomicnum = (*atoms)[indexmap[i]].atomicnum;
00819     strncpy((*atoms)[i].type, (*atoms)[indexmap[i]].type, 10);
00820   }
00821 
00822   free(expanded);
00823   free(indexmap);
00824   return TRUE;
00825 }
00826 
00827 
00828 /******************************************************
00829  *
00830  * string functions and other parsing helpers
00831  *
00832  ******************************************************/
00833 
00834 /* skip n lines at a time */
00835 static void eatline(FILE * fd, int n)
00836 {
00837   int i;
00838   for (i=0; i<n; i++) {
00839     char readbuf[1025];
00840     fgets(readbuf, 1024, fd);
00841   }
00842 }
00843 
00844 
00845 /* Skip all following consecutive white lines. */
00846 static void eatwhitelines(FILE *file) {
00847   char buffer[BUFSIZ];
00848   long filepos;
00849   filepos = ftell(file);
00850   while (fgets(buffer, sizeof(buffer), file)) {
00851     if (strlen(trimright(buffer))) {
00852       fseek(file, filepos, SEEK_SET);
00853       break;
00854     }
00855     filepos = ftell(file);
00856   }
00857 }
00858 
00859 
00860 /* Returns a pointer to the first non-whitespace
00861  * character in a string.
00862  * The input string s must be null-terminated. 
00863  */
00864 static char* trimleft(char* the_string)
00865 {
00866   char *new_string = the_string;
00867   while ( (*new_string=='\n' || *new_string==' ' || *new_string=='\t') && 
00868           (*new_string != '\0'))
00869   {
00870     new_string++;
00871   }
00872 
00873   return new_string;
00874 }
00875 
00876 
00877 /* Places a NULL character after the last non-whitespace
00878  * character in a string and returns the pointer to the
00879  * beginning of the string. 
00880  * The input string s must be null-terminated.
00881  */
00882 static char* trimright(char* s)
00883 {
00884   int i;
00885   for (i=strlen(s)-1; i>=0; i--) {
00886     if (!isspace(s[i])) break;
00887   }
00888   s[i+1] = '\0';
00889 
00890   return s;
00891 }
00892 
00893 
00894 /* Return 1 if the string contains only whitespace. */
00895 static int iswhiteline(char *s) {
00896   return (!strlen(trimleft(s)));
00897 }
00898 
00899 
00900 /* Convert a string to upper case */
00901 static char *strtoupper(char *s) {
00902   int i;
00903   int sz = strlen(s);
00904 
00905   if (s != NULL) {
00906     for(i=0; i<sz; i++)
00907       s[i] = toupper(s[i]);
00908   }
00909 
00910   return s;
00911 }
00912 
00913 
00914 
00915 /* Places file pointer AFTER the line containing one of the
00916  * keystrings. If keystring2 is NULL then it will be ignored. */
00917 static int pass_keyline(FILE *file, const char *keystring,
00918         const char *keystring2) {
00919   char buffer[BUFSIZ];
00920   char *line;
00921   int found = 0;
00922   long filepos;
00923   filepos = ftell(file);
00924 
00925   do {
00926     if (!fgets(buffer, sizeof(buffer), file)) {
00927       fseek(file, filepos, SEEK_SET);
00928       return 0;
00929     }
00930     line = trimleft(buffer);
00931     if (strstr(line, keystring)) {
00932       found = 1;
00933       break;
00934     }
00935     else if (keystring2 && strstr(line, keystring2)) {
00936       found = 2;
00937       break;
00938     }
00939   } while (1);
00940     
00941   if (!found) {
00942     fseek(file, filepos, SEEK_SET);
00943   }
00944 
00945   return found;
00946 }
00947 
00948 
00949 /* Advances the file pointer until the first appearance
00950  * of a line containing one of the given keystrings.
00951  * The keystrings are specified as a list of const char*
00952  * function arguments. The last argument must be NULL
00953  * signifying the end of the keystring argument list.
00954  * Returns the 1-based index of the keystring from the
00955  * argument list for which a match was found.
00956  * The file pointer will be placed AT THE BEGINNING of
00957  * the line containing the matched keystring.
00958  * If no keystring was found before EOF then the file
00959  * is rewound to the position where the search started
00960  * and 0 is returned.
00961  */
00962 static int goto_keyline(FILE *file, ...) {
00963   char buffer[BUFSIZ];
00964   const char *keystring;
00965   int found=0, loop=1;
00966   long filepos, curline;
00967   va_list argptr;
00968   filepos = ftell(file);
00969 
00970   /* loop over lines */
00971   while (loop) {
00972     int narg = 0;
00973     curline = ftell(file);
00974     if (!fgets(buffer, sizeof(buffer), file)) {
00975       fseek(file, filepos, SEEK_SET);
00976 
00977       return 0;
00978     }
00979 
00980     /* loop over the list of search strings */
00981     va_start(argptr, file);
00982     while ((keystring = va_arg(argptr, const char*)) != NULL) {
00983       /* search for keystring in line buffer */
00984       if (strstr(buffer, keystring)) {
00985         found = narg+1;
00986         /* rewind to beginning of current line */
00987         fseek(file, curline, SEEK_SET);
00988         loop = 0;
00989         break;
00990       }
00991       narg++;
00992     }
00993     va_end (argptr);
00994   };
00995     
00996   if (!found) {
00997     /* no match, rewind to beginning of search */
00998     fseek(file, filepos, SEEK_SET);
00999   }
01000 
01001   return found;
01002 }
01003 
01004 /* Check wether keystring1 occurs before keystring2 and
01005  * jumps back to beginning of search */
01006 static int have_keyline(FILE *file, const char *keystring1,
01007         const char *keystring2) {
01008   int found;
01009   long filepos;
01010   filepos = ftell(file);
01011 
01012   found = pass_keyline(file, keystring1, keystring2);
01013 
01014   fseek(file, filepos, SEEK_SET);
01015   return found;
01016 }
01017 
01018 
01019 /* Some MOLDEN files use Fortran-style notation where 'D' is used instead of 
01020  * 'E' as exponential character. 
01021  * If you pass this function a C string that contains floating point numbers
01022  * in Fortran-style scientific notation, it will find the "D" characters
01023  * that correspond to just the floating point numbers and fix them by
01024  * changing them to "E".  Once converted, one should be able to
01025  * call scanf() on the string to parse it normally.  The code modifies
01026  * the string in-place, which should work well for the molden plugin. */
01027 
01028 static int fpexpftoc(char *ftocstr) {
01029   int convcnt = 0;
01030   int len = strlen(ftocstr);
01031 
01032   // Compute string length, minus three chars, since any floating point 
01033   // number in scientific notation is going to have at least a "-" or "+",
01034   // and one or more integer digits following the "E" or "D" character for
01035   // exponent notation.
01036   int lenm2 = len - 2;
01037 
01038   // Replace "D" exponential characters with "E", so we can parse with 
01039   // ANSI stdio routines.
01040   // Our loop starts at the second character in the string since the shortest
01041   // possible FP number in scientific notation would be "1D+01" or something
01042   // similar.
01043   int i;
01044   for (i=1; i<lenm2; i++) {
01045     // Check for a 'D' character.  If we find one, then we look to see that the
01046     // preceding character is numeric, that the following character is +/-,
01047     // and that the character following +/- is also numeric.  If all of the
01048     // necessary conditions are true then we replace the 'D' with an 'E'.
01049     // The strict checking should allow us to safely convert entire strings
01050     // of characters where there are also integers, and other kinds of strings
01051     // on the same line.
01052     if (ftocstr[i] == 'D' &&
01053         ((ftocstr[i+1] == '-') || (ftocstr[i+1] == '+')) &&
01054         ((ftocstr[i-1] >= '0') && (ftocstr[i-1] <= '9')) &&
01055         ((ftocstr[i+2] >= '0') && (ftocstr[i+2] <= '9'))) {
01056       ftocstr[i] = 'E';
01057       convcnt++;
01058     }
01059   }
01060 
01061   return convcnt;
01062 }
01063 
01064 
01065 
01066 
01067 /*****************************************************
01068  *
01069  *   Debugging functions
01070  *
01071  *****************************************************/
01072 
01073 /* Print the current line in the format
01074  * "HERE) <contents of line>".
01075  * Doesn't advance the file pointer.
01076  */
01077 static void thisline(FILE *file) {
01078   char buffer[BUFSIZ];
01079   long filepos;
01080   filepos = ftell(file);
01081   if (!fgets(buffer, sizeof(buffer), file)) {
01082     if (feof(file)) printf("HERE) EOF\n");
01083     else printf("HERE) ????\n");
01084     return;
01085   }
01086   printf("HERE) %s\n", buffer);
01087   fseek(file, filepos, SEEK_SET);
01088 }
01089 
01090 
01091 /* Print all lines up to and including the next line that contains
01092  * non-whitespace characters. The file pointer will be restored to
01093  * its previous position. 
01094  * Output format: "HERE) <contents of line>".
01095  * If the end of file has been reached print "HERE) EOF".
01096  */
01097 static void whereami(FILE *file) {
01098   char buffer[BUFSIZ];
01099   char *line;
01100   long filepos;
01101   filepos = ftell(file);
01102   do {
01103     if (!fgets(buffer, sizeof(buffer), file)) {
01104       if (feof(file)) printf("HERE) EOF\n");
01105       else printf("HERE) ????\n");
01106       return;
01107     }
01108     line = trimleft(buffer);
01109     printf("HERE) %s\n", buffer);
01110   } while (!strlen(line));
01111 
01112   fseek(file, filepos, SEEK_SET);
01113 }
01114 
01115 
01116 #endif

Generated on Thu Mar 28 03:08:18 2024 for VMD Plugins (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002