Heat transfer through covalent versus non-covalent bonding: a case study on crystalline π-conjugated P3HT polymer using approach-to-equilibrium molecular dynamics

Abstract

Heat transport in a crystalline polymer, poly(3-hexylthiophene) (P3HT), is studied using molecular dynamics. The potential energy is described by density functional theory (DFT) with an additional correction to account for van der Waals dispersion forces. The approach-to-equilibrium molecular dynamics (AEMD) method is used to create and analyse thermal transients and target transport along or between polymer chains. As expected, the approach leads to a length-dependent thermal conductivity of the system along the covalent polymer backbone. Less expectedly, we observe heat transfer between two polymers, i.e. along a non-covalent heat path, and quantify this transfer using a thermal conductance.

---

1 Introduction

Thermal transport by conduction is fundamentally associated with the concept of phonons, which are derived from the force constants between atoms. In cases where bond strengths are not as substantial as in covalent bonds, the mechanisms underlying heat transfer are less well-established. Recent theoretical work has concentrated on the case of two-dimensional materials and the transfer associated with van der Waals forces between graphene and hexagonal boron nitride, for example, as cited in ref. 1. The impact of van der Waals forces has also been considered in the total thermal conductivity calculated for a disordered system of small molecules.² In the same vein, the heat path was experimentally adjusted to pass through more or fewer hydrogen bonds between molecules, thereby increasing the thermal conductivity of the system.³

In order to provide a valuable addition to the existing literature on this subject and enhance the understanding of the relationship between the nature of bonds and the intensity of heat transfer, we study in the present work a π-conjugated crystalline system, namely, poly(3-hexylthiophene) (P3HT). In this system (see Fig. 1), polymers are stacked by π–π interactions, thereby combining covalent character along the polymer chains and non-covalent character between them. In contrast to the amorphous blends of polymers or molecules, where entities may be more or less distant from each other, crystallinity goes hand in hand with coupling along the entire polymer.

Fig. 1 Unit cells in the xy-, xz- and yz-planes of models: (A) reference model containing two P3HT strands, with each strand consisting of eight thiophenes; (B) model A duplicated along the stacking direction (z); (C) two P3HT strands of twelve thiophenes; (D) two P3HT strands of sixteen thiophenes, i.e. model A duplicated along the strand axis (x); (E) two strands of eight thiophenes in which the hexyl side chains have been replaced by methyl groups (P3MT). Carbon atoms are represented in brown, sulfur in yellow, and hydrogen in white. The configurations are captured during a trajectory at T = 300 K.

Thermal studies have already been conducted on crystalline polymer systems, in particular to identify a hydrodynamic transport regime.⁴ In the present study, we expand upon existing models limited to an ideal crystal at T = 0 K by considering strands of varying lengths to probe the length dependence and at finite temperature to account for thermal fluctuations, but still maintaining a π-stacked arrangement. This configuration aligns with the aligned polymer systems (P3HT in particular) employed for the development of organic thermoelectric generators,⁵ for which a comprehensive understanding of thermal transport is crucial in the optimisation of the overall device performances.

The P3HT crystal is studied at room temperature using first-principles molecular dynamics (FPMD). The description of interatomic forces in this scheme, as determined by the calculation of the electronic structure, is independent of empirical assumptions as done in classical molecular dynamics. The density functional theory (DFT), augmented here with the consideration of van der Waals forces, is selected to facilitate the incorporation of interactions that are weaker than covalent or ionic forces, as is the case in the direction of π-stacking. We resort to FPMD to create thermal transients and exploit them in the specific framework of approach-to-equilibrium molecular dynamics (AEMD).

The study is organised as follows: in Section 2, crystalline P3HT systems and the simulation method are described. In Section 3, thermal conductivities characterising heat transport along polymer chains are calculated and analysed. Section 4 focuses on heat transport between π-conjugated polymers and the calculation of the corresponding thermal conductance. Conclusions are drawn in the final section.

2 Methods and system models

