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