Problem with lost atom

Dear All,
I am simulating the thermal conductance along the z axis for bilayer graphene sandwiched by Al2O3 system. This is my input setting.

variable x equal 47.59
variable y equal 41.214148966101
variable z equal 136.636
variable t equal 300.0

#Setup parameters
units real
atom_style atomic
read_data atomic_structure_atomic
mass 1 16.0
mass 2 27.0
mass 3 12.0
mass 4 12.0
velocity all create $t 87287

pair_style hybrid/overlay lj/cut 12.0 morse 6.287 tersoff shift -0.0462
pair_coeff 1 1 lj/cut 0.194633996 3.541
pair_coeff 1 2 lj/cut 0.31351524159 4.02
pair_coeff 1 3 lj/cut 0.0 1
pair_coeff 1 4 lj/cut 0.0 1
pair_coeff 2 2 lj/cut 0.50500834301 4.499
pair_coeff 2 3 morse 10.81786819 1.738 2.246
pair_coeff 2 4 morse 10.81786819 1.738 2.246
pair_coeff * * tersoff ./Copt.tersoff NULL NULL C C

neighbor 0.3 bin
neigh_modify delay 0 every 1

#layers for heat flux
region hot block INF INF INF INF 28.69356 41.810616
region cold block INF INF INF INF 96.32838 97.01156
compute Thot all temp/region hot
compute Tcold all temp/region cold

#1st equilibirum run
fix 1 all nvt temp $t $t 100.0
thermo 1
run 3000000
velocity all scale $t
unfix 1

#2nd equilibrium run
fix 1 all nve
fix hot all heat 1 23.0609 region hot
fix cold all heat 1 -23.0609 region cold
thermo_style custom step temp c_Thot c_Tcold
thermo_modify colname c_Thot Temp_hot colname c_Tcold Temp_cold
thermo 10
run 10000

#thermal conductivity calculation
compute ke all ke/atom
variable temp atom c_ke/23.0609/1.5/8.617333262145e-5
compute layers all chunk/atom bin/1d z lower 1.0 units reduced
fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.heat
variable tdiff equal f_2[20][3]-f_2[36][3]
fix ave all ave/time 1 1 1000 v_tdiff ave running start 13000

variable kappa equal 23.0609/(lx*ly)*lz/f_ave
thermo_style custom step temp c_Thot c_Tcold v_tdiff f_ave
thermo_modify colname c_Thot Temp_hot colname c_Tcold Temp_cold &
              colname v_tdiff dTemp_step colname f_ave dTemp

run 3000000
print "Running average thermal conductivity: $(v_kappa:%.4f)"

This is the error message I received after I launched the calculation.

LAMMPS (2 Aug 2023 - Update 1)
variable x equal 47.59
variable y equal 41.214148966101
variable z equal 136.636
variable t equal 300.0

#Setup parameters
units real
atom_style atomic
read_data atomic_structure_atomic
Reading data file ...
  triclinic box = (0 0 0) to (47.59 41.214149 136.636) with tilt (-23.795 0 0)
WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (../domain.cpp:220)
  2 by 2 by 8 MPI processor grid
  reading atoms ...
  19600 atoms
  read_data CPU = 0.131 seconds
mass 1 16.0
mass 2 27.0
mass 3 12.0
mass 4 12.0
velocity all create $t 87287
velocity all create 300 87287

#pair_style hybrid/overlay lj/cut 10.0 morse 6.287 tersoff shift -0.0462 kolmogorov/crespi/z 20.0
pair_style hybrid/overlay lj/cut 12.0 morse 6.287 tersoff shift -0.0462
pair_coeff 1 1 lj/cut 0.194633996 3.541
pair_coeff 1 2 lj/cut 0.31351524159 4.02
pair_coeff 1 3 lj/cut 0.0 1
pair_coeff 1 4 lj/cut 0.0 1
pair_coeff 2 2 lj/cut 0.50500834301 4.499
pair_coeff 2 3 morse 10.81786819 1.738 2.246
pair_coeff 2 4 morse 10.81786819 1.738 2.246
pair_coeff * * tersoff ./Copt.tersoff NULL NULL C C

neighbor 0.3 bin
neigh_modify delay 0 every 1

#layers for heat flux
region hot block INF INF INF INF 28.69356 41.810616
region cold block INF INF INF INF 96.32838 97.01156
compute Thot all temp/region hot
compute Tcold all temp/region cold

#1st equilibirum run
fix 1 all nvt temp $t $t 100.0
fix 1 all nvt temp 300 $t 100.0
fix 1 all nvt temp 300 300 100.0
thermo 1
run 3000000
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12.3
  ghost atom cutoff = 12.3
  binsize = 6.15, bins = 12 7 23
  6 neighbor lists, perpetual/occasional/extra = 6 0 0
  (1) pair lj/cut, perpetual, skip from (4)
      attributes: half, newton on
      pair build: skip
      stencil: none
      bin: none
  (2) pair morse, perpetual, skip from (5)
      attributes: half, newton on, cut 6.587
      pair build: skip
      stencil: none
      bin: none
  (3) pair tersoff, perpetual, skip from (6)
      attributes: full, newton on, cut 2.4
      pair build: skip
      stencil: none
      bin: none
  (4) neighbor class addition, perpetual
      attributes: half, newton on
      pair build: half/bin/newton/tri
      stencil: half/bin/3d/tri
      bin: standard
  (5) neighbor class addition, perpetual, trim from (4)
      attributes: half, newton on, cut 6.587
      pair build: trim
      stencil: none
      bin: none
  (6) neighbor class addition, perpetual
      attributes: full, newton on, cut 2.4
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.496 | 8.707 | 10.81 Mbytes
   Step          Temp          E_pair         E_mol          TotEng         Press     
         0   300            4.3067738e+08  0              4.306949e+08   4.4405907e+08
         1   9.6397316e-73  2.5266396e+09  0              2.5266396e+09  2.5914761e+09
ERROR: Lost atoms: original 19600 current 19595 (../thermo.cpp:488)
Last command: run 3000000

I checked the unit for the potential and energy settings. It is said that the metal units would be converted into real unit by LAMMPS for the tersoff potential automatically. I did not changed the tersoff potential and only change the unit for energy from eV to Kcal/mol in the lj and mose potential.

I think that the force was not that large to blow up some atoms out of the unit cell. I do not understand why I still receive this error message.

Would anyone please give me some suggestions/comments on the parameter settings in my input file?

Thank you in advance.

I think differently. You have rather large energies and a pressure of 4.4405907e+08. Such a large pressure requires extremely large forces.

That will only happen, if the tersoff potential file contains the suitable metadata to identify which units they data is stored. If that metadata tag is missing, there is no conversion. Whether parameters are converted will be logged in the log file.

Thank you for the explanation.

I think that I will switch to metal unit.

As an interesting heuristic, one mole of TNT weighs 227 grams, and its explosion releases about 227 kcal of energy, so that the molecular free energy density of TNT is about 227 kcal/mol over 21 atoms.

A system of about twenty thousand particles with a potential energy of 430 million kcal/mol would be … twitchy.

1 Like

Thank you for the reply and hint.

I like the way you use to estimate the problem. It is very inspiring.