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, 4 months ago) by tjod
File size: 31443 byte(s)
Log Message:
test

Line File contents
1 #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 }