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

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