ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/ebufchrg.c
Revision: 93
Committed: Sat Jan 17 23:24:55 2009 UTC (13 years, 5 months ago) by wdelano
File size: 9245 byte(s)
Log Message:
branch created
Line User Rev File contents
1 tjod 3 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5    
6     #include "energies.h"
7     #include "angles.h"
8     #include "derivs.h"
9     #include "hess.h"
10     #include "utility.h"
11     #include "field.h"
12     #include "attached.h"
13     #include "cutoffs.h"
14     #include "pot.h"
15     #include "solv.h"
16    
17     EXTERN struct t_minim_values {
18     int iprint, ndc, nconst;
19     float dielc;
20     } minim_values;
21    
22    
23     void ebufcharge()
24     {
25     int i, j, k,temp;
26     double xi, yi, zi, xk,yk,zk;
27     double xr, yr, zr;
28     double e, rik, cc, rik2;
29     double rdielc;
30     double rdn;
31     double cutoff;
32    
33     cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
34    
35     energies.eu = 0.0;
36     cc = 332.0716;
37     rdielc = 1.0/units.dielec;
38    
39     rdn = 1.0;
40     if (field.type == MMX || field.type == MM2 )
41     rdn = .915;
42     else if (field.type == MM3)
43     rdn = .923;
44    
45     if (minim_values.iprint)
46     {
47 wdelano 58 fprintf(pcmlogfile,"\nCharge-Charge Interactions\n");
48     fprintf(pcmlogfile," At1 At2 q1 q2 Rik Eqq\n");
49 tjod 3 }
50    
51     for (i=1; i < natom; i++)
52     {
53     if (atom[i].charge != 0.0 )
54     {
55    
56     for(j=i+1; j <= natom; j++)
57     {
58     if (atom[i].use || atom[j].use)
59     {
60 gilbertke 89 if ( skip[i][j] != i)
61 tjod 3 {
62     if (atom[j].charge != 0.0)
63     {
64     if (atom[i].type == 5 || atom[i].type == 36)
65     {
66     // compute reduced distance
67     xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
68     yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
69     zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
70     } else
71     {
72     xi = atom[i].x;
73     yi = atom[i].y;
74     zi = atom[i].z;
75     }
76     if (atom[j].type == 5 || atom[j].type == 36)
77     {
78     // compute reduced distance
79     xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
80     yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
81     zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
82     } else
83     {
84     xk = atom[j].x;
85     yk = atom[j].y;
86     zk = atom[j].z;
87     }
88     xr = xi - xk;
89     yr = yi - yk;
90     zr = zi - zk;
91     rik2 = xr*xr + yr*yr + zr*zr;
92     if (rik2 < cutoff)
93     {
94     rik = sqrt( xr*xr + yr*yr + zr*zr);
95     e = cc*atom[i].charge*atom[j].charge/(units.dielec*(rik+0.05));
96 gilbertke 89 if (skip[i][j] == -i)
97 tjod 3 e *= units.chgscale;
98     energies.eu += e;
99     atom[i].energy += e;
100     atom[j].energy += e;
101     if (minim_values.iprint)
102 wdelano 58 fprintf(pcmlogfile,"QQ: %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.3f = %8.4f\n",
103 tjod 3 atom[i].name, i, atom[j].name, j,atom[i].charge, atom[j].charge,rik,e);
104     }
105     }
106     }
107     }
108     }
109     }
110     }
111     }
112    
113     void ebufcharge1()
114     {
115     int i,j, k;
116     double xi, yi, zi;
117     double xr, yr, zr;
118     double xk,yk,zk,rdn;
119     double e, rik, cc, rik2, st;
120     double rdielc;
121     double cutoff;
122    
123     cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
124    
125     energies.eu = 0.0;
126     cc = 332.0716;
127     rdielc = 1.0/units.dielec;
128     rdn = 1.0;
129     if (field.type == MMX || field.type == MM2 )
130     rdn = .915;
131     else if (field.type == MM3)
132     rdn = .923;
133    
134     for (i=0; i <= natom; i++)
135     {
136     deriv.deqq[i][0] = 0.0;
137     deriv.deqq[i][1] = 0.0;
138     deriv.deqq[i][2] = 0.0;
139     }
140    
141    
142     for (i=1; i < natom; i++)
143     {
144     if (atom[i].charge != 0.0 )
145     {
146     for(j=i+1; j <= natom; j++)
147     {
148     if (atom[i].use || atom[j].use)
149     {
150 gilbertke 89 if (atom[j].charge != 0.0 && skip[i][j] != i )
151 tjod 3 {
152     if (atom[i].type == 5 || atom[i].type == 36)
153     {
154     // compute reduced distance
155     xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
156     yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
157     zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
158     } else
159     {
160     xi = atom[i].x;
161     yi = atom[i].y;
162     zi = atom[i].z;
163     }
164     if (atom[j].type == 5 || atom[j].type == 36)
165     {
166     // compute reduced distance
167     xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
168     yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
169     zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
170     } else
171     {
172     xk = atom[j].x;
173     yk = atom[j].y;
174     zk = atom[j].z;
175     }
176     xr = xi - xk;
177     yr = yi - yk;
178     zr = zi - zk;
179     rik2 = xr*xr +yr*yr + zr*zr;
180     if (rik2 < cutoff)
181     {
182     rik = sqrt( rik2);
183     st = -cc*rdielc*atom[i].charge*atom[j].charge/(rik*(rik+0.05)*(rik+0.05));
184 gilbertke 89 if (skip[i][j] == -i)
185 tjod 3 st *= units.chgscale;
186     e = cc*rdielc*atom[i].charge*atom[j].charge/(rik+0.05);
187 gilbertke 89 if (skip[i][j] == -i)
188 tjod 3 e *= units.chgscale;
189    
190     deriv.deqq[i][0] += st*xr;
191     deriv.deqq[i][1] += st*yr;
192     deriv.deqq[i][2] += st*zr;
193    
194     deriv.deqq[j][0] -= st*xr;
195     deriv.deqq[j][1] -= st*yr;
196     deriv.deqq[j][2] -= st*zr;
197 gilbertke 89
198 tjod 3 energies.eu += e;
199     }
200     }
201     }
202     }
203     }
204     }
205     }
206    
207     void ebufcharge2(int i)
208     {
209     int j,k, jj,temp;
210     double fik,de,d2e,d2edx,d2edy,d2edz;
211     double xi,yi,zi,xr,yr,zr,term[3][3];
212     double xk,yk,zk,rdn;
213     double r,r2;
214     double cc, rdielc;
215     double cutoff;
216    
217     if (atom[i].charge == 0.0)
218     return;
219    
220     cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
221    
222     cc = 332.0716;
223     rdielc = 1.0/units.dielec;
224     rdn = 1.0;
225     if (field.type == MMX || field.type == MM2 )
226     rdn = .915;
227     else if (field.type == MM3)
228     rdn = .923;
229    
230    
231     if (atom[i].type == 5 || atom[i].type == 36)
232     {
233     // compute reduced distance
234     xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
235     yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
236     zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
237     } else
238     {
239     xi = atom[i].x;
240     yi = atom[i].y;
241     zi = atom[i].z;
242     }
243    
244     for (k=1; k <= natom; k++)
245     {
246 gilbertke 89 if (atom[k].charge != 0.0 && (skip[i][k] != i) )
247 tjod 3 {
248     if (atom[j].type == 5 || atom[j].type == 36)
249     {
250     // compute reduced distance
251     xk = rdn*(atom[k].x-atom[atom[k].iat[0]].x) + atom[atom[k].iat[0]].x;
252     yk = rdn*(atom[k].y-atom[atom[k].iat[0]].y) + atom[atom[k].iat[0]].y;
253     zk = rdn*(atom[k].z-atom[atom[k].iat[0]].z) + atom[atom[k].iat[0]].z;
254     } else
255     {
256     xk = atom[k].x;
257     yk = atom[k].y;
258     zk = atom[k].z;
259     }
260     xr = xi - xk;
261     yr = yi - yk;
262     zr = zi - zk;
263     r2 = xr*xr + yr*yr + zr*zr;
264     if (r2 < cutoff)
265     {
266     r = sqrt(r2);
267     fik = cc*rdielc*atom[i].charge*atom[k].charge;
268 gilbertke 89 if (skip[i][k] == -i)
269 tjod 3 fik *= units.chgscale;
270     de = -fik/((r+0.05)*(r+0.05));
271     d2e = -2.0*de/r;
272     de = de / r;
273     d2e = (d2e-de) / r2;
274     d2edx = d2e * xr;
275     d2edy = d2e * yr;
276     d2edz = d2e * zr;
277     term[0][0]= d2edx*xr + de;
278     term[1][0] = d2edx*yr;
279     term[2][0] = d2edx*zr;
280     term[0][1] = term[1][0];
281     term[1][1] = d2edy*yr + de;
282     term[2][1] = d2edy*zr;
283     term[0][2] = term[2][0];
284     term[1][2] = term[2][1];
285     term[2][2] = d2edz*zr + de;
286    
287     for (j=0; j < 3; j++)
288     {
289     hess.hessx[i][j] += term[j][0];
290     hess.hessy[i][j] += term[j][1];
291     hess.hessz[i][j] += term[j][2];
292     hess.hessx[k][j] -= term[j][0];
293     hess.hessy[k][j] -= term[j][1];
294     hess.hessz[k][j] -= term[j][2];
295     }
296     }
297     }
298     }
299     }