Molecular dynamics (MD) simulations are performed using the Car–Parrinello (CP) method⁶ as implemented in developers version 4.3 of the CPMD code.⁷ For the exchange–correlation component of the total Kohn–Sham energy, the formulation proposed by Perdew, Burke and Ernzerhof (PBE)⁸ is employed. The valence–core interactions are described by norm-conserving Troullier–Martins⁹ pseudopotentials. The valence states are developed on a plane wave basis with a 70 Ry cutoff, and the Brillouin zone sampling is limited to the Γ point. The fictitious electronic mass is taken to be equal to 400 a.u. and the time step is set to Δt = 0.1 fs (4 a.u.) to ensure that constants of motion are conserved. The temperature of the ions is controlled using a Nosé–Hoover thermostat,^10,11 while the fictitious electronic kinetic energy is regulated using a Blöchl–Parrinello thermostat¹² to ensure that adiabaticity is preserved. The masses of thermostats (in units of corresponding frequencies) are equal to 1100 cm⁻¹ for the ions and 5000 cm⁻¹ for the fictitious electronic degrees of freedom.

The study of thermal properties is conducted through the implementation of temperature transients within the system, employing the AEMD method.¹³ This approach involves the application of local thermostats to establish hot and cold blocks, constituting “Phase 1” of the process. Subsequently, the thermostats are released, and the temperature of the blocks is monitored throughout the approach-to-equilibrium process, during “Phase 2”. The characteristic time of the transient is then utilised to determine the thermal conductivity or the interface conductance depending on the configuration, as detailed in Sections 3 and 4.

In order to estimate size effects, supercells of P3HT crystals of varying dimensions were constructed and the configurations are captured at room temperature and are shown in Fig. 1. Periodic boundary conditions are applied in the three directions. The size and number of atoms are specified in Table 1. The size along the x direction is systematically adjusted so that the leftmost atom of the chain forms a covalent bond with the rightmost one. This means that an infinitely long polymer chain is obtained. The reference system (A) contains two strands of eight thiophenes each, with the hexyl chains included to allow the structure to naturally adopt the experimentally observed ‘Form I’ configuration. In this configuration, the backbone formed by the thiophenes is inclined with respect to the stacking direction (z). This system is obtained in accordance with the protocol that was employed in our previous work.¹⁴ To summarise, this entails commencing from a planar configuration (not tilted) and subsequently introducing temperature through gradual heating to ambient temperature. System B is duplicated along the stacking directions and contains four strands per unit cell. Systems C and D are analogous to system A, but with larger sizes along the polymer axis x. Finally, in system E, hexyl chains have been substituted by methyl groups.

Table 1 Characteristics of the five models (number of atoms N and size of the orthorombic supercell along x (L_x), y (L_y) and z (L_z))

Label	N	Lx (Å)	Ly (Å)	Lz (Å)
A	400	31.0	15.7	7.6
B	800	31.0	15.7	15.2
C	400	46.6	15.7	7.6
D	800	62.0	15.7	7.6
E	160	31.0	7.3	7.1

---

For each of the aforementioned systems, the electronic structure is calculated using a convergence criterion of 10⁻⁶ a.u. The geometry is then relaxed in the corresponding potential energy, and the system is heated to T = 10 K, followed by T = 100 K, prior to being stabilised for approximately 10 ps at T = 300 K.

3 Heat transport along polymers and thermal conductivity

3.1 Assessment of the thermal conductivity

A regime of thermal transport along the polymer chains is initiated by subjecting the systems to an initial non-equilibrium state, as depicted in the upper section of Fig. 2. The x position of atoms within the skeleton serves as the basis for delineating two distinct blocks, which are then subjected to thermostats at hot and cold temperatures, respectively. The alkyl chain atoms are kept, via the use of a thermostat, at the temperature of the backbone atom to which they are attached. The number of atoms in each block thus formed is identical. Local thermostats are maintained for approximately 10 ps to stabilise the hot and cold temperatures, as illustrated in Fig. 3. The corresponding temperature profile is depicted in ochre in Fig. 2. The resulting shape, which is not uniform in each block, is due to the small size of the system. In the second phase, thermostats are removed, leading to a decline in the temperature of the hot block and an increase in the temperature of the cold block (Fig. 3). This process continues until the temperature of the system is uniform. The duration of the approach-to-equilibrium is system-dependent, being specifically influenced by its thermal conductivity and size. This relationship can be substantiated by solving the heat equation under the same conditions as employed in the molecular dynamics simulation, i.e. a 1D heat transport and initial periodic hot/cold conditions.¹³ The result is a Fourier series, with a predominant term as follows:where T_av is the average temperature reached at the end of phase 2, and ΔT₀ is the initial temperature difference between the hot and cold block, equal to 200 K in the present case.

