Thermal Conductivities of Molecular Liquids by Reverse

Jul 20, 2005 - Note that, in principle, in multicomponent mixtures, alternative definitions of the heat flux are possible.19 They contain diffusional ...
30 downloads 8 Views 112KB Size
Download PDF
15060

J. Phys. Chem. B 2005, 109, 15060-15067

Thermal Conductivities of Molecular Liquids by Reverse Nonequilibrium Molecular Dynamics Meimei Zhang,*,† Enrico Lussetti,† Luı´s E. S. de Souza,‡ and Florian Mu1 ller-Plathe† International UniVersity Bremen, P.O. Box 750 561, D-28725 Bremen, Germany, and Banco Central do Brasil, Desup/Diesp, AVenida Paulista 1804, 01310-922 Sa˜ o Paulo-SP, Brazil ReceiVed: March 9, 2005

The reverse nonequilibrium molecular dynamics method for thermal conductivities is adapted to the investigation of molecular fluids. The method generates a heat flux through the system by suitably exchanging velocities of particles located in different regions. From the resulting temperature gradient, the thermal conductivity is then calculated. Different variants of the algorithm and their combinations with other system parameters are tested: exchange of atomic velocities versus exchange of molecular center-of-mass velocities, different exchange frequencies, molecular models with bond constraints versus models with flexible bonds, united-atom versus all-atom models, and presence versus absence of a thermostat. To help establish the range of applicability, the algorithm is tested on different models of benzene, cyclohexane, water, and n-hexane. We find that the algorithm is robust and that the calculated thermal conductivities are insensitive to variations in its control parameters. The force field, in contrast, has a major influence on the value of the thermal conductivity. While calculated and experimental thermal conductivities fall into the same order of magnitude, in most cases the calculated values are systematically larger. United-atom force fields seem to do better than all-atom force fields, possibly because they remove high-frequency degrees of freedom from the simulation, which, in nature, are quantum-mechanical oscillators in their ground state and do not contribute to heat conduction.

1. Introduction Within linear-response theory, the thermal conductivity λ relates a heat flux j to a temperature gradient ∇T

j ) -λ∇T

(1)

In an isotropic medium, such as a liquid, j and ∇T are collinear and λ is scalar. It may be calculated using standard molecular dynamics techniques, be they equilibrium (Einstein or GreenKubo) or nonequilibrium. For technical as well as scientific reasons, it is however advantageous to use so-called reverse nonequilibrium methods. They are called reverse, because they reverse cause and effect with respect to experiment and conventional nonequilibrium simulation methods. In experiment, one applies a temperature gradient to the system, waits for steady state, and then measures the heat flux. In reverse nonequilibrium molecular dynamics (RNEMD), the heat flux is imposed on the system, and the resulting temperature gradient is measured. The advantages and disadvantages of RNEMD have been discussed in detail elsewhere,1 and this discussion is not to be repeated here. The RNEMD method originally had been developed for the calculation of thermal conductivities in monatomic LennardJones liquids.1 Since then, its scope has been extended to the calculation of shear viscosities of atomic fluids2,3 and Soret coefficients of mixtures of Lennard-Jones particles.4,5 The algorithm version for thermal conductivities has been used for Yukawa fluids,6 thin crystalline films,7 carbon nanotubes,8 and some molecular systems: water,9 octahydro-1,3,5,7-tetranitro† ‡

International University Bremen. Banco Central do Brasil.

1,3,5,7-tetrazocine,10 and DMSO/water mixtures.11 In the molecular calculations, some of the modifications described below have already been implemented. The algorithm version for shear viscosities has been extended and validated for atomistic models of molecular liquids12 and simple models of block copolymers,13 amphiphiles,14 and liquid crystals.15 It has also been successfully combined with the dissipatitive particle dynamics equation of motion.16 The method has recently been reviewed in detail.17 Worries about convective heat transport as a source of error are unnecessary, because of the small dimensions. Estimates from instability theory18 indicate that temperature differences of the order of 1017 K would be needed to drive convection on the scale of nanometers. In addition, for the single-component fluids, gravitation would be required to generate a RayleighBenard instability. This paper extends our RNEMD method to the thermal conductivity of molecular liquids and validates it for a series of molecular liquids and mixtures, testing algorithm variants and simulation parameters. The purpose is to help establish the range of applicability and the limitations of the method. A similar detailed investigation was previously carried out for the shear variant of RNEMD.12 As the goals are mainly technical, the reproduction of known experimental thermal conductivities is a secondary issue here. We have mostly used tested and trustable force fields for the different liquids, but made no attempt to tune them to give experimental thermal conductivities. 2. Reverse Nonequilibrium Molecular Dynamics for Atomic Fluids In this section, we briefly review the original RNEMD method for calculating the thermal conductivity of atomic liquids. It

10.1021/jp0512255 CCC: $30.25 © 2005 American Chemical Society Published on Web 07/20/2005

Thermal Conductivities of Molecular Liquids

J. Phys. Chem. B, Vol. 109, No. 31, 2005 15061

