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 |
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 |
}
|