ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/kcharge.c
Revision: 103
Committed: Thu Feb 19 01:37:38 2009 UTC (12 years, 4 months ago) by gilbertke
File size: 16623 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 #include "atom_k.h"
7 #include "utility.h"
8 #include "job_control.h"
9
10 #define MAXOPLS 600
11
12 void numeral(int, char *, int);
13 void compute_fcharge(void);
14 int find_rsize(int,int);
15 void search_rings(int);
16 void message_alert(char *,char *);
17 void get_rsize(int,int,int, int *);
18 int get_field(void);
19
20 EXTERN struct ElementType {
21 char symbol[3];
22 int atomnum;
23 float weight, covradius, vdwradius;
24 int s,p,d,f, type;
25 } Elements[] ;
26
27
28 EXTERN struct t_charge_k {
29 int ncharge, nbndchrg, nbndchrgdel;
30 int type[MAXATOMTYPE], btype[MAXBONDCONST], btypedel[MAXBONDCONST];
31 float charge[MAXATOMTYPE], bcharge[MAXBONDCONST], formchrg[MAXATOMTYPE], bchargedel[MAXBONDCONST];
32 float typechrg[MAXATOMTYPE];
33 } charge_k;
34
35 EXTERN struct t_dipole_k {
36 int ndipole,ndipole3,ndipole4,ndipole5;
37 char kb[MAXBONDCONST][7], kb3[MAXBOND3CONST][7], kb4[MAXBOND4CONST][7], kb5[MAXBOND5CONST][7];
38 float bmom[MAXBONDCONST],bmom3[MAXBOND3CONST],bmom4[MAXBOND4CONST],bmom5[MAXBOND5CONST];
39 } dipole_k;
40
41 void kcharge()
42 {
43 int i, j, ia, ib, it, iit, kit, jji,field;
44 long int mask;
45 char pa[4],pb[4],pt[7];
46 float bmomj;
47
48 mask = 1L << 0;
49 field = get_field();
50 for (i=1; i <= natom; i++)
51 {
52 atom.sigma_charge[i] = 0.0;
53 if (field == MMX && (atom.type[i] < 300))
54 atom.charge[i] = 0.0;
55 else if (field != MMX && atom.type[i] < 300)
56 atom.charge[i] = 0.0;
57 }
58
59 if (field == MMFF94)
60 {
61 compute_fcharge();
62 for(i=0; i < bonds_ff.nbnd; i++)
63 {
64 ia = bonds_ff.i12[i][0];
65 ib = bonds_ff.i12[i][1];
66 iit = atom.type[ia];
67 kit = atom.type[ib];
68 if (iit < kit)
69 it = iit*100+kit;
70 else
71 it = kit*100+iit;
72 if (bonds_ff.index[i] == 1)
73 {
74 for(j=0; j < charge_k.nbndchrgdel; j++)
75 {
76 if (it == charge_k.btypedel[j])
77 {
78 if (iit < kit)
79 {
80 atom.charge[ia] -= charge_k.bchargedel[j];
81 atom.charge[ib] += charge_k.bchargedel[j];
82 }else
83 {
84 atom.charge[ib] -= charge_k.bchargedel[j];
85 atom.charge[ia] += charge_k.bchargedel[j];
86 }
87 break;
88 }
89 }
90 } else
91 {
92 for(j=0; j < charge_k.nbndchrg; j++)
93 {
94 if (it == charge_k.btype[j])
95 {
96 if (iit < kit)
97 {
98 atom.charge[ia] -= charge_k.bcharge[j];
99 atom.charge[ib] += charge_k.bcharge[j];
100 }else
101 {
102 atom.charge[ib] -= charge_k.bcharge[j];
103 atom.charge[ia] += charge_k.bcharge[j];
104 }
105 break;
106 }
107 }
108 }
109 }
110 // add in formal charges
111 for (i=1; i <= natom; i++)
112 {
113 if (atom.atomnum[i] == 15)
114 atom.formal_charge[i] = 0.0;
115 if (atom.atomnum[i] == 16 && atom.type[i] != 72)
116 atom.formal_charge[i] = 0.0;
117 }
118 for (i=1; i <= natom; i++)
119 {
120 atom.charge[i] += (1.0 - atom_k.ligands[atom.type[i]]*charge_k.formchrg[atom.type[i]])*atom.formal_charge[i] ;
121 for(j=0; j < MAXIAT; j++)
122 {
123 if (atom.iat[i][j] != 0)
124 {
125 if (atom.formal_charge[atom.iat[i][j]] < 0.0)
126 atom.charge[i] += charge_k.formchrg[atom.type[atom.iat[i][j]]]*atom.formal_charge[atom.iat[i][j]];
127 }
128 }
129 }
130 } else if (field == MMX || field == MM3)
131 {
132 // compute charges from bond moments
133 for (i=0; i < bonds_ff.nbnd; i++)
134 {
135 ia = bonds_ff.i12[i][0];
136 ib = bonds_ff.i12[i][1];
137 iit = atom.type[ia];
138 kit = atom.type[ib];
139 if( field == MMX)
140 {
141 if (iit >= 300)
142 iit = 300;
143 if( kit >= 300)
144 kit = 300;
145 if (iit == 40 && ( kit != 2 && kit != 3 && kit != 4 && kit != 40) )
146 iit = 2;
147 if (kit == 40 && ( iit != 2 && iit != 3 && iit != 4 && iit != 40) )
148 kit = 2;
149 }
150 numeral(iit,pa,3);
151 numeral(kit,pb,3);
152
153 if (iit < kit)
154 {
155 strcpy(pt,pa);
156 strcat(pt,pb);
157 } else
158 {
159 strcpy(pt,pb);
160 strcat(pt,pa);
161 }
162
163 bmomj = 0.0;
164 for (j=0; j < dipole_k.ndipole; j++)
165 {
166 if (strcmp(pt,dipole_k.kb[j]) == 0)
167 {
168 bmomj = dipole_k.bmom[j];
169 break;
170 }
171 }
172
173 if (bmomj != 0.0)
174 {
175 bmomj = 0.42*bmomj/bonds_ff.bl[i];
176 if (iit <= kit)
177 {
178 atom.charge[ia] += 0.5*bmomj;
179 atom.charge[ib] -= 0.5*bmomj;
180 } else
181 {
182 atom.charge[ia] -= 0.5*bmomj;
183 atom.charge[ib] += 0.5*bmomj;
184 }
185 }
186 if (kit == 66)
187 {
188 if (iit == 3)
189 atom.charge[ib] -= 0.5;
190 if (iit == 18)
191 atom.charge[ib] -= 0.33;
192 } else if (iit == 66)
193 {
194 if (kit == 3)
195 atom.charge[ia] -= 0.5;
196 if (kit == 18)
197 atom.charge[ia] -= 0.33;
198 }
199 }
200 // check for isolated charged atoms
201 for (i=1; i <= natom; i++)
202 {
203 if ( !(atom.flags[i] & mask))
204 {
205 if (atom.type[i] == 16 || atom.type[i] == 30)
206 atom.charge[i] += 1.0;
207 if (atom.type[i] == 41) // N+
208 {
209 for(j=0; j < MAXIAT; j++)
210 {
211 if (atom.iat[i][j] != 0)
212 {
213 if (atom.type[atom.iat[i][j]] == 42 || atom.type[atom.iat[i][j]] == 4
214 || atom.type[atom.iat[i][j]] == 66)
215 goto L_10;
216 }
217 }
218 atom.charge[i] += 1.0;
219 }
220 if (atom.type[i] == 46) // O+
221 {
222 for (j=0; j < MAXIAT; j++)
223 {
224 if (atom.bo[i][j] == 3) // co in metal carbonyl
225 goto L_10;
226 }
227 atom.charge[i] += 1.0;
228 }
229 if (atom.type[i] == 27 || atom.type[i] == 42 || atom.type[i] == 48)
230 atom.charge[i] -= 1.0;
231 if (atom.type[i] >= 11 && atom.type[i] <= 14)
232 {
233 jji = 0;
234 for (j=0; j < MAXIAT; j++)
235 {
236 if (atom.iat[i][j] != 0)
237 jji++;
238 }
239 if (jji == 0)
240 atom.charge[i] -= 1.0;
241 }
242 if (atom.type[i] == 25 || atom.type[i] == 47)
243 {
244 jji = 0;
245 for (j=0; j < MAXIAT; j++)
246 {
247 if (atom.iat[i][j] != 0)
248 jji++;
249 }
250 if (jji == 6)
251 atom.charge[i] -= 1.0;
252 }
253 if (atom.type[i] == 58) // Al
254 atom.charge[i] += 1.0;
255
256 }
257 L_10:
258 continue;
259 }
260 // search list of charged atoms read in
261 for (i=1; i <= natom; i++)
262 {
263 for(j=0; j < charge_k.ncharge; j++)
264 {
265 if (atom.type[i] == charge_k.type[j])
266 {
267 atom.charge[i] = charge_k.charge[j];
268 break;
269 }
270 }
271 }
272 } else
273 {
274 // search list of charges
275 for (i=1; i <= natom; i++)
276 {
277 for(j=0; j < charge_k.ncharge; j++)
278 {
279 if (atom.type[i] == charge_k.type[j])
280 {
281 atom.charge[i] = charge_k.charge[j];
282 break;
283 }
284 }
285 }
286 }
287 for (i=1;i <= natom; i++)
288 atom.sigma_charge[i] = atom.charge[i];
289
290 // scale or turn off charges based on control from PYMOL
291 if (job_control.use_charge)
292 {
293 for (i=1; i <= natom; i++)
294 atom.charge[i] *= job_control.scale;
295 }
296 }
297 /* -------------------------------------------------------- */
298 void compute_fcharge()
299 {
300 int i, j, jjbo;
301 int jatm, k, jji, attached[4];
302 int *used;
303 float charge;
304
305 used = ivector(1,natom+1);
306 for (i=1; i <= natom; i++)
307 used[i] = FALSE;
308
309 for (i=1; i <= natom; i++)
310 {
311 atom.formal_charge[i] = charge_k.typechrg[atom.type[i]];
312 }
313 // find variable charges: types 32, 72, 76, 81
314 for (i=1; i <= natom; i++)
315 {
316 if (used[i] == FALSE)
317 {
318 if (atom.formal_charge[i] != 0.0)
319 {
320 if (atom.type[i] == 32)
321 {
322 jatm = atom.iat[i][0];
323 jji = 0;
324 jjbo = 0;
325 for (k=0; k < MAXIAT; k++)
326 {
327 if (atom.iat[jatm][k] != 0 && atom.bo[jatm][k] != 9)
328 {
329 if (atom.type[atom.iat[jatm][k]] == 32)
330 {
331 if (atom.bo[jatm][k] == 1)
332 {
333 jjbo++;
334 attached[jji] = atom.iat[jatm][k];
335 jji++;
336 } else
337 {
338 attached[jji] = atom.iat[jatm][k];
339 jji++;
340 }
341 }
342 }
343 }
344 charge = -(double)jjbo/jji;
345 if (jji == 1) // only 1 type 32
346 atom.formal_charge[i] = 0.0;
347 else if (jji == 2)
348 {
349 if (atom.atomnum[jatm] == 15) // p-o treat as single negative charge;
350 charge = -1.00/jji;
351 if (atom.type[jatm] == 73)
352 {
353 charge = -0.500;
354 atom.formal_charge[attached[0]] = charge;
355 atom.formal_charge[attached[1]] = charge;
356 used[attached[0]] = TRUE;
357 used[attached[1]] = TRUE;
358 } else if (atom.atomnum[jatm] != 7)
359 {
360 atom.formal_charge[attached[0]] = charge;
361 atom.formal_charge[attached[1]] = charge;
362 used[attached[0]] = TRUE;
363 used[attached[1]] = TRUE;
364 }else
365 {
366 atom.formal_charge[attached[0]] = 0.0;
367 atom.formal_charge[attached[1]] = 0.0;
368 used[attached[0]] = TRUE;
369 used[attached[1]] = TRUE;
370 }
371 } else if (jji == 3)
372 {
373 atom.formal_charge[attached[0]] = charge;
374 atom.formal_charge[attached[1]] = charge;
375 atom.formal_charge[attached[2]] = charge;
376 used[attached[0]] = TRUE;
377 used[attached[1]] = TRUE;
378 used[attached[2]] = TRUE;
379 } else if (jji == 4)
380 {
381 atom.formal_charge[attached[0]] = charge;
382 atom.formal_charge[attached[1]] = charge;
383 atom.formal_charge[attached[2]] = charge;
384 atom.formal_charge[attached[3]] = charge;
385 used[attached[0]] = TRUE;
386 used[attached[1]] = TRUE;
387 used[attached[2]] = TRUE;
388 used[attached[3]] = TRUE;
389 }
390 }else if (atom.type[i] == 72)
391 {
392 jatm = atom.iat[i][0];
393 jji = 0;
394 jjbo = 0;
395 for (k=0; k < MAXIAT; k++)
396 {
397 if (atom.iat[jatm][k] != 0 && atom.bo[jatm][k] != 9)
398 {
399 if (atom.type[atom.iat[jatm][k]] == 72)
400 {
401 if (atom.bo[jatm][k] == 1)
402 {
403 jjbo++;
404 attached[jji] = atom.iat[jatm][k];
405 jji++;
406 } else
407 {
408 attached[jji] = atom.iat[jatm][k];
409 jji++;
410 }
411 }
412 }
413 }
414 if (jji == 1)
415 {
416 if (atom.atomnum[jatm] == 15)
417 atom.formal_charge[i] = 0.0;
418 else
419 atom.formal_charge[i] = -1.00;
420 }else
421 {
422 atom.formal_charge[attached[0]] = -0.5;
423 atom.formal_charge[attached[1]] = -0.5;
424 used[attached[0]] = TRUE;
425 used[attached[1]] = TRUE;
426 }
427 } else if (atom.type[i] == 81)
428 {
429 jji = 0;
430 for (k=0; k < MAXIAT ; k++)
431 {
432 if (atom.iat[i][k] != 0 && atom.type[atom.iat[i][k]] == 80)
433 {
434 jatm = atom.iat[i][k];
435 jji = 0;
436 for (j=0; j < MAXIAT; j++)
437 {
438 if (atom.iat[jatm][j] != 0 && atom.bo[jatm][j] != 9)
439 {
440 if (atom.atomnum[atom.iat[jatm][j]] == 7)
441 {
442 attached[jji] = atom.iat[jatm][j];
443 jji++;
444 }
445 }
446 }
447 }
448 }
449 if (jji == 0)
450 {
451 charge = 1.00;
452 atom.formal_charge[i] = charge;
453 used[i] = TRUE;
454 } else
455 {
456 charge = 1.00/jji;
457 for (j=0; j < jji; j++)
458 {
459 atom.formal_charge[attached[j]] = charge;
460 used[attached[j]] = TRUE;
461 }
462 }
463 }
464 }
465 }
466 }
467 free_ivector(used,1, natom+1);
468 }