Fig. 2 Temperature profile along the polymer axis during phases 1 and 2 and fit by a sinusoidal curve. Each point in the profile is obtained by averaging the kinetic energies of the atoms belonging to the interval centred at x. P3HT, 4 polymers (system B).

Fig. 3 Temperatures in the hot (red curve) and cold (turquoise curve) blocks during phase 1 and phase 2. System B is shown in the inset with the delimitations of the two blocks.

As illustrated in Fig. 2, the temperature profile during phase 2 assumes a sinusoidal form, as indicated by the perfect superposition of the pink curve with the averaged profile over the duration of the phase 2 (shown in blue). By calculating the temperature difference ΔT(t) between the hot and cold blocks analytically using the temperature profile from eqn (1), we obtain

As demonstrated in Fig. 4, the present simulations yield the same exponential decay of ΔT(t). It follows from the aforementioned discussion that the relationship obtained from the heat equation between the thermal conductivity κ and the transient decay time τ will be employed in order to calculate the thermal conductivity of the system studied using molecular dynamics. The relationship in question is as follows:

with the heat capacity equal to C_v = 3γNk_B. N denotes the number of atoms in the supercell of volume V, and the factor γ = 1.01 ± 0.01 is calculated from the dependence of the total energy on the temperature during the trajectories simulated for P3HT (0.86 ± 0.11 for P3MT).

Fig. 4 Semi-log plot of the temperature difference ΔT between the hot and cold blocks during phase 2 and exponential fit to extract the decay time τ, here equal to 5.17 ± 0.02 ps, which corresponds to a thermal conductivity equal to κ = 0.21 ± 0.01 W K⁻¹ m⁻¹ (system B).

In order to evaluate the statistical error in the value of τ and therefore in the corresponding κ, one has to take the average over different implementations of phase 2. This amounts to extending phase 1, selecting a new starting point for phase 2 fully decorrelated from the first and extracting a new value of τ. This procedure can be repeated at will for each system under consideration to obtain distinct values for the thermal conductivity, with each one of them equally contributing to the average. In our case, three uncorrelated phase 2 processes have been created.

3.2 Length dependence of the thermal conductivity

Fig. 5 summarises the calculated thermal conductivities of systems A, C, D and E. The results are presented as a function of the length of the polymer chains, which also corresponds to the size of the simulation cell according to x and to the period of the temperature profile. The figure shows the results of successive AEMDs as small dots, with the large dots and error bars representing the mean and standard deviation of the calculations, respectively.

Fig. 5 Thermal conductivity of P3HT versus length (systems A, C and D), without hexyl chains (P3MT, system E), and range of the experimental data.

A comparison of the thermal conductivities of systems A and E, which are identical in length but differ in the presence of hexyl (A) or methyl (E) chains, reveals that the mean values are in perfect agreement, with the error bar being higher in the case of P3MT (system E). This last point can be explained by the much smaller number of atoms in system E (160 compared with 400 for P3HT). The observed agreement can be interpreted as indicating a negligible role of hexyl chains in heat transport along polymers, suggesting that covalent interactions between atoms in the polymer backbone are predominantly responsible for this process.

In contrast, for the diverse P3HT systems with lengths ranging from 3 to 6 nm, a substantial variation in thermal conductivity is observed, which increases with the length of the simulation cell. This size dependence of thermal conductivity has been observed in all the other systems studied by AEMD. To comprehend the physical origin of this effect, it is necessary to revert to the origin of thermal conduction, namely, heat carriers, or phonons, in crystalline materials. According to the principles of kinetic theory, thermal equilibrium is established by collisions. In this case, these collisions involve phonon–phonon scattering and possibly scattering of phonons by defects in the material or surfaces and interfaces. Consequently, heat transport in a material is attributable, on the one hand, to the contribution of the ballistic transport mode of phonons (between two collisions) and, on the other hand, to the diffusive transport of phonons (where collisions dominate). In ballistic transport, phonons evolve without dispersion, and this evolution is limited only by the size of the crystal. Consequently, ballistic thermal conductivity increases with the crystal size in the direction of heat propagation.¹⁵

