ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/kcharge.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 14826 byte(s)
Log Message:
test

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 #include "pdb.h"
10
11 #define MAXOPLS 600
12
13 void numeral(int, char *, int);
14 void compute_fcharge(void);
15 void assign_amber_chrg(int);
16 void assign_gast_charge();
17 void compute_opls(void);
18 int find_rsize(int,int);
19 void search_rings(int);
20 void message_alert(char *,char *);
21 void find_residues(void);
22 void bbchk(void);
23 void pireset(void);
24 int get_oplsconst(int,int,char *);
25 void ksortr(int, int *,int *,int *);
26 void get_rsize(int,int,int, int *);
27 void ksort_ring(int,int,int *);
28 void get_firstlevel(int,int *,int *,int *,int *);
29 void get_secondlevel(int,int,int *,int *,int *,int *);
30 void get_residue_type(int,int*);
31
32 EXTERN struct ElementType {
33 char symbol[3];
34 int atomnum;
35 float weight, covradius, vdwradius;
36 int s,p,d,f, type;
37 } Elements[] ;
38
39 EXTERN struct BioType {
40 char name_3[4], name_1[2];
41 int type, natom;
42 } BioGroups[];
43
44 typedef struct Gastype {
45 int type;
46 float a,b,c,d;
47 } Gastype;
48
49 EXTERN Gastype Gastchrg[];
50
51 EXTERN struct t_bondpk {
52 int npibond,npibond44i,npibond45i,npibond46i,npibond55i,npibond56i,npibond66i;
53 int npibond44o,npibond45o,npibond46o,npibond55o,npibond56o,npibond66o;
54 int npibond40,npibond50,npibond60;
55 int use_pibond,use_pibond44i,use_pibond45i,use_pibond46i,use_pibond55i,use_pibond56i,use_pibond66i;
56 int use_pibond44o,use_pibond45o,use_pibond46o,use_pibond55o,use_pibond56o,use_pibond66o;
57 int use_pibond40,use_pibond50,use_pibond60;
58 char kb[MAXPIBONDCONST][7],kb44i[20][7],kb45i[20][7],kb46i[20][7],kb55i[20][7],kb56i[20][7],kb66i[20][7];
59 char kb44o[20][7],kb45o[20][7],kb46o[20][7],kb55o[20][7],kb56o[20][7],kb66o[20][7];
60 char kb40[20][7],kb50[20][7],kb60[20][7];
61 float bk[MAXPIBONDCONST], bl[MAXPIBONDCONST], bmom[MAXPIBONDCONST],sslop[MAXPIBONDCONST], tslop[MAXPIBONDCONST], tslop2[MAXPIBONDCONST];
62 float bk44i[20], bl44i[20], bmom44i[20],sslop44i[20], tslop44i[20], tslop244i[20];
63 float bk45i[20], bl45i[20], bmom45i[20],sslop45i[20], tslop45i[20], tslop245i[20];
64 float bk46i[20], bl46i[20], bmom46i[20],sslop46i[20], tslop46i[20], tslop246i[20];
65 float bk55i[20], bl55i[20], bmom55i[20],sslop55i[20], tslop55i[20], tslop255i[20];
66 float bk56i[20], bl56i[20], bmom56i[20],sslop56i[20], tslop56i[20], tslop256i[20];
67 float bk66i[20], bl66i[20], bmom66i[20],sslop66i[20], tslop66i[20], tslop266i[20];
68 float bk44o[20], bl44o[20], bmom44o[20],sslop44o[20], tslop44o[20], tslop244o[20];
69 float bk45o[20], bl45o[20], bmom45o[20],sslop45o[20], tslop45o[20], tslop245o[20];
70 float bk46o[20], bl46o[20], bmom46o[20],sslop46o[20], tslop46o[20], tslop246o[20];
71 float bk55o[20], bl55o[20], bmom55o[20],sslop55o[20], tslop55o[20], tslop255o[20];
72 float bk56o[20], bl56o[20], bmom56o[20],sslop56o[20], tslop56o[20], tslop256o[20];
73 float bk66o[20], bl66o[20], bmom66o[20],sslop66o[20], tslop66o[20], tslop266o[20];
74 float bk40[20], bl40[20], bmom40[20],sslop40[20], tslop40[20], tslop240[20];
75 float bk50[20], bl50[20], bmom50[20],sslop50[20], tslop50[20], tslop250[20];
76 float bk60[20], bl60[20], bmom60[20],sslop60[20], tslop60[20], tslop260[20];
77 } bondpk;
78
79 EXTERN struct t_charge_k {
80 int ncharge, nbndchrg, nbndchrgdel;
81 int type[MAXATOMTYPE], btype[MAXBONDCONST], btypedel[MAXBONDCONST];
82 float charge[MAXATOMTYPE], bcharge[MAXBONDCONST], formchrg[MAXATOMTYPE], bchargedel[MAXBONDCONST];
83 float typechrg[MAXATOMTYPE];
84 } charge_k;
85
86 EXTERN struct t_amberchrg_k {
87 int ncharge, type[MAXATOMTYPE], res_type[MAXATOMTYPE];
88 char symbol[MAXATOMTYPE][4];
89 float chrg[MAXATOMTYPE];
90 } amberchrg_k;
91
92 EXTERN struct t_oplschrg_k {
93 int ncharge, type[MAXOPLS];
94 char chrgstring[MAXOPLS][25];
95 float chrg[MAXOPLS];
96 } oplschrg_k;
97
98 EXTERN struct t_dipole_k {
99 int ndipole,ndipole3,ndipole4,ndipole5;
100 char kb[MAXBONDCONST][7], kb3[MAXBOND3CONST][7], kb4[MAXBOND4CONST][7], kb5[MAXBOND5CONST][7];
101 float bmom[MAXBONDCONST],bmom3[MAXBOND3CONST],bmom4[MAXBOND4CONST],bmom5[MAXBOND5CONST];
102 } dipole_k;
103
104 EXTERN struct t_residues {
105 int nchain, ichainstart[10];
106 int ngroup, iresnum[200], irestype[200], istartatom[200];
107 } residues;
108 EXTERN struct t_ringdata {
109 int nring3, nring4, nring5, nring6;
110 int tot3, tot4, tot5, tot6;
111 int **iring3,**iring4,**iring5,**iring6;
112 } ringdata;
113
114
115 void kcharge()
116 {
117 int i, j, ia, ib, it, iit, kit, jji;
118 long int mask;
119 char pa[4],pb[4],pt[7];
120 float bmomj;
121
122 mask = 1L << 0;
123
124 for (i=1; i <= natom; i++)
125 {
126 atom[i].sigma_charge = 0.0;
127 if (field.type == MMX && (atom[i].type < 300))
128 atom[i].charge = 0.0;
129 else if (field.type != MMX && atom[i].type < 300)
130 atom[i].charge = 0.0;
131 }
132
133 if (field.type == MMFF94)
134 {
135 compute_fcharge();
136 for(i=0; i < bonds_ff.nbnd; i++)
137 {
138 ia = bonds_ff.i12[i][0];
139 ib = bonds_ff.i12[i][1];
140 iit = atom[ia].type;
141 kit = atom[ib].type;
142 if (iit < kit)
143 it = iit*100+kit;
144 else
145 it = kit*100+iit;
146 if (bonds_ff.index[i] == 1)
147 {
148 for(j=0; j < charge_k.nbndchrgdel; j++)
149 {
150 if (it == charge_k.btypedel[j])
151 {
152 if (iit < kit)
153 {
154 atom[ia].charge -= charge_k.bchargedel[j];
155 atom[ib].charge += charge_k.bchargedel[j];
156 }else
157 {
158 atom[ib].charge -= charge_k.bchargedel[j];
159 atom[ia].charge += charge_k.bchargedel[j];
160 }
161 break;
162 }
163 }
164 } else
165 {
166 for(j=0; j < charge_k.nbndchrg; j++)
167 {
168 if (it == charge_k.btype[j])
169 {
170 if (iit < kit)
171 {
172 atom[ia].charge -= charge_k.bcharge[j];
173 atom[ib].charge += charge_k.bcharge[j];
174 }else
175 {
176 atom[ib].charge -= charge_k.bcharge[j];
177 atom[ia].charge += charge_k.bcharge[j];
178 }
179 break;
180 }
181 }
182 }
183 }
184 // add in formal charges
185 for (i=1; i <= natom; i++)
186 {
187 if (atom[i].atomnum == 15)
188 atom[i].formal_charge = 0.0;
189 if (atom[i].atomnum == 16 && atom[i].type != 72)
190 atom[i].formal_charge = 0.0;
191 }
192 for (i=1; i <= natom; i++)
193 {
194 atom[i].charge += (1.0 - atom_k.ligands[atom[i].type]*charge_k.formchrg[atom[i].type])*atom[i].formal_charge ;
195 for(j=0; j < MAXIAT; j++)
196 {
197 if (atom[i].iat[j] != 0)
198 {
199 if (atom[atom[i].iat[j]].formal_charge < 0.0)
200 atom[i].charge += charge_k.formchrg[atom[atom[i].iat[j]].type]*atom[atom[i].iat[j]].formal_charge;
201 }
202 }
203 }
204 }
205 }
206 /* -------------------------------------------------------- */
207 void compute_fcharge()
208 {
209 int i, j, jjbo;
210 int jatm, k, jji, attached[4];
211 int *used;
212 float charge;
213
214 used = ivector(1,natom+1);
215 for (i=1; i <= natom; i++)
216 used[i] = FALSE;
217
218 for (i=1; i <= natom; i++)
219 {
220 atom[i].formal_charge = charge_k.typechrg[atom[i].type];
221 }
222 // find variable charges: types 32, 72, 76, 81
223 for (i=1; i <= natom; i++)
224 {
225 if (used[i] == FALSE)
226 {
227 if (atom[i].formal_charge != 0.0)
228 {
229 if (atom[i].type == 32)
230 {
231 jatm = atom[i].iat[0];
232 jji = 0;
233 jjbo = 0;
234 for (k=0; k < MAXIAT; k++)
235 {
236 if (atom[jatm].iat[k] != 0 && atom[jatm].bo[k] != 9)
237 {
238 if (atom[atom[jatm].iat[k]].type == 32)
239 {
240 if (atom[jatm].bo[k] == 1)
241 {
242 jjbo++;
243 attached[jji] = atom[jatm].iat[k];
244 jji++;
245 } else
246 {
247 attached[jji] = atom[jatm].iat[k];
248 jji++;
249 }
250 }
251 }
252 }
253 charge = -(double)jjbo/jji;
254 if (jji == 1) // only 1 type 32
255 atom[i].formal_charge = 0.0;
256 else if (jji == 2)
257 {
258 if (atom[jatm].atomnum == 15) // p-o treat as single negative charge;
259 charge = -1.00/jji;
260 if (atom[jatm].type == 73)
261 {
262 charge = -0.500;
263 atom[attached[0]].formal_charge = charge;
264 atom[attached[1]].formal_charge = charge;
265 used[attached[0]] = TRUE;
266 used[attached[1]] = TRUE;
267 } else if (atom[jatm].atomnum != 7)
268 {
269 atom[attached[0]].formal_charge = charge;
270 atom[attached[1]].formal_charge = charge;
271 used[attached[0]] = TRUE;
272 used[attached[1]] = TRUE;
273 }else
274 {
275 atom[attached[0]].formal_charge = 0.0;
276 atom[attached[1]].formal_charge = 0.0;
277 used[attached[0]] = TRUE;
278 used[attached[1]] = TRUE;
279 }
280 } else if (jji == 3)
281 {
282 atom[attached[0]].formal_charge = charge;
283 atom[attached[1]].formal_charge = charge;
284 atom[attached[2]].formal_charge = charge;
285 used[attached[0]] = TRUE;
286 used[attached[1]] = TRUE;
287 used[attached[2]] = TRUE;
288 } else if (jji == 4)
289 {
290 atom[attached[0]].formal_charge = charge;
291 atom[attached[1]].formal_charge = charge;
292 atom[attached[2]].formal_charge = charge;
293 atom[attached[3]].formal_charge = charge;
294 used[attached[0]] = TRUE;
295 used[attached[1]] = TRUE;
296 used[attached[2]] = TRUE;
297 used[attached[3]] = TRUE;
298 }
299 }else if (atom[i].type == 72)
300 {
301 jatm = atom[i].iat[0];
302 jji = 0;
303 jjbo = 0;
304 for (k=0; k < MAXIAT; k++)
305 {
306 if (atom[jatm].iat[k] != 0 && atom[jatm].bo[k] != 9)
307 {
308 if (atom[atom[jatm].iat[k]].type == 72)
309 {
310 if (atom[jatm].bo[k] == 1)
311 {
312 jjbo++;
313 attached[jji] = atom[jatm].iat[k];
314 jji++;
315 } else
316 {
317 attached[jji] = atom[jatm].iat[k];
318 jji++;
319 }
320 }
321 }
322 }
323 if (jji == 1)
324 {
325 if (atom[jatm].atomnum == 15)
326 atom[i].formal_charge = 0.0;
327 else
328 atom[i].formal_charge = -1.00;
329 }else
330 {
331 atom[attached[0]].formal_charge = -0.5;
332 atom[attached[1]].formal_charge = -0.5;
333 used[attached[0]] = TRUE;
334 used[attached[1]] = TRUE;
335 }
336 } else if (atom[i].type == 81)
337 {
338 jji = 0;
339 for (k=0; k < MAXIAT ; k++)
340 {
341 if (atom[i].iat[k] != 0 && atom[atom[i].iat[k]].type == 80)
342 {
343 jatm = atom[i].iat[k];
344 jji = 0;
345 for (j=0; j < MAXIAT; j++)
346 {
347 if (atom[jatm].iat[j] != 0 && atom[jatm].bo[j] != 9)
348 {
349 if (atom[atom[jatm].iat[j]].atomnum == 7)
350 {
351 attached[jji] = atom[jatm].iat[j];
352 jji++;
353 }
354 }
355 }
356 }
357 }
358 if (jji == 0)
359 {
360 charge = 1.00;
361 atom[i].formal_charge = charge;
362 used[i] = TRUE;
363 } else
364 {
365 charge = 1.00/jji;
366 for (j=0; j < jji; j++)
367 {
368 atom[attached[j]].formal_charge = charge;
369 used[attached[j]] = TRUE;
370 }
371 }
372 }
373 }
374 }
375 }
376 free_ivector(used,1, natom+1);
377 }