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?
All usable features in LAMMPS are documented. This is a requirement for contributing to the distribution. So if it is not mentioned in the manual it is not available. You would need to modify the source to implement it.
From my memory, the Drude pairs in LAMMPS can be connected with any (suitable) bond_style available in LAMMPS. So you can look through that list and see if any suits (not bond_style quartic, funny enough) and use those to link the Drude pairs. Maybe bond_style table will suit your use.
Bond style table works as a charm! Didn’t know about that one. (It is not so easy to know all the features of LAMMPS, thankfully there is always help to get in this forum :)) Thanks a lot!
Hello, I too am encountering the issue of polarization catastrophe. I’m curious, were you able to fully resolve this bond-breaking problem caused by the polarization catastrophe using the Bond Style Table?
Please note that the paper he’s referencing is one particular example of polarisation catastrophe in Drude oscillators. The more usual approach is to use Thole damping (Google it), which is implemented in LAMMPS with the thole styles: pair_style thole command — LAMMPS documentation
I was not able to fully resolve it, no. But I did manage to stabilize the system sufficiently to enable running long enough simulations to obtain the properties I sought.
I am using the Thole damping function and a small timestep around 0.25 fs, but still encountering the polarization catastrophe. A more recent treatment is to introduce a “hard-wall” https://pubs.acs.org/doi/10.1021/jp402860e . Unfortunately, this feature is only implemented in Gromacs and OpenMM. I am still trying to figure out whether there is an alternative approach.
In my (extremely personal) opinion this “hard wall” approach (of “reflecting” the Drude charge off a containing sphere if it is detected to leave it mid-timestep) is not very sound. It makes simulations non-reversible, and it is not documented in the official GROMACS documentation (it seems to be maintained on a separate lab’s fork).
The class2 and nonlinear bond types are worth experimenting with – they are approximately quadratic near the equilibrium distance (of 0 for this use case) but give much stronger forces at the periphery. You can also use the lepton bond style to implement your own expression for testing purposes.