In the diffusive mode, heat transport is predominantly constrained by collisions. Consequently, diffusive thermal conductivity is a function of the range of phonon mean free paths existing in the system. In a given material, these free paths can span a broad range of values (from 10 nm to 10 μm in silicon, for example¹⁶). Within the AEMD framework, for a simulation cell of dimension L_x in the direction of heat flow, phonons responsible for transporting heat from the hot block to the cold block can be categorised into two distinct groups. The first category comprises phonons with mean free paths shorter than L_x, which undergo collisions as if they were in an infinite crystal. Conversely, phonons with mean free paths greater than L_x will not undergo collisions within the cell and their trajectory will be ballistic. As the size of the simulation cell increases, the proportion of diffusive phonons increases until the size exceeds the largest phonon mean free path of the material. At this juncture, the thermal conductivity ceases to be contingent on the dimensions of the simulation cell, a phenomenon that has been observed, for instance, in crystalline Ge₂Sb₂Te₅.¹⁷ In such a scenario, the thermal conductivity of the material attains a value observed at the macroscopic scale. These phenomena were modelled analytically by Alvarez and Jou,¹⁸ and the resulting law was used in our previous works to extrapolate the AEMD calculations to infinite size. However, extrapolation would not make sense in the present case because the sizes studied are too small compared to the saturation value.

The thermal conductivity of P3HT in the direction of the polymer chains has been determined through experimental methods and documented in a variety of works in the existing literature. For instance, as illustrated in Fig. 5, the horizontal grey band corresponds to the value of κ as determined by Degousée et al.¹⁹ in a mechanically aligned P3HT film. The experimental value is attained by our calculations for system D of length equal to L_x = 6 nm. In addition, the size of the crystalline domain is equal to 8–9 nm for P3HT aligned by this method.²⁰ However, owing to the computational workload of the present calculations, the study of systems of this length was not possible and the observation of the saturation of the thermal conductivity at a constant value was not feasible. In any case, it should be borne in mind that comparison with experimental results should take into account the contribution of the amorphous domains of smaller thermal conductivity, which adds to that of the crystalline domains to form the average macroscopic thermal conductivity of P3HT.

3.3 Other dependencies

In addition to the variation in thermal conductivity as a function of system length in the direction of thermal transport previously discussed, other dependencies were investigated using system B. This system is twice the size in the z direction compared with system A. Initially, the hot and cold blocks are delineated in the x direction, as was previously employed to obtain the results presented in Fig. 5. The thermal conductivities of these systems are presented in Fig. 6, along with those of system D, which is twice as large in the x direction of heat transport. It is evident that the thermal conductivity exhibits a substantially lesser variation when the size is increased along the cross-section than when the size is increased in the direction of thermal transport (B vs. D). A further comparison of the thermal conductivities of systems A and B reveals that the error bar decreases significantly as the cross-section increases. This is due to the fact that there are more atoms in each of the hot/cold blocks and therefore fewer fluctuations in their average temperature. In summary, the effect of the size of the cross-section can be summarised as an impact on the error bar, but not significantly on the thermal conductivity itself.

Fig. 6 Thermal conductivity of P3HT: reference (system A) and impact of the section (system B) and length (system D), and the heat transport direction (system B).

System B was also employed to ascertain the thermal conductivity in the direction perpendicular to the polymer axis. In order to do this, the top two polymers were thermalised to T = 400 K, and the bottom two were thermalised to T = 200 K. In this case, heat propagation during the approach to equilibrium (phase 2) can no longer take place via paths including covalent bonds. Instead, it must pass through π–π interactions between chains. We therefore expect the thermal conductivity to be lower. This hypothesis is confirmed by the results presented in Fig. 6, whichs show a significant decrease in thermal conductivity (from 0.28 to 0.028 W K⁻¹ m⁻¹) in the system under investigation when compared with the reference, system A. It is important to note that size effects may also be present in this direction. Ideally, the thermal conductivity values should be compared when the κ(L) curves saturate. However, the scope of the present study does not permit a full exploration of this issue. Nevertheless, the conclusion remains unambiguous: the presence of a gap between polymers does not impede heat transport, although it is less efficient, which provides the first indication of the role of non-covalent interactions, in this case π–π interactions. This topic is more specifically addressed in the following section.

