Negative dielectric, it's a feature, not a bug!

I don’t think that your change is doing what you expect it to do. The only thing that the variable q2 is used for is for error estimation and to get a guess of g_ewald. It doesn’t directly affect the dynamics in any way.

force->qqrd2e is already being divided by the dielectric constant in force.cpp:114 in Force::init(). Previously, force->qqrd2e was being divided by the dielectric constant twice, which is not correct (hence previously your negative sign on the dielectric was canceled out). This directly affects only the accuracy metric, in that for a dielectric constant less than one, one would get better than requested accuracy. If the dielectric constant was greater than one though, then a user would get less accuracy than requested.

My guess is that LAMMPS doesn’t know to estimate the accuracy or g_ewald when the dielectric constant is negative, and that is why it is stalling. Just try setting the dielectric constant equal to 1.0 (not -1.0) in the current version, and I think you will get similar behavior as before (or maybe even better accuracy because you could get a better estimation g_ewald).


Yes, that makes sense. Perhaps the variable q2 should use the absolute value of the dielectric constant, may have to check that.

I should warn you that the Ewald sum is on shaky theoretical grounds when the system is not charge neutral (basically you are assuming a neutralizing background plasma at the infinite boundary). I hope you have thought hard about how this change may affect Ewald sum convergence, as presumably all your masses (and hence charges) add up to a large positive value.

Best regards,