Thermal conductivity of carbon dioxide from non-equilibrium molecular dynamics : A systematic study of several common force fields

We report a systematic investigation of the thermal conductivity of various three-site models of carbon dioxide (CO2) using nonequilibrium molecular dynamics in the temperature range 300-1000 K and for pressures up to 200 MPa. A direct comparison with experimental data is made. Three popular CO2 force fields (MSM, EPM2, and TraPPE) and two flexible models (based on EPM2) were investigated. All rigid force fields accurately predict the equation of state for carbon dioxide for the given range of variables. They can also reproduce the thermal conductivity of CO2 at room temperature and predict a decrease of the thermal conductivity with increasing temperature. At high temperatures, the rigid models underestimate the thermal conductivity.


I. INTRODUCTION
Carbon dioxide (CO 2 ) has an important impact on the climate and is therefore widely studied.Huge efforts are being made, for instance, to reduce emissions of CO 2 to the atmosphere, by capture-and sequestration techniques. 1,2 n that context, membrane separation techniques are needed, at high as well as low temperatures. 3,4 ossil-fueled power systems, natural gas processes, or production of hydrogen gas include all high-temperature separation technologies. 2,5 he thermal conductivity of CO 2 is needed for process modelling in these processes.
Molecular simulation is a popular technique for the prediction of thermal conductivities of fluids. 6][9] Nieto-Draghi et al. used semi or fully flexible models of CO 2 to predict the thermal conductivity in the temperature range 300-400 K. 7 More recently Liang et al. showed that a CO 2 model with site-site interactions fitted from ab initio calculations provides a good prediction of the thermal conductivity at low density. 9Most of these studies focused on a narrow temperature and pressure range (around room temperature and up to 10 MPa), while experimental data of the thermal conductivity of CO 2 are also available at elevated temperatures, and for pressures up to 1000 K and 200 MPa. 10 Hence, it is important to develop further a molecular dynamics simulation model for the thermal conductivity of this important molecule.2][13][14][15][16] The quality of these rigid models of CO 2 for the prediction of thermal conductivity has not been reported, however.a) Author to whom correspondence should be addressed.Electronic mail: signe.kjelstrup@ntnu.no In this work, we report thermal conductivities of CO 2 for the temperature range 300-1000 K and for pressures up to 200 MPa.We will use three common models of CO 2 , namely, MSM, 17 EPM2, 18 and TraPPE, 19 and compare the computed conductivities with those of the National Institute of Standards and Technology database (NIST). 10The models differ in their bonds lengths of C-O and in their values of the Lennard-Jones (LJ) potential parameters, as well as in the partial charges that are used.We will show that all rigid models can correctly predict the thermal conductivity of CO 2 .The TraPPE model is slightly superior.However, at high temperatures all models underestimate the thermal conductivity of CO 2 .In order to understand why, two flexible models, based on EPM2, will also be examined.We shall see that flexibility may partially explain the discrepancy observed.
The paper is structured as follows.In Sec.II, an overview of the technical and simulation details is provided.The predictions of the equation of state are presented in Sec.III.The thermal conductivity of CO 2 at various pressures and temperatures from the various models are next compared with the experimental data from NIST.We will close the paper in Sec.IV with comments and conclusions.

