ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/freemol/branches/sync4pymol12/src/mengine/src/kcharge.c

#	Line 1 | Line 1
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	<	}
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 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	>	}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines