File: | gromacs/gmxlib/ewald_util.c |
Location: | line 274, column 29 |
Description: | Value stored to 'fscal' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2004, The GROMACS development team. |
6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
8 | * and including many others, as listed in the AUTHORS file in the |
9 | * top-level source directory and at http://www.gromacs.org. |
10 | * |
11 | * GROMACS is free software; you can redistribute it and/or |
12 | * modify it under the terms of the GNU Lesser General Public License |
13 | * as published by the Free Software Foundation; either version 2.1 |
14 | * of the License, or (at your option) any later version. |
15 | * |
16 | * GROMACS is distributed in the hope that it will be useful, |
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
19 | * Lesser General Public License for more details. |
20 | * |
21 | * You should have received a copy of the GNU Lesser General Public |
22 | * License along with GROMACS; if not, see |
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
25 | * |
26 | * If you want to redistribute modifications to GROMACS, please |
27 | * consider that scientific software is very special. Version |
28 | * control is crucial - bugs must be traceable. We will be happy to |
29 | * consider code for inclusion in the official distribution, but |
30 | * derived work must not be called official GROMACS. Details are found |
31 | * in the README & COPYING files - if they are missing, get the |
32 | * official version at http://www.gromacs.org. |
33 | * |
34 | * To help us fund GROMACS development, we humbly ask that you cite |
35 | * the research papers on the package. Check out http://www.gromacs.org. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <stdio.h> |
42 | #include <math.h> |
43 | #include "gromacs/math/utilities.h" |
44 | #include "typedefs.h" |
45 | #include "types/commrec.h" |
46 | #include "gromacs/math/vec.h" |
47 | #include "coulomb.h" |
48 | #include "gromacs/utility/smalloc.h" |
49 | #include "physics.h" |
50 | #include "txtdump.h" |
51 | #include "gromacs/utility/futil.h" |
52 | #include "names.h" |
53 | #include "macros.h" |
54 | |
55 | real calc_ewaldcoeff_q(real rc, real dtol) |
56 | { |
57 | real x = 5, low, high; |
58 | int n, i = 0; |
59 | |
60 | |
61 | do |
62 | { |
63 | i++; |
64 | x *= 2; |
65 | } |
66 | while (gmx_erfc(x*rc)gmx_erfcf(x*rc) > dtol); |
67 | |
68 | n = i+60; /* search tolerance is 2^-60 */ |
69 | low = 0; |
70 | high = x; |
71 | for (i = 0; i < n; i++) |
72 | { |
73 | x = (low+high)/2; |
74 | if (gmx_erfc(x*rc)gmx_erfcf(x*rc) > dtol) |
75 | { |
76 | low = x; |
77 | } |
78 | else |
79 | { |
80 | high = x; |
81 | } |
82 | } |
83 | return x; |
84 | } |
85 | |
86 | static real ewald_function_lj(real x, real rc) |
87 | { |
88 | real xrc, xrc2, xrc4, factor; |
89 | xrc = x*rc; |
90 | xrc2 = xrc*xrc; |
91 | xrc4 = xrc2*xrc2; |
92 | #ifdef GMX_DOUBLE |
93 | factor = exp(-xrc2)*(1 + xrc2 + xrc4/2.0); |
94 | #else |
95 | factor = expf(-xrc2)*(1 + xrc2 + xrc4/2.0); |
96 | #endif |
97 | |
98 | return factor; |
99 | } |
100 | |
101 | real calc_ewaldcoeff_lj(real rc, real dtol) |
102 | { |
103 | real x = 5, low, high; |
104 | int n, i = 0; |
105 | |
106 | do |
107 | { |
108 | i++; |
109 | x *= 2.0; |
110 | } |
111 | while (ewald_function_lj(x, rc) > dtol); |
112 | |
113 | n = i + 60; /* search tolerance is 2^-60 */ |
114 | low = 0; |
115 | high = x; |
116 | for (i = 0; i < n; ++i) |
117 | { |
118 | x = (low + high) / 2.0; |
119 | if (ewald_function_lj(x, rc) > dtol) |
120 | { |
121 | low = x; |
122 | } |
123 | else |
124 | { |
125 | high = x; |
126 | } |
127 | } |
128 | return x; |
129 | } |
130 | |
131 | void ewald_LRcorrection(int start, int end, |
132 | t_commrec *cr, int thread, t_forcerec *fr, |
133 | real *chargeA, real *chargeB, |
134 | real *C6A, real *C6B, |
135 | real *sigmaA, real *sigmaB, |
136 | real *sigma3A, real *sigma3B, |
137 | gmx_bool calc_excl_corr, |
138 | t_blocka *excl, rvec x[], |
139 | matrix box, rvec mu_tot[], |
140 | int ewald_geometry, real epsilon_surface, |
141 | rvec *f, tensor vir_q, tensor vir_lj, |
142 | real *Vcorr_q, real *Vcorr_lj, |
143 | real lambda_q, real lambda_lj, |
144 | real *dvdlambda_q, real *dvdlambda_lj) |
145 | { |
146 | int i, i1, i2, j, k, m, iv, jv, q; |
147 | atom_id *AA; |
148 | double Vexcl_q, dvdl_excl_q, dvdl_excl_lj; /* Necessary for precision */ |
149 | double Vexcl_lj; |
150 | real one_4pi_eps; |
151 | real v, vc, qiA, qiB, dr2, rinv, enercorr; |
152 | real Vself_q[2], Vself_lj[2], Vdipole[2], rinv2, ewc_q = fr->ewaldcoeff_q, ewcdr; |
153 | real ewc_lj = fr->ewaldcoeff_lj, ewc_lj2 = ewc_lj * ewc_lj; |
154 | real c6Ai = 0, c6Bi = 0, c6A = 0, c6B = 0, ewcdr2, ewcdr4, c6L = 0, rinv6; |
155 | rvec df, dx, mutot[2], dipcorrA, dipcorrB; |
156 | tensor dxdf_q, dxdf_lj; |
157 | real vol = box[XX0][XX0]*box[YY1][YY1]*box[ZZ2][ZZ2]; |
158 | real L1_q, L1_lj, dipole_coeff, qqA, qqB, qqL, vr0_q, vr0_lj = 0; |
159 | gmx_bool bFreeEnergy = (chargeB != NULL((void*)0)); |
160 | gmx_bool bMolPBC = fr->bMolPBC; |
161 | gmx_bool bDoingLBRule = (fr->ljpme_combination_rule == eljpmeLB); |
162 | |
163 | /* This routine can be made faster by using tables instead of analytical interactions |
164 | * However, that requires a thorough verification that they are correct in all cases. |
165 | */ |
166 | |
167 | one_4pi_eps = ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/fr->epsilon_r; |
168 | vr0_q = ewc_q*M_2_SQRTPI1.12837916709551257390; |
169 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
170 | { |
171 | vr0_lj = -pow(ewc_lj, 6)/6.0; |
172 | } |
173 | |
174 | AA = excl->a; |
175 | Vexcl_q = 0; |
176 | Vexcl_lj = 0; |
177 | dvdl_excl_q = 0; |
178 | dvdl_excl_lj = 0; |
179 | Vdipole[0] = 0; |
180 | Vdipole[1] = 0; |
181 | L1_q = 1.0-lambda_q; |
182 | L1_lj = 1.0-lambda_lj; |
183 | /* Note that we have to transform back to gromacs units, since |
184 | * mu_tot contains the dipole in debye units (for output). |
185 | */ |
186 | for (i = 0; (i < DIM3); i++) |
187 | { |
188 | mutot[0][i] = mu_tot[0][i]*DEBYE2ENM0.02081941; |
189 | mutot[1][i] = mu_tot[1][i]*DEBYE2ENM0.02081941; |
190 | dipcorrA[i] = 0; |
191 | dipcorrB[i] = 0; |
192 | } |
193 | dipole_coeff = 0; |
194 | switch (ewald_geometry) |
195 | { |
196 | case eewg3D: |
197 | if (epsilon_surface != 0) |
198 | { |
199 | dipole_coeff = |
200 | 2*M_PI3.14159265358979323846*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/((2*epsilon_surface + fr->epsilon_r)*vol); |
201 | for (i = 0; (i < DIM3); i++) |
202 | { |
203 | dipcorrA[i] = 2*dipole_coeff*mutot[0][i]; |
204 | dipcorrB[i] = 2*dipole_coeff*mutot[1][i]; |
205 | } |
206 | } |
207 | break; |
208 | case eewg3DC: |
209 | dipole_coeff = 2*M_PI3.14159265358979323846*one_4pi_eps/vol; |
210 | dipcorrA[ZZ2] = 2*dipole_coeff*mutot[0][ZZ2]; |
211 | dipcorrB[ZZ2] = 2*dipole_coeff*mutot[1][ZZ2]; |
212 | break; |
213 | default: |
214 | gmx_incons("Unsupported Ewald geometry")_gmx_error("incons", "Unsupported Ewald geometry", "/home/alexxy/Develop/gromacs/src/gromacs/gmxlib/ewald_util.c" , 214); |
215 | break; |
216 | } |
217 | if (debug) |
218 | { |
219 | fprintf(debug, "dipcorr = %8.3f %8.3f %8.3f\n", |
220 | dipcorrA[XX0], dipcorrA[YY1], dipcorrA[ZZ2]); |
221 | fprintf(debug, "mutot = %8.3f %8.3f %8.3f\n", |
222 | mutot[0][XX0], mutot[0][YY1], mutot[0][ZZ2]); |
223 | } |
224 | clear_mat(dxdf_q); |
225 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
226 | { |
227 | clear_mat(dxdf_lj); |
228 | } |
229 | if ((calc_excl_corr || dipole_coeff != 0 || EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) && !bFreeEnergy) |
230 | { |
231 | for (i = start; (i < end); i++) |
232 | { |
233 | /* Initiate local variables (for this i-particle) to 0 */ |
234 | qiA = chargeA[i]*one_4pi_eps; |
235 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
236 | { |
237 | c6Ai = C6A[i]; |
238 | if (bDoingLBRule) |
239 | { |
240 | c6Ai *= sigma3A[i]; |
241 | } |
242 | } |
243 | if (calc_excl_corr) |
244 | { |
245 | i1 = excl->index[i]; |
246 | i2 = excl->index[i+1]; |
247 | |
248 | /* Loop over excluded neighbours */ |
249 | for (j = i1; (j < i2); j++) |
250 | { |
251 | k = AA[j]; |
252 | /* |
253 | * First we must test whether k <> i, and then, |
254 | * because the exclusions are all listed twice i->k |
255 | * and k->i we must select just one of the two. As |
256 | * a minor optimization we only compute forces when |
257 | * the charges are non-zero. |
258 | */ |
259 | if (k > i) |
260 | { |
261 | qqA = qiA*chargeA[k]; |
262 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
263 | { |
264 | c6A = c6Ai * C6A[k]; |
265 | if (bDoingLBRule) |
266 | { |
267 | c6A *= pow(0.5*(sigmaA[i]+sigmaA[k]), 6)*sigma3A[k]; |
268 | } |
269 | } |
270 | if (qqA != 0.0 || c6A != 0.0) |
271 | { |
272 | real fscal; |
273 | |
274 | fscal = 0; |
Value stored to 'fscal' is never read | |
275 | rvec_sub(x[i], x[k], dx); |
276 | if (bMolPBC) |
277 | { |
278 | /* Cheap pbc_dx, assume excluded pairs are at short distance. */ |
279 | for (m = DIM3-1; (m >= 0); m--) |
280 | { |
281 | if (dx[m] > 0.5*box[m][m]) |
282 | { |
283 | rvec_dec(dx, box[m]); |
284 | } |
285 | else if (dx[m] < -0.5*box[m][m]) |
286 | { |
287 | rvec_inc(dx, box[m]); |
288 | } |
289 | } |
290 | } |
291 | dr2 = norm2(dx); |
292 | /* Distance between two excluded particles |
293 | * may be zero in the case of shells |
294 | */ |
295 | if (dr2 != 0) |
296 | { |
297 | rinv = gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); |
298 | rinv2 = rinv*rinv; |
299 | if (qqA != 0.0) |
300 | { |
301 | real dr; |
302 | |
303 | dr = 1.0/rinv; |
304 | ewcdr = ewc_q*dr; |
305 | vc = qqA*gmx_erf(ewcdr)gmx_erff(ewcdr)*rinv; |
306 | Vexcl_q += vc; |
307 | #ifdef GMX_DOUBLE |
308 | /* Relative accuracy at R_ERF_R_INACC of 3e-10 */ |
309 | #define R_ERF_R_INACC0.1 0.006 |
310 | #else |
311 | /* Relative accuracy at R_ERF_R_INACC of 2e-5 */ |
312 | #define R_ERF_R_INACC0.1 0.1 |
313 | #endif |
314 | /* fscal is the scalar force pre-multiplied by rinv, |
315 | * to normalise the relative position vector dx */ |
316 | if (ewcdr > R_ERF_R_INACC0.1) |
317 | { |
318 | fscal = rinv2*(vc - qqA*ewc_q*M_2_SQRTPI1.12837916709551257390*exp(-ewcdr*ewcdr)); |
319 | } |
320 | else |
321 | { |
322 | /* Use a fourth order series expansion for small ewcdr */ |
323 | fscal = ewc_q*ewc_q*qqA*vr0_q*(2.0/3.0 - 0.4*ewcdr*ewcdr); |
324 | } |
325 | |
326 | /* The force vector is obtained by multiplication with |
327 | * the relative position vector |
328 | */ |
329 | svmul(fscal, dx, df); |
330 | rvec_inc(f[k], df); |
331 | rvec_dec(f[i], df); |
332 | for (iv = 0; (iv < DIM3); iv++) |
333 | { |
334 | for (jv = 0; (jv < DIM3); jv++) |
335 | { |
336 | dxdf_q[iv][jv] += dx[iv]*df[jv]; |
337 | } |
338 | } |
339 | } |
340 | |
341 | if (c6A != 0.0) |
342 | { |
343 | rinv6 = rinv2*rinv2*rinv2; |
344 | ewcdr2 = ewc_lj2*dr2; |
345 | ewcdr4 = ewcdr2*ewcdr2; |
346 | /* We get the excluded long-range contribution from -C6*(1-g(r)) |
347 | * g(r) is also defined in the manual under LJ-PME |
348 | */ |
349 | vc = -c6A*rinv6*(1.0 - exp(-ewcdr2)*(1 + ewcdr2 + 0.5*ewcdr4)); |
350 | Vexcl_lj += vc; |
351 | /* The force is the derivative of the potential vc. |
352 | * fscal is the scalar force pre-multiplied by rinv, |
353 | * to normalise the relative position vector dx */ |
354 | fscal = 6.0*vc*rinv2 + c6A*rinv6*exp(-ewcdr2)*ewc_lj2*ewcdr4; |
355 | |
356 | /* The force vector is obtained by multiplication with |
357 | * the relative position vector |
358 | */ |
359 | svmul(fscal, dx, df); |
360 | rvec_inc(f[k], df); |
361 | rvec_dec(f[i], df); |
362 | for (iv = 0; (iv < DIM3); iv++) |
363 | { |
364 | for (jv = 0; (jv < DIM3); jv++) |
365 | { |
366 | dxdf_lj[iv][jv] += dx[iv]*df[jv]; |
367 | } |
368 | } |
369 | } |
370 | } |
371 | else |
372 | { |
373 | Vexcl_q += qqA*vr0_q; |
374 | Vexcl_lj += c6A*vr0_lj; |
375 | } |
376 | } |
377 | } |
378 | } |
379 | } |
380 | /* Dipole correction on force */ |
381 | if (dipole_coeff != 0) |
382 | { |
383 | for (j = 0; (j < DIM3); j++) |
384 | { |
385 | f[i][j] -= dipcorrA[j]*chargeA[i]; |
386 | } |
387 | } |
388 | } |
389 | } |
390 | else if (calc_excl_corr || dipole_coeff != 0 || EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
391 | { |
392 | for (i = start; (i < end); i++) |
393 | { |
394 | /* Initiate local variables (for this i-particle) to 0 */ |
395 | qiA = chargeA[i]*one_4pi_eps; |
396 | qiB = chargeB[i]*one_4pi_eps; |
397 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
398 | { |
399 | c6Ai = C6A[i]; |
400 | c6Bi = C6B[i]; |
401 | if (bDoingLBRule) |
402 | { |
403 | c6Ai *= sigma3A[i]; |
404 | c6Bi *= sigma3B[i]; |
405 | } |
406 | } |
407 | if (calc_excl_corr) |
408 | { |
409 | i1 = excl->index[i]; |
410 | i2 = excl->index[i+1]; |
411 | |
412 | /* Loop over excluded neighbours */ |
413 | for (j = i1; (j < i2); j++) |
414 | { |
415 | k = AA[j]; |
416 | if (k > i) |
417 | { |
418 | qqA = qiA*chargeA[k]; |
419 | qqB = qiB*chargeB[k]; |
420 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
421 | { |
422 | c6A = c6Ai*C6A[k]; |
423 | c6B = c6Bi*C6B[k]; |
424 | if (bDoingLBRule) |
425 | { |
426 | c6A *= pow(0.5*(sigmaA[i]+sigmaA[k]), 6)*sigma3A[k]; |
427 | c6B *= pow(0.5*(sigmaB[i]+sigmaB[k]), 6)*sigma3B[k]; |
428 | } |
429 | } |
430 | if (qqA != 0.0 || qqB != 0.0 || c6A != 0.0 || c6B != 0.0) |
431 | { |
432 | real fscal; |
433 | |
434 | fscal = 0; |
435 | qqL = L1_q*qqA + lambda_q*qqB; |
436 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
437 | { |
438 | c6L = L1_lj*c6A + lambda_lj*c6B; |
439 | } |
440 | rvec_sub(x[i], x[k], dx); |
441 | if (bMolPBC) |
442 | { |
443 | /* Cheap pbc_dx, assume excluded pairs are at short distance. */ |
444 | for (m = DIM3-1; (m >= 0); m--) |
445 | { |
446 | if (dx[m] > 0.5*box[m][m]) |
447 | { |
448 | rvec_dec(dx, box[m]); |
449 | } |
450 | else if (dx[m] < -0.5*box[m][m]) |
451 | { |
452 | rvec_inc(dx, box[m]); |
453 | } |
454 | } |
455 | } |
456 | dr2 = norm2(dx); |
457 | if (dr2 != 0) |
458 | { |
459 | rinv = gmx_invsqrt(dr2)gmx_software_invsqrt(dr2); |
460 | rinv2 = rinv*rinv; |
461 | if (qqA != 0.0 || qqB != 0.0) |
462 | { |
463 | real dr; |
464 | |
465 | dr = 1.0/rinv; |
466 | v = gmx_erf(ewc_q*dr)gmx_erff(ewc_q*dr)*rinv; |
467 | vc = qqL*v; |
468 | Vexcl_q += vc; |
469 | /* fscal is the scalar force pre-multiplied by rinv, |
470 | * to normalise the relative position vector dx */ |
471 | fscal = rinv2*(vc-qqL*ewc_q*M_2_SQRTPI1.12837916709551257390*exp(-ewc_q*ewc_q*dr2)); |
472 | dvdl_excl_q += (qqB - qqA)*v; |
473 | |
474 | /* The force vector is obtained by multiplication with |
475 | * the relative position vector |
476 | */ |
477 | svmul(fscal, dx, df); |
478 | rvec_inc(f[k], df); |
479 | rvec_dec(f[i], df); |
480 | for (iv = 0; (iv < DIM3); iv++) |
481 | { |
482 | for (jv = 0; (jv < DIM3); jv++) |
483 | { |
484 | dxdf_q[iv][jv] += dx[iv]*df[jv]; |
485 | } |
486 | } |
487 | } |
488 | |
489 | if ((c6A != 0.0 || c6B != 0.0) && EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
490 | { |
491 | rinv6 = rinv2*rinv2*rinv2; |
492 | ewcdr2 = ewc_lj2*dr2; |
493 | ewcdr4 = ewcdr2*ewcdr2; |
494 | v = -rinv6*(1.0 - exp(-ewcdr2)*(1 + ewcdr2 + 0.5*ewcdr4)); |
495 | vc = c6L*v; |
496 | Vexcl_lj += vc; |
497 | /* fscal is the scalar force pre-multiplied by rinv, |
498 | * to normalise the relative position vector dx */ |
499 | fscal = 6.0*vc*rinv2 + c6L*rinv6*exp(-ewcdr2)*ewc_lj2*ewcdr4; |
500 | dvdl_excl_lj += (c6B - c6A)*v; |
501 | |
502 | /* The force vector is obtained by multiplication with |
503 | * the relative position vector |
504 | */ |
505 | svmul(fscal, dx, df); |
506 | rvec_inc(f[k], df); |
507 | rvec_dec(f[i], df); |
508 | for (iv = 0; (iv < DIM3); iv++) |
509 | { |
510 | for (jv = 0; (jv < DIM3); jv++) |
511 | { |
512 | dxdf_lj[iv][jv] += dx[iv]*df[jv]; |
513 | } |
514 | } |
515 | } |
516 | } |
517 | else |
518 | { |
519 | Vexcl_q += qqL*vr0_q; |
520 | dvdl_excl_q += (qqB - qqA)*vr0_q; |
521 | Vexcl_lj += c6L*vr0_lj; |
522 | dvdl_excl_lj += (c6B - c6A)*vr0_lj; |
523 | } |
524 | } |
525 | } |
526 | } |
527 | } |
528 | /* Dipole correction on force */ |
529 | if (dipole_coeff != 0) |
530 | { |
531 | for (j = 0; (j < DIM3); j++) |
532 | { |
533 | f[i][j] -= L1_q*dipcorrA[j]*chargeA[i] |
534 | + lambda_q*dipcorrB[j]*chargeB[i]; |
535 | } |
536 | } |
537 | } |
538 | } |
539 | for (iv = 0; (iv < DIM3); iv++) |
540 | { |
541 | for (jv = 0; (jv < DIM3); jv++) |
542 | { |
543 | vir_q[iv][jv] += 0.5*dxdf_q[iv][jv]; |
544 | vir_lj[iv][jv] += 0.5*dxdf_lj[iv][jv]; |
545 | } |
546 | } |
547 | |
548 | Vself_q[0] = 0; |
549 | Vself_q[1] = 0; |
550 | Vself_lj[0] = 0; |
551 | Vself_lj[1] = 0; |
552 | |
553 | /* Global corrections only on master process */ |
554 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && thread == 0) |
555 | { |
556 | for (q = 0; q < (bFreeEnergy ? 2 : 1); q++) |
557 | { |
558 | if (calc_excl_corr) |
559 | { |
560 | /* Self-energy correction */ |
561 | Vself_q[q] = ewc_q*one_4pi_eps*fr->q2sum[q]*M_1_SQRTPI0.564189583547756; |
562 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
563 | { |
564 | Vself_lj[q] = fr->c6sum[q]*0.5*vr0_lj; |
565 | } |
566 | } |
567 | |
568 | /* Apply surface dipole correction: |
569 | * correction = dipole_coeff * (dipole)^2 |
570 | */ |
571 | if (dipole_coeff != 0) |
572 | { |
573 | if (ewald_geometry == eewg3D) |
574 | { |
575 | Vdipole[q] = dipole_coeff*iprod(mutot[q], mutot[q]); |
576 | } |
577 | else if (ewald_geometry == eewg3DC) |
578 | { |
579 | Vdipole[q] = dipole_coeff*mutot[q][ZZ2]*mutot[q][ZZ2]; |
580 | } |
581 | } |
582 | } |
583 | } |
584 | if (!bFreeEnergy) |
585 | { |
586 | *Vcorr_q = Vdipole[0] - Vself_q[0] - Vexcl_q; |
587 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
588 | { |
589 | *Vcorr_lj = -Vself_lj[0] - Vexcl_lj; |
590 | } |
591 | } |
592 | else |
593 | { |
594 | *Vcorr_q = L1_q*(Vdipole[0] - Vself_q[0]) |
595 | + lambda_q*(Vdipole[1] - Vself_q[1]) |
596 | - Vexcl_q; |
597 | *dvdlambda_q += Vdipole[1] - Vself_q[1] |
598 | - (Vdipole[0] - Vself_q[0]) - dvdl_excl_q; |
599 | if (EVDW_PME(fr->vdwtype)((fr->vdwtype) == evdwPME)) |
600 | { |
601 | *Vcorr_lj = -(L1_lj*Vself_lj[0] + lambda_lj*Vself_lj[1]) - Vexcl_lj; |
602 | *dvdlambda_lj += -Vself_lj[1] + Vself_lj[0] - dvdl_excl_lj; |
603 | } |
604 | } |
605 | |
606 | if (debug) |
607 | { |
608 | fprintf(debug, "Long Range corrections for Ewald interactions:\n"); |
609 | fprintf(debug, "start=%d,natoms=%d\n", start, end-start); |
610 | fprintf(debug, "q2sum = %g, Vself_q=%g c6sum = %g, Vself_lj=%g\n", |
611 | L1_q*fr->q2sum[0]+lambda_q*fr->q2sum[1], L1_q*Vself_q[0]+lambda_q*Vself_q[1], L1_lj*fr->c6sum[0]+lambda_lj*fr->c6sum[1], L1_lj*Vself_lj[0]+lambda_lj*Vself_lj[1]); |
612 | fprintf(debug, "Electrostatic Long Range correction: Vexcl=%g\n", Vexcl_q); |
613 | fprintf(debug, "Lennard-Jones Long Range correction: Vexcl=%g\n", Vexcl_lj); |
614 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && thread == 0) |
615 | { |
616 | if (epsilon_surface > 0 || ewald_geometry == eewg3DC) |
617 | { |
618 | fprintf(debug, "Total dipole correction: Vdipole=%g\n", |
619 | L1_q*Vdipole[0]+lambda_q*Vdipole[1]); |
620 | } |
621 | } |
622 | } |
623 | } |
624 | |
625 | real ewald_charge_correction(t_commrec *cr, t_forcerec *fr, real lambda, |
626 | matrix box, |
627 | real *dvdlambda, tensor vir) |
628 | |
629 | { |
630 | real vol, fac, qs2A, qs2B, vc, enercorr; |
631 | int d; |
632 | |
633 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) |
634 | { |
635 | /* Apply charge correction */ |
636 | vol = box[XX0][XX0]*box[YY1][YY1]*box[ZZ2][ZZ2]; |
637 | |
638 | fac = M_PI3.14159265358979323846*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/(fr->epsilon_r*2.0*vol*vol*sqr(fr->ewaldcoeff_q)); |
639 | |
640 | qs2A = fr->qsum[0]*fr->qsum[0]; |
641 | qs2B = fr->qsum[1]*fr->qsum[1]; |
642 | |
643 | vc = (qs2A*(1 - lambda) + qs2B*lambda)*fac; |
644 | |
645 | enercorr = -vol*vc; |
646 | |
647 | *dvdlambda += -vol*(qs2B - qs2A)*fac; |
648 | |
649 | for (d = 0; d < DIM3; d++) |
650 | { |
651 | vir[d][d] += vc; |
652 | } |
653 | |
654 | if (debug) |
655 | { |
656 | fprintf(debug, "Total charge correction: Vcharge=%g\n", enercorr); |
657 | } |
658 | } |
659 | else |
660 | { |
661 | enercorr = 0; |
662 | } |
663 | |
664 | return enercorr; |
665 | } |