Using lammps's compute heat flux

Dear All,

I am using lammps to confirm the result of another simulation. I need to create a temperature gradient in x direction and measure the heat flux. I have attached my lammps code. Am I correct in only measuring q over the two regions that I am thermostatting and then finding the difference? Lammps results seem to be an order of magnitude higher than what I get from my other code or Fourier’s law. I would really appreciate your help.

Regards,

Rambod

units real variable T equal 80
variable V equal vol
variable dt equal 1.0
variable p equal 1000
variable s equal 1
variable d equal $p*$s
variable NA equal 6.0221418*1e23
#–real unit conversion for viscosity–
variable kB equal 1.3806504e-23
variable atm2Pa equal 101325.0
variable A2m equal 1.0e-10
variable fs2s equal 1.0e-15
dimension 3

boundary p p p
neigh_modify check yes
neighbor 3.0 bin
lattice fcc 5.72 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region box block -20.0 20.0 -2.5 2.5 -2.5 2.5 units lattice
region left block -15.0 -13.0 INF INF INF INF units lattice
region right block 13.0 15.0 INF INF INF INF units lattice
region center block -1.25 1.25 INF INF INF INF units lattice
create_box 2 box

create_atoms 1 box

group ar type 1
group cu type 2
group hot region left
group cold region right
mass 1 39.948
mass 2 127.092
pair_style lj/cut 9.534
pair_modify mix arithmetic
pair_coeff 1 1 0.2381 3.405
pair_coeff 2 2 9.4500 2.3377

min_style cg

minimize 0.0 0.0 10000 10000
compute Thot all temp/region left
compute Tcold all temp/region right
fix NVT all nvt temp $T T 10 drag 0.2 timestep {dt}
thermo_style custom step temp c_Thot c_Tcold
thermo 1000
run 100000
velocity all scale $T
unfix NVT
reset_timestep 0
fix 1 all nve
fix hot_rescale1 all temp/rescale 1 $T 80 0.0001 1.0
fix cold_rescale1 all temp/rescale 1 T 80 0.0001 1.0 compute hot_temp1 all temp/region left compute cold_temp1 all temp/region right fix_modify hot_rescale1 temp hot_temp1 fix_modify cold_rescale1 temp cold_temp1 run 200000 compute myKE all ke/atom compute myPE all pe/atom compute myStress all stress/atom NULL virial compute fluxhot hot heat/flux myKE myPE myStress compute fluxcold cold heat/flux myKE myPE myStress variable Jxhot equal c_fluxhot[1] variable Jxcold equal c_fluxcold[1] variable fluxdiff equal v_Jxhot-v_Jxcold fix hot_rescale all temp/rescale 1 80 80 0.0001 1.0 fix cold_rescale all temp/rescale 1 80 80 0.0001 1.0 compute hot_temp all temp/region left compute cold_temp all temp/region right fix_modify hot_rescale temp hot_temp fix_modify cold_rescale temp cold_temp fix avediff all ave/time 1 1 1000 v_fluxdiff ave running start 10000 thermo_style custom step temp c_hot_temp c_cold_temp v_fluxdiff f_avediff fix fluxprint all print 1000 "{Jxhot} ${Jxcold}" file data.txt screen no
run 1000000