LAMMPS input script for thermal conductivity of bismuth telluride
# thermostatting 2 regions via fix langevin
# LAMMPS input script for thermal conductivity of bismuth telluride
# thermostatting 2 regions via fix langevin
# Settings
dimension 3
boundary p p p
units metal
atom_style atomic
# Variables
variable kB equal 0.00008617
variable t equal 10 #150k=14.79360; 300k=5.00961; 450k=1.72856; 600k=0.28392; 750k=-0.82723 #900k=0.23078
variable tlo equal 9
variable thi equal 11
# Setup problem
read_data Bi2Te36x6x1.lammps-data
mass 1 208.98
mass 2 127.60
velocity all create $t 12312
# Interatomic potential
pair_style morse 5.5
pair_coeff * * 0.975 1.285 3.089
pair_coeff * * 0.085 2.212 4.203
pair_coeff * * 0.076 1.675 3.642
# Adjusted timestep value
timestep 0.000000005
# Heat layers
region hot block 0 1 0 1 7 9
region cold block 5 6 5 6 23 25
compute Thot all temp/region hot
compute Tcold all temp/region cold
# 1st equilibration run
fix 1 all nvt temp $t $t 0.1
thermo 10
run 100
velocity all scale $t
unfix 1
# 2nd equilibration run
fix 1 all nve
fix hot all langevin ${thi} ${thi} 0.1 45623 tally yes
fix cold all langevin ${tlo} ${tlo} 0.1 78934 tally yes
fix_modify hot temp Thot
fix_modify cold temp Tcold
variable tdiff equal c_Thot-c_Tcold
thermo_style custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff
thermo 1000
run 10000
# Thermal conductivity calculation
# Reset Langevin thermostats to zero energy accumulation
compute ke all ke/atom
variable temp atom c_ke/1.5
fix hot all langevin ${thi} ${thi} 0.1 45623 tally yes
fix cold all langevin ${tlo} ${tlo} 0.1 78934 tally yes
fix_modify hot temp Thot
fix_modify cold temp Tcold
fix ave all ave/time 10 100 1000 v_tdiff ave running
thermo_style custom step temp c_Thot c_Tcold f_hot f_cold v_tdiff f_ave
thermo_modify lost ignore flush yes
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.langevin
variable start_time equal time
# Calculate average value of f_ave
compute ave_f_ave all reduce ave f_ave
# Calculate thermal conductivity with divide-by-zero handling
variable kappa equal (0.5*(abs(f_hot)+abs(f_cold)))/((time-${start_time}) != 0 ? (time-${start_time})*(lx*ly) : 1e-30)*(lz/c_ave_f_ave)
run 20000
print "Running average thermal conductivity metal unit style: $(v_kappa:%.5f)"
print "Running average thermal conductivity W/mK: $(v_kappa*1.60218e+3:%.5f)"
The timestep is way too small. I have no direct experience of simulating inorganic solids, but something in the range 1e-3 to 5e-3 seems more reasonable.
I am not sure if the fix_modify is enough to achieve thermalisation only on the atoms inside each region:
I would apply the fix langevin not to the group all, but only to the atoms inside these groups:
group Ahot region hot
group Acold region cold
fix hot Ahot langevin ${thi} ${thi} 0.1 45623 tally yes
fix cold Acold langevin ${tlo} ${tlo} 0.1 78934 tally yes
and keep the rest the same. I assume your sample is solid, and no dynamic option is needed for the group.
It would also be helpful if you could elaborate on what " wrong " means. If the syntax is fine, no LAMMPS script is inherently wrong: it could implement an incorrect physical model, but that’s up to you to explain.
Yes, it is. Please look in the manual for applying a bias to a thermostat. 8.2.4. Thermostats — LAMMPS documentation
It is preferred over using a group since it is dynamic.
This seems quite suspicious to me. The region command defaults to using lattice units but since you are loading data instead of defining a lattice and creating atoms, I think it switches to using box units. Are you intending for your thermostat region to be quite so small?
I think there might be more issues with this but I’m not so good at visualizing these regions in my head so I am not confident in my interpretation. Have you visualized your simulation and highlighted the thermostat regions to confirm they are where you expect them to be?
And to echo @hothello, knowing what you think should happen / what is happening would likely increase the quality of advice you receive over just “it’s wrong”.
I haven’t had my coffee yet, so forgive me if I am misinterpreting your input.
It doesn’t switch, but without a lattice command, the lattice spacings are set by default to 1.0 1.0 1.0, so that there is no difference between lattice and box units.