II. SIMULATION DETAILS
All models for CO 2 studied here are 3-site models with fixed C=O bond length and a fixed angle of 180 • .Only intermolecular interactions are needed to describe the system.The intermolecular potential consists of long-range Coulombic interactions, and a shifted and truncated 12-6 LJ potential 20  1).
where r ij is the distance between atoms i and j, ε ij and σ ij are LJ potential parameters, and r c is the cutoff radius.The LJ interaction parameters between different types of atoms were calculated from the Lorentz-Berthlot mixing rules 20 The Coulombic interactions equal where q i , q j are the partial charges on atoms i, j, and ε 0 is the dielectric constant of vacuum.In our work, we used the particle-particle particle-mesh solver implemented in LAMMPS 21 for electrostatic interactions, see Ref. 22 for more details.The force field parameters for the MSM, 17 EPM2, 18 and TraPPE 19 models are listed in Table I.They were determined at room temperature to reproduce liquid vapor equilibrium of CO 2 .
For the fully flexible models (EPM2_Flex1, EPM2_flex2), additional functions were used to describe bond stretching (harmonic potential Eq. (7) or Morse potential Eq. ( 8)) and angle bending of CO 2 (Eq.( 9)) where r ij is the distance between atom i and j; θ ijk is the angle between atoms i,j,k; k S and k B are the force constant.The nonbonding parameters for the MSM, 17 EPM2, 18 and TraPPE 19 models force fields are also listed in Table I.For flexible models k B was taken from Harris and Yung. 18In EMP2_flex1, k S was taken from Nieto-Draghi et al. 7 using harmonic potential for C-O bond stretching; however, the harmonic equation (7) was not able to describe the separation of atoms at longer distance. 6Hence, a Morse potential for EPM2_flex2 was fitted with Eq. ( 8) using quantum chemistry data to overcome this limitation of harmonic potential.We used Gaussian 09 package 23 to perform a Density Functional Theory (DFT) calculation with B3LYP functional 24,25 and a full-electron 6-311+G(d,p) basis set. 26Geometry optimization and potential energy scan was made to fit with a Morse potential (Figure 1).The cut-off radius was 12 Å, which corresponds to r cut ∼ 4σ O .The tail correction was not used for this cut-off.An increase of the cut-off value in our simulations, did not improve the accuracy of the computed thermal conductivity of CO 2 , but significantly increased the computational cost. 27It was observed by Bugel and Galliero for Lennard-Jones fluids that r cut = 2.5σ was sufficient for reliable results for thermal conductivity. 28he thermal conductivity can be either obtained from equilibrium molecular dynamics (EMD) or from nonequilibrium molecular dynamics (NEMD) simulations. 20The Green-Kubo formulations are commonly used in EMD to access the thermal conductivity.There are several NEMD techniques to obtain a heat flux and the corresponding temperature gradient, from which the thermal conductivity is computed.The non-equilibrium situation can be obtained by swapping particle momenta, 29,30 by using a heat exchange algorithm (HEX), 31 or by thermostating the boundaries. 32,33 n this work, we applied the thermostating technique, which was used earlier to successfully calculate the thermal conductivity of hydrocarbons 32,34 in zeolites and water. 33Here, we summarize the essentials of the simulation technique, previously FIG. 1. Energy potential of CO 2 bond stretching calculated by DFT B3LYP/6-311+G(d,p).The continuous line represents a harmonic potential Eq. ( 7) by Nieto-Draghi et al. 7 and Morse potential Eq. ( 8) (this work).The parameters are in Table I.
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: 131.180.130.205On: Wed, 05 Nov 2014 13:07:26 described in Refs.32 and 33.We took a simulation box size of 20 × 3 × 3 (nm 3 ) with periodic boundary conditions in three dimensions.The total number of CO 2 molecules was in the range 200-3000 to cover a wide range of densities, from 80 to 1200 kg/m 3 .In total, 300 simulations were performed to compute the thermal conductivity of CO 2 in the temperature range 300 K-1000 K for the four different CO 2 models.The time step for integration of the equations of motion was 1 fs.The initial configuration was constructed by randomly distributing the CO 2 .The system was stabilized during 1 ns by NVT runs with the Nosé-Hoover thermostat. 35When the system was in thermal equilibrium, we performed NEMD simulations of 5 ns runs in the microcanonical ensemble (NVE) with thermostats 6 at the hot and cold regions.A longer run of 7.5 ns was performed at high temperature (600-1000 K).The heat flux and temperature gradient were checked for convergence.The average values of temperature and pressure in NEMD simulations were within 1% of expected values.The last 2.5 ns of the run was used to determine the temperature gradient and the heat flux through the system.This is sufficiently long to obtain sufficient statistics and consistent trajectories.The simulation temperature was set as T = 1 2 (T cold + T hot ).The simulation box was divided into 20 small equal slabs.The cold and hot regions were chosen in slabs 1-2 and slabs 11-12, respectively.The volumes of the cold and hot regions are the same along the x-direction of the simulation box (Figure 1).In the NEMD simulation, the average temperature of each slab was recorded.The temperature gradient was obtained by fitting all average temperatures along the x-axis to a straight line, excluding the thermostat regions.
The temperature was maintained in each thermostat by supply or withdrawal of kinetic energy.The total energy of the system is unchanged, meaning that the energy withdrawn is the energy supplied.The heat (energy) flux through the system can therefore be computed from the change in kinetic energy K in any of the two thermostats during a single time step where J q is the heat flux through the simulation box, δt is the time step, and A is the cross sectional area.The factor 2 arises from the fact that due to the periodic boundary conditions that there are two temperature gradients in the symmetric simulation box.The thermal conductivity is obtained from Fourier's law ( 11) In order to compare with experimental data, we used the root mean squared error of a model, defined as where λ sim , λ exp , and k are the calculated thermal conductivity, experimental thermal conductivity, and number of simulations, respectively.The minimum error (%) is defined as FIG. 2. The NEMD simulations to apply temperature gradient in the simulation box.See text for details.
The relative error is defined as