will then be easier to discuss the necessary modifications for molecular systems in the following section. For detailed discussions, see refs 1 and 17. The periodic system is partitioned into ∼20 slabs along, say, the z direction. Of these slabs, the one at the beginning of the simulation box is designated as the “hot” slab and the one at the center of the box as the “cold” slab. The nonequilibrium algorithm moves, in an artificial manner, energy from the cold slab to the hot slab, i.e., against the temperature gradient (dT/dz). Since energy is conserved, it flows back through the liquid by physical transport mechanisms, resulting in a heat flux jz. When the steady state has been reached, jz is equal in magnitude but opposite in direction to the known imposed artificial energy flow. Equation 1 then reads in this geometry

(dTdz )

jz ) -λ

(2)

Equation 2 shows that, at steady state, a temperature gradient results: A constant jz * 0 is imposed, λ is a positive constant, so dT/dz must be non-zero, too. The temperature gradient is thus the primary data, which must be extracted. Note that, in principle, in multicomponent mixtures, alternative definitions of the heat flux are possible.19 They contain diffusional mass fluxes of the components weighted with their partial molar enthalpies. In our simulations, the steady state is also reached for the concentrations of components. The concentration profiles become stationary, and there are no net fluxes of components. Hence, the additional terms vanish, and the definition of the heat flux remains unique. The RNEMD algorithm comprises two ingredients, both of which require modification when going from atomic to molecular fluids: a means of moving energy from the cold to the hot slab and a way of calculating the temperature gradient. For atomic liquids, the energy transfer mechanism uses a Maxwell Daemon, whichsevery so many MD time stepssdetermines the atomic velocities (or momenta) of all atoms in the two slabs. It then selects the hottest atom of the cold slab and the coldest atom of the hot slab and swaps their Cartesian velocity vectors. Other schemes are possible, which do not use the hottest and the coldest atoms, respectively, but two other atoms with a smaller difference in kinetic energy, and which exchange velocities more often. The details are relatively unimportant, as long as energy is moved from the cold to the hot slab. This is always possible, as the Maxwell-Boltzmann distribution of atomic velocities is much wider than the difference between the average velocities in the two slabs. It can be shown that, if velocities are exchanged between atoms of like mass, the algorithm conserves both total energy and total linear momentum up to a negligible error coming from the discretization of the equations of motion. Note that there is also a recent variant of the RNEMD schemes which allows velocity exchange between objects of different masses.20 The heat flux is given by the sum of all energy transfers per time per area

jz )

1

m

(Vhot2 - Vcold2) ∑ 2tAtransfers 2

(3)

where t is the length of the simulation, A is the cross-sectional area of the simulation box perpendicular to flow direction z, Vhot and Vcold denote the velocities of the hot and cold atoms of like mass m, respectively, whose velocities are exchanged, and the factor 2 arises from the periodicity. In our implementation, the strength of the nonequilibrium perturbation is governed by

the period W, with which velocity swaps are executed, a larger W meaning a smaller perturbation. The second task is the calculation of the temperature gradient. To this end, one analyzes the intervening slabs, where no velocity exchange is performed and where unperturbed Newtonian dynamics takes place. One evaluates the temperature Tslab within each slab as

3NkBTslab 2

1 N miVi2〉 ) 〈 2 atoms i

∑

(4)

in slab

where the number of atoms in the slab is N and kB is Boltzmann’s constant. The angle brackets denote averaging over the atoms in the slab and over time. For a sufficiently small perturbation (large W), the linear response holds and the temperature profile is linear. Its slope (dT/dz) is extracted by a linear regression. With jz known from eq 3, the thermal conductivity can be calculated. A critical issue for the calculation of transport coefficients by any equilibrium or nonequilibrium molecular dynamics method is thermostating. In principle, any thermostat generates or destroys energy in the system. Energy is no longer a conserved property, and calculated heat fluxes may lose their significance. As a result, thermostats should be avoided altogether in the calculation of thermal conductivities. (As most thermostats do not conserve momentum either, analogous considerations apply in the calculation of shear viscosities.17) However, calculations on realistic molecular fluids may necessitate a thermostat, be it because NVT conditions are desired or be it to counteract slow temperature drifts due to round-off or cutoff noise. There is evidence from calculations on LennardJones fluids that a gently coupled Berendsen thermostat does not significantly alter the calculated λ.1 This may be due to the uniform velocity scaling employed by this thermostat, which mainly shifts the entire temperature profile in a simulation, but does not destroy it. As the temperature variation within the simulation is usually much smaller than the average temperature, the influence of the thermostat on the temperature gradient and, hence, the thermal conductivity is small, too. In this respect, the otherwise widely criticized Berendsen thermostat might have an advantage over alternatives. Thermostats, which apply individual velocity scaling based on local friction, are definitely not suited. In this work, we will therefore also investigate the effect of a Berendsen thermostat on calculated thermal conductivities for realistic molecular liquids. 3. Modifications for Molecular Fluids The RNEMD method as outlined in section 2 can be used for molecular fluids, too, if the molecules are fully flexible and no constraints are present. In this case, the appropriate velocity exchanges are applied to individual atoms of the molecules. Energy and momentum conservation demand, as before, that the two atoms have the same mass. The analysis of the velocity or temperature profiles is also done on an atom-by-atom basis. This variant of the method will be called atomic exchange. For molecular models with holonomic constraints, the method has to be amended if energy and momentum conservation are to be maintained. This is described, for example, in refs 9 and 12. The description that follows is our implementation in the atomistic MD program YASP.21 The velocity exchange is now done between entire molecules of like mass, rather than between atoms. Otherwise, the atomic velocities after exchange would not comply with the constraints. The selection of the two