4 Heat transport between polymers and thermal conductance

4.1 Assessment of the thermal conductance

A regime of thermal transport between the polymer chains is initiated by subjecting the systems to an initial non-equilibrium state, as depicted in the right section of Fig. 7. Local thermostats are maintained for approximately 6 ps to form a hot and a cold block, as illustrated in Fig. 8. The corresponding temperature profile is depicted in ochre in Fig. 7, and the evolution of blocks temperature with time is shown in Fig. 8. In phase 2, thermostats are removed. The temperature profile maintains a similar shape, with one polymer at a higher temperature than the other (Fig. 7), but a decline in the temperature of the hottest one and an increase in the temperature of the coldest one are observed (Fig. 8). Once again, this is a sign that there is heat transfer across the gap separating the two polymers, which is therefore via π–π interactions. The transfer continues until the temperature of the system is uniform. In what follows, we will examine how this transient can be used to characterise this transfer by means of a thermal conductance G.

Fig. 7 Temperature profiles along the z axis during phases 1 and 2 (left) and illustration of the hot and cold blocks, i.e. the upper polymer is hot, and the lowest one is cold (right) (system D).

Fig. 8 Temperatures in the hot (red curve) and cold (turquoise curve) polymers during phase 1 and phase 2. System D is shown in the inset with the delimitations of the two blocks.

Firstly, it should be noted that even in the event of transfer between the two polymers, this process is likely to be less efficient than heat diffusion along the polymer. Consequently, the assumption that the temperature remains consistent throughout the polymers is substantiated and moreover corroborated by our simulations. This approximation is commonly referred to as the lumped capacitance approximation.²¹ It has been employed in first-principles molecular dynamics to analyze heat transfer between two silicon blocks connected by a molecular layer.²² The difference in the present work is that there are two channels of exchange between the hot and cold blocks, i.e. with the polymer above and below in the direction of the π stacking. In such conditions, transient phase 2 is governed by the equilibrium between the internal energy variation of each polymer (the left-hand term in the below equation) and the exchanges with the two neighbouring polymers (the right-hand term), i.e. for the hot polymer:

and for the cold polymer:where G is the thermal conductance, in units of W K⁻¹, and T_H and T_C are the temperature of the hot and cold polymers, respectively. By subtracting eqn (4) and (5), we obtain an equation governing the evolution of the temperature difference ΔT(t) = T_H(t) − T_C(t) between the hot and cold blocks:

The temperature difference therefore decreases exponentially with time according to

ΔT(t) = ΔT₀ exp(−t/τ), (7)

with the decay time constant τ being related to the conductance G by

As demonstrated in Fig. 9, the present molecular dynamics simulations yield the same exponential decay of ΔT(t). It follows from the aforementioned discussion that the above relationship (eqn (8)) will be employed to calculate the thermal conductance G between the polymers from the transient decay time τ obtained by molecular dynamics.

Fig. 9 Semi-log plot of the temperature difference ΔT between the hot and cold polymers during phase 2 and exponential fit to extract the decay time τ, here equal to 6.04 ± 0.03 ps, which corresponds to a thermal conductance equal to G = 1.37 ± 0.01 nW K⁻¹ (system D).

4.2 Impact of lateral chains and the length dependence

The thermal conductances are presented in Fig. 10. The results are presented as a function of the length of the polymer chains, which also corresponds to the size of the simulation cell according to x. The figure shows the results of successive AEMDs as small ochre dots, with the large black dots and error bars representing the mean and standard deviation of the calculations, respectively. The magnitude of interchain thermal conductance is found to be lower in P3MT (system E) in comparison to P3HT (system A). However, the distance d defined as the shortest distance between the polymer chains is equivalent in both systems.¹⁴ This observation indicates that the hexyl chains enhance heat transfer between the polymers, likely by introducing steric interactions that strengthen the coupling between the chains.

