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