Nonequilibrium Molecular Dynamics Calculation of the Thermal

11516

J. Phys. Chem. B 2007, 111, 11516-11523

Nonequilibrium Molecular Dynamics Calculation of the Thermal Conductivity of Amorphous Polyamide-6,6 Enrico Lussetti,† Takamichi Terao,‡ and Florian Mu1 ller-Plathe*,† Eduard-Zintl-Institut fu¨r Anorganische und Physikalische Chemie, Technische UniVersita¨t Darmstadt, Petersenstrasse 20, D-64287 Darmstadt, Germany, and Department of Mathematical and Design Engineering, Gifu UniVersity, Yanagido 1-1, Gifu 501-1193, Japan ReceiVed: May 16, 2007; In Final Form: July 17, 2007

The thermal conductivity of the amorphous phase of polyamide-6,6 is investigated by nonequilibrium molecular dynamics simulations. Two different algorithms are used, reverse nonequilibrium molecular dynamics and the dual-thermostat method. Particular attention is paid to the force field used. Four different models are tested, flexible and rigid bonds and all-atom and united-atom descriptions. They mainly differ in the number of high-frequency degrees of freedom retained. The calculated thermal conductivity depends systematically on the number of degrees of freedom of the model. This dependence is traced to the quantum nature of the fast molecular vibrations, which are incorrectly described by classical mechanics. The best agreement with experiment is achieved for a united-atom model with all bonds kept rigid. It could be shown that both the thermal conductivity and the heat capacity of a model show a similar (but not equal) dependence on its number of degrees of freedom. Hence, the computationally more convenient heat capacity can be used for force field optimization. A side result is the anisotropy of the thermal conductivity in stretched polyamide; heat conduction is faster parallel to the drawing direction than perpendicular to it.

1. Introduction The thermal conductivity of a polymer is a central quantity for engineering applications. It is of importance for the finished polymer products and devices. In thermal insulation applications, the thermal conductivity of the material should be as low as possible, while in the encapsulation of electronic devices, a high thermal conductivity is desirable. The role of the thermal conductivity is, however, not limited to the final product. In the manufacturing and processing of polymers, thermal conductivities may be design parameters. For these cases, they may be hard to obtain experimentally because the process may be carried out at temperatures and pressures far from ambient conditions, they may take place in closed vessels, or they may involve mixtures, in which chemical reactions are still going on. In such exotic conditions, simulations often offer a convenient way to estimate missing materials’ properties. In contrast to other transport properties, the calculation of thermal conductivities by molecular dynamics simulation has, for a long time, been impaired by the lack of a reliable and robust simulation method.1 This has been mainly due to the necessity to define and calculate a heat flux vector on the molecular scale for both equilibrium (Green-Kubo) and traditional nonequilibrium molecular dynamics methods. With the arrival of nonequilibrium methods, which bypass this need, such as the heat-exchange method,2 the reverse nonequilibrium molecular dynamics (RNEMD) method,3 the dual-thermostat (DT). and heat-injection methods,4 the calculation of thermal conductivities has become much more practical. In addition to model fluids, molecular fluids and mixtures, modeled by realistic force fields, have been treated.5,6 † ‡

Technische Universita¨t Darmstadt. Gifu University.

In this paper, we extend the nonequilibrium methods to the calculation of the thermal conductivity of polymers. As a first example, we use polyamide-6,6, which is of technical importance in high-temperature applications and, therefore, well characterized. We focus on the amorphous phase of this polymer, although the commercial material usually is semicrystalline. (The study of purely crystalline and purely amorphous polymer phases will be the subject of forthcoming papers.) With several robust numerical schemes for nonequilibrium molecular dynamics in hand, we turn to the requirements on the molecular dynamics model or force field for a faithful reproduction or prediction of thermal conductivities of polymer materials. There have been strong hints from our simulations of molecular fluids6 that the thermal conductivity is predicted too high by conventional force fields and that this is due to the classical treatment of fast quantum degrees of freedom in molecular dynamics. This hypothesis is to be checked systematically with four different models of polyamide-6,6, with the aim of developing a protocol for practical molecular dynamics calculations of the thermal conductivity with existing force fields. A second objective of the present contribution is a first study of the anisotropy of the thermal conductivity in oriented polymer samples. As stiff degrees are known to transport heat more efficiently (by phonons) than soft degrees (by collisions), stretching an amorphous polymer should increase the thermal conductivity parallel to the stretching direction at the expense of the perpendicular directions. The effect of stiff versus soft interactions has recently been demonstrated by both experiment and simulation on a crystalline but disordered inorganic system, namely, WSe2.7 2. Nonequilibrium Methods for Polymers 2.1. Reverse Nonequilibrium Molecular Dynamics. RNEMD has been discussed in detail several times; see, for example, ref

10.1021/jp0737956 CCC: $37.00 © 2007 American Chemical Society Published on Web 09/07/2007

Thermal Conductivity of Amorphous Polyamide-6,6

J. Phys. Chem. B, Vol. 111, No. 39, 2007 11517

