ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/kbond.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 8 months ago) by wdelano
File size: 10614 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5 #include "bonds_ff.h"
6 #include "rings.h"
7 #include "field.h"
8
9 int is_ring31(int);
10 int is_ring41(int);
11 int is_ring51(int);
12 int is_ring61(int);
13 int is_ring42(int, int);
14 int is_ring52(int, int);
15 int is_ring62(int,int);
16 void numeral(int, char *, int);
17 int is_delocalbond(int ia, int ib);
18 void message_alert(char *, char *);
19 double SIGN(double ,double );
20 int is_bond(int,int);
21 float get_bond(int, int);
22
23
24 EXTERN struct t_bondk1 {
25 int use_ring3, use_ring4, use_ring5;
26 int nbnd, nbnd3, nbnd4, nbnd5, ndeloc;
27 char kb[MAXBONDCONST][7], kb3[MAXBOND3CONST][7],kb4[MAXBOND4CONST][7],kb5[MAXBOND5CONST][7];
28 char kbdel[MAXBONDDELOC][7];
29 float s[MAXBONDCONST], t[MAXBONDCONST];
30 float s3[MAXBOND3CONST], t3[MAXBOND3CONST];
31 float s4[MAXBOND4CONST], t4[MAXBOND4CONST];
32 float s5[MAXBOND5CONST], t5[MAXBOND5CONST];
33 float sdel[MAXBONDDELOC], tdel[MAXBONDDELOC];
34 } bondk1;
35
36 EXTERN int Missing_constants;
37
38 int kbond()
39 {
40 long int pi_mask;
41 int ia4,ia5,ia6, ib4,ib5,ib6, iab4, iab5, iab6;
42 int i, j,k, iit, kit, it, ierr,nbpi, ii, kk, i1, i2;
43 int igeni, igenk, ifound;
44 int ia, ib, iend, kend, iy;
45 int numi, numj;
46 float bconst, blength;
47 float sslope, tslope, tslope2;
48 double dbl,dks,bl0;
49 char hatext[80];
50 char pa[4],pb[4], pt[7], pt_gen[7];
51
52 pi_mask = 1L << PI_MASK; /* this mask is for pi atoms - now uses bit 0 */
53 nbpi = 0;
54 if( natom != 0 )
55 {
56
57 for( i = 0; i < bonds_ff.nbnd; i++ )
58 {
59 ia = bonds_ff.i12[i][0];
60 ib = bonds_ff.i12[i][1];
61 iit = atom[bonds_ff.i12[i][0]].type;
62 kit = atom[bonds_ff.i12[i][1]].type;
63
64 igeni = atom[bonds_ff.i12[i][0]].tclass;
65 igenk = atom[bonds_ff.i12[i][1]].tclass;
66
67 // metal general class
68 if (iit > 300)
69 igeni = 300;
70 if (kit > 300)
71 igenk = 300;
72
73 if( field.type == MMX)
74 {
75 if (iit == 40 && ( kit != 2 && kit != 3 && kit != 4 && kit != 40) )
76 iit = 2;
77 if (kit == 40 && ( iit != 2 && iit != 3 && iit != 4 && iit != 40) )
78 kit = 2;
79 }
80
81 numeral(iit, pa, 3);
82 numeral(kit, pb, 3);
83 strcpy(pt,pa);
84 strcat(pt,pb);
85 if( iit <= kit )
86 {
87 strcpy(pt,pa);
88 strcat(pt,pb);
89 }else
90 {
91 strcpy(pt,pb);
92 strcat(pt,pa);
93 }
94 numeral(igeni, pa, 3);
95 numeral(igenk, pb, 3);
96 if( igeni <= igenk )
97 {
98 strcpy(pt_gen,pa);
99 strcat(pt_gen,pb);
100 }else
101 {
102 strcpy(pt_gen,pb);
103 strcat(pt_gen,pa);
104 }
105
106 ierr = FALSE;
107 // check 3 membered ring
108 if (rings.nring3 > 0 && bondk1.nbnd3 > 0 && ( is_ring31(ia) && is_ring31(ib) ) )
109 {
110 for (j=0; j < bondk1.nbnd3; j++)
111 {
112 if ( (strcmp(bondk1.kb3[j], pt) == 0) || (strcmp(bondk1.kb3[j],pt_gen) == 0) )
113 {
114 bonds_ff.bk[i] = bondk1.s3[j];
115 bonds_ff.bl[i] = bondk1.t3[j];
116 bonds_ff.index[i] = 0;
117 ierr = TRUE;
118 break;
119 }
120 }
121 if (ierr != TRUE)
122 {
123 Missing_constants = TRUE;
124 sprintf(hatext,"Bond constants missing for cyclopropane bond %d - %d of type %d-%d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
125 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
126 fprintf(pcmlogfile,"%s\n",hatext);
127 //fprintf(stderr,"%s\n",hatext);
128 }
129 }
130 if (ierr == TRUE)
131 goto L_10;
132 // check 4 membered ring
133 if (rings.nring4 > 0 && bondk1.nbnd4 > 0 && ( is_ring41(ia) && is_ring41(ib) ) )
134 {
135 for (j=0; j < bondk1.nbnd4; j++)
136 {
137 if ( (strcmp(bondk1.kb4[j],pt) == 0) || (strcmp(bondk1.kb4[j],pt_gen) == 0) )
138 {
139 bonds_ff.bk[i] = bondk1.s4[j];
140 bonds_ff.bl[i] = bondk1.t4[j];
141 bonds_ff.index[i] = 0;
142 ierr = TRUE;
143 break;
144 }
145 }
146 if (ierr != TRUE)
147 {
148 Missing_constants = TRUE;
149 sprintf(hatext,"Bond constants missing for cyclobutane bond %d - %d of type %d-%d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
150 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
151 fprintf(pcmlogfile,"%s\n",hatext);
152 // fprintf(stderr,"%s\n",hatext);
153 }
154 }
155 if (ierr == TRUE)
156 goto L_10;
157 // check 5 membered ring
158 if (rings.nring5 > 0 && bondk1.nbnd5 > 0 && (is_ring51(ia) != -1) && (is_ring51(ib) != -1))
159 {
160 for (j=0; j < bondk1.nbnd5; j++)
161 {
162 if (strcmp(bondk1.kb5[j],pt) == 0 )
163 {
164 bonds_ff.bk[i] = bondk1.s5[j];
165 bonds_ff.bl[i] = bondk1.t5[j];
166 bonds_ff.index[i] = 0;
167 ierr = TRUE;
168 break;
169 }
170 }
171 // don't fail on missing param - check for regular param first
172 }
173 if (ierr == TRUE)
174 goto L_10;
175 // delocalized bonds in MMFF94
176 if (field.type == MMFF94 && bondk1.ndeloc > 0)
177 {
178 if ( is_delocalbond(ia, ib) )
179 {
180 for (k = 0; k < MAXIAT; k++)
181 {
182 if (atom[ia].iat[k] == ib)
183 {
184 if (atom[ia].bo[k] == 1) // single bond
185 {
186 for (j=0; j < bondk1.ndeloc; j++)
187 {
188 if (strcmp(bondk1.kbdel[j],pt) == 0)
189 {
190 bonds_ff.bk[i] = bondk1.sdel[j];
191 bonds_ff.bl[i] = bondk1.tdel[j];
192 bonds_ff.index[i] = 1;
193 ierr = TRUE;
194 break;
195 }
196 }
197 }
198 }
199 }
200 }
201 if (ierr == TRUE)
202 goto L_10;
203 }
204 // regular parameters
205 for (j=0; j < bondk1.nbnd; j++)
206 {
207 if ( (strcmp(bondk1.kb[j], pt) == 0) )
208 {
209 bonds_ff.bk[i] = bondk1.s[j];
210 bonds_ff.bl[i] = bondk1.t[j];
211 bonds_ff.index[i] = 0;
212 ierr = TRUE;
213 if ( ((strcmp(pt," 2 2") == 0) || (strcmp(pt," 2 3") == 0))&& field.type == MMX)
214 {
215 for (k=0; k < MAXIAT; k++)
216 {
217 if (atom[ia].iat[k] == ib)
218 {
219 if (atom[ia].bo[k] == 1) // butadiene
220 {
221 bonds_ff.bl[i] = 1.48;
222 break;
223 }
224 }
225 }
226 }
227 break;
228 }
229 }
230 if (ierr == TRUE)
231 goto L_10;
232 //
233 // generalized parameters
234 for (j=0; j < bondk1.nbnd; j++)
235 {
236 if ( (strcmp(bondk1.kb[j], pt_gen) == 0) )
237 {
238 bonds_ff.bk[i] = bondk1.s[j];
239 bonds_ff.bl[i] = bondk1.t[j];
240 bonds_ff.index[i] = 0;
241 ierr = TRUE;
242 break;
243 }
244 }
245 if (ierr == TRUE)
246 goto L_10;
247 if (ierr != TRUE)
248 {
249 Missing_constants = TRUE;
250 sprintf(hatext,"Bond constants missing for bond %d - %d of type %d - %d\n",bonds_ff.i12[i][0],bonds_ff.i12[i][1],
251 atom[bonds_ff.i12[i][0]].type, atom[bonds_ff.i12[i][1]].type);
252 fprintf(pcmlogfile,"%s\n",hatext);
253 }
254 L_10:
255 continue;
256
257 }
258 }
259 if (Missing_constants == TRUE)
260 {
261 fprintf(pcmlogfile,"Error assigning constants in: %s\n",Struct_Title);
262 fprintf(pcmlogfile,"%s\n",hatext);
263 return FALSE;
264 }
265
266 return TRUE;
267 }
268 /* ------------------------------------------ */
269 int find_bond(int ia, int ib)
270 {
271 int i;
272 int iz;
273
274 iz = -1;
275 for (i=0; i < bonds_ff.nbnd; i++)
276 {
277 if (ia == bonds_ff.i12[i][0] && ib == bonds_ff.i12[i][1])
278 return(i);
279 if (ib == bonds_ff.i12[i][0] && ia == bonds_ff.i12[i][1])
280 return(i);
281 }
282 return(iz);
283 }
284
285 /* ------------------------------------------ */
286 int is_delocalbond(int ia, int ib)
287 {
288 // test for delocalized bond
289 // if ia & ib are aromatic and part of 6 membered ring
290 // if not aromatic is bond order 1
291
292 long int aromatic_mask;
293 int i, j;
294 int jdbl, kdbl;
295
296 aromatic_mask = (1L << AROMATIC_MASK);
297
298 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
299 && (is_ring62(ia,ib) == TRUE) )
300 {
301 return(FALSE);
302 }
303
304 if ( (atom[ia].flags & aromatic_mask) && (atom[ib].flags & aromatic_mask)
305 && (is_ring52(ia,ib) == TRUE) )
306 {
307 return(FALSE);
308 }
309 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring62(ia,ib) != FALSE))
310 return FALSE;
311 if (atom[ia].type == 37 && atom[ib].type == 37 && (is_ring52(ia,ib) != FALSE))
312 return FALSE;
313
314 for (i=0; i < MAXIAT; i++)
315 {
316 if (atom[ia].iat[i] == ib)
317 {
318 if (atom[ia].bo[i] == 1)
319 {
320 jdbl = FALSE;
321 kdbl = FALSE;
322 for (j=0; j < MAXIAT; j++)
323 {
324 if (atom[ia].bo[j] >= 2 && atom[ia].bo[j] != 9)
325 jdbl = TRUE;
326 if (atom[ib].bo[j] >= 2 && atom[ib].bo[j] != 9)
327 kdbl = TRUE;
328 }
329 if (jdbl == TRUE && kdbl == TRUE)
330 return (TRUE);
331 else if (jdbl == TRUE && (atom[ib].flags & aromatic_mask) )
332 return (TRUE);
333 else if (kdbl == TRUE && (atom[ia].flags & aromatic_mask) )
334 return (TRUE);
335 else
336 return(FALSE);
337 }
338 }
339 }
340 return(FALSE);
341 }
342 #define IS_ODD(j) ((j) & 1 )
343
344 double SIGN(double a,double b ) /* floating polong int SIGN transfer */
345 {
346 return( b < 0.0 ? -fabs(a) : fabs(a) );
347 }
348