[lammps-users] High pressure and atoms lost in CAS glass simulation

Dear LAMMPS community,

I am trying to equilibrate a calcium aluminosilicate glass system using the force field developed by Bouhadja et al. (J. Chem. Phys. 138, 224510, 2013), and also used in J. Chem. Phys. 141, 024507 (2014) (see table III in that paper for the force field parameters). I have copied to the bottom of the email the input file that I’m using.

Briefly, I generate the initial system at random (at a slightly lower density than reported in the papers: 3000 atoms in a cubic box 40 A side vs. 2995 atoms in a box of 35 A side) and use nve/limit to remove close contacts, which is something that has worked well for me in the past. However, when I run the simulation I get extremely high pressures and when I switch to nvt the system blows up (ERROR: Lost atoms…). After the 60,000 steps of nve/limit at high temperature any initial bad geometry should have been dealt with.

I suspect that the value of some parameter somewhere is wrong, but I have double checked the values of the force field parameters and everything else looks fine to me. I would really appreciate any suggestions from a fresh pair of eyes.

Thanks in advance!

#-------------------------- Initial Setup -------------------------------------#

units real

atom_style charge

dimension 3

boundary p p p

timestep 1.0

Hi Luis,

Hi Gideon,

Thanks for the reply! The nvt is actually at 5000 K (that was a typo). I used a very low Langevin damping constant to try to keep the temperature stable during the nve/limit stage. The 60 ps (or even less) are usually enough to remove close contacts, before starting a true equilibration stage. I wanted to avoid it, but I may have to use a smaller time step at least during the equilibration.

Thanks again!




have you applied the same procedure to the same family of potentials before?
they are notorious for not being sufficiently repulsive and thus when you initialize a system the way you do, you will have close contacts and then the coulomb interaction will become dominant and those close contacts can never separate.

what could work in that case would be start initially with some other potential and since you are far, far away from equilibrium you don’t really need good parameters, just an educated guess that will put you in the ballpark should do, e.g. like this:

pair_style lj/cut/coul/cut 8.0 8.0
pair_coeff * * 0.2000 3.5

with that you can then run some minimization and a first thermalization and then switch to your actually desired potential style and parameters and minimize and thermalize again.
the steep repulsion will avoid any “coulomb catastrophe” as you are currently seeing.


Thanks Axel,

That actually worked fine. I have no prior experience working with these potentials, I should have suspected that electrostatics were lurking in there somewhere.

Thanks again!