Cooling using fix langevin

Dear all
I am trying to perform thermal quasistatic shear deformation over a 2D glass system. To do so, I refer to the following paper: Avalanches in the Athermal Quasistatic Limit of Sheared Amorphous.
Solids: An Atomistic Perspective(https://doi.org/10.1007/s11249-021-01439-5). I am trying to validate the results of the paper mentioned.
I am using July 2021, version of LAMMPS.
To prepare the glass, I am using the following code.
I am equilibrating the system on temperature 0.65, using langevin thermostat as mentioned in the paper. Then, the system was cooled to temperature T = 0 at a rate of 0.002.

# Big colloid particles and small LJ particles
#../../src/lmp_mpi -i in.colloid

units		lj   # lj or real or metal or si or cgs or electron or micro or nano (unitless
atom_style	sphere
boundary p p p 
dimension	2

region      box prism 0.0 53.0 0.0 53.0 -0.5 0.5 0.0 0.0 0.0 


create_box	2 box
create_atoms 1 random 1282 3000 box units box
create_atoms 2 random 1584  4000 box units box

set type 1 mass 1.0 
set type 2 mass 1.0
set type 1 diameter 1.175
set type 2 diameter 0.618

comm_modify vel yes

pair_style lj/cut 2.5
pair_coeff 1 1 0.5 1.175
pair_coeff 2 2 0.5 0.618
pair_coeff 1 2 1.0 1.0
 
neigh_modify every 1 delay 0 check yes
timestep 0.0001
minimize 0.0 1.0e-20 100000 1000000

fix bd all langevin 0.65 0.65 1.0 15111
fix 20 all nve/sphere
fix 5 all enforce2d
thermo 10000
dump    mydump all atom 10000 glass.dump
dump             1 all custom 10000 *.data_equilibration id type diameter mass x y z vx vy fx fy
dump_modify 1 sort id
restart 10000 restart.*.glass_0.65
run 2000000
unfix bd
undump mydump
undump 1 
minimize 0.0 0.0 100000 1000000 
thermo 1000
fix bd all langevin 0.65 0.0 1.0 15111
dump    mydump all atom 10000 glass_cool.dump
dump             1 all custom 10000 *.data_cool id type diameter mass x y z vx vy fx fy
dump_modify 1 sort id
restart 10000 restart.*.glass_cool
run 3250000

Now while cooling the system from temperature 0.65 to 0.0 , temperature of the system starts to increase , ultimately leading to the loss of particles in the system.
I am getting the following error:

fix bd all langevin 0.65 0.0 1.0 15111
dump    mydump all atom 10000 glass_cool.dump
dump             1 all custom 10000 *.data_cool id type diameter mass x y z vx vy fx fy
dump_modify 1 sort id
restart 10000 restart.*.glass_cool
run 32500000
Per MPI rank memory allocation (min/avg/max) = 4.331 | 4.333 | 4.335 Mbytes
Step Temp E_pair E_mol TotEng Press
 2001390   0.63792997   -2.4150167            0   -1.7773093    1.6039234
 2002000   0.34016477   -1.5372197            0   -1.1971736    3.3175798
 2003000   0.45384006   -1.5060123            0   -1.0523306    3.1868157
ERROR: Lost atoms: original 2866 current 2301 (../thermo.cpp:439)
Last command: run 32500000

Can you please help me to resolve this issue. Any suggestion will be helpful.

Your particles have rotational degrees of freedom (atom style sphere), but I see no interactions that affect them. What is the point of this?

Have you checked your center of mass momentum?
or visualized your trajectory?
You may be victim of the so-called “flying ice cube effect”.

1 Like