ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/kcharge.c
(Generate patch)
# 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