ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/ktorsion.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 11 months ago) by tjod
File size: 17281 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 "torsions.h"
6     #include "rings.h"
7     #include "field.h"
8     #include "atom_k.h"
9     #include "fix.h"
10    
11     int isbond(int,int);
12     int is_delocalbond(int,int);
13     void numeral(int,char *,int);
14     void angle(int,int,int,float *);
15     void four(char *,char *,char *,char *, char *);
16     int is_ring42(int,int);
17     int is_ring54(int, int, int, int);
18     void message_alert(char *, char *);
19     void transomeg(float *, float *, float *, int *, int *, char *, char *, char *, char *,
20     char *, char *, char *, char *, char *, char *, char *);
21    
22     EXTERN struct t_allene {
23     int nallene, ntor[10];
24     } allene;
25    
26    
27     EXTERN struct t_torkn1 {
28     int use_tor4, use_tor5;
29     int ntor, ntor4, ntor5, ntordel, torindex[MAXTORDEL];
30     char kv[MAXTORCONST][13], kv4[MAXTOR4CONST][13], kv5[MAXTOR5CONST][13];
31     char kvdel[MAXTORDEL][13];
32     float tv1[MAXTORCONST], tv2[MAXTORCONST], tv3[MAXTORCONST];
33     float tv4[MAXTORCONST], tv5[MAXTORCONST], tv6[MAXTORCONST];
34     int phase1[MAXTORCONST], phase2[MAXTORCONST], phase3[MAXTORCONST];
35     int phase4[MAXTORCONST], phase5[MAXTORCONST], phase6[MAXTORCONST];
36     float tv41[MAXTOR4CONST], tv42[MAXTOR4CONST], tv43[MAXTOR4CONST];
37     int phase41[MAXTOR4CONST], phase42[MAXTOR4CONST], phase43[MAXTOR4CONST];
38     float tv51[MAXTOR5CONST], tv52[MAXTOR5CONST], tv53[MAXTOR5CONST];
39     int phase51[MAXTOR5CONST], phase52[MAXTOR5CONST], phase53[MAXTOR5CONST];
40     float tvdel1[MAXTORDEL], tvdel2[MAXTORDEL], tvdel3[MAXTORDEL];
41     int phasedel1[MAXTORDEL], phasedel2[MAXTORDEL], phasedel3[MAXTORDEL];
42     } torkn1;
43    
44     EXTERN struct t_torknp {
45     int npitor;
46     char kp[MAXPITORCONST][13];
47     int ph1[MAXPITORCONST], ph2[MAXPITORCONST], ph3[MAXPITORCONST];
48     float tv1[MAXPITORCONST], tv2[MAXPITORCONST], tv3[MAXPITORCONST];
49     } torknp;
50    
51    
52     EXTERN struct t_ts_bondorder {
53     float fbnd[15];
54     } ts_bondorder;
55     EXTERN struct t_high_coord {
56     int ncoord, i13[400][3];
57     } high_coord;
58    
59     EXTERN int Missing_constants;
60     //EXTERN FILE *errfile;
61    
62     void ktorsion()
63     {
64     int i, j, ierr, itor, ii, ki, k;
65     int ia, ib, ic, id;
66     int ita, itb, itc, itd;
67     int cl_a, cl_b, cl_c, cl_d;
68     int use_ring4, use_ring5;
69     long int mask;
70     char izero[4];
71     char pa[4],pb[4],pc[4],pd[4],pt[13], kv1[13];
72     char k1[13], k2[13], k3[13], k4[13], k5[13], k6[13], k7[13], k8[13], k9[13], k10[13],
73     k11[13], k12[13], k13[13], k14[13], k15[13];
74     float v1, v2, v3, rangle;
75    
76     mask = (1L << 0);
77     itor = 0;
78     use_ring4 = FALSE;
79     if (rings.nring4 > 0 && torkn1.ntor4 > 0)
80     use_ring4 = TRUE;
81    
82     use_ring5 = FALSE;
83     if (rings.nring5 > 0 && torkn1.ntor5 > 0)
84     use_ring5 = TRUE;
85    
86     for (i=0; i < torsions.ntor; i++)
87     {
88     ia = torsions.i14[i][0];
89     ib = torsions.i14[i][1];
90     ic = torsions.i14[i][2];
91     id = torsions.i14[i][3];
92    
93     ita = atom[ia].type;
94     itb = atom[ib].type;
95     itc = atom[ic].type;
96     itd = atom[id].type;
97    
98     cl_a = atom[ia].tclass;
99     cl_b = atom[ib].tclass;
100     cl_c = atom[ic].tclass;
101     cl_d = atom[id].tclass;
102    
103     if (field.type == MMX)
104     {
105     if (ita == 40)
106     ita = 2;
107     if (itb == 40)
108     itb = 2;
109     if (itc == 40)
110     itc = 2;
111     if (itd == 40)
112     itd = 2;
113    
114     if (ita == 56 && itb == 56 && itc == 56 && itd == 56)
115     {
116     ita = 1; itb = 1; itc = 1; itd = 1;
117     }
118     }
119     numeral(ita,pa,3);
120     numeral(itb,pb,3);
121     numeral(itc,pc,3);
122     numeral(itd,pd,3);
123    
124     if (itb < itc )
125     four(pa,pb, pc, pd, pt);
126     else if (itb > itc)
127     four(pd,pc, pb, pa, pt);
128     else if (ita < itd)
129     four(pa,pb, pc, pd, pt);
130     else
131     four(pd,pc, pb, pa, pt);
132     strcpy(izero," 0");
133    
134     if (field.type == MMFF94)
135     {
136    
137     numeral(atom_k.tclass1[ita],pa,3);
138     numeral(atom_k.tclass[itb],pb,3);
139     numeral(atom_k.tclass[itc],pc,3);
140     numeral(atom_k.tclass1[itd],pd,3);
141    
142    
143     if (atom_k.tclass[itb] < atom_k.tclass[itc] )
144     {
145     four(pa, pb, pc, izero,k2);
146     four(izero ,pb, pc, pd,k3);
147     } else if (atom_k.tclass[itb] > atom_k.tclass[itc])
148     {
149     four(izero,pc, pb, pa, k2);
150     four(pd,pc, pb, izero, k3);
151     } else if (atom_k.tclass1[ita] < atom_k.tclass1[itd])
152     {
153     four(pa,pb, pc, izero, k2);
154     four(izero,pb, pc, pd, k3);
155     } else
156     {
157     four(izero,pc, pb, pa, k2);
158     four(pd,pc, pb, izero, k3);
159     }
160     } else
161     {
162     four( izero, pc, pb, pa, k1 );
163     four( pd, pc, pb, izero, k2 );
164     four( izero, pb, pc, pd, k3 );
165     }
166    
167     numeral(ita,pa,3);
168     numeral(itb,pb,3);
169     numeral(itc,pc,3);
170     numeral(itd,pd,3);
171    
172     four( pa, pb, pc, izero, k4 );
173     four( izero, pc, pb, izero, k5 );
174     four( izero, pb, pc, izero, k6 );
175     four( pd, pc, izero, izero, k7 );
176     four( izero, izero, pc,pd, k8 );
177    
178     four( izero, izero, pb, pa, k9 );
179     four( pa, pb, izero, izero, k10 );
180     four( izero, izero, pc, izero, k11 );
181     four( izero, pc, izero, izero, k12 );
182     four( izero, izero, pb, izero, k13 );
183     four( izero, pb, izero, izero, k14 );
184    
185     numeral(cl_a,pa,3);
186     numeral(cl_b,pb,3);
187     numeral(cl_c,pc,3);
188     numeral(cl_d,pd,3);
189     if (itb < itc )
190     four(pa,pb, pc, pd, k15);
191     else if (itb > itc)
192     four(pd,pc, pb, pa, k15);
193     else if (ita < itd)
194     four(pa,pb, pc, pd, k15);
195     else
196     four(pd,pc, pb, pa, k15);
197    
198     ierr = FALSE;
199     // check for linear angles in high coordinate atoms
200     if (high_coord.ncoord > 0)
201     {
202     for (k = 0; k < high_coord.ncoord; k++)
203     {
204     if (high_coord.i13[k][1] == ib) // high coordinate atom in center
205     {
206     if (high_coord.i13[k][0] == ia && high_coord.i13[k][2] == ic)
207     {
208     angle(ia,ib,ic,&rangle);
209     if (rangle > 175 || rangle < -175)
210     {
211     torsions.v1[i] = 0.0;
212     torsions.v2[i] = 0.0;
213     torsions.v3[i] = 0.0;
214     ierr = TRUE;
215     goto L_10;
216     }
217     } else if (high_coord.i13[k][2] == ia && high_coord.i13[k][0] == ic)
218     {
219     angle(ia,ib,ic,&rangle);
220     if (rangle > 175 || rangle < -175)
221     {
222     torsions.v1[i] = 0.0;
223     torsions.v2[i] = 0.0;
224     torsions.v3[i] = 0.0;
225     ierr = TRUE;
226     goto L_10;
227     }
228     }
229     } else if ( high_coord.i13[k][1] == ic)
230     {
231     if (high_coord.i13[k][0] == ib && high_coord.i13[k][2] == id)
232     {
233     angle(ib,ic,id,&rangle);
234     if (rangle > 175 || rangle < -175)
235     {
236     torsions.v1[i] = 0.0;
237     torsions.v2[i] = 0.0;
238     torsions.v3[i] = 0.0;
239     ierr = TRUE;
240     goto L_10;
241     }
242     } else if (high_coord.i13[k][2] == ib&& high_coord.i13[k][0] == id)
243     {
244     angle(ib,ic,id,&rangle);
245     if (rangle > 175 || rangle < -175)
246     {
247     torsions.v1[i] = 0.0;
248     torsions.v2[i] = 0.0;
249     torsions.v3[i] = 0.0;
250     ierr = TRUE;
251     goto L_10;
252     }
253     }
254     }
255     }
256     }
257    
258     // check for four membered rings
259     if ( isbond(ia,id) )
260     {
261     for(j=0; j < torkn1.ntor4; j++)
262     {
263     strcpy(kv1,torkn1.kv4[j]);
264     if (strcmp(pt,kv1) == 0 || strcmp(k1,kv1) == 0 || strcmp(k2,kv1) == 0 || strcmp(k3,kv1) == 0 ||
265     strcmp(k4,kv1) == 0 || strcmp(k5,kv1) == 0 || strcmp(k6,kv1) == 0 || strcmp(k7,kv1) == 0 ||
266     strcmp(k8,kv1) == 0 || strcmp(k9,kv1) == 0 || strcmp(k10,kv1) == 0 || strcmp(k11,kv1) == 0 ||
267     strcmp(k12,kv1) == 0 || strcmp(k13,kv1) == 0 || strcmp(k14,kv1) == 0 )
268     {
269     torsions.v1[i] = torkn1.tv41[j];
270     torsions.ph1[i] = torkn1.phase41[j];
271     torsions.v2[i] = torkn1.tv42[j];
272     torsions.ph2[i] = torkn1.phase42[j];
273     torsions.v3[i] = torkn1.tv43[j];
274     torsions.ph3[i] = torkn1.phase43[j];
275     ierr = TRUE;
276     goto L_10;
277     }
278     }
279     }
280     // delocalized torsions
281     if ( is_delocalbond(ib,ic) )
282     {
283     for(j=0; j < torkn1.ntordel; j++)
284     {
285     strcpy(kv1,torkn1.kvdel[j]);
286     if (strcmp(pt,kv1) == 0 || strcmp(k1,kv1) == 0 || strcmp(k2,kv1) == 0 || strcmp(k3,kv1) == 0 ||
287     strcmp(k4,kv1) == 0 || strcmp(k5,kv1) == 0 || strcmp(k6,kv1) == 0 || strcmp(k7,kv1) == 0 ||
288     strcmp(k8,kv1) == 0 || strcmp(k9,kv1) == 0 || strcmp(k10,kv1) == 0 || strcmp(k11,kv1) == 0 ||
289     strcmp(k12,kv1) == 0 || strcmp(k13,kv1) == 0 || strcmp(k14,kv1) == 0 )
290     {
291     torsions.v1[i] = torkn1.tvdel1[j];
292     torsions.ph1[i] = torkn1.phasedel1[j];
293     torsions.v2[i] = torkn1.tvdel2[j];
294     torsions.ph2[i] = torkn1.phasedel2[j];
295     torsions.v3[i] = torkn1.tvdel3[j];
296     torsions.ph3[i] = torkn1.phasedel3[j];
297     ierr = TRUE;
298     goto L_10;
299     }
300     }
301     }
302     // check for five membered rings
303     if (is_ring54(ia,ib,ic,id) )
304     {
305     for(j=0; j < torkn1.ntor5; j++)
306     {
307     strcpy(kv1,torkn1.kv5[j]);
308     if (strcmp(pt,kv1) == 0 || strcmp(k1,kv1) == 0 || strcmp(k2,kv1) == 0 || strcmp(k3,kv1) == 0 ||
309     strcmp(k4,kv1) == 0 || strcmp(k5,kv1) == 0 || strcmp(k6,kv1) == 0 || strcmp(k7,kv1) == 0 ||
310     strcmp(k8,kv1) == 0 || strcmp(k9,kv1) == 0 || strcmp(k10,kv1) == 0 || strcmp(k11,kv1) == 0 ||
311     strcmp(k12,kv1) == 0 || strcmp(k13,kv1) == 0 || strcmp(k14,kv1) == 0 )
312     {
313     torsions.v1[i] = torkn1.tv51[j];
314     torsions.ph1[i] = torkn1.phase51[j];
315     torsions.v2[i] = torkn1.tv52[j];
316     torsions.ph2[i] = torkn1.phase52[j];
317     torsions.v3[i] = torkn1.tv53[j];
318     torsions.ph3[i] = torkn1.phase53[j];
319     ierr = TRUE;
320     goto L_10;
321     }
322     }
323     }
324     // check regular parameters
325     if (field.type == MMFF94)
326     {
327     for(j=0; j < torkn1.ntor; j++) // specific constants
328     {
329     strcpy(kv1,torkn1.kv[j]);
330     if (strcmp(pt,kv1) == 0)
331     {
332     torsions.v1[i] = torkn1.tv1[j];
333     torsions.ph1[i] = torkn1.phase1[j];
334     torsions.v2[i] = torkn1.tv2[j];
335     torsions.ph2[i] = torkn1.phase2[j];
336     torsions.v3[i] = torkn1.tv3[j];
337     torsions.ph3[i] = torkn1.phase3[j];
338     torsions.v4[i] = torkn1.tv4[j];
339     torsions.ph4[i] = torkn1.phase4[j];
340     torsions.v5[i] = torkn1.tv5[j];
341     torsions.ph5[i] = torkn1.phase5[j];
342     torsions.v6[i] = torkn1.tv6[j];
343     torsions.ph6[i] = torkn1.phase6[j];
344     goto L_10;
345     }
346     }
347     for(j=0; j < torkn1.ntor; j++) // class 2
348     {
349     strcpy(kv1,torkn1.kv[j]);
350     if (strcmp(k2,kv1) == 0 || strcmp(k3,kv1) == 0 )
351     {
352     torsions.v1[i] = torkn1.tv1[j];
353     torsions.ph1[i] = torkn1.phase1[j];
354     torsions.v2[i] = torkn1.tv2[j];
355     torsions.ph2[i] = torkn1.phase2[j];
356     torsions.v3[i] = torkn1.tv3[j];
357     torsions.ph3[i] = torkn1.phase3[j];
358     torsions.v4[i] = torkn1.tv4[j];
359     torsions.ph4[i] = torkn1.phase4[j];
360     torsions.v5[i] = torkn1.tv5[j];
361     torsions.ph5[i] = torkn1.phase5[j];
362     torsions.v6[i] = torkn1.tv6[j];
363     torsions.ph6[i] = torkn1.phase6[j];
364     goto L_10;
365     }
366     }
367     for(j=0; j < torkn1.ntor; j++) // wild cards
368     {
369     strcpy(kv1,torkn1.kv[j]);
370     if (strcmp(k5,kv1) == 0 || strcmp(k6,kv1) == 0 || strcmp(k7,kv1) == 0 || strcmp(k8,kv1) == 0 )
371     {
372     torsions.v1[i] = torkn1.tv1[j];
373     torsions.ph1[i] = torkn1.phase1[j];
374     torsions.v2[i] = torkn1.tv2[j];
375     torsions.ph2[i] = torkn1.phase2[j];
376     torsions.v3[i] = torkn1.tv3[j];
377     torsions.ph3[i] = torkn1.phase3[j];
378     torsions.v4[i] = torkn1.tv4[j];
379     torsions.ph4[i] = torkn1.phase4[j];
380     torsions.v5[i] = torkn1.tv5[j];
381     torsions.ph5[i] = torkn1.phase5[j];
382     torsions.v6[i] = torkn1.tv6[j];
383     torsions.ph6[i] = torkn1.phase6[j];
384     goto L_10;
385     }
386     }
387     }
388     // check regular specific parameters
389     for(j=0; j < torkn1.ntor; j++)
390     {
391     strcpy(kv1,torkn1.kv[j]);
392     if (strcmp(pt,kv1) == 0)
393     {
394     torsions.v1[i] = torkn1.tv1[j];
395     torsions.ph1[i] = torkn1.phase1[j];
396     torsions.v2[i] = torkn1.tv2[j];
397     torsions.ph2[i] = torkn1.phase2[j];
398     torsions.v3[i] = torkn1.tv3[j];
399     torsions.ph3[i] = torkn1.phase3[j];
400     torsions.v4[i] = torkn1.tv4[j];
401     torsions.ph4[i] = torkn1.phase4[j];
402     torsions.v5[i] = torkn1.tv5[j];
403     torsions.ph5[i] = torkn1.phase5[j];
404     torsions.v6[i] = torkn1.tv6[j];
405     torsions.ph6[i] = torkn1.phase6[j];
406     goto L_10;
407     break;
408     }
409     }
410     // check regular generalized parameters
411     for(j=0; j < torkn1.ntor; j++)
412     {
413     strcpy(kv1,torkn1.kv[j]);
414     if (strcmp(pt,kv1) == 0 || strcmp(k1,kv1) == 0 || strcmp(k2,kv1) == 0 || strcmp(k3,kv1) == 0 ||
415     strcmp(k4,kv1) == 0 || strcmp(k5,kv1) == 0 || strcmp(k6,kv1) == 0 || strcmp(k7,kv1) == 0 ||
416     strcmp(k8,kv1) == 0 || strcmp(k9,kv1) == 0 || strcmp(k10,kv1) == 0 || strcmp(k11,kv1) == 0 ||
417     strcmp(k12,kv1) == 0 || strcmp(k13,kv1) == 0 || strcmp(k14,kv1) == 0 || strcmp(k15,kv1) == 0)
418     {
419     torsions.v1[i] = torkn1.tv1[j];
420     torsions.ph1[i] = torkn1.phase1[j];
421     torsions.v2[i] = torkn1.tv2[j];
422     torsions.ph2[i] = torkn1.phase2[j];
423     torsions.v3[i] = torkn1.tv3[j];
424     torsions.ph3[i] = torkn1.phase3[j];
425     torsions.v4[i] = torkn1.tv4[j];
426     torsions.ph4[i] = torkn1.phase4[j];
427     torsions.v5[i] = torkn1.tv5[j];
428     torsions.ph5[i] = torkn1.phase5[j];
429     torsions.v6[i] = torkn1.tv6[j];
430     torsions.ph6[i] = torkn1.phase6[j];
431     goto L_10;
432     break;
433     }
434     }
435     // missing parameters
436     if (ierr == FALSE)
437     {
438     // Missing_constants = TRUE;
439     }
440     L_10:
441     continue;
442     }
443    
444     if (allene.nallene > 0)
445     {
446     for (i = 0; i < allene.nallene; i++)
447     {
448     torsions.v1[allene.ntor[i]] = 0;
449     torsions.v2[allene.ntor[i]] = -11.5;
450     torsions.v3[allene.ntor[i]] = 0;
451     }
452     }
453     }
454     /* -------------------------------------------- */
455     void four(char *pa, char *pb, char *pc, char *pd, char *pt)
456     {
457     strcpy(pt,pa);
458     strcat(pt,pb);
459     strcat(pt,pc);
460     strcat(pt,pd);
461     }