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