ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/kcharge.c
Revision: 93
Committed: Sat Jan 17 23:24:55 2009 UTC (13 years, 4 months ago) by wdelano
File size: 19386 byte(s)
Log Message:
branch created
Line File contents
1 #define EXTERN extern
2
3 #include "pcwin.h"
4 #include "pcmod.h"
5 #include "bonds_ff.h"
6 #include "field.h"
7 #include "atom_k.h"
8 #include "utility.h"
9 #include "pdb.h"
10
11 #define MAXOPLS 600
12
13 void numeral(int, char *, int);
14 void compute_fcharge(void);
15 int find_rsize(int,int);
16 void search_rings(int);
17 void message_alert(char *,char *);
18 void get_rsize(int,int,int, int *);
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_bondpk {
29 int npibond,npibond44i,npibond45i,npibond46i,npibond55i,npibond56i,npibond66i;
30 int npibond44o,npibond45o,npibond46o,npibond55o,npibond56o,npibond66o;
31 int npibond40,npibond50,npibond60;
32 int use_pibond,use_pibond44i,use_pibond45i,use_pibond46i,use_pibond55i,use_pibond56i,use_pibond66i;
33 int use_pibond44o,use_pibond45o,use_pibond46o,use_pibond55o,use_pibond56o,use_pibond66o;
34 int use_pibond40,use_pibond50,use_pibond60;
35 char kb[MAXPIBONDCONST][7],kb44i[20][7],kb45i[20][7],kb46i[20][7],kb55i[20][7],kb56i[20][7],kb66i[20][7];
36 char kb44o[20][7],kb45o[20][7],kb46o[20][7],kb55o[20][7],kb56o[20][7],kb66o[20][7];
37 char kb40[20][7],kb50[20][7],kb60[20][7];
38 float bk[MAXPIBONDCONST], bl[MAXPIBONDCONST], bmom[MAXPIBONDCONST],sslop[MAXPIBONDCONST], tslop[MAXPIBONDCONST], tslop2[MAXPIBONDCONST];
39 float bk44i[20], bl44i[20], bmom44i[20],sslop44i[20], tslop44i[20], tslop244i[20];
40 float bk45i[20], bl45i[20], bmom45i[20],sslop45i[20], tslop45i[20], tslop245i[20];
41 float bk46i[20], bl46i[20], bmom46i[20],sslop46i[20], tslop46i[20], tslop246i[20];
42 float bk55i[20], bl55i[20], bmom55i[20],sslop55i[20], tslop55i[20], tslop255i[20];
43 float bk56i[20], bl56i[20], bmom56i[20],sslop56i[20], tslop56i[20], tslop256i[20];
44 float bk66i[20], bl66i[20], bmom66i[20],sslop66i[20], tslop66i[20], tslop266i[20];
45 float bk44o[20], bl44o[20], bmom44o[20],sslop44o[20], tslop44o[20], tslop244o[20];
46 float bk45o[20], bl45o[20], bmom45o[20],sslop45o[20], tslop45o[20], tslop245o[20];
47 float bk46o[20], bl46o[20], bmom46o[20],sslop46o[20], tslop46o[20], tslop246o[20];
48 float bk55o[20], bl55o[20], bmom55o[20],sslop55o[20], tslop55o[20], tslop255o[20];
49 float bk56o[20], bl56o[20], bmom56o[20],sslop56o[20], tslop56o[20], tslop256o[20];
50 float bk66o[20], bl66o[20], bmom66o[20],sslop66o[20], tslop66o[20], tslop266o[20];
51 float bk40[20], bl40[20], bmom40[20],sslop40[20], tslop40[20], tslop240[20];
52 float bk50[20], bl50[20], bmom50[20],sslop50[20], tslop50[20], tslop250[20];
53 float bk60[20], bl60[20], bmom60[20],sslop60[20], tslop60[20], tslop260[20];
54 } bondpk;
55
56 EXTERN struct t_charge_k {
57 int ncharge, nbndchrg, nbndchrgdel;
58 int type[MAXATOMTYPE], btype[MAXBONDCONST], btypedel[MAXBONDCONST];
59 float charge[MAXATOMTYPE], bcharge[MAXBONDCONST], formchrg[MAXATOMTYPE], bchargedel[MAXBONDCONST];
60 float typechrg[MAXATOMTYPE];
61 } charge_k;
62
63 EXTERN struct t_dipole_k {
64 int ndipole,ndipole3,ndipole4,ndipole5;
65 char kb[MAXBONDCONST][7], kb3[MAXBOND3CONST][7], kb4[MAXBOND4CONST][7], kb5[MAXBOND5CONST][7];
66 float bmom[MAXBONDCONST],bmom3[MAXBOND3CONST],bmom4[MAXBOND4CONST],bmom5[MAXBOND5CONST];
67 } dipole_k;
68
69 EXTERN struct t_ringdata {
70 int nring3, nring4, nring5, nring6;
71 int tot3, tot4, tot5, tot6;
72 int **iring3,**iring4,**iring5,**iring6;
73 } ringdata;
74
75
76 void kcharge()
77 {
78 int i, j, ia, ib, it, iit, kit, jji;
79 long int mask;
80 char pa[4],pb[4],pt[7];
81 float bmomj;
82
83 mask = 1L << 0;
84
85 for (i=1; i <= natom; i++)
86 {
87 atom[i].sigma_charge = 0.0;
88 if (field.type == MMX && (atom[i].type < 300))
89 atom[i].charge = 0.0;
90 else if (field.type != MMX && atom[i].type < 300)
91 atom[i].charge = 0.0;
92 }
93
94 if (field.type == MMFF94)
95 {
96 compute_fcharge();
97 for(i=0; i < bonds_ff.nbnd; i++)
98 {
99 ia = bonds_ff.i12[i][0];
100 ib = bonds_ff.i12[i][1];
101 iit = atom[ia].type;
102 kit = atom[ib].type;
103 if (iit < kit)
104 it = iit*100+kit;
105 else
106 it = kit*100+iit;
107 if (bonds_ff.index[i] == 1)
108 {
109 for(j=0; j < charge_k.nbndchrgdel; j++)
110 {
111 if (it == charge_k.btypedel[j])
112 {
113 if (iit < kit)
114 {
115 atom[ia].charge -= charge_k.bchargedel[j];
116 atom[ib].charge += charge_k.bchargedel[j];
117 }else
118 {
119 atom[ib].charge -= charge_k.bchargedel[j];
120 atom[ia].charge += charge_k.bchargedel[j];
121 }
122 break;
123 }
124 }
125 } else
126 {
127 for(j=0; j < charge_k.nbndchrg; j++)
128 {
129 if (it == charge_k.btype[j])
130 {
131 if (iit < kit)
132 {
133 atom[ia].charge -= charge_k.bcharge[j];
134 atom[ib].charge += charge_k.bcharge[j];
135 }else
136 {
137 atom[ib].charge -= charge_k.bcharge[j];
138 atom[ia].charge += charge_k.bcharge[j];
139 }
140 break;
141 }
142 }
143 }
144 }
145 // add in formal charges
146 for (i=1; i <= natom; i++)
147 {
148 if (atom[i].atomnum == 15)
149 atom[i].formal_charge = 0.0;
150 if (atom[i].atomnum == 16 && atom[i].type != 72)
151 atom[i].formal_charge = 0.0;
152 }
153 for (i=1; i <= natom; i++)
154 {
155 atom[i].charge += (1.0 - atom_k.ligands[atom[i].type]*charge_k.formchrg[atom[i].type])*atom[i].formal_charge ;
156 for(j=0; j < MAXIAT; j++)
157 {
158 if (atom[i].iat[j] != 0)
159 {
160 if (atom[atom[i].iat[j]].formal_charge < 0.0)
161 atom[i].charge += charge_k.formchrg[atom[atom[i].iat[j]].type]*atom[atom[i].iat[j]].formal_charge;
162 }
163 }
164 }
165 } else if (field.type == MMX || field.type == MM3)
166 {
167 // compute charges from bond moments
168 for (i=0; i < bonds_ff.nbnd; i++)
169 {
170 ia = bonds_ff.i12[i][0];
171 ib = bonds_ff.i12[i][1];
172 iit = atom[ia].type;
173 kit = atom[ib].type;
174 if( field.type == MMX)
175 {
176 if (iit >= 300)
177 iit = 300;
178 if( kit >= 300)
179 kit = 300;
180 if (iit == 40 && ( kit != 2 && kit != 3 && kit != 4 && kit != 40) )
181 iit = 2;
182 if (kit == 40 && ( iit != 2 && iit != 3 && iit != 4 && iit != 40) )
183 kit = 2;
184 }
185 numeral(iit,pa,3);
186 numeral(kit,pb,3);
187
188 if (iit < kit)
189 {
190 strcpy(pt,pa);
191 strcat(pt,pb);
192 } else
193 {
194 strcpy(pt,pb);
195 strcat(pt,pa);
196 }
197 if( (atom[ia].flags & mask) && (atom[ib].flags & mask) )
198 { // search through pi constants
199 bmomj = 0.0;
200 for (j=0; j < bondpk.npibond; j++)
201 {
202 if (strcmp(pt, bondpk.kb[j]) == 0)
203 {
204 bmomj = bondpk.bmom[j];
205 break;
206 }
207 }
208 } else
209 {
210 bmomj = 0.0;
211 for (j=0; j < dipole_k.ndipole; j++)
212 {
213 if (strcmp(pt,dipole_k.kb[j]) == 0)
214 {
215 bmomj = dipole_k.bmom[j];
216 break;
217 }
218 }
219 }
220 if (bmomj != 0.0)
221 {
222 bmomj = 0.42*bmomj/bonds_ff.bl[i];
223 if (iit <= kit)
224 {
225 atom[ia].charge += 0.5*bmomj;
226 atom[ib].charge -= 0.5*bmomj;
227 } else
228 {
229 atom[ia].charge -= 0.5*bmomj;
230 atom[ib].charge += 0.5*bmomj;
231 }
232 }
233 if (kit == 66)
234 {
235 if (iit == 3)
236 atom[ib].charge -= 0.5;
237 if (iit == 18)
238 atom[ib].charge -= 0.33;
239 } else if (iit == 66)
240 {
241 if (kit == 3)
242 atom[ia].charge -= 0.5;
243 if (kit == 18)
244 atom[ia].charge -= 0.33;
245 }
246 }
247 // check for isolated charged atoms
248 for (i=1; i <= natom; i++)
249 {
250 if ( !(atom[i].flags & mask))
251 {
252 if (atom[i].type == 16 || atom[i].type == 30)
253 atom[i].charge += 1.0;
254 if (atom[i].type == 41) // N+
255 {
256 for(j=0; j < MAXIAT; j++)
257 {
258 if (atom[i].iat[j] != 0)
259 {
260 if (atom[atom[i].iat[j]].type == 42 || atom[atom[i].iat[j]].type == 4
261 || atom[atom[i].iat[j]].type == 66)
262 goto L_10;
263 }
264 }
265 atom[i].charge += 1.0;
266 }
267 if (atom[i].type == 46) // O+
268 {
269 for (j=0; j < MAXIAT; j++)
270 {
271 if (atom[i].bo[j] == 3) // co in metal carbonyl
272 goto L_10;
273 }
274 atom[i].charge += 1.0;
275 }
276 if (atom[i].type == 27 || atom[i].type == 42 || atom[i].type == 48)
277 atom[i].charge -= 1.0;
278 if (atom[i].type >= 11 && atom[i].type <= 14)
279 {
280 jji = 0;
281 for (j=0; j < MAXIAT; j++)
282 {
283 if (atom[i].iat[j] != 0)
284 jji++;
285 }
286 if (jji == 0)
287 atom[i].charge -= 1.0;
288 }
289 if (atom[i].type == 25 || atom[i].type == 47)
290 {
291 jji = 0;
292 for (j=0; j < MAXIAT; j++)
293 {
294 if (atom[i].iat[j] != 0)
295 jji++;
296 }
297 if (jji == 6)
298 atom[i].charge -= 1.0;
299 }
300 if (atom[i].type == 58) // Al
301 atom[i].charge += 1.0;
302
303 }
304 L_10:
305 continue;
306 }
307 // search list of charged atoms read in
308 for (i=1; i <= natom; i++)
309 {
310 for(j=0; j < charge_k.ncharge; j++)
311 {
312 if (atom[i].type == charge_k.type[j])
313 {
314 atom[i].charge = charge_k.charge[j];
315 break;
316 }
317 }
318 }
319 } else
320 {
321 // search list of charges
322 for (i=1; i <= natom; i++)
323 {
324 for(j=0; j < charge_k.ncharge; j++)
325 {
326 if (atom[i].type == charge_k.type[j])
327 {
328 atom[i].charge = charge_k.charge[j];
329 break;
330 }
331 }
332 }
333 }
334 for (i=1;i <= natom; i++)
335 atom[i].sigma_charge = atom[i].charge;
336 }
337 /* -------------------------------------------------------- */
338 void compute_fcharge()
339 {
340 int i, j, jjbo;
341 int jatm, k, jji, attached[4];
342 int *used;
343 float charge;
344
345 used = ivector(1,natom+1);
346 for (i=1; i <= natom; i++)
347 used[i] = FALSE;
348
349 for (i=1; i <= natom; i++)
350 {
351 atom[i].formal_charge = charge_k.typechrg[atom[i].type];
352 }
353 // find variable charges: types 32, 72, 76, 81
354 for (i=1; i <= natom; i++)
355 {
356 if (used[i] == FALSE)
357 {
358 if (atom[i].formal_charge != 0.0)
359 {
360 if (atom[i].type == 32)
361 {
362 jatm = atom[i].iat[0];
363 jji = 0;
364 jjbo = 0;
365 for (k=0; k < MAXIAT; k++)
366 {
367 if (atom[jatm].iat[k] != 0 && atom[jatm].bo[k] != 9)
368 {
369 if (atom[atom[jatm].iat[k]].type == 32)
370 {
371 if (atom[jatm].bo[k] == 1)
372 {
373 jjbo++;
374 attached[jji] = atom[jatm].iat[k];
375 jji++;
376 } else
377 {
378 attached[jji] = atom[jatm].iat[k];
379 jji++;
380 }
381 }
382 }
383 }
384 charge = -(double)jjbo/jji;
385 if (jji == 1) // only 1 type 32
386 atom[i].formal_charge = 0.0;
387 else if (jji == 2)
388 {
389 if (atom[jatm].atomnum == 15) // p-o treat as single negative charge;
390 charge = -1.00/jji;
391 if (atom[jatm].type == 73)
392 {
393 charge = -0.500;
394 atom[attached[0]].formal_charge = charge;
395 atom[attached[1]].formal_charge = charge;
396 used[attached[0]] = TRUE;
397 used[attached[1]] = TRUE;
398 } else if (atom[jatm].atomnum != 7)
399 {
400 atom[attached[0]].formal_charge = charge;
401 atom[attached[1]].formal_charge = charge;
402 used[attached[0]] = TRUE;
403 used[attached[1]] = TRUE;
404 }else
405 {
406 atom[attached[0]].formal_charge = 0.0;
407 atom[attached[1]].formal_charge = 0.0;
408 used[attached[0]] = TRUE;
409 used[attached[1]] = TRUE;
410 }
411 } else if (jji == 3)
412 {
413 atom[attached[0]].formal_charge = charge;
414 atom[attached[1]].formal_charge = charge;
415 atom[attached[2]].formal_charge = charge;
416 used[attached[0]] = TRUE;
417 used[attached[1]] = TRUE;
418 used[attached[2]] = TRUE;
419 } else if (jji == 4)
420 {
421 atom[attached[0]].formal_charge = charge;
422 atom[attached[1]].formal_charge = charge;
423 atom[attached[2]].formal_charge = charge;
424 atom[attached[3]].formal_charge = charge;
425 used[attached[0]] = TRUE;
426 used[attached[1]] = TRUE;
427 used[attached[2]] = TRUE;
428 used[attached[3]] = TRUE;
429 }
430 }else if (atom[i].type == 72)
431 {
432 jatm = atom[i].iat[0];
433 jji = 0;
434 jjbo = 0;
435 for (k=0; k < MAXIAT; k++)
436 {
437 if (atom[jatm].iat[k] != 0 && atom[jatm].bo[k] != 9)
438 {
439 if (atom[atom[jatm].iat[k]].type == 72)
440 {
441 if (atom[jatm].bo[k] == 1)
442 {
443 jjbo++;
444 attached[jji] = atom[jatm].iat[k];
445 jji++;
446 } else
447 {
448 attached[jji] = atom[jatm].iat[k];
449 jji++;
450 }
451 }
452 }
453 }
454 if (jji == 1)
455 {
456 if (atom[jatm].atomnum == 15)
457 atom[i].formal_charge = 0.0;
458 else
459 atom[i].formal_charge = -1.00;
460 }else
461 {
462 atom[attached[0]].formal_charge = -0.5;
463 atom[attached[1]].formal_charge = -0.5;
464 used[attached[0]] = TRUE;
465 used[attached[1]] = TRUE;
466 }
467 } else if (atom[i].type == 81)
468 {
469 jji = 0;
470 for (k=0; k < MAXIAT ; k++)
471 {
472 if (atom[i].iat[k] != 0 && atom[atom[i].iat[k]].type == 80)
473 {
474 jatm = atom[i].iat[k];
475 jji = 0;
476 for (j=0; j < MAXIAT; j++)
477 {
478 if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] != 9)
479 {
480 if (atom[atom[jatm].iat[j]].atomnum == 7)
481 {
482 attached[jji] = atom[jatm].iat[j];
483 jji++;
484 }
485 }
486 }
487 }
488 }
489 if (jji == 0)
490 {
491 charge = 1.00;
492 atom[i].formal_charge = charge;
493 used[i] = TRUE;
494 } else
495 {
496 charge = 1.00/jji;
497 for (j=0; j < jji; j++)
498 {
499 atom[attached[j]].formal_charge = charge;
500 used[attached[j]] = TRUE;
501 }
502 }
503 }
504 }
505 }
506 }
507 free_ivector(used,1, natom+1);
508 }