Fig. 10 Thermal conductances between chains of P3MT (system E) and P3HT (systems A and D).

Furthermore, in the case where the length of the polymers is doubled (system D), the thermal conductance doubles compared with the reference system (A). This increase is simply an effect of scale: doubling the length L_x results in a corresponding doubling of the effective cross-section of the interchain interactions, thereby increasing the thermal conductance. Moreover, if we normalise the conductance by the cross-section, which is generally done experimentally, we obtain values of 138 ± 20 MW K⁻¹ m⁻² for a length L_x = 3 nm and 135 ± 3 MW K⁻¹ m⁻² for L_x = 6 nm. The other effect of increasing length is that the uncertainty in the conductance value is reduced, which can be attributed to the reduction in temperature fluctuations in systems containing more atoms.

According to our calculations, heat transfer is not limited to the backbone of the polymer where atoms are covalently bonded. In this scenario, it can be hypothesised that the transfer is enabled by π–π interactions, with a potential contribution from alkyl chains. Here we can assume that it is the π–π interactions that allow the transfer, with possibly a complementary effect provided by the alkyl chains. It is anticipated that the thermal conductances thus obtained will be lower than those reported in the existing literature, which pertain to interfaces between covalently bonded materials. For instance, a thermal conductance of 700 MW K⁻¹ m⁻² was measured at the TiN/Al₂O₃ and TiN/MgO interfaces by the time-domain thermoreflectance (TDTR) measurement.²³ At the Al/Si interface, the conductance varies between 125 and 200 MW K⁻¹ m⁻² depending on the roughness of the interface and the preparation conditions.²⁴ A value of 300 MW K⁻¹ m⁻² has been measured at the PMMA/Si interface.²⁵ At the CuPc/C₆₀ interface, the thermal conductance is 450 ± 130 MW K⁻¹ m⁻².²⁶ Furthermore, Losego et al.²⁷ have shown that the thermal conductance increases as the chemical bond becomes stronger, as in the case of the Au/quartz interface. These values remain lower than the thermal conductance at the metal–metal interface, which is, for example, 4 GW K⁻¹ m⁻² at the Al/Cu interface.²⁸

5 Conclusions

In this article, we have presented the study of heat transfer along and between P3HT polymer chains. The calculations were performed using thermal transients and the AEMD method. Regarding thermal conductivity along the chain axis, we observed a size effect, whereby conductivity increased with the chain length. The value obtained for the longest system (6 nm) aligns well with experimental data obtained for a mechanically aligned P3HT film. The effects of dimensions perpendicular to thermal transport, studied by stacking four polymers, are negligible. We also demonstrated that thermal conductivity remains unaffected when hexyl chains are replaced by methyls. For thermoelectric applications, the calculated thermal conductivity for the smallest size considered (3 nm) suggests that reducing the size of the crystalline domains could halve the thermal conductivity and increase the figure of merit.

We have also demonstrated that heat transfer does take place even in the absence of covalent bonding. This was first demonstrated by calculating the thermal conductivity in the π–π stacking direction and confirmed in the second part of the study by tackling heat transfer between two polymer chains thermalised at different temperatures. We evidenced that the temperature of the hot chain decreased and that of the cold chain increased, indicating heat exchange. We characterised these exchanges by calculating the thermal conductance, which, as expected from an exchange not involving strong bonds, has a low value compared with measurements in the literature. Furthermore, we demonstrated that hexyl chains significantly increase conductance in P3HT compared to P3MT. There is no size effect when the polymer chain length is increased, apart from an increase in the exchange surface area.

Author contributions

C. O. Diarra: investigation, software, validation, and visualization. C. Massobrio: writing – review and editing. E. Martin: supervision, methodology, and writing – original draft.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the SI.

Supplementary information available: Atomic configurations of systems A to E. See DOI: https://doi.org/10.1039/d5cp02084k

Acknowledgements

