[lammps-users] Temperature not stabilizing in thermal conductivity simulation

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,

I am new to lammps and I am simply trying for the moment to calculate the thermal conductivity of crystalline silicon at 300K using the Muller-Plathe method (input script below).
I use a VELOCITY CREATE to set initial velocities and then FIX NVT to establish the temperature. Then I use FIX NVE for a few timesteps before starting the Muller plate heat flux.
I have two issues:

  1. Temperature is divided by 2 between the first and second step before it slowly goes back up to the NVT thermostat set temperature (300K). I do not understand why? Should I set my initial velocities differently?
  2. Once temperature gets back to 300K, it keeps oscillating between 285 and 315K at least, no matter how long I wait for. And then I never get a clean linear temperature profile once I impose the heat flux.

What is it that I am doing wrong? Any piece of advice or information would be greatly appreciated.

Thank you in advance

3d bulk silicon

#Initialization############################################################

units metal
dimension 3
boundary p p p
atom_style atomic

Atom definition##########################################################

lattice diamond 5.43
region simbox block 0 10 0 4 0 4 units lattice

create_box 1 simbox
create_atoms 1 box
mass 1 28.0855

Atoms interactions settings##################################

pair_style sw
pair_coeff * * si.sw Si
#kspace_style pppm 1.0e-4

neighbor 2.0 bin
neigh_modify every 20 delay 0 check no

Equilibrate temperature in NVT

fix 1 all nvt 300.0 300.0 100.0
velocity all create 300.0 5674845

timestep 0.001
thermo_style custom step temp press &
ke vol lx ly lz epair emol etotal
thermo 1000
dump 1 all atom 1000 dump1.silicon
run 1000000

undump 1
unfix 1

#equilibrate in NVE mode ####################################

fix 3 all nve
thermo 1000
dump 2 all atom 1000 dump2.silicon
run 50000

undump 2
reset_timestep 0

#impose heat flux #########################

fix 4 all thermal/conductivity 180 x 20 swap 1

Output temperature profile##############################################

compute ke all ke/atom
variable temp atom c_ke[]/(1.5*8.617343e-5)
fix 5 all ave/spatial 1 1000 1000 x lower 0.5 v_temp &
file tmp.profile units lattice

Output heat flux value##################################################

thermo_style custom step temp press f_4 &
ke vol lx ly lz epair emol etotal

timestep 0.001
dump 3 all atom 1000 dump3.silicon
thermo 900
run 100000

Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
thomas.coquil@…24…

2

There is some value for the total energy of your system where the kinetic part will come out to be 300K when running NVE. I bet the potential energy for the initial configuration is what is off. Nose-Hoover will make the temperature fluctuate for a very long time. You can damp, or simply skip NVT and find the right initial velocities that give you the right long time temperature in NVE. The temperature profile will take a while to converge.
Matt

Quoting Thomas Coquil <[email protected]...>:

3

Thomas, if you start your lattice in an equilibrium configuration, your potential energy will be lower than a thermalized system. You will obtain faster convergence if you start your velocities with 2X the desired temperature which means the system will have the correct total energy. As far as fluctuations, that is what you should expect from MD.
Jeremy

4

Hello,

Thank you for your quick help.
So is there a way that would be better for me to stabilize my temperature (no NVT? velocity + NVE only? whatever?..)?
How am I ever supposed to get a meaningful temperature profile with Muller Plathe if my overall temperature and that of each individual slice keeps fluctuating so much (290 to 310K at least).
Note that I output the thermo data and the temperature of the slices every 1000 steps (1fs steps).
I have already tried waiting for a very long time (10M steps = much longer than people usually do).

Thank you very much

input script:

3d bulk silicon

#Initialization############################################################

units metal
dimension 3
boundary p p p
atom_style atomic

Atom definition##########################################################

lattice diamond 5.43
region simbox block 0 10 0 4 0 4 units lattice

create_box 1 simbox
create_atoms 1 box
mass 1 28.0855

Atoms interactions settings##################################

pair_style sw
pair_coeff * * si.sw Si
#kspace_style pppm 1.0e-4

neighbor 2.0 bin
neigh_modify every 20 delay 0 check no

Equilibrate temperature in NVT

fix 1 all nvt 300.0 300.0 100.0
velocity all create 300.0 5674845

timestep 0.001
thermo_style custom step temp press &
ke vol lx ly lz epair emol etotal
thermo 1000
dump 1 all atom 1000 dump1.silicon
run 1000000

undump 1
unfix 1

#equilibrate in NVE mode ####################################

fix 3 all nve
thermo 1000
dump 2 all atom 1000 dump2.silicon
run 50000

undump 2
reset_timestep 0

#impose heat flux #########################

fix 4 all thermal/conductivity 180 x 20 swap 1

Output temperature profile##############################################

compute ke all ke/atom
variable temp atom c_ke[]/(1.5*8.617343e-5)
fix 5 all ave/spatial 1 1000 1000 x lower 0.5 v_temp &
file tmp.profile units lattice

Output heat flux value##################################################

thermo_style custom step temp press f_4 &
ke vol lx ly lz epair emol etotal

timestep 0.001
dump 3 all atom 1000 dump3.silicon
thermo 900
run 10000000

Thomas Coquil
Research Assistant
UCLA Mechanical and Aerospace Engineering Department
420 Westwood Plaza
43-132 ENGINEERING IV
Los Angeles CA, 90095
Tel: 1 626 429 2110
thomas.coquil@…24…

5

Quoting Thomas Coquil <[email protected]...>:

Hello,

Thank you for your quick help.
So is there a way that would be better for me to stabilize my temperature
(no NVT? velocity + NVE only? whatever?...)?

If you run NVE at some total energy (kinetic+potential), there will be some constant temperature. That is, some constant *average* temperature - assuming you conserve energy. There will be fluctuations. Sometimes very large fluctuations depending on your system size, potential, and so forth. But the long time average will be constant. My guess is that you are starting from a poorly chosen configuration, i.e. a far from equilibrium potential energy for your desired temperature. Imagine you start a water simulation from an ice configuration, and chose the kinetic energy from a Boltzmann at 300K. The temperature would become constant, but not 300K. You can choose whatever velocities you want to get the correct total energy. That's what I do.

How am I ever supposed to get a meaningful temperature profile with Muller
Plathe if my overall temperature and that of each individual slice keeps
fluctuating so much (290 to 310K at least).

Why do you think the fluctuations are unreasonable? I think you have 64 atoms per bin? I would expect large fluctuations, so you'll need to average more.

Note that I output the thermo data and the temperature of the slices every
1000 steps (1fs steps).
I have already tried waiting for a very long time (10M steps = much longer
than people usually do).

I don't know how quickly it should converge for your system. Maybe you need to swap more frequently. No harm in averaging every step, too.