Note: Local thermal conductivities from boundary driven non-equilibrium molecular dynamics simulations F. Bresme and J. Armstrong Citation: The Journal of Chemical Physics 140, 016102 (2014); doi: 10.1063/1.4858434 View online: http://dx.doi.org/10.1063/1.4858434 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/1?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Thermal conductivity of carbon dioxide from non-equilibrium molecular dynamics: A systematic study of several common force fields J. Chem. Phys. 141, 134504 (2014); 10.1063/1.4896965 Thermal conductivity of carbon nanotube—polyamide-6,6 nanocomposites: Reverse non-equilibrium molecular dynamics simulations J. Chem. Phys. 135, 184905 (2011); 10.1063/1.3660348 The Thermal Conductivity of Amorphous Polymers Calculated by NonEquilibrium Molecular Dynamics Simulation AIP Conf. Proc. 982, 486 (2008); 10.1063/1.2897842 Equilibrium and nonequilibrium molecular dynamics simulations of the thermal conductivity of molten alkali halides J. Chem. Phys. 126, 204511 (2007); 10.1063/1.2734965 A critical comparison of equilibrium, non-equilibrium and boundary-driven molecular dynamics techniques for studying transport in microporous materials J. Chem. Phys. 115, 8112 (2001); 10.1063/1.1407002
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 141.212.109.170 On: Thu, 18 Dec 2014 18:32:30
THE JOURNAL OF CHEMICAL PHYSICS 140, 016102 (2014)
Note: Local thermal conductivities from boundary driven non-equilibrium molecular dynamics simulations F. Bresme1,2,a) and J. Armstrong1,b) 1
Department of Chemistry, Chemical Physics Section, Imperial College London, London SW7 2AZ, United Kingdom 2 Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway
(Received 4 November 2013; accepted 12 December 2013; published online 7 January 2014) [http://dx.doi.org/10.1063/1.4858434] The quantification of the thermal transport properties of fluids and materials is a problem of scientific and technological interest. In particular, the optimization of heat management in microelectronics1 or the design of efficient energy conversion thermoelectrics,2 requires a good understanding of the variables controlling the thermal conductivity (TC). Computer simulations, particularly non-equilibrium molecular dynamics (NEMD),3–8 offers an excellent approach to quantify the TCs of solids and liquids from a microscopic approach by using accurately parametrized force-fields. NEMD simulations rely on the use of very large thermal gradients ∼1010 K/m. Although this has raised doubts about the linearity of the system’s response, we note that in most NEMD simulations the particles move in an environment that is isotropic to a large extent and in the diffusive regime, as can be judged from the relation ∇Ta/T 1,9 where ∇T is the thermal gradient, a the characteristic molecular dimensions, and T the average temperature of the system. The Green-Kubo (GK) approach on the other hand, does not require explicit gradients or other perturbations in the equations of motion, and the TC can be computed from the heat flux autocorrelation function. Despite the differences between NEMD and GK approaches, good agreement between them has been reported in simple liquids,10 solids,11 and water.12, 13 We note that, often, the NEMD TCs are obtained by performing a linear regression of the thermal gradient over a wide range of temperatures. Because the temperature gradient induces a density gradient, the resulting TCs will be an average over different thermodynamic states, which depending on the strength of the gradient can cover a wide range of densities and temperatures. Possible discrepancies between NEMD and GK calculations may have been masked by the large uncertainties in the equilibrium and non-equilibrium computations, ∼10%.10 In this Note, we make use of the nonlinear dependence of the temperature profile to quantify the “local” TCs of a molecular fluid. The idea of locality or local equilibrium (LE), is one of the main hypotheses of non-equilibrium thermodynamics.14 This hypothesis has been corroborated in NEMD computations of thermodynamic properties and the equation of state of simple liquids, molten salts, and molecular liquids.6, 13, 15 Bridgman16 established a connection a) [email protected] b) [email protected]
0021-9606/2014/140(1)/016102/2/$30.00
between the TC and the thermodynamic properties. Hence, if LE is fulfilled by the thermodynamic properties, so too should the TC. We thus expect that the thermal gradients obtained in NEMD simulations will change locally, following the thermal conductivity corresponding to specific density/temperature states. We explore this idea in the following using NEMD simulations, and we show how to extract local thermal conductivities from the analysis of the thermal gradients. We compare our NEMD results with equilibrium GK computations. We have chosen as a model system a supercritical molecular fluid consisting of diatomic molecules, where the two atomic sites have the same size and mass. We use the truncated and shifted Lennard-Jones (LJ) potential to model the interactions, using as cut-off rc = 3.2σ , where σ is the atomic diameter. The interatomic interaction strength is defined by ε. In addition, we used a harmonic potential to model the intramolecular interactions between the atoms in the same molecule, ub (r) = k/2(r − r0 )2 , with a force constant k = 400 ε/σ 2 and molecular bond length, r0 /σ = 0.5. The NEMD simulations were performed by using thermostatting layers located at the edges and center of the simulation box. The width of the layers was set to 3σ . The simulation cell was fully periodic and had dimensions (18 × 18 × 140)σ 3 . A typical simulation consisted of 22 680 molecules. The thermostats were applied every time step, δt∗ = (ε/(mσ 2 ))1/2 δt = 0.0025, on any atoms lying inside the regions, and the momentum was reset after each thermostatting event. A typical simulation involved 1.8×107 steps of production run. The GK equilibrium simulations were performed in cubic boxes containing 2500 molecules. A careful pre-equilibration in the NVE ensemble ensured that the run temperatures where within 0.5% of the target temperatures. The time step in this case was set to δt∗ = 0.001, to enable an accurate integration of the autocorrelation functions. A typical GK-NVE simulation involved 5 × 107 steps to ensure good statistics and convergence of the heat flux autocorrelation functions. All the simulations, NEMD and GK, were performed with LAMMPS.17 We show in Figure 1(top) a representative temperature profile, as well as the corresponding equation of state, constructed using the pairs of density and temperature obtained from a single NEMD simulation. These results are compared with simulation data obtained from equilibrium
140, 016102-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 141.212.109.170 On: Thu, 18 Dec 2014 18:32:30
016102-2
F. Bresme and J. Armstrong
J. Chem. Phys. 140, 016102 (2014)
4
T (ε/KB)
4 3
T (ε/KB)
2
3 1
0
25
50
75
100
125
z-axis (σ)
2
1 0.6
1
0.8
ρ
1.2
*
10 125
10
100
8
9 8
6
50
4
25
2
*
λ
λ
*
7
75
6
1
0.5
1.5
2
*
t
5 4 3 1
1.5
2
2.5
T
3
3.5
4
*
FIG. 1. (Top) The main plot shows the equation of state obtained from NEMD simulation (blue circles) and data points of the NpT equilibrium simulations (red triangles). The inset plot shows the temperature profile for the NEMD simulation. (Bottom) Local TC obtained from NEMD (solid black line, with dashed black lines representing the error bars). The TC from GK are shown as red circles and the inset plot shows the energy flux autocorrelation function (solid red line) for ρ ∗ = (N/V )σ 3 = 1.205 (N being the number of atoms) and T∗ = kB T/ε = 1.5, and its running integral (dashed blue line). The shadowed region indicates the time interval used to obtain the thermal conductivity. For the time discretization of the correlation functions, we employed 2000 points with a step of δt∗ . Both the TC and the heat flux are represented in LJ reduced units. The errors for the NEMD and GK data are
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 141.212.109.170 On: Thu, 18 Dec 2014 18:32:30
THE JOURNAL OF CHEMICAL PHYSICS 140, 016102 (2014)
Note: Local thermal conductivities from boundary driven non-equilibrium molecular dynamics simulations F. Bresme1,2,a) and J. Armstrong1,b) 1
Department of Chemistry, Chemical Physics Section, Imperial College London, London SW7 2AZ, United Kingdom 2 Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway
(Received 4 November 2013; accepted 12 December 2013; published online 7 January 2014) [http://dx.doi.org/10.1063/1.4858434] The quantification of the thermal transport properties of fluids and materials is a problem of scientific and technological interest. In particular, the optimization of heat management in microelectronics1 or the design of efficient energy conversion thermoelectrics,2 requires a good understanding of the variables controlling the thermal conductivity (TC). Computer simulations, particularly non-equilibrium molecular dynamics (NEMD),3–8 offers an excellent approach to quantify the TCs of solids and liquids from a microscopic approach by using accurately parametrized force-fields. NEMD simulations rely on the use of very large thermal gradients ∼1010 K/m. Although this has raised doubts about the linearity of the system’s response, we note that in most NEMD simulations the particles move in an environment that is isotropic to a large extent and in the diffusive regime, as can be judged from the relation ∇Ta/T 1,9 where ∇T is the thermal gradient, a the characteristic molecular dimensions, and T the average temperature of the system. The Green-Kubo (GK) approach on the other hand, does not require explicit gradients or other perturbations in the equations of motion, and the TC can be computed from the heat flux autocorrelation function. Despite the differences between NEMD and GK approaches, good agreement between them has been reported in simple liquids,10 solids,11 and water.12, 13 We note that, often, the NEMD TCs are obtained by performing a linear regression of the thermal gradient over a wide range of temperatures. Because the temperature gradient induces a density gradient, the resulting TCs will be an average over different thermodynamic states, which depending on the strength of the gradient can cover a wide range of densities and temperatures. Possible discrepancies between NEMD and GK calculations may have been masked by the large uncertainties in the equilibrium and non-equilibrium computations, ∼10%.10 In this Note, we make use of the nonlinear dependence of the temperature profile to quantify the “local” TCs of a molecular fluid. The idea of locality or local equilibrium (LE), is one of the main hypotheses of non-equilibrium thermodynamics.14 This hypothesis has been corroborated in NEMD computations of thermodynamic properties and the equation of state of simple liquids, molten salts, and molecular liquids.6, 13, 15 Bridgman16 established a connection a) [email protected] b) [email protected]
0021-9606/2014/140(1)/016102/2/$30.00
between the TC and the thermodynamic properties. Hence, if LE is fulfilled by the thermodynamic properties, so too should the TC. We thus expect that the thermal gradients obtained in NEMD simulations will change locally, following the thermal conductivity corresponding to specific density/temperature states. We explore this idea in the following using NEMD simulations, and we show how to extract local thermal conductivities from the analysis of the thermal gradients. We compare our NEMD results with equilibrium GK computations. We have chosen as a model system a supercritical molecular fluid consisting of diatomic molecules, where the two atomic sites have the same size and mass. We use the truncated and shifted Lennard-Jones (LJ) potential to model the interactions, using as cut-off rc = 3.2σ , where σ is the atomic diameter. The interatomic interaction strength is defined by ε. In addition, we used a harmonic potential to model the intramolecular interactions between the atoms in the same molecule, ub (r) = k/2(r − r0 )2 , with a force constant k = 400 ε/σ 2 and molecular bond length, r0 /σ = 0.5. The NEMD simulations were performed by using thermostatting layers located at the edges and center of the simulation box. The width of the layers was set to 3σ . The simulation cell was fully periodic and had dimensions (18 × 18 × 140)σ 3 . A typical simulation consisted of 22 680 molecules. The thermostats were applied every time step, δt∗ = (ε/(mσ 2 ))1/2 δt = 0.0025, on any atoms lying inside the regions, and the momentum was reset after each thermostatting event. A typical simulation involved 1.8×107 steps of production run. The GK equilibrium simulations were performed in cubic boxes containing 2500 molecules. A careful pre-equilibration in the NVE ensemble ensured that the run temperatures where within 0.5% of the target temperatures. The time step in this case was set to δt∗ = 0.001, to enable an accurate integration of the autocorrelation functions. A typical GK-NVE simulation involved 5 × 107 steps to ensure good statistics and convergence of the heat flux autocorrelation functions. All the simulations, NEMD and GK, were performed with LAMMPS.17 We show in Figure 1(top) a representative temperature profile, as well as the corresponding equation of state, constructed using the pairs of density and temperature obtained from a single NEMD simulation. These results are compared with simulation data obtained from equilibrium
140, 016102-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 141.212.109.170 On: Thu, 18 Dec 2014 18:32:30
016102-2
F. Bresme and J. Armstrong
J. Chem. Phys. 140, 016102 (2014)
4
T (ε/KB)
4 3
T (ε/KB)
2
3 1
0
25
50
75
100
125
z-axis (σ)
2
1 0.6
1
0.8
ρ
1.2
*
10 125
10
100
8
9 8
6
50
4
25
2
*
λ
λ
*
7
75
6
1
0.5
1.5
2
*
t
5 4 3 1
1.5
2
2.5
T
3
3.5
4
*
FIG. 1. (Top) The main plot shows the equation of state obtained from NEMD simulation (blue circles) and data points of the NpT equilibrium simulations (red triangles). The inset plot shows the temperature profile for the NEMD simulation. (Bottom) Local TC obtained from NEMD (solid black line, with dashed black lines representing the error bars). The TC from GK are shown as red circles and the inset plot shows the energy flux autocorrelation function (solid red line) for ρ ∗ = (N/V )σ 3 = 1.205 (N being the number of atoms) and T∗ = kB T/ε = 1.5, and its running integral (dashed blue line). The shadowed region indicates the time interval used to obtain the thermal conductivity. For the time discretization of the correlation functions, we employed 2000 points with a step of δt∗ . Both the TC and the heat flux are represented in LJ reduced units. The errors for the NEMD and GK data are