15062 J. Phys. Chem. B, Vol. 109, No. 31, 2005

Zhang et al.

molecules, whose velocities are to be swapped, is based on their center-of-mass momenta

∑

P)

mivi ) V

atoms i in molecule

∑

mi

(5)

atoms i in molecule

To establish a temperature gradient, the Maxwell Daemon selects the translationally hottest and coldest molecules in the cold and hot slabs, respectively. The translational kinetic energy of a molecule is defined as 1

/2V2

∑

mi

atoms i in molecule

with V its center-of-mass velocity from eq 5. The entire center-of-mass velocity vectors V of the two selected molecules are exchanged. Thus, the velocity v′i of atom i of molecule 1 after the exchange is given in terms of velocities before the exchange.

v′i ) vi - V1 + V2

(6)

As only center-of-mass velocities are exchanged, the relative velocities of atoms within each molecule remain unchanged. Hence, the velocities after the exchange are compatible with the application of a constraint algorithm. The exchange of molecular center-of-mass velocities instead of atomic velocities is called molecular exchange. It can, of course, also be applied to flexible molecules, even if it is not necessary in this case. The presence of constraints complicates the calculation of the temperature profile, too. The temperature in a molecular dynamics simulation with constraints is given by the equipartition theorem

(

)

3N - C 2

1 kBT ) 〈 mivi2〉 2 atoms i

∑

(7)

where C is the number of constraints in the system. If a local (slab) temperature is to be determined, the angle brackets denote time averaging only over the N atoms within the slab, and C refers to the number of constraints in the slab. This number can be estimated analytically only if the constraints are assumed to be distributed uniformly throughout the system. This is only the case if (i) the density is uniform (small perturbation) and (ii) the composition is uniform (only one molecular species or Soret effect negligible). In all other cases, the local C’s must be evaluated in the simulation. In YASP, bond constraints are solved by a modified SHAKE procedure.22 The constraints are stored as a list of pairs of atoms whose distance is to be kept constant. The program uses this list to determine the slab(s) in which the two atoms defining a constraint reside. For every atom, the constraint counter in its slab is incremented by 1/2. If an atom is part of c constraints, its contribution to the constraints in the slab is c/2. In this way, constraints spanning two slabs are split between the slabs. As a consequence, the scheme also allows for molecules extending over several slabs. 4. Computational Details All molecular dynamics simulations were carried out with the YASP package,21 which uses the leapfrog algorithm23 and orthorhombic periodic boundary conditions. The simulation cells were usually elongated in the z direction, in which the heat flow was imposed (Lx ) Ly ) Lz/3). Bond constraints, if present,

were solved by the SHAKE method22 to a relative tolerance of 10-6. When the temperature was kept constant, this was done by Berendsen’s method,24 with the target temperature being 300 K, unless specified otherwise. Temperatures and densities were chosen at or close to ambient conditions, where experimental data are most abundant. The intramolecular force field contained harmonic bond stretching or constraints, harmonic bond angle bending, periodic cosine-type torsional potentials (n-hexane and cyclohexane), and harmonic dihedral potentials (benzene outof-plane). The nonbonded potential includes Lennard-Jones terms, with Lorentz-Berthelot mixing rules for unlike interactions, and electrostatic interactions for the models with atomic partial charges. The latter were treated using the reaction-field method. For details of the functional form, see ref 21. Nonbonded interactions were evaluated from a Verlet neighbor list, which was updated using a link-cell method. Nonbonded interactions within a molecule were excluded, except for n-hexane (see section 5.5). Specific settings of simulation parameters are reported below for the individual systems. We sample the temperature profiles only in time steps, in which no velocity swaps are performed. Therefore, the sampling period is typically slightly different from the exchange period W or its multiples, for example W + 1. This is not specifically mentioned for the individual cases in section 5. 5. Results and Discussion 5.1. Benzene. Benzene (645 molecules, density F ) 874.9 kg/m3, T ) 307.5 ( 2.4 K) was simulated at constant NVE conditions for 900 ps using a time step of 1 fs. A modified OPLS force field (bond constraints, but flexible bond angles and dihedral angles) has been used.25,26 (Note a misprint in Table 2 of ref 25: σC should read 0.355 nm.) Nonbonded interactions were truncated at 1.1 nm, and the neighbor list (cutoff of 1.2 nm) was updated every 15 time steps. The reaction-field approximation was used for Coulombic interactions; the dielectric constant was 2.5. Molecular velocity swaps were performed every 150, 300, and 500 steps (W), and sampling of the temperature profile was done every 150, 300, or 500 time steps for the last 500 ps (W ) 150) or 600 ps (W ) 300 and 500) of the simulation. Benzene is a simple, relatively stiff molecule with small atomic partial charges. It can be run without coupling to a thermostat, if the reaction-field approximation is used. In the 900 ps simulation, only a small heating from 303 to ∼310 K has been observed, which is due to truncation of the nonbonded interactions. The observed temperature profiles are linear, as are the concomitant density profiles. An example (W ) 150) is shown in Figure 1. The fact that the temperature profile is linear means that dT/dz is constant across the system. As, by construction, the average heat flux is also constant across the system, it follows that the thermal conductivity is also independent of position (cf. eq 2). Therefore, the perturbations chosen here are small enough for λ not to depend on the distance from the exchange layers. This would be different if the perturbation were large enough to drive the system into nonlinear response. It may be noted in Figure 1 that the lowest temperature that is reached is ∼260 K, which falls below the experimental freezing point of benzene (279 K). The benzene model used here is known to be liquid down to 240 K,27 so there are no complications from a two-phase situation. In addition, we have verified that the diffusion tensor is isotropic. Diffusion in the direction of heat flow would be significantly lower than perpendicular to it, if the cold slab formed a frozen or supercooled barrier. The convergence of the thermal conductiv-

