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