compute heat/flux metal units

Dear Lammps-users,
Sorry if this topic has been visited before, but its not clear to me what
I'm missing. I want to reproduce the results of the compute heat/flux
example script as described in the documentation, only using metal units
instead of real units.

The documentation gives the correct result as ~0.29 W/mK. The result I get
for the script in real units is 0.286 W/mK which is in agreement. When I
modify the script with metal units I get 0.229 W/mK which looks wrong.

Any advice would be greatly appreciated!
Thanks,
John

Modified Example Script for Metals:

# Sample LAMMPS input script for thermal conductivity of solid Ar

units metal
variable T equal 70
variable V equal vol
variable dt equal 0.004
variable p equal 200 # correlation length
variable s equal 10 # sample interval
variable d equal $p*$s # dump interval

# convert from LAMMPS real units to SI

variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable eV2J equal 1.60217646e-19
variable A2m equal 1.0e-10
variable ps2s equal 1.0e-12
variable convert equal \{eV2J\}\*{eV2J}/\{ps2s\}/{A2m}

# setup problem

dimension 3
boundary p p p
lattice fcc 5.376 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1

region box block 0 4 0 4 0 4
create_box 1 box
create_atoms 1 box
mass 1 39.948
pair_style lj/cut 13.0
pair_coeff * * 0.010325 3.405
timestep ${dt}
thermo $d

# equilibration and thermalization

velocity all create $T 102486 mom yes rot yes dist gaussian
fix NVT all nvt temp $T $T 10 drag 0.2
run 8000

# thermal conductivity calculation, switch to NVE if desired

unfix NVT
fix NVE all nve

reset_timestep 0
compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom virial
compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol
fix JJ all ave/correlate $s $p d &              c\_flux\[1\] c\_flux\[2\] c\_flux\[3\] type auto file J0Jt\.dat &              ave running variable scale equal {convert}/${kB}/$T/$T/$V*s\*{dt}
variable k11 equal trap(f_JJ[3])*\{scale\} variable k22 equal trap\(f\_JJ\[4\]\)\*{scale}
variable k33 equal trap(f_JJ[5])*${scale}
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33
run 100000
variable k equal (v_k11+v_k22+v_k33)/3.0
variable ndens equal count(all)/vol
print "average conductivity: $k[W/mK] @ T K, {ndens} /A^3"

There are many, many ways to screw this kind of conversion
up. If you do it right, then the simulation you run will
be identical (step by step, for any quantity you monitor)
to the original simulation. I suggest you start simple
and add new commands one at a time and see when
you mess up.

Steve