← Back to branch summary

~ubuntu-branches/ubuntu/quantal/psicode/quantal

« back to all changes in this revision

Viewing changes to src/bin/input/read_zmat.cc

  • Committer: Bazaar Package Importer
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2009-02-23 00:12:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20090223001202-rutldoy3dimfpesc
Tags: 3.4.0-1
* New upstream release.

[ Michael Banck ]
* debian/patches/01_DESTDIR.dpatch: Refreshed.
* debian/patches/02_FHS.dpatch: Removed, applied upstream.
* debian/patches/03_debian_docdir: Likewise.
* debian/patches/04_man.dpatch: Likewise.
* debian/patches/06_466828_fix_gcc_43_ftbfs.dpatch: Likewise.
* debian/patches/07_464867_move_executables: Fixed and refreshed.
* debian/patches/00list: Adjusted.
* debian/control: Improved description.
* debian/patches-held: Removed.
* debian/rules (install/psi3): Do not ship the ruby bindings for now.

[ Daniel Leidert ]
* debian/rules: Fix txtdir via DEB_MAKE_INSTALL_TARGET.
* debian/patches/01_DESTDIR.dpatch: Refreshed.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/bin/input/read_zmat.cc
 
1
/*! \file
 
2
    \ingroup INPUT
 
3
    \brief Enter brief description of file here 
 
4
*/
 
5
#define EXTERN
 
6
#include <cstdio>
 
7
#include <cstdlib>
 
8
#include <cctype>
 
9
#include <cmath>
 
10
#include <cstring>
 
11
#include <libciomr/libciomr.h>
 
12
#include <libipv1/ip_lib.h>
 
13
#include <libqt/qt.h>
 
14
#include "input.h"
 
15
#include "global.h"
 
16
#include "defines.h"
 