Thermal Conductivities of Molecular Liquids

Figure 1. Temperature (a) and density (b) profiles in the RNEMD simulation of benzene for the strongest perturbation (W ) 150). The two symmetric sides of the simulation cells have been averaged. Linear least-squares fits to the data points are shown, too.

Figure 2. Cumulative average of the thermal conductivity of benzene at the smallest perturbation (W ) 500). It shows the slowest convergence of all simulations in this contribution.

ity to less than 10% is shown in Figure 2 for a W of 500 by a cumulative average, the convergence being even faster for W values of 150 and 300. The thermal conductivities are as follows: λ ) 0.193, 0.197, and 0.208 W m-1 K-1 for W values of 150, 300, and 500, respectively. These values agree well among each other, indicating that all chosen perturbation strengths (W) are in the linear-response domain and eq 1 holds. The fact that the smallest perturbation (W ) 500) shows a slightly larger deviation from the other values is in line with previous RNEMD studies of

J. Phys. Chem. B, Vol. 109, No. 31, 2005 15063 both the thermal diffusion1 and the shear viscosity.2,12 It is attributed to a decreasing statistical accuracy due to a lower signal-to-noise ratio. The calculated thermal conductivity is larger than the experimental one (0.141 W m-1 K-1 at 298 K),28 the disagreement being small for a transport coefficient (the simulated self-diffusion coefficient of this model is a factor of 2 lower than experiment).26,29 Still, there could be other reasons for the disagreement apart from inaccuracies of the model or the experimental measurements. The model is semiflexible, and treats classically degrees of freedom like bond angles involving hydrogens. These are stiff enough to be quantum oscillators in their ground states in reality. As quantum oscillators, they would not be available to transport energy, causing a lower λ. To check for possible size effects on λ, we have repeated the calculations with smaller versions of the benzene system. They had the same extension in directions perpendicular to the heat flux (Lx and Ly), but in the direction of the flux (Lz), they were shortened to 1/3 and 2/3 of the original value, respectively. The thermal conductivity of the smallest system (Lz/3) was 0.270.29 W m-1 K-1, i.e., significantly larger than the value of ∼0.195 W m-1 K-1 found above. The intermediate system (2Lz/ 3), however, exhibited a value (∼0.20 W m-1 K-1) not significantly different from that of the full system. We, therefore, conclude that our system size of approximately 3.17 nm × 3.17 nm × 9.51 nm is large enough to observe a converged thermal conductivity. This should also hold for the other fluids studied here, as system dimensions and molecule sizes are similar (cf. section 5.4). For safety, convergence with system size should be re-evaluated when much larger molecules are treated. 5.2. Cyclohexane. Cyclohexane (750 molecules, F ) 771.1 kg/m3, T ) 303 ( 1.8 K) was simulated at constant NVE conditions using a time step of 1 fs. An automatically optimized all-atom force field with bond constraints has been used,25,30,31 which has no atomic partial charges. Nonbonded interactions were truncated at 1.1 nm, and the neighbor list (cutoff of 1.2 nm) was updated every 15 time steps. Molecular velocity swaps were performed every 150, 300, and 500 steps (W) (simulation times of 500, 900, and 600 ps and sampling times of 300, 600, and 400 ps, respectively), and sampling of the temperature profile was done every 150-500 time steps. Cyclohexane has been chosen because it can be described by a force field without electrostatic interactions; therefore, no heating due to cutoff noise is to be expected, and thermostating can be completely avoided. At the same time, it is slightly more flexible than benzene. Temperature and density profiles are linear for all perturbation strengths (Figure 3); λ is uniform across the system, and the cumulative average of λ shows convergence within a few hundred picoseconds (Figure 4). The calculated thermal conductivities are 0.235 (W ) 150), 0.253 (W ) 300), and 0.238 W m-1 K-1 (W ) 500), the mutual agreement showing that at all perturbation strengths the linear response holds. All values are approximately twice as large as the experimental value (0.123 W m-1 K-1 at 298 K).28 The same considerations as with benzene (section 5.1) apply. The deviation between experimental and calculated thermal conductivity is not too big for a transport property (the diffusion coefficient of this cyclohexane model is 36% below the experimental value),31 but there is the possibility that the classical description predicts a systematically too high λ. The fact that the overestimation is larger for cyclohexane than for benzene could then be due to the larger number of internal quantum degrees of freedom. For cyclohexane, the action of a thermostat (NVT conditions) has also been investigated. We have used a Berendsen thermo-

15064 J. Phys. Chem. B, Vol. 109, No. 31, 2005

Zhang et al.

Figure 5. Calculated thermal conductivities of benzene/cyclohexane mixtures at around 300 K: W ) 150 (pure fluids) and W ) 100 (mixtures).