8. The method was originally designed for atomic liquids3 and later extended to molecular fluids and mixtures.6 In this section, we only recall what is necessary to understand the modifications for polymers proposed in this contribution. The idea of the method is to drive a heat flux through the simulated system as the primary perturbation. Energy is continuously transferred from one region of the system, denoted as the “cold” region, to another one, the “hot” region. Let us assume the two regions to be separated in the z direction. The transfer is done in the artificial, nonphysical manner discussed below. When the steady state has been reached, an equal amount of energy per time and area jz flows in the opposite direction through the system by a physical mechanism (thermal conduction) since energy is conserved. As a result, a temperature gradient 〈dT/dz〉 develops in the system. The ways of calculating the temperature profile and its gradient have been described for fluid models with and without constraints.8 They are used unchanged for polymers in this work. From the known imposed jz and the calculated 〈dT/dz〉, the thermal conductivity λ is calculated as λ ) -jz/〈dT/dz〉, assuming linear response (Fourier’s law) to hold. The second ingredient needs a modification to be applicable to polymers. The nonphysical energy transfer is implemented as a Maxwell’s daemon, which, at certain intervals, selects two atoms for their velocities, the fastest (hottest) atom in the cold region and the slowest (coldest) atom in the hot region.3 By exchanging their velocity vectors, energy is transferred from cold to hot. The average energy flux in the steady state is evaluated as

jz )

1

m

∑ (V2hot - V2cold) 2tA transfers 2

(1)

where t is the duration of the simulation, A is the cross-sectional area of the simulation box perpendicular to the flow direction z, Vhot and Vcold denote the velocities of the hot and the cold atom of equal mass m, and the factor 2 arises from the periodicity. This so-called atomic exchange algorithm works well for atomic fluids and for molecular systems without constraints. Exchanging atom velocities, however, is not possible if the atom is fixed by constraints, for example, by rigid bond constraints. Changing the velocity of one atom in a constraint is not compatible with the constraint conditions, which dictate also the velocities of the atoms involved, and it would lead to nonconservation of energy and ultimately to the failure of the constraints’ solver. On the other hand, one often would like to use force fields with constraints for physical as well as technical reasons. For systems of small molecules with constraints, there is a simple solution; one selects entire molecules based on their center-of-mass velocities and exchanges the latter.6 As centerof-mass velocities do not interfere with molecular conformations or the velocities of the member atoms with respect to each other, exchanging them does not affect constraints either. The development of the analysis proceeds analogously to the atomic case. This exchange is called molecular. Exchanging the center-of-mass velocities is not an option if the molecules are polymers. First, a single molecule may extend through large parts of the system; therefore, it is difficult to direct the energy transfer to the hot and cold regions. Second, there may be very few macromolecules in a simulation, and it may be difficult to find a suitable molecule (right region, right velocity) for exchange. Third, the velocity swap would perturb the Newtonian dynamics also outside the exchange slabs, where

later analysis takes place. If there are no constraints in the model, the atomic exchange can be used without modifications also for polymers. If the model involves constraints, we resort to a partial treatment, which we call semimolecular exchange. The polymer is broken down into groups of neighboring atoms. Within a group, there may be constraints. The connection between different, including adjacent, groups must not involve constraints but must be fully flexible. One could for example constrain all C-H bonds in polyamide. This would lead to CH2 moieties being groups in the sense of our definition. Then, the C-C bonds between neighboring methylene groups would have to be treated by a flexible potential, such as harmonic. In polystyrene, a phenyl ring may be completely rigid. However, the connection between styrene units must contain one flexible backbone C-C bond. With the polymer partitioned into groups like this, one can perform selection and energy transfer based on the center-of-mass velocities of suitable groups. As before, if the swap of center-of-mass velocity vectors is performed between groups of equal mass, the total energy and total linear momentum are conserved. This remains true also in the presence of constraints internal to the semimolecular groups. 2.2. Dual-Thermostat Method. Like in RNEMD, in the DT method, there is a hot region and a cold region. The regions are maintained at their respective temperatures by two locally acting Berendsen thermostats.9 Note that the intervening space, where the analysis is performed, is not thermostatted at all, but free Newtonian dynamics takes place. As in RNEMD, the calculation of the heat flux is straightforward because the amount of energy added or removed at each time step by the Berendsen thermostat can be easily calculated, and it is not necessary to define a heat flux on the molecular scale. The DT method, which is fully described in ref 4, can be applied if the system is in the linear-response regime. The DT method has been found suitable for any force field with or without constraints, although, strictly speaking, the velocity rescaling action of the Berendsen thermostat violates the constraint conditions on the atomic velocities. In practice, however, the velocities are scaled by a factor very close to one; therefore, the perturbation is numerically unimportant. 3. Computational Details 3.1. Simulation Details. All molecular dynamics simulations were carried out with the YASP package,10 which uses the leapfrog algorithm11 and orthorhombic periodic boundary conditions. The simulation cells were elongated in the z direction, in which the energy flow was imposed (Lx ) Ly ) Lz/2). Bond constraints, if present, were solved by the SHAKE method12 to a relative tolerance of 10-6. Temperatures and densities were chosen at or close to ambient conditions, where experimental data are most abundant. When the temperature was kept constant, this was done by the Berendsen method,9 with the target temperature being 300 K, unless specified otherwise. The intramolecular force field contains harmonic bond stretching or constraints, harmonic bond angle bending, and periodic cosine-type torsional potentials. 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, with a dielectric constant of 5. For details of the functional form, see ref 10. Nonbonded interactions were evaluated from a Verlet neighbor list, which was updated every 30 time steps using a link-cell method. Unless specified otherwise, all simulations were carried out at constant volume and temperature. Normally, the time step

