What you report is very common when you have densely charged ions (as in your case with lithium) or strong H-bonds involving induced dipoles (such as intramolecular H-bonds in ethylene glycol).
The physical reason is that at very short range the point charge representation is not robust and we need to represent smeared electron densities. This is done using a “damping” or smearing function for the electrostatic interactions at short range.
I’d recommend using a Tang-Toennis damping function between your lithium cations and the polarizable atoms of your system (at least negative ones), as in Goloviznina et al, J Chem Theory Comput 17, 1606 (2021). We took this from simulations of high-temperature molten salts (papers by Paul Madden).
We parameterised this damping function and implemented it in pair_style coul/tt. We also wrote a script (coul_tt, which you will find in gothub/paduagroup/clandpol) which adds these terms to the input files.
We have very stable trajectories and the Drude degrees of freedom are kept around 1 K with no spikes.
If you change the Drude force constant without changing the charges you are changing the polarizability of the atoms.
Thank you very much for your helpful answer.
Using the TT damping function between Li and polarizable atoms, my simulation is now going well.
Should the damping function be used for larger cations such as Na and K as well?
What about H atoms of water molecules and the additional site M in the case of four-site water models? Can they also cause a polarization catastrophe and is it required to apply the same damping to their interactions with polarizable atoms?
If you see spikes in the temperature of the Drude degrees of freedom then I’d say use the TT damping function.
Maybe for Na it is necessary. I don’t know for K.
We don’t use it on 4-site water models, where the M site is well inside the O LJ site. Only for small ions and some strong intramolecular H bonds.
So, is it wrong to use the damping function between H atoms of water molecules and Drude cores/Drude particles?
In the case of the systems containing SPC/E water molecules and polarizable ions, the simulations go well when the TT damping function is used only between cations and polarizable atoms, as you recommended before (however, adding the same damping function between H of water and polarizable atoms decreases the Drude temperature fluctuations). In the case of polarizable water models based on Drude oscillators, however, the simulations go well only if the TT damping function is used for the reactions between polarizable atoms and both cations and H atoms of water molecules (if I don’t use the damping function for H, I receive “Non-numeric box dimensions - simulation unstable” error).