Effect of Temperature and Electric Field

Hello.
I am simulating a system where I have deposited H atoms on Si surface. I am trying to see the effect of temperature and electric field on H atoms and dumped Kinectic Energy of H. When there’s no Electric field in the system, H shows random hopping and KE of H shows random values due to Temperature (I have Vx, Vy, and Vz values). Again, when I am applying Efield, I observed H hopped on Si atoms making bonds (bonds observed from bond record files and OVITO) in the direction of Efield. However, the KE of H fall down to nearly zero and increased over time, which represents the characteristics of free particles (I have only Vx values, Vy, Vz are zero).

  1. If H is making bonds while hopping, why it’s KE and Velocity is always increasing? Why Vy, Vz becomes 0? Why is there no effect of temperature?

  2. As per my understanding, the KE or Velocity should fall at the time of bonding. Is there any explanation from LAMMPS perspective?

  3. When Efield is applied, it seems there is no effect of temperature. Any possible way to see both the effects of temp and Efield?

Following is the code:

units real
dimension 3
boundary f f f
atom_style full

lattice diamond 5.4307
region box block 0 13 0 7 0 7
create_box 2 box

region silicon block 0 13 0 7 0 2
group silicon region silicon

create_atoms 1 region silicon

group hydrogen type 2

pair_style reax/c NULL
pair_coeff * * ffield.txt Si H

mass 1 28.0855 # Si
mass 2 1.008 # H

neighbor 2.0 bin
neigh_modify every 1 delay 0 check no

timestep 0.05

fix 1 all nve
fix 2 all langevin 380 380 2 587283

region slab block 0 13 0 7 6 7
group slab region slab

fix 5 hydrogen deposit 1 2 100 17895 region slab vz -2.2 -2.2
fix 6 hydrogen deposit 1 2 100 14938 region slab vz -2.2 -2.2
fix 7 hydrogen deposit 1 2 100 83659 region slab vz -2.2 -2.2
fix 8 hydrogen deposit 1 2 100 54320 region slab vz -2.2 -2.2
fix 9 hydrogen deposit 1 2 100 19283 region slab vz -2.4 -2.4
fix 10 hydrogen deposit 1 2 100 30045 region slab vz -2.2 -2.2
fix 11 hydrogen deposit 1 2 100 23190 region slab vz -2.1 -2.1

compute eng all pe/atom
compute eatoms all reduce sum c_eng
compute reax all pair reax/c
variable eb equal c_reax[1]

fix parameter all qeq/reax 1 0.0 12.0 1e-6 param.qeq16

thermo_style custom step atoms temp c_eatoms v_eb
thermo 100
thermo_modify lost warn

fix forty all reax/c/bonds 100 x_bondrec.txt #pg 1066
fix fifty all reax/c/species 1 2 100 deposition.out element Si H cutoff 1 2 0.3

compute ke_all all ke/atom

dump 2 all custom 20 velocity.lammpstrj id type q x y z vx vy vz c_ke_all c_eng

min_style cg
minimize 1e-25 1e-25 5000 10000

run 500

###If i don’t use these commands of velocity and setforce, Efiled seems not working
velocity silicon set 0 0 0
velocity hydrogen set 0 0 0
fix 15 silicon setforce 0.0 0.0 0.0
fix 16 hydrogen setforce 0.0 0.0 0.0

#---------Electri Field-‘X’-----------------------------

fix electric all efield 0.00003 0.0 0.0
run 1400000

variable natoms equal “count(all)”
variable teng equal “c_eatoms”
variable length equal “lx”
variable ecoh equal “v_teng/v_natoms”