How to rectify the error "Lost atoms: original 28437 current 3659" in LAMMPS?

Question

I have modeled a complex structure of a protein and a Carbon nano-cone (used as an indenter to the protein) using PACKMOL. The LAMMPS data file of the complex structure consists of 4 atom types (C O N S). But I constantly get an error saying 'Atoms lost' (I have tried with different potentials like ReaxFF, LJ, etc...). I know this means a few atoms have been blown off so far from the simulation box that they have been lost. LAMMPS mailing list archives all suggest for a bad geometry, but there is nothing I can do to modify the data file. Hence I really need some help with correcting the code(seen below). The main motive of the code is for the Carbon nano cone to indent the protein. The code below is the code I have written for LJ potential.

# -------------- Initialisation -------------------#

units           metal
boundary        p p p
atom_style      atomic
timestep        0.001

read_data       SE_atomic.dat
mass    1 14.0067   #Nitrogen
mass    2 12.0107   #Carbon
mass    3 15.9940   #Oxygen
mass    4 32.0650   #Sulphur


#Substrate

region           sars block INF INF INF INF INF 300 units box
group            sars region sars

#Indenter

region           top block INF INF INF INF 300 INF units box
group            top region top

#Potentials


pair_style       lj/cut 2.5

pair_coeff  1  1  0.0163859696   2.9565
pair_coeff  1  2  0.07146236771  3.1365
pair_coeff  1  3  0.07168744659  2.9165
pair_coeff  1  4  0.02856991424  2.8232
pair_coeff  2  2  0.12219124354  3.2615
pair_coeff  2  3  0.12257609881  3.0415
pair_coeff  2  4  0.04885079323  2.9482
pair_coeff  3  3  0.07146236771  3.1365
pair_coeff  3  4  0.07168744659  2.9165
pair_coeff  4  4  0.09181993247  3.0815


compute         peratom all pe/atom
neighbor        2.0 bin
neigh_modify delay 20 every 1 check yes page 500000 one 50000

#Initialisation

velocity all create 353 511124 rot yes mom yes
variable        t_step equal "step"
variable        t_temp equal "temp"
thermo          100
thermo_style    custom step press temp pe

#Relaxation

fix rigid sars addforce 0.0 0.0 0.0
fix relax top nvt temp 250 250 0.01
variable t equal 0.005
timestep ${t}
dump 1 all custom 500 equil.*.dump id type x y z fx fy fz

run 3000


#Indentation

fix 1 top addforce 0.0 -0.05 0.0
fix fix1 all nvt temp 250 250 0.01
dump 2 all custom 100 sars-Indent.*.dump id type x y z fx fy fz

run 1000


print "All done"

This code is written for ReaxFF potential

# -------------- Initialisation -------------------#

units           metal
boundary        p p p
atom_style      charge
timestep        0.001

read_data       SE_charge.dat
mass    1 14.0067   #Nitrogen
mass    2 12.0107   #Carbon
mass    3 15.9940   #Oxygen
mass    4 32.0650   #Sulphur


#Substrate

region           sars block INF INF INF INF INF 300 units box
group            sars region sars

#Indenter

region           top block INF INF INF INF 300 INF units box
group            top region top

#Potentials

pair_style      reax/c NULL
pair_coeff      * * SARS.ff N C O S
compute         peratom all pe/atom
neighbor        0.5 bin
neigh_modify delay 20 every 1 check yes page 500000 one 50000

#Initialisation

velocity all create 353 511124 rot yes mom yes
variable        t_step equal "step"
variable        t_temp equal "temp"
thermo          100
thermo_style    custom step press temp pe

#Relaxation

fix rigid sars addforce 0.0 0.0 0.0
fix relax top nvt temp 250 250 0.01
fix charge all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
variable t equal 0.005
timestep ${t}
dump 1 all custom 500 equil.*.dump id type x y z fx fy fz

run 2000


#Indentation

fix 1 top addforce 0.0 -0.05 0.0
fix fix1 all nvt temp 250 250 0.01
dump 2 all custom 100 sars-Indent.*.dump id type x y z fx fy fz

run 1000


print "All done"

If anyone in the community can help regarding the same, it would be really helpful for me! Cheers!

I have uploaded my LAMMPS datafile for the LJ case and my datafile for the ReaxFF case, and all data files, input files, and log files are located in this folder.

Answer 1

The "Lost atoms" error typically happens when huge forces blow up an MD calculation. This is also the problem here:

The main issue are the units. The input parameters for reax force fields are given in "real" units by default. Unless they are converted real units have to be used in the lammps input file. Setting the units appropriately eliminates the "Lost atoms" error because the forces now have reasonable values.

Two more (not directly related) issue are standing out: The "fix nvt" and the "region" commands. Unless explicitly unset a fix persists until the end of the simulation. Therefore,

fix relax top nvt temp 250 250 0.01

and

fix fix1 all nvt temp 250 250 0.01

each apply a time integrator to all atoms in "top". This will certainly lead to undesired results.

The simulation box defined in "SE_charge.dat" only extends to 216, so the region top is always empty.

region           top block INF INF INF INF 300 INF units box

Therefore, the entire Relaxation block does nothing because the nvt time integrator is applies to 0 atoms.

TL;DR With a few tweaks to the input file the simulation with the reaxFF is running. All I changed is

# The units for the parameters in reaxFF input decks are real
units           real                 # time is given in fs

# Your simulation box only extends to 216.506012, setting 300 in the commands below leaves the region top empty. Set this value to something usefull
region           sars block INF INF INF INF INF 200 units box
region           top block INF INF INF INF 200 INF units box

# adjust the time step to the new units
fix relax top nvt temp 250 250 100
variable t equal 1

# use unfix to remove the nvt fix, otherwise you have two time integrators active!
unfix rigid
unfix relax

# again adjust the time step
fix fix1 all nvt temp 250 250 100

EDIT: The entire input file with changes is here.