I am not sure if this is a bug or I am missing something.
I am using LAMMPS to simulate collisional cascades. In such simulations, there is a box of atoms, which is thermalised first, and then one atom called the primary knock-on atom (PKA) is boosted with some kinetic energy usually given in eV or keV.
This is the basic LAMMPS input I use for a Fe PKA of 1 keV:
# mpiexec -np 8 lmp -in .\fecrhe.in
# 1 keV > 587.92 A/ps
units metal
atom_style atomic
boundary p p p
atom_modify map hash
# # Initialisation --------------------------------------------------------------
lattice bcc 2.8650
region R1 block 0 40 0 40 0 40
region R2 block 2 38 2 38 2 38
create_box 1 R1 # 1 atom types
create_atoms 1 box # Fe atoms
group Region1 region R1
group Region2 region R2
group thermoboundary subtract Region1 Region2
# # Load restart ---------------------------------------------------------------
# read_restart initial.restart
# Calculations ---------------------------------------------------------------
fix integrate all nve
fix reset all dt/reset 10 0.00001 0.001 0.026 units box
compute atom_pe all pe/atom
compute atom_ke all ke/atom
compute atom_st all stress/atom NULL
pair_style eam/alloy
pair_coeff * * fe_mendelev_eam_alloy.potential Fe
# Save -------------------------------------------------------------------------
write_restart initial.restart
# Saving data ----------------------------------------------------------------
dump simulation all custom 100 simulation.xyz id type x y z c_atom_pe c_atom_ke c_atom_st[1] c_atom_st[2] c_atom_st[3]
thermo_style custom step time temp ke pe etotal pxx pyy pzz press lx ly lz
thermo 200
# Thermalisation --------------------------------------------------------------
velocity all create 200.0 376847 loop geom # 100 K
run 1000
# Collisional cascade ---------------------------------------------------------
group PKA id 112001
velocity PKA set 587.92 0.0 0.0
fix control thermoboundary langevin 100.0 100.0 1.0 699483
run 1000
And this is the energy evolution (kinetic energy and potential energy minus initial potential energy) from log.lammps:
You can see that the thermalisation takes place correctly, but when the PKA is boosted, the energy increases around 8 keV.
Now if I execute the same code but from the restart file (comment initialisation and uncomment read_restart initial.restart), I get:
The jump in the kinetic energy is 1 keV as expected.
Why are these results different? Maybe it is related to units?