I am trying to recreate the results from this paper. Specifically, the nanowire I am testing is cylindrical with a cross-sectional area of 24 nm2 and a length of 10 nm. I am expecting a thermal conductivity of around 1.9 W/mK, however, I am only getting around 1.0 W/mK. The paper used Müller-Plathe’s formalism.
I appropriated the Müller-Plathe sample script (KAPPA/in.mp). I am not sure if my code is proper.
# Settings
units metal
atom_style full
# Variables
variable V equal vol
variable dt equal 0.0005
variable p equal 200 # correlations length
variable s equal 10 # sample interval
variable d equal $p*$s # dump interval
variable seed equal 190511
## Conversion
variable kB equal 1.3806504e-23 # [J/K] Boltzmann
variable eV2J equal 1.6021e-19
variable A2m equal 1.0e-10
variable nm2A equal 10
variable ps2s equal 1.0e-12
variable convert equal ${eV2J}*${eV2J}/${ps2s}/${A2m}
# Configuration
## NW settings
variable area equal 24 # [nm]
variable nwrad equal sqrt(${area}/PI)*${nm2A}
variable nwlen equal 10 # [nm]
## Ge Ratio: 0.1-0.9
variable ratio equal 0.0
## Temperature: 100K-1100K
variable T equal 1100
## Simbox
variable simxl equal ${nwrad}*6 # simbox x length
variable simyl equal ${nwrad}*6 # simbox y length
variable simzl equal ${nwlen}*${nm2A} # simbox z length
variable simx1 equal -${simxl}/2
variable simx2 equal ${simxl}/2
variable simy1 equal -${simyl}/2
variable simy2 equal ${simyl}/2
variable simz1 equal -${simzl}/2
variable simz2 equal ${simzl}/2
# Problem Setup
dimension 3
boundary f f p
lattice diamond 5.43
region simbox &
block ${simx1} ${simx2} &
${simy1} ${simy2} &
${simz1} ${simz2} &
units box
create_box 2 simbox
region nanowire &
cylinder z 0 0 ${nwrad} INF INF &
units box
create_atoms 1 region nanowire
set region nanowire &
type/ratio 2 ${ratio} ${seed}
mass 1 28.0855 # Silicon
mass 2 72.6400 # Germanium
velocity all create $T ${seed} &
mom yes &
rot yes &
dist gaussian
pair_style tersoff
pair_coeff * * SiCGe.tersoff Si(D) Ge
neighbor 0.3 bin
neigh_modify delay 0 every 1
# 1st Equilibriation Run
timestep ${dt}
fix 1 all nvt temp $T $T 0.5
thermo 100
run 100000
dump 1 all custom 100 ./output/Ge0.0_T1100.lammpstrj id type mass xs ys zs fx fy fz
velocity all scale $T
unfix 1
# 2nd equilibration run
compute ke all ke/atom
variable temp atom c_ke/1.5
fix 1 all nve
compute layers all chunk/atom bin/1d z lower 0.05 units reduced
fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.mp
fix 3 all thermal/conductivity 10 z 20
variable tdiff equal f_2[11][3]-f_2[1][3]
thermo_style custom step temp epair etotal f_3 v_tdiff
thermo_modify colname f_3 E_delta colname v_tdiff dTemp_step
thermo 1000
run 100000
dump 2 all custom 100 ./output/la_Ge0.0_T1100.lammpstrj id type mass xs ys zs fx fy fz
# thermal conductivity calculation
# reset fix thermal/conductivity to zero energy accumulation
fix 3 all thermal/conductivity 10 z 20
variable start_time equal time
variable kappa equal (f_3/(time-${start_time})/(lx*ly)/2.0)*(lz/2.0)/f_ave
fix ave all ave/time 1 1 1000 v_tdiff ave running
thermo_style custom step temp epair etotal f_3 v_tdiff f_ave
thermo_modify colname f_3 E_delta colname v_tdiff dTemp_step colname f_ave dTemp
run 100000
print 'Running average thermal conductivity: $(v_kappa:%.2f)'
Are there other settings I can fine-tune (e.g., temperature difference between the cold and hot slabs, number of slabs, thickness of slabs) that might achieve the desirable values?
