# Inconsisten kinetic energy calculation

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
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.
Please also note that the allocated data will be freed during a `lammps.close()` call, so the print statements accessing the extracted data must come before the lmp.close() or else you are likely to get a segmentation fault.