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