Figure 3. Temperature (a) and density (b) profiles in the RNEMD simulation of cyclohexane for an intermediate perturbation (W ) 300). The two symmetric sides of the simulation cells have been averaged. Linear least-squares fits to the data points are shown, too.

Figure 4. Cumulative average of the thermal conductivity of cyclohexane at the strongest perturbation (W ) 150).

stat24 with a short coupling time of 0.2 ps, to accentuate its possible effects. We have repeated the W ) 300 calculation (simulation time of 600 ps and sampling time of 400 ps). Although forbidden in principle (see section 2), the thermostat has little effect on λ, which is 0.234 W m-1 K-1. The small difference compared to the NVE calculation above (0.253 W m-1 K-1) is probably more because of the shorter sampling time and the slightly different temperature (300.0 ( 1.7 K) than because of the action of the thermostat. We have also investigated the sensitivity to the force field. We have repeated the W ) 300 calculation (simulation time of

540 ps, sampling time of 360 ps, NVE conditions at T ) 300.5 ( 1.8 K) with nonbonded parameters from the OPLS model of cyclohexane (C ) 0.276 kJ/mol, σC ) 0.35 nm, H ) 0.125 kJ/mol, σH ) 0.25 nm, qC ) -2qH ) -0.12 e, reaction-field dielectric constant of 2), with the same bonded terms as before. This changes the thermal conductivity to 0.182 W m-1 K-1, in the direction of the experimental value, showing that the force field has an influence far greater than that of either the exchange frequency W or the presence or absence of a thermostat. It should be mentioned that both force fields describe the same state of cyclohexane. Both pressures can be considered essentially zero, given the large pressure fluctuations. 5.3. Mixtures of Benzene and Cyclohexane. Using the force fields and settings of sections 5.1 and 5.2, mixtures of benzene and cyclohexane were simulated at 300 K and three compositions: xbenzene ) 0.25 (324 benzenes, 972 cyclohexanes, and F ) 793.7 kg/m3), xbenzene ) 0.5 (645 benzenes, 645 cyclohexanes, and F ) 815.0 kg/m3), and xbenzene ) 0.75 (645 benzenes, 216 cyclohexanes, and F ) 845.5 kg/m3). A very gentle thermostat (coupling time of 50 ps and time step of 2 fs) was used, and W was 100; sampling times were 1.2 ns (xbenzene ) 0.25 and 0.5) and 1.4 ns (xbenzene ) 0.75), and molecular velocity swaps were used. The force field was specifically selected for reproduction of thermodynamic properties of benzene/cyclohexane mixtures,25 but not with thermal conductivities in mind. The thermal conductivity varies in a regular way and exhibits a shallow minimum when xbenzene ) 0.75 (Figure 5). However, the variation of λ with composition is small; between xbenzene values of 0.5 and 1, it is approximately the same as the variations due to different W values, small differences in T, and whether a thermostat is used. Moreover, when compared to the absolute errors for the mixtures,32 the composition influence on λ is minimal. For a start, the experimental λ of benzene is higher than that of cyclohexane, whereas the situation is reversed in the simulations. The only similarity between the two data sets is the negative deviation of λ for mixtures from a linear mixing rule. This points to a clear limitation of the calculation of thermal conductivities by molecular dynamics: The λ values of different systems have to differ significantly to be able to be distinguished by simulation. Our case of a mixture of two liquids with very similar λ values is certainly not a problem usefully studied by MD. 5.4. Water. Water (900 molecules, F ) 1001.6 kg/m3, T ) 301.3 ( 0.06 K) has been simulated using the rigid SPC/E model33 at constant NVT (coupling time of 1.3 ps) conditions for 1200 ps using a time step of 1 fs and constraining all atom-

Thermal Conductivities of Molecular Liquids

Figure 6. Temperature profiles for SPC/E water [303 K, molecular exchange, W ) 150 (9)] and flexible SPC/E water [300 K, molecular exchange, W ) 300 (O)]. The two symmetric sides of the simulation cells have been averaged.

