infinite relative permittivity is incompletely implemented
0 ir->epsilon_r[f] is supposed to signal infinite relative permittivity, i.e. fully conducting, so no electrostatic potential.As partly reported on gmx-developers today by Mark Tianwu Zang:
- ewald_LRcorrection, ewald_charge_correction, do_ewald and calc_verlet_buffer_size use ir->epsilon_r without checking for zero
- solve_pme_yzx does likewise (indirectly)
- generalized Born code does likewise (genborn.c, and all the way into the group kernels)
So various code will divide by zero.
I think the PME, GB and ic data structures should store a multiplicative constant (as fr and ic try to do), that their init routines should all apply the same logic when using inputrec->epsilon_r[f], and that only inputrec should have members epsilon_r[f].
Because this goes everywhere, it's probably easiest to do in master branch with some tests that (0 epsilon_r[f]) => (no crash and no Coulomb) for all coulombtype and GB. I have some as-yet unpublished toys that will support doing that. Then we cherry-pick the result back.
#1 Updated by Mark Abraham almost 6 years ago
Of course, it doesn't make much sense to trigger these code paths with 0 == epsilon_r[f] - all you do is waste time. tpbconv -zeroq supports doing a rerun that happens to have no/few electrostatic interactions, because the searching notices atoms have no charge and the kernels become an empty loop.
We can probably extend coulombtype to include None in 5.0 - a null electrostatics model might be useful for testing code, or with rerun.
The question becomes "why do we support infinite relative permittivity?"
#3 Updated by Berk Hess almost 5 years ago
Infinite epsilon_rf is fully and properly supported.
Infinite epsilon_r simply completely turns off all electrostatic interactions. I don't see any need to support this with PME, GB, etc. grompp should either give a fatal_error and suggest to switch to plain cut-off with epsilon_r=0 for optimal performance.