ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/ebufchrg.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 10908 byte(s)
Log Message:
test

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_units {
18 double bndunit, cbnd, qbnd;
19 double angunit, cang, qang, pang, sang, aaunit;
20 double stbnunit, ureyunit, torsunit, storunit, v14scale;
21 double aterm, bterm, cterm, dielec, chgscale;
22 } units;
23
24 EXTERN struct t_minim_values {
25 int iprint, ndc, nconst;
26 float dielc;
27 } minim_values;
28
29 EXTERN int *skip;
30
31 void image(double *, double *, double *, int mode);
32
33 void ebufcharge()
34 {
35 int i, j, k,temp;
36 double xi, yi, zi, xk,yk,zk;
37 double xr, yr, zr;
38 double e, rik, cc, rik2;
39 double rdielc;
40 double rdn;
41 double cutoff;
42
43 cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
44
45 energies.eu = 0.0;
46 cc = 332.0716;
47 rdielc = 1.0/units.dielec;
48
49 rdn = 1.0;
50 if (field.type == MMX || field.type == MM2 )
51 rdn = .915;
52 else if (field.type == MM3)
53 rdn = .923;
54
55 if (minim_values.iprint)
56 {
57 fprintf(pcmoutfile,"\nCharge-Charge Interactions\n");
58 fprintf(pcmoutfile," At1 At2 q1 q2 Rik Eqq\n");
59 }
60 for (i=1; i <= natom; i++)
61 skip[i] = 0;
62
63 for (i=1; i < natom; i++)
64 {
65 if (atom[i].charge != 0.0 )
66 {
67 for (k=0; k < MAXIAT; k++)
68 {
69 if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
70 skip[atom[i].iat[k]] = i;
71 }
72 for(k=0; k < attached.n13[i]; k++)
73 skip[attached.i13[k][i]] = i;
74 for(k=0; k < attached.n14[i]; k++)
75 {
76 temp = attached.i14[k][i];
77 skip[attached.i14[k][i]] = -i;
78 }
79
80 for(j=i+1; j <= natom; j++)
81 {
82 if (atom[i].use || atom[j].use)
83 {
84 if ( skip[j] != i)
85 {
86 if (atom[j].charge != 0.0)
87 {
88 if (atom[i].type == 5 || atom[i].type == 36)
89 {
90 // compute reduced distance
91 xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
92 yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
93 zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
94 } else
95 {
96 xi = atom[i].x;
97 yi = atom[i].y;
98 zi = atom[i].z;
99 }
100 if (atom[j].type == 5 || atom[j].type == 36)
101 {
102 // compute reduced distance
103 xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
104 yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
105 zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
106 } else
107 {
108 xk = atom[j].x;
109 yk = atom[j].y;
110 zk = atom[j].z;
111 }
112 xr = xi - xk;
113 yr = yi - yk;
114 zr = zi - zk;
115 rik2 = xr*xr + yr*yr + zr*zr;
116 if (rik2 < cutoff)
117 {
118 rik = sqrt( xr*xr + yr*yr + zr*zr);
119 e = cc*atom[i].charge*atom[j].charge/(units.dielec*(rik+0.05));
120 if (skip[j] == -i)
121 e *= units.chgscale;
122 energies.eu += e;
123 atom[i].energy += e;
124 atom[j].energy += e;
125 if (minim_values.iprint)
126 fprintf(pcmoutfile,"QQ: %2s(%-3d)- %2s(%-3d) %-8.3f %-8.3f %-8.3f = %8.4f\n",
127 atom[i].name, i, atom[j].name, j,atom[i].charge, atom[j].charge,rik,e);
128 }
129 }
130 }
131 }
132 }
133 }
134 }
135 }
136
137 void ebufcharge1()
138 {
139 int i,j, k;
140 double xi, yi, zi;
141 double xr, yr, zr;
142 double xk,yk,zk,rdn;
143 double e, rik, cc, rik2, st;
144 double rdielc;
145 double cutoff;
146
147 cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
148
149 energies.eu = 0.0;
150 cc = 332.0716;
151 rdielc = 1.0/units.dielec;
152 rdn = 1.0;
153 if (field.type == MMX || field.type == MM2 )
154 rdn = .915;
155 else if (field.type == MM3)
156 rdn = .923;
157
158 for (i=0; i <= natom; i++)
159 {
160 deriv.deqq[i][0] = 0.0;
161 deriv.deqq[i][1] = 0.0;
162 deriv.deqq[i][2] = 0.0;
163 }
164
165 for (i=1; i <= natom; i++)
166 skip[i] = 0;
167
168 for (i=1; i < natom; i++)
169 {
170 if (atom[i].charge != 0.0 )
171 {
172 for (k=0; k < MAXIAT; k++)
173 {
174 if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
175 skip[atom[i].iat[k]] = i;
176 }
177 for(k=0; k < attached.n13[i]; k++)
178 skip[attached.i13[k][i]] = i;
179 for(k=0; k < attached.n14[i]; k++)
180 skip[attached.i14[k][i]] = -i;
181
182 for(j=i+1; j <= natom; j++)
183 {
184 if (atom[i].use || atom[j].use)
185 {
186 if (atom[j].charge != 0.0 && skip[j] != i )
187 {
188 if (atom[i].type == 5 || atom[i].type == 36)
189 {
190 // compute reduced distance
191 xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
192 yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
193 zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
194 } else
195 {
196 xi = atom[i].x;
197 yi = atom[i].y;
198 zi = atom[i].z;
199 }
200 if (atom[j].type == 5 || atom[j].type == 36)
201 {
202 // compute reduced distance
203 xk = rdn*(atom[j].x-atom[atom[j].iat[0]].x) + atom[atom[j].iat[0]].x;
204 yk = rdn*(atom[j].y-atom[atom[j].iat[0]].y) + atom[atom[j].iat[0]].y;
205 zk = rdn*(atom[j].z-atom[atom[j].iat[0]].z) + atom[atom[j].iat[0]].z;
206 } else
207 {
208 xk = atom[j].x;
209 yk = atom[j].y;
210 zk = atom[j].z;
211 }
212 xr = xi - xk;
213 yr = yi - yk;
214 zr = zi - zk;
215 rik2 = xr*xr +yr*yr + zr*zr;
216 if (rik2 < cutoff)
217 {
218 rik = sqrt( rik2);
219 st = -cc*rdielc*atom[i].charge*atom[j].charge/(rik*(rik+0.05)*(rik+0.05));
220 if (skip[j] == -i)
221 st *= units.chgscale;
222 e = cc*rdielc*atom[i].charge*atom[j].charge/(rik+0.05);
223 if (skip[j] == -i)
224 e *= units.chgscale;
225
226 deriv.deqq[i][0] += st*xr;
227 deriv.deqq[i][1] += st*yr;
228 deriv.deqq[i][2] += st*zr;
229
230 deriv.deqq[j][0] -= st*xr;
231 deriv.deqq[j][1] -= st*yr;
232 deriv.deqq[j][2] -= st*zr;
233
234 virial.virx += xr*st*xr;
235 virial.viry += yr*st*yr;
236 virial.virz += zr*st*zr;
237
238 energies.eu += e;
239 }
240 }
241 }
242 }
243 }
244 }
245 }
246
247 void ebufcharge2(int i)
248 {
249 int j,k, jj,temp;
250 double fik,de,d2e,d2edx,d2edy,d2edz;
251 double xi,yi,zi,xr,yr,zr,term[3][3];
252 double xk,yk,zk,rdn;
253 double r,r2;
254 double cc, rdielc;
255 double cutoff;
256
257 if (atom[i].charge == 0.0)
258 return;
259
260 cutoff = cutoffs.chrgcut*cutoffs.chrgcut;
261
262 cc = 332.0716;
263 rdielc = 1.0/units.dielec;
264 rdn = 1.0;
265 if (field.type == MMX || field.type == MM2 )
266 rdn = .915;
267 else if (field.type == MM3)
268 rdn = .923;
269
270 for (j=1; j <= natom; j++)
271 skip[j] = 0;
272
273 skip[i] = i;
274 for (k=0; k < MAXIAT; k++)
275 {
276 if (atom[i].iat[k] != 0 && atom[i].bo[k] != 9)
277 skip[atom[i].iat[k]] = i;
278 }
279 for(k=0; k < attached.n13[i]; k++)
280 skip[attached.i13[k][i]] = i;
281 for(k=0; k < attached.n14[i]; k++)
282 {
283 temp = attached.i14[k][i];
284 skip[attached.i14[k][i]] = -i;
285 }
286
287 if (atom[i].type == 5 || atom[i].type == 36)
288 {
289 // compute reduced distance
290 xi = rdn*(atom[i].x-atom[atom[i].iat[0]].x) + atom[atom[i].iat[0]].x;
291 yi = rdn*(atom[i].y-atom[atom[i].iat[0]].y) + atom[atom[i].iat[0]].y;
292 zi = rdn*(atom[i].z-atom[atom[i].iat[0]].z) + atom[atom[i].iat[0]].z;
293 } else
294 {
295 xi = atom[i].x;
296 yi = atom[i].y;
297 zi = atom[i].z;
298 }
299
300 for (k=1; k <= natom; k++)
301 {
302 if (atom[k].charge != 0.0 && (skip[k] != i) )
303 {
304 if (atom[j].type == 5 || atom[j].type == 36)
305 {
306 // compute reduced distance
307 xk = rdn*(atom[k].x-atom[atom[k].iat[0]].x) + atom[atom[k].iat[0]].x;
308 yk = rdn*(atom[k].y-atom[atom[k].iat[0]].y) + atom[atom[k].iat[0]].y;
309 zk = rdn*(atom[k].z-atom[atom[k].iat[0]].z) + atom[atom[k].iat[0]].z;
310 } else
311 {
312 xk = atom[k].x;
313 yk = atom[k].y;
314 zk = atom[k].z;
315 }
316 xr = xi - xk;
317 yr = yi - yk;
318 zr = zi - zk;
319 r2 = xr*xr + yr*yr + zr*zr;
320 if (r2 < cutoff)
321 {
322 r = sqrt(r2);
323 fik = cc*rdielc*atom[i].charge*atom[k].charge;
324 if (skip[k] == -i)
325 fik *= units.chgscale;
326 de = -fik/((r+0.05)*(r+0.05));
327 d2e = -2.0*de/r;
328 de = de / r;
329 d2e = (d2e-de) / r2;
330 d2edx = d2e * xr;
331 d2edy = d2e * yr;
332 d2edz = d2e * zr;
333 term[0][0]= d2edx*xr + de;
334 term[1][0] = d2edx*yr;
335 term[2][0] = d2edx*zr;
336 term[0][1] = term[1][0];
337 term[1][1] = d2edy*yr + de;
338 term[2][1] = d2edy*zr;
339 term[0][2] = term[2][0];
340 term[1][2] = term[2][1];
341 term[2][2] = d2edz*zr + de;
342
343 for (j=0; j < 3; j++)
344 {
345 hess.hessx[i][j] += term[j][0];
346 hess.hessy[i][j] += term[j][1];
347 hess.hessz[i][j] += term[j][2];
348 hess.hessx[k][j] -= term[j][0];
349 hess.hessy[k][j] -= term[j][1];
350 hess.hessz[k][j] -= term[j][2];
351 }
352 }
353 }
354 }
355 }