Gota Kikugawa1,a), Tapan G. Desai2, Pawel Keblinski3, and Taku Ohara1
Journal of Applied Physics 114, published online July 16, 2013
1Institute of Fluid Science, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan
2Advanced Cooling Technologies, Inc., 1046 New Holland Ave., Lancaster, Pennsylvania 17601, USA
3Department of Materials Science and Engineering, Rensselaer Polytechnic Institute, 110 8th Street, Troy, New York 12180, USA
(Received 18 April 2013; accepted 23 June 2013; published online 16 July 2013)
We performed molecular dynamics (MD) simulations on amorphous polyethylene (PE) and polystyrene (PS) in order to elucidate the effect of crosslinks between polymer chains on heat conduction. In each polymer system, thermal conductivities were measured for a range of crosslink concentration by using nonequilibrium MD techniques. PE comprised of 50 carbon atom long chains exhibited slightly higher conductivity than that of 250 carbon atom long chains at the standard state. In both cases for PE, crosslinking significantly increased conductivity and the increase was more or less proportional to the crosslink density. On the other hand, in the PS case, although the thermal conductivity increased with the crosslinking, the magnitude of change in
thermal conductivity was relatively small. We attribute this difference to highly heterogeneous PS based network including phenyl side groups. In order to elucidate the mechanism for the increase of thermal conductivity with the crosslink concentration, we decomposed energy transfer into modes associated with various bonded and non-bonded interactions. VC 2013 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4813505]
I. Introduction
Crosslinking in polymer materials has been extensively studied in order to modify mechanical, chemical and thermal properties1–6 and applied not only to the industrial products like rubber or thermosetting resins but also to the biological materials7,8. For instance, crosslink formation in polyethylene (PE) enables to overcome the disadvantages such as low melting point, stress cracking tendency, and softening and flowing at elevated temperature. Although abundant studies concerning the mechanical behavior of crosslinked polymers exist, there seems relatively little research on the effect of crosslinking on thermal transport properties6,9–11. So far, thermal conductivity in crosslinked polymers was mostly examined focusing on the engineering aspects via experimental measurements and phenomenological models in amorphous polymers12. However, these investigations have not been based on the microscopic factors and mechanisms such as the contribution of crosslink bond between polymer chains to macroscopic thermal conductivity. Therefore, in the present paper, we aim to elucidate how a crosslink formation affects the thermal transport from the molecular viewpoint.
In terms of thermal conduction of the polymeric materials, it is of critical importance to understand the thermal transport properties in an alkyl chain, which is a basic component of a backbone chain in most polymers. Bulk polymer materials comprised of alkyl chains typically have low thermal conductivity on the order of ~0.1 Wm-1K-1. This low value is associated with strong phonon scattering along curved chains and weak interactions between non-bonded chain segments. However, for straight chains such as present in highly drawn PE, conductivity can be over two orders of magnitude higher13,14. Recently, using MD simulations, Henry et al. demonstrated that the thermal conductivity of a straight single PE chain even exceeds ~100Wm-1K–113,14. Wang et al. also showed extremely high thermal conductance of an alkyl chain of SAM15. In these studies, ballistic thermal transport inside an alkyl chain plays a critical role in determining the high thermal conductivity. By using MD simulations, Nakano et al. evaluated high thermal conductance in a single alkyl chain of a lipid molecule inside bilayer lipid membranes16. Sasikumar and Keblinski showed a result of MD simulations that introducing only a small fraction of a kink conformation in a PE chain can drastically reduce thermal conductance due to the scattering of phonons that were previously ballistic17. As for thermal conductivity in bulk polymers, Shenogin et al. examined the effect of anharmonicity of atomic interactions in amorphous polystyrene (PS)18. Ni et al. reported thermal conductivity of crystalline PE (completely aligned PE chains) by using MD simulations19. They also addressed the effect of crosslink on thermal conductivity and confirmed the decrease in thermal conductivity even with a low crosslink concentration due to the scattering of ballistic phone propagation.
In our study, to gain a quantitative understanding of the relationship between crosslinking and thermal conduction, MD simulations on model amorphous polymers including PE and atactic PS were performed. We found that thermal conductivity increases more or less linearly with crosslink concentration in the amorphous PE. On the other hand, in the PS case, the variation in thermal conductivity is small for all the crosslink concentrations studied. Our analysis suggests that this striking difference in thermal conductivity change is partly attributed to the fact that side chains (phenyl groups) in PS introduce heterogeneity in terms of the covalent bond interaction into the polymer network and disrupt the efficient thermal paths through the alkyl backbone network.
Moreover, in the PE systems, we analyzed microscopic energy transfer modes obtained from the decomposition of macroscopic heat flow20,21 in order to elucidate the underlying mechanism behind the crosslink bond contribution to the energy transfer. We decomposed the total heat flux into components associated with molecular translation, the energy exchange by a nonbonded interaction, and the energy exchange by a covalent bonded interaction through polymer chains and crosslink bonds.
II. Computational Details
A. Simulation Systems and Molecular Models
To obtain a well-relaxed structure of amorphous polymers as an initial configuration, the following procedure, which is similar to the proposed one in the literature22, was performed. Polymeric chains were arranged into the simulation cell with three-dimensional periodic boundary condition such that the torsion angles reproduce the RIS (rotational isomeric state) probability23 using XENOVIEW software package24. Then, the system was subject to two successive annealing cycles between 300 K and 600 K with a 1 ns run for each cycle. After that, the system was equilibrated in the NPT ensemble at 1 atm and 300 K.
In the present study, we examined two types of PE models with different hydrocarbon chain lengths comprised of 50 and 250 carbons in a single chain (referred to as C50 and C250 hereafter), and the total number of chains in the system was 200 and 40, respectively. In the atactic PS model, a single chain contains 50 carbon atoms in an alkyl backbone and the total number of chains is 90. The system sizes were 41.33 x 41.331 x 65.31 Å3 and 40.92 x 40.92 x 163.68 Å3 for C50 and C250 PE, respectively, and 50.4 x 50.4 x 151.1 Å3 for PS. These system sizes were decided by confirming that the size was large enough that the size effect on thermal conductivity can be neglected. For example, for C50 PE without crosslinking, thermal conductivities in the 41.33 x 41.33 x 82.66, 41.33 x 41.33 x 165.31, and 41.33 x 41.33 x 247.97 Å3 systems were 0.13, 0.18, and 0.18 Wm-1K-1, respectively. Therefore, the system size of 41.33 x 41.33 x 165.31 Å3 was adopted in the present study. The monomer molecular structure of both PE and PS is depicted in Fig. 1. Also, a snapshot of an equilibrated configuration for the C50 PE system is shown in Fig. 2.
To configure crosslinked polymers, based on noncrosslinked equilibrated systems, atomic pairs with an interatomic distance within 4 Å for PE and 5 Å for PS were randomly selected and linked. In both the PE and PS polymer, the crosslinks were formed only between the different polymer backbone carbons. In particular, the side chains were not connected by the crosslink for PS. Thus, the crosslink in this study are bonded via single sp3 bonds between carbon atoms. The crosslinked polymers were also subject to a 1 ns annealing cycle and equilibrated in the NPT ensemble at 1 atm and 300 K.
The number of crosslink bonds introduced into the system was determined using a degree of crosslinking (DC).
FIG. 1. Chemical structure of a monomer unit for (a) polyethylene and (b) polystyrene.
The definition of DC is the ratio of the number of crosslinked monomers to the total number of monomers, which is described as:
DC = 2NCL/Nmono,
where the NCL and Nmono denote the number of crosslinks and the total number of monomers, respectively. A factor of two in numerator of the right hand side means that each crosslink is formed by two monomers in the different chains. All simulation conditions including equilibrated mass densities are listed in Table I. The equilibrated densities in the non-crosslinked C50 and C250 PE systems were found to be 0.83 and 0.85 g/cm3, respectively, which are slightly lower than other MD results22,25 ranging from 0.87 to 0.91 g/cm3, and experimental data, 0.91 to 0.94 g/cm3 for low density polyethylene. Non-crosslinked polystyrene exhibits the equilibrated density of 0.99 g/cm3 which is fairly close to the experimental value, 1.04 g/cm326. The density in all the polymers increases with a degree of crosslinking, but a dependency on a degree of crosslinking is much weaker for the PS system. This is partly because in the PE systems the short crosslink bonds between polymer chains enable interchain distances to be closer. This packing effect is not clearly seen in the PS systems because originally the large side chain of phenyl rings would be the dominant factor to determine the system density.
FIG. 2. A snapshot of the C50 polyethylene system and a NEMD setup to impose constant heat flux. Each molecule (chain) is differently colored for clarity.
TABLE I. Simulation conditions and equilibrated densities for (a) C50 polyethylene (b) C250 polyethylene, and (c) polystyrene.
The force field for PE molecules was modeled by the united-atom NERD potential27. For the crosslinking bonds between the PE chains, the NERD potential for branched alkanes was also applied28,29. For the PS molecules including crosslink bonds, we used the PCFF all-atom force field30. Both the cutoff distances of 12-6 Lennard-Jones (LJ) interaction for the NERD and 9–6 LJ intermolecular interaction for the PCFF were set to 14 Å.
The time integration of equations of motion for all the systems were performed by multiple time scale r-RESPA scheme31 with a global time step of 1 fs and an inner time step for an intramolecular motion of 0.2 fs. The MD simulations were carried out using the in-house MD program for the PE systems and the LAMMPS MD simulator32 for the PS systems.
B. NEMD Heat Transfer Simulations for Thermal Conductivity Estimation
In the nonequilibrium MD (NEMD) simulations, the constant heat flux was imposed along a longer side of the simulation cell (z axis). To this end, the same amount of kinetic energy was added at constant rate to the heat source region with a slab thickness of 5.0 Å located at the end of the cell and was subtracted from the heat sink region located at the center of the cell33 (see also Fig. 2). The magnitude of heat flux in the PE systems was chosen from 450 MW/m2 to 900 MW/m2, so that a temperature difference between the heat sink and heat source region was sufficiently large (over around 10 K). For the PS systems, a constant heat flux of 600 MW/m2 was used in all simulations. After the system reached a steady state, a temperature profile was measured from a total kinetic energy of all atoms inside the slabs with a thickness of around 2.0 Å generated along the z axis. Thermal conductivity in each system can be evaluated by this temperature profiles, using the Fourier’s law.
In order to investigate the heat conduction mechanism associated with the effect of crosslink formation on thermal conductivity, we monitored the microscopic building blocks of the macroscopic heat flow for PE in the NEMD simulation systems, i.e., the total heat flux is decomposed into the microscopic heat transfer modes. The detailed description for this decomposition can be found in Ref. 20, but here we briefly summarize the methodology. The formula for multibody potentials, including angle bending potential and torsional potential, is given by:
The first term in the right hand side of Eq. (2) represents the transport of internal energy of sites carried by their motion. The quantities concerning site s, vz,s and Es, are the z component of velocity and the sum of potential and kinetic energies, respectively. The symbol Σ is the summation of all sites in a measurement volume. In this study, the measurement volume is assigned as a slab region, which is located between the energy control volumes. The second term in the right hand side of Eq. (2) represents the contributions of the work performed by inter- and intramolecular forces on the sites, which are defined by n-body potentials including bond stretching, angle bending, and dihedral potential. The vector Fs,U is the intra- or intermolecular forces acting on the site s due to the interaction of an n-body potential defined for a set of sites U={s1, s2,…, sn}, and zs is the position z of site s. The double summation ΣΣ in the bracket [], which is summed over all the n-body potentials, is taken over all pairs of sites under the condition that either the line segment connecting the two sites is contained in the measurement volume or it crosses one or both boundaries of the measurement volume. The quantity (zsα – zsβ) represents the portion of (zsα – zsβ) involved in the measurement volume. Jz and VCV are the total heat flux and measurement volume, respectively. In summary, the total energy flux is mainly decomposed into three parts, i.e., the contributions of energy transfer associated with molecular motion, energy exchange by bonded interactions through covalent bond inside an alkyl chain, and energy exchange by non-bonded interactions including the van der Waals (vdW) interaction expressed by LJ potential. Taking into account the percentages of each contribution to the total heat flux, the thermal conductivity can also be similarly decomposed into these three contributions that corresponds to each energy transfer mode.
III. Results and Discussion
A. Thermal Conductivity of Linear and Crosslinked Polymers
Temperature profiles obtained from the NEMD simulations in the non-crosslinked PE and PS systems are shown in Fig. 3. The temperatures exhibit more or less linear profiles between the heat source and sink regions. The thermal conductivities in both the PE and PS as a function of a degree of crosslinking are plotted in Fig. 4.
Thermal conductivities in the non-crosslinked polymers (DC 0% in Fig. 4) were found to be 0.18Wm-1K-1 for C50 PE, 0.27Wm-1K-1 for C250 PE, and 0.16Wm-1K-1 for PS. Although PE has a wide variety of forms and properties depending on a synthesis method, thermal conductivity of PE typically exhibits ~0.3Wm-1K-1 (e.g., see Ref. 34). Therefore, our result in the longer PE case shows reasonable agreement with experimental data. The reason for larger thermal conductivity in the PE with a longer alkyl chain is that the proportion of covalent bond interaction to take part in thermal energy transfer is relatively larger in the C250 PE system as described in Sec.III B. The thermal conductivity in the non-crosslinked PS also is in good agreement with experimental data, ~0.16Wm-1K-134, and the value is nearly half of that for PE.
In the crosslinked PE systems, thermal conductivity more or less linearly increases as the DC increases. However, in the low DC range, up to 10%, there is no visible change in thermal conductivity in both the PE systems. We speculate that this change in thermal conductivity originates from how the increase of number of covalent bonds affects thermal energy transport in the crosslinked PE. The crosslinks could influence thermal conductivity in both positive and negative ways, i.e., on one hand the addition of strong covalent bonds increases thermal conductivity between prior non-bonded chain segments, but on the other hand, crosslinks can decrease thermal conductivity due to the scattering of phonon propagation along the chain backbone19. With small DC, these two contributions more or less cancel each other, while with larger DC positive contribution dominates.
FIG. 3. Temperature profiles obtained from NEMD simulations on amorphous non-crosslinked (a) C50 PE, (b) C250 PE, and (c) PS.
FIG. 4. Thermal conductivities in the C50 PE (blue marked line), C250 PE (green marked line), PS (red marked line) as a function of a degree of crosslinking. Filled marks denote thermal conductivity in the 80% DC crosslinked PE with mass heterogeneity.
In the PS system, the variation of thermal conductivity with increasing DC is insignificant as compared with that of the PE systems, although thermal conductivity slightly increases up to around 0.2 Wm-1K-1. We speculate that this difference is attributed to highly heterogeneous PS structure including phenyl groups and a mix of sp3 and sp2 carbon bonds. Such “disordered” structures with heterogeneous interaction strength or molecular conformation could exhibit low thermal conductivity, which apparently does not increase significantly with increasing DC.
In order to give more credence of the above explanation on the different effects of DC on thermal conductivity of PE and PS, we examined the effect of mass heterogeneity on thermal conductivity using the PE systems. The heterogeneity of intramolecular interactions is present in the PS system due to its phenyl side chain (sp3 and sp2 mixed structure of covalent bonding) in contrast to the PE system. Instead of introducing interaction heterogeneity into the system which could significantly affect the polymer conformation, mass heterogeneity was imposed in the present study. To this end, a certain amount of atomic mass is randomly added or subtracted to/from every carbon atom so as to maintain total mass of the whole system. In our simulations, 50% of carbon atoms were 6C and the rest were 18C. Both the C50 and C250 PE in the 80% DC condition were examined. The results of thermal conductivity were plotted in Fig. 4. It was revealed that the thermal conductivity decreases by ~0.1 Wm-1K-1 from that in the uniform mass case. While this trend is consistent with our reasoning, the understanding of the very low thermal conductivity of crosslinked PS requires further study.
B. Microscopic Contributions to Thermal Conductivity
As described in Sec. II B, we evaluated the microscopic heat transfer mode to understand the underlying mechanism of the increase in thermal conductivity with a crosslink concentration in the PE systems.
The decomposed contributions in the non-crosslinked PE systems are shown in Fig. 5. The contribution are divided into 3 categories: (i) The energy transfer associated with molecular motion (first term in Eq. (2)), (ii) the energy transfer by nonbonded interactions including the contribution of intra- and intermolecular vdW interaction, and (iii) the energy transfer by bonded interactions. They are labeled as “translation”, “nonbonded,” and “bonded,” respectively. The contribution of nonbonded interactions is similar for C50 and C250 chains, while the difference in the contribution of the bonded interactions is significant. Therefore, the difference in thermal conductivity between non-crosslinked C50 and C250 PE systems mainly comes from the bonded interactions. The longer chain PE can transfer more thermal energy through its longer bonded alkyl chain without heat conduction interrupted by chain ends.
FIG. 5. Contributions of energy transfer modes to thermal conductivity in amorphous C50 PE and C250 PE. The red, green, and blue boxes denote the contributions of molecular translation (first term in Eq. (2)), nonbonded interactions, and bonded interaction, respectively.
In the cross-linked PE cases, the energy transfer modes were evaluated in all the DC conditions (Fig. 6). In both systems, the change in the contribution of nonbonded interactions is more or less insignificant except that in the 80% DC. On the other hand, the bonded contribution increases with increasing crosslink concentration. This means that introducing crosslink bonds opens new heat paths and directly contributes to the increase in thermal conductivity. Therefore, we conclude that the main contribution to the larger thermal conductivity in highly cross-linked PE is attributed to the bonded interactions of crosslink bonds. In the PS system, increasing DC might not change the percentage of each contribution enough.
FIG. 6. Contributions of energy transfer modes to thermal conductivity in (a) C50 and (b) C250 PE systems with varying the degree of crosslinking. The colors of boxes denote the same contributions as those in Fig. 5.
IV. Conclusion
In the present study, we performed molecular dynamics simulations on amorphous polyethylene and polystyrene with varying crosslink concentrations to investigate the effect of crosslink bonds on heat conduction. For amorphous PE system, two different hydrocarbon chain lengths comprised of 50 and 250 carbons in a single chain was examined. It was found that thermal conductivity in the system with a longer PE chain is larger and similar to the experimental value. The thermal conductivity was found to more or less linearly increase with crosslink concentration.
The decomposition of total heat flux into each contribution of microscopic energy transfer mode reveals that the increase in thermal conductivity is dominated by the contribution of energy exchange by bonded interaction. Thus, introducing crosslink bonds open new heat channels, which directly contribute to the increase in thermal conductivity. The variation of thermal conductivity in atactic PS with a crosslink concentration is insignificant as compared with that in PE. This is likely due to highly heterogeneous PS structure in terms of the covalent bond interaction strength and presence of the phenol rings.
Acknowledgments
We want to thank Dr. William J. Evans for the technical supports and useful discussions. This work was supported by the Young Researcher Overseas Visits Program for Vitalizing Brain Circulation by the Japan Society for the Promotion of Science (JSPS). Numerical simulations were performed on the SGI Altix UV1000 at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University. TGD was supported by the NASA Small Business Innovation Research Grant (SBIR), under Contract No. NNX11CA05C.
References
1I. Chodak, “Properties of crosslinked polyolefin-based materials,” Prog. Polym. Sci. 20, 1165–1199 (1995).
2L. Calcagno and G. Foti, “Ion irradiation of polymers,” Nucl. Instrum. Methods B59/60, 1153–1158 (1991).
3H. A. Khonakdara, J. Morshediana, U. Wagenknechtb, and S. H. Jafari, “An investigation of chemical crosslinking effect on properties of highdensity polyethylene,” Polymer 44, 4301–4309 (2003).
4H. Lin, E. Van Wagner, J. S. Swinnea, B. D. Freeman, S. J. Pas, A. J. Hill, S. Kalakkunnath, and D. S. Kalika, “Transport and structural characteristics of crosslinked poly(ethylene oxide) rubbers,” J. Memb. Sci. 276, 145–161 (2006).
5I. Krupa and A. S. Luyt, “Mechanical properties of uncrosslinked and crosslinked linear low-density polyethylene/wax blends,” J. Appl. Poly. Sci. 81, 973–980 (2001).
6B. Likozar and M. Krajnc, “Cross-linking of polymers: Kinetics and transport phenomena,” Ind. Eng. Chem. Res. 50, 1558–1570 (2011).
7G. Lewis, “Properties of crosslinked ultra-high-molecular-weight polyethylene,”Biomater 22, 371–401 (2001).
8S. M. Kurtz, M. L. Villarraga, M. P. Herr, J. S. Bergstr€om, C. M. Rimnac, and A. A. Edidin, “Thermomechanical behavior of virgin and highly crosslinked ultra-high molecular weight polyethylene used in total joint replacements,” Biomater 23, 3681–3697 (2002).
9D. N. Simavilla, J. D. Schieber, and D. C. Venerus, “Anisotropic thermal transport in a crosslinked polyisoprene rubber subjected to uniaxial elongation,” J. Poly. Sci. B50, 1638–1644 (2012).
10V. T. Morgan, “The thermal conductivity of crosslinked polyethylene insulation in aerial bundled cables,” IEEE Trans. Elect. Insul. 26, 1153–1158 (1991).
11B. Likozar and M. Krajnc, “A study of heat transfer during molding of elastomers,”Chem. Eng. Sci. 63, 3181–3192 (2008).
12C. Zhong, Q. Yang, and W. Wang, “Correlation and prediction of the thermal conductivity of amorphous polymers,” Fluid Phase Equilibria 181, 195–202 (2001).
13A. Henry and G. Chen, “High thermal conductivity of single polyethylene chains using molecular dynamics simulations,” Phys. Rev. Lett. 101, 235502 (2008).
14A. Henry and G. Chen, “Anomalous heat conduction in polyethylene chains: Theory and molecular dynamics simulations,” Phys. Rev. B 79, 144305 (2009).
15Z. Wang, J. A. Carter, A. Lagutchev, Y. K. Koh, N.-H. Seong, D. G. Cahill, and D. D. Dlott, “Ultrafast flash thermal conductance of molecular chains,” Science 317, 787–790 (2007).
16T. Nakano, G. Kikugawa, and T. Ohara, “Molecular heat transfer in lipid bilayers with symmetric and asymmetric tail chains,” J. Heat Trans. 135, 061301 (2013).
17K. Sasikumar and P. Keblinski, “Effect of chain conformation in the phonon transport across a Si-polyethylene single-molecule covalent junction,” J. Appl. Phys. 109, 114307 (2011).
18S. Shenogin, A. Bodapati, P. Keblinski, and A. J. H. McGaughey, “Predicting the thermal conductivity of inorganic and polymeric glasses: The role of anharmonicity,” J. Appl. Phys. 105, 034906 (2009).
19B. Ni, T. Watanabe, and S. R. Phillpot, “Thermal transport in polyethylene and at polyethylene–diamond interfaces investigated using molecular dynamics simulation,” J. Phys.: Condens. Matter 21, 084219 (2009).
20D. Torii, T. Nakano, and T. Ohara, “Contribution of inter- and intramolecular energy transfers to heat conduction in liquids,” J. Chem. Phys. 128, 044504 (2008).
21T. Ohara, T. C. Yuan, D. Torii, G. Kikugawa, and N. Kosugi, “Heat conduction in chain polymer liquids: Molecular dynamics study on the contributions of inter- and intramolecular energy transfer,” J. Chem. Phys. 135, 034507 (2011).
22A. P. Awasthi, D. C. Lagoudas, and D. C. Hammerand, “Modeling of graphene–polymer interfacial mechanical behavior using molecular dynamics,” Model. Simul. Mater. Sci. Eng. 17, 015002 (2009).
23P. J. Flory, Statistical Mechanics of Chain Molecules(Hanser, 1989).
24S. Shenogin and R. Ozisik, see http://xenoview.mat.rpi.edu/ For Xenoview Ver. 3.7 (2009).
25D. Hossain, M. A. Tschopp, D. K. Ward, J. L. Bouvard, P. Wang, and M. F. Horstemeyer, “Molecular dynamics simulations of deformation mechanisms of amorphous polyethylene,” Polymer 51, 6071–6083 (2010).
26A. Quach and R. Simha, “Pressure-volume-temperature properties and transitions of amorphous polymers; polystyrene and poly (orthomethylstyrene),” J. Appl. Phys. 42, 4592–4606 (1971).
27S. K. Nath, F. A. Escobedo, and J. J. de Pablo, “On the simulation of vapor–liquid equilibria for alkanes,” J. Chem. Phys. 108, 9905–9911 (1998).
28S. K. Nath and J. J. de Pablo, “Simulation of vapour-liquid equilibria for branched alkanes,” Mol. Phys. 98, 231–238 (2000).
29S. K. Nath and R. Khare, “New forcefield parameters for branched hydrocarbons,” J. Chem. Phys. 115, 10837–10844 (2001).
30H. Sun, S. J. Mumby, J. R. Maple, and A. T. Hagler, “An ab initio CFF93 all-atom force field for polycarbonates,” J. Am. Chem. Soc. 116, 2978–2987 (1994).
31M. Tuckerman, B. J. Berne, and G. J. Martyna, “Reversible multiple time scale molecular dynamics,” J. Chem. Phys. 97, 1990–2001 (1992).
32S. Plimpton, “Fast parallel algorithms for short-range molecular dynamics,” J. Comput. Phys. 117, 1–19 (1995); see also http://lammps.sandia.gov.
33P. Jund and R. Jullien, “Molecular-dynamics calculation of the thermal conductivity of vitreous silica,” Phys. Rev. B59, 13707–13711 (1999).
34X. Zhang and M. Fujii, “Measurements of the thermal conductivity and thermal diffusivity of polymers,” Poly. Eng. Sci. 43, 1755–1764 (2003).