III. RESULTS AND DISCUSSIONS
We investigated first 120 state points, in order to generally assess the CO 2 -models.The average temperature in the simulation box varied from 300 K to 1000 K, while the pressure was varied up to 200 MPa.The models were first tested for their accuracy in the prediction of the equation of state for CO 2 .Figure 2 shows the pressure and density obtained for various temperatures together with the NIST experimental data. 10All models reproduced very well the thermodynamic equation of state, at low as well as high temperatures and pressures.Calculated points fell almost on top of each other.In spite of the CO 2 models being developed to fit experimental data below the critical temperature (for vapor-liquid phase equilibria of pure CO 2 and CO 2 in mixtures with hydrocarbons), a very good prediction was found for the equation of state of CO 2 , also at higher temperatures.The flexible model created here, gives an equally good prediction at low pressure, but slightly better at high pressure (Figure 3).An average error (defined by Eq. ( 13)) is 6.3% and 3.6% for rigid and flexible models, respectively.The NEMD simulations to calculate the thermal conductivity of CO 2 were therefore done with all models.The microstructure of CO 2 was previously studied by ab initio molecular dynamics and flexible model. 36The radial distribution functions and distribution of bond and angle of CO 2 obtained by flexible models (see the supplementary material 27 ) are in agreement with literature. 36At high temperature 1000 K, the bond angles are more flexible than that at 300 K.However, both bond length and angle distribution obtained by flexible classical force field are more rigid than that obtained by ab initio results. 36igure 4 shows the average temperature and density profile of a typical run at T = 300 K with the TraPPE force field.We obtain similar profiles for the MSM, EPM2, and flexible models (not shown).The temperature profile was fitted to a straight line positioned in the analysis layers (discarding the hot and cold regions).A good accuracy fit was obtained for all temperatures (regression coefficient R 2 > 0.90).The temperature gradient was typical chosen in range 3-5 K nm −1 .Larger gradients have been selected earlier (e.g., 15-20 K(nm) −1 for water 33 ).We verified that the magnitude of the temperature gradient does not influence the value of the thermal conductivity. 28The intermolecular potential parameter of the model is the most important factor in its determination.Typical results from the TraPPE model at 300 K and various densities are listed in Table II.The error bar (the estimate is based on the error bar of heat flux and temperature gradient) of the calculated thermal conductivity is maximum 5 mW m −1 K −1 .The simulated value predicted the TABLE II.Thermal conductivity of carbon dioxide from the TraPPE model at 300 K as a function of density.The symbols λ sim , λ exp denote the thermal conductivity from simulations and from the NIST data of Ref. 10, respectively.T is the temperature gradient.J q is the heat flux through the system.experimental result within this accuracy at medium densities.The deviation between the two was relatively larger at low densities.This has also been observed by others. 33However, the thermal conductivity of CO 2 is small at low densities, so the relative deviation becomes large for this reason.
The simulated thermal conductivities of CO 2 at 300 and 400 K are presented as a function of the pressure in Figure 5. Figure 5(a) for 300 K shows that all models predict well the experimental values.At 400 K (Fig. 5(b)), the computed values underestimate the experimental thermal conductivity of CO 2 .This tendency is strengthened as the temperature rises, see Figures 6 and 7.All rigid models underestimate the thermal conductivity in these figures.This effect is stronger at high temperatures.
This tendency is quantified by the root mean square error (RMSE) computed from Eq. ( 12) and listed in Table III.For example, for the TraPPE model at 300 K the RMSE is only 6.8 mW m −1 K −1 .However, at 1000 K the RMSE has increased 6 times for the same model.As the consequence, the minimum error increases 9 times, up to 32.1%, at 1000 K.
Two main factors may contribute to this discrepancy.In the first place, the force fields of the rigid models that we have used were obtained at temperatures below the critical temperature of carbon dioxide (304 K).In spite of the good prediction, these models of the equation of state (Figure 3) may not have the wanted effect on a transport property like the thermal conductivity.Also, a rigid model may be too limited.Even CO 2 is a small molecule and the vibrational energy of CO 2 could be neglected up to 1000 K, 9, 37 a possibility to bend or vibrate may have an impact on the simulated thermal conductivity.
To investigate the last factor further, we included two flexible models EPM2_Flex1 and EPM2_Flex2 as described  in Sec.II.The results are plotted along with the rigid model results in Figures 5-7.We see from Figs. 5-7 that the added possibility for bending and vibration in the molecule increases the thermal conductivity in all cases.On the whole, however, it leads to an overestimation of the experimental value.
The root mean square errors of all models for the whole temperature range are listed in Table III.In general, the TraPPE model is superior to the MSM or EPM2 models.At 300 K, the TraPPE model can forecast the thermal conductivity within an average error of 7 mW m −1 K −1 , while the EPM2 and MSM models show larger deviations, 15.4 and 14.3 mW m −1 K −1 , respectively.Using a semi-flexible or a full-flexible EMP2 model did not improve significantly the accuracy of the prediction. 7The results predicted by these models were still 22%-30% from the experimental value.This suggests that the force fields need be optimized at higher temperatures, and possibly also for transport properties.
The observation is to some degree supported by Liang et al. who showed that fully flexible models with site-site interaction based on ab initio potential did not capture the thermal conductivity of CO 2 at density higher than 135 kg/m 3 . 9Also, a simple rigid model was able to reproduce the complex thermal conductivity of water, even at high temperature of 700 K. 33 At low temperatures, the rigid model of water overestimated thermal conductivity. 33n general, we have observed that the TraPPE model can better predict the thermal conductivity of CO 2 than the EPM2 and MSM-models do, and that the predictions are satisfactory at 300 K and low densities.All models yield an excellent prediction of the equation of state for CO 2 , but they fail to predict experimental results at temperatures above 400 K.A clear effect of adding bond stretching and bondbending is present in Figure 8.The rigid models are superior in low temperature region (below 600 K).However, at high temperature range (above 600 K), flexible models seem to be better.Using a better description of CO 2 bonding with Morse potential (EPM2_Flex2) already shows a slightly better results than a harmonic potential (EMP2_Flex1).However, a renewed evaluation of the force fields may also be needed to conclude on how to best reproduce the conductivity of CO 2 at high temperature.Adding flexibility to a rigid model (e.g., TraPPE) may require fully re-optimize the parameters of the force field.This will be included in our future work.
TABLE III.The root mean squared error (RMSE) (in mW m −1 K −1 ) of the thermal conductivity of carbon dioxide, obtained by different rigid force field models in the temperature interval 300-1000 K.The value in the parentheses is the error (in %) from Eq. ( 13).The TraPPE model is better than the EPM2 and MSM models.

IV. CONCLUSIONS
We have employed nonequilibrium molecular dynamics simulation to explore the thermal conductivity of CO 2 of three popular rigid three-site models (TraPPE, EPM2, MSM) and two fully flexible models based on EPM2.It is remarkable that the models provide an excellent equation of state for CO 2 for low as well as high temperatures.The rigid models, especially the TraPPE model of CO 2 , can predict the thermal conductivity of CO 2 within 5% error in a wide range of pressure for 300-400 K.This model can well be used to study heat and mass transfer in mixtures of CO 2 with CH 4 or H 2 for membrane separation processes near room temperature.The underestimation of the experimental results by all rigid models at high temperatures is probably due to a lack of optimization of the force fields at these conditions.In view of the importance of the properties of carbon dioxide, this deficiency should be mended.

FIG. 3 .
FIG.3.Equation of state of CO 2 for various temperatures in the interval 300 K-1000 K. Four different CO 2 models are reported (see text for explanation).The solid lines represent experimental data taken from the NIST database.10

TABLE I .
Parameters of several CO 2 potential models used in simulations.The flexible models were based on EPM2 model with additional bond stretching and angle bending terms.The EPM2_Flex2 was fitted in this work based on quantum calculation (Figure