self-diffusion and thermal conductivity of h by molecular

SELF-DIFFUSION AND THERMAL
CONDUCTIVITY OF H2 BY MOLECULAR
DYNAMICS SIMULATIONS
EIRIK KJØNSTAD
[email protected]
SARAI DERY FOLKESTAD
[email protected]
TKJ4200 IRREVERSIBLE THERMODYNAMICS PROJECT
REPORT
Abstract
Equilibrium molecular dynamics (EMD) and non-equilibrium molecular dynamics (NEMD) simulations were performed to calculate the
thermal conductivity of hydrogen at temperatures of 100 K, 150 K,
and 200 K and at densities ranging from 0.0135 g/cm3 to 0.0314 g/cm3 .
The temperature-density curves from the NEMD simulations were compared with experimental equilibrium data. The self-diffusion coefficient
was also determined by EMD. Total simulation time was 2 ns. The
canonical ensemble (NVT) was used in the EMD simulations, and the
microcanonical ensemble (NVE) for NEMD simulations.
The thermal conductivity was found to increase with increasing temperature, and increase with increasing density. Comparison of NEMD
temperature-density curves with experimental equilibrium data show,
if the assumption of local equilibrium holds, a systematic error in the
force field. A comparison of the NEMD results with EMD simulations
of the same force field, is needed to confirm local equilibrium.
The self-diffusion coefficient was found to increase with increasing temperature, and decrease with increasing density. The self-diffusion coefficient calculated by EMD, from the mean square displacement, converges quickly. The EMD method is considered to be efficient, also for
short simulation times, and the errors will be due to the force field.
The results, in thermal conductivity and self-diffusion, showed a dependence on temperature and density that was expected from kinetic
gas theory.
The EMD, Green-Kubo, thermal conductivity calculations introduce
a non-systematic error not present in the NEMD simulations. The
heat-flux autocorrelation function converges slowly, and thus the error
is significant for short run times. NEMD is therefore thought to be
a superior method for determining the thermal conductivity, and at
temperatures of 150 K and 200 K, this was confirmed by the simulation
results.
Contents
1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
2. Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2.1. The Thermodynamics of the System . . . . . . . . . . . . . . . . . . .
2
2.2. Compressibility Factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.3. Atomic and Molecular Interactions . . . . . . . . . . . . . . . . . . . .
7
2.4. Diffusion Coefficient and Thermal Conductivity . . . . . . . .
8
3. Procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3.1. System Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3.2. Simulation Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3.3. Comparison of the Force Field to Experimental Data . . 11
3.4. Compressibility Factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4. Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.1. Diffusion Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
4.2. Thermal Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.3. Comparison of the Force Field to Experimental Data . . 18
4.4. Compressibility Factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5. Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.1. Diffusion Coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.2. Thermal Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
5.3. Comparison of the Force Field to Experimental Data . . 22
5.4. Compressibility Factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
6. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
7. Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
List of Symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
Appendix A. Compressibility factor . . . . . . . . . . . . . . . . . . . . . . . . . . 31
Appendix B. MATLAB scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
Appendix C. LAMMPS scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
1
1. Introduction
Carbon materials have many applications in hydrogen powered fuel
cells[1], for example as catalytic support[2]. One of the main challenges of fuel cells, particularly for vehicular use, is on-board storage of
hydrogen [3]. None of the methods of hydrogen storage existing today,
compression, liquefaction at cryogenic conditions or storage as metal
hydrides, meet the requirements for efficiency, capacity and safety[4].
The low molecular weight and production costs of some carbon allotropes make adsorption of hydrogen on carbon materials a viable
option for on-board storage[4], although such systems do not, as yet,
meet the requirements of storage capacity[3]. On these counts, systems of hydrogen gas on solid carbon materials, such as graphite, are
of interest.
Molecular dynamics (MD) simulations can be used to successfully
model such systems, as done by Haas et al.[2]. Haas et al. modeled
the hydrogen-carbon using Equilibrium Molecular Dynamics (EMD)
simulations. Non-Equilibrium Molecular Dynamics (NEMD) is analogous to the experimental situation, and may yield better results for
properties that are known to be diffucult to obtain by EMD, such as
the thermal conductivity[5].
As a first step in modelling the carbon-hydrogen system, we model
hydrogen alone to determine transport properties. The thermal conductivity and self diffusion of hydrogen is determined for 100 K, 150
K and 200 K at five different densities. The density of the system is
varied by keeping the box volume constant at 400 ˚
A × 100 ˚
A × 100 ˚
A
and increasing the number of molecules by 5000 from 15000 to 35000.
The thermal conductivity will be found both by EMD simulations using
the Green-Kubo relations, and with NEMD simulations using Fourier’s
law. The values obtained from the two methods will be compared.
Temperature-density curves obtained from the NEMD simulations will
also be compared with experimental equilibrium data.
Figure 1. Snapshot of the simulation.
2
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
2. Theory
While the fluxes of the hydrogen system can be described by simple transport laws, non-equilibrium thermodynamics is needed to describe the coupled transport of heat and mass in the extended carbonhydrogen system. Since the objective of this project is to model a
carbon-hydrogen system, the required theory is presented below.
2.1. The Thermodynamics of the System. Linear non-equilibrium
thermodynamics can be used to describe the coupled transport and the
entropy production of a system.
The entropy production of a system with multiple transport phenomena, is the product sum of conjugate fluxes Ji and driving forces Xi ,
σ=
n
X
Ji Xi .
(2.1)
i=1
For a system of two chemical components and a temperature gradient, where one component is used as a frame of reference, the entropy
production is
σ=
Jq′
1 ∂T
1 ∂µ1,T
− 2
+ J1 −
,
T ∂x
T ∂x
(2.2)
where Jq′ is the measurable heat flux, J1 is the mass flux of component
1, T is the temperature and µ1,T is the chemical potential of component
1 at constant temperature T .
2.1.1. Derivation of Entropy Production. Assuming local equilibrium,
Gibbs’ law is valid locally. The entropy production can then be derived
by combining the local form of Gibbs’ law with mass balances and the
first law of thermodynamics. For a one-dimensional two-component
system, in which only mass and heat is transported, the derivation is
as follows.
For an infinitesimal volume element between x and x + dx, the rate of
change in entropy density, s, is given by
∂s
∂
= − Js + σ,
∂t
∂x
(2.3)
where Js is the entropy flux through the volume element and σ is the
entropy production of the element.
3
The volume element must contain enough particles to give a statistical
basis for thermodynamic calculations[6]. The state of the volume element is defined by the temperature, T (x), pressure, p(x), and chemical
potential of component i, µi (x).
The general Gibbs’ equation is given by
dU = T dS − pdV +
n
X
µi dNi ,
(2.4)
i=1
where U is the internal energy of the system, T the temperature, S the
entropy, p the pressure, V the volume, and Ni is the amount of moles
of component i. Integrating (2.4) yields
U = T S − pV +
n
X
µi Ni .
(2.5)
i=1
Dividing equation (2.5) by the volume V and performing the variable
changes u = U/V , s = S/V , yields
u = Ts − p +
n
X
µ i ci ,
(2.6)
i=1
where ci = Ni /V is the concentration of component i. On differential
form, this is the local Gibbs’ equation
du = T ds +
n
X
µi dci .
(2.7)
i=1
By differentiating with respect to time, t, the rate of change of the
entropy density can be expressed as
n
∂s
1 ∂u
1 X ∂ci
µi
=
−
.
∂t
T ∂t
T i=1 ∂t
(2.8)
The mass balance for component i, when no chemical reaction takes
place, is
∂ci
∂
= − Ji .
∂t
∂x
(2.9)
4
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
The first law of thermodynamics for a system, when only heat is exchanged with the system, yields
∂u
∂
= − Jq ,
∂t
∂x
(2.10)
where Jq is the heat flux. Inserting (2.9) and (2.10) into (2.7) yields
∂s
1 ∂
µ1 ∂
=−
Jq +
J1 .
∂t
T ∂x
T ∂x
(2.11)
By using the product rule of derivation, (2.11) becomes
∂
∂s
=−
∂t
∂x
Jq J1 µ 1
−
T
T
+ Jq
∂ 1
∂ µ1
− J1
.
∂x T
∂x T
(2.12)
The entropy flux is defined by
∂
Js = −
∂x
Jq J1 µ 1
−
T
T
.
(2.13)
By setting (2.2) equal to (2.12), the entropy production,
σ = Jq
∂ µ1
∂ 1
− J1
,
∂x T
∂x T
(2.14)
is obtained. It is practical to express the entropy production in terms
of the measurable heat flux, Jq′ . The relation between the heat flux,
Jq , and Jq′ is given by
Jq = Jq′ +
n
X
H i Ji ,
(2.15)
i=1
where Hi is the molar enthalpy of component i.
Inserting (2.15) into (2.14), and using the product rule, yields
σ=
Jq′
∂ 1
∂ 1
+ H 1 J1
− J1
∂x T
∂x T
∂ 1
1 ∂
µ1 + µ1
T ∂x
∂x T
.
(2.16)
By using the relation
∂
∂
∂T
µ1 =
µ1,T + S1
,
∂x
∂x
∂x
(2.17)
5
where S1 is the molar entropy of component 1, it is possible to rearrange
(2.16) into
σ = Jq′
∂ 1
1 ∂
∂ 1
− J1
µ1,T + J1 (H1 − µ1 − T S1 )
.
∂x T
T ∂x
∂x T
(2.18)
The relation between µi , Hi and Si is
H i = µ i + T Si ,
(2.19)
thus the entropy production for the system is therefore given by
σ = Jq′
1 ∂
∂ 1
− J1
µ1,T .
∂x T
T ∂x
(2.20)
This is equivalent to equation (2.2).
2.1.2. Flux Equations. In linear non-equilibrium thermodynamics, the
fluxes is expressed as the weighted sum of all driving forces within the
system,
Ji =
X
Lij Xj ,
(2.21)
j
the coefficients Lij being independent of the driving forces Xj . Onsager
showed that, for a microscopically reversible system[6],
Lij = Lji ,
(2.22)
which are known as the Onsager reciprocal relations.
2.1.3. Carbon-Hydrogen System. For the carbon-hydrogen system, a
graphite layer is placed under a three-dimensional box containing hydrogen and there is a temperature gradient is along the x direction of
the graphite surface, the coupled transport of heat and mass is given
by
1 ∂T
1 ∂µH2 ,T
+ Lµq − 2
JH2 = Lµµ −
T ∂x
T ∂x
1 ∂µH2 ,T
1 ∂T
′
Jq = Lqµ −
+ Lqq − 2
T ∂x
T ∂x
when the graphite surface is used as a frame of reference.
(2.23)
(2.24)
6
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
The simple transport laws for heat and mass may be derived as a special
case of the flux equations. If there is no temperature gradient, then
∂T /∂x = 0, and (2.23) is reduced to
JH 2
1 ∂µH2 ,T
= Lµµ −
.
T ∂x
(2.25)
Assuming that the activity is 1, one finds Fick’s law,
JH 2 = −
Lµµ ∂c
,
T ∂x
(2.26)
where the diffusion coefficient is defined as D = Lµµ /T .
Similarly, if there is no concentration gradient, equation (2.24) reduces
to Fourier’s law,
Jq′ = −
Lqq ∂T
,
T 2 ∂x
(2.27)
where the thermal conductivity is defined as λ = Lqq /T 2 .
For the system with only hydrogen gas, there is no diffusion flux. The
heat flux is described by Fourier’s law. It can be written in terms of
averages as[5]
λ=
lim
lim −
∂T /∂x→0 t→∞
hJq (t)i
,
h∂T /∂xi
(2.28)
where hJq (t)i and h∂T /∂xi are the ensemble average of the heat flux
and temperature gradient, respectively. For this system the entropy
production is
σ = Jq′
∂ 1
.
∂x T
(2.29)
2.2. Compressibility Factor. An ideal gas is defined as a gas for
which the equation of state is
pV = nRT,
(2.30)
where R is the gas constant. The deviation of a real gas from the ideal
can then be measured by defining the compressibility factor, Z, by
Z=
Preal
.
Pideal
(2.31)
7
An interpretation[7] of the compressibility factor is that Z < 1 indicates
that attractive interactions are dominating in the gas, while Z > 1
indicates that repulsive interactions are dominating.
2.3. Atomic and Molecular Interactions. In molecular mechanics,
when modelling intramolecular bonds, such as the H−H bond in the
H2 molecule, the Morse potential VM has the appropriate shape. It is
given by
VM (r) = De 1 − e−α(r−r0 )
2
,
(2.32)
where r is the distance between two H-atoms in the hydrogen molecule,
r0 is the “equilibrium” bond length, α determines the rigidity of the
bond, and De determines the depth of the potential.
Interatomic interactions may be modelled by a Lennard-Jones (LJ)
potential. This potential contains a repulsive term, varying as r −12 ,
accounting for exchange effects, and an attractive term, varying as
r −6 , accounting for dispersive interactions[8]. The 12-6 LJ potential
VLJ , with one atom type, is given by
VLJ (rij ) = 4ǫLJ
σ
rij
12
−
σ
rij
6 !
,
(2.33)
where rij is the distance from atom i to atom j, ǫLJ determines the
potential well depth, and σ is the separation diameter, the distance for
which VLJ = 0. A truncated LJ potential may be used to lessen the
computional expense. Truncation means neglecting the interaction at
large distances. To find an explicit form of the truncated LJ potential,
one can define the cut-off distance rC as the distance at which, for
larger distances, the interaction between the molecules is zero. The
truncated LJ potential Vnb can therefore be defined as
(
VLJ (r) − VLJ (rC ) if r ≤ rC ,
Vnb (r) =
0
if r > rC .
(2.34)
8
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
9
8
7
VM [eV]
6
5
4
3
2
1
0
−1
0.2
0.4
0.6
0.8
1
˚]1.2
x[A
1.4
1.6
1.8
2
Figure 2. The Morse potential, VM , plotted as a function of interatomic distance, x.
160
140
120
VLJ [K]
100
80
60
40
20
0
−20
−40
2.5
3
3.5
x[˚
A]
4
4.5
5
Figure 3. The Lennard-Jones potential, VLJ , plotted as
a function of intermolecular distance, x.
2.4. Diffusion Coefficient and Thermal Conductivity.
2.4.1. Self-diffusion Coefficient. The self diffusion coefficient is defined
by Fick’s first law, see equation (2.26), and was shown by Einstein[8]
for a d-dimensional system to be given by
D=
h|r(t) − r(0)|2i
1
lim
,
d t→∞
2t
(2.35)
where h|r(t) − r(0)|2 i is the mean square displacement at time t. If
the fraction, h|r(t) − r(0)|2 i /2t, is constant for some t > t0 , the selfdiffusion coefficient is proportional to the mean square displacement,
h|r(t) − r(0)|2 i = (2dD)t t > t0 .
(2.36)
9
2.4.2. Thermal Conductivity. The Green-Kubo (GK) relations can be
used to determine the transport properties for a system in equilibrium.
This is based on local fluctuations occuring in the gas[8]. A description
of the GK relation for determining the thermal conductivity λ is given
below.
The thermal conductivity can be found from the auto-correlation function, hJ(τ ) · J (0)i, of the heat flux, J (t). The auto-correlation function
is the ensemble average of the dot product J (τ ) · J (0). The explicit
relation is[9]
1
λ=
3V kB T 2
Z
∞
hJ (τ ) · J (0)i dτ,
(2.37)
0
where kB is the Boltzmann constant.
3. Procedure
The LAMMPS Molecular Dynamics Simulator[10] was used for running
the simulations.
3.1. System Geometry. The simulations were performed in a cuboid
box of dimensions 400 ˚
A × 100 ˚
A × 100 ˚
A, with periodic boundary conditions in the x, y, and z direction. The x direction is the direction for
which the length is 400 ˚
A. The simulations were done at temperatures
of 100 K, 150 K, and 200 K, and at different densities, the number of
H2 molecules ranging from NH2 = 15000 to NH2 = 35000 at intervals
of 5000.
3.2. Simulation Details. The time evolution of the system was found
by solving Newton’s equations of motion with a standard velocityVerlet algorithm[11]. The timestep was set to 1 fs, and the number
of timesteps used per simulation was 106 , corresponding to a total run
time of 1 ns. For each temperature and density, two consecutive runs
were performed. One equilibration run, and one run where the desired
properties were calculated, giving a total simulation time of 2 ns.
For the equilibration run, the initial configuration was set by choosing
velocities from a gaussian distribution with the given mean temperature, 100 K, 150 K, or 200 K. The canonical ensemble (NVT) was used
for equilibration.
3.2.1. EMD Simulations. For the EMD simulations, NVT was also
used for the calculation run. The mean square displacement and thermal conductivity was averaged and stored every 2000 steps.
10
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
Figure 4. The simulation box.
For the calculation of the self-diffusion coefficient, the mean square
displacement was plotted against the time, see equation (2.36). A
linear regression line was adapted to the data. The regression line was
used to confirm the linear relation and used to calculate the diffusion
coefficient.
3.2.2. NEMD Simulations. For the NEMD simulations NVE was used
for the calculation run.
The box was divided into 40 slabs in the x direction. In each segment,
the density and temperature were averaged spatially and stored every
100 steps. The cumulative kinetic energy transferred between the first
and the middle slab was stored every 2000 steps.
To impose a temperature gradient along the x axis, the Muller-Plathe
method[5] was used. First, the molecule of highest velocity in the first
slab, and the molecule of lowest velocity in the middle slab, was identified. Then their velocities were swapped. This generated a transfer
of kinetic energy, resulting in a temperature gradient.
Fourier’s law in terms of averages, see equation (2.28), was used to determine the thermal conductivity. Specifically, if h∆Ek i is the average
kinetic energy transfer per swap, δt is the time between swaps, and Ω
11
is the cross-sectional area, then
λ=
h∆Ek i
.
i
2δtΩ h ∂T
∂x
(3.1)
The factor 2 in the denominator is added because of the periodic boundary conditions.
A linear regression line was fitted to the temperature profile in both
halves of the box. The average of the temperature gradients from the
two linear regression lines was used to calculate the thermal conductivity by Fourier’s law. The temperatures of the first and middle slabs,
where heat was added or taken out, was not included in the regression.
3.2.3. Simulation Parameters. To model the intermolecular interactions between hydrogen molecules, a truncated LJ-potential was used,
see equation (2.34), for which the cut-off distance was set to rC =
12.0 ˚
A.
A Morse potential was chosen to model the H−H bond, see equation
(2.32). The parameters used are listed in Table 3.1. The parameters in
the Morse potential is from the MM2 force field[12], and the parameters
of the Lennard-Jones potential is from the Delft Molecular Mechanics
(DMM) force field [13].
Table 3.1. Simulation parameters for intra- and intermolecular interactions[12][13].
Morse potential Lennard-Jones potential
De = 4.747 eV
α = 1.946 ˚
A−1
r0 = 0.7414 ˚
A
ǫLJ = 27.655 K
σ = 2.63984 ˚
A
3.3. Comparison of the Force Field to Experimental Data. To
check how well the potentials of the force field describe the real system,
the average density of each of the forty slabs was plotted against the
average temperature of that same slab. This was then compared to
experimental equilibrium data from the NIST database[14].
3.4. Compressibility Factor. The pressures from the EMD simulations, and for the equlibration of the NEMD simulations were compared
to the pressure of the ideal gas of the given temperature and density.
12
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
4. Results
4.1. Diffusion Coefficient. The diffusion coefficient of hydrogen gas
was determined and the values obtained are listed in Table 4.1 and
presented in Figure 6. The mean square displacement for 100 K and
15000 molecules is shown in Figure 5.
5
5
x 10
4.5
h|r(t) − r(0)|2i [˚
A2 ]
4
3.5
3
2.5
2
1.5
1
0.5
0
0
0.2
0.4
0.6
0.8
1
t[fs]
1.2
1.4
1.6
1.8
2
6
x 10
Figure 5. Mean square displacement, h|r(t) − r(0)|2 i,
plotted as a function of time, t, for 15000 molecules and
100 K.
13
−3
8
x 10
100K
150K
200K
7
D[cm2 /s]
6
5
4
3
2
1
6
7
8
9
10
11
ρ[g/cm3 ]
12
13
14
15
Figure 6. The self-diffusion coefficient, D, of hydrogen
gas from EMD simulations at density, ρ, for the T = 100
K, T = 150 K, and T = 200 K isoterms.
Table 4.1. The self-diffusion coefficient, D, of hydrogen
gas by EMD simulations at temperature, T , density, ρ.
T [K] ρ[g/cm3 ]
D[cm2 /s]
100.0
100.0
100.0
100.0
100.0
0.0135
0.0179
0.0224
0.0269
0.0314
4.06 · 10−3
3.07 · 10−3
2.46 · 10−3
2.04 · 10−3
1.78 · 10−3
150.0
150.0
150.0
150.0
150.0
0.0135
0.0179
0.0224
0.0269
0.0314
5.69 · 10−3
4.21 · 10−3
3.37 · 10−3
2.75 · 10−3
2.37 · 10−3
201.0
200.0
200.0
200.0
200.0
0.0135
0.0179
0.0224
0.0269
0.0314
7.01 · 10−3
5.24 · 10−3
4.12 · 10−3
3.39 · 10−3
2.86 · 10−3
14
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
4.2. Thermal Conductivity. The thermal conductivity of hydrogen
gas was determined by EMD and NEMD simulations and are listed in
Table 4.2 and Table 4.3, respectively. The values obtained by both
methods are presented as isotherms, and compared with literature
values[15], see Figure 7–9.
0.24
Literature
EMD simulation
NEMD simulation
0.22
0.2
λ[W/mK]
0.18
0.16
0.14
0.12
0.1
0.08
0.06
0
0.01
0.02
0.03
0.04
ρ[g/cm3 ]
0.05
0.06
0.07
0.08
Figure 7. Thermal conductivity of hydrogen gas, λ,
from literature[15], EMD-, and NEMD simulations at
100 K for densities ρ.
0.24
Literature
EMD simulation
NEMD simulation
0.22
λ[W/mK]
0.2
0.18
0.16
0.14
0.12
0.1
0
0.01
0.02
0.03
0.04
ρ[g/cm3 ]
0.05
0.06
0.07
Figure 8. Thermal conductivity of hydrogen gas, λ,
from literature[15], EMD-, and NEMD simulations at
150 K for densities ρ.
15
0.3
Literature
EMD simulation
NEMD simulation
0.28
0.26
λ[W/mK]
0.24
0.22
0.2
0.18
0.16
0.14
0.12
0
0.01
0.02
0.03
ρ[g/cm3 ]
0.04
0.05
0.06
Figure 9. Thermal conductivity of hydrogen gas, λ,
from literature[15], EMD-, and NEMD simulations at
200 K for densities ρ.
NEMD
Regression curve left, R2 = 0.9987
Regression curve right, R2 = 0.9990
115
110
T [K]
105
100
95
90
85
0
50
100
150
200
x[˚
A]
250
300
350
400
Figure 10. Temperature profile of the simulation box
for 15000 molecules and 100 K. The temperature, T , is
plotted as a function of box position, x. Fitted regression
curves are shown for the left and right side of the box.
16
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
0.017
0.016
ρ[g/cm3 ]
0.015
0.014
0.013
0.012
0.011
0
50
100
150
200
x[˚
A]
250
300
350
400
Figure 11. Density profile of the simulation box for
15000 molecules and 100 K. The density, ρ, is plotted as
a function of box position, x.
Table 4.2. Thermal conductivity, λ, by EMD simulations. The temperature and density is given by T and ρ,
respectively.
T [K] ρ[g/cm3 ] λ[W/mK]
100.0
100.0
100.0
100.0
100.0
0.0135
0.0179
0.0224
0.0269
0.0314
0.113
0.135
0.152
0.164
0.193
150.0
150.0
150.0
150.0
150.0
0.0135
0.0179
0.0224
0.0269
0.0314
0.175
0.168
0.198
0.206
0.207
201.0
200.0
200.0
200.0
200.0
0.0135
0.0179
0.0224
0.0269
0.0314
0.213
0.223
0.238
0.251
0.260
17
Table 4.3. Thermal conductivity, λ, by NEMD simulations. The average temperature and density is given by
Tav and ρav , respectively. R¯2 is the average coefficient of
determination from the two linear regression lines. The
average temperature difference, obtained by linear regression, is ∆av T .
Tav [K] ρav [g/cm3 ] ∆av T [K] λ[W/mK]
R¯2
98.80
99.99
100.00
100.20
100.20
0.0135
0.0179
0.0224
0.0269
0.0314
18.9
14.6
6.3
6.7
10.2
0.131
0.148
0.184
0.178
0.238
0.999
0.997
0.992
0.990
0.986
147.90
148.50
148.40
148.60
149.40
0.0135
0.0179
0.0224
0.0269
0.0314
23.8
23.0
21.4
20.1
18.2
0.156
0.171
0.190
0.209
0.237
0.997
0.999
0.999
0.998
0.998
196.40
197.00
197.60
198.00
198.70
0.0135
0.0179
0.0224
0.0269
0.0314
19.3
25.0
25.3
23.9
21.5
0.200
0.209
0.213
0.230
0.266
0.996
0.998
0.999
0.999
0.998
18
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
4.3. Comparison of the Force Field to Experimental Data. The
performance of the potentials compared to the real system, were investigated by comparing the NEMD temperature-density curves with
experimental equilibrium data from NIST[14]. The results are listed in
Table 4.4 and shown in Figure 12 and 13.
0.035
15k NEMD
15k NIST
20k NEMD
20k NIST
25k NEMD
25k NIST
30k NEMD
30k NIST
35k NEMD
35k NIST
ρ[g/cm3 ]
0.03
0.025
0.02
0.015
0.01
120
130
140
150
T [K]
160
170
180
190
Figure 12. Isobars in NEMD simulations for Tav = 150
K, density ρ, and temperature T , compared to experimental data from the NIST database[14].
0.035
15k NEMD
15k NIST
20k NEMD
20k NIST
25k NEMD
25k NIST
30k NEMD
30k NIST
35k NEMD
35k NIST
ρ[g/cm3 ]
0.03
0.025
0.02
0.015
0.01
170
180
190
200
T [K]
210
220
230
240
Figure 13. Isobars in NEMD simulations for Tav = 200
K, density ρ, and temperature T , compared to experimental data from the NIST database[14].
19
Table 4.4. Average temperature, Tav , density ρ, pressure, p, and average relative error between NEMD calculations and experimental equilibrium data from the NIST
database[14], ǫ = (NEMD − NIST )/NIST .
Tav [K] ρav [g/cm3 ] p[bar]
ǫ
100.0
100.0
100.0
100.0
100.0
0.0135
0.0179
0.0224
0.0269
0.0314
50.6
63.1
73.9
83.4
91.8
0.0831
0.184
0.28231
0.3783
0.4704
150.0
150.0
150.0
150.0
150.0
0.0135
0.0179
0.0224
0.0269
0.0314
87.7 -0.00540
115.4 0.0389
142.6 0.0787
170.0 0.1184
198.1 0.1594
201.0
200.0
200.0
200.0
200.0
0.0135
0.0179
0.0224
0.0269
0.0314
123.6
164.1
206.6
249.8
294.4
-0.0308
0.0035
0.0322
0.0642
0.0951
0.8
Error 100K
0.7
Error 150K
Error 200K
0.6
Regression 100K, R2=0.9997
Regression 150K, R2=0.9996
ρ[g/cm3 ]
0.5
Regression 200K, R2=0.9994
0.4
0.3
0.2
0.1
0
−0.1
0.01
0.015
0.02
ǫ[−]
0.025
0.03
0.035
Figure 14. Average relative error, ǫ, plotted as a function of density for average temperatures of 100 K, 150 K,
and 200 K with linear regression lines.
20
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
4.4. Compressibility Factor. The compressibility factor for the
EMD simulations are listed in Table A.1 and shown in Figure 15.
1.25
100K
150K
200K
1.2
1.15
1.1
Z
1.05
1
0.95
0.9
0.85
0.8
0.75
50
100
150
p[bar]
200
250
300
Figure 15. Compressibility factor, Z, at pressure, p,
and temperatures 100 K, 150 K, and 200 K.
21
5. Discussion
5.1. Diffusion Coefficient. The mean square displacement was plotted as a function of the simulation time, see Figure 5, and a linear
regression line was adapted to the data. The coefficient of determination, R2 , confirmed a clear linear relation between the time and mean
square displacement throughout the entire simulation. This indicates
that the fraction in equation (2.35) converges quickly, and therefore
that the method for calculating the self-diffusion coefficient is accurate. Possible errors will most likely be systematic and due to the force
field.
There are two definite trends in the values of the diffusion coefficients
obtained from the simulation, see Figure 6. Firstly, the diffusion coefficient, for a constant density, increases with increasing temperature.
Increasing temperature is equivalent to increasing kinetic energy. As
the mean square displacement increases with the kinetic energy a rise
in diffusion coefficient is expected. Secondly, for all three isotherms,
the coefficient values decrease when the density of the system increases.
For a real gas, the mean free path, and therefore also the mean square
displacement, is reduced when the density of the box increases. This
is due to the increasing crowdedness of the box.
5.2. Thermal Conductivity. The thermal conductivity obtained
from both EMD and NEMD simulations correspond well with, but
overestimate, the experimental literature values, see Figure 7–9. The
simulation results, EMD and NEMD, indicate that there is a systematic
error.
The thermal conductivity at 100 K, for the five different densities, see
Figure 7, show that EMD and NEMD both methods determine the
thermal conductivity with the same accuracy. For 100 K and 35000
molecules, the value from the NEMD simulation is an outlier. Modelling hydrogen gas at temperatures below 100 K may be unreliable,
due to quantum effects, which are not taken into account in the interaction potentials we have chosen.
For 150 K the NEMD values are slightly better, compared to the experimental values, see Figure 8. For the 200 K simulations, see Figure
9, NEMD yields significantly better values than EMD.
For the EMD approach, using the Green-Kubo relations, the truncation
of the integral (2.37) from t = 0 to t = ∞ may introduce some error.
In addition, the heat-flux autocorrelation function converges slowly[9].
This error is non-systematic, and may be remedied by increasing total
simulation time.
22
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
The Muller-Plathe NEMD method depends on the convergence of
h∂T /∂xi, which is expected to be fast[5]. Therefore, from Fourier’s
law, see equation (3.1), the method is expected to have good accuracy.
Since the temperature and densities vary over the system, the calculated thermal conductivity will be an average over these temperatures
and densities. This can introduce errors in the determined values if the
temperature difference is large[5].
At steady state, the heat flux, hJ(t)i, and temperature gradient,
h∂T /∂xi, will be constant by Fourier’s law. This means that for a
well-behaved system, the temperature profile is linear. The temperature gradient was determined as the average of the absolute value of
the slopes from the linear regression, see Figure 10. The average coefficients of determination, R¯2 , from the two regression lines, were all
higher than 0.98. Hence, the temperature profiles are linear.
5.3. Comparison of the Force Field to Experimental Data. In
Figure 12-13, the densities and temperatures of each of the 40 slabs,
at the given average temperature of 150 K and 200 K, are plotted as
isobars and compared to experimental data. For each comparison, the
average relative error compared to the experimental data is given in
Table 4.4. For the three lowest densities, at 150 K and 200 K, the
average relative error is below 8%. For 100 K, and at high densities for
150 K, the average relative error compared with experiment is much
greater. The high deviation for 100 K may be due to quantum effects.
The average relative error, ǫ, is plotted against the density in Figure
14. Since R2 > 0.999 for all three isotherms, the average relative error
varies linearly with respect to density. The simulation data follows the
trend of the experimental data, see Figure 12-13, and thus the average
relative error, ǫ, is approximately equal to the actual relative error at a
given temperature. Based on the assumption of local equilibrium, this
error is systematic. It would be of interest to check if EMD simulations,
with the force field in question, would show this same relative error
trend. If so, the NEMD data would be in good agreement with EMD,
validating the assumption of local equilibrium, and hence, of systematic
error. Further simulations are needed to confirm this.
5.4. Compressibility Factor. The compressibility factor, Z, was
used as a first check that the simulation system was well-behaved. For
the EMD simulations, Z is shown in Figure 15. The compressibility
factor increases with increasing temperature. For the temperatures
and pressures in these simulations, hydrogen gas is expected not to
deviate largely from the ideal gas. For 100 K, Z < 1, and attractive
interactions dominate. For 150 K and 200 K, Z > 1, and repelling
23
interactions dominate. The calculated Z is consistent with the LJ potential, see Figure 3.
When the temperature increases, and thus has higher kinetic energy,
the molecules are, classically, allowed to travel closer to each other.
At a certain intermolecular distance the potential becomes positive,
VLJ > 0. For higher temperatures, more particles will be in this positive
region of the LJ potential. Thus, by increasing the temperature, the
intermolecular interactions will shift from being dominantly atractive
to dominantly repelling.
For increasing density, Z will increase in some cases and decrease in
other cases. This can also be explained by the LJ potential. For a
given density there is a mean distance between molecules. Consider a
hydrogen molecule and its interactions with other hydrogen molecules
at two densities, ρ1 > ρ2 . It is then more likely that the neighbours
are on average closer to each other for the gas at density ρ1 than for
the gas at density ρ2 . If an increase in density leads to more molecules
occupying a region of the LJ potential where Vnb is lower, Z will be
reduced. Equivalently, more molecules occupying a region of the LJ
potential where Vnb is higher, will cause Z to increase.
24
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
6. Conclusion
The Green-Kubo method and the NEMD method will in theory yield
exact values for the thermal conductivity. Since the the temperature gradient converges more quickly than the heat-flux autocorrelation function, the NEMD method is superior. This is confirmed by our
results at 150 K and 200 K.
The EMD simulation results indicate that the mean square displacement converges quickly, and that EMD simulations are sufficient for
calculating the self-diffusion coefficient. Such values were calculated
and compared with experimental data by Haas et al.[2].
Comparison of the NEMD results with experimental data indicate a
systematic error, if the assumption of local equilibrium holds. To confirm local equilibrium, a comparison of the NEMD results with EMD
simulations for the given force field should be performed.
25
7. Acknowledgements
Thanks to Dr. Thuat T. Trinh, who supervised and helped us in our
project work.
Sarai Dery Folkestad
Eirik Kjønstad
27
References
[1] Andrew L. D. The role of carbon in fuel cells. Journal of Power Sources,
156(2):128 – 141, 2006.
[2] O. Haas, J. M. Simon, and S. Kjelstrup. Surface self-diffusion and mean displacement of hydrogen on graphite and a pem fuel cell catalyst support. The
Journal of Physical Chemistry C, 113(47):20281–20289, 2009.
[3] R. Str¨
obel, J. Garche, P. T. Moseley, L. J¨orissen, and G. Wolf. Hydrogen
storage by carbon materials. Journal of Power Sources, 159(2):781–801, 9 2006.
[4] J. Burress, M. Kraus, M. Beckner, R. Cepel, G. Suppes, C. Wexler, and
P. Pfeifer. Hydrogen storage in engineered carbon nanospaces. Nanotechnology,
20(20):204026, 2009.
[5] F. M¨
uller-Plathe. A simple nonequilibrium molecular dynamics method for calculating the thermal conductivity. The Journal of chemical physics, 106:6082,
1997.
[6] S. Kjelstrup. Non-equilibrium thermodynamics for engineers. World Scientific,
Hackensack, N.J., 2010.
[7] S. Skogestad. Chemical and energy process engineering. CRC Press, Boca Raton, Fla., 2009.
[8] A. R. Leach. Molecular modelling: principles and applications. Prentice Hall,
Harlow, 2nd edition, 2001.
[9] P. K. Schelling, S. R. Phillpot, and P. Keblinski. Comparison of atomic-level
simulation methods for computing thermal conductivity. Physical Review B,
65(14):144306–, 04 2002.
[10] S. Plimpton, P. Crozier, and A. Thompson. Lammps-large-scale
atomic/molecular massively parallel simulator. Sandia National Laboratories, 2007.
[11] M. Tuckerman, B. J. Berne, and G. J. Martyna. Reversible multiple time scale
molecular dynamics. The Journal of chemical physics, 97(3):1990, 1992.
[12] J. B. Nicholas, F. R. Trouw, J. E. Mertz, L. E. Iton, and A. J Hopfinger. Molecular dynamics simulation of propane and methane in silicalite. The Journal of
Physical Chemistry, 97(16):4149–4163, 1993.
[13] A. C. T. Van Duin, J. M. A. Baas, and B. Van de Graaf. Delft molecular mechanics: A new approach to hydrocarbon force fields: Inclusion of a geometrydependent charge calculation. Journal of the Chemical Society. Faraday transactions, 90(19):2881–2895, 1994.
[14] P. J. Linstrom and W. G. Mallard. Nist chemistry webbook, nist standard
reference database number 69, national institute of standards and technology,
gaithersburg md, 20899. 2003.
[15] H. M. Roder. Thermal conductivity of hydrogen for temperatures between 78
and 310 k with pressures to 70 mpa. International journal of thermophysics,
5(4):323–350, 1984.
29
List of Symbols
Notation
Description
α
˚
A−1
Parameter in the Morse potential,
see equation ref.
Concentration of component i.
mol/m3
Diffusion coefficient.
cm2 /s
Parameter in the Morse potential,
eV
see equation ref.
Dimension.
−
Transfer of kinetic energy
J
Time between swaps of velocities
s
NEMD simulations
Relative error in density from
comparison of NEMD simulations, and experimental data from
the NIST database.
Parameter in the Lennard-Jones
eV
potential.
Molar enthalpy of component i.
J/mol
Mass flux of component i.
kg/sm2
Heat flux.
W/m2
Measurable heat flux.
W/m2
Entropy flux.
J/Km2
Boltzmann’s constant.
J/K
Phenomenological coefficients
Thermal conductivity of H2 .
W/mK
Mass.
kg
Chemical potential of component
J/mol
i
Chemical potential of component
J/mol
i at constant temperature.
Moles of component i.
mol
moles
mol
Cross-sectional area.
m2
Pressure.
Pa
Gas constant.
l bar/K mol
Position vector.
−
˚
A
Distance
between
hydrogen
atoms in the hydrogen molecule.
˚
A
Parameter in the Morse potential.
˚
Cut-off distance.
A
ci
D
De
d
∆Ek
δt
ǫ
ǫLJ
Hi
Ji
Jq
Jq′
Js
kB
Lij
λ
m
µi
µi,T
Ni
n
Ω
p
R
r
r
r0
rC
Unit
30
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
s
S
Si
σ
σLJ
T
t
t0
u
U
V
VLJ
VM
Vnb
v
vh , vc
x
Z
Entropy density.
J/Km3
Entropy.
J/K
Molar entropy of component i.
J/K mol
Entropy production.
J/Km3 s
˚
A
Parameter in the Lennard-Jones
potential.
Temperature.
K
Time.
s
Time to reach convergency.
s
Internal energy density.
J/m3
Internal energy.
J
Volume.
m3
Lennard-Jones potential.
eV
Morse potential.
eV
Truncated Lennard-Jones poteneV
tial.
Velocity vector.
−
Velocity of particle in hot
m/s.
slab/cold slab.
Distance along the x-axis
m
Compressibility factor.
-
31
Appendix A. Compressibility factor
Table A.1. Compressibility factor, Z, from the EMD
simulations, listed with the average temperature Tav ,
pressure, p, and density, ρ.
Tav [K] ρ[g/cm3 ] p[bar]
Z
100.0
100.0
100.0
100.0
100.0
0.0135
0.0179
0.0224
0.0269
0.0314
50.6
63.1
73.9
83.5
91.8
0.966
0.905
0.847
0.797
0.752
150.0
150.0
150.0
150.0
150.0
0.0135
0.0179
0.0224
0.0269
0.0314
87.7
115.4
142.6
170.0
198.1
1.118
1.103
1.090
1.083
1.081
201.0
200.0
200.0
200.0
200.0
0.0135
0.0179
0.0224
0.0269
0.0314
123.6
164.1
206.6
249.8
294.4
1.181
1.176
1.184
1.193
1.205
32
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
Appendix B. MATLAB scripts
% EMD script
% Conductivity
kondData = importdata(’nvt.out’);
conductivityVector = kondData(501:length(kondData(:,1)),5);
averageConductivity = sum(conductivityVector)/length(conductivityVector);
fprintf(’Average conductivity: %.8f \n’, averageConductivity);
plot(kondData(501:length(kondData(:,1)),1),conductivityVector,’-k’);
% Pressure
pressureVector = kondData(501:length(kondData(:,1)),3);
averagePressure = sum(pressureVector)/length(pressureVector);
fprintf(’Average pressure: %.5f bar \n’, averagePressure);
figure(2)
plot(kondData(501:length(kondData(:,1)),1),pressureVector);
% Diffusion
timeSteps = kondData(:,1); % Time in femtoseconds
meanSquaredDisplacement = kondData(:,6); % mean squared displacement in
Regresjon = polyfitn(timeSteps, meanSquaredDisplacement,1);
figure(4);
plot(timeSteps, meanSquaredDisplacement,’-k’);
figure(3);
plot(X,Y);
% Temperature
temperatureVector = kondData(501:length(kondData(:,1)),2);
averageT = sum(temperatureVector)/length(temperatureVector);
33
fprintf(’Diffusion coefficient: %.6f \n’, Regresjon.Coefficients(1)/6);
fprintf(’Average temperature: %.2f \n \n’, averageT);
34
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
clear all;
% NEMD script
% Import the temperature profile file.
data = importdata(’temp_h2.profile’);
% The temperature column is collected
temperatur = data(:,4);
% Allocation of average temperature-vector
averageT = zeros(40,1);
%
i
j
k
Set counters for reading through the data.
= 1;
= 1;
= 1;
count = 1;
h = 1;
while 41*h < length(temperatur)
count = count + 1;
h = h + 1;
end
n = count;
bin = zeros(40,n);
while i < 41
j = 1;
l = k;
while l < length(temperatur)
bin(i,j) = temperatur(l);
j = j + 1;
l = l + 41;
end
k = k + 1;
averageT(i) = sum(bin(i,:))/length(bin(i,:));
i = i + 1;
end
% Calculate total average temperature
numberOfT = 0;
for c=1:40
numberOfT = numberOfT + length(bin(c,:));
end
35
allTvector = zeros(numberOfT,1);
index = 1;
for c=1:40
for n=1:length(bin(c,:))
allTvector(index) = bin(c,n);
index = index + 1;
end
end
averageOverallTemperature = sum(allTvector)/numberOfT;
position = zeros(40,1);
teller = 5;
for g=1:40
position(g) = teller;
teller = teller + 10;
end
% -------------------------------% Density, similar formatting of file, so a repeat of the above code.
densityData = importdata(’density.data’);
density = densityData(:,4);
averageP = zeros(40,1);
i = 1;
j = 1;
k = 1;
count = 1;
h = 1;
while 41*h < length(density)
count = count + 1;
h = h + 1;
end
n = count;
binD = zeros(40,n);
while i < 41
j = 1;
l = k;
while l < length(density)
binD(i,j) = density(l);
36
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
j = j + 1;
l = l + 41;
end
k = k +1;
averageP(i) = sum(binD(i,:))/length(binD(i,:));
i = i + 1;
end
% ------------------% Import the cumulative kinetic energy transfer and set it equal to cumKin.
kineticData = importdata(’nvt.out’);
cumKinVector = kineticData(:,2);
cumKin = cumKinVector(length(cumKinVector))*1.602*10^(-19); %J
% Linear regression.
%
Temperature, not including the regions where swapping occurs.
regr1temp = polyfitn(position(2:20),averageT(2:20),1);
regr2temp = polyfitn(position(22:40), averageT(22:40),1);
% Average slop/temperature gradient
tempGradient = (regr1temp.Coefficients(1) - regr2temp.Coefficients(1))/2;
% Average temperature difference
tempDifference = tempGradient*180;
% Coefficients of determination
coeffOfDetermination = [regr1temp.R2 regr2temp.R2];
% Calculate the thermal conductivity
time = 10^(-15)*10^6; %[s], 10^(-15) ps, 10^6 steps, 10^(-12) s/ps
area = 100*100*10^(-20); %[m^2]
distance = 180*10^(-10); %[m]
thermalConductivity = -(cumKin/(2*time))/(area*tempGradient*10^(10));
% -----------------------------% -----------------------------% Plots of temeperature profile (Figure 1), density profile (Figure 2) and
% temperature-density plot (Figure 3).
figure(1);
37
hold on
plot(position,averageT,’.k’);
plot(position(2:20),regr1temp.Coefficients(1)*position(2:20) + regr1temp.Coeff
plot(position(22:40), regr2temp.Coefficients(1)*position(22:40) + regr2temp.Co
hold off
figure(2);
plot(position, averageP, ’.k’);
% Local equilibrium asssumption valid? Plot vs EXPERIMENTAL DATA
% Read data file for density-temperature experimental data.
RhoTData = importdata(’densityTemperatureData.data’); %NIST
% INPUT FILE MUST BE CHANGED DEPENDING ON THE TEMPERATURE AND PRESSURE
% Collect the important data
Temperature = RhoTData(:,1); % K
Density = RhoTData(:,3); % g/cm^3
figure(3);
hold on
plot(Temperature, Density, ’-k’);
plot(averageT, averageP,’.k’);
hold off
% Calculation of average deviation from experimental data.
% For each temperature, find the closest data point in the experimental
% data.
% Left-hand side.
deviationVector = zeros(40,1);
for dataPoint = 2:20
done = 0;
closestIndex = 1;
while (done == 0)
if (averageT(dataPoint) > Temperature(closestIndex))
closestIndex = closestIndex + 1;
else
deviationVector(dataPoint) = (averageP(dataPoint)-Density(closestI
done = 1;
end
end
end
38
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
for dataPoint = 22:40
done = 0;
closestIndex = 1;
while (done == 0)
if (averageT(dataPoint) > Temperature(closestIndex))
closestIndex = closestIndex + 1;
else
deviationVector(dataPoint) = (averageP(dataPoint)-Density(closestI
done = 1;
end
end
end
% Average relative error
averageDeviation = sum(deviationVector)/38;
% -------------------
% Printing the relevant data to the console
fprintf(’Thermal conductivity: %.8f \n Temperature difference: %.1f \n Average
fprintf(’Coefficient of determination (left): %.4f \n’, coeffOfDetermination(1
fprintf(’Coefficient of determination (right): %.4f \n’, coeffOfDetermination(
fprintf(’Average deviation: %.4f \n \n’, averageDeviation);
% --------------------
39
Appendix C. LAMMPS scripts
# Input file for H2, EMD
boundary p p p
restart 100000 restart
#Simulation Units and Simulation Type
units metal
atom_style full
#Reads configuration from file
#read_restart restart.last
read_data init.data
#Mass relative mass of atoms
mass 1 1.08000
bond_style morse
bond_coeff 1 4.747 1.946 0.7414
pair_style
lj/cut 12.0
#kspace_style
pppm 5.0e-5
#kspace_modify
mesh 50 50 50 order 5
#kspace_modify slab 3.0
pair_coeff
1 1 0.002383109 2.63984
thermo_style
multi
neigh_modify delay 0 every 1 check yes
#Applies or "fixes" NVT MD
#velocity
all create 100.0 234324324 dist gaussian mom yes rot yes
group h2 type 1
variable T equal 100.0
###set NVT
fix
1 all nvt temp $T $T 1.0
##compute temp over x axe
#dump
1 all xyz 500 h2.xyz
#dump 1 all custom 500 co2.xyz element x y z
#dump_modify 1 element H
#dump_modify
1 unwrap yes
40
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
#msd calculation
compute
mymsd h2 msd com yes
variable
msdx equal "c_mymsd[1]"
variable
msdy equal "c_mymsd[2]"
variable
msdz equal "c_mymsd[3]"
variable
msdtot equal "c_mymsd[4]"
variable
timestep equal "step"
fix
mymsd h2 ave/time 1 1 1 c_mymsd[4] file h2msd.dat
#fix
2 all print 10 "${timestep} ${msdx} ${msdy} ${msdz} ${msdtot}
# -------------- Flux calculation in nve --------------variable s equal 10
variable p equal 200
variable d equal $p*$s
variable
dt equal 0.001
variable V equal vol
variable
variable
variable
variable
variable
kB equal 1.3806504e-23
# [J/K] Boltzmann
kCal2J equal 1.602177e-19
A2m equal 1.0e-10
fs2s equal 1.0e-12
convert equal ${kCal2J}*${kCal2J}/${fs2s}/${A2m}
reset_timestep 0
compute
myKE all ke/atom
compute
myPE all pe/atom
compute
myStress all stress/atom virial
compute
flux all heat/flux myKE myPE myStress
variable
Jx equal c_flux[1]/vol
variable
Jy equal c_flux[2]/vol
variable
Jz equal c_flux[3]/vol
fix
JJ all ave/correlate $s $p $d &
c_flux[1] c_flux[2] c_flux[3] type auto file J0Jt.dat ave running
variable
scale equal ${convert}/${kB}/$T/$T/$V*$s*${dt}
variable
k11 equal trap(f_JJ[3])*${scale}
variable
k22 equal trap(f_JJ[4])*${scale}
variable
k33 equal trap(f_JJ[5])*${scale}
fix
variable
density h2 ave/spatial $s $p $d x lower 5.0 density/mass file densit
k equal (v_k11+v_k22+v_k33)/3.0
#Dumps output files and thermodynamic properties
thermo_style custom step temp press pe v_k c_mymsd[4]
thermo
2000
#Timestep is 1 fs (0.001 ps)
41
timestep 0.001
#Number of time steps
run
1000000
42
SARAI DERY FOLKESTAD, EIRIK KJØNSTAD
#Input file H2, NEMD
boundary p p p
restart 100000 restart
units metal
atom_style full
read_restart restart.last
#read_data
init.data
#Mass relative mass of atoms
mass
1 1.08000
bond_style
bond_coeff
pair_style
pair_coeff
morse
1 4.747 1.946 0.7414
lj/cut 12.0
1 1 0.002383109 2.63984
thermo_style multi
neigh_modify delay 0 every 1 check yes
group h2 type 1
# Fixes
variable T equal 100
fix 1 h2 nve
#fix 1 all nvt temp $T $T 1.0
#velocity
all create $T 234324324 dist gaussian mom yes rot yes
# NEMD setup
# Fix thermal/conductivity swap method (Muller-Plathe)
fix heatSwap h2 thermal/conductivity 400 x 40 swap 1
# Compute temperature in regions
compute ke h2 ke/atom
variable temp atom c_ke/(1.5*8.61e-5)
fix temp_profile h2 ave/spatial 5 10 100 x lower 10.0 v_temp file temp_h2.pro
#dump 1 all xyz 1000 h2.xyz
#dump_modify 1 element H
fix
density h2 ave/spatial 5 10 100 x lower 10.0 density/mass file de
thermo_style custom step f_heatSwap
43
thermo 2000
timestep 0.001
run 1000000

Download Report