Hello

So I want to simulate a system composed of water and Na+ and Cl- ions using drude oscillators. I follow the method described here: 8.5.7. Tutorial for Thermalized Drude oscillators in LAMMPS — LAMMPS documentation, using the SWM4-NDP water model. This paper describes how to model polarizable ions: https://doi.org/10.1021/ct900576a. The paper describes the “polarization catastrophe” which I think I have encountered regarding the Cl- ions. I guess the problem is that when the drude and core particle have high negative and positive atomic charges (which is the case for Cl-, Br-, I-), then the drude particle can get far away from the core particle and the simulation becomes unstable. As explained in the paper: “The simplest treatment is to introduce an additional anharmonic restoring force to prevent excessively large excursions of the Drude particle away from the atom”. For highly polarizable anions such as Br- and I- (and Cl- I think), they implement an anharmonic restoring force active only beyond a certain stretching distance:

U_hyp = K_hyp · (ΔR - ΔR_cut)^4 , where ΔR_cut is the cutoff distance and ΔR is the distance.

K_hyp is chosen to be a very high value (40 000 kcal/mol/Å^4) so that this potential effectively reduces the distance between the drude and core particle, also reducing the induced dipole and stabilizes the system. This has been implemented in CHARMM apparently.

So the question is: is this implemented in LAMMPS as well, (I haven’t found anything about it though), or is there an alternative way of doing this?