2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
43 #include "gromacs/math/utilities.h"
45 #include "types/commrec.h"
48 #include "gromacs/utility/smalloc.h"
51 #include "gromacs/fileio/futil.h"
55 real
calc_ewaldcoeff_q(real rc
, real dtol
)
57 real x
= 5, low
, high
;
66 while (gmx_erfc(x
*rc
) > dtol
);
68 n
= i
+60; /* search tolerance is 2^-60 */
71 for (i
= 0; i
< n
; i
++)
74 if (gmx_erfc(x
*rc
) > dtol
)
86 static real
ewald_function_lj(real x
, real rc
)
88 real xrc
, xrc2
, xrc4
, factor
;
93 factor
= exp(-xrc2
)*(1 + xrc2
+ xrc4
/2.0);
95 factor
= expf(-xrc2
)*(1 + xrc2
+ xrc4
/2.0);
101 real
calc_ewaldcoeff_lj(real rc
, real dtol
)
103 real x
= 5, low
, high
;
111 while (ewald_function_lj(x
, rc
) > dtol
);
113 n
= i
+ 60; /* search tolerance is 2^-60 */
116 for (i
= 0; i
< n
; ++i
)
118 x
= (low
+ high
) / 2.0;
119 if (ewald_function_lj(x
, rc
) > dtol
)
131 /* There's nothing special to do here if just masses are perturbed,
132 * but if either charge or type is perturbed then the implementation
133 * requires that B states are defined for both charge and type, and
134 * does not optimize for the cases where only one changes.
136 * The parameter vectors for B states are left undefined in atoms2md()
137 * when either FEP is inactive, or when there are no mass/charge/type
138 * perturbations. The parameter vectors for LJ-PME are likewise
139 * undefined when LJ-PME is not active. This works because
140 * bHaveChargeOrTypePerturbed handles the control flow. */
141 void ewald_LRcorrection(int start
, int end
,
142 t_commrec
*cr
, int thread
, t_forcerec
*fr
,
143 real
*chargeA
, real
*chargeB
,
144 real
*C6A
, real
*C6B
,
145 real
*sigmaA
, real
*sigmaB
,
146 real
*sigma3A
, real
*sigma3B
,
147 gmx_bool bHaveChargeOrTypePerturbed
,
148 gmx_bool calc_excl_corr
,
149 t_blocka
*excl
, rvec x
[],
150 matrix box
, rvec mu_tot
[],
151 int ewald_geometry
, real epsilon_surface
,
152 rvec
*f
, tensor vir_q
, tensor vir_lj
,
153 real
*Vcorr_q
, real
*Vcorr_lj
,
154 real lambda_q
, real lambda_lj
,
155 real
*dvdlambda_q
, real
*dvdlambda_lj
)
157 int i
, i1
, i2
, j
, k
, m
, iv
, jv
, q
;
159 double Vexcl_q
, dvdl_excl_q
, dvdl_excl_lj
; /* Necessary for precision */
162 real v
, vc
, qiA
, qiB
, dr2
, rinv
, enercorr
;
163 real Vself_q
[2], Vself_lj
[2], Vdipole
[2], rinv2
, ewc_q
= fr
->ewaldcoeff_q
, ewcdr
;
164 real ewc_lj
= fr
->ewaldcoeff_lj
, ewc_lj2
= ewc_lj
* ewc_lj
;
165 real c6Ai
= 0, c6Bi
= 0, c6A
= 0, c6B
= 0, ewcdr2
, ewcdr4
, c6L
= 0, rinv6
;
166 rvec df
, dx
, mutot
[2], dipcorrA
, dipcorrB
;
167 tensor dxdf_q
, dxdf_lj
;
168 real vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
169 real L1_q
, L1_lj
, dipole_coeff
, qqA
, qqB
, qqL
, vr0_q
, vr0_lj
= 0;
170 gmx_bool bMolPBC
= fr
->bMolPBC
;
171 gmx_bool bDoingLBRule
= (fr
->ljpme_combination_rule
== eljpmeLB
);
173 /* This routine can be made faster by using tables instead of analytical interactions
174 * However, that requires a thorough verification that they are correct in all cases.
177 one_4pi_eps
= ONE_4PI_EPS0
/fr
->epsilon_r
;
178 vr0_q
= ewc_q
*M_2_SQRTPI
;
179 if (EVDW_PME(fr
->vdwtype
))
181 vr0_lj
= -pow(ewc_lj
, 6)/6.0;
192 L1_lj
= 1.0-lambda_lj
;
193 /* Note that we have to transform back to gromacs units, since
194 * mu_tot contains the dipole in debye units (for output).
196 for (i
= 0; (i
< DIM
); i
++)
198 mutot
[0][i
] = mu_tot
[0][i
]*DEBYE2ENM
;
199 mutot
[1][i
] = mu_tot
[1][i
]*DEBYE2ENM
;
204 switch (ewald_geometry
)
207 if (epsilon_surface
!= 0)
210 2*M_PI
*ONE_4PI_EPS0
/((2*epsilon_surface
+ fr
->epsilon_r
)*vol
);
211 for (i
= 0; (i
< DIM
); i
++)
213 dipcorrA
[i
] = 2*dipole_coeff
*mutot
[0][i
];
214 dipcorrB
[i
] = 2*dipole_coeff
*mutot
[1][i
];
219 dipole_coeff
= 2*M_PI
*one_4pi_eps
/vol
;
220 dipcorrA
[ZZ
] = 2*dipole_coeff
*mutot
[0][ZZ
];
221 dipcorrB
[ZZ
] = 2*dipole_coeff
*mutot
[1][ZZ
];
224 gmx_incons("Unsupported Ewald geometry");
229 fprintf(debug
, "dipcorr = %8.3f %8.3f %8.3f\n",
230 dipcorrA
[XX
], dipcorrA
[YY
], dipcorrA
[ZZ
]);
231 fprintf(debug
, "mutot = %8.3f %8.3f %8.3f\n",
232 mutot
[0][XX
], mutot
[0][YY
], mutot
[0][ZZ
]);
235 if (EVDW_PME(fr
->vdwtype
))
239 if ((calc_excl_corr
|| dipole_coeff
!= 0) && !bHaveChargeOrTypePerturbed
)
241 for (i
= start
; (i
< end
); i
++)
243 /* Initiate local variables (for this i-particle) to 0 */
244 qiA
= chargeA
[i
]*one_4pi_eps
;
245 if (EVDW_PME(fr
->vdwtype
))
256 i2
= excl
->index
[i
+1];
258 /* Loop over excluded neighbours */
259 for (j
= i1
; (j
< i2
); j
++)
263 * First we must test whether k <> i, and then,
264 * because the exclusions are all listed twice i->k
265 * and k->i we must select just one of the two. As
266 * a minor optimization we only compute forces when
267 * the charges are non-zero.
271 qqA
= qiA
*chargeA
[k
];
272 if (EVDW_PME(fr
->vdwtype
))
277 c6A
*= pow(0.5*(sigmaA
[i
]+sigmaA
[k
]), 6)*sigma3A
[k
];
280 if (qqA
!= 0.0 || c6A
!= 0.0)
285 rvec_sub(x
[i
], x
[k
], dx
);
288 /* Cheap pbc_dx, assume excluded pairs are at short distance. */
289 for (m
= DIM
-1; (m
>= 0); m
--)
291 if (dx
[m
] > 0.5*box
[m
][m
])
293 rvec_dec(dx
, box
[m
]);
295 else if (dx
[m
] < -0.5*box
[m
][m
])
297 rvec_inc(dx
, box
[m
]);
302 /* Distance between two excluded particles
303 * may be zero in the case of shells
307 rinv
= gmx_invsqrt(dr2
);
315 vc
= qqA
*gmx_erf(ewcdr
)*rinv
;
318 /* Relative accuracy at R_ERF_R_INACC of 3e-10 */
319 #define R_ERF_R_INACC 0.006
321 /* Relative accuracy at R_ERF_R_INACC of 2e-5 */
322 #define R_ERF_R_INACC 0.1
324 /* fscal is the scalar force pre-multiplied by rinv,
325 * to normalise the relative position vector dx */
326 if (ewcdr
> R_ERF_R_INACC
)
328 fscal
= rinv2
*(vc
- qqA
*ewc_q
*M_2_SQRTPI
*exp(-ewcdr
*ewcdr
));
332 /* Use a fourth order series expansion for small ewcdr */
333 fscal
= ewc_q
*ewc_q
*qqA
*vr0_q
*(2.0/3.0 - 0.4*ewcdr
*ewcdr
);
336 /* The force vector is obtained by multiplication with
337 * the relative position vector
339 svmul(fscal
, dx
, df
);
342 for (iv
= 0; (iv
< DIM
); iv
++)
344 for (jv
= 0; (jv
< DIM
); jv
++)
346 dxdf_q
[iv
][jv
] += dx
[iv
]*df
[jv
];
353 rinv6
= rinv2
*rinv2
*rinv2
;
354 ewcdr2
= ewc_lj2
*dr2
;
355 ewcdr4
= ewcdr2
*ewcdr2
;
356 /* We get the excluded long-range contribution from -C6*(1-g(r))
357 * g(r) is also defined in the manual under LJ-PME
359 vc
= -c6A
*rinv6
*(1.0 - exp(-ewcdr2
)*(1 + ewcdr2
+ 0.5*ewcdr4
));
361 /* The force is the derivative of the potential vc.
362 * fscal is the scalar force pre-multiplied by rinv,
363 * to normalise the relative position vector dx */
364 fscal
= 6.0*vc
*rinv2
+ c6A
*rinv6
*exp(-ewcdr2
)*ewc_lj2
*ewcdr4
;
366 /* The force vector is obtained by multiplication with
367 * the relative position vector
369 svmul(fscal
, dx
, df
);
372 for (iv
= 0; (iv
< DIM
); iv
++)
374 for (jv
= 0; (jv
< DIM
); jv
++)
376 dxdf_lj
[iv
][jv
] += dx
[iv
]*df
[jv
];
383 Vexcl_q
+= qqA
*vr0_q
;
384 Vexcl_lj
+= c6A
*vr0_lj
;
390 /* Dipole correction on force */
391 if (dipole_coeff
!= 0)
393 for (j
= 0; (j
< DIM
); j
++)
395 f
[i
][j
] -= dipcorrA
[j
]*chargeA
[i
];
400 else if (calc_excl_corr
|| dipole_coeff
!= 0)
402 for (i
= start
; (i
< end
); i
++)
404 /* Initiate local variables (for this i-particle) to 0 */
405 qiA
= chargeA
[i
]*one_4pi_eps
;
406 qiB
= chargeB
[i
]*one_4pi_eps
;
407 if (EVDW_PME(fr
->vdwtype
))
420 i2
= excl
->index
[i
+1];
422 /* Loop over excluded neighbours */
423 for (j
= i1
; (j
< i2
); j
++)
428 qqA
= qiA
*chargeA
[k
];
429 qqB
= qiB
*chargeB
[k
];
430 if (EVDW_PME(fr
->vdwtype
))
436 c6A
*= pow(0.5*(sigmaA
[i
]+sigmaA
[k
]), 6)*sigma3A
[k
];
437 c6B
*= pow(0.5*(sigmaB
[i
]+sigmaB
[k
]), 6)*sigma3B
[k
];
440 if (qqA
!= 0.0 || qqB
!= 0.0 || c6A
!= 0.0 || c6B
!= 0.0)
445 qqL
= L1_q
*qqA
+ lambda_q
*qqB
;
446 if (EVDW_PME(fr
->vdwtype
))
448 c6L
= L1_lj
*c6A
+ lambda_lj
*c6B
;
450 rvec_sub(x
[i
], x
[k
], dx
);
453 /* Cheap pbc_dx, assume excluded pairs are at short distance. */
454 for (m
= DIM
-1; (m
>= 0); m
--)
456 if (dx
[m
] > 0.5*box
[m
][m
])
458 rvec_dec(dx
, box
[m
]);
460 else if (dx
[m
] < -0.5*box
[m
][m
])
462 rvec_inc(dx
, box
[m
]);
469 rinv
= gmx_invsqrt(dr2
);
471 if (qqA
!= 0.0 || qqB
!= 0.0)
476 v
= gmx_erf(ewc_q
*dr
)*rinv
;
479 /* fscal is the scalar force pre-multiplied by rinv,
480 * to normalise the relative position vector dx */
481 fscal
= rinv2
*(vc
-qqL
*ewc_q
*M_2_SQRTPI
*exp(-ewc_q
*ewc_q
*dr2
));
482 dvdl_excl_q
+= (qqB
- qqA
)*v
;
484 /* The force vector is obtained by multiplication with
485 * the relative position vector
487 svmul(fscal
, dx
, df
);
490 for (iv
= 0; (iv
< DIM
); iv
++)
492 for (jv
= 0; (jv
< DIM
); jv
++)
494 dxdf_q
[iv
][jv
] += dx
[iv
]*df
[jv
];
499 if ((c6A
!= 0.0 || c6B
!= 0.0) && EVDW_PME(fr
->vdwtype
))
501 rinv6
= rinv2
*rinv2
*rinv2
;
502 ewcdr2
= ewc_lj2
*dr2
;
503 ewcdr4
= ewcdr2
*ewcdr2
;
504 v
= -rinv6
*(1.0 - exp(-ewcdr2
)*(1 + ewcdr2
+ 0.5*ewcdr4
));
507 /* fscal is the scalar force pre-multiplied by rinv,
508 * to normalise the relative position vector dx */
509 fscal
= 6.0*vc
*rinv2
+ c6L
*rinv6
*exp(-ewcdr2
)*ewc_lj2
*ewcdr4
;
510 dvdl_excl_lj
+= (c6B
- c6A
)*v
;
512 /* The force vector is obtained by multiplication with
513 * the relative position vector
515 svmul(fscal
, dx
, df
);
518 for (iv
= 0; (iv
< DIM
); iv
++)
520 for (jv
= 0; (jv
< DIM
); jv
++)
522 dxdf_lj
[iv
][jv
] += dx
[iv
]*df
[jv
];
529 Vexcl_q
+= qqL
*vr0_q
;
530 dvdl_excl_q
+= (qqB
- qqA
)*vr0_q
;
531 Vexcl_lj
+= c6L
*vr0_lj
;
532 dvdl_excl_lj
+= (c6B
- c6A
)*vr0_lj
;
538 /* Dipole correction on force */
539 if (dipole_coeff
!= 0)
541 for (j
= 0; (j
< DIM
); j
++)
543 f
[i
][j
] -= L1_q
*dipcorrA
[j
]*chargeA
[i
]
544 + lambda_q
*dipcorrB
[j
]*chargeB
[i
];
549 for (iv
= 0; (iv
< DIM
); iv
++)
551 for (jv
= 0; (jv
< DIM
); jv
++)
553 vir_q
[iv
][jv
] += 0.5*dxdf_q
[iv
][jv
];
554 vir_lj
[iv
][jv
] += 0.5*dxdf_lj
[iv
][jv
];
563 /* Global corrections only on master process */
564 if (MASTER(cr
) && thread
== 0)
566 for (q
= 0; q
< (bHaveChargeOrTypePerturbed
? 2 : 1); q
++)
570 /* Self-energy correction */
571 Vself_q
[q
] = ewc_q
*one_4pi_eps
*fr
->q2sum
[q
]*M_1_SQRTPI
;
572 if (EVDW_PME(fr
->vdwtype
))
574 Vself_lj
[q
] = fr
->c6sum
[q
]*0.5*vr0_lj
;
578 /* Apply surface dipole correction:
579 * correction = dipole_coeff * (dipole)^2
581 if (dipole_coeff
!= 0)
583 if (ewald_geometry
== eewg3D
)
585 Vdipole
[q
] = dipole_coeff
*iprod(mutot
[q
], mutot
[q
]);
587 else if (ewald_geometry
== eewg3DC
)
589 Vdipole
[q
] = dipole_coeff
*mutot
[q
][ZZ
]*mutot
[q
][ZZ
];
594 if (!bHaveChargeOrTypePerturbed
)
596 *Vcorr_q
= Vdipole
[0] - Vself_q
[0] - Vexcl_q
;
597 if (EVDW_PME(fr
->vdwtype
))
599 *Vcorr_lj
= -Vself_lj
[0] - Vexcl_lj
;
604 *Vcorr_q
= L1_q
*(Vdipole
[0] - Vself_q
[0])
605 + lambda_q
*(Vdipole
[1] - Vself_q
[1])
607 *dvdlambda_q
+= Vdipole
[1] - Vself_q
[1]
608 - (Vdipole
[0] - Vself_q
[0]) - dvdl_excl_q
;
609 if (EVDW_PME(fr
->vdwtype
))
611 *Vcorr_lj
= -(L1_lj
*Vself_lj
[0] + lambda_lj
*Vself_lj
[1]) - Vexcl_lj
;
612 *dvdlambda_lj
+= -Vself_lj
[1] + Vself_lj
[0] - dvdl_excl_lj
;
618 fprintf(debug
, "Long Range corrections for Ewald interactions:\n");
619 fprintf(debug
, "start=%d,natoms=%d\n", start
, end
-start
);
620 fprintf(debug
, "q2sum = %g, Vself_q=%g c6sum = %g, Vself_lj=%g\n",
621 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]);
622 fprintf(debug
, "Electrostatic Long Range correction: Vexcl=%g\n", Vexcl_q
);
623 fprintf(debug
, "Lennard-Jones Long Range correction: Vexcl=%g\n", Vexcl_lj
);
624 if (MASTER(cr
) && thread
== 0)
626 if (epsilon_surface
> 0 || ewald_geometry
== eewg3DC
)
628 fprintf(debug
, "Total dipole correction: Vdipole=%g\n",
629 L1_q
*Vdipole
[0]+lambda_q
*Vdipole
[1]);
635 real
ewald_charge_correction(t_commrec
*cr
, t_forcerec
*fr
, real lambda
,
637 real
*dvdlambda
, tensor vir
)
640 real vol
, fac
, qs2A
, qs2B
, vc
, enercorr
;
645 /* Apply charge correction */
646 vol
= box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
648 fac
= M_PI
*ONE_4PI_EPS0
/(fr
->epsilon_r
*2.0*vol
*vol
*sqr(fr
->ewaldcoeff_q
));
650 qs2A
= fr
->qsum
[0]*fr
->qsum
[0];
651 qs2B
= fr
->qsum
[1]*fr
->qsum
[1];
653 vc
= (qs2A
*(1 - lambda
) + qs2B
*lambda
)*fac
;
657 *dvdlambda
+= -vol
*(qs2B
- qs2A
)*fac
;
659 for (d
= 0; d
< DIM
; d
++)
666 fprintf(debug
, "Total charge correction: Vcharge=%g\n", enercorr
);