[lammps-users] Problem with Drude particles

Dear Lammps users,

I want to simulate an aqueous electrolyte using CL&Pol polarizable force field (https://doi.org/10.1002/wcms.1572), but I receive errors such as “out of range atoms”, “bond/angle/shake atoms missing”, and “Non-numeric pressure” after a few hundred to a few thousand time steps, which I think indicate bad dynamics. It might be interesting that at the time step when I get the error, the Drude temperature goes to zero.
I have attached my input files here. To create the initial configuration, Drude particles are added to the atoms of anions and cations using the repository “https://github.com/agiliopadua/clandpol” while water is modeled using SPC/E model. It seems that there is no problem with the initial configuration: I am using the box dimensions obtained from my non-polarizable NPT simulation and there are no overlapping atoms in my initial geometry (I made sure about this by using the command “delete_atoms overlap”). I tried to solve the problem by decreasing the time step, but I get errors even when the time step is set to very small values. The only reason I can think of for the errors I receive is that due to the high polarizability of the atoms forming anions, Drude particles move too far from their cores and make the simulation unstable! Indeed, the simulation goes well when the polarizability of these atoms is decreased compared to what was reported in the literature, but using the values reported in the literature I get the mentioned errors. Just to check whether the problem is somehow caused by the excessive motion of Drude particles, I simply increased the force constants of the Drude-core harmonic bonds (I mean K_drude=k_D/2 that is specified as the first coefficient for DP-DC bonds in data file) while keeping the Drude charges as calculated from “alpha=q_D^2/k_D” considering k_D=1000 kcal/molA2 and alphas from the literature. I found that when K_drude is increased to around 1000kcal/molA2 (twice as much as it should be), the errors disappear and the simulation goes well! To me, it means that the motion of the Drude particles should be somehow controlled (but how?).

Does anyone have an idea to fix the problem?


in.lmp (1.6 KB)

pair-sc.lmp (4.13 KB)

pair-drude.lmp (5.98 KB)

data-p.lmp (431 KB)