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, 10 months ago) by gilbertke
File size: 16414 byte(s)
Log Message:
combine rings and srings to eliminate duplication
Line User Rev File contents
1 wdelano 58 #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 gilbertke 98
162 wdelano 58 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 gilbertke 98
172 wdelano 58 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     }