atom distances with SHAKE. Nonbonded interactions were truncated at 0.98 nm, and the neighbor list (cutoff of 1.035 nm) was updated every 15 time steps. The reaction-field dielectric constant was 72. Molecular velocity swaps were performed every 150 steps (W), and the temperature profile was sampled every 150 time steps for the last 530 ps of the simulation. The thermal conductivity calculated for SPC/E water is 0.81 ( 0.01 W m-1 K-1. The same value has been obtained using a completely different nonequilibrium technique.34 It also compares well with previous simulations of SPC/E water9 at 318 K and using Ewald summation, which gave values of 0.730 (W ) 250) and 0.784 W m-1 K-1 (W ) 500), respectively. The rigid SPC/E model overestimates the experimental value of 0.607 W m-1 K-1 at 298 K.28 Details of the long-range interactions do not seem to matter much. A set of simulations using a nonbonded cutoff of 0.75 nm (neighbor list cutoff of 0.80 nm), which for other properties would be considered inadequate, yielded λ values of 0.82 ( 0.01 (W ) 150), 0.779 (W ) 200), and 0.782 W m-1 K-1 (W ) 300). (For the latter two simulations, the conditions were slightly different: F ) 1004.8 kg/m3, T ) 304.2 ( 4.7 K, temperature coupling time of 2 ps, and sampling for 600 ps.) To test the effect of atomic versus molecular velocity exchange, we constructed also a flexible SPC/E water, by replacing ad hoc the three distance constraints with harmonic springs of the same length and force constants of 200 000 kJ mol-1 nm-2. This model used a shorter time step (0.5 fs) and temperature control (0.5 ps). Three simulations with a duration of 1.1-2 ns at 300.4 K using a molecular exchange with a W of 300 (corresponding to a W of 150 of the rigid model because of the shorter time step) yielded a λ of 0.98 ( 0.01 W m-1 K-1. Using atomic exchange instead (four simulations between 1.1 and 2 ns in duration), resulted in a λ of 0.95 ( 0.01 W m-1 K-1. Therefore, the difference between the two exchange algorithms is small. The thermal conductivity also seems to be relatively insensitive to small variations in the harmonic force constants. Changing them to 282 000 kJ mol-1 nm-2 left λ at 0.98 W m-1 K-1. All water models produce linear temperature gradients (Figure 6) for the same perturbation, the one for rigid SPC/E being larger. Thus, λ is again uniform across the system. The flexible SPC/E exhibits a larger thermal conductivity. This is in line with our speculation above that additional degrees of freedom are available for conducting heat. It is interesting to note that the rigid SPC/E water model also has a linear density profile

J. Phys. Chem. B, Vol. 109, No. 31, 2005 15065

Figure 7. Density profiles for rigid SPC/E water [303 K (9)] and flexible SPC/E water [300 K, molecular exchange (O)].

(Figure 7), whereas the density of the flexible SPC/E shows slow convergence and more erratic behavior. This is under further investigation. To study another example for possible size effects, we have repeated the calculation of rigid SPC/E water with identical settings (T ) 301 ( 3 K, 1 ns of sampling, W ) 150), but with a simulation box twice as long in the direction of the heat flux. The thermal conductivity was 0.79 ( 0.01 W m-1 K-1, which is practically unchanged from the value of 0.81 ( 0.01 W m-1 K-1 of the smaller system. As in the case of benzene (section 5.1), the system size of the simulation can, therefore, be considered sufficiently large to achieve a converged value of λ. 5.5. n-Hexane. This molecule was studied as an example of a flexible molecule with soft internal degrees of freedom. Four runs with 300 molecules were performed (F ) 659.7 kg/m3, T ) 302.8 ( 0.1 K) at constant NVE conditions for 1620 ps (after equilibration for 1500 ps) using a time step of 2 fs. A unitedatom modified Ryckaert-Bellemans model35 with nonbonded interactions from the GROMOS96 force field36 has been used; the force field parameters are given in Table 1. Bond lengths were constrained using SHAKE. Nonbonded interactions were truncated at 1.17 nm, and the neighbor list (cutoff of 1.3 nm) was updated every 15 time steps. Molecular velocity swaps were performed every 80 (W) steps, and sampling of the temperature profile was done every 80 time steps. A second simulation 4100 ps in duration (after equilibration for 2400 ps) was performed, where the constraints were replaced by harmonic bonds. The exchange period W was 150, and the temperature was 301.55 ( 0.55 K. The thermal conductivity of n-hexane is 0.107 ( 0.002 W m-1 K-1 for the model with bond constraints and 0.134 ( 0.004 W m-1 K-1 without them. Both values are within 10% of the experimental value of 0.12 W m-1 K-1 at 298 K.28 For both models, the temperature and density profiles are perfectly linear (not shown). As for water, the unconstrained model with its larger number of degrees of freedom has a higher thermal conductivity. In contrast with all other systems discussed so far in this article, the simulated thermal conductivities for n-hexane scatter around the experimental value, rather than being systematically too high. We believe this has to do with the use of united-atom models versus all-atom models. This aspect is discussed further in the Conclusions. The unconstrained model has a slightly lower average pressure (-4.2 MPa) than the constrained model (-0.3 MPa), although both are within the pressure fluctuations (27-50 MPa). To examine the influence of pressure, we have also modified

15066 J. Phys. Chem. B, Vol. 109, No. 31, 2005

Zhang et al.

TABLE 1: United-Atom Force Field for n-Hexanea nonbondedb

mass (amu)

 (kJ/mol)

σ (nm)

The exchange slabs are excluded from the calculation of the temperature gradient. Similar equipartitioning has been found also for other molecular fluids (data not shown).

CH3 CH2

15 14

0.73 0.48

0.393 0.393

6. Conclusions

bond/constraint

length (nm)

force constant (kJ mol-1 nm-2) c

CH2/3-CH2

0.153

350000

bond angle

angle (deg)

force constant (kJ mol-1 rad-2)

CH2/3-CH2-CH2/3

109.47

520.0

dihedral angle CH2/3-CH2-CH2-CH2/3

dihedral angle (deg)

force constant (kJ/mol)

periodicity

180 180 180 180 180

19.68 -13.12 21.23 -6.56 3.94

1 2 3 4 5

a For the functional form of the force field in YASP, see ref 21. First-, second-, and third-neighbor nonbonded interactions within one molecule are excluded. The Lorentz-Berthelot mixing rules are used for interactions between unlike atoms. c Only for harmonic bond terms.

b

Figure 8. Temperature profiles for the different degrees of freedom for n-hexane [model with bond constraints, 59 ps run, molecular exchange every 80 steps (W)]. The two symmetric sides of the simulation cells have been averaged. Equipartition between the degrees of freedom (translation, rotation, and vibration) is evident. Exceptions are the exchange slabs, where the translational temperature deviates from the others, because the exchange of center-of-mass velocities selectively influences this degree of freedom. The exchange slabs are excluded from the calculation of the temperature profile.

the nonbonded parameters (increased the Lennard-Jones σ from 0.393 to 0.395 nm), which brings the pressure of the unconstrained model up to -0.4 MPa, i.e., very close to that of the constrained model. The thermal conductivity changes marginally from 0.134 to 0.139 W m-1 K-1. We therefore conclude that the difference in λ between the constrained and unconstrained n-hexane model is mainly due to the presence of the constraints and not to different pressures. Figure 8 shows the profiles of temperatures calculated from the average kinetic energies of the different degrees of freedom (translation, rotation, and vibration). It shows that the kinetic energy is equipartitioned between them, not only for the system as a whole, but also locally in the individual slabs. Exceptions are the exchange slabs, where the translational temperature deviates from the others, because the exchange of center-ofmass velocities selectively influences this degree of freedom.

The reverse nonequilibrium molecular dynamics method is suitable also for the calculation of thermal conductivities of molecular liquids. The calculated thermal conductivities are relatively insensitive to variations in the RNEMD parameters. In benzene, cyclohexane, and water, λ changed by less than 7% when the exchange period W was varied. This is in line with results on the Lennard-Jones fluid1 and shows that it is easy to keep the simulation in the linear-response regime. For all the cases that have been investigated, the linearity of the temperature profiles showed that the calculated λ is uniform across the system. Similarly, it does not seem to matter much if atomic or molecular exchange is being used for flexible molecules. For water, the results differ by no more than 3%. The simulations also show that, for dense fluids of small to medium size molecules, system sizes of several hundred to some thousand molecules and several nanometers in length are sufficient to avoid size effects. Despite the theoretical misgivings about thermostating simulations aimed at the thermal conductivity, it has turned out that in calculations of practical interest a Berendsen thermostat is harmless. Even with an unrealistically short coupling time of 0.2 ps (in the other simulations, the coupling times are between 1.3 and 50 ps, except one at 0.5 ps), the λ values calculated in the NVT ensemble differed less than 10% from that of an NVE simulation for cyclohexane. Normally, temperature coupling times are at least 1 order of magnitude longer. Thus, any effect of the thermostat should be smaller. This finding, which is in keeping with previous results for the Lennard-Jones fluid,1 has important implications for future applications, as simulations often have to be thermostated for a variety of technical or physical reasons. From the few examples studied here, it appears that the force field potentially has the largest influence on the thermal conductivity. The results for the two cyclohexane nonbonded parameter sets, one with atomic partial charges and the other without, differ by ∼30%. In contrast, changing the harmonic force constant of the flexible water model from one reasonable value to another reasonable value has no effect at all. More examples are needed to determine the force field parameters to which the thermal conductivity is most sensitive. This study has not been undertaken with the purpose of reproducing experimental thermal conductivities or of finding optimized force field parameters. Nonetheless, the results accumulated here allow some speculation about how a computer model of a molecular liquid must be designed, to reproduce or predict its thermal conductivity. It has already been noted (section 5) that the thermal conductivity increases with the available degrees of freedom. All all-atom models (water, benzene, and cyclohexane) produce thermal conductivities larger than experiment. The excess is 40% for benzene, 100 and 50% for the two cyclohexane models, and 30 and 60% for the rigid and flexible water model, respectively. The united-atom nhexane, on the other hand, produces values of -10% (semiflexible with bond constraints) and 10% (fully flexible) around the experiment. Degrees of freedom are removed both by removing explicit atoms (united-atom model) and by using constraints. Removing aliphatic and aromatic hydrogens and their bond and angle vibrations appears to be particularly useful. Their high frequen-

Thermal Conductivities of Molecular Liquids cies make them quantum oscillators at room temperature, which are not available for heat conduction. This view is supported by previous simulations of Hoheisel et al.37-39 They report equilibrium calculations of the thermal conductivities of cyclohexane, benzene, and mixtures thereof using rigid united-atom Lennard-Jones models. They achieve better agreement with experiment than our work, and their deviations are always negative: -15 and -6% for the two cyclohexane models37,38 and -1% for benzene.38 The encouraging implication of this speculation is that, for λ, cheap united-atom models may be the models of choice. On the other hand, there would be not one model that describes the thermal conductivity as well as other properties. A united-atom benzene without an electric quadrupole moment, for example, does not reproduce the experimental liquid structure. The need for constraining further degrees of freedom is not clear from the results. Our unitedatom n-hexane model seems to have struck a good balance about which and how many motions to disallow. While the foregoing discussion of deviations between simulation and experiment is important for the selection of force fields, it must not be forgotten that the agreement is very good. The largest deviation found is a factor of 2 (cyclohexane). Smaller deviations, between 30 and 50%, are much more typical. This makes molecular dynamics a useful predictive tool for thermal conductivities. Acknowledgment. Helpful discussions with Giuseppe Milano and Sylvain Goudeau are gratefully acknowledged. Financial support was provided by the SPP1155 “Molecular Modelling and Simulation in Chemical Engineering” of the Deutsche Forschungsgemeinschaft and by the Fonds der Chemischen Industrie. The John von Neumann Institute for Computing at the Forschungszentrum Ju¨lich made available computer time. References and Notes (1) Mu¨ller-Plathe, F. J. Chem. Phys. 1997, 106, 6082. (2) Mu¨ller-Plathe, F. Phys. ReV. E 1999, 59, 4894. (3) Mu¨ller-Plathe, F.; Reith, D. Comput. Theor. Polym. Sci. 1999, 9, 203. (4) Reith, D.; Mu¨ller-Plathe, F. J. Chem. Phys. 2000, 112, 2436. (5) Bordat, P.; Reith, D.; Mu¨ller-Plathe, F. J. Chem. Phys. 2001, 115, 8978. (6) Donko´, Z.; Hartmann, P. Phys. ReV. E 2004, 69, 016405. (7) Choi, S.-H.; Maruyama, S. J. Korean Phys. Soc. 2004, 45, 897. (8) Padgett, C. W.; Brenner, D. W. Nano Lett. 2004, 4, 1051. (9) Bedrov, D.; Smith, G. D. J. Chem. Phys. 2000, 113, 8080. (10) Bedrov, D.; Smith, G. D.; Sewell, T. D. Chem. Phys. Lett. 2000, 324, 64.

J. Phys. Chem. B, Vol. 109, No. 31, 2005 15067 (11) Nieto-Draghi, C.; Bonet A Ä valos, J.; Rousseau, B. J. Chem. Phys. 2003, 119, 4782. (12) Bordat, P.; Mu¨ller-Plathe, F. J. Chem. Phys. 2002, 116, 3362. (13) Soddemann, T.; Du¨nweg, B.; Kremer, K. Eur. Phys. J. E 2001, 6, 409. (14) Guo, H. X.; Kremer, K.; Soddemann, T. Phys. ReV. E 2002, 66, 061503. Loison, C.; Mareschal, M.; Kremer, K.; Schmid, F. J. Chem. Phys. 2004, 119, 13138. Guo, H. X.; Kremer, K. J. Chem. Phys. 2003, 118, 7714. Guo, H. X.; Kremer, K. J. Chem. Phys. 2003, 119, 9308. (15) Soddemann, T.; Auernhammer, G. K.; Guo, H. X.; Du¨nweg, B.; Kremer, K. Eur. Phys. J. E 2004, 13, 141. (16) Soddemann, T.; Du¨nweg, B.; Kremer, K. Phys. ReV. E 2003, 68, 046702. (17) Mu¨ller-Plathe, F.; Bordat, P. in NoVel Methods in Soft Matter Simulations; Karttunen, M., Vattulainen, I., Lukkarinen, A., Eds.; Lecture Notes in Physics Vol. 640; Springer: Heidelberg, Germany, 2004; p 310. (18) Normand, C.; Pomeau, Y.; Velarde, M. G. ReV. Mod. Phys. 1993, 49, 581. (19) de Groot, S. R.; Mazur, P. Non-Equilibrium Thermodynamics, Dover: New York, 1984. (20) Nieto-Draghi, C.; Bonet A Ä valos, J. Mol. Phys. 2003, 101, 2303. (21) Mu¨ller-Plathe, F. Comput. Phys. Commun. 1993, 78, 77. (22) Mu¨ller-Plathe, F.; Brown, D. Comput. Phys. Commun. 1991, 64, 7. (23) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids, Oxford University Press: Oxford, U.K., 1987. (24) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Di Nola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (25) Milano, G.; Mu¨ller-Plathe, F. J. Phys. Chem. B 2004, 108, 7415. (26) Mu¨ller-Plathe, F. Macromolecules 1996, 29, 4782. (27) Witt, R.; Sturz, L.; Do¨lle, A.; Mu¨ller-Plathe, F. J. Phys. Chem. A 2000, 104, 5716. (28) Lide, D. R., Ed. CRC Handbook of Chemistry and Physics, 82nd ed., CRC Press: Boca Raton, FL, 2001. (29) Mu¨ller-Plathe, F. Chem. Phys. Lett. 1996, 252, 419. (30) Faller, R.; Schmitz, H.; Biermann, O.; Mu¨ller-Plathe, F. J. Comput. Chem. 1999, 20, 1009. (31) Schmitz, H.; Faller, R.; Mu¨ller-Plathe, F. J. Phys. Chem. B 1999, 103, 9731. (32) Landolt-Bo¨rnstein, Zahlenwerte und Funktionen; Springer: Berlin, 1969; Vol. II, Chapter 5. (33) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (34) Terao, T.; Mu¨ller-Plathe, F. J. Chem. Phys. 2005, 122, 081103. (35) van der Spoel, D.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Thieleman, D. P.; Sijbers, A. L. T. M.; Hess, B.; Feenstra, K. A.; Lindahl, E.; van Drunen, R.; Berendsen, H. J. C. Gromacs User Manual Version 3.0 (http://diamond.kist.re.kr/SMS/md-code/GROMACS/manual-ah-3.0.pdf). (36) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, I. G. Biomolecular Simulations: The GROMOS96 Manual and User Guide; VDF: Zu¨rich, 1996. (37) Hoheisel, C.; Wu¨rflinger, A. J. Chem. Phys. 1989, 91, 473. (38) Luo, H.; Hoheisel, C. J. Chem. Phys. 1992, 96, 3173. (39) Schaink, H. M.; Luo, H.; Hoheisel, C. J. Chem. Phys. 1993, 99, 9912.