Hello, everyone.

I have been using LAMMPS these last months for my bachelor’s degree final project about particle irradiation on tungsten atoms. For that purpose, I use the interatomic potential designed by Bonny and his colleagues in its EAM1 variant.

At one point, I wanted to calculate the kinetic energy needed to displace an atom, but my results happened to be ten times lower than expected. After some testing, I found that LAMMPS energy calculations happened to be ten times larger the input.

I will give a brief explanation of my code before copying it here. In order to calculate the displacement energy in several directions, I used a python code that, given a Kinetic energy in (eV) and a direction (angles in spherical coordinates), it calculate the velocity in each cartesian direction. Then I execute a LAMMPS code where I set that velocity to a single atom (PKA).

Since I am using the “units metal” command, I need to do the following conversion:

```
E(J)=E(eV)*1.602e-19J/eV
m(kg)=183.84uma*1.660e-27kg/uma #atomic mass of W
v(m/s) = sqrt(2*E/m)
v(angstrom/ps)=v(m/s)/100
```

With this velocity, I treat it as the module of my velocity vector.

Now, my python code is as follows:

```
#Imports//////////////////
from lammps import lammps
import numpy as np
#spherical coordinates ////////////////////
def esfericas (E,theta,fi):
E=E*1.602e-19#de ev a J
m=3.0527348e-25#masa W de uma a Kg
v=np.sqrt(2*E/m)#de J a m/s
v=v/100
vx=v*np.sin(np.radians(theta))*np.cos(np.radians(fi))
vy=v*np.sin(np.radians(theta))*np.sin(np.radians(fi))
vz=v*np.cos(np.radians(theta))
return vx, vy, vz
#angles corresponding to the main directions I want to study
th100,fi100=90,0
th110,fi110=90,45
th111,fi111=54.74,45
#Function to extract the LAMMPS kinetic energy
def desp_obt (e,theta,fi):
vx,vy,vz=esfericas(e,theta,fi)
lmp=lammps()
lmp.command("variable vx equal %.3f"%vx)
lmp.command("variable vy equal %.3f"%vy)
lmp.command("variable vz equal %.3f"%vz)
lmp.file("in.desplazamientoprueba")
ke_exp=lmp.numpy.extract_compute("KE",0,0)
ke_atom=lmp.numpy.extract_compute("KEATOM",1,1)
print("Velocity according to:\nMe LAMMPS")
print("vx",vx,lmp.numpy.extract_variable("vx",vartype=0))
print("vy",vy,lmp.numpy.extract_variable("vy",vartype=0))
print("vz",vz,lmp.numpy.extract_variable("vz",vartype=0))
lmp.close()
print("Kinetic energy")
print("My calculation:",e)
print("LAMMPS (compute ke):",ke_exp)
print("LAMMPS (compute ke/atom):",np.unique(ke_atom))
desp_obt(20,th111,fi100)
```

My input script is as follows:

```
clear
#--------
#variables
#--------
variable a0 equal 3.13400196 #W lattice parameter
variable t_i equal 300 #temperature
variable dim equal 5 #System size
# ---------------
# settings
# ---------------
dimension 3
boundary p p p
atom_style atomic
neighbor 0 bin
neigh_modify delay 5
units metal
timestep 0.001
# ---------------
# create simulation box
# ---------------
lattice bcc ${a0}
region bigcube block -${dim} ${dim} -${dim} ${dim} -${dim} ${dim}
region box block -${dim} ${dim} -${dim} ${dim} -${dim} ${dim}
create_box 1 bigcube
mass 1 183.84
create_atoms 1 region box
# ---------------
# Potential
# ---------------
pair_style eam/alloy
pair_coeff * * Bonny_WHHe.txt W
# ---------------
# groups
# ---------------
group celda region box
group celda type 1
region centro sphere 0 0 0 0.1
group proyectil region centro
compute tcen all temp
velocity celda create ${t_i} 123456 temp tcen
#Minimization
fix boxrelax all box/relax iso 0.0 vmax 0.001
minimize 1e-18 1e-4 10000 10000
unfix boxrelax
#First NPT relaxation
fix 1 all npt temp ${t_i} ${t_i} 0.01 iso 0.1 0.1 1
run 100
unfix 1
#NVE relaxation
fix nvenergy all nve
fix lang all langevin ${t_i} ${t_i} 0.01 12345
run 100
#Kinetic energy calculation and velocity set
compute KE proyectil ke
compute KEATOM proyectil ke/atom
velocity proyectil set ${vx} ${vy} ${vz}
run 0
```

For what I read in the documetation, the calculation of the kinetic energy the summation of 1/2 mv^2, but given a single atom in the group I should obtain the same value as I calculated and I do not know where my calculation could have failed (I cheked several times).

Any insight would be really appreciated.

Thank you for reading.