11518 J. Phys. Chem. B, Vol. 111, No. 39, 2007

Lussetti et al.

TABLE 1: Thermal Conductivity and Heat Capacity for the Different Models of Amorphous Polyamide-6,6 at 300 K model fl-AA fl-UA pc-UA fc-UA experiment

flexible bonds per repeat unit

degrees of freedom per repeat unit

λ (W m-1 K-1)

CV (J mol-1 K-1)a

38 18 5 0

114 54 41 36

0.45 ( 0.04 0.32 ( 0.02 0.29 ( 0.02 0.27 ( 0.03 0.15-0.30b

992 ( 45 469 ( 35 358 ( 20 310 ( 20 319c; 310d; 362e

a Heat capacity per 1 mole of repeat units. b Our reference value is 0.24 W m-1 K-1;14 the span reflects typical discrepancies between experimental values for different polymer samples (see, e.g., refs 15 and 16). The experimental values typically are extrapolations of thermal conductivities of semicrystalline samples to zero crystallinity. c Ref 17. d Ref 18. e Ref 19.

was 1 fs. The cutoff radius, unless specified otherwise, was 0.9 nm, and the neighbor list cutoff was 1.0 nm. For the Berendsen thermostat, coupling times of 0.8 and 8 ps were used (unless specified otherwise) in order to verify the independence of the result from this difference. Both coupling times were sufficient to keep the measured average temperature within 1 K of the target temperature. Specific settings of simulation parameters are reported below for the individual simulations. A nonequilibrium run typically lasted 1000 ps, of which the last 800 ps were used for the production. 3.2. Sample Preparation. The sample contained 48 chains of molecular weight 4539 g/mol, which corresponds to 20 repeat units and 40 amide groups. The initial atomic coordinates and velocities were generated as described in ref 13. First conformations were obtained via a pivot Monte Carlo calculation of single polymer chains in vacuum with a suitable bond-based interaction cutoff to generate melt-like structures. The chains were inserted into the periodic simulation cell and then equilibrated for 4 ns by equilibrium MD. For further equilibration, the system underwent an annealing cycle, consisting of a quick warm-up to 500 K, 1 ns of equilibration at 500 K, and a cooling down to 300 K over 300 ps, followed by a further 300 ps of equilibration at 300 K. The experimental density of amorphous polyamide6,6 is 1.07-1.1 g/cm3; in the simulations, the values 1.07 and 1.11 were used in order to compare the differences in the results and to ascertain the independence of the conclusions from this difference. 3.3. Force Fields. In order to assess the influence of different force fields on the thermal conductivity, with special attention to the number of degrees of freedom per chain, four different force fields were used. Flexible All-Atom Model (denoted as fl-AA), Without Any Constraints. This is a flexible version of the fully bondconstrained all-atom model described in ref 13, where all of the constraints have been replaced ad hoc by harmonic terms, all with a force constant of 200000 kJ mol-1 nm-2. Flexible United-Atom Model (fl-UA), Without Constraints, DeriVed from the PreVious Model. All CH2 and CH3 groups are treated as united atoms, which incorporate the aliphatic hydrogens. All intramolecular bonds are modeled by harmonic terms, all with a force constant of 335000 kJ mol-1 nm-2. Partially Constrained United-Atom Model (pc-UA), DeriVed from the PreVious (fl-UA) Model. The force field parameters are the same; all bonds are constrained except (a) the bonds between amide groups and adjacent CH2 groups and (b) the central CH2-CH2 bonds of the hexamethylene units. In this way, the chain is partitioned into internally bond-constrained groups (amide, (CH2)3, (CH2)4), which are connected to each other by harmonic bonds with a force constant of 335000 kJ mol-1 nm-2. Fully Bond-Constrained United-Atom Model (fc-UA), DeriVed from the fl-UA Model. All bonds are constrained. From the first to the last model, the number of degrees of freedom per repeat unit decreases (Table 1). As the mass

densities are the same for all models, the number of degrees of freedom per volume also decreases. The force field parameters are listed in detail in Table 2. All models are compatible with at least one of the nonequilibrium methods described in section 2. The RNEMD method was used for fl-AA, fl-UA, and pcUA, and the DT method was used for fc-UA. It has been checked elsewhere4 that RNEMD and DT give the same thermal conductivity for a given system. 3.4. Simulations: Isotropic Sample. For the fl-AA model, RNEMD with atomic velocity exchange was used. Velocity swaps were performed every 110 time steps and in one case every 75 time steps in order to assess the possible effect of this difference. The density was normally 1.11 g/cm3; in one case, a density of 1.07 g/cm3 was tested. Coupling times of 0.8 and 8 ps for the Berendsen thermostat were used in order to evaluate the possible effect of this difference. Two different choices for the spring constants were made (section 4.1), adjusting the time step accordingly, in order to make the comparison with the other simulations more transparent. For the fl-UA model, the atomic velocity exchange was used, and velocity swaps were performed every 75 time steps; further tests with swaps every 55 and 130 time steps were done. The density was 1.07 g/cm3. The coupling time for the Berendsen thermostat was 8 ps. A bigger cutoff of 1.17 nm (and 1.30 nm for the neighbor list cutoff) were tested in order to verify that the results do dot depend on the choice of the cutoff distance. For the pc-UA model, the semimolecular velocity exchange was used, and velocity swaps were performed every 75 time steps; further tests with swaps every 55 and 165 time steps were done. Results for densities of 1.07 g/cm3 and 1.11 g/cm3 were compared. Coupling times of 0.8 and 8 ps for the Berendsen thermostat were used in order to evaluate the possible effect of this difference. For the fc-UA model, the DT method was used. The thermostatted slabs were kept at 290 and 310 K using a coupling time of 0.1 ps. The intervening region was not thermostatted. The density was 1.07 g/cm3. 3.5. Simulations: Stretched Sample. An amorphous polyamide-6,6 sample was stretched in one direction, and the thermal conductivity parallel and perpendicular to the stretching direction was calculated. The fl-UA model was used. Where NPT conditions were needed, the Berendsen method was used also for controlling the pressure. The preparation of a stretched sample of 24 polymer chains initially at 1.07 g/cm3 and 300 K took the following steps: (1) Heat to 500 K over 400 ps at constant volume. The coupling time of the thermostat was 0.1 ps, and the bath temperature was increased at a rate of 0.5 K/ps. (2) Equilibrate at 500 K for 400 ps (with a coupling time of 0.1 ps) at constant volume. (3) Stretch the sample at 500 K in a NPT simulation. The temperature coupling time was 1 ps. The pressure coupling was anisotropic; the initial target pressure was 200000 kPa (the average pressure of step 2), and the coupling times were 150,

Bond Stretching: Vbonds ) (kb/2)|r - r0|2 or Bond Constraints fl-AA

pc-UAb

fl-UA

fc-UA

constraint length/nm

type

ro/nm

kb/kJ mo1-1 nm-2

type

ro

kb

type

Ch-H Ch-Ch Ch-Ca Ch-N Ca-N N-Ha Ca-O

0.109 0.1522 0.152 0.1449 0.1335 0.101 0.1229

200000 200000 200000 200000 200000 200000 200000

CH2-CH2/3 CH2/3-Ca CH2-N Ca-N N-Ha Ca-O

0.1522 0.152 0.1449 0.1335 0.101 0.1229

350000 350000 350000 350000 350000 350000


ro

kb

type

constraint length

0.1522 0.152 0.1449

350000 350000 350000


0.152 0.152 0.1449 0.1335 0.101 0.1229

0.1335 0.101 0.1229

Angle Bending: Vangles ) (kφ/2)(φ - φ0)2 fl-AA type

φo/deg

Ch-Ch-Ch/a H-Ch-H Ch-Ch-H Ca-Ch-H Ch-Ca-O Ch-Ca-N Ca-N-Ha Ch-N-Ha

111° 109.5° 109.5° 109.5° 120.4° 116.6° 120° 120°

fl-UA kφ/kJ mo1

-1

rad

-2

250 145 150 210 335 293 125 209

pc-UA

type

fc-UA

kφ

type

kφ

type

CH2/3-CH2-CH2/Ca

111°

250

CH2/3-CH2-CH2/Ca

111°

250

CH2/3-CH2-CH2/Ca

111°

250

CH2/3-Ca-O CH2/3-Ca-N Ca-N-Ha CH2-N-Ha

120.4° 116.6° 120° 120°

335 293 125 209


120.4° 116.6° 120° 120°

335 293 125 209


120.4° 116.6° 120° 120°

335 293 125 209

φo

φo

φo

kφ


TABLE 2: Force Field Parametersa

Torsions: Vtors ) ∑pktors p [1 + cos p(τ - τ0)] (+ symbols below separate terms of the sum with different periodicity p) fl-AA type

p

fl-UA

-1 τo/deg ktors p /kJ mol

2.0 + 3.6 8.1 + 3.6 25 + 40 9 35

fc-UA

p

τo

ktors p

type

p

τo

ktors p

type

p

τo

ktors p

CH2-CH2-CH2-CH2/3 CH2-CH2-CH2-Ca CH2-N-Ca-CH2/3 Ca-N-CH2-CH2 O-Ca-N-Ha

1+3 1+3 1+2 1 1

180° 180° 180° 180° 180°

2.0 + 3.6 8.1 + 3.6 25 + 40 9 35


1+3 1+3 1+2 1 1

180° 180° 180° 180° 180°

2.0 + 3.6 8.1 + 3.6 25 + 40 9 35


1+3 1+3 1+2 1 1

180° 180° 180° 180° 180°

2.0 + 3.6 8.1 + 3.6 25 + 40 9 35

Nonbonded Interactions c: Vnb ) (qiqi/4π0)[1/rij + ((rf - 1)/(2rf + 1)) r2ij/r3cutoff] + [4[(σ/rij)12 - (σ/rij)]] fl-AA

fl-UA

pc-UA

fc-UA

type

/kJ mol-1

σ/nm

q/e

type

σ

q

type

σ

q

type

σ

q

H Ha Ch Ca N O

0.0657 0.0657 0.458 0.36 0.7118 0.8792

0.24 0.1 0.35 0.35 0.33 0.32

0.08229 0.2819 -0.1506 0.5973 -0.4057 -0.5479

Ha CH2/3 Ca N O

0.0657 0.48/0.73 0.36 0.7118 0.8792

0.1 0.393 0.35 0.33 0.32

0.292 0 0.625 -0.39 -0.527

Ha CH2/3 Ca N O

0.0657 0.48/0.73 0.36 0.7118 0.8792

0.1 0.393 0.35 0.33 0.32

0.292 0 0.625 -0.39 -0.527

Ha CH2/3 Ca N O

0.0657 0.48/0.73 0.36 0.7118 0.8792

0.1 0.393 0.35 0.33 0.32

0.292 0 0.625 -0.39 -0.527

a C stands for a carbon involved in CH or CH groups, and C is for an amide carbon. H is a CH or CH hydrogen, and H is an amide hydrogen. In the united-atom models, CH and CH groups are h 2 3 a 2 3 a 2 3 treated as single atoms. The chains are terminated at each end by a CH3 group; at one end, the CH3 group replaces the CH2 group belonging to a tetramethylene group; at the other end, the CH3 group is bonded to an amide carbon. b In the pc-UA model, all bonds are constrained except (a) the bonds between amide groups and adjacent CH2 groups and (b) the central CH2-CH2 bonds of the hexamethylene units. c First- and second-neighbor nonbonded interactions within one chain are excluded but not third-neighbor interactions. The cutoff is 0.9 nm; the reaction field dielectric constant rf for electrostatic interactions is 5. The Lorentz-Berthelot combination rule is used to determine the σ and parameters for unlike pairs: ij ) (iijj)0.5 and σij ) (σii + σjj)/2.


Ch-Ch-Ch-Ch 1+3 180° Ch-Ch-Ch-Ca 1+3 180° Ch-N-Ca-Ch 1+2 180° Ca-N-Ch-Ch 1 180° O-Ca-N-Ha 1 180°

pc-UA

type

11520 J. Phys. Chem. B, Vol. 111, No. 39, 2007 150, and 1 ps for the x, y, and z directions, respectively. The target pressure was kept constant for the x and y directions; for the z direction, it was decreased by -800 kPa/ps. As a consequence, the z dimension of the box increased from 5.527 to 7.0 nm, and the x and y dimensions shrank from 5.527 to 5.1 nm. The total density was 0.99 g/cm3 at the end of this step, which is lower than the density of the isotropic sample. (4) Cool to 300 K for 400 ps in a NPT simulation. The initial bath temperature of 500 K was decreased at a rate of -0.5 K/ps, and the target pressure was 50000 kPa in all directions. The coupling times were 1 ps for the temperature and 35 ps for the pressure. At the end of this step, the density was 1.08 g/cm3. (5) Isotropic compression of the sample in a NPT simulation at 300 K for 280 ps, with a target pressure of 250000 kPa and coupling times like those in step 4. At the end of this step, the density was 1.12 g/cm3. (6) Stretch in a NPT simulation at 300 K for 360 ps. The purpose of this stage was to adjust the final density to the same value as that of the isotropic sample. The coupling times were 1 ps for the temperature and 107, 107, and 1 ps for the pressure in the x, y, z directions, respectively. The initial target pressure of 50000 kPa was decreased at a rate of -3000 kPa/ps. The z dimension of the box increased from 6.7 to 7.0 nm; the other dimensions remained unchanged at 4.9 nm. Summarizing, most of the stretching, which totalled ∼27%, took place at 500 K in step 3, while only minor adjustments of the box dimensions resulted from the other NPT simulations. As a measure of the anisotropy of the sample, we considered the gyration tensor R and calculated the quantity

2〈Rzz〉/〈Rxx + Ryy〉 as an order parameter. Its value was 1.5 ( 0.5 for the unstretched sample and 3.0 ( 0.9 for the stretched sample. For calculation of the thermal conductivity in the z direction (in the stretch direction), the prepared sample was duplicated in the z direction so that the elongated side of the box was 14 nm and the short sides were 4.9 nm. For calculation of the thermal conductivity in the x direction (perpendicular to the stretch direction), the initial sample was triplicated in the x direction so that the elongated side of the box was 14.7 nm and the short sides were 4.9 and 7 nm. 4. Results and Discussion 4.1. Influence of Simulation Parameters. The influence of a number of technical simulation parameters on the thermal conductivity must be clarified in order to establish confidence in the accuracy of the results; for details, see, for example, refs 3, 6, and 8. For both nonequilibrium methods, the system must be in the linear response regime, that is, the heat flux jz and the temperature gradient 〈dT/dz〉 must be proportional. This has been checked for all systems. As an example, we show the slope of the temperature profile as a function of the imposed heat flux for the pc-UA model (Figure 1), which was simulated by the RNEMD method. For the fl-UA and fc-AA models, the dependence was equally linear. To quote a number, for the flAA force field, the thermal conductivity changed from 0.452 ( 0.035 W m-1 K-1 for a velocity swap interval of W ) 110 to 0.446 ( 0.035 W m-1 K-1 for W ) 75, both values agreeing well to within their error bars. In the RNEMD simulations, the temperatures of the hottest and coldest slabs typically were ∼330 and ∼270 K, respectively. In the DT calculation on the fc-UA model, the temperature spread was much smaller (between 310 and 290 K). Therefore, linear response holds here even more. The temperature coupling time of the Berendsen thermostat had no influence on the calculated thermal conductivity.

Lussetti et al.

Figure 1. Linear response of the reverse nonequilibrium method for the partially constrained united-atom (pc-UA) model of polyamide6,6. The measured temperature gradient is a linear function of the imposed heat flux. The three perturbation strengths (as a number of velocity swaps per picosecond or (W∆t)-1) are indicated in the figure. They correspond to velocity swaps every W ) 165, 75, and 55 time steps of length ∆t ) 1 ps.

Coupling times of 0.8 and 8 ps were tested for both the fl-AA model and the pc-UA model. For the fl-AA force field, the thermal conductivities are 0.452 ( 0.035 and 0.445 ( 0.035 W m-1 K-1 for coupling times of 0.8 and 8 ps, respectively. For the pc-UA model, the corresponding thermal conductivities are 0.287 ( 0.015 and 0.290 ( 0.015 W m-1 K-1. For both models, the coupling time does not influence thermal conductivity. The values agree to within their error bars. The force constants for the harmonic bonds of the fl-AA model were increased from 200000 to 300000 kJ mol-1 nm2. At the same time, the time step was reduced from 1 to 0.4 fs, with the velocity exchange frequency being held constant at 9.1 ps-1. The value of the thermal conductivity changed from 0.45 ( 0.04 to 0.46 ( 0.04 W m-1 K-1, again by less than the error bar. For the other models, all of which have fewer flexible bonds, the differences were even smaller, as expected. Similarly, the cutoff distance for the nonbonded interactions was increased from 0.9 to 1.17 nm for the fl-UA model. Again, the resulting thermal conductivity, 0.33 ( 0.02 W m-1 K-1, was the same to within the error bars as that for the shorter cutoff (0.32 ( 0.02 W m-1 K-1). Thus, the influence of force field details for the degrees of freedom retained is much less than the actual number of degrees of freedom (section 4.2). 4.2. Influence of the Number of Degrees of Freedom of the Model. The most important result of this work is the dependence of the calculated thermal conductivity of polyamide6,6 on the force field. The four different force fields span a range of thermal conductivities from 0.26 to 0.45 W m-1 K-1 (Table 1). The dependence on the number of degrees of freedom of the model is striking (Figure 2). The model with the most degrees of freedom is the fl-AA model, which has no constraints and all atoms present, including aliphatic hydrogens. It has the highest thermal conductivity. At the other extreme is the fcUA model, with all bonds constrained and no aliphatic hydrogens. The other two models fall in between. The experimental thermal conductivity of amorphous polyamide-6,6 is estimated to be around 0.24 W m-1 K-1; see the footnote of Table 1. This is close to the calculated value for the most constrained model in our series, the fc-UA model (0.27 W m-1 K-1). The following hypothesis attempts to explain the increasing deviation from experiment with the degrees of freedom: in reality, the stiff degrees of freedom are quantum oscillators in their vibrational ground state. Thermally, they


Figure 2. The calculated thermal conductivity of polyamide-6,6 of a force field depends on the number of degrees of freedom of the force field. The line is a linear fit to the data.

cannot be excited, as their vibrational energy (in wavenumbers) is 1000-2500 cm-1, whereas kT at room temperature is 207 cm-1. Hence, they are not available for the transport of energy. In contrast, in classical molecular dynamics simulations, all degrees of freedom can store and transport energy, as they are all populated according to the equipartition theorem. Eliminating the hydrogen atoms with their fast stretching and bending vibrations and freezing all remaining bond vibrations leaves just those soft degrees of freedom, which do contribute to heat transport at room temperature. While this concept explains the correlation between observed thermal conductivity and the force field ansatz, it is certainly oversimplified. The relation appears to be linear in the limited range presented in Figure 2. If the assumption of every degree of freedom contributing equally to the thermal conductivity were true, however, the line would have to pass through the origin, which it evidently does not. The reason is that vibrational modes contribute differently to the thermal conductivity; high-frequency modes such as C-H bond vibrations are often localized and do not contribute as much to the thermal conductivity as delocalized low-frequency vibrations, even if they are treated classically.20,21 Therefore, the erroneous inclusion of a highfrequency mode into the classical treatment is not as damaging as the erroneous inclusion of a low-frequency mode. It is known that in the calculation of heat capacities CV, one faces similar problems and simulated heat capacities are much closer to experiment if the models are highly constrained. As an example, a united-atom fully rigid model of dimethyl sulfoxide, that is, all bonds and all angles kept fixed, gave a CV of 101 J mol-1 K-1, compared to an experimental value of 118 J mol-1 K-1.22 We have therefore calculated the heat capacities at constant volume for all of our polyamide-6,6 models from a line fitted to the total energy at three different temperatures, 290, 300, and 310 K. They are listed in Table 1, together with several estimates of the heat capacity at constant pressure Cp derived from experiment. (The difference between Cp and CV can be calculated using the thermodynamic relation Cp - CV ) R2TV/KT, where the thermal expansion coefficient R is 2.43 K-1 for polyamide-6,6 and the isothermal compressibility KT is (3300 MPa-1).23 A density of 1.07 g/cm3 was used. The difference between Cp and CV was 12 J/K per mole of repeat units and, therefore, insignificant.) The data have been collected in Figure 3. It shows that there is indeed a correlation between the disagreements of the


Figure 3. Comparison of the deviation of calculated and experimental thermal conductivities λ and heat capacities CV. The data for the various polyamide-6,6 models are from this work. The data for low molecular fluids from ref 6 are included for comparison. The line is a linear fit to all data points.

Figure 4. Thermal conductivity of the pc-UA model of amorphous polyamide-6,6 at 300 K as a function of the density. The line is a linear fit to the data.

calculated thermal conductivities or the heat capacities and their respective experimental values. It also holds for the corresponding data of low-molecular-weight fluids.6 The scatter from the least-squares line is much lower for the polymer values of this work than those for the molecular liquids. One should also note that all classical degrees of freedom contribute equally to the heat capacity, irrespective of their frequency or localization, whereas their contribution to the thermal conductivity is more varied and complex, as outlined above. Therefore, the graph of Figure 3 is no more than a correlation, not a systematic mathematical relationship. This manifests itself, for example, in the slope of the least-squares line being close to 1/2 rather than 1, which would be expected if both heat capacity and thermal conductivity depended solely and linearly on the number of degrees of freedom. Yet, it seems that the thermal conductivity and the heat capacity, calculated from classical molecular dynamics, are effected fundamentally by the same type of error due to the neglect of the quantum nature of some fast degrees of freedom, whereas the details are different. 4.3. Influence of Isotropic Compression/Expansion and Uniaxial Stretch. Isotropically compressed and expanded samples have been prepared for the pc-UA model of amorphous polyamide-6,6. The densities range from 0.99 to 1.19 g/cm3 and include the experimentally estimated range of 1.07-1.11 g/cm3. Figure 4 shows that the thermal conductivity increases linearly with the density. This is to be expected of a classical mechanical simulation since increasing the density also increases the degrees

11522 J. Phys. Chem. B, Vol. 111, No. 39, 2007 of freedom per volume of the sample. A similar increase is observed also for the fl-AA model, where the thermal conductivity increases from 0.41 ( 0.04 to 0.45 ( 0.04 W m-1 K-1 between 1.07 and 1.11 g/cm3. For the stretched sample prepared as described in section 3.5, the thermal conductivities parallel to the stretching direction λ|| and perpendicular to it λ⊥ have been evaluated with the flUA model. The resulting λ|| was 0.44 ( 0.02 W m-1 K-1 (37% higher than that for the isotropic sample at the same density), and the λ⊥ was 0.28 ( 0.02 W m-1 K-1 (12% lower than that for the isotropic sample). The anisotropy of the thermal conductivity λ||/λ⊥ was approximately 1.5. As demonstrated by the anisotropy of the gyration tensor (section 3.5), stretching causes polymer chains to orient in the stretching direction. Therefore, in a stretched sample, there are more backbone bonds oriented in the stretching direction than perpendicular to it. A similar increase of the thermal conductivity in the stretching direction and a concomitant decrease perpendicular to it have also been found experimentally by forced Rayleigh scattering on samples of liquid polyisobutylene (MW ∼ 130000) after shearing.24,25 The fact that there is also an increased thermal conductivity in the same direction supports the hypothesis that there are two types of heat transport in a polymer. The faster is solid-like, through-bond, and progresses along the chain via skeletal vibrations or phonons. The slower is liquid-like, through-space, and transfers energy between different chains by collisions. 5. Conclusions The nonequilibrium algorithms employed (the reverse nonequilibrium molecular dynamics method and the dual-thermostat method) have yielded accurate results for the thermal conductivity of glassy polyamide-6,6 at room temperature. For the first time, to our knowledge, the thermal conductivity of an industrial polymer has been calculated by molecular dynamics simulations. Both algorithms are easy to implement and fast; the simulations described in this paper typically required 1 ns of simulation time. This corresponds to a few days of real time, which makes the approach attractive for practical applications. Further studies on different systems (other polymers, other thermodynamic states) are now needed to extend this conclusion and to determine the range of applicability of the algorithms. We tested four different force fields, which differ in the number of degrees of freedom retained. The degrees of freedom were successively eliminated by removing aliphatic hydrogens (united-atom models) and/or constraining bonds. The thermal conductivity is very sensitive to the number of degrees of freedom of the models; a fully atomistic and flexible model overestimates experiment by about 80%, whereas a fully constrained united-atom model yields a value in excellent agreement with the experimental data. This result is consistent with previous findings for molecular liquids like n-hexane6 and with the interpretation suggested earlier.4 At room temperature, the stiff quantum oscillators cannot be excited, and they basically do not contribute to the heat transport; therefore, they should be excluded or frozen in a classical model. Only a full quantum mechanical treatment could unambiguously clarify which degrees of freedom to keep and which ones to exclude. Such a treatment is not available, and if it was, it would be computationally very expensive. On the basis of the experience of this work, we confirm the earlier suspicion that using a united-atom force field (remove all hydrogen vibrations) with constraints (remove the remaining bond vibrations) works well, in practice, for a solid polymer near room temperature. Recent work on

Lussetti et al. some molecular liquids at ambient conditions hinted that the same heuristics performs well also in those cases.6 There appears, however, to be a guideline which can be followed to select and parametrize a force field suitable for thermal conductivities; with all of the caveats of section 4.2, the heat capacity follows a similar trend. A correct value of the heat capacity is a strong hint that the thermal conductivity of the model will also reproduce experiment. This conclusion is still limited, as it has only been confirmed for a small number of systems and only near ambient conditions, where, however, it seems to hold. If confirmed for other systems, it would be of great practical importance because the calculation of the heat capacity through computer simulations, as well as its experimental measurement, is cheap and simple, unlike the calculation or the measurement of the thermal conductivity. It should be kept in mind that force fields, which yield a correct result for the thermal conductivity, will not necessarily describe other properties of the materials under investigation as correctly. This possible limitation, however, does not preclude a better qualitative understanding of the physics underlying thermal conduction. An example is the calculation of the thermal conductivity in a stretched, anisotropic polymer sample; the thermal conductivity along the stretching direction, mainly due to heat transfer by intrachain mechanisms, is bigger than that in the direction perpendicular to it, in which interchain collisions play a larger role. The thermal conductivity of the isotropic sample, as expected, lies in the middle because both mechanisms contribute to it. This result is consistent with experimental results for other drawn polymers26 and is of interest for its physical meaning as well as for applications to fibers. Acknowledgment. This work has been supported by the Priority Programme 1155 “Molecular Simulation in Chemical Engineering” of the Deutsche Forschungsgemeinschaft. Takamichi Terao thanks the Research Center for Computational Science, Okazaki, Japan, for the use of their facilities. The help of Sylvain Goudeau with the preparation of the initial polyamide structures is gratefully acknowledged. References and Notes (1) Evans, D. J.; Morriss, D. J. Statistical Mechanics of Nonequilibrium Liquids; Academic: London, 1990. (2) Hafskjold, B.; Ikeshoji, T.; Ratkje, S. Mol. Phys. 1993, 80, 13891412. (3) Mu¨ller-Plathe, F. J. Chem. Phys. 1997, 106, 6082-6085. (4) Terao, T.; Lussetti, E.; Mu¨ller-Plathe, F. Phys. ReV. E 2007, 75, 057701/1-057701/4. (5) Bedrov, D.; Smith, G. D. J. Chem. Phys. 2000, 113, 8080-8084. (6) Zhang, M.; Lussetti, E.; de Souza, L. E. S.; Mu¨ller-Plathe, F. J. Phys. Chem. B 2005, 109, 15060-15067. (7) Chiritescu, C.; Cahill, D. G.; Nguyen, N.; Johnson, D.; Bodapati, A.; Keblinski, P.; Zschack, P. Science 2007, 315, 351-353. (8) Mu¨ller-Plathe, F.; Bordat, P. In NoVel Methods in Soft Matter Simulations; Karttunen, M., Vattulainen, I., Lukkarinen, A., Eds.; Lecture Notes in Physics; Springer: Heidelberg, Germany, 2004; Vol. 640, p 310. (9) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Di, Nola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (10) (a) Mu¨ller-Plathe, F. Comput. Phys. Commun. 1993, 78, 77-94. (b) Tarmyshov, K.; Mu¨ller-Plathe, F. J. Chem. Inf. Model. 2005, 45, 19431952. (11) Allen M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (12) Mu¨ller-Plathe, F.; Brown, D. Comput. Phys. Commun. 1991, 64, 7-14. (13) Goudeau, S.; Charlot, M.; Vergelati, C.; Mu¨ller-Plathe, F. Macromolecules 2004, 37, 8072-8081. (14) Crawford, R. J. Plastics Engineering; Butterworth-Heinemann: Oxford, U.K., 1998. (15) dos Santos, W. N.; Gregorio, R., Jr. J. Appl. Polym. Sci. 2002, 85, 1779-1786.

Thermal Conductivity of Amorphous Polyamide-6,6 (16) Abu Isa, I. A.; Jodeh, S. W. Mat. Res. InnoVations 2001, 4, 135143. (17) Gaur, U.; Wunderlich, B. B.; Wunderlich, B. J. Phys. Chem. Ref. Data 1983, 12, 65-89. (18) Wilhoit, R. C.; Dole, M. J. Phys. Chem. 1953, 57, 14-21. (19) Warfield, R. W.; Kayser, E. G.; Hartmann, B. Makromol. Chem. 1983, 184, 1927-1935. (20) Allen, P. B.; Feldman, J. L.; Fabian, J.; Wooten, F. Philos. Mag. B 1999, 79, 1715-1731. (21) Orbach, R. Philos. Mag. B 1992, 65, 289-301.

J. Phys. Chem. B, Vol. 111, No. 39, 2007 11523 (22) Bordat, P.; Sacristan, J.; Reith, D.; Girard, S.; Gla¨ttli, A.; Mu¨llerPlathe, F. Chem. Phys. Lett. 2003, 374, 201-205. (23) Mark, J. E. Polymer Data Handbook; Oxford University Press: New York, 1999. (24) Balasubramanian, V.; Bush, K.; Smoukov, S.; Venerus, D. C.; Schieber, J. D. Macromolecules 2005, 38, 6210-6215. (25) Venerus, D. C.; Schieber, J. D.; Balasubramanian, V.; Bush, K.; Smoukov, S. Phys. ReV. Lett. 2004, 93, 098301. (26) Choy, C. L.; Wong, Y. W.; Yang, G. W.; Kanamoto, T. J. Polym. Sci., Part B: Polym. Phys. 1999, 37, 3359-3367.

Nonequilibrium Molecular Dynamics Calculation of the Thermal

Recommend Documents