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

Line User Rev File contents
1 tjod 3 #define EXTERN extern
2     #include "pcwin.h"
3     #include "pcmod.h"
4    
5     #include "angles.h"
6     #include "bonds_ff.h"
7     #include "field.h"
8    
9     int find_bond(int, int);
10     void numeral(int, char *, int);
11    
12     EXTERN struct t_crossterm_k {
13     int nangang, nstrbnd, nstrtor;
14     int ang_ang[MAXAA], stbnindex[MAXSTBN] ;
15     char str_tor[MAXSTRTOR][7], stbn[MAXSTBN][10];
16     float aacon[MAXAA][3], stbncon[MAXSTBN][3], str_torcon[MAXSTRTOR];
17     } crossterm_k;
18    
19     struct t_strbnd {
20     int nstrbnd, isb[MAXANG][3];
21     float ksb1[MAXANG], ksb2[MAXANG];
22     } strbnd;
23     struct stbndat {
24     int irow, krow, lrow;
25     float con1, con2;
26     } stbn[] =
27     {
28     0, 1, 0, 0.15, 0.15,
29     0, 1, 1, 0.10, 0.30,
30     0, 1, 2, 0.05, 0.35,
31     0, 1, 3, 0.05, 0.35,
32     0, 1, 4, 0.05, 0.35,
33     0, 2, 0, 0.00, 0.00,
34     0, 2, 1, 0.00, 0.15,
35     0, 2, 2, 0.00, 0.15,
36     0, 2, 3, 0.00, 0.15,
37     0, 2, 4, 0.00, 0.15,
38     1, 1, 1, 0.30, 0.30,
39     1, 1, 2, 0.30, 0.50,
40     1, 1, 3, 0.30, 0.50,
41     1, 1, 4, 0.30, 0.50,
42     2, 1, 2, 0.50, 0.50,
43     2, 1, 3, 0.50, 0.50,
44     2, 1, 4, 0.50, 0.50,
45     3, 1, 3, 0.50, 0.50,
46     3, 1, 4, 0.50, 0.50,
47     4, 1, 4, 0.50, 0.50,
48     1, 2, 1, 0.30, 0.30,
49     1, 2, 2, 0.25, 0.25,
50     1, 2, 3, 0.25, 0.25,
51     1, 2, 4, 0.25, 0.25,
52     2, 2, 2, 0.25, 0.25,
53     2, 2, 3, 0.25, 0.25,
54     2, 2, 4, 0.25, 0.25,
55     3, 2, 3, 0.25, 0.25,
56     3, 2, 4, 0.25, 0.25,
57     4, 2, 4, 0.25, 0.25 };
58    
59     EXTERN int Missing_constants;
60     //EXTERN FILE *errfile;
61    
62     void kstrbnd()
63     {
64     int i, j, iz, nb1, nb2, irow, krow, lrow;
65     int ia, ib, ic, itb, ita, itc, ierr;
66     int k, nstbn, index;
67     char pa[4],pb[4],pc[4],pt[10],pt1[10];
68    
69     strbnd.nstrbnd = 0;
70     nb1 = -1;
71     nb2 = -1;
72    
73     for (i=0; i < angles.nang; i++)
74     {
75     ierr = FALSE;
76     ia = angles.i13[i][0];
77     ib = angles.i13[i][1];
78     ic = angles.i13[i][2];
79     ita = atom[ia].type;
80     itb = atom[ib].type;
81     itc = atom[ic].type;
82     if (field.type == MMX && (atom[ia].type >= 300 || atom[ib].type >= 300 || atom[ic].type >= 300))
83     goto L_10;
84     if (field.type == MMX && (atom[ia].type >= 20 || atom[ic].type >= 20))
85     goto L_10;
86     if (field.type == MM3 && (atom[ia].type >= 300 || atom[ib].type >= 300 || atom[ic].type >= 300) )
87     goto L_10;
88     nstbn = -1;
89     numeral(ita,pa,3);
90     numeral(itb,pb,3);
91     numeral(itc,pc,3);
92     if (ita <= itc)
93     {
94     strcpy(pt,pa); strcat(pt,pb); strcat(pt,pc);
95     }else
96     {
97     strcpy(pt,pc); strcat(pt,pb); strcat(pt,pa);
98     }
99     itb = atom[ib].tclass;
100     numeral(itb,pb,3);
101     strcpy(pt1," 0"); strcat(pt1,pb); strcat(pt1," 0");
102    
103     for (j=0; j < crossterm_k.nstrbnd; j++)
104     {
105     if (strcmp(crossterm_k.stbn[j],pt) == 0 || strcmp(crossterm_k.stbn[j],pt1) == 0 )
106     {
107     nstbn = j;
108     ierr = TRUE;
109     break;
110     }
111     }
112     for (j=0; j < bonds_ff.nbnd; j++)
113     {
114     iz = find_bond(ia, ib);
115     if (iz != -1)
116     {
117     nb1 = iz;
118     break;
119     }
120     }
121     for (j=0; j < bonds_ff.nbnd; j++)
122     {
123     iz = find_bond(ib, ic);
124     if (iz != -1)
125     {
126     nb2 = iz;
127     break;
128     }
129     }
130     if (angles.index[i] == 0)
131     index = 0;
132     else if (angles.index[i] == 1)
133     {
134     if (bonds_ff.index[nb1] == 1)
135     index = 1;
136     else
137     index = 2;
138     } else if (angles.index[i] == 2)
139     index = 3;
140     else if (angles.index[i] == 3)
141     index = 5;
142     else if (angles.index[i] == 4)
143     index = 4;
144     else if (angles.index[i] == 5)
145     {
146     if (bonds_ff.index[nb1] == 1)
147     index = 6;
148     else
149     index = 7;
150     } else if (angles.index[i] == 6)
151     index = 8;
152     else if (angles.index[i] == 7)
153     {
154     if (bonds_ff.index[nb1] == 1)
155     index = 9;
156     else
157     index = 10;
158     } else if (angles.index[i] == 8)
159     index = 11;
160    
161     if (crossterm_k.stbnindex[nstbn] != index)
162     nstbn = -1;
163    
164     if (field.type == MMFF94)
165     {
166     ierr = TRUE;
167     if (angles.anat[i] < 175.0)
168     {
169     if (nstbn == -1) // constant not found assign constants base on periodic table
170     {
171     if (atom[ia].atomnum <= 1)
172     irow = 0;
173     else if (atom[ia].atomnum <= 10)
174     irow = 1;
175     else if (atom[ia].atomnum <= 18)
176     irow = 2;
177     else if (atom[ia].atomnum <= 36)
178     irow = 3;
179     else if (atom[ia].atomnum <= 54)
180     irow = 4;
181    
182     if (atom[ib].atomnum <= 1)
183     krow = 0;
184     else if (atom[ib].atomnum <= 10)
185     krow = 1;
186     else if (atom[ib].atomnum <= 18)
187     krow = 2;
188     else if (atom[ib].atomnum <= 36)
189     krow = 3;
190     else if (atom[ib].atomnum <= 54)
191     krow = 4;
192    
193     if (atom[ic].atomnum <= 1)
194     lrow = 0;
195     else if (atom[ic].atomnum <= 10)
196     lrow = 1;
197     else if (atom[ic].atomnum <= 18)
198     lrow = 2;
199     else if (atom[ic].atomnum <= 36)
200     lrow = 3;
201     else if (atom[ic].atomnum <= 54)
202     lrow = 4;
203    
204     for (j=0; j < 30; j++)
205     {
206     if ( irow == stbn[j].irow && krow == stbn[j].krow && lrow == stbn[j].lrow )
207     {
208     strbnd.isb[strbnd.nstrbnd][0] = i;
209     strbnd.isb[strbnd.nstrbnd][1] = nb1;
210     strbnd.isb[strbnd.nstrbnd][2] = nb2;
211     strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con1;
212     strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con2;
213     strbnd.nstrbnd++;
214     break;
215     } else if (lrow == stbn[j].irow && krow == stbn[j].krow && irow == stbn[j].lrow)
216     {
217     strbnd.isb[strbnd.nstrbnd][0] = i;
218     strbnd.isb[strbnd.nstrbnd][1] = nb1;
219     strbnd.isb[strbnd.nstrbnd][2] = nb2;
220     strbnd.ksb1[strbnd.nstrbnd] = stbn[j].con2;
221     strbnd.ksb2[strbnd.nstrbnd] = stbn[j].con1;
222     strbnd.nstrbnd++;
223     break;
224     }
225     }
226     }else
227     {
228     strbnd.isb[strbnd.nstrbnd][0] = i;
229     strbnd.isb[strbnd.nstrbnd][1] = nb1;
230     strbnd.isb[strbnd.nstrbnd][2] = nb2;
231     if (angles.anat[i] > 175.0)
232     {
233     strbnd.ksb1[strbnd.nstrbnd] = 0.0;
234     strbnd.ksb2[strbnd.nstrbnd] = 0.0;
235     } else
236     {
237     if (ita < itc)
238     {
239     strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
240     strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
241     } else
242     {
243     strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][1];
244     strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][0];
245     }
246     }
247     strbnd.nstrbnd++;
248     }
249     }
250     } else
251     {
252     k = 0;
253     if ( atom[ia].atomnum <= 1)
254     {
255     k++;
256     if (crossterm_k.stbncon[nstbn][2] == 0.0)
257     nb1 = -1;
258     }
259     if (atom[ic].atomnum <= 1)
260     {
261     k++;
262     if (crossterm_k.stbncon[nstbn][2] == 0.0)
263     nb2 = -1;
264     }
265     if (crossterm_k.stbncon[nstbn][k] != 0.0)
266     {
267     strbnd.isb[strbnd.nstrbnd][0] = i;
268     strbnd.isb[strbnd.nstrbnd][1] = nb1;
269     strbnd.isb[strbnd.nstrbnd][2] = nb2;
270     strbnd.ksb1[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];
271     strbnd.ksb2[strbnd.nstrbnd] = crossterm_k.stbncon[nstbn][k];
272     strbnd.nstrbnd++;
273     }
274     }
275     // if (ierr == FALSE)
276     // {
277     // Missing_constants = TRUE;
278     // fprintf(errfile,"Missing strbnd constants for %d %d %d \n",ia,ib,ic);
279     // }
280     L_10:
281     continue;
282     }
283     }
284    
285