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