Langevin thermostat in thermal conduction problem

Dear lammps users,

I’m trying to use langevin thermostats as hot/cold reservoir to obtain thermal conducivity. I can obtain a value of thermal conductivity from the script. However, I find that the temperature value is not the one I set. (I set 330K for the hot end, but it fluctuates around 315K instead. Similarly for the cold end). I wonder if I can fix it at 330K by changing setting of langevin. My script is in the attachment.
Any suggestions will be appreciated

variable T equal 300

units metal
boundary f p f
atom_style full
read_data data.2G
neighbor 2.0 bin
neigh_modify delay 5

#Potential function

pair_style airebo 2.5
pair_coeff * * E:/lammps/CH.airebo C

mass 1 12.081

#groups & heat layers

group HOT molecule 1
group COLD molecule 2
group Leftbound molecule 3
group Rightbound molecule 4

group total subtract all Leftbound Rightbound

compute Thot HOT temp
compute Tcold COLD temp

compute Ttotal total temp
compute kinetic all ke
compute potential all pe

variable thot equal c_Thot
variable tcold equal c_Tcold
variable totalenergy equal c_kinetic+c_potential

timestep 0.001

fix LB Leftbound setforce 0 0 0
fix RB Rightbound setforce 0 0 0

#Initial condition: System is at 300K
velocity total create $T 798456 dist gaussian

#1st Equilibrium, NH
fix 1 total nvt temp $T $T 50

thermo_style custom step c_Tcold c_Thot c_Ttotal
thermo_modify lost warn
thermo 1000
dump trajectory all atom 5000 3c60tra.dump

run 400000

unfix 1

#2nd equilibrium nve
variable Th equal 330
variable Tc equal 270

fix 1 total nve

run 200000
write_restart restart.equil
fix hot HOT langevin {Th} {Th} 1 6453 tally yes
fix cold COLD langevin {Tc} {Tc} 1 12354 tally yes


run 500000
compute kine all ke/atom

fix spatialoutput total ave/spatial 1 20000 20000 x lower 0.025 c_kine units reduced &
file spatialT.dat
fix Delta_out all ave/time 1 1 10000 f_hot f_cold file fluxfix.dat

run 2000000

in.langevin (1.65 KB)

How many atoms in the hot/cold groups. Are those

atoms constrained in any way? If you run Thot = Tcold

do the 2 groups have the desired temperature?



Thank you Steve for replying. There are 432 atoms in hot reservoir, and the same for cold end. There are 72 fixed atoms near hot/cold end. I just tried that if Thot=Tcold=330K the system can hold a temperature at the prescribed value.



72 atoms is not a lot to maintain a temperature, esp if they

are losing/gaining heat with their neighbors. I’m guessing

if you increased the model size, your 2 temps would stay

closer to what you want.

Another issue is how you are measuring the temp once atoms

start moving. If some of the 72 atoms move closer into the

center, that is another complicating factor. You can use fix heat

with a region (at each end) and likewise measure the temp of

the atoms in the same regions (see compute temp/region),

which might be more accurate than what you are doing.




Thanks for your advice