[lammps-users] About thermal conductivity simulation

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users:

I desperately need some help. I intend to calculate thermal conductivity of uranium dioxide by “fix thermal/conductivity” commend. I went through the history of this Mailing-List and found some useful posters and used some of the strategies in my input file, which is pasted in the end of this mail. Please be advised the potential is customized and has been tested through thermal expansion simulation.

My questions are listed below:

1, the uranium dioxide is partially ionic system that has long range columbic interaction. And I only let “fix thermal/conductivity” to exchange the kinetic energies between oxygen particles. I wonder whether this method is valid or I should exchange the kinetic energies of uranium particles as well.

2,the parameters for “fix thermal/conductivity” and corresponding time averages are shown below.

fix heat_swap_o oxygen thermal/conductivity 10 z 10

fix e_exchange all ave/time 10 1000 10000 f_heat_swap_o file e_exchange.dat ave one

The result I got in the e_exchange.data is
Timestep Energy

20000 71.3858

30000 213.342

40000 355.356

50000 498.597

60000 644.172

70000 787.423

Interesting thing is that the outputted value increase by approximately the same interval which is two folds of the very first one. I though the output should be equilibrate at a constant value. I don’t know how to interpret this result. Does anyone can help me with this?

3, About the temperature profile, I used a linear function to fit the result. And the R is about 0.98, I wonder whether this linearity is good enough to infer the temperature gradiant.

4, The last question is about the equation to calculate the thermal conductivity. According to the manual, the thermal conductivity should be

Thermal conductivity=(energy exchange from “fix thermal/conductivity”)/{(cross section of the simulation box)(Time which corresponds to the Nevery in the “fix thermal/conductivity ”)(Slope of temperature gradient) }

I wonder whether I got this right?

Any comment or suggestion would be highly appreciated.
Thank you.

Below: Input file

# Thermal conductivity of solid UO2 @ 300K

units metal

boundary p p p

atom_style charge

variable T equal 300

lattice fcc 5.4704

region box block -3 3 -3 3 -3 3

create_box 2 box

lattice fcc 5.4704

create_atoms 1 box

lattice custom 5.4704 a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0.25 0.25 0.25 basis 0.75 0.25 0.25 basis 0.25 0.75 0.25 basis 0.75 0.75 0.25 basis 0.25 0.25 0.75 basis 0.75 0.25 0.75 basis 0.25 0.75 0.75 basis 0.75 0.75 0.75

create_atoms 2 box

timestep 0.001

mass 1 238

mass 2 16

group uranium type 1

group oxygen type 2

set group uranium charge 2.2208

set group oxygen charge -1.1104

kspace_style pppm 1e-4

pair_style hybrid/overlay table linear 88001 coul/long 10

pair_coeff 1 1 table Yakub.table Uranium 10

pair_coeff 2 2 table Yakub.table Oxygen 10

pair_coeff 1 2 table Yakub.table Dioxide 10

pair_coeff * * coul/long

neighbor 0.4 bin

neigh_modify every 5 delay 0 check yes

velocity all create 300 825577 dist gaussian

fix temp all temp/berendsen 300 300 0.1

fix nve all nve

compute ke all ke/atom

variable temp atom c_ke/(1.5*8.617343e-5)

fix temp_profile all ave/spatial 1 10000 10000 z lower 0.6 v_temp file temp.profile units lattice

thermo_style custom step temp etotal vol enthalpy

thermo_modify lost warn

thermo 1000

run 10001

# Impose heat flux

unfix temp

fix heat_swap_o oxygen thermal/conductivity 10 z 10

fix e_exchange all ave/time 10 1000 10000 f_heat_swap_o file e_exchange.dat

# Output heat flux value

run 150001

2

  1. I think that exchanging KE between just the O’s or between both the U’s and O’s are both OK. Probably better to do both. If in doubt, try both approaches and see if you get the same answer. I don’t think that the following applies to your system, but here’s a note from the docs just in case: “LAMMPS does not check, but you should not use this fix to swap the kinetic energy of atoms that are in constrained molecules, e.g. via fix shake or fix rigid. This is because application of the constraints will alter the amount of transferred momentum. You should, however, be able to use flexible molecules. See the Zhang paper for a discussion and results of this idea. When running a simulation with large, massive particles or molecules in a background solvent, you may want to only exchange kinetic energy bewteen solvent particles.”

  2. The scalar value that fix thermal/conductivity produces is cumulative, so as you’ve discovered, it continually grows. See: http://lammps.sandia.gov/doc/fix_thermal_conductivity.html Specifically: "This fix computes a global scalar which can be accessed by various output commands. The scalar is the cummulative kinetic energy transferred between the bottom and middle of the simulation box (in the edim direction) is stored as a scalar quantity by this fix. This quantity is zeroed when the fix is defined and accumlates thereafter, once every N steps. The units of the quantity are energy; see the units command for details. The scalar value calculated by this fix is “intensive”. "

  3. Sounds pretty good to me.

  4. Also looks fine. Make sure the units work out.

Paul