ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/ebufchrg.c
Revision: 99
Committed: Mon Jan 19 04:31:29 2009 UTC (10 years, 8 months ago) by wdelano
File size: 9245 byte(s)
Log Message:
synchronized with trunk, less openmp lib/include
Line File contents
1 #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 fprintf(pcmlogfile,"\nCharge-Charge Interactions\n");
48 fprintf(pcmlogfile," At1 At2 q1 q2 Rik Eqq\n");
49 }
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 if ( skip[i][j] != i)
61 {
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 if (skip[i][j] == -i)
97 e *= units.chgscale;
98 energies.eu += e;
99 atom[i].energy += e;
100 atom[j].energy += e;
101 if (minim_values.iprint)
102 fprintf(pcmlogfile,"QQ: %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.3f = %8.4f\n",
103 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 if (atom[j].charge != 0.0 && skip[i][j] != i )
151 {
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 if (skip[i][j] == -i)
185 st *= units.chgscale;
186 e = cc*rdielc*atom[i].charge*atom[j].charge/(rik+0.05);
187 if (skip[i][j] == -i)
188 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
198 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 if (atom[k].charge != 0.0 && (skip[i][k] != i) )
247 {
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 if (skip[i][k] == -i)
269 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 }