Skip to content

Commit fc8f8fe

Browse files
authored
Make ELC tuning follow literature error-formula (#5280)
The ELC error formula now matches with equation 18 in Arnold 2002 (https://doi.org/10.1063/1.1491955). ELC and ELCIC can now tune to slightly larger summation limits.
1 parent b4a5b03 commit fc8f8fe

1 file changed

Lines changed: 8 additions & 7 deletions

File tree

src/core/electrostatics/elc.cpp

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -961,22 +961,23 @@ double ElectrostaticLayerCorrection::tune_far_cut() const {
961961
auto const box_l_y_inv = box_geo.length_inv()[1];
962962
auto const min_inv_boxl = std::min(box_l_x_inv, box_l_y_inv);
963963
auto const box_l_z = box_geo.length()[2];
964+
auto const h = elc.box_h;
964965
// adjust lz according to dielectric layer method
965-
auto const lz =
966-
(elc.dielectric_contrast_on) ? elc.box_h + elc.space_layer : box_l_z;
966+
auto const lz = (elc.dielectric_contrast_on) ? h + elc.space_layer : box_l_z;
967967

968968
auto tuned_far_cut = min_inv_boxl;
969969
double err;
970970
do {
971+
// following equation 18 in arnold02d
971972
auto const pref = 2. * std::numbers::pi * tuned_far_cut;
972973
auto const sum = pref + 2. * (box_l_x_inv + box_l_y_inv);
973-
auto const den = -expm1(-pref * lz);
974-
auto const num1 = exp(pref * (elc.box_h - lz));
975-
auto const num2 = exp(-pref * (elc.box_h + lz));
974+
auto const den = expm1(pref * lz);
975+
auto const num1 = exp(pref * h);
976+
auto const num2 = 1. / num1; // exp(-pref * h);
976977

977978
err = 0.5 / den *
978-
(num1 * (sum + 1. / (lz - elc.box_h)) / (lz - elc.box_h) +
979-
num2 * (sum + 1. / (lz + elc.box_h)) / (lz + elc.box_h));
979+
(num1 / (lz - h) * (sum + 1. / (lz - h)) +
980+
num2 / (lz + h) * (sum + 1. / (lz + h)));
980981

981982
tuned_far_cut += min_inv_boxl;
982983
} while (err > elc.maxPWerror and tuned_far_cut < maximal_far_cut);

0 commit comments

Comments
 (0)