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