ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/readz.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 11 months ago) by tjod
File size: 31443 byte(s)
Log Message:
test

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5     #include "substr.h"
6     #include "utility.h"
7    
8     #define MAIN 0L
9     #define SUBSTRUCTURE 1L
10    
11     struct t_mopdata {
12     int iaopt, ibopt, idopt;
13     int scf1, ef, force,precise, grad, vectors;
14     int Hamiltonian;
15     char keywords[90],title[80],title1[80];
16     } mopdata;
17     EXTERN struct ElementType { char symbol[3];
18     int atomnum;
19     float weight, covradius, vdwradius;
20     int s,p,d,f, type;
21     } Elements[];
22    
23     void distjm(int,int, float *);
24     void angljm(int, int, int, float *);
25     void type(void);
26     void hdel(int);
27     void hadd(void);
28     double dihdrl(int, int, int, int);
29     void generate_bonds(void);
30     void readz(int, int);
31     void read_arcfile(int, int);
32     void mopaco(int,int);
33     int make_atom(int,float,float,float,char *);
34     void message_alert(char *, char *);
35     void z2xyz(int, float *, float *, float *,int *, int *, int *,float *, float *, float *);
36     int is_metal(char *);
37     FILE * fopen_path ( char * , char * , char * ) ;
38     int FetchRecord(FILE *, char *);
39    
40     float *bl, *alpha, *beta, *mcharge;
41     float *xtmp, *ytmp, *ztmp;
42    
43     void readz(int mndotype, int isubred)
44     {
45     FILE *infile;
46    
47     char inputline[251], dummy[3];
48     int i,j,ibotptr,niatom,idummy,newatom;
49     int ia1, ia2, ia3;
50     int inumarg, iline, ntemp;
51     int icart, isym;
52     float wtt,blen, bangle, dihed;
53     int *na, *nb, *nc, *mtype;
54     int metalcount;
55    
56     if (isubred == MAIN)
57     {
58     ibotptr = 0;
59     }else
60     {
61     ibotptr = natom;
62     substr.istract[substr.icurstr] = TRUE;
63     }
64    
65     infile = fopen_path(Openbox.path,Openbox.fname,"r");
66    
67     FetchRecord(infile,inputline); // keyword
68     FetchRecord(infile,inputline); // title
69     strncpy(Struct_Title,inputline,60);
70     FetchRecord(infile,inputline); // comment
71    
72     ntemp= 0;
73     metalcount = 0;
74     icart = FALSE;
75     isym = TRUE;
76     while ( FetchRecord(infile,inputline))
77     {
78     inumarg = sscanf(inputline,"%s %f %d %f %d %f %d %d %d %d %f ",
79     dummy, &blen, &j, &bangle, &j, &dihed, &j, &ia1, &ia2,
80     &ia3, &wtt);
81     if (inumarg == 1 && ntemp > 1)
82     break;
83     if (inumarg <= 0)
84     break;
85     if (inumarg > 0)
86     ntemp++;
87     if (ntemp == 1 && ia1 == 99)
88     icart = TRUE;
89     if (ntemp == 1)
90     {
91     if (isdigit(dummy[0]) != 0)
92     isym = FALSE;
93     }
94     }
95    
96     fclose(infile);
97    
98     bl = vector(0, ntemp);
99     alpha = vector(0,ntemp);
100     beta = vector(0,ntemp);
101     mcharge = vector(0,ntemp);
102     mtype = ivector(0,ntemp);
103     na = ivector(0,ntemp);
104     nb = ivector(0,ntemp);
105     nc = ivector(0,ntemp);
106     xtmp = vector(0,ntemp);
107     ytmp = vector(0,ntemp);
108     ztmp = vector(0,ntemp);
109    
110     for (i = 0; i < ntemp; i++)
111     {
112     bl[i] = 0.0F;
113     alpha[i] = 0.0F;
114     beta[i] = 0.0F;
115     mtype[i] = 0;
116     na[i] = 0;
117     nb[i] = 0;
118     nc[i] = 0;
119     }
120    
121     infile = fopen_path(Openbox.path,Openbox.fname,"r");
122     FetchRecord(infile,inputline);
123     FetchRecord(infile,inputline);
124     FetchRecord(infile,inputline);
125     niatom = 0;
126     iline = 0;
127     while ( FetchRecord(infile,inputline))
128     {
129     iline++;
130     blen = 0.0F;
131     bangle = 0.0F;
132     dihed = 0.0F;
133     ia1 = 0;
134     ia2 = 0;
135     ia3 = 0;
136     inumarg = sscanf(inputline,"%s %f %d %f %d %f %d %d %d %d %f ",
137     dummy, &blen, &j, &bangle, &j, &dihed, &j, &ia1, &ia2,
138     &ia3, &mcharge[niatom]);
139    
140     if ( inumarg >= 10 && (strcmp(dummy,"0") == 0)) // termination with line of 0's
141     goto L_10;
142     else if ( inumarg >= 10)
143     {
144     bl[niatom] = blen;
145     alpha[niatom] = bangle;
146     beta[niatom] = dihed;
147     na[niatom] = ia1;
148     nb[niatom] = ia2;
149     nc[niatom] = ia3;
150     } else if (inumarg == 1 && iline > 4)
151     {
152     goto L_10;
153     } else if (iline == 1)
154     {
155     bl[niatom] = 0.0F;
156     alpha[niatom] = 0.0F;
157     beta[niatom] = 0.0F;
158     na[niatom] = 0;
159     nb[niatom] = 0;
160     nc[niatom] = 0;
161     } else if (iline == 2) // line 2
162     {
163     sscanf(inputline,"%s %f %d %d %d %d", dummy, &blen, &j, &ia1, &ia2, &ia3);
164     bl[niatom] = blen;
165     na[niatom] = ia1;
166     nb[niatom] = 0;
167     nc[niatom] = 0;
168     } else if (iline == 3) // line 3
169     {
170     sscanf(inputline,"%s %f %d %f %d %d %d %d", dummy, &blen, &j, &bangle,&j,&ia1, &ia2, &ia3);
171     bl[niatom] = blen;
172     alpha[niatom] = bangle;
173     na[niatom] = ia1;
174     nb[niatom] = ia2;
175     nc[niatom] = 0;
176     } else if (inumarg < 10 && iline > 4)
177     goto L_10;
178    
179     if (isdigit(dummy[0]) != 0) // we have atomic numbers
180     mtype[niatom] = atoi(dummy);
181     else if (strcasecmp(dummy,"XX") == 0 || dummy[0] == 'X')
182     mtype[niatom] = 99;
183     else
184     {
185     for (j=0; j < 103; j++)
186     {
187     if (strcasecmp(dummy,Elements[j].symbol) == 0)
188     {
189     mtype[niatom] = Elements[j].atomnum;
190     break;
191     }
192     }
193     }
194     idummy = mtype[niatom];
195     niatom++;
196     }
197     L_10:
198     fclose(infile);
199     if (icart == FALSE)
200     z2xyz(niatom,bl,alpha,beta,na,nb,nc,xtmp,ytmp,ztmp);
201     else
202     {
203     for (i=0; i < niatom; i++)
204     {
205     xtmp[i] = bl[i];
206     ytmp[i] = alpha[i];
207     ztmp[i] = beta[i];
208     }
209     }
210    
211     /* need to remove dummy atoms */
212     while( TRUE )
213     {
214     idummy = 0;
215     for (i=0; i < niatom; i++)
216     {
217     if (mtype[i] == 99)
218     {
219     idummy = 1;
220     for ( j = i; j < niatom-1; j++)
221     {
222     mtype[j] = mtype[j+1];
223     xtmp[j] = xtmp[j+1];
224     ytmp[j] = ytmp[j+1];
225     ztmp[j] = ztmp[j+1];
226     }
227     mtype[niatom] = 0;
228     xtmp[niatom] = 0.0F;
229     ytmp[niatom] = 0.0F;
230     ztmp[niatom] = 0.0F;
231     }
232     }
233     if (idummy != 1)
234     break;
235     }
236    
237     /* generate atoms */
238     idummy = 0;
239     for (i= 0; i < niatom; i++)
240     {
241     if (mtype[i] > 0)
242     {
243     newatom = make_atom(0,xtmp[i],ytmp[i],ztmp[i],Elements[mtype[i]-1].symbol);
244     if (isubred == 1)
245     {
246     atom[newatom].substr[0] = 0;
247     atom[newatom].substr[0] |= (1L << substr.icurstr);
248     }
249     }
250    
251     }
252    
253    
254     /* need to generate bonds using mopaco */
255     mopaco(ibotptr+1,natom);
256    
257     free_vector(bl, 0,ntemp);
258     free_vector(alpha, 0,ntemp);
259     free_vector(beta, 0,ntemp);
260     free_vector(mcharge, 0,ntemp);
261     free_ivector(na, 0,ntemp);
262     free_ivector(nb, 0,ntemp);
263     free_ivector(nc, 0,ntemp);
264     free_ivector(mtype, 0,ntemp);
265     free_vector(xtmp, 0,ntemp);
266     free_vector(ytmp, 0,ntemp);
267     free_vector(ztmp, 0,ntemp);
268     }
269     /* ========================================= */
270     void read_arcfile(int isel, int isubred)
271     {
272     FILE *infile;
273    
274     char inputline[251],atomchar[3];
275     char *p;
276     int i,j,ibotptr,niatom,idummy,newatom;
277     int ia1, ia2, ia3;
278     int inumarg, iline, ntemp, ipoint;
279     float blen, bangle, dihed;
280     int *na, *nb, *nc, *mtype;
281     float *bl, *alpha, *beta, *mcharge;
282     float *xtmp, *ytmp, *ztmp;
283    
284     if (isubred == MAIN)
285     ibotptr = 0;
286     else
287     ibotptr = natom;
288    
289     infile = fopen_path(Openbox.path,Openbox.fname,"r");
290     if (infile == NULL)
291     {
292     message_alert("Error reading arc file","Arcfile Setup");
293     return;
294     }
295     iline = 0;
296     while ( FetchRecord(infile,inputline) )
297     iline++;
298     fclose(infile);
299     infile = fopen_path(Openbox.path,Openbox.fname,"r");
300    
301     ntemp = iline-3;
302     bl = vector(0, ntemp);
303     alpha = vector(0,ntemp);
304     beta = vector(0,ntemp);
305     mcharge = vector(0,ntemp);
306     mtype = ivector(0,ntemp);
307     na = ivector(0,ntemp);
308     nb = ivector(0,ntemp);
309     nc = ivector(0,ntemp);
310     xtmp = vector(0,ntemp);
311     ytmp = vector(0,ntemp);
312     ztmp = vector(0,ntemp);
313    
314     for (i = 0; i <= ntemp; i++)
315     {
316     bl[i] = 0.0F;
317     alpha[i] = 0.0F;
318     beta[i] = 0.0F;
319     mtype[i] = 0L;
320     na[i] = 0L;
321     nb[i] = 0L;
322     nc[i] = 0L;
323     }
324    
325     ipoint = 0;
326     while ( FetchRecord(infile,inputline) )
327     {
328     p = strtok(inputline," ");
329     if (p != NULL)
330     {
331     if (strcasecmp(p,"FINAL") == 0)
332     ipoint++;
333     if (ipoint == isel)
334     goto L_20;
335     }
336     }
337    
338     L_20:
339     // start file read
340     FetchRecord(infile,inputline); // keyword
341     FetchRecord(infile,inputline); // title
342     strncpy(Struct_Title,inputline,60);
343     FetchRecord(infile,inputline); // comment
344    
345     niatom = 0;
346     iline = 0;
347     while ( FetchRecord(infile,inputline))
348     {
349     iline++;
350     blen = 0.0F;
351     bangle = 0.0F;
352     dihed = 0.0F;
353     ia1 = 0;
354     ia2 = 0;
355     ia3 = 0;
356     inumarg = sscanf(inputline,"%s %f %d %f %d %f %d %d %d %d %f ",
357     atomchar, &blen, &j, &bangle, &j, &dihed, &j, &ia1, &ia2,
358     &ia3, &mcharge[niatom]);
359    
360     if (strcmp(atomchar,"0") == 0) // end of file in Ampac
361     goto L_10;
362    
363     if ( inumarg >= 10)
364     {
365     bl[niatom] = blen;
366     alpha[niatom] = bangle;
367     beta[niatom] = dihed;
368     na[niatom] = ia1+ibotptr;
369     nb[niatom] = ia2+ibotptr;
370     nc[niatom] = ia3+ibotptr;
371     } else if (inumarg == 1 && iline > 4)
372     {
373     goto L_10;
374     } else if (inumarg == 1 && iline == 1)
375     {
376     bl[niatom] = 0.0F;
377     alpha[niatom] = 0.0F;
378     beta[niatom] = 0.0F;
379     na[niatom] = 0+ibotptr;
380     nb[niatom] = 0+ibotptr;
381     nc[niatom] = 0+ibotptr;
382    
383     } else if (inumarg == 4 && iline == 2) // line 2
384     {
385     bl[niatom] = blen;
386     na[niatom] = (int) bangle+ibotptr;
387     } else if (inumarg == 7 && iline == 3) // line 3
388     {
389     bl[niatom] = blen;
390     alpha[niatom] = bangle;
391     na[niatom] = (int) dihed+ibotptr;
392     nb[niatom] = j+ibotptr;
393     } else if (inumarg < 10 && iline > 4)
394     goto L_10;
395    
396     for (j=0; j < 103; j++)
397     {
398     if (strcasecmp(atomchar,Elements[j].symbol) == 0)
399     {
400     mtype[niatom] = Elements[j].atomnum;
401     break;
402     }
403     }
404     niatom++;
405     }
406     L_10:
407     fclose(infile);
408    
409     z2xyz(niatom,bl,alpha,beta,na,nb,nc,xtmp,ytmp,ztmp);
410    
411     /* need to remove dummy atoms */
412     while( TRUE )
413     {
414     idummy = 0;
415     for (i=0; i < niatom; i++)
416     {
417     if (mtype[i] == 99)
418     {
419     idummy = 1;
420     for ( j = i; j < niatom-1; j++)
421     {
422     mtype[j] = mtype[j+1];
423     xtmp[j] = xtmp[j+1];
424     ytmp[j] = ytmp[j+1];
425     ztmp[j] = ztmp[j+1];
426     }
427     mtype[niatom] = 0;
428     xtmp[niatom] = 0.0F;
429     ytmp[niatom] = 0.0F;
430     ztmp[niatom] = 0.0F;
431     }
432     }
433     if (idummy != 1)
434     break;
435     }
436    
437     /* generate atoms */
438     for (i= 0; i < niatom; i++)
439     newatom = make_atom(0,xtmp[i],ytmp[i],ztmp[i],Elements[mtype[i]-1].symbol);
440    
441    
442     /* need to set atom types and charges */
443     for (i = ibotptr; i < ibotptr + niatom; i++)
444     atom[i+1].charge= mcharge[i];
445    
446     /* need to generate bonds using mopaco */
447     mopaco(ibotptr+1,natom);
448    
449     free_vector(bl, 0,ntemp);
450     free_vector(alpha, 0,ntemp);
451     free_vector(beta, 0,ntemp);
452     free_vector(mcharge, 0,ntemp);
453     free_ivector(na, 0,ntemp);
454     free_ivector(nb, 0,ntemp);
455     free_ivector(nc, 0,ntemp);
456     free_ivector(mtype, 0,ntemp);
457     free_vector(xtmp, 0,ntemp);
458     free_vector(ytmp, 0,ntemp);
459     free_vector(ztmp, 0,ntemp);
460     }
461     /* ======================================================= */
462     void z2xyz(int nat, float *bl, float *ang, float *dihed,int *na, int *nb, int *nc,
463     float *xtmp, float *ytmp, float *ztmp)
464     {
465     // nat = number of atoms
466     // bl, ang, dihed, na, nb, nc = internal coord and connectivity
467     // xtmp, ytmp, ztmp = new cartesian coordinates
468    
469     int i, k, ma, mb, mc;
470     float ccos, cosa, cosd, coskh, cosph, costh, dpr, pi, rbc, sina,
471     sind, sinkh, sinph, sinth, xa, xb, xd, xpa, xpb, xpd, xqa, xqd,
472     xrd, xyb, ya, yb, yd, ypa, ypd, yqd, yza, za, zb, zd, zpd, zqa,
473     zqd;
474    
475     /* nz = end number of atoms
476     * nat = start number of atoms = ibotptr
477     * converts z matrix to cartesians.obtain as many cartesian sets
478     * as entries in z matrix, so dummy atoms are included. */
479    
480     pi = 4*atan( 1.0F );
481     dpr = 180.0F/pi;
482     // atom 1
483     xtmp[0] = 0.0F;
484     ytmp[0] = 0.0F;
485     ztmp[0] = 0.0F;
486     // atom 2
487     xtmp[1] = bl[1];
488     ytmp[1] = 0.0F;
489     ztmp[1] = 0.0F;
490    
491     if (nat > 2)
492     {
493     // atom 3
494     ccos = cos(ang[2]/dpr);
495     if ( na[2] == 1)
496     xtmp[2] = xtmp[0] + bl[2]*ccos;
497     else
498     xtmp[2] = xtmp[1] - bl[2]*ccos;
499    
500     ytmp[2] = bl[2]*sin(ang[2]/dpr);
501     ztmp[2] = 0.0;
502    
503     // atom 4 to end
504     for (i=3; i < nat; i++)
505     {
506     cosa = cos(ang[i]/dpr);
507     mb = nb[i]-1;
508     mc = na[i]-1;
509     xb = xtmp[mb] - xtmp[mc];
510     yb = ytmp[mb] - ytmp[mc];
511     zb = ztmp[mb] - ztmp[mc];
512     rbc = 1.0F/sqrt(xb*xb + yb*yb + zb*zb);
513     if ( fabs(cosa) < (1.0 - 1.0e-10))
514     {
515     // atoms are not collinear
516     ma = nc[i]-1;
517     xa = xtmp[ma]-xtmp[mc];
518     ya = ytmp[ma]-ytmp[mc];
519     za = ztmp[ma]-ztmp[mc];
520     xyb = sqrt(xb*xb + yb*yb);
521     k = -1;
522     if (xyb <= 0.1F)
523     {
524     xpa = za;
525     za = -xa;
526     xa = xpa;
527     xpb = zb;
528     zb = -xb;
529     xb = xpb;
530     xyb = sqrt( xb*xb + yb*yb );
531     k = 1;
532     }
533     /* rotate about the y-axis to make zb vanish */
534     costh = xb/xyb;
535     sinth = yb/xyb;
536     xpa = xa*costh + ya*sinth;
537     ypa = ya*costh - xa*sinth;
538     sinph = zb*rbc;
539     cosph = sqrt( fabs( 1.e00 - sinph*sinph ) );
540     xqa = xpa*cosph + za*sinph;
541     zqa = za*cosph - xpa*sinph;
542     /* rotate about the x-axis to make za=0, and ya positive. */
543     yza = sqrt( ypa*ypa + zqa*zqa );
544     if( yza < 1.e-4 )
545     {
546     /* angle too small to be important */
547     coskh = 1.0F;
548     sinkh = 0.0F;
549     } else
550     {
551     coskh = ypa/yza;
552     sinkh = zqa/yza;
553     }
554     /* coordinates :- a=(xqa,yza,0), b=(rbc,0,0), c=(0,0,0)
555     * none are negative.
556     * the coordinates of i are evaluated in the new frame. */
557     sina = sin( ang[i]/dpr );
558     sind = -sin( dihed[i]/dpr );
559     cosd = cos( dihed[i]/dpr );
560     xd = bl[i]*cosa;
561     yd = bl[i]*sina*cosd;
562     zd = bl[i]*sina*sind;
563     /* transform the coordinates back to the original system. */
564     ypd = yd*coskh - zd*sinkh;
565     zpd = zd*coskh + yd*sinkh;
566     xpd = xd*cosph - zpd*sinph;
567     zqd = zpd*cosph + xd*sinph;
568     xqd = xpd*costh - ypd*sinth;
569     yqd = ypd*costh + xpd*sinth;
570     if( k >= 1 )
571     {
572     xrd = -zqd;
573     zqd = xqd;
574     xqd = xrd;
575     }
576     xtmp[i] = xqd + xtmp[mc];
577     ytmp[i] = yqd + ytmp[mc];
578     ztmp[i] = zqd + ztmp[mc];
579     }else
580     {
581     /* atom mc, mb and i are collinear */
582     rbc *= bl[i]*cosa;
583     xtmp[i] = xtmp[mc] + xb*rbc;
584     ytmp[i] = ytmp[mc] + yb*rbc;
585     ztmp[i] = ztmp[mc] + zb*rbc;
586     }
587     }
588     }
589     }
590     // =======================================
591     void write_mopac()
592     {
593     int **bvec, **dvec;
594     int i, i1, ndx, mckel,
595     iatk1, ihvy, ij, katm1, katm2, katm3 ,
596     iorn, ipass, j, k, jatm, katm, katm4,
597     lps, niatm, nrec, *order, *maxn;
598     float xangl, xdist;
599     double xdihe;
600     FILE *mopfile;
601    
602     /* allocate space for bvec and dvec */
603    
604     bvec = imatrix(0,MAXIAT,1,natom+1);
605     dvec = imatrix(0L,3L,1L,natom+1);
606     order = ivector(0, natom+1);
607     maxn = ivector(0,natom+1);
608    
609     if (selatom[1] != 0 && selatom[2] != 0 && selatom[3] != 0)
610     {
611     katm1 = selatom[1];
612     katm2 = selatom[2];
613     katm3 = selatom[3];
614     mckel = 1;
615     }else
616     {
617     katm1 = 1;
618     katm2 = 2;
619     katm3 = 3;
620     mckel = 0;
621     }
622    
623     lps = 0;
624     /* see if we have lone pairs or metals */
625     for( i = 1; i <= natom; i++ )
626     {
627     if( atom[i].type == 20 )
628     lps = 1;
629     }
630    
631     if( lps == 1 )
632     hdel(1);
633     type();
634     mopfile = fopen_path(Savebox.path,Savebox.fname,"w");
635    
636    
637     /* MOPAC */
638     fprintf(mopfile,"%s\n",mopdata.keywords);
639     fprintf(mopfile,"%s\n",mopdata.title);
640     fprintf(mopfile,"%s\n",mopdata.title1);
641    
642     for( ij = 1; ij <= natom; ij++ )
643     {
644     order[ij] = 0;
645     maxn[ij] = 0;
646     ndx = 0;
647     for (i=0; i < MAXIAT; i++)
648     {
649     if (atom[ij].iat[i] != 0)
650     {
651     if (atom[ij].iat[i] > ndx && atom[ij].iat[i] < ij)
652     ndx = atom[ij].iat[i];
653     }
654     }
655     maxn[ij] = ndx;
656     if (ndx == 0 && ij == 2)
657     maxn[ij] = 1;
658     else if (ndx == 0 && ij == 3)
659     maxn[ij] = 2;
660     if (ndx == 0 && ij > 3) // check on separate structures
661     {
662     maxn[ij] = maxn[ij-1];
663     }
664     }
665     for( i = 1; i <= natom; i++ )
666     {
667     for( k = 0; k < 4; k++ )
668     {
669     dvec[k][i] = 0.;
670     }
671     for( k = 0; k < MAXIAT; k++ )
672     {
673     bvec[k][i] = atom[i].iat[k];
674     }
675     }
676    
677    
678    
679     i1 = katm1;
680     order[katm1] = 1;
681     dvec[0][i1] = i1;
682     dvec[1][i1] = 0;
683     dvec[2][i1] = 0;
684     dvec[3][i1] = 0;
685    
686     fprintf( mopfile, "%2s%12.6f%3d%12.6f%3d%12.6f%3d %3d %3d %3d \n", Elements[atom[i1].atomnum-1].symbol,
687     0.0, 0, 0.0, 0, 0.0, 0, 0, 0, 0);
688    
689     i1 = katm2;
690     order[katm2] = 2;
691     dvec[0][i1] = katm2;
692     dvec[1][i1] = katm1;
693     dvec[2][i1] = 0;
694     dvec[3][i1] = 0;
695     distjm( dvec[0][i1], dvec[1][i1], &xdist );
696     fprintf( mopfile, "%2s%12.6f%3d%12.6f%3d%12.6f%3d %3d %3d %3d \n",
697     Elements[atom[i1].atomnum-1].symbol, xdist, mopdata.ibopt,0.0, 0, 0.0, 0, order[dvec[1][i1]], 0, 0 );
698    
699     i1 = katm3;
700     order[katm3] = 3;
701     dvec[0][i1] = katm3;
702     dvec[1][i1] = katm2;
703     dvec[2][i1] = katm1;
704     dvec[3][i1] = 0;
705     for (ij=0; ij < MAXIAT; ij++) // check to see if 3 is bonded to 1
706     {
707     if (atom[katm3].iat[ij] != 0 && atom[katm3].iat[ij] == katm1)
708     {
709     dvec[1][i1] = katm1;
710     dvec[2][i1] = katm2;
711     break;
712     }
713     }
714     distjm( dvec[0][i1], dvec[1][i1], &xdist );
715     angljm( dvec[0][i1], dvec[1][i1], dvec[2][i1], &xangl );
716     fprintf( mopfile, "%2s%12.6f%3d%12.6f%3d%12.6f%3d %3d %3d %3d \n",
717     Elements[atom[i1].atomnum-1].symbol, xdist, mopdata.ibopt, xangl, mopdata.iaopt, 0.0, 0, order[dvec[1][i1]],
718     order[dvec[2][i1]], 0 );
719    
720     nrec = 3;
721     niatm = katm3;
722     ipass = 0;
723     ihvy = 0;
724    
725     if (mckel == 0)
726     {
727     for (ij = 4; ij <= natom; ij++)
728     {
729     i1 = ij;
730     order[ij] = ij;
731     dvec[0][i1] = ij;
732     dvec[1][i1] = maxn[ij];
733     if (maxn[ij] == 1)
734     {
735     dvec[2][i1] = 2;
736     dvec[3][i1] = 3;
737     } else if (maxn[ij] == 2)
738     {
739     dvec[3][i1] = 3;
740     }
741     if (dvec[2][i1] == 0)
742     dvec[2][i1] = maxn[maxn[ij]];
743     if (dvec[3][i1] == 0)
744     dvec[3][i1] = maxn[dvec[2][i1]];
745     if (dvec[3][i1] == 0)
746     {
747     jatm = dvec[2][i1];
748     for (i=0; i < MAXIAT; i++)
749     {
750     if (atom[jatm].iat[i] != 0 && atom[jatm].iat[i] != dvec[1][i1] && atom[jatm].iat[i] != dvec[0][i1])
751     {
752     dvec[3][i1] = atom[jatm].iat[i];
753     break;
754     }
755     }
756     }
757    
758     if (dvec[3][i1]*dvec[2][i1]*dvec[1][i1] == 0)
759     message_alert("bad vector in mopac setup","Error");
760    
761     distjm( dvec[0][i1], dvec[1][i1], &xdist );
762     angljm( dvec[0][i1], dvec[1][i1], dvec[2][i1], &xangl );
763     xdihe = dihdrl( dvec[0][i1], dvec[1][i1], dvec[2][i1], dvec[3][i1] );
764     fprintf( mopfile, "%2s%12.6f%3d%12.6f%3d%12.6f%3d %3d %3d %3d \n",
765     Elements[atom[i1].atomnum-1].symbol, xdist, mopdata.ibopt, xangl, mopdata.iaopt,
766     xdihe, mopdata.idopt, order[dvec[1][i1]],
767     order[dvec[2][i1]], order[dvec[3][i1]] );
768     }
769     } else
770     {
771    
772     while( TRUE )
773     {
774     if( niatm > natom )
775     niatm = 1;
776     if( order[niatm] && atom[niatm].iat[1] > 0 )
777     {
778     for( j = 0; j < MAXIAT; j++ )
779     {
780     katm = bvec[j][niatm];
781     if( katm != 0 )
782     {
783     if( (atom[katm].iat[1] != 0 && ihvy == 0) || (atom[katm].iat[1] == 0 && ihvy == 1) )
784     {
785     if( katm && !order[katm] )
786     {
787     nrec += 1;
788     if( nrec == 4 )
789     katm4 = katm;
790     if( nrec > natom )
791     goto L_200;
792     order[katm] = nrec;
793     dvec[0][katm] = katm;
794     iorn = order[niatm];
795     if( iorn == 1 )
796     {
797     dvec[1][katm] = katm1;
798     dvec[2][katm] = katm2;
799     dvec[3][katm] = katm3;
800     }else if( iorn == 2 )
801     {
802     dvec[1][katm] = katm2;
803     if( !atom[katm3].iat[1] )
804     {
805     dvec[2][katm] = katm1;
806     for( ij = 0; ij < MAXIAT; ij++ )
807     {
808     iatk1 = atom[katm1].iat[ij];
809     if( iatk1 )
810     {
811     if( iatk1 != katm2 && order[iatk1] > 0 )
812     goto L_98766;
813     }
814     }
815     dvec[3][katm] = katm3;
816     goto L_110;
817     L_98766:
818     dvec[3][katm] = iatk1;
819     }else
820     {
821     dvec[2][katm] = katm3;
822     dvec[3][katm] = katm1;
823     }
824     }else
825     {
826     dvec[1][katm] = dvec[0][niatm];
827     dvec[2][katm] = dvec[1][niatm];
828     dvec[3][katm] = dvec[2][niatm];
829     }
830     L_110:
831     distjm( dvec[0][katm], dvec[1][katm], &xdist );
832     angljm( dvec[0][katm], dvec[1][katm], dvec[2][katm], &xangl );
833     xdihe = dihdrl( dvec[0][katm], dvec[1][katm], dvec[2][katm], dvec[3][katm] );
834     fprintf( mopfile, "%2s%12.6f%3d%12.6f%3d%12.6f%3d %3d %3d %3d \n",
835     Elements[atom[katm].atomnum-1].symbol, xdist, mopdata.ibopt, xangl, mopdata.iaopt,
836     xdihe, mopdata.idopt, order[dvec[1][katm]],
837     order[dvec[2][katm]], order[dvec[3][katm]] );
838     }
839     }
840     }
841     }
842     }
843    
844     if( nrec == natom )
845     goto L_200;
846     if( niatm == natom )
847     ipass += 1;
848     if( ipass <= 5 )
849     {
850     niatm += 1;
851     }
852     else if( nrec < natom && !ihvy )
853     {
854     ipass = 0;
855     ihvy = 1;
856     }
857     else
858     {
859     break;
860     }
861     }
862     goto L_98765;
863     }
864     L_200:
865     ;
866     fprintf(mopfile,"\n");
867    
868     L_98765:
869     fclose(mopfile);
870     free_imatrix(bvec,0,MAXIAT,1,natom+1);
871     free_imatrix(dvec,0,3,1,natom+1);
872     free_ivector(order,0, natom+1);
873     free_ivector(maxn,0, natom+1);
874    
875     if( lps == 1 )
876     hadd();
877     generate_bonds();
878     return;
879     }
880    
881    
882     /*---------------------------------------------------------------- */
883     void distjm(int id1,int id2, float *xdist)
884     {
885     float xdf, xx, ydf, zdf;
886    
887    
888     xdf = atom[id1].x - atom[id2].x;
889     ydf = atom[id1].y - atom[id2].y;
890     zdf = atom[id1].z - atom[id2].z;
891    
892     xx = xdf*xdf + ydf*ydf + zdf*zdf;
893     if( xx < 0.0 )
894     xx = 0.0;
895     *xdist = sqrt( xx );
896     return;
897     }
898     /*--------------------------------------------------------------- */
899     void angljm(int ia1, int ia2, int ia3, float *xangl)
900     {
901     int i, ia[3], j;
902     float a[3], aa123, b[3], bb123, fxcos, p[3][3], xa, xb, xcos,
903     xxcos;
904    
905     ia[0] = ia1;
906     ia[1] = ia2;
907     ia[2] = ia3;
908    
909     for( j = 0; j < 3; j++ )
910     {
911     p[0][j] = atom[ia[j]].x;
912     p[1][j] = atom[ia[j]].y;
913     p[2][j] = atom[ia[j]].z;
914     }
915    
916     for( i = 0; i < 3; i++ )
917     {
918     a[i] = p[i][0] - p[i][1];
919     b[i] = p[i][2] - p[i][1];
920     }
921    
922     aa123 = a[0]*a[0] + a[1]*a[1] + a[2]*a[2];
923     xa = sqrt( aa123 );
924     bb123 = b[0]*b[0] + b[1]*b[1] + b[2]*b[2];
925     xb = sqrt( bb123 );
926    
927     xxcos = a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
928     fxcos = xxcos/xa;
929     xcos = fxcos/xb;
930     if( xcos > 1.0 )
931     xcos = 1.0;
932     if( xcos < -1.0 )
933     xcos = -1.0;
934     *xangl = acos( xcos );
935     *xangl *= 57.29578;
936    
937     return;
938     }