17
 
 
18
namespace psi { namespace input {
 
19
 
 
20
/* read_zmat(): parse the zmat = () keyword in the input file.
 
21
**
 
22
** Written by EFV, original date unknown
 
23
**
 
24
** Modified 8/03 by TDC to allow a simpler z-matrix format without all the 
 
25
** "extra" parentheses.  This should allow users to cut-and-paste z-matrices from
 
26
** other programs for use in PSI.
 
27
**
 
28
** The following H2O examples should both be valid:
 
29
**
 
30
** %original format
 
31
** zmat = (
 
32
**   (o)
 
33
**   (h 1 r)
 
34
**   (h 1 r 2 a)
 
35
** )
 
36
**
 
37
** %new format (8/03)
 
38
** zmat = (
 
39
**   o
 
40
**   h 1 r
 
41
**   h 1 r 2 a
 
42
** )
 
43
*/
 
44
 
 
45
/*this is used to hold ZVARS info*/
 
46
struct definition{
 
47
      char* variable;
 
48
      double value;
 
49
};
 
50
 
 
51
/*zmat parser, code is found at bottom of this file*/
 
52
void parse_zmat(int i, int position, double *value,struct definition *array, 
 
53
                int num_vals, int zval_exist, char *zmat_lbl);
 
54
void parse_zmat_simple(int i, int position, double *value,struct definition *array, 
 
55
                       int num_vals, int zval_exist, char *zmat_lbl);
 
56
 
 
57
void read_zmat()
 
58
{
 
59
  int i, j, k, m, a, b, c, errcod, value_set, zvar_exist;
 
60
  char *buffer;
 
61
  int num_vals, entry_length, atomcount, all_atomcount;
 
62
  int A, B, C, D, f;
 
63
  int linearOn = 1;     /* Flag indicating a linear fragment in the beginning of the Z-matrix */
 
64
  double rAB, rBC, rCD, thetaABC, thetaBCD, phiABCD, val, norm1, norm2;
 
65
  double cosABC, sinABC, cosBCD, sinBCD, cosABCD, sinABCD;
 
66
  double eAB[3], eBC[3], ex[3], ey[3], ez[3];
 
67
  double Z = 0.0;
 
68
  char *name, **zvar_lbl, **zmat_lbl;
 
69
  int zmat_len, simple_zmat, simple_zvars;
 
70
 
 
71
  /*these are passed to parsing function to tell it what position in the row to look at*/
 
72
  int bond=2, angle=4, tors=6;
 
73
  struct definition *def_arr;
 
74
 
 
75
  frag_num_atoms = (int *) malloc(nfragments*sizeof(int));
 
76
  frag_num_allatoms = (int *) malloc(nfragments*sizeof(int));
 
77
  frag_atom = (int *) malloc(nfragments*sizeof(int));
 
78
  frag_allatom = (int *) malloc(nfragments*sizeof(int));
 
79
  nref_per_fragment = (int *) malloc(nfragments*sizeof(int));
 
80
  ref_pts_lc = (double ***) malloc(nfragments*sizeof(double **));
 
81
  num_atoms = 0;
 
82
  num_allatoms = 0;
 
83
      
 
84
  zvar_lbl = (char **) malloc(nfragments*sizeof(char *));
 
85
  zmat_lbl = (char **) malloc(nfragments*sizeof(char *));
 
86
 
 
87
  for (f=0; f<nfragments; ++f) { /* loop over fragments */
 
88
    zvar_lbl[f] = (char *) malloc(10*sizeof(char)); 
 
89
    zmat_lbl[f] = (char *) malloc(10*sizeof(char)); 
 
90
    if (f == 0) {
 
91
      sprintf(zvar_lbl[f],"ZVARS");
 
92
      sprintf(zmat_lbl[f],"ZMAT");
 
93
    }
 
94
    else { 
 
95
      sprintf(zvar_lbl[f],"ZVARS%d",f+1);
 
96
      sprintf(zmat_lbl[f],"ZMAT%d",f+1);
 
97
    }
 
98
 
 
99
    /* must first determine "type" of z-matrix input: nested or simple */
 
100
    simple_zmat = 1;
 
101
    ip_count(zmat_lbl[f], &zmat_len,0);
 
102
    entry_length = 0;
 
103
    for(i=0; i < zmat_len; i++) {  
 
104
      ip_count(zmat_lbl[f], &entry_length,1,i);
 
105
      if(entry_length > 1) {
 
106
        simple_zmat = 0;
 
107
      }
 
108
    }
 
109
 
 
110
    /* Read number of lines and count atoms in ZMAT */
 
111
    frag_num_atoms[f] = frag_num_allatoms[f] = 0;
 
112
    if(simple_zmat) {
 
113
      if(zmat_len >= 1) {
 
114
        frag_num_allatoms[f] = 1;
 
115
        errcod = ip_string(zmat_lbl[f],&buffer,1,0);
 
116
        if(errcod != IPE_OK) punt("Problem with the Z-matrix");
 
117
        if(strcmp(buffer,"X")) { free(buffer); frag_num_atoms[f]++; }
 
118
      }
 
119
      if(zmat_len >= 4) {
 
120
        frag_num_allatoms[f] = 2;
 
121
        errcod = ip_string(zmat_lbl[f],&buffer,1,1);
 
122
        if(errcod != IPE_OK) punt("Problem with the Z-matrix");
 
123
        if(strcmp(buffer,"X")) { free(buffer); frag_num_atoms[f]++; }
 
124
      }
 
125
      if(zmat_len >= 9) {
 
126
        frag_num_allatoms[f] = 3;
 
127
        errcod = ip_string(zmat_lbl[f],&buffer,1,4);
 
128
        if(errcod != IPE_OK) punt("Problem with the Z-matrix");
 
129
        if(strcmp(buffer,"X")) { free(buffer); frag_num_atoms[f]++; }
 
130
      }
 
131
      if(zmat_len > 9) { 
 
132
        frag_num_allatoms[f] = 3 + (zmat_len-9)/7;
 
133
        if((zmat_len-9)%7 != 0) punt("Error in z-matrix input!");
 
134
  
 
135
        for(i=0; i < frag_num_allatoms[f]-3; i++) {
 
136
          errcod = ip_string(zmat_lbl[f],&buffer,1,9+(i*7));
 
137
          if(errcod != IPE_OK) punt("Problem with the Z-matrix");
 
138
          if(strcmp(buffer,"X")) { free(buffer); frag_num_atoms[f]++; }
 
139
        }
 
140
      }
 
141
    }
 
142
    else { /* nested ZMAT */
 
143
      ip_count(zmat_lbl[f],&frag_num_allatoms[f],0);
 
144
      if (frag_num_allatoms[f] == 0)
 
145
        punt("A ZMAT is empty!");
 
146
      for(i=0;i<frag_num_allatoms[f];i++){
 
147
        errcod = ip_string(zmat_lbl[f],&buffer,2,i,0);
 
148
        if (errcod != IPE_OK)
 
149
          punt("Problem with the Z-matrix");
 
150
        if (strcmp(buffer,"X"))
 
151
          frag_num_atoms[f]++;
 
152
        free(buffer);
 
153
      }
 
154
    }
 
155
    if (frag_num_atoms[f] == 0) punt("A Z-matrix has no non-dummy atoms!");
 
156
    num_atoms += frag_num_atoms[f];
 
157
    num_allatoms += frag_num_allatoms[f];
 
158
  }
 
159
 
 
160
  if(num_atoms > MAXATOM) punt("There are more atoms than allowed!");
 
161
 
 
162
  frag_allatom[0] = frag_atom[0] = 0;
 
163
  for (f=1; f<nfragments; ++f) {
 
164
    frag_atom[f] = frag_atom[f-1] + frag_num_atoms[f-1];
 
165
    frag_allatom[f] = frag_allatom[f-1] + frag_num_allatoms[f-1];
 
166
  }
 
167
 
 
168
  /*-----------------------
 
169
    Allocate global arrays
 
170
    -----------------------*/
 
171
  full_geom = block_matrix(num_allatoms,3);
 
172
  geometry = (double **) malloc(num_atoms*sizeof(double *));
 
173
  atom_dummy = (int *) malloc(sizeof(int)*num_allatoms);
 
174
  /* see chkpt.h for info about z_entry structure */
 
175
  z_geom = (struct z_entry *) malloc(sizeof(struct z_entry)*num_allatoms); 
 
176
  element = (char **) malloc(sizeof(char *)*num_atoms);
 
177
  full_element = (char **) malloc(sizeof(char *)*num_allatoms);
 
178
  elemsymb_charges = init_array(num_atoms);
 
179
 
 
180
  atomcount = 0;
 
181
  all_atomcount = 0;
 
182
  for (f=0; f<nfragments; ++f) {
 
183
 
 
184
    /*need this to avoid seg fault if character string used but no zvar given*/
 
185
    zvar_exist = ip_exist(zvar_lbl[f],0);    
 
186
    errcod = 0;
 
187
    if (zvar_exist) { /* read in ZVARS */
 
188
      errcod += ip_count(zvar_lbl[f],&num_vals,0);
 
189
      simple_zvars = 1;
 
190
      entry_length = 0;
 
191
      for (i=0; i<num_vals; i++) {
 
192
        ip_count(zvar_lbl[f], &entry_length,1,i);
 
193
        if(entry_length>1) simple_zvars = 0;
 
194
      }
 
195
      if(simple_zvars) {
 
196
        if(num_vals%2) punt("Problem with number of simple ZVARS entries.");
 
197
        num_vals /= 2;
 
198
        def_arr = (struct definition *) malloc( num_vals * sizeof(struct definition) );
 
199
        for(i=0; i < num_vals; i++) {
 
200
          errcod += ip_string(zvar_lbl[f], &def_arr[i].variable, 1, 2*i);
 
201
          errcod += ip_data(zvar_lbl[f], "%lf", &def_arr[i].value, 1, 2*i+1);
 
202
        }
 
203
      }
 
204
      else { /* nested ZVARS */
 
205
        def_arr = (struct definition *) malloc( num_vals * sizeof(struct definition) );
 
206
        for(i=0;i<num_vals;++i) {
 
207
          errcod += ip_string(zvar_lbl[f],&def_arr[i].variable,2,i,0);
 
208
          errcod += ip_data(zvar_lbl[f],"%lf",&def_arr[i].value,2,i,1);
 
209
        }
 
210
      }
 
211
      if(errcod > 0)
 
212
        punt("Problem parsing ZVARS");
 
213
    }
 
214
 
 
215
    for (i=0; i<frag_num_allatoms[f]; ++i) {
 
216
      /* Process a line of ZMAT and write to z_geom array */
 
217
      if (!simple_zmat) {
 
218
        entry_length = 0;
 
219
        ip_count(zmat_lbl[f],&entry_length,1,i);
 
220
        if ( ((i == 0) && (entry_length != 1)) ||
 
221
             ((i == 1) && (entry_length != 3)) ||
 
222
             ((i == 2) && (entry_length != 5)) ||
 
223
             ((i  > 2) && (entry_length != 7)) ) {
 
224
          fprintf(outfile,"  Line %d of ZMAT has a wrong number of entries.\n",i+1);
 
225
          punt("Invalid ZMAT");
 
226
        }
 
227
      }
 
228
 
 
229
      if (i == 0) {                                     /*      1st atom */
 
230
        full_geom[all_atomcount][0] = 0.0;
 
231
        full_geom[all_atomcount][1] = 0.0;
 
232
        full_geom[all_atomcount][2] = 0.0;
 
233
 
 
234
z_geom[all_atomcount].bond_atom= z_geom[all_atomcount].angle_atom= z_geom[all_atomcount].tors_atom= -1;
 
235
z_geom[all_atomcount].bond_opt = z_geom[all_atomcount].angle_opt = z_geom[all_atomcount].tors_opt = -1;
 
236
z_geom[all_atomcount].bond_val = z_geom[all_atomcount].angle_val = z_geom[all_atomcount].tors_val = -999.9;
 
237
z_geom[all_atomcount].bond_label[0]= z_geom[all_atomcount].angle_label[0]= z_geom[all_atomcount].tors_label[0]=' ';
 
238
      }
 
239
 
 
240
      else if (i == 1) {                                        /*      2nd atom */
 
241
        if(!simple_zmat) {
 
242
          ip_data(zmat_lbl[f],"%d",&a,2,i,1);
 
243
          parse_zmat(i,bond,&rAB,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
244
        }
 
245
        else {
 
246
          ip_data(zmat_lbl[f], "%d", &a, 1, 2);
 
247
          parse_zmat_simple(i,bond,&rAB,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
248
        }
 
249
        
 
250
        z_geom[all_atomcount].bond_atom = a + frag_allatom[f];
 
251
        z_geom[all_atomcount].angle_atom= z_geom[all_atomcount].tors_atom = -1;
 
252
        z_geom[all_atomcount].angle_opt = z_geom[all_atomcount].tors_opt  = -1; 
 
253
        z_geom[all_atomcount].bond_val  = rAB * conv_factor;
 
254
        z_geom[all_atomcount].angle_val = z_geom[all_atomcount].tors_val  = -999.9;
 
255
 
 
256
        if (rAB < ZERO_BOND_DISTANCE)
 
257
          punt("Invalid bond length in atom 2.");
 
258
        full_geom[all_atomcount][0] = 0.0;
 
259
        full_geom[all_atomcount][1] = 0.0;
 
260
        full_geom[all_atomcount][2] = rAB;
 
261
      }
 
262
 
 
263
      else if (i == 2) {                                        /*      3rd atom */
 
264
        if(!simple_zmat) {
 
265
          ip_data(zmat_lbl[f],"%d",&a,2,i,1);
 
266
          ip_data(zmat_lbl[f],"%d",&b,2,i,3);
 
267
        }
 
268
        else {
 
269
          ip_data(zmat_lbl[f], "%d", &a, 1, 5);
 
270
          ip_data(zmat_lbl[f], "%d", &b, 1, 7);
 
271
        }
 
272
        if ( ((a == 2) && (b == 1)) ||
 
273
             ((a == 1) && (b == 2)) ) {
 
274
 
 
275
          if(!simple_zmat) parse_zmat(i,bond,&rBC,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
276
          else parse_zmat_simple(i,bond,&rBC,def_arr,num_vals,zvar_exist, zmat_lbl[f]);
 
277
            
 
278
          if (rBC <= ZERO_BOND_DISTANCE) {
 
279
            fprintf(outfile,"  Invalid bond length in atom 3.\n");
 
280
            punt("Invalid ZMAT");
 
281
          }
 
282
 
 
283
          if(!simple_zmat) parse_zmat(i,angle,&thetaABC,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
284
          else parse_zmat_simple(i,angle,&thetaABC,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
285
        
 
286
          //if (thetaABC <= ZERO_BOND_ANGLE) { // does not permit 0 bond angle with atoms 1-2-3
 
287
          if (thetaABC < 0.0) {
 
288
            fprintf(outfile,"  Invalid bond angle in atom 3.\n");
 
289
            punt("Invalid ZMAT");
 
290
          }
 
291
          z_geom[all_atomcount].bond_atom = a + frag_allatom[f];
 
292
          z_geom[all_atomcount].bond_val  = rBC * conv_factor;
 
293
          z_geom[all_atomcount].angle_atom= b + frag_allatom[f];
 
294
          z_geom[all_atomcount].angle_val = thetaABC;
 
295
          z_geom[all_atomcount].tors_atom = -1;
 
296
          z_geom[all_atomcount].tors_val  = -999.9;
 
297
          z_geom[all_atomcount].tors_opt  = -1;
 
298
 
 
299
          thetaABC = thetaABC*M_PI/180.0;
 
300
           
 
301
              if (a == 2) {                             /*      ABC case */
 
302
            full_geom[all_atomcount][0] = full_geom[a+frag_allatom[f]-1][0] + rBC*sin(thetaABC);
 
303
            full_geom[all_atomcount][1] = full_geom[a+frag_allatom[f]-1][1];
 
304
            full_geom[all_atomcount][2] = full_geom[a+frag_allatom[f]-1][2] - rBC*cos(thetaABC);
 
305
              }
 
306
          else {                                        /*      BAC case */
 
307
            full_geom[all_atomcount][0] = full_geom[a+frag_allatom[f]-1][0] + rBC*sin(thetaABC);
 
308
            full_geom[all_atomcount][1] = full_geom[a+frag_allatom[f]-1][1];
 
309
            full_geom[all_atomcount][2] = full_geom[a+frag_allatom[f]-1][2] + rBC*cos(thetaABC);
 
310
          }
 
311
        }
 
312
        else {
 
313
              fprintf(outfile,"  Problem in atom 3 in zmat.\n");
 
314
              punt("Invalid ZMAT");;
 
315
        }
 
316
      }
 
317
 
 
318
      else { 
 
319
        if(!simple_zmat) {
 
320
          ip_data(zmat_lbl[f],"%d",&c,2,i,1);
 
321
          ip_data(zmat_lbl[f],"%d",&b,2,i,3);
 
322
          ip_data(zmat_lbl[f],"%d",&a,2,i,5);
 
323
        }
 
324
        else {
 
325
          ip_data(zmat_lbl[f],"%d",&c,1,9+((i-3)*7)+1);
 
326
          ip_data(zmat_lbl[f],"%d",&b,1,9+((i-3)*7)+3);
 
327
          ip_data(zmat_lbl[f],"%d",&a,1,9+((i-3)*7)+5);
 
328
        }
 
329
        a -= 1; b -= 1; c -= 1;
 
330
  
 
331
        if ( (a == b) || (b == c) || (a == c) ||
 
332
             (a >= i) || (b >= i) || (c >= i) ) {
 
333
          fprintf(outfile,"  Problem in atom %d of zmat.\n",i);
 
334
          punt("Invalid ZMAT");
 
335
        }
 
336
  
 
337
        if(!simple_zmat) parse_zmat(i,bond,&rCD,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
338
        else parse_zmat_simple(i,bond,&rCD,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
339
          
 
340
        if (rCD <= ZERO_BOND_DISTANCE) {
 
341
          fprintf(outfile,"  Invalid bond length in atom %d.\n",i+1);
 
342
          punt("Invalid ZMAT");
 
343
        }
 
344
  
 
345
        if(!simple_zmat) parse_zmat(i,angle,&thetaBCD,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
346
        else parse_zmat_simple(i,angle,&thetaBCD,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
347
          
 
348
        if (thetaBCD <= ZERO_BOND_ANGLE) {
 
349
          fprintf(outfile,"  Invalid bond angle in atom %d.\n",i+1);
 
350
          punt("Invalid ZMAT");
 
351
        }
 
352
  
 
353
        if(!simple_zmat) parse_zmat(i,tors,&phiABCD,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
354
        else parse_zmat_simple(i,tors,&phiABCD,def_arr,num_vals,zvar_exist,zmat_lbl[f]);
 
355
             
 
356
        z_geom[all_atomcount].bond_atom  = c+1+frag_allatom[f];
 
357
        z_geom[all_atomcount].angle_atom = b+1+frag_allatom[f];
 
358
        z_geom[all_atomcount].tors_atom  = a+1+frag_allatom[f];
 
359
        z_geom[all_atomcount].bond_val   = rCD * conv_factor;
 
360
        z_geom[all_atomcount].angle_val  = thetaBCD;
 
361
        z_geom[all_atomcount].tors_val   = phiABCD;
 
362
        
 
363
        thetaBCD = thetaBCD * M_PI/180.0; 
 
364
        phiABCD = phiABCD * M_PI/180.0;
 
365
 
 
366
        /* If you want to have linear fragment defined in 
 
367
            the beginning of the Z-matrix - fine, but you still have to
 
368
            supply the third "dummy" atom and a value for the dihedral angle*/
 
369
        a += frag_allatom[f]; /* switch to full index */
 
370
        b += frag_allatom[f];
 
371
        c += frag_allatom[f];
 
372
 
 
373
        if ( ((full_geom[all_atomcount-1][0] - LINEAR_CUTOFF) < 0.0) && linearOn ) {
 
374
          if ((full_geom[c][2] - full_geom[b][2]) > 0.0) {
 
375
            full_geom[all_atomcount][0] = rCD*sin(thetaBCD);
 
376
            full_geom[all_atomcount][1] = 0.0;
 
377
            full_geom[all_atomcount][2] = full_geom[c][2] - rCD*cos(thetaBCD);
 
378
          }
 
379
          else {
 
380
            full_geom[all_atomcount][0] = rCD*sin(thetaBCD);
 
381
            full_geom[all_atomcount][1] = 0.0;
 
382
            full_geom[all_atomcount][2] = full_geom[c][2] + rCD*cos(thetaBCD);
 
383
          }
 
384
        }
 
385
 
 
386
        else { /* code for nonlinear ABC fragment */
 
387
          /* Set "linearity" Flag to false */
 
388
          linearOn = 0;
 
389
           
 
390
              /* Note : x,y,z - unit vectors defining a coordinate 
 
391
              system associated with atom C; BC defines z, ABC defines the xz-plane */
 
392
 
 
393
          unit_vec(full_geom[b],full_geom[a],eAB);      /* U1 in the old zmat */
 
394
          unit_vec(full_geom[c],full_geom[b],ez);       /* U2 in the old zmat */
 
395
          cosABC = -dot_prod(ez,eAB);   /* ez is essentially eBC */
 
396
          sinABC = sqrt(1 - (cosABC * cosABC) );
 
397
           
 
398
              if ( (sinABC - LINEAR_CUTOFF) < 0.0 ) {
 
399
                fprintf(outfile,"  Dihedral angle in line %d is defined with respect to a linear fragment.\n",i+1);
 
400
                punt("Invalid ZMAT");
 
401
              }
 
402
           
 
403
              cross_prod(eAB,ez,ey);            /* ey is U3 in the old zmat */
 
404
              for(j=0;j<3;j++)                  /* normalization of ey */
 
405
                ey[j] /= sinABC;
 
406
              cross_prod(ey,ez,ex);             /* ex is U4 in the old zmat */
 
407
 
 
408
              /* Intermediates - to avoid calling 
 
409
                 trig. functions multiple times */
 
410
              cosBCD = cos(thetaBCD);
 
411
              sinBCD = sin(thetaBCD);
 
412
              cosABCD = cos(phiABCD);
 
413
              sinABCD = sin(phiABCD);
 
414
           
 
415
              for (j=0;j<3;j++)
 
416
                full_geom[all_atomcount][j] = full_geom[c][j] + rCD * 
 
417
                  ( - ez[j] * cosBCD + ex[j] * sinBCD * cosABCD + ey[j] * sinBCD * sinABCD );
 
418
        }
 
419
      }
 
420
 
 
421
      if(!simple_zmat) errcod = ip_string(zmat_lbl[f],&buffer,2,i,0);
 
422
      else {
 
423
        if(i==0) errcod = ip_string(zmat_lbl[f],&buffer,1,0);
 
424
        else if(i==1) errcod = ip_string(zmat_lbl[f],&buffer,1,1);
 
425
        else if(i==2) errcod = ip_string(zmat_lbl[f],&buffer,1,4);
 
426
        else if(i>2) errcod = ip_string(zmat_lbl[f],&buffer,1,9+((i-3)*7));
 
427
      }
 
428
 
 
429
      if (strcmp(buffer,"X")) {
 
430
        atom_num(buffer, &Z);
 
431
        free(buffer);
 
432
        elemsymb_charges[atomcount] = Z;
 
433
        element[atomcount] = elem_name[(int)Z];
 
434
        full_element[all_atomcount] = elem_name[(int)Z];
 
435
        geometry[atomcount] = full_geom[all_atomcount];
 
436
        atom_dummy[all_atomcount] = 0;
 
437
        ++atomcount;
 
438
      }
 
439
      else {
 
440
        static const char* dummy_elem_name = "X";
 
441
        full_element[all_atomcount] = const_cast<char*>(dummy_elem_name);
 
442
        free(buffer);
 
443
        atom_dummy[all_atomcount] = 1;
 
444
      }
 
445
      ++all_atomcount;
 
446
    }
 
447
  }
 
448
 
 
449
  for(i=0;i<num_allatoms;++i) {
 
450
    full_geom[i][0] *=conv_factor;
 
451
    full_geom[i][1] *=conv_factor;
 
452
    full_geom[i][2] *=conv_factor;
 
453
  }
 
454
 
 
455
  read_charges();
 
456
 
 
457
  fprintf(outfile,"Coordinates after reading z-matrices\n");
 
458
  print_mat(geometry, num_atoms, 3, outfile);
 
459
 
 
460
  return;
 
461
}
 
462
 
 
463
 
 
464
 
 
465
/*******************************************************************************************************************
 
466
 
 
467
  parse_zmat
 
468
 
 
469
  function used to parse z-matrix input
 
470
 
 
471
  position should be 2 for bond, 4 for angle, 6 for torsion (corresponds to postion of value/variable in zmat row)
 
472
  since I didn't want to make more stuff global this is pretty ugly, but it cleans things up in read_zmat();
 
473
  ******************************************************************************************************************/
 
474
 
 
475
void parse_zmat(int i, int position, double *value, struct definition
 
476
*array, int num_vals, int zvar_exist, char *zmat_lbl) {
 
477
 
 
478
   int j, a, value_set;
 
479
   double r;
 
480
   char *temp_string, *dollar;
 
481
 
 
482
   value_set = 0;
 
483
   ip_string(zmat_lbl,&temp_string,2,i,position);
 
484
   if( !isalnum( temp_string[0] ) &&
 
485
       temp_string[0] != '-' &&
 
486
       temp_string[0] != '+' &&
 
487
       temp_string[0] != '.')
 
488
     punt("z-matrix entry doesn't begin with letter or number");
 
489
   dollar = strchr(temp_string,'$');
 
490
   if( dollar != NULL ) {
 
491
     if(position == 2) 
 
492
       z_geom[i].bond_opt = 1;
 
493
     else if(position == 4)
 
494
       z_geom[i].angle_opt = 1;
 
495
     else if(position == 6)
 
496
       z_geom[i].tors_opt = 1;
 
497
     *dollar = '\0';
 
498
   }
 
499
   else {
 
500
     if(position == 2) 
 
501
       z_geom[i].bond_opt = 0;
 
502
     else if(position == 4)
 
503
       z_geom[i].angle_opt = 0;
 
504
     else if(position == 6)
 
505
       z_geom[i].tors_opt = 0;
 
506
   }
 
507
   
 
508
   if( isdigit( temp_string[0] ) ||
 
509
       ( (temp_string[0] == '-') && isdigit(temp_string[1]) )  ||
 
510
       ( (temp_string[0] == '+') && isdigit(temp_string[1]) )  ||
 
511
       ( (temp_string[0] == '.') && isdigit(temp_string[1]) ) ) { 
 
512
     *value = atof( temp_string );
 
513
     ++value_set;
 
514
   }
 
515
   if( isalpha( temp_string[0] ) && zvar_exist ) {
 
516
     for(j=0;j<num_vals;++j) 
 
517
       if( !strcmp(temp_string,array[j].variable) ) {
 
518
         *value = array[j].value;
 
519
         ++value_set;
 
520
       }
 
521
   }
 
522
   if(value_set != 1 ) {
 
523
     fprintf(outfile,"  Problem with variable definition from line %i\n",i+1);
 
524
     punt("Invalid ZMAT");
 
525
   }
 
526
   
 
527
   {
 
528
     const char* tstring = isalpha(temp_string[0]) ? temp_string : "";
 
529
     if(position == 2) strncpy( z_geom[i].bond_label, tstring, CHKPT_ZMAT_LABEL_LEN);
 
530
     if(position == 4) strncpy( z_geom[i].angle_label, tstring, CHKPT_ZMAT_LABEL_LEN);
 
531
     if(position == 6) strncpy( z_geom[i].tors_label, tstring, CHKPT_ZMAT_LABEL_LEN);
 
532
   }
 
533
   
 
534
   free(temp_string);   
 
535
   
 
536
   return;
 
537
}
 
538
 
 
539
void parse_zmat_simple(int i, int position, double *value, struct definition *array,
 
540
                       int num_vals, int zvar_exist, char *zmat_lbl)
 
541
{
 
542
  int j, value_set;
 
543
  char *temp_string, *dollar;
 
544
 
 
545
  value_set=0;
 
546
 
 
547
  if(i==1) ip_string(zmat_lbl, &temp_string, 1, 3);
 
548
  else if(i==2) ip_string(zmat_lbl, &temp_string, 1, 4+position);
 
549
  else ip_string(zmat_lbl, &temp_string, 1, 9+((i-3)*7)+position);
 
550
 
 
551
  if(!isalnum(temp_string[0]) && temp_string[0] != '-' && 
 
552
     temp_string[0] != '+' && temp_string[0] != '.')
 
553
    punt("Z-matrix entry doesn't begin with a letter or number.");
 
554
 
 
555
  dollar = strchr(temp_string, '$');
 
556
  if(dollar != NULL) {
 
557
    if(position == 2) z_geom[i].bond_opt = 1;
 
558
    else if(position == 4) z_geom[i].angle_opt = 1;
 
559
    else if(position == 6) z_geom[i].tors_opt = 1;
 
560
    *dollar = '\0';
 
561
  }
 
562
  else {
 
563
    if(position == 2) z_geom[i].bond_opt = 0;
 
564
    else if(position == 4) z_geom[i].angle_opt = 0;
 
565
    else if(position == 6) z_geom[i].tors_opt = 0;
 
566
   }
 
567
 
 
568
  if ( isdigit(temp_string[0]) || 
 
569
       ( (temp_string[0] == '-') && isdigit(temp_string[1]) )  ||
 
570
       ( (temp_string[0] == '+') && isdigit(temp_string[1]) )  ||
 
571
       ( (temp_string[0] == '.') && isdigit(temp_string[1]) ) ) { 
 
572
    *value = atof(temp_string);
 
573
    ++value_set;
 
574
  }
 
575
 
 
576
  if(isalpha(temp_string[0]) && zvar_exist) {
 
577
    for(j=0; j < num_vals; j++)
 
578
      if(!strcmp(temp_string, array[j].variable)) {
 
579
        *value = array[j].value;
 
580
        ++value_set;
 
581
      }
 
582
  }
 
583
  if(!value_set) {
 
584
    fprintf(outfile, "  Problem with variable definition for atom %d.\n", i+1);
 
585
    punt("Invalid Z-matrix.");
 
586
  }
 
587
 
 
588
  {
 
589
    const char* tstring = isalpha(temp_string[0]) ? temp_string : "";
 
590
    if(position == 2) strncpy( z_geom[i].bond_label, tstring, CHKPT_ZMAT_LABEL_LEN);
 
591
    if(position == 4) strncpy( z_geom[i].angle_label, tstring, CHKPT_ZMAT_LABEL_LEN);
 
592
    if(position == 6) strncpy( z_geom[i].tors_label, tstring, CHKPT_ZMAT_LABEL_LEN);
 
593
  }
 
594
 
 
595
  free(temp_string);
 
596
 
 
597
  return;
 
598
}
 
599
 
 
600
}} // namespace psi::input