In NVE irradiation simulations, the temperature rises too much, and excessive defect generation occurs

LAMMPS LAMMPS General Discussion
1

Hi Everybody!
I’ve encountered a minor issue. I supplied 1 keV of energy to a carbon atom within a diamond crystal to observe the resulting damage in the crystal. Initially, I set up a small crystal model.

Theoretical predictions suggested that the number of defects would be around a few dozen (the displacement energy of a carbon atom in diamond is 45 eV). However, in the actual simulation, I observed the system’s temperature skyrocketing to the order of 10^8 K (thermo_style custom temp). Under these conditions, all atoms in the system dispersed. Eventually, the lattice structure ceased to exist, and the number of defects significantly exceeded expectations. What could be the cause of this?
thank you

this is my in.lmp

1. Initialization

units metal # Use metal units
dimension 3 # Dimension of simulation system
boundary p p p # Periodic boundary conditions

2. Atom/molecule definition

atom_style atomic # Type of atoms
lattice diamond 3.567 # FCC lattice parameter of diamond, using metal units
region box block 0 4 0 4 0 10 units lattice # Size of simulation box
create_box 1 box # Create box with 1 type of atoms
create_atoms 1 box # Create single atom in box
region rpka sphere 2 2 1 0.2 units lattice
group pka region rpka #pka

mass 1 12.011

3. Interaction potential

pair_style tersoff/zbl # Use Tersoff potential
pair_coeff * * SiC.tersoff.zbl C # Import appropriate Tersoff/ZBL potential parameter file

4. Set initial conditions

fix 1 all nvt temp 300 300 0.05 # Set NVT for relaxation
timestep 0.001
thermo 5 # Output the thermo information
thermo_style custom step temp etotal time
thermo_modify lost warn norm yes
run 20000
unfix 1 # cancel NVT

fix 1 all nve # Set NVE for radiation
#fix 1 all nvt temp 300.0 300.0 100.0
velocity all set 0.0 0.0 0.0 units box # Initial velocity of all atoms is zero
velocity pka set 0 0 1266.97
#fix 2 all dt/reset 1 0.00000001 0.0005 0.005 units lattice
timestep 0.001
thermo_modify lost ignore
thermo 5 # Output the thermo information
thermo_style custom step temp etotal time

dump 1 all atom 10 10keV__box50_10ps_1fs.dump # Record atomic positions

Run simulation for 10 steps

thermo_modify lost warn norm yes
run 20000

2

Two suggestions:

  • you may need to use a larger system to better dissipate the kinetic energy. You are putting a large amount of energy into a small system without any way of removing it.
  • you may need to use a shorter time step. You are accelerating a single atom by a very large amount of kinetic energy. This may break the numerical integration since the discretization must be adapted to the fastest motion in the system. If the time step is too large, the velocity Verlet integrator will diverge and the time integration will produce garbage. As the energy dissipates through collisions, the time step may be increased again. This is what fix dt/reset is for if you cannot afford to run with the short time step throughout.
1 Like
4

Hi akohlmey,

Thank you for your advice, it seems these two reasons are the main factors. I think the setting of the time step is more important. After using the option dt/reset NULL 0.0005 , it looks much better. Additionally, I hadn’t realized that the default unit for velocity is units lattice. I had always assumed it was Å/ps, which resulted in the velocities I set being slightly more than three times what I intended.