The authors acknowledge the High Performance Computing Center of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. Part of the computing resources were funded by the Equipex Equip@Meso project (Programme Investissements dAvenir) and the CPER Alsacalcul/Big Data. Calculations on the larger systems were performed by using resources from GENCI (Grand Equipement National de Calcul Intensif) (grant no. 0910296 and 0905071). This work of the Interdisciplinary Institute HiFunMat, as part of the ITI 2021–2028 program of the University of Strasbourg, CNRS and Inserm, was supported by IdExUnistra (ANR-10-IDEX-0002) and SFRI (STRATUS project, ANR-20-SFRI-0012) under the framework of the French Investments for the Future Program.

References

1. H. M. Dong, H. P. Liang, Z. H. Tao, Y. F. Duan, M. V. Milosevic and K. Chang, Phys. Chem. Chem. Phys., 2024, 26, 4047–4051.
2. T. Zhou, Z. Li, Y. Cheng, Y. Ni, S. Volz, D. Donadio, S. Xiong, W. Zhang and X. Zhang, Phys. Chem. Chem. Phys., 2020, 22, 3058–3065.
3. G.-H. Kim, D. Lee, A. Shanker, L. Shao, M. S. Kwon, D. Gidley, J. Kim and K. P. Pipe, Nat. Mater., 2015, 14, 295.
4. Z. Zhang, Y. Ouyang, Y. Guo, T. Nakayama, M. Nomura, S. Volz and J. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2020, 102, 195302.
5. N. J. Pataki, S. Guchait, B. Jismy, N. Leclerc, A. Kyndiah, M. Brinkmann and M. Caironi, Adv. Energy Mater., 2025, 2404656.
6. R. Car and M. Parrinello, Phys. Rev. Lett., 1985, 55, 2471–2474.
7. CPMD, Copyright IBM Corp. 1990–2022, copyright MPI für Festkörperforschung Stuttgart 1997–2001. Available on GitHub under MIT License. https://github.com/OpenCPMD.
8. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
9. N. Troullier and J. L. Martins, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993–2006.
10. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
11. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
12. P. E. Blöchl and M. Parrinello, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 9413–9416.
13. E. Lampin, P. L. Palla, P.-A. Francioso and F. Cleri, J. Appl. Phys., 2013, 114, 033525.
14. C. O. Diarra, M. Boero, E. Steveler, T. Heiser and E. Martin, Phys. Chem. Chem. Phys., 2023, 25, 15539–15546.
15. R. Saito, M. Mizuno and M. S. Dresselhaus, Phys. Rev. Appl., 2018, 9, 024017.
16. A. S. Henry and G. Chen, J. Comput. Theor. Nanosci., 2008, 5, 1–12.
17. I. Bel-Hadj, M. Guerboub, A. Lambrecht, G. Ori, C. Massobrio, E. Martin and ADynMat consortium, J. Phys. D: Appl. Phys., 2024, 57, 235303.
18. F. X. Alvarez and D. Jou, Appl. Phys. Lett., 2007, 90, 083109.
19. T. Degousée, V. Untilova, V. Vijayakumar, X. Xu, Y. Sun, M. Palma, M. Brinkmann, L. Biniek and O. Fenwick, J. Mater. Chem. A, 2021, 9, 16065–16075.
20. M. Brinkmann, J. Polym. Sci., Part B: Polym. Phys., 2011, 49, 1218–1233.
21. T. L. Bergman, Fundamentals of heat and mass transfer, John Wiley & Sons, 2011.
22. T.-Q. Duong, C. Massobrio, G. Ori, M. Boero and E. Martin, J. Chem. Phys., 2020, 153, 074704.
23. R. M. Costescu, M. A. Wall and D. G. Cahill, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 67, 054302.
24. P. E. Hopkins, L. M. Phinney, J. R. Serrano and T. E. Beechem, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 085307.
25. M. D. Losego, L. Moh, K. A. Arpin, D. G. Cahill and P. V. Braun, Appl. Phys. Lett., 2010, 97, 011908.
26. Y. Jin, C. Shao, J. Kieffer, K. P. Pipe and M. Shtein, J. Appl. Phys., 2012, 112, 093503.
27. M. Losego, M. Grady, N. Sottos, D. G. Cahill and P. V. Braun, Nat. Mater., 2012, 11, 502–506.
28. B. C. Gundrum, D. G. Cahill and R. S. Averback, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 245426.

---