Lost atoms during NVT heating of BKS silica (LAMMPS)

LAMMPS LAMMPS Beginners
, ,
1

Hi,

I’m trying to generate amorphous SiO₂ using the BKS potential in LAMMPS, but I consistently encounter a “Lost atoms” error during the heating stage.

Simulation details

  • LAMMPS (recent version)

  • Units: metal

  • Initial structure: β-cristobalite (288 atoms)

  • Boundary: p p p

Force field

pair_style hybrid/overlay buck/coul/long 5.5 lj/cut 1.2
kspace_style ewald 1.0e-4

Charges:

  • Si = +2.4

  • O = –1.2

The Buckingham parameters are the standard BKS values.

To avoid the Buckingham short-range divergence, I added the repulsive Lennard-Jones overlay suggested in the literature:

  • Si–O: ε = 2.0 eV, σ = 1.2 Å

  • O–O: ε = 2.6 eV, σ = 1.6 Å

Simulation protocol

  1. Energy minimization (CG)

  2. Assign velocities at 298 K

  3. fix nve

  4. fix temp/berendsen 298 6000 1.0

  5. Heat from 298 K to 6000 K over 1000 ps (≈5.7 K/ps)

The simulation eventually stops with

ERROR: Lost atoms: original 288 current 286

I observe the same issue even when the target temperature is reduced (e.g., 4000 K).

What I’ve already tried

To rule out the common causes, I have already tested several alternatives:

  • Reduced the timestep from 1.0 fs to 0.5 fs.

  • Performed the heating in the NPT ensemble (instead of NVT), allowing the simulation cell to expand.

  • Successfully heated the system to 6000 K and even 8000 K under NPT without immediate instability.

  • The same problem still appears when switching to the NVT heating stage.

Because of this, I’m beginning to suspect that the issue may not simply be due to the timestep or the inability of the simulation box to expand.

I’m wondering whether the instability could instead be related to:

  • insufficient relaxation of the initial β-cristobalite structure,

  • an issue with my BKS + LJ overlay implementation,

  • the thermostat choice or damping parameter,

  • an incorrect equilibrium density before switching to NVT,

  • or some other aspect of the simulation setup that I’m overlooking.

Has anyone encountered similar behavior when generating amorphous silica using the BKS potential? Any suggestions on obtaining a stable melt–quench simulation would be greatly appreciated.

Thank you!

2

These are all items that you can test for yourself. I doubt that anybody will do it for you.

There should be several discussions from the past about generating amorphous SiO2 archived. I suggest you dig in and systematically evaluate what people have tried and what was suggested as remedy for similar problems. This can be tedious, but it is not likely that somebody will do this for you. This kind of legwork is just the normal crunch of doing research.

3

Thank you for your suggestion. I have already gone through several previous discussions and tried many of the commonly recommended approaches. For example, I reduced the timestep from 1.0 fs to 0.5 fs, performed NPT heating up to both 6000 K and 8000 K, and experimented with different thermostat/barostat settings.

However, I’m still encountering instability. During NPT heating, I observe large pressure oscillations and significant density changes, and in some cases the simulation box contracts before eventually leading to a “Lost atoms” error.

I’m therefore trying to determine whether this behavior indicates a more fundamental issue with my simulation setup (e.g., the force-field implementation or initial structure) rather than simply an integration or thermostat problem. If anyone has encountered similar behavior with BKS silica, I’d appreciate any further insights.

4

The magnitude of pressure fluctuations are usually correlated with the system size.
For larger systems, the fluctuations are smaller.

Please also note that a kspace convergence of 1.0e-4 is sufficient for converged forces (thanks to error cancellation) but not for a converged pressure. You may need 1.0e-6.

Please also note that using kspace style ewald is suggesting your system is very small since its performance for larger systems is very bad.

5

Thank you for your comments.

One thing I’ve observed is that, during NPT heating, the pressure exhibits large oscillations and the density initially increases. Eventually, the simulation box contracts significantly, after which atoms are lost. This behavior persists even after reducing the timestep to 0.5 fs.

Since there are published BKS melt–quench studies with systems of comparable size (around 300 atoms), I’m wondering whether the issue is more likely related to my force-field implementation or initial structure rather than the system size itself.

Does the sequence of pressure oscillations → density increase/box contraction → lost atoms suggest a particular type of instability

bks.txt (175.8 KB)

6

That is for you to figure out. You have all the published information, you just need to figure you
what you are doing different.

That is a question about your research and not about LAMMPS. I an neither your adviser nor your collaborator and thus this is not for me to advise you on.

Unless you can produce some proof and a usable simple reproduce input deck (see for example the discussion on fix gcmc misbehaving that we just resolved for a case that was on-topic) showing that LAMMPS does not behave like its documentation says (i.e. is producing provably incorrect results and not just results you don’t want) this is off-topic and something you need ask your collaborators or adviser or tutor for assistance with.

7

Hello, I have also been working on modeling in this area recently. Regarding the Lennard-Jones overlay suggested in the literature that you mentioned, could you please tell me which specific paper it is? I would appreciate it if you could provide the details. In addition, I am using the NPT ensemble, and I think there are some aspects where we could have more in‑depth discussions.

Sincerely,
Young

8

This reply is also relevant. Just cite the matsci forum for the protocol :wink: