Next Article in Journal
Inhibition of Cutaneous TRPV3 Channels by Natural Caffeic Acid for the Alleviation of Skin Inflammation
Next Article in Special Issue
Design of pH/Redox Co-Triggered Degradable Diselenide-Containing Polyprodrug via a Facile One-Pot Two-Step Approach for Tumor-Specific Chemotherapy
Previous Article in Journal
Modification of the Protein Amino Acid Content in Hen Eggs as a Consequence of Different Concentrations of Lupine and Soy in Feed
Previous Article in Special Issue
Atomistic Origins of Resurrection of Aged Acetylcholinesterase by Quinone Methide Precursors
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Designing Accurate Moment Tensor Potentials for Phonon-Related Properties of Crystalline Polymers

1
Institute of Solid State Physics, NAWI Graz, Graz University of Technology, 8010 Graz, Austria
2
Institute of Materials Chemistry, TU Wien, 1060 Vienna, Austria
*
Author to whom correspondence should be addressed.
Molecules 2024, 29(16), 3724; https://doi.org/10.3390/molecules29163724
Submission received: 13 June 2024 / Revised: 16 July 2024 / Accepted: 29 July 2024 / Published: 6 August 2024
(This article belongs to the Special Issue Exclusive Feature Papers on Molecular Structure)

Abstract

:
The phonon-related properties of crystalline polymers are highly relevant for various applications. Their simulation is, however, particularly challenging, as the systems that need to be modeled are often too extended to be treated by ab initio methods, while classical force fields are too inaccurate. Machine-learned potentials parametrized against material-specific ab initio data hold the promise of being extremely accurate and also highly efficient. Still, for their successful application, protocols for their parametrization need to be established to ensure an optimal performance, and the resulting potentials need to be thoroughly benchmarked. These tasks are tackled in the current manuscript, where we devise a protocol for parametrizing moment tensor potentials (MTPs) to describe the structural properties, phonon band structures, elastic constants, and forces in molecular dynamics simulations for three prototypical crystalline polymers: polyethylene (PE), polythiophene (PT), and poly-3-hexylthiophene (P3HT). For PE, the thermal conductivity and thermal expansion are also simulated and compared to experiments. A central element of the approach is to choose training data in view of the considered use case of the MTPs. This not only yields a massive speedup for complex calculations while essentially maintaining DFT accuracy, but also enables the reliable simulation of properties that, so far, have been entirely out of reach.

Graphical Abstract

1. Introduction

Crystalline polymers have a peculiar structure: they are characterized by strong covalent bonds along their polymer chains and weak van der Waals (vdW) bonds in directions perpendicular to them. This leads to highly anisotropic characteristics, especially for phonon-related properties and heat transport, which are in the focus of our interest. A relevant reference material in this context is polyethylene (PE), which is also the structurally simplest polymer. In its common, amorphous form, PE is a thermal insulator with a very low thermal conductivity, κ, from around 0.3 to 0.5 Wm−1K−1 [1,2]. When its crystallinity is increased, its thermal conductivity increases dramatically. In fact, in (partly) crystalline PE thin films, thermal conductivities ranging from 22.5 Wm−1K−1 [3] and ~40 Wm−1K−1 [4] to even 62 Wm−1K−1 [5] and 65 Wm−1K−1 [6] have been measured, with the actual values depending on the quality of the samples. The highest thermal conductivity values have been observed in stretched fibers. For these, κ reaches values between 90 Wm−1K−1 [7] and 104 Wm−1K−1 [8].
To better understand the underlying mechanisms of thermal conduction and to know its theoretical limits, a reliable simulation of κ for a perfectly ordered single crystal is extremely helpful. To achieve this, a reliable simulation of harmonic and anharmonic phonon properties is required. The knowledge of phonon characteristics also provides insight into vibrational spectra, elastic material properties, and temperature-related quantities like thermal expansion. The precondition for correctly describing phonon properties is an accurate simulation of the forces acting on the polymer constituents upon displacing specific atoms. A complication in this context is that, when calculating anharmonic phonon properties, (tens of) thousands of individual force calculations need to be performed, especially when dealing with extended, complex systems containing many atoms in the considered supercells. Provided that such force calculations are possible at a sufficiently high level of accuracy, the material’s thermal conductivity can, for example, be calculated via the Boltzmann transport equation (BTE) by simulating phonon heat capacities, group velocities, and lifetimes [9]. Alternatively, with high-quality forces at hand, the heat transport properties can be calculated from particle trajectories, employing approaches such as non-equilibrium molecular dynamics (NEMD) [10], the Green–Kubo method [11,12], or by performing approach-to-equilibrium molecular dynamics (AEMD) [13] simulations. The number of required simulations is, however, even larger in such molecular dynamics (MD)-based simulations.
In either case, the repeated calculation of forces on atoms becomes so computationally costly that it is often infeasible when employing ab initio methods like density functional theory (DFT). Therefore, the community often resorts to empirical force fields, which can, however, be very inaccurate, especially when the phonons of interest describe intra-molecular vibrations. This can, for example, be concluded from the large spread in the thermal conductivities along polyethylene chains calculated with such “off-the-shelf” empirical force fields. Depending on the nature of the force field, the thermal conductivity along PE chains has been predicted to amount to ~50 Wm−1K−1 [14] (using the COMPASS force field), ~45 Wm−1K−1 [15] (employing the AIREBO potential) or 310 ± 190 Wm−1K−1 [16] (with the REBO potential). Notably, the former values are even significantly below the highest measured ones, although in experiments, the measured thermal conductivity is inevitably reduced by sample imperfections, which extrinsically enhances phonon scattering. Some of us have also recently shown that empirical force fields like COMPASS and GAFF fail at accurately reproducing the phonon band structure of molecular crystals like naphthalene, especially in regions in which intra-molecular vibrations dominate [17]. This is in sharp contrast to DFT (in particular PBE with the Grimme D3 or the many-body dispersion (MBD) correction [18,19,20]), which yields phonon band structures [17] and low-frequency Raman spectra [21] in excellent agreement with experiments. Thus, overall, van-der-Waals-corrected DFT accurately describes phonon-related properties, but is slow; conversely, transferrable, empirical force fields are fast, but often suffer from inaccuracies.
Recently developed machine-learned (force field) potentials (MLPs) promise to resolve this dilemma, especially when they are parametrized specifically for the material in question [22,23,24]. The latter can be achieved by training them on accurate DFT data. The type of MLPs in the focus of this paper are moment tensor potentials (MTPs), as implemented in the Machine Learning Interatomic Potentials (MLIP) code [25,26]. Their conceptual structure will be briefly discussed below. Since, to the best of our knowledge, MTPs have not been systematically used for polymers, we focus here on devising optimized strategies for their parametrization and benchmarking. This serves as a crucial prerequisite for using them to systematically understand the heat transport and related phenomena in polymers in follow-up studies.
To benchmark the obtained potentials, we will compare a range of observables calculated with the MTPs to the results obtained when using the same DFT approach that is used in the MTPs’ parametrization. When experimental data for certain properties are available as well, these will also be considered in the benchmarking. In the latter case, however, one has to keep in mind that a comparison to experiments not only tests the quality of the parametrization of the MTPs, but also the suitability of the underlying ab initio method (e.g., DFT).
In the context of thermal and phonon-related properties, MTPs have been shown to correctly predict phonon band structures and thermal conductivities of simple crystals (diamond, silicon, InAs, and BAs) and 2D materials [24]. Also in another study, MTPs have correctly predicted the phonons of 2D materials and the MTP-calculated thermal expansion of graphene agreed reasonably well with the DFT results [27]. Moreover, the MTP-derived phonon band structures of BC2N monolayers [28], diamond [23], and BAs [23] agree very well with the corresponding DFT results. Moreover, in recent studies, some of us showed that MTPs also provide an exceptionally accurate description of phonon-related properties for metal–organic frameworks (MOFs) [29] and, when suitably parametrized, can even reproduce the spin-dependent vibrational properties of MOFs [30]. Very recently, MLPs have been used to accelerate the thermal conductivity calculations of 103 inorganic materials [31]. However, these results cannot be directly transferred to crystalline polymers, as, due to their extreme anisotropy, a successfully parametrized MTP must be able to correctly describe particularly strong covalent interactions along the polymer chains and very weak van der Waals and multipole interactions perpendicular to them, potentially in conjunction with highly flexible side chains attached to the polymer backbones. In the present study, we will discuss under which circumstances MTPs can live up to this challenge.
In particular, we will focus on developing strategies for the parametrization of accurate, system-specific MTPs for various use cases. For this, different sets of reference structures for the parametrization are calculated using dispersion-corrected density functional theory, which allows for simultaneously describing covalent, electrostatic, and van der Waals interactions (for details see Section 3). On the one hand, we will benchmark the performance of the obtained MTPs for predicting unit cells, the elastic stiffness tensor, as well as harmonic and anharmonic phonon properties. All these quantities can be calculated in a quasi-static manner considering comparably small atomic displacements. On the other hand, we will also test the accuracy of the suitably parametrized MTPs for performing molecular dynamics simulations at room temperature, at which much larger atomic displacements occur. In this context, we will also assess to what extent it can be beneficial to parametrize distinct MTPs for such complementary use cases.
Structures of Polyethylene, Polythiophene, and Poly(3-hexylthiophene-2,5-diyl) (P3HT)
In the present study, we will investigate three polymers with complementary properties: polyethylene (PE), polythiophene (PT), and poly-3-hexylthiophene (P3HT), which are shown in Figure 1. Throughout the manuscript, these polymers are considered in their crystalline form, rather than in their amorphous state or as a melt. Therefore, in the simulations, periodic boundary conditions are used, which have the practical consequence of the polymer chains not having a beginning or an end. The unit cells periodically repeated in all directions of space are indicated by the black rectangles in Figure 1. PE is characterized by fully saturated single bonds along its chain. This results in a certain degree of flexibility of the polymer backbone. Moreover, it has a rather simple structure, with only 12 atoms in the unit cell, which also allows us to perform properly converged DFT calculations of phonon lifetimes, at least when restricting the description to three-phonon scattering processes (i.e., third-order force constants). Relying on state-of-the-art computational resources, this is no longer feasible for PT (28 atoms in the unit cell) and even less so for P3HT (100 atoms in the unit cell). In that sense, PT and P3HT serve as examples for what becomes possible when using machine-learned potentials. Here, PT is a classical reference system for a semiconducting polymer with a rigid, π-conjugated backbone. P3HT combines the features of PT and PE, as it consists of a rigid PT backbone surrounded by flexible hexyl sidechains that make it soluble in common organic solvents. From a practical point of view, P3HT is also one of the most applied conjugated polymers, as, for a long time, it served as a classical benchmark material for polymer-based organic thin-film transistors [32,33,34,35].
PE is typically observed as one of two polymorphs, being either orthorhombic or monoclinic [36]. The monoclinic polymorph is metastable and reverts back to the orthorhombic form when subject to temperatures in excess of 60–70 °C [36]. Thus, in the following, we will consider only the orthorhombic polymorph of PE, where, as a starting geometry for our geometry optimizations, we picked the known structure from Huan et al. [37]. To the best of our knowledge, the crystal structure of PT has not been unambiguously determined yet. Therefore, we constructed it assuming an orthorhombic unit cell with two chains, in accordance with previous studies [38,39]. The structure was then relaxed with DFT. More details on how the structure was created can be found in Section S1 of the Supplementary Materials. For the starting structure of P3HT, we relied on the work of Zhugayevych et al., which deals with a theory-based identification of the lowest-energy polymorphs [40]. In that work, multiple nearly isoenergetic polymorphs were found, whose energies at room temperature differed by less than kBT. We opted for the lowest-energy polymorph labelled “CAaac” in the paper [40], whose structure and monoclinic unit cell were then relaxed with DFT. The structures were aligned with the chain directions parallel to the Cartesian z-direction, which is parallel to the lattice vector c, in accordance with the literature [36,41]. In addition to the real space structures in Figure 1, Figure 2 illustrates the first Brillouin zones of PE and P3HT. The first Brillouin zone of PT is not shown, as it is analogous to that of PE.
Figure 1. Optimized real space structures of orthorhombic polyethylene (PE) (a,b), of polythiophene (PT) (c,d), and of poly(3-hexyl-thiophene (P3HT) (e,f) with the viewing directions along the direction of the polymer chains (a,c,e) and perpendicular to them (b,d,f). The tripods indicate the cartesian directions x, y, and z. (Color code—brown: carbon; white: hydrogen; yellow: sulfur; plots produced using VESTA [42]).
Figure 1. Optimized real space structures of orthorhombic polyethylene (PE) (a,b), of polythiophene (PT) (c,d), and of poly(3-hexyl-thiophene (P3HT) (e,f) with the viewing directions along the direction of the polymer chains (a,c,e) and perpendicular to them (b,d,f). The tripods indicate the cartesian directions x, y, and z. (Color code—brown: carbon; white: hydrogen; yellow: sulfur; plots produced using VESTA [42]).
Molecules 29 03724 g001
For polyethylene (PE), all of the quantities considered in this work are calculated, because the corresponding DFT calculations are feasible and/or experimental data exist. For polythiophene (PT) and P3HT, selected quantities like the thermal conductivity and thermal expansion will not be benchmarked, as for these, no suitable reference data exist, and their calculation with DFT is beyond current computational possibilities.

2. Results and Discussion

In the following, we will first describe the details of the used MTPs together with the strategy used for their parametrization. Here, we will address the question of how the parametrization strategy should be adapted to the intended use case of the MTP. Subsequently, we will benchmark the performance of the “best” MTPs for predicting structural parameters, phonon bands, elastic constants, and anharmonic properties (like phonon lifetimes and thermal expansion). Finally, their performance in MD simulations and the achievable speedup of the computations will be assessed.

2.1. Generating Moment Tensor Potentials

2.1.1. The Tested Moment Tensor Potentials (MTPs) and Their Training Process

In the following, a general overview of the employed computational approach will be given, which is necessary for assessing the results discussed below. The more technical details will be compiled at the end of the manuscript in Section 3. As mentioned in Section 1, the machine-learned potentials of choice for the present study are moment tensor potentials (MTPs), as implemented in the MLIP code (Version 2) and as described in detail in Refs. [25,26]. In short, in MTPs, the local environments around atoms are identified within spheres with a chosen radius. For all systems discussed below, this radius is set to the default value of 5 Å. We typically find a comparably weak dependence of the performance of the MTPs on this parameter, with the best results typically obtained around the default value. The energy contribution of an atom is then determined by the neighbors within its environment. Summing the contributions of all atoms yields the total energy of the system. The atomic contributions due to the local environments are expanded in terms of basis functions, which are built from moment tensor descriptors. These consist of a radial part represented by a series of polynomials up to a given maximum order (defined by the “radial basis set size” parameter), and an angular part, which is represented by a tensor up to a certain rank. The multitude of combinations of moment tensors included in the potential is defined by a “level” parameter. Increasing the level and the radial basis size increases the number of trainable parameters, which typically increases the accuracy of the MTP, but at the same time, reduces its speed. According to our experience, increasing the level has a stronger impact than increasing the radial basis size, both in terms of speed and accuracy. Therefore, throughout this work, we kept the radial basis size fixed to 10. The influence of the level is investigated by setting it to 18, 22, and 26 (see below).
In this context, it should be mentioned that, for the quantities of interest in the present study, even the highest-level MTPs are many orders of magnitude faster than any DFT simulation. Therefore, whether employing a computationally less efficient MTP could become a problem depends on the situation. This is best discussed in the context of specific use cases: one set of target applications of the MTPs discussed in this work is the simulation of harmonic and anharmonic force constants to obtain phonon band structures and phonon lifetimes. In the case of PE, only 12 force calculations on structures comprising a single, displaced atom are needed when calculating phonon band structures. In contrast, around 8000 force calculations on structures with pairs of displaced atoms are required for determining the phonon lifetimes from third-order force constants (for details see Section 3). The latter constitutes a considerable computational workload when performing the simulations with DFT. In contrast, such calculations are straightforward on any PC, even for the most complex MTPs considered here. Therefore, for such calculations, the speed of the MTP is not really an issue, and very high levels (like level 26) can be chosen. Conversely, when performing molecular dynamics (MD) simulations (here, for evaluating the thermal expansion of PE), the computational demand is usually significantly higher. This is particularly true when studying thermal transport with NEMD and AEMD approaches, where forces for systems comprising tens of thousands of atoms need to be calculated many millions of times. Thus, for MD simulations, it will typically be advisable to use a lower level, like 22, 18, or even 16, which makes the form of the MTPs mathematically less complex with fewer adjustable parameters. This results in a significant speedup of the MTPs (see Section 2.9, Section S27, and especially Section S28, as well as Ref. [29]).
Another relevant aspect of MTPs is that, by default, they treat all atoms of one chemical element equally, employing the same set of parameters. For PE, this is perfectly fine, as all carbons have identical chemical environments. However, if one considers more complex polymers, carbon atoms can appear in significantly different configurations (i.e., hybridizations, chemical environments, etc.). In such cases, it is beneficial to treat them as different atom types. Therefore, for PT, we treat carbon atoms as two atom types, distinguishing between carbon atoms bonded to two other carbon atoms and one hydrogen and carbon atoms bonded to two carbons and one sulfur. In P3HT, the chemical environments of the carbons are even more diverse, as illustrated in Figure 3. A major aspect is that now one is dealing with sp2-hybridized carbon atoms in the backbone and sp3-hybridized ones in the side chains. Describing the differently hybridized carbons by different atom types turns out to be crucial for the accuracy of the MTPs. Additionally, we will also assess to what extent further differentiation between carbon atoms in non-identical chemical environments would be useful. These different atom types are illustrated by the different shadings of red and blue in Figure 3 (see also figure caption). In terms of speed, distinguishing between atom types significantly increases the required computational effort to parametrize the MTPs due to the larger number of parameters that need to be fitted (a quantitative assessment will be provided below). However, a larger number of atom types has no significant impact on the runtime of simulations when using an already parametrized MTP.
Training MTPs means adjusting their variable parameters, such that a cost function is minimized. In MLIP, the cost function comprises the weighted square deviations between the energies, stresses, and forces calculated by the respective MTP and a reference methodology for certain training structures. In our case, that reference methodology is dispersion-corrected DFT employing periodic boundary conditions. In short, the PBE functional [43] is combined with Grimme’s D3 correction [18,19] (with Becke–Johnson damping [20]) to account for long-range van der Waals interactions. For these simulations, converged basis sets and k-point grids are employed and the Vienna Ab initio Simulation Package (VASP) [44] is used (further technical details are provided in Section 3). The generation of the training structures will be described in Section 2.1.2. The different contributions to the cost function are scaled, such that one can give more weight to either energies, stresses, or forces. Here, we use the default weights of 1 (eV)−1, 0.01 Å(eV)−1, and 0.001 (GPa)−1 for energies, forces, and stresses, respectively. Training ends when the cost function drops by less than a factor of 10−3 over the previous 50 iterations.
The random initialization of the fit parameters at the start of the training results in a stochastic nature of the obtained MTPs [26]. Therefore, it is useful to train ensembles of differently initialized MTPs. We typically train five MTPs for any system and level of MTP complexity. These five MTPs are used to calculate specific quantities of interest (i.e., they are benchmarked on a validation set [45]), and then the “best” MTP is picked for further use. For MTPs meant for primarily modeling phonon properties, the selection criterion could, for example, be their ability to reproduce Γ-point frequencies or phonons in the whole Brillouin zone, where the latter is adopted here. Thus, phonon frequencies homogeneously sampled over the entire Brillouin zone are used to calculate the root mean square deviation (RMSD) up to 12.5 THz for choosing the “best” MTP that is meant for phonon calculations. The reason for restricting the frequency range will be explained in Section 2.2, where also the deviations between different parametrization runs will be discussed. A summary of various possible quality metrics for the MTPs is provided in Section S2.
Conversely, for ranking the MTPs meant for MD simulations, a possible criterion that is also chosen here is their performance when describing the forces acting on atoms in a set of validation structures not used in the parametrization (see Section 3 for details on the generation of the validation structures). As an alternative to picking the “best” MTP, one can calculate the quantities of interest with all five MTPs, and then take the mean and standard deviation of the different predictions. This approach was used by the MLIP developers in their original paper [26] and it has the advantage of obtaining an uncertainty estimate by the means of the standard deviation, which can, for example, be compared to the difference between the MTP- and DFT-calculated quantities. Thus, for selected quantities, in the following, mean values will also be specified.

2.1.2. Generating the Training Data for the MTP Parametrization

When generating reference data for the parametrization of MTPs, it is crucial to perform an efficient sampling of the configuration space. For that, different strategies have been employed in the past. To benchmark these approaches in the context of heat transport modeling, Choi et al. [46] generated data for training neural network interatomic potential (NNIP) by ab initio molecular dynamics, by random displacements of atoms, and by generating displacements as a superposition of phonon eigenmodes. They concluded that the best-performing NNIPs were those generated based on the ab initio molecular-dynamics-derived training data. Here, we employ a similar, yet significantly more computationally efficient approach: in short, for all our systems, we perform molecular dynamics simulations within VASP relying on DFT accelerated by on-the-fly parametrized kernel-based force fields [47]. In this approach, VASP uses a number of initial ab initio MD steps to generate a VASP machine-learned potential (VMLP), which is used to predict further MD steps. Depending on whether the Bayesian error estimation lies below or above a threshold, that predicted step is retained or a step using DFT-calculated forces is performed. The threshold of the error can be automatically determined or user-defined (as further discussed in Section S7). Any new DFT calculation then becomes part of the pool of DFT reference data that is available for the training of the VMLP. This process continues, and the pool of reference data grows over time until enough training data have been collected. In our cases, such an active learning run typically occurs for 10,000 to 100,000 MD steps, comprising several hundred DFT simulations (with the actual numbers for each parametrization listed below). Due to the nature and parametrization of the VMLP, this strategy is particularly computationally efficient (compared, e.g., to MLIP’s internal active learning approach), as argued in more detail in [29].
The DFT-calculated energies, forces, and stresses from the VASP on-the-fly learning process are then used as training data for fitting the MTPs. This two-step training approach has, for example, been implemented by the VASP group for modeling interfaces [48]. A detailed comparison of the performances of VMLPs and MTPs in terms of speed and accuracy can be found in [29], suggesting similar accuracies of both potential families and, at the time of the study, somewhat faster simulations with MTPs. Most importantly, the MTPs can be straightforwardly interfaced with popular molecular dynamics codes like LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) [49] using the LAMMPS–MLIP interface. This boosts the versatility of MTPs and makes them particularly suited for dealing with all use cases in the focus of the present study on an equal footing. Therefore, in the following, we will exclusively discuss the performance of the MTPs.
A central question for any MD-based generation of training data is at what temperature the MD simulations should be performed, or whether employing a temperature ramp during the run would be more appropriate. We opt for the latter in line with most previous active-learning endeavors employing VASP [29,48,50]. But also, in that case, a maximum temperature for the heat ramp needs to be set, where it again appears advisable to consider the envisioned use case of the MTPs: in single-point calculations geared at obtaining harmonic phonon properties, all atoms except for one are at their equilibrium positions, while the one atom is typically displaced by 0.01 Å (a value that we checked for convergence, and which is known to provide well-converged results for organic materials [17]). Similarly, when calculating phonon lifetimes via finite displacements of pairs of atoms, the considered atomic displacements are very small. For example, the default displacement value in phono3py [9] (which is also applied here; for convergence tests, see Section S29) amounts to only 0.03 Å. As a consequence, in such calculations, the MTPs are required to accurately describe the very weak forces on atoms in structures that are particularly close to equilibrium. Consequently, for calculating harmonic and anharmonic phonon properties, training data at small displacements (i.e., at low temperatures) ought to be sufficient or even preferable. Thus, for calculating phonon properties, we test the impact of generating training data only up to 100 K. In contrast, when using MTPs for performing MD simulations around room temperature, much larger displacements are encountered. For example, for the polymers considered here, the average atomic displacements in MD runs amount to roughly ~0.13 Å at 15 K and ~0.6 Å at 300 K, as described in Section S8. Thus, when performing MD simulations, especially at room temperature, the MTP must accurately describe forces at comparably larger displacements. This suggests that the use of training data at elevated temperatures—possibly even significantly above the envisioned target temperature—would be advisable. In fact, including high-temperature reference data is a commonly applied strategy for improving the force fields eventually used for molecular dynamics runs [29,51]. Also, the VASP manual suggests overshooting the application temperature in the training by about 30% [51]. Thus, for the MD use case of the present study, training data are sampled up to a temperature of 500 K, which is significantly higher than the maximum temperature of 300 K used when simulating thermal expansion. Notably, at such high temperatures, polymers like PE or P3HT would be melts in their equilibrium configurations. However, the simulation times applied are much too short for the polymers to melt, such that, in the simulations, they remain in their crystalline structure. Moreover, the high-temperature simulations are exclusively performed for sampling the configuration space for parametrizing machine-learned potentials.

2.2. Impact of the Choice of the Reference Data, the Level of the MTP, and the Number of Considered Atom Types

Prior to providing an extensive, quantitative comparison between MTP- and DFT-calculated observables, it is useful to assess the impacts of several factors influencing the parametrization process. This, in particular, applies to the reference data, the MTP level, and the number of considered atom types. Adjusting these settings provides handles for improving the accuracy of MTPs for the two use-cases outlined above, namely for calculating phonon-related properties (using small atomic displacements) and for performing MD runs at finite temperatures (associated with much larger atomic displacements). For assessing the various, differently parametrized potentials, two metrics are defined: to quantify the accuracy of MTPs in predicting phonons, the root mean square deviation (RMSD) between the MTP phonon frequencies and the DFT reference data, RMSDphonon, is calculated on a dense mesh in the whole Brillouin zone up to 12.5 THz throughout the paper (further details, e.g., on the mesh size are given in Section 3). The choice of a cutoff is motivated by the fact that low-frequency phonons are insofar more relevant, as their thermal occupation is higher, and they typically have higher group velocities and phonon lifetimes. This makes them the main contributors to the thermal conductivity. As shown in Section S2, the exact value of the cutoff is not significant. As a second parameter complementary to RMSDphonon, a metric to assess the accuracy of the MTPs in MD simulations is defined. It is referred to as RMSDMD and is calculated as the RMSD between the MTP-calculated forces and DFT forces for a validation set generated via a 300 K active learning run (see Section 3 for details on the generation of the validation set).
With these two metrics (RMSDphonon and RMSDMD) in hand, three strategies for improving the accuracy of the MTP will be assessed. The first is optimizing the generation of the training data, where, for the reasons discussed in Section 2.1.2, the focus is on the impact of the temperature range over which the training data are generated. The “default” training sets for PE and PT originate from active learning MD runs from 15 K to 500 K in an NPT ensemble with 150 MD steps per Kelvin temperature increase. The “default” MTPs are then parametrized at level 22 with one carbon type for PE and two carbon types for PT. To test the impact of emphasizing training data at low temperatures, additionally smaller reference data sets are defined by including only the subsets of DFT-calculated structures obtained in temperature ranges from 15 K to 100 K, 15 K to 200 K, 15 K to 300 K, and 15 K to 400 K. For PE (PT), this yields five training data sets with 129 (173) configurations generated at temperatures ranging from 15 K to 100 K, 260 (260) configurations from 15 K to 200 K, 368 (344) configurations from 15 K to 300 K, 507 (456) configurations from 15 K to 400 K, and the original 550 (490) configurations from 15 K to 500 K. The resulting RMSDphonon are then calculated for five independently parametrized MTPs for PE and PT employing a DFT-relaxed unit cell. The resulting deviations between the MTP and DFT results are shown in Figure 4a for PE (blue) and PT (orange).
As will be discussed in more detail below, the frequency deviations are generally very small. Nevertheless, the random initialization of the MTP fitting procedure leads to a non-negligible spread within each set of five MTPs parametrized with identical settings. The majority of the MTPs within each set result in rather similar and small values of RMSDphonon, with typically one or, at most, two outliers. The RMSDphonon values for the best-performing MTPs are shown in Figure 4a by large, open symbols connected by dashed lines. The values obtained for all other MTPs are shown by small, solid symbols, with the values for the median MTPs connected by solid lines. In view of the outliers obtained in each temperature range, stressing the median values appears more useful than calculating the average values of RMSDphonon. The trends for both the median and the best-performing MTPs show that using only low-temperature training data clearly improves the RMSDphonon for PE as well as for PT. For the best-performing MTPs, in the case of PE, the value of RMSDphonon drops by more than a factor of two from 0.15 THz (for 15–500 K reference date) to 0.07 THz (for 15–100 K reference data). Similarly, for PT, it drops from 0.08 THz to 0.04 THz. This is remarkable insofar as the number of reference configurations used in the parametrization is reduced by a factor of more than two for the restricted temperature range. The improved accuracy arises from the fact that the forces in the phonon calculations are more similar to the forces encountered in the training data sampled at a low temperature than to those occurring during training at higher temperatures (as illustrated in more detail in Section S9). To test whether the “low-temperature” MTPs can be further improved by increasing the number of reference configurations, we repeat the active learning run for PE between 15 K and 100 K, starting with the already obtained data from the run between 15 K and 100 K (to avoid adding very similar configurations). This increases the number of training configurations from 129 to 397. We then reparametrize five new MTPs on the expanded “low-temperature” reference configurations. The resulting RMSDphonon remains virtually unchanged at 0.72 THz for the respective best-performing MTPs, while the median RMSDphonon decreases somewhat from 0.10 THz to 0.08 THz (which appears to be within the statistical fluctuations). Also, fixing the unit cell size in the training (employing an NVT instead of an NPT ensemble) does not yield statistically relevant changes in RMSDphonon (with hardly any change for the best-performing MTPs and some decrease for the median, see also Figure S1). A possible reason for this is the very small thermal expansion of PE below 100 K.
As a next step, the potentials, assessed in terms of their accuracy when evaluating the phonon properties (Figure 4a), are also tested in terms of their accuracy in predicting the forces for the structures encountered in MD runs at 300 K (Figure 4b). The values of RMSDMD follow an opposite trend compared to the values of RMSDphonon: the larger the temperature range considered in the generation of reference configurations, the more accurate the force predictions become. This is not unexpected, considering that the MTPs trained on the data obtained from active learning runs between 15 K and 100 K have never encountered the large atomic displacements occurring in the validation structures at 300 K. In view of this, it might appear surprising that the performance of the “low-temperature” MTPs is still reasonable, but for MD simulations, at elevated temperatures, they fail in a different context: the simulations become unstable, resulting in a disintegration of the polymer chains. This is explicitly tested for the best-performing “low-temperature” MTP of PE in 800 ps NPT runs using supercells containing 3456 atoms for fixed temperatures set to values between 100 K and 340 K in 40 K steps. While the simulations are stable up to 220 K, the polymer chain disintegrates at all tested higher temperatures.
To assess the impact of the level on the MTPs’ accuracy, MTPs trained with data generated between 15 K and 100 K are parametrized at level 18, level 22, and level 26. Figure 4c,d clearly show that both RMSDphonon and RMSDMD decrease with increasing level. In fact, RMSDphonon for the best-performing MTP of this series decreases by more than a factor of two between level 18 and level 26, and also, for RMSDMD, this decrease amounts to 33%. Moreover, while for the level 22 MTPs, one of the stochastically initialized MTPs is still an outlier, generating comparably large deviations, in the present test, this does not occur at level 26. This is consistent with our observation that higher MTP levels tend to decrease the stochastic spread of the MTPs. A higher MTP level, however, comes with higher computational costs both in the parametrization and in production calculations, as discussed in Section 2.1.1. As a rule of thumb, in our experience, increasing the level by two increases the computational demands by roughly a factor of 1.5, as illustrated in Section S27 of the Supplementary Materials.
As a third strategy for increasing the accuracy of the MTPs (especially in complex materials), different atom types can be assigned to identical elements in different chemical environments. P3HT is the ideal test system for the following comparison, considering the significant number of chemically non-identical carbon atoms (see Figure 3). In Figure 5, the RMSDphonon and RMSDMD values for sets of level 22 MTPs trained with reference data between 15 K and 500 K are compared for different degrees of atom typing of the carbons of P3HT: this includes treating all carbons in P3HT equally, introducing three atom types (distinguishing between carbon in the side chains, carbon in the backbone bonded to sulfur, and carbon in the backbone not bonded to sulfur), and using six atom types (distinguishing between all chemically inequivalent carbons, as illustrated by the different shadings of red and blue in Figure 3). Splitting into three types gives a clear improvement compared to a single atom type for both RMSDphonon and RMSDMD. Splitting into six types improves the RMSDphonon and RMSDMD of the best and median MTPs even further. Notably, the arithmetic mean RMSDphonon actually does not improve, but this is the consequence of a statistical fluctuation caused by a particularly poorly performing outlier for the MTP with six atom types (see Table S9). In passing, we note that, also for PT, a significant improvement in the MTP accuracy is observed, when the chemically inequivalent carbon atoms are treated as separate atom types (see Table S8). Speed wise, as already discussed in Section 2.1.1, splitting into atom types does not decrease the speed of the MTPs, but it somewhat increases their training time (by ~60% for three atom types and by ~170% for six types, when using the computing hardware described in Section S28). Thus, maximizing the number of atom types in a material appears advisable, as long as this does not have any adverse consequences on the stability of the MTPs. In fact, when using MTPs for P3HT in MD, the simulations tend to be more stable for three atom types for carbon compared to six.
For designing the MTPs used in the following benchmarking (and also in further calculations), all three discussed strategies can be combined to boost their overall accuracy. For example, for P3HT, when using training data generated between 15 K and 100 K, setting the MTP level to 26, and splitting the carbon atoms into six types, the RMSDphonon for the best-performing MTP (out of two) amounts to 0.04 THz (i.e., 0.8 cm−1). This is remarkably low, considering that P3HT is a rather complex system with 100 atoms in the primitive unit cell. As the above settings (level 26, training between 15K and 100 K, and a maximum number of atom types for carbons) are ideal for phonon calculations, they are suggested for that task. Thus, MTPs parametrized using the said settings in the following will be referred to as MTPphonon. Furthermore, we also define a suggested MTP for MD simulations (abbreviated by MTPMD), which needs to be considerably faster and must also work at elevated temperatures. It is based on training data from 15 to 500 K, has a level of 22, and also maximally split atom types. Amongst all parametrized MTPphonon for a specific material, the one with the smallest RMSDphonon in the following will be referred to as the “best” MTPphonon. Likewise, the MTPMD with the smallest RMSDMD will be considered to be the “best” MTPMD. In the following, the predictive performance of these “best” MTPs will be assessed for a variety of physical observables.

2.3. Predicting Unit Cell Parameters

Phonon frequencies (and other phonon-related properties) strongly depend on the lattice parameters used in their calculation. Therefore, as a first assessment of the accuracy of the machine-learned potentials, static lattice parameters calculated with MTPs and DFT via minimizing the potential energy are compared. Where available, experimental values are also reported. The results can be found in Table 1.
For all materials, an excellent agreement between the lattice parameters obtained with DFT and those with the MTPs is observed. In the vdW bonding directions (lattice vectors a and b), the deviations are around ~1%, while in the chain direction (lattice vector c), the agreement is even better, with deviations below 0.1%. Generally, the MTPsphonon perform slightly better than the MTPsMD, which is not surprising, considering that the former set of MTPs is designed for describing structures close to equilibrium. In Table 1, we also report an uncertainty estimate that we obtain by parametrizing five MTPs and calculating the standard deviations of the predicted lattice parameters. Interestingly, for MTPMD, this uncertainty estimate is typically of a comparable magnitude as the deviation of the mean values of the MTP calculations from the DFT results, whereas for MTPphonon, the uncertainty estimate is smaller.
Importantly, DFT and MTP relaxations not only yield equivalent lattice constants, but also yield the same symmetries for all materials. PE and PT have an orthogonal unit cell (all angles, including α, are 90°). P3HT has a monoclinic symmetry with α, the angle between the b and c lattice vector, being slightly below 90°. When comparing the MTP-simulated unit cell parameters to experiments, the overall agreement is still convincing. For most values, the deviations are, however, clearly larger than the deviations from the DFT structures. This is due to three reasons: (i) there is a non-negligible experimental error manifested in the deviations between different experiments that are sometimes even larger than the deviations in individual experiments from the DFT results, (ii) several of the experiments are performed at elevated temperatures, while, for the simulated values reported in Table 1, thermal expansion is not considered (compare Section 2.7), and (iii) even though PBE + D3 used in the DFT simulations is known to provide a very good description of the structures of molecular crystals, it is certainly not perfect [56,57], and one cannot expect a machine-learned potential to perform better than the ab initio methodology used for its parametrization. Overall, the MTPs are still capable of excellently reproducing the densities of the crystalline polymers.

2.4. Elastic Constants

The elastic constants of polymers are of practical and fundamental relevance, as illustrated, for example, by the exceptionally high measured values of ~200–300 GPa for the Young’s modulus of highly oriented PE in the chain direction [58,59,60]. Moreover, the elements of the elastic stiffness tensor can be regarded as phonon-related properties, since they are intimately coupled to the group velocities of acoustic phonons via the Christoffel equations [61]. For calculating the elastic constants using the MTPs, we applied two approaches: the “clamped ion method”, in which the atomic positions are not relaxed directly for strained unit cells. Instead, the effect of the relaxation is estimated via the Hessian matrix [62] to save computational time. This approach is also the one used internally in VASP for calculating elastic constants. Additionally, for the MTP-based calculations, we also employ the “relaxed ion method”, in which the atomic positions for strained unit cells are relaxed directly. For PE and PT, the resulting elastic constants, calculated with identical MTPs, are very similar, irrespective of the used method, but for P3HT, the two methods give rather different results for all elastic constants apart from Czz (see Section S10). This discrepancy shows that there is an issue with the “clamped ion method” for P3HT. We attribute this to a numerical instability and difficulties in estimating the atomic relaxation via the Hessian matrix discussed in Section S11 of the Supplementary Materials. As the “relaxed ion method”, which does not suffer from this effect, is beyond reach in the DFT calculations of P3HT, the following benchmark for P3HT will only be concerned with Czz. Correspondingly, Figure 6 compares all elements of the elasticity tensor calculated with DFT and with MTPs as a parity plot only for PE and PT. Overall the agreement between DFT and the MTPs is quite satisfactory, especially for MTPphonon, which provides a better agreement than MTPMD. On the one hand, we attribute this to the higher level of MTPphonon. On the other hand, the small unit cell changes in the training data for MTPphonon are closer to the small unit cell changes of 1% that are applied in the calculations of the elastic constants.
To provide a more quantitative comparison, we focus on the Czz element, which is by far the largest element and corresponds to how much stress along the chain is required to deform the material in that direction. Naturally, as there are strong covalent bonds along the chain, Czz reaches large values of 328 GPa, 383 GPa, and 151 GPa for PE, PT, and P3HT when calculated with DFT (with the full elastic stiffness tensor given in Section S12). To understand why Czz is 2.5 times larger in PT than in P3HT, one has to consider that, for the same number of polymer chains in the unit cell, the area of the x-y plane is 2.7 times larger in P3HT than in PT (see Figure 1c,d). Comparing the “best” MTPphonon and DFT, Czz deviates by 1.6%, 0.8%, and 0.1% (for PE, PT, and P3HT), respectively. This can be regarded as an excellent agreement between the two approaches. The “best” MTPMD performs somewhat worse (albeit still very good) with deviations of 3.1%, 2.4%, and 2.7% (for PE, PT, and P3HT, respectively). To put the agreement between MTPs and DFT into perspective, DFT calculations are also performed with the default VASP energy cutoff of 400 eV rather than the converged cutoff of 900 eV (PE) or 700 eV (PT and P3HT). These calculations with the low-energy cutoff yield deviations that are larger than the deviations between the MTPs and DFT calculations with converged cutoffs they were trained on (see Figure S7).
For PE, a comparison to a large number of experiments is also possible, where the list of experimental data in Table 2 is even truncated, omitting X-ray diffraction and Raman spectroscopy results from before 1980. A more complete list is provided in Ref. [63]. Experiments typically only report the elastic stiffness tensor elements Czz and show a significant spread, which could be due to imperfect PE crystals [64]. In fact, especially for Czz, the spread in experimental data is much larger than the deviations between the MTP and DFT results. Notably, however, if one attributes too small experimental values to crystal imperfections and, thus, considers the largest values to be closest to the ideal, intrinsic situation, the agreement of the experiments with both MTP and DFT data is indeed excellent. Regarding Cxx and Cyy, the agreement is clearly already worse for the DFT data.
For PT and P3HT, the literature is sparser than that for PE. In fact, for these materials, values varying by orders of magnitude have been reported for different samples, and the directions for which the elastic constants were determined were often not clear [65,66]. Thus, we shifted the discussion of these data to Section S16 of the Supplementary Materials.
Table 2. Elastic constants Cxx, Cyy, and Czz (also called Young’s moduli) of PE calculated in this study and from the literature for temperatures T (“RT” denotes room temperature). They are calculated with five MTPs (each for MTPphonon and MTPMD), for which the arithmetic mean and standard deviation are calculated. Corresponding values are denoted as “mean” in the table. Additionally, the results with the “best” MTPs are listed.
Table 2. Elastic constants Cxx, Cyy, and Czz (also called Young’s moduli) of PE calculated in this study and from the literature for temperatures T (“RT” denotes room temperature). They are calculated with five MTPs (each for MTPphonon and MTPMD), for which the arithmetic mean and standard deviation are calculated. Corresponding values are denoted as “mean” in the table. Additionally, the results with the “best” MTPs are listed.
First AuthorYearMethodT [K] C x x [GPa] C y y [GPa] C z z [GPa]
Theoretical
This study2024DFT012.211.4328
This study2024MTPphonon, “best”013.212.2322
This study2024MTPphonon, mean013.0 ± 0.612.4 ± 0.5322 ± 1
This study2024MTPMD, “best”013.913.8318
This study2024MTPMD, mean015.8 ± 2.112.8 ± 2.0316 ± 4
Kurita [63]2018DFT010.97.8333
Experimental
Matsuo [67]1986X-ray293 213–229
Nakamae [58]1991X-rayRT 235
117 254
Kobayashi [68]1983RamanRT 281
Tashiro [69]1988RamanRT 260
Pietralla [59]1997RamanRT 315
Holliday [60]1971neutron29866329
Twisleton [70]1982neutron7698326

2.5. Phonon Band Structure

As a first directly phonon-specific property, the MTP and DFT calculated phonon band structures are compared. In Figure 7, we show the results for the best-case scenario, i.e., using MTPphonon for all three systems. This yields a truly amazing agreement between the MTP and DFT results, even for P3HT, despite its highly complex band structure. In this context, it should be noted that the deviations for P3HT in Figure 7c are amplified, as, for the sake of clarity, the plotted frequency range is only half that of the other polymers (with a plot of the full frequency in Figure S9). Importantly, the calculated band structures also agree very well with the experimental data on the dispersion of acoustic phonons in PE along the Γ → Z direction deduced from the overtones of Raman spectra, as discussed in Section S18 in the Supplementary Materials. Another aspect that needs to be stressed for phonon band structures is that excellent agreement is observed not only along the chain direction (which is parallel to Γ → Z), but also for the directions in which inter-chain interactions are dominated by vdW and multipole interactions (parallel to Γ → X and Γ → Y). This illustrates that the MTPs used here are able to simultaneously mimic a highly anisotropic bonding situation that leads to large and small forces occurring in different directions.
In view of the amazing agreement between the DFT- and MTP-calculated phonon band structures, the questions arise of whether this is a direct consequence of choosing the “best” potentials for that task (the MTPphonon potentials with the minimum RMSDphonon) and how the accuracy would change for the MTP with the median RMSDphonon. To quantify the agreement, again, the values of RMSDphonon are used, which are particularly well suited for that task considering that they represent the deviations between the MTP- and DFT-calculated phonon frequencies for a homogeneously sampled Brillouin zone up to 12.5 THz. The corresponding values are summarized in Table 3 for various MTPs. For the MTPphonon used to plot the band structures in Figure 7, they amount to 0.043 THz (1.43 cm−1), 0.029 THz (0.97 cm−1), and 0.036 THz (1.20 cm−1) for PE, PT, and P3HT, respectively, i.e., they are extremely small, even for the structurally highly complex P3HT crystal. This agreement only very mildly deteriorates for other MTPphonon variants, such that, for the MTPs with the median RMSDphonon values, the deviations increase by only ~0.01 THz (also see Figures S2 and S3 in the Supplementary Materials for the actual band structures).
When using MTPMD variants, the agreement becomes clearly worse, with the values of RMSDphonon increasing by roughly a factor of three (see Table 3). Notably, this increase is much larger than that seen between the “best” and the median potentials for all variants. Even worse, the “best” MTPMD of PE produces two modes with imaginary frequencies at the Y-point (shown in Figure S13). This further supports the notion that, when using MTPs for calculating phonon properties, the “phonon” parametrization strategy defined in the last paragraph of Section 2.2 is advisable.
As a last note, concerning low-frequency phonon band structures, it should be mentioned that, for the above-discussed simulations, DFT-optimized unit cells were used in all cases in order to not mix aspects regarding the accuracy of calculating phonon bands with the accuracy of predicting unit cells. We also test the use of an MTP-optimized unit cell for calculating phonon bands with the “best” MTPphonon. Especially for PE and P3HT, this causes only a minor increase in RMSDphonon to 0.053 THz and 0.046 THz, respectively (see also Table S15). For PT, the increase in RMSDphonon is larger (to a value of 0.094 THz, which is similar to the value obtained for the MTPMD with the DFT unit cell). This is insofar not surprising, as for PT, the deviations between the DFT- and MTPphonon-optimized unit cell parameters are also larger (see Table 2). Still, even in this case, the changes in the band structure are rather minor, as shown in Figure S14.
For benchmarking phonons over the entire frequency range, plotting band structures is not ideal, as the resulting plots would become far too busy. Thus, Figures S11–S13 in the Supplementary Materials compare the phonon densities of states for the “best” versions of MTPphonon. Again, the agreement is extraordinary. To illustrate that the excellent agreement between the phonon band structures calculated with DFT and classical force fields cannot be taken for granted, Figure 7d compares the DFT-generated bands with those obtained using a traditional, transferrable force field for PE. Here, we choose the AIREBO potential [71], which has been commonly used for PE, even as recently as 2020 [8,15,72]. While the AIREBO potential does reproduce the overall shape of the phonon band structure of PE, it severely (in some instances, even by a factor of two) overestimates the frequencies of all phonons. This also results in significantly too large group velocities, which makes this force field unsuitable for calculating, for example, thermal conductivities, which will be the topic of Section 2.6.

2.6. Calculating the Thermal Conductivity of PE Using the Boltzmann Transport Equation

After demonstrating that the MTPs provide an excellent performance when calculating harmonic phonon properties, the question arises whether this also applies to anharmonic quantities, like phonon lifetimes. These are, for example, relevant for thermal conductivities. For PE, (with considerable computational efforts) we were able to use DFT with converged numerical settings to calculate the phonon lifetimes from third-order force constants. These phonon lifetimes were then used to calculate the thermal conductivity via the Boltzmann transport equation (BTE). For more complex materials like PT and especially P3HT, performing such calculations is clearly beyond current computational possibilities. Thus, the benchmarking in this section focuses solely on PE. In this context, it also needs to be stressed that the present section deals with a comparison of DFT and MTP results, employing an approach that is consistent with what has been used in the literature [72,73,74]. There are several conceptual shortcomings of this approach that need to be considered for a quantitatively accurate description of thermal conductivity. As they have no direct impact on the benchmarking performed here, they will only be discussed briefly at the end of this section.
The linearized form of the BTE can be solved directly using, for example, phono3py [75]. We refer to this as the “full BTE” solution. Another approach is to solve the BTE within the relaxation time approximation (RTA) [75]. This requires less computational effort and is commonly conducted in the literature, as it typically has comparably little impact on the results for materials with low thermal conductivities [76]. For PE along the chain direction, where thermal conductivities are large, the RTA has been shown to underestimate the thermal conductivity by ~20% as compared to the full BTE [72,73]. Nevertheless, we will also report results obtained within the relaxation time approximation, as this allows for a phonon-resolved analysis, which is beneficial for benchmarking purposes. The BTE expression for the thermal conductivity within the RTA is
κ R T A = 1 N q V c λ C λ v λ v λ   τ λ .
Here, κ R T A is the thermal conductivity tensor, N q is the number of q-points in the sampling of reciprocal space, V c is the volume of the unit cell, C λ is the heat capacity associated with the phonon mode λ , v λ is the group velocity vector of that mode, and τ λ is its lifetime. The composite index, λ , contains the wave vector q and the band index n of each phonon, and the sum is formed over all phonon modes, λ . The mode heat capacities and the group velocities are harmonic phonon properties that can be obtained from the phonon bands, for which we already know that the MTPphonon results perfectly match the DFT data. Phonon lifetimes, however, require the calculation of at least third-order force constants, whose simulation is described in the Section 3.
Figure 8 compares the DFT- and MTPphonon-calculated phonon band structures of PE with the sizes of symbols scaled linearly and their coloring logarithmically with the respective phonon lifetimes. Again, the agreement between the two methods is excellent. Interestingly, for both methodologies, one observes an unusual behavior, as phonon lifetimes do not decrease steadily with frequency (as observed usually), but become particularly large for longitudinal acoustic phonons between 12 THz and 15 THz. Following the arguments of Wang et al. [73], who observed this unusual behavior for a single PE chain, it can be explained in the following way: longitudinal acoustic (LA) phonon modes between 5.5 THz and 11 THz have only a low contribution to the thermal conductivity, as the corresponding phonon lifetimes are small. This is the consequence of efficient Umklapp scattering processes involving twisting and transverse modes, e.g., LA → TWA + TA1 (LA = longitudinal acoustic, TWA = twisting acoustic, TA1 = transverse acoustic) [73]. Energy conservation, however, forbids scattering processes of the form LA → TWA + TA1 for phonons with frequencies above ~12 THz, as then, the LA phonons have more than twice the maximum energy of the TA and TWA modes. As a consequence, when only considering three-phonon processes a steep increase in phonon lifetime between 12 THz and 15 THz is obtained. The phonons of a single PE chain deviate in some respects from that of the PE crystal studied here. For example, there are no optical modes below 20 THz [73]. Moreover, for a single PE chain, the twisting mode has zero frequency at the Γ-point, and is, thus, an acoustic mode, whereas, in a PE crystal, the twisting mode has a nonzero frequency at the Γ-point and is, thus, an optical mode. Despite these differences between the single PE chain and the PE crystal, the core argument of Wang et al. concerning energy conservation can also be applied in the case of PE crystals.
With all ingredients in hand, the anisotropic thermal conductivity of PE can be calculated using Equation (1). The results for the full BTE and for the BTE in the RTA are listed in Table 4. This table also compares the results for the “best” MTPphonon with the mean values for all five parametrized versions of MTPphonon. A comparison to the results obtained with the “best” MTPMD is prevented by the above-mentioned imaginary frequency phonons calculated with this MTP. The overall agreement between the “best” MTPphonon and the DFT result is again excellent, especially in the chain direction. The observation that κzz differs by only 1 Wm−1K−1 must be considered as a fortunate coincidence. Still, the values for the arithmetic mean of the thermal conductivities for all five parametrized variants of MTPphonon also agree very well with the DFT results. Interestingly, as observed for the lattice constants, the standard deviations of the thermal conductivities calculated with the five MTPs are again similar to the deviations of the MTP means from the DFT results.
To provide more quantitative analyses of the excellent agreement between the DFT and MTPphonon calculated thermal conductivities (beyond the phonon lifetimes indicated in the phonon band structures in Figure 8), Figure 9 compares the contributions of phonons in specific frequency ranges. A complication in this context is that the apparent deviations, to some extent, depend on the way the frequency ranges are chosen: too large frequency ranges allow for an efficient cancellation of positive and negative deviations, while too small frequency ranges suffer from the finite sampling of reciprocal space in conjunction with the possibility that minor shifts of phonon frequencies can change the frequency ranges that specific phonons are associated with. Nevertheless, Figure 9 clearly demonstrates that the overall trends in the contributions of phonons in specific frequency ranges agree rather well between MTPphonon and DFT. A general trend is that, for phonons below ~10 THz, the MTPphonon-calculated contributions are somewhat larger, while the opposite is the case for higher-lying phonons. In view of the near perfect agreement of the phonon bands in Figure 7a, we attribute this to minor quantitative deviations in the calculated phonon lifetimes. Phonons above the band gap starting at ~16 THz do not play a role in the thermal conductivity.
The values of the thermal conductivity reported here compare reasonably well to the results described in the literature for applying very similar DFT-based theoretical approaches [73,74]. This is sufficient for the purpose of benchmarking the performance of the MTPs, but the used approach for obtaining κ suffers from several conceptual shortcomings. For one, phonon tunneling effects [77] are neglected (both for the DFT and the MTPphonon calculations). These primarily affect small thermal conductivities and, thus, in the present context, only the values for transport perpendicular to the polymer chains. Moreover, thermal expansion effects are not considered, i.e., for calculating thermal conductivities at room temperature, molecular-dynamics-calculated unit cells could be used (see also Section 2.7). Furthermore, effects like frequency renormalization and especially higher-order phonon scattering processes are neglected, which can distinctly reduce the phonon lifetimes above 12 THz due to a substantial increase in the allowed scattering processes. Also, the application of the relaxation time approximation within the BTE affects the results. All these aspects can be relevant for the final values of the thermal conductivity obtained from phonon properties, but they do not directly influence the benchmarking of the MTPs. Another aspect, especially when high-frequency phonons with appreciable lifetimes exist, is that molecular dynamics simulations would suffer from the classical statistics typically assumed for phonon occupations. This affects NEMD and AEMD simulations, as well as the phonon lifetimes obtained in Fourier transforms of equilibrium molecular dynamics runs. A further complication in NEMD and AEMD studies is that the mathematical form of the finite-size extrapolation is not uniquely defined. For systems with thermal conductivities as large as those obtained for PE and PT, this poses a problem, even when extending the simulations to supercells containing more than 100,000 atoms and considering several million time steps. In passing, we note that we performed such simulations with the MTPs discussed here on a contemporary supercomputer. A detailed discussion of all these aspects becomes rather involved and, thus, would seriously distract from the scope of the present paper. Therefore, we will address them in a forthcoming manuscript. In short, there it will be shown that, when considering all aforementioned effects, consistent values of the thermal conductivities for BTE and MD-based simulations are obtained with deviations in the range of a few percent for the best simulations of a specific type.

2.7. Studying the Thermal Expansion of PE

To calculate thermal expansion, MD simulations in an NPT ensemble at specific temperatures are performed, as described in Section 3 (with convergence tests provided in Section S29.8). As such simulations are not feasible over sufficiently long times when employing converged ab initio methods, we will benchmark the performance of MTPs only against temperature-dependent scattering experiments, which, to the best of our knowledge, only exist for PE. Thus, we will (again) focus on this material in the present section. Stability issues for MD at elevated temperatures prevent simulations with MTPphonon at relevant temperatures (see Section 2.2). Therefore, the focus here is on simulations employing MTPMD. Our results are compared to experiments by Davis et al. [78], Shen et al. [55], and Takahashi [53] in Figure 10. To better assess the following comparison, the peculiarities of the experimentally determined thermal expansion of PE are discussed first. The thermal expansion of PE is highly anisotropic, due to the dominance of vdW bonds perpendicular to the chain axis (a and b vectors) and covalent bonds along the chain (c vector). Lattice vector a (c.f., Figure 1a) significantly expands with temperature by around 0.3 Å from 100 K to 300 K [78], such that its evolution dominates the change in the unit cell volume. The expansion of b is about an order of magnitude smaller than that of a. In stark contrast to a and b, PE’s lattice vector c (along the chain) slightly shrinks with an increasing temperature. This means that, in the experiment by Davis et al. [78], PE has a negative expansion coefficient in that direction with the overall change between 100 K and 300 K being as small as 0.006 Å. In the experiment by Shen et al., the change in the length of c between 10 K and 297 K could not even be resolved [55]. Each of the five variants of MTPMDs that we parametrize captures these trends (see Figure S18). That is, they all show positive expansions of a and b that fit the experiments well, and a small negative expansion of c.
In Figure 10, our calculations are reported in two ways: for one, we ran the MD simulations for each of the five MTPMD and took the mean and standard deviation for the lattice constants. These are shown as blue data points with error bars. Additionally, the violet datapoints represent the results for the “best” variant of MTPMD. The calculated thermal expansion of lattice vectors a and b (vdW-bonding directions) agree well with experiments. For c, the agreement appears to be worse, with a seemingly larger offset in the absolute value. This is, however, mostly a consequence of the very small magnitudes of the changes in c. This means that, also in the c-direction, the MTPMD correctly predicts that there is virtually no change with temperature. For example, the calculated 100 K lattice constant |c| differs by only 0.003 Å (equaling 0.1%) compared to the “32 days” sample from Davis et al. [79]. For the sake of comparison, in Figure S20, a plot of the thermal expansion calculated with the AIREBO potential is contained. Not unexpectedly, that potential fares much worse than the MTPs: it yields lattice constants strongly deviating from experiments and predicts a wrong trend for the evolution of the lattice constant c with temperature. Finally, it should be noted that the agreement between theory and the experiments for the thermal expansion of PE is better than what we observed for metal–organic frameworks, potentially due to the ionic character and porous nature of these systems [29]. One also must not forget that, even when using ab initio methods, accurately simulating thermal expansion is far from trivial and prone to serious errors, with even minor changes in the used DFT functionals potentially changing the signs of the expansion coefficients [79]. This suggests that an unknown, but potentially significant fraction of the remaining deviation between theory and experiments in Figure 10 is due to shortcomings of the DFT methodology used for the learning of the MTPs, rather than due to the nature of the MTPs and the employed parametrization process.

2.8. Energy, Force, and Stress in Molecular Dynamics

As a final step, to test the accuracy of the MTPs in a context intrinsically relevant for MD simulations, we benchmark the performance of the “best” variant of MTPMD in predicting the total energies, forces, and stresses. This is conducted for the DFT-calculated structures generated during the validation active learning run described in Section 2.2. Notably, none of these data have been used in the parametrization process of the MTP. Figure 11 contains bar plots of the differences between the MTP- and DFT-calculated energies, forces, and stresses. The errors of the MTP are generally very low: for total energies significantly below 1%, for forces mostly below 5%, and for stresses also usually below 5%, albeit with a second maximum in the histogram at a deviation around 10% (as shown in Figure 11a,c,e). To put the MTP accuracy into perspective, we compared the deviations between the MTP and DFT results to those between the AIREBO and DFT data in the right panels of Figure 11b,d,f. There, much larger scales for plotting ranges of energy, force, and stress errors are necessary to even display the AIREBO data. These data show that the MTPMD vastly outperforms the AIREBO potential: for the energies, the RMSD of MTPMD is more than two orders of magnitude (factor of 308) smaller than for the AIREBO potential (0.072 meV/atom vs. 22 meV/atom), for the forces, the difference amounts to nearly two orders of magnitude (factor of 64) with RMSDMD values of 0.011 eV/Å vs. 0.72 eV/Å, and for the stress errors, the deviation is similar (factor of 32) with 0.26 kbar vs. 8.3 kbar.

2.9. Speed Gain Due to Use of Moment Tensor Potentials

A comprehensive assessment of the speedup by using MTPs instead of employing DFT is futile, as it depends on a variety of aspects, including the nature of the studied system, the number of atoms in the unit cell, the employed computer architecture, and the chosen degree of parallelization—e.g., the number of cores used in the simulations, the details of the DFT approach (including aspects like the chosen basis sets, the k-point sampling, convergence criteria, etc.). Also, the codes used for the specific calculations play an important role. Still, to qualitatively assess the overall situation, a specific example shall be discussed here (with details on the used computer architecture and on the MTP, as well as on the DFT simulations provided in Section S27 and in the Section 3): to calculate the phonon band structure of PE, the converged supercell contains 432 atoms. Starting from an already parametrized MTP, the MTPphonon calculation to obtain the bands (shown in Figure 7) was roughly 105 times faster than the DFT simulation (see Section S27 on how we arrived at this estimate). While the computational complexity of DFT calculations scales cubically with the number of electrons in a naive DFT implementation, MTP simulation times scale linearly with the number of atoms. Therefore, MTPs have an even larger advantage for larger cell sizes, e.g., for converging NEMD and AEMD simulations to obtain thermal conductivities. The exact number of atoms that need to be considered in such calculations strongly depends on the magnitude of the thermal conductivity, but it typically reaches several ten thousand atoms; for materials with high thermal conductivities like PE, in the chain direction, this is potentially even more. For such a case, one can estimate a speedup on the order of nearly 1010 (needless to say that the exact value cannot be determined, as converged DFT simulations are entirely infeasible for such a number of atoms).
What still needs to be accounted for is the computational effort for parametrizing the MTP, which has to be performed separately for every considered material. For generating the training structures for MTPphonon, we used the same supercell as that for the phonon calculations. For PE, this resulted in 129 DFT-based force calculations in the active learning run between 15 K and 100 K. Training one MTPphonon using these reference data requires a similar amount of resources due to the very high level of the MTPphonon (see Section S27). Thus, in total, the parametrization of a single variant of MTPphonon costs the equivalent of ~250 DFT calculations (here, with 432 atoms). For parametrizing a single MTPMD, 550 ab initio MD steps were performed, but since the effort of parametrizing the lower-level MTPMD is much lower, the overall effort is ~2.5 times as large as that for the MTPphonon. This means that, for obtaining quantities that are relatively easy to calculate, parametrizing an MTP is not necessarily advisable. For example, due to the high symmetry of the material and the small primitive unit cell, for calculating the phonon band structure of PE, only 12 DFT calculations are needed. This means that, for PE, calculating the phonon band structure with an MTPphonon is approximately twenty times as expensive as calculating the phonon band structure with DFT in the first place due to the MTP parametrization. For more complex materials requiring a much larger number of DFT calculations for obtaining phonon band structures (but a similar number of DFT reference calculations for the MTP parametrization), this gap diminishes or even closes. This, we have, for example, observed for a variety of metal–organic frameworks.
MTPs become much more efficient than DFT when phonon lifetimes need to be calculated, because, in that case, even for PE, 8000 DFT calculations are necessary. This becomes more of an issue for larger systems, where the number of necessary DFT calculations becomes insurmountable (around 50,000 and 360,000 for PT and P3HT, respectively). Even more significantly, performing NEMD calculations over 5 ns with 0.5 fs time steps requires calculating forces ten million times. This typically involves a much larger unit cell than what is needed for converging phonon properties. These numbers illustrate the main benefit of the presented approach of parametrizing MTPs for polymer materials, namely, that it opens up the possibility for performing new types of calculations, for which DFT would be far too expensive and off-the-shelf, transferrable force fields would be far too inaccurate.

3. Methods

3.1. Details of the Applied Computational Approach

DFT calculations were performed with VASP version 6.3.0 [44]. Its primary task was to generate training data for parametrizing the MTPs. VASP was employed to calculate elastic constants and (in conjunction with phonopy and phono3py [9,80,81]) to calculate phonon band structures and the thermal conductivity of PE (employing the BTE, as described in the Section 2.6 and Section 3.2). The PBE functional [43], complemented by Grimme’s D3 correction [18,19] with Becke–Johnson damping [20], was employed, as this combination is well established for accurately describing the phonons of organic semiconductors [17,21,82,83]. Following convergence tests for the phonon band-structures shown in Figure S21, the energy cutoff determining the size of the plane-wave basis set was set to 900 eV (PE) and to 700 eV (PT and P3HT). When relaxing cells, the energy cutoff was increased to 1350 eV to mitigate Pulay stress. Geometry relaxations (with a convergence criterion for maximum forces of 0.001 eV/Å) were performed in two steps: first, all degrees of freedom, including the unit cell parameters, were relaxed with the 1350 eV cutoff. Subsequently, the unit cell parameters were kept fixed, and only the atomic positions were relaxed with an energy cutoff of 900 eV (PE) or 700 eV (PT and P3HT). The second step was necessary to relax the atomic positions with the same energy cutoff that was later used in the phonon calculations.
K-meshes for describing the electronic wavefunctions were set such that the phonon band structures were converged. We found that phonon calculations with supercells (see below) converged with 1 × 1 × 1 (PE and P3HT) and 1 × 1 × 2 (PT) k-meshes. This defined the k-meshes used for the primitive unit cells as 2 × 3 × 6 (PE), 2 × 3 × 4 (PT), and 2 × 1 × 2 (P3HT). The convergence tests of the k-mesh can be found in Figure S22. The break condition for the electronic self-consistency loop was set to 10−8 eV. Partial occupancies of orbitals were treated with a Gaussian smearing with a width of 0.05 eV.
In the VASP active learning runs, mostly the same settings as those described above were used. Also, in active learning runs, unit cells are distorted, again causing Pulay stress. To mitigate this, an increase in the plane-wave cutoff of at least 30% compared to fixed volume calculations is recommended on the VASP wiki [51]. As large volume fluctuations indeed occurred for P3HT, we increased the energy cutoff in the active learning from 700 eV to 900 eV. As no large volume fluctuations occurred for PE and PT, for these systems, the cutoffs were kept at 900 eV and 700 eV, respectively. As an additional measure to mitigate Pulay stress, the active learning runs were split into smaller sections to reinitialize the plane wave basis sets. For P3HT, the runs were performed from 15 K to 300 K, 300 K to 400 K, and 400 K to 500 K. For PE, the learning runs were split into sections from 15 K to 30 K, 30 K to 150 K, 150 K to 300 K, 300 K to 400 K, and 400 K to 500 K (see discussion in Section S7). For PT, the run was not split into multiple segments, as the change in the unit cell size was particularly small for that polymer. Supercells for the training data were chosen such that they extended approximately 15 Å long in each direction. This choice generally also yields converged phonon band structures, as detailed below and shown in Figure S26. Correspondingly, a 2 × 3 × 6 supercell was chosen for PE, a 2 × 3 × 2 supercell for PT, and a 2 × 1 × 2 supercell for P3HT. For the active learning MD runs, a Langevin thermostat with a friction of 10 ps−1 for all atoms and the unit cell was used. VASP active learning runs are generally performed in an NPT ensemble (except for Section S5, where we describe the influence of using an NVT ensemble). For the active learning runs, the fictitious mass of the lattice degrees-of-freedom was set to 1000 amu. The time step was set to 0.5 fs. For training the MTPs geared towards phonon calculations (MTPphonon), we performed an NPT run from 15 to 100 K for 12,750 time steps (i.e., 150 steps/K). This yielded 129 DFT-calculated training structures for PE, 173 training structures for PT, and 139 structures for P3HT. For the MTPMD training data, the run was continued to 500 K, yielding a total of 72,750 steps (i.e., again with 150 steps/K). This resulted in 550, 490, and 549 DFT-calculated training structures for PE, PT, and P3HT, respectively.
To assess the accuracy of the MTPs, we created a validation set that consisted of independent DFT reference calculations from an MD run at a fixed temperature of 300 K starting from scratch (i.e., disregarding any data from the above-mentioned training). This run lasted for 10,000 time steps in an NPT ensemble for PE and P3HT. For PT, because the structure became unstable during the run (i.e., atoms detached from the polymer chains), we had to enforce more DFT steps to make it stable. This was achieved by setting the tag ML_CX to −0.15 (see VASP manual [84]). This tag controlled the threshold that determined whether a DFT calculation was performed and helped to avoid too few training data for runs at constant temperature [51]. As a result of this procedure, we ended up with more DFT reference structures for PT as compared to the other two materials. Thus, the validation sets consisted of 136 structures for PE, 296 for PT, and 182 for P3HT. Some initial configurations are associated with a lower temperature, because it takes some time for the thermostat to reach the desired 300 K. Considering the example of PE, it took around 250 time steps to reach 300 K. We found that this was irrelevant for the calculated RMSDMD, whether we cut out these poorly equilibrated structures or not. Therefore, we kept them in the validation set.
The AIREBO potential [71] was used as implemented in LAMMPS, with its functional form given in Section S23. The scaling factor for the Lennard–Jones term, which defines the cutoff, was set to 3.

3.2. Modeling Physical Observables

Harmonic phonon properties were obtained using the phonopy package [9,85]. It uses a finite differences scheme with a displacement distance of 0.01 Å. This requires the calculation of forces on all atoms upon displacing individual atoms. The number of displacements that eventually needs to be obtained can be significantly reduced by considering symmetries. For example, in the case of PE (with 12 atoms in the primitive unit cell), only 12 calculations are necessary. Notably, in such calculations, sufficiently large supercells must be used with details on the supercell convergence contained in Figure S26. For our materials, we typically found well-converged phonon band structures for supercells, with an extent of about 15 Å in each direction. An exception was PT, for which a supercell of 2 × 3 × 2 would be around 15 Å large in each direction. However, for this supercell size, the phonon band structure had a small dip in imaginary values along the Γ-X path close to Γ. To mitigate this dip, a 2 × 3 × 4 supercell was used. Correspondingly, 2 × 3 × 6, 2 × 3 × 4, and 2 × 1 × 2 supercells were used for the phonon calculations of PE, PT, and P3HT. The same supercells were used in the DFT and MTP simulations. As described above, to quantitatively measure the agreement between MTP and DFT, we defined the root mean square deviation (RMSD) between the phonon frequencies of the MTP and DFT, RMSDphonon. We calculated the RMSDphonon over the whole Brillouin zone, which we sampled with dense q-point meshes of 30 × 44 × 83, 30 × 42 × 30, and 10 × 5 × 9 q-points for PE, PT, and P3HT, respectively.
For calculating the anharmonic phonon properties, we employed the phono3py code [9,80], which considers third-order force constants and, consequently, accounts for three-phonon scattering processes. Phono3py was used to calculate the phonon lifetimes and the thermal conductivity via the Boltzmann transport equation (BTE), either directly (keyword “lbte”) or within the relaxation time approximation. Displacements of 0.03 Å were used for calculating these third-order force constants. This is the default in phono3py and was checked for convergence (see Section S29.6). The supercell for the third-order force constant calculation, which was performed only for PE, was a 2 × 2 × 3 cell for both the DFT and the MTP calculations. Convergence tests for choosing that supercell size are provided in Figure S27. Brillouin zone integration in phono3py is performed with the default tetrahedron method [86], which has the advantage of being a parameter-free approach. Brillouin zones are sampled with a 10 × 15 × 160 mesh for the RTA calculation. When employing the full BTE, a 10 × 15 × 60 mesh is used for κxx and κyy, while a 4 × 6 × 160 mesh is used for κzz. These different mesh sizes for the respective elements of κ were used to reduce computational and memory costs, as a 10 × 15 × 160 mesh for the full BTE calculation would be infeasible due to the excessive demand for memory. Convergence tests for these meshes are provided in Figure S28. Phonopy and phono3py require forces, agnostic to whether they are calculated with VASP or MLIP, therefore, the same phonopy and phono3py settings can be used for both force calculators.
To obtain the elastic stiffness tensors, a finite differences approach was employed [87]. For the DFT calculations, we used an in-house developed code, in which the cell was strained by 1%. This strain distance is converged, as elastic constants do not change significantly whether strains of 1% or 0.1% are used, and change only slightly for strains of 0.01%, as is shown in Section S29.3. The contribution from relaxing the ions in the strained cells is approximated from second-order force constants, and the elastic tensor is obtained from the strain–stress relationship [62,87]. This approximation is performed to save computational resources, and we refer to this approach as the “clamped ion method” [88]. This method is identical to the VASP-internal routine (setting IBRION = 6 and ISIF = 3), with the important difference that, in the VASP-internal approach, the POTIM parameter controls both the atomic displacement distance to calculate second-order force constants as well as the strain distance, whereas, in our in-house developed code, these two settings are disentangled. This allowed us to keep the atomic displacement at 0.01 Å for the second-order force constants (calculated via phonopy), while converging the strain distance independently. In fact, we used the force constants for which we already obtained converged phonon band structures. The energy cutoffs were also kept identical to those used in the phonon band structure calculation. A calculation of P3HT with an energy cutoff of 900 eV is given in Table S15. It shows only small deviations from the calculation with an energy cutoff of 700 eV, thereby suggesting that the energy cutoff is converged for calculating the elastic tensor. For the elastic constants calculations with MTPs in the main paper, we employed the same method (“clamped ion method”) as that with DFT for consistency reasons. However, since the MTPs are orders of magnitude faster than DFT, it is also substantially cheaper to calculate elastic constants with a less approximative method. In this method, after straining the cell by 0.1%, the ionic positions were relaxed and the stress tensor for the strain–stress relationship was directly calculated with the MTPs. We refer to this less approximative method as the “relaxed ion method”, and the results are shown in Section S15.
Thermal expansion was calculated in conventional MD simulations employing the MTPs. MD was performed in LAMMPS [49] (Version 2nd July 2021) in combination with the LAMMPS–MLIP interface (Version 2) [26]. A Nose–Hoover thermostat with a damping parameter of 0.1 ps and a Nose–Hoover barostat with damping parameter of 1 ps were used. A time step of 0.5 fs was chosen, which is consistent with previous works, where time steps between 0.25 fs and 2 fs were employed for (non-equilibrium) molecular dynamics simulations of PE [14,89,90]. For calculating thermal expansion, an orthorhombic 4 × 6 × 12 supercell was used, which has an extent of approximately 30 Å in all dimensions. An NPT run was performed for 800 ps at each set temperature, during which the lattice parameters were tracked (see Figure S29 for a convergence test). Ignoring the first 10 ps to account for equilibration effects, the lattice parameters were then averaged over the remaining time of the run.

4. Conclusions

Machine-learned potentials are a powerful tool in materials’ modeling, as they allow simulating structural and dynamic properties with close to DFT accuracy at computational costs that are reduced by many orders of magnitude. This is illustrated here for moment tensor potentials (MTP), which are parametrized for specific materials, with a focus on phonon-related properties and molecular dynamics (MD) simulations. In particular, we discuss how to optimally parametrize MTPs and how to benchmark their performance.
For the parametrization, it turned out to be useful to distinguish between two distinct use cases: as the first use case, MTPs were applied for finite difference calculations, where all atoms remained close to their equilibrium positions. These included calculations of phonon band structures and the thermal conductivity via the Boltzmann transport equation. The second use case was molecular dynamics calculations at elevated temperatures (e.g., room temperature), which can be employed to, for example, simulate thermal expansion. In such simulations, significantly larger atomic displacements occur. This calls for establishing separate protocols for these two scenarios. In that regard, a crucial aspect is the suitable selection of training data, which, in our case, were generated during active learning molecular dynamics simulations. The magnitudes of atomic displacements covered by the training data were determined by the considered temperature ranges in the active learning runs. In this context, we found it to be advantageous to parametrize the MTPs for the phonon calculations using reference structures generated at low temperatures (e.g., between 15 K and 100 K), while a force field for MD calculations had to be parametrized against reference data, whose generation involved higher-temperature reference structures (generated, e.g., between 15 K and 500 K). In particular, we found that MTPs trained on low-temperature data were twice as accurate in phonon calculations as their counterparts trained also on structures generated at higher temperatures. The opposite was observed for MD simulations, where MTPs containing only low-temperature data in their parametrization even turned out to be unstable.
Further strategies for improving the accuracy are increasing the level of the MTP and splitting chemically inequivalent atoms into multiple atom types. Both of these strategies increase the number of parameters of the MTPs, which, in the case of the level, comes at the cost of considerably increased computational costs. For the phonon calculation, computation time is not such a serious concern, as the size of the considered supercells and the number of force calculations are not excessively large. Therefore, a higher level can be used than for MTPs that are meant for molecular dynamics calculations. The situation is different when increasing the number of atom types. This only increases the time needed for the parametrization, which leads us to suggest to, by default, treat each atom in a chemically different environment as a separate species. These considerations yield two optimum styles of MTPs for the two described use cases: MTPphonon, which is a level 26 MTP trained against active learning data generated between 10 K and 100 K, and MTPMD, as a level 22 MTP trained on data generated between 10 K and 500 K. The two types of MTPs were benchmarked in terms of their ability to describe unit cells, phonon band structures, and elastic constants for polyethylene (PE), polythiophene (PT), and poly-3-hexylthiophene (P3HT). For PE, also the thermal conductivity (via the Boltzmann transport equation), the thermal expansion (via MD), and forces, energies, and strains during an MD run were assessed. All these tests (i) testify to the exceptional accuracy of system-specific MTPs (relative to DFT data and experiments) and (ii) show that the MTP variants parametrized for the specific use cases display a distinctly improved performance.
All this was achieved with the MTP simulations being between five and (an estimated) ten orders of magnitude faster than the DFT calculations using the methodology employed for their parametrization. For calculating relatively simple quantities (like structural parameters and phonon band structures) of relatively simple materials (like PE), this increase in speed is, however, more than offset by the computational efforts needed to parametrize the MTPs. This is clearly not the case when calculating computationally more demanding quantities, like phonon lifetimes or MD trajectories (for obtaining thermal expansion coefficients or thermal conductivities).
Many of these calculations only become possible by machine-learned potentials, as ab initio methods like DFT would be far too expensive, and off-the-shelf force fields would be too inaccurate. Thus, we are confident that the current manuscript will have a profound impact on the future simulation of the structural and dynamic properties of polymers and will, for example, lay the basis for the quantitatively reliable simulation of their heat transport properties employing molecular dynamics and/or lattice dynamics simulations.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/molecules29163724/s1. References [91,92,93,94,95,96,97,98,99,100] are cited in the Supplementary Materials.

Author Contributions

Conceptualization, E.Z. and L.R.; methodology, L.R. and S.W.; software, L.R., S.W. and L.L.; validation, L.R.; formal analysis, L.R. and L.L.; investigation, L.R.; resources, E.Z.; data curation, L.R.; writing—original draft preparation, L.R.; writing—review and editing, E.Z., L.R., L.L. and S.W.; visualization, L.R.; supervision, E.Z.; project administration, E.Z.; funding acquisition, E.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Austrian Science Fund (FWF) [primarily Grant-DOI: 10.55776/P33903 and in part also 10.55776/P36129]. For the purpose of open access, the authors have applied a CC-BY public copyright license to any author accepted manuscript version arising from this submission. We also acknowledge the Graz University of Technology for support through the Lead Project Porous Materials @ Work for Sustainability (LP-03). Computational results have been obtained using the Vienna Scientific Cluster, VSC-5. Open Access Funding by the Graz University of Technology.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All relevant data can be downloaded from the TU Graz data repository: https://doi.org/10.3217/7s1ce-ss195.

Acknowledgments

We acknowledge the work of Tomas Kamencek, who with his thorough studies on how to best perform phonon calculations, laid the groundwork for this paper. We thank Enrico Di Lucente for stimulating discussions.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. He, L.; Ye, Z.; Zeng, J.; Yang, X.; Li, D.; Yang, X.; Chen, Y.; Huang, Y. Enhancement in Electrical and Thermal Properties of LDPE with Al2O3 and H-BN as Nanofiller. Materials 2022, 15, 2844. [Google Scholar] [CrossRef]
  2. Yang, C.; Navarro, M.E.; Zhao, B.; Leng, G.; Xu, G.; Wang, L.; Jin, Y.; Ding, Y. Thermal Conductivity Enhancement of Recycled High Density Polyethylene as a Storage Media for Latent Heat Thermal Energy Storage. Sol. Energy Mater. Sol. Cells 2016, 152, 103–110. [Google Scholar] [CrossRef]
  3. Liu, J.; Xu, Z.; Cheng, Z.; Xu, S.; Wang, X. Thermal Conductivity of Ultrahigh Molecular Weight Polyethylene Crystal: Defect Effect Uncovered by 0 K Limit Phonon Diffusion. ACS Appl. Mater. Interfaces 2015, 7, 27279–27288. [Google Scholar] [CrossRef]
  4. Kim, T.; Drakopoulos, S.X.; Ronca, S.; Minnich, A.J. Origin of High Thermal Conductivity in Disentangled Ultra-High Molecular Weight Polyethylene Films: Ballistic Phonons within Enlarged Crystals. Nat. Commun. 2022, 13, 2452. [Google Scholar] [CrossRef] [PubMed]
  5. Xu, Y.; Kraemer, D.; Song, B.; Jiang, Z.; Zhou, J.; Loomis, J.; Wang, J.; Li, M.; Ghasemi, H.; Huang, X.; et al. Nanostructured Polymer Films with Metal-like Thermal Conductivity. Nat. Commun. 2019, 10, 1771. [Google Scholar] [CrossRef]
  6. Ronca, S.; Igarashi, T.; Forte, G.; Rastogi, S. Metallic-like Thermal Conductivity in a Lightweight Insulator: Solid-State Processed Ultra High Molecular Weight Polyethylene Tapes and Films. Polymer 2017, 123, 203–210. [Google Scholar] [CrossRef]
  7. Shrestha, R.; Li, P.; Chatterjee, B.; Zheng, T.; Wu, X.; Liu, Z.; Luo, T.; Choi, S.; Hippalgaonkar, K.; De Boer, M.P.; et al. Crystalline Polymer Nanofibers with Ultra-High Strength and Thermal Conductivity. Nat. Commun. 2018, 9, 1664. [Google Scholar] [CrossRef] [PubMed]
  8. Shen, S.; Henry, A.; Tong, J.; Zheng, R.; Chen, G. Polyethylene Nanofibres with Very High Thermal Conductivities. Nat. Nanotechnol. 2010, 5, 251–255. [Google Scholar] [CrossRef] [PubMed]
  9. Togo, A. First-Principles Phonon Calculations with Phonopy and Phono3py. J. Phys. Soc. Jpn. 2023, 92, 012001. [Google Scholar] [CrossRef]
  10. Müller-Plathe, F. A Simple Nonequilibrium Molecular Dynamics Method for Calculating the Thermal Conductivity. J. Chem. Phys. 1997, 106, 6082–6085. [Google Scholar] [CrossRef]
  11. Green, M.S. Markoff Random Processes and the Statistical Mechanics of Time-Dependent Phenomena. II. Irreversible Processes in Fluids. J. Chem. Phys. 1954, 22, 398–413. [Google Scholar] [CrossRef]
  12. Kubo, R. Statistical-Mechanical Theory of Irreversible Processes. I. General Theory and Simple Applications to Magnetic and Conduction Problems. J. Phys. Soc. Jpn. 1957, 12, 570–586. [Google Scholar] [CrossRef]
  13. Lampin, E.; Palla, P.L.; Francioso, P.A.; Cleri, F. Thermal Conductivity from Approach-to-Equilibrium Molecular Dynamics. J. Appl. Phys. 2013, 114, 33525. [Google Scholar] [CrossRef]
  14. Zhang, T.; Wu, X.; Luo, T. Polymer Nanofibers with Outstanding Thermal Conductivity and Thermal Stability: Fundamental Linkage between Molecular Characteristics and Macroscopic Thermal Properties. J. Phys. Chem. C 2014, 118, 21148–21159. [Google Scholar] [CrossRef]
  15. Henry, A.; Chen, G.; Plimpton, S.J.; Thompson, A. 1D-to-3D Transition of Phonon Heat Conduction in Polyethylene Using Molecular Dynamics Simulations. Phys. Rev. B Condens. Matter Mater. Phys. 2010, 82, 144308. [Google Scholar] [CrossRef]
  16. Ni, B.; Watanabe, T.; Phillpot, S.R. Thermal Transport in Polyethylene and at Polyethylene-Diamond Interfaces Investigated Using Molecular Dynamics Simulation. J. Phys. Condens. Matter 2009, 21, 084219. [Google Scholar] [CrossRef] [PubMed]
  17. Kamencek, T.; Wieser, S.; Kojima, H.; Bedoya-Martínez, N.; Dürholt, J.P.; Schmid, R.; Zojer, E. Evaluating Computational Shortcuts in Supercell-Based Phonon Calculations of Molecular Crystals: The Instructive Case of Naphthalene. J. Chem. Theory Comput. 2020, 16, 2716–2735. [Google Scholar] [CrossRef] [PubMed]
  18. Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the Damping Function in Dispersion Corrected Density Functional Theory. J. Comput. Chem. 2010, 32, 1456. [Google Scholar] [CrossRef] [PubMed]
  19. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104–154119. [Google Scholar] [CrossRef]
  20. Johnson, E.R.; Becke, A.D. A Post-Hartree-Fock Model of Intermolecular Interactions: Inclusion of Higher-Order Corrections. J. Chem. Phys. 2006, 124, 174104. [Google Scholar] [CrossRef]
  21. Bedoya-Martínez, N.; Giunchi, A.; Salzillo, T.; Venuti, E.; Della Valle, R.G.; Zojer, E. Toward a Reliable Description of the Lattice Vibrations in Organic Molecular Crystals: The Impact of van Der Waals Interactions. J. Chem. Theory Comput. 2018, 14, 4380–4390. [Google Scholar] [CrossRef] [PubMed]
  22. Mortazavi, B.; Zhuang, X.; Rabczuk, T.; Shapeev, A.V. Atomistic Modeling of the Mechanical Properties: The Rise of Machine Learning Interatomic Potentials. Mater. Horiz. 2023, 10, 1956–1968. [Google Scholar] [CrossRef] [PubMed]
  23. Ouyang, Y.; Yu, C.; He, J.; Jiang, P.; Ren, W.; Chen, J. Accurate Description of High-Order Phonon Anharmonicity and Lattice Thermal Conductivity from Molecular Dynamics Simulations with Machine Learning Potential. Phys. Rev. B 2022, 105, 115202. [Google Scholar] [CrossRef]
  24. Mortazavi, B.; Podryabinkin, E.V.; Novikov, I.S.; Rabczuk, T.; Zhuang, X.; Shapeev, A.V. Accelerating First-Principles Estimation of Thermal Conductivity by Machine-Learning Interatomic Potentials: A MTP/ShengBTE Solution. Comput. Phys. Commun. 2021, 258, 107583. [Google Scholar] [CrossRef]
  25. Shapeev, A.V. Moment Tensor Potentials: A Class of Systematically Improvable Interatomic Potentials. Multiscale Model. Simul. 2016, 14, 1153–1173. [Google Scholar] [CrossRef]
  26. Novikov, I.S.; Gubaev, K.; Podryabinkin, E.V.; Shapeev, A.V. The MLIP Package: Moment Tensor Potentials with MPI and Active Learning. Mach. Learn. Sci. Technol. 2021, 2, 025002. [Google Scholar] [CrossRef]
  27. Mortazavi, B.; Rajabpour, A.; Zhuang, X.; Rabczuk, T.; Shapeev, A.V. Exploring Thermal Expansion of Carbon-Based Nanosheets by Machine-Learning Interatomic Potentials. Carbon 2022, 186, 501–508. [Google Scholar] [CrossRef]
  28. Mortazavi, B.; Novikov, I.S.; Shapeev, A.V. A Machine-Learning-Based Investigation on the Mechanical/Failure Response and Thermal Conductivity of Semiconducting BC2N Monolayers. Carbon 2022, 188, 431–441. [Google Scholar] [CrossRef]
  29. Wieser, S.; Zojer, E. Machine Learned Force-Fields for an Ab-Initio Quality Description of Metal-Organic Frameworks. NPJ Comput. Mater. 2024, 10, 18. [Google Scholar] [CrossRef]
  30. Strasser, N.; Wieser, S.; Zojer, E. Predicting Spin-Dependent Phonon Band Structures of HKUST-1 Using Density Functional Theory and Machine-Learned Interatomic Potentials. Int. J. Mol. Sci. 2024, 25, 3023. [Google Scholar] [CrossRef]
  31. Togo, A.; Seko, A. On-the-Fly Training of Polynomial Machine Learning Potentials in Computing Lattice Thermal Conductivity. J. Chem. Phys. 2024, 160, 211001. [Google Scholar] [CrossRef] [PubMed]
  32. Sirringhaus, H. Device Physics of Solution-Processed Organic Field-Effect Transistors. Adv. Mater. 2005, 17, 2411–2425. [Google Scholar] [CrossRef]
  33. Knobloch, A.; Manuelli, A.; Bernds, A.; Clemens, W. Fully Printed Integrated Circuits from Solution Processable Polymers. J. Appl. Phys. 2004, 96, 2286–2291. [Google Scholar] [CrossRef]
  34. Lim, J.A.; Liu, F.; Ferdous, S.; Muthukumar, M.; Briseno, A.L. Polymer Semiconductor Crystals. Mater. Today 2010, 13, 14–24. [Google Scholar] [CrossRef]
  35. Pacher, P.; Lex, A.; Proschek, V.; Etschmaier, H.; Tchernychova, E.; Sezen, M.; Scherf, U.; Grogger, W.; Trimmel, G.; Slugovc, C.; et al. Chemical Control of Local Doping in Organic Thin-Film Transistors: From Depletion to Enhancement. Adv. Mater. 2008, 20, 3143–3148. [Google Scholar] [CrossRef]
  36. Peacock, A. Handbook of Polyethylene; Taylor & Francis: London, UK, 2000. [Google Scholar]
  37. Huan, T.D.; Ramprasad, R. Polymer Structure Prediction from First Principles. J. Phys. Chem. Lett. 2020, 11, 5823–5829. [Google Scholar] [CrossRef] [PubMed]
  38. Corish, J.; Morton-Blake, D.A.; Veluri, K.; Bénière, F. Atomistic Simulations of the Structures of the Pristine and Doped Lattices of Polypyrrole and Polythiophene. J. Mol. Struct. THEOCHEM 1993, 283, 121–134. [Google Scholar] [CrossRef]
  39. Cheng, P.; Shulumba, N.; Minnich, A.J. Thermal Transport and Phonon Focusing in Complex Molecular Crystals: Ab Initio Study of Polythiophene. Phys. Rev. B 2019, 100, 94306. [Google Scholar] [CrossRef]
  40. Zhugayevych, A.; Mazaleva, O.; Naumov, A.; Tretiak, S. Lowest-Energy Crystalline Polymorphs of P3HT. J. Phys. Chem. C 2018, 122, 9141–9151. [Google Scholar] [CrossRef]
  41. Brückner, S.; Porzio, W. The Structure of Neutral Polythiophene. An Application of the Rietveld Method. Die Makromol. Chem. 1988, 189, 961–967. [Google Scholar] [CrossRef]
  42. Momma, K.; Izumi, F. VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Cryst. 2011, 44, 1272–1276. [Google Scholar] [CrossRef]
  43. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. [Google Scholar] [CrossRef] [PubMed]
  44. Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558–561. [Google Scholar] [CrossRef] [PubMed]
  45. Morrow, J.D.; Gardner, J.L.A.; Deringer, V.L. How to Validate Machine-Learned Interatomic Potentials. J. Chem. Phys. 2023, 158, 121501. [Google Scholar] [CrossRef] [PubMed]
  46. Choi, J.M.; Lee, K.; Kim, S.; Moon, M.; Jeong, W.; Han, S. Accelerated Computation of Lattice Thermal Conductivity Using Neural Network Interatomic Potentials. Comput. Mater. Sci. 2022, 211, 111472. [Google Scholar] [CrossRef]
  47. Jinnouchi, R.; Karsai, F.; Verdi, C.; Asahi, R.; Kresse, G. Descriptors Representing Two-and Three-Body Atomic Distributions and Their Effects on the Accuracy of Machine-Learned Inter-Atomic Potentials. J. Chem. Phys. 2020, 152, 234102. [Google Scholar] [CrossRef]
  48. Liu, P.; Wang, J.; Avargues, N.; Verdi, C.; Singraber, A.; Karsai, F.; Chen, X.-Q.; Kresse, G. Combining Machine Learning and Many-Body Calculations: Coverage-Dependent Adsorption of CO on Rh(111). Phys. Rev. Lett. 2022, 130, 78001. [Google Scholar] [CrossRef]
  49. Thompson, A.P.; Aktulga, H.M.; Berger, R.; Bolintineanu, D.S.; Brown, W.M.; Crozier, P.S.; in ’t Veld, P.J.; Kohlmeyer, A.; Moore, S.G.; Nguyen, T.D.; et al. LAMMPS—A Flexible Simulation Tool for Particle-Based Materials Modeling at the Atomic, Meso, and Continuum Scales. Comput. Phys. Commun. 2022, 271, 108171. [Google Scholar] [CrossRef]
  50. Verdi, C.; Karsai, F.; Liu, P.; Jinnouchi, R.; Kresse, G. Thermal Transport and Phase Transitions of Zirconia by On-the-Fly Machine-Learned Interatomic Potentials. NPJ Comput. Mater. 2021, 7, 156. [Google Scholar] [CrossRef]
  51. VASP Best Practices for Machine-Learned Force Fields. Available online: https://www.vasp.at/wiki/index.php/Best_practices_for_machine-learned_force_fields (accessed on 4 June 2024).
  52. Avitabile, G.; Napolitano, R.; Pirozzi, B.; Rouse, K.D.; Thomas, M.W.; Willis, B.T.M. Low Temperature Crystal Structure of Polyehtylene: Results from a Neutron Diffraction Study and from Potential Energy Calculations. J. Polym. Sci. Polym. Lett. Ed. 1975, 13, 351–355. [Google Scholar] [CrossRef]
  53. Takahashi, Y. Neutron Structure Analysis of Polyethylene-D4. Macromolecules 1998, 31, 3868–3871. [Google Scholar] [CrossRef]
  54. Mo, Z.; Lee, K.B.; Moon, Y.B.; Kobayashi, M.; Heeger, A.J.; Wudl, F. X-Ray Scattering from Polythiophene: Crystallinity and Crystallographic Structure. Macromolecules 1985, 18, 1972–1977. [Google Scholar] [CrossRef]
  55. Shen, M.; Hansen, W.N.; Romo, P.C. Thermal Expansion of the Polyethylene Unit Cell. J. Chem. Phys. 1969, 51, 425–430. [Google Scholar] [CrossRef]
  56. Xavier, N.F.; Da Silva, A.M.; Bauerfeldt, G.F. What Rules the Relative Stability of α-, β-, and γ-Glycine Polymorphs? Cryst. Growth Des. 2020, 20, 4695–4706. [Google Scholar] [CrossRef]
  57. Červinka, C.; Beran, G.J.O. Ab Initio Thermodynamic Properties and Their Uncertainties for Crystalline α-Methanol. Phys. Chem. Chem. Phys. 2017, 19, 29940–29953. [Google Scholar] [CrossRef] [PubMed]
  58. Nakamae, K.; Nishino, T.; Ohkubo, H. Elastic Modulus of Crystalline Regions of Polyethylene with Different Microstructures: Experimental Proof of Homogeneous Stress Distribution. J. Macromol. Sci. Part B 1991, 30, 1–23. [Google Scholar] [CrossRef]
  59. Pietralla, M.; Hotz, R.; Engst, T.; Siems, R. Chain Direction Elastic Modulus of PE Crystal and Interlamellar Force Constant of N-Alkane Crystals from RAMAN Measurements. J. Polym. Sci. B Polym. Phys. 1997, 35, 47–57. [Google Scholar] [CrossRef]
  60. Holliday, L.; White, J. The Stiffness of Polymers in Relation to Their Structure. Pure Appl. Chem. 1971, 26, 545–582. [Google Scholar] [CrossRef]
  61. Jaeken, J.W.; Cottenier, S. Solving the Christoffel Equation: Phase and Group Velocities. Comput. Phys. Commun. 2016, 207, 445–451. [Google Scholar] [CrossRef]
  62. Wu, X.; Vanderbilt, D.; Hamann, D.R. Systematic Treatment of Displacements, Strains, and Electric Fields in Density-Functional Perturbation Theory. Phys. Rev. B Condens. Matter Mater. Phys. 2005, 72, 035105. [Google Scholar] [CrossRef]
  63. Kurita, T.; Fukuda, Y.; Takahashi, M.; Sasanuma, Y. Crystalline Moduli of Polymers, Evaluated from Density Functional Theory Calculations under Periodic Boundary Conditions. ACS Omega 2018, 3, 4824–4835. [Google Scholar] [CrossRef] [PubMed]
  64. Barham, P.J.; Keller, A. Achievement of High-Modulus Polyethylene Fibres and the Modulus of Polyethylene Crystals. J. Polym. Sci. Polym. Lett. Ed. 1979, 17, 591–593. [Google Scholar] [CrossRef]
  65. Root, S.E.; Savagatrup, S.; Printz, A.D.; Rodriquez, D.; Lipomi, D.J. Mechanical Properties of Organic Semiconductors for Stretchable, Highly Flexible, and Mechanically Robust Electronics. Chem. Rev. 2017, 117, 6467–6499. [Google Scholar] [CrossRef] [PubMed]
  66. Yang, Y.; Jiang, H.; Zhou, Z.; Yang, J.; Wang, Y.; Bi, K. The Experimental Young’s Modulus of Polythiophene Nanofibers. Mater. Sci. Eng. B 2024, 299, 117014. [Google Scholar] [CrossRef]
  67. Matsuo, M.; Sawatari, C. Elastic Modulus of Polyethylene in the Crystal Chain Direction As Measured by X-Ray Diffraction. Macromolecules 1986, 19, 2036–2040. [Google Scholar] [CrossRef]
  68. Kobayashi, M.; Sakagami, K.; Tadokoro, H. Effects of Interlamellar Forces on Longitudinal Acoustic Modes of N-alkanes. J. Chem. Phys. 1983, 78, 6391–6398. [Google Scholar] [CrossRef]
  69. Tashiro, K.; Wu, G.; Kobayashi, M. Morphological Effect on the Raman Frequency Shift Induced by Tensile Stress Applied to Crystalline Polyoxymethylene and Polyethylene: Spectroscopic Support for the Idea of an Inhomogeneous Stress Distribution in Polymer Material. Polymer 1988, 29, 1768–1778. [Google Scholar] [CrossRef]
  70. Twisleton, J.F.; White, J.W.; Reynolds, P.A. Dynamical Studies of Fully Oriented Deuteropolyethylene by Inelastic Neturon Scattering. Polymer 1982, 23, 578–588. [Google Scholar] [CrossRef]
  71. Stuart, S.J.; Tutein, A.B.; Harrison, J.A. A Reactive Potential for Hydrocarbons with Intermolecular Interactions. J. Chem. Phys. 2000, 112, 6472–6486. [Google Scholar] [CrossRef]
  72. Zhang, Z.; Ouyang, Y.; Guo, Y.; Nakayama, T.; Nomura, M.; Volz, S.; Chen, J. Hydrodynamic Phonon Transport in Bulk Crystalline Polymers. Phys. Rev. B 2020, 102, 195302. [Google Scholar] [CrossRef]
  73. Wang, X.; Kaviany, M.; Huang, B. Phonon Coupling and Transport in Individual Polyethylene Chains: A Comparison Study with the Bulk Crystal. Nanoscale 2017, 9, 18022–18031. [Google Scholar] [CrossRef] [PubMed]
  74. Shulumba, N.; Hellman, O.; Minnich, A.J. Lattice Thermal Conductivity of Polyethylene Molecular Crystals from First-Principles Including Nuclear Quantum Effects. Phys. Rev. Lett. 2017, 119, 185901. [Google Scholar] [CrossRef] [PubMed]
  75. Chaput, L. Direct Solution to the Linearized Phonon Boltzmann Equation. Phys. Rev. Lett. 2013, 110, 265506. [Google Scholar] [CrossRef] [PubMed]
  76. Lindsay, L.; Broido, D.A.; Reinecke, T.L. Ab Initio Thermal Transport in Compound Semiconductors. Phys. Rev. B Condens. Matter Mater. Phys. 2013, 87, 165201. [Google Scholar] [CrossRef]
  77. Simoncelli, M.; Marzari, N.; Mauri, F. Unified Theory of Thermal Transport in Crystals and Glasses. Nat. Phys. 2019, 15, 809–813. [Google Scholar] [CrossRef]
  78. Davis, G.T.; Eby, R.K.; Colson, J.P. Thermal Expansion of Polyethylene Unit Cell: Effect of Lamella Thickness. J. Appl. Phys. 1970, 41, 4316–4326. [Google Scholar] [CrossRef]
  79. Kamencek, T.; Schrode, B.; Resel, R.; Ricco, R.; Zojer, E. Understanding the Origin of the Particularly Small and Anisotropic Thermal Expansion of MOF-74. Adv. Theory Simul. 2022, 5, 2200031. [Google Scholar] [CrossRef]
  80. Togo, A.; Chaput, L.; Tanaka, I. Distributions of Phonon Lifetimes in Brillouin Zones. Phys. Rev. B Condens. Matter Mater. Phys. 2015, 91, 094306. [Google Scholar] [CrossRef]
  81. Togo, A.; Chaput, L.; Tadano, T.; Tanaka, I. Implementation Strategies in Phonopy and Phono3py. J. Phys. Condens. Matter 2023, 35, 353001. [Google Scholar] [CrossRef]
  82. Bedoya-Martínez, N.; Schrode, B.; Jones, A.O.F.; Salzillo, T.; Ruzié, C.; Demitri, N.; Geerts, Y.H.; Venuti, E.; Della Valle, R.G.; Zojer, E.; et al. DFT-Assisted Polymorph Identification from Lattice Raman Fingerprinting. J. Phys. Chem. Lett. 2017, 8, 3690–3695. [Google Scholar] [CrossRef]
  83. George, J.; Wang, R.; Englert, U.; Dronskowski, R. Lattice Thermal Expansion and Anisotropic Displacements in Urea, Bromomalonic Aldehyde, Pentachloropyridine, and Naphthalene. J. Chem. Phys. 2017, 147, 074112. [Google Scholar] [CrossRef] [PubMed]
  84. VASP Manual ML_CX. Available online: https://www.vasp.at/wiki/index.php/ML_CX (accessed on 4 June 2024).
  85. Togo, A.; Tanaka, I. First Principles Phonon Calculations in Materials Science. Scr. Mater. 2015, 108, 1–5. [Google Scholar] [CrossRef]
  86. Togo, A. Phono3py Documentation. Available online: https://phonopy.github.io/phono3py/tips.html (accessed on 19 December 2023).
  87. Le Page, Y.; Saxe, P. Symmetry-General Least-Squares Extraction of Elastic Data for Strained Materials from Ab Initio Calculations of Stress. Phys. Rev. B Condens. Matter Mater. Phys. 2002, 65, 104104. [Google Scholar] [CrossRef]
  88. Baroni, S.; De Gironcoli, S.; Dal Corso, A.; Giannozzi, P. Phonons and Related Crystal Properties from Density-Functional Perturbation Theory. Rev. Mod. Phys. 2001, 73, 515. [Google Scholar] [CrossRef]
  89. Zhang, T.; Luo, T. Morphology-Influenced Thermal Conductivity of Polyethylene Single Chains and Crystalline Fibers. J. Appl. Phys. 2012, 112, 094304. [Google Scholar] [CrossRef]
  90. Sæther, S.; Falck, M.; Zhang, Z.; Lervik, A.; He, J. Thermal Transport in Polyethylene: The Effect of Force Fields and Crystallinity. Macromolecules 2021, 54, 6563–6574. [Google Scholar] [CrossRef]
  91. Blum, V.; Gehrke, R.; Hanke, F.; Havu, P.; Havu, V.; Ren, X.; Reuter, K.; Scheffler, M. Ab Initio Molecular Simulations with Numeric Atom-Centered Orbitals. Comput. Phys. Commun. 2009, 180, 2175–2196. [Google Scholar] [CrossRef]
  92. Hermann, J.; Tkatchenko, A. Density Functional Model for van Der Waals Interactions: Unifying Many-Body Atomic Approaches with Nonlocal Functionals. Phys. Rev. Lett. 2020, 124, 146401. [Google Scholar] [CrossRef]
  93. Samdal, S.; Samuelsen, E.J.; Volden, H.V. Molecular Conformation of 2,2′-Bithiophene Determined by Gas Phase Electron Diffraction and Ab Initio Calculations. Synth. Met. 1993, 59, 259–265. [Google Scholar] [CrossRef]
  94. Kuhn, H.W. The Hungarian Method for the Assignment Problem. Nav. Res. Logist. Q. 1955, 2, 83–97. [Google Scholar] [CrossRef]
  95. Stukowski, A. Visualization and Analysis of Atomistic Simulation Data with OVITO—The Open Visualization Tool. Model. Simul. Mater. Sci. Eng. 2009, 18, 015012. [Google Scholar] [CrossRef]
  96. Schaufele, R.F.; Shimanouchi, T. Longitudinal Acoustical Vibrations of Finite Polymethylene Chains. J. Chem. Phys. 1967, 47, 3605–3610. [Google Scholar] [CrossRef]
  97. Feldkamp, L.A.; Venkataraman, G.; King, J.S. Dispersion Relation for Skeletal Vibrations in Deuterated Polyethylene. In Proceedings of the IAEA Symposium on Neutron Inelastic Scattering, Copenhagen, Denmark, 20–25 May 1968. [Google Scholar]
  98. Brenner, D.W.; Shenderova, O.A.; Harrison, J.A.; Stuart, S.J.; Ni, B.; Sinnott, S.B. A Second-Generation Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. J. Phys. Condens. Matter 2002, 14, 783–802. [Google Scholar] [CrossRef]
  99. Case, D.A.; Cheatham, T.E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K.M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R.J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26, 1668–1688. [Google Scholar] [CrossRef]
  100. Guo, Z.; Han, Z.; Feng, D.; Lin, G.; Ruan, X. Sampling-Accelerated Prediction of Phonon Scattering Rates for Converged Thermal Conductivity and Radiative Properties. NPJ Comput. Mater. 2024, 10, 31. [Google Scholar] [CrossRef]
Figure 2. First Brillouin zones of PE (a) and P3HT (b). The high symmetry points X, Y, Z, and Γ are labelled. The paths Γ-X, Γ-Y, and Γ-Z are shown with red arrows and correspond to the phonon band structure paths shown in Section 2.5. In all systems, Γ-Z is parallel to the polymer axis. The Bravais lattice is primitive orthorhombic for PE (and also for PT—not shown here) and primitive monoclinic for P3HT. Thus, in P3HT, the Γ-X direction (red arrow) is slightly inclined relative to the real space direction x (purple arrow).
Figure 2. First Brillouin zones of PE (a) and P3HT (b). The high symmetry points X, Y, Z, and Γ are labelled. The paths Γ-X, Γ-Y, and Γ-Z are shown with red arrows and correspond to the phonon band structure paths shown in Section 2.5. In all systems, Γ-Z is parallel to the polymer axis. The Bravais lattice is primitive orthorhombic for PE (and also for PT—not shown here) and primitive monoclinic for P3HT. Thus, in P3HT, the Γ-X direction (red arrow) is slightly inclined relative to the real space direction x (purple arrow).
Molecules 29 03724 g002
Figure 3. Illustration of how the carbon atoms of P3HT can be split into at most six atom types. Sulfur and hydrogen atoms are white. Carbon atoms with up to six different chemical environments can be distinguished: sp2-hybridized carbon atoms bonded to two other carbon atoms and one sulfur (orange), bonded to two carbon atoms and one hydrogen (dark red), and bonded to three carbon atoms (light red); furthermore, sp3-hybridized carbon atoms bonded to the sp2-hybridized backbone (medium bright blue), sp3-hybridized carbons in the middle of the sidechain (light blue), and at the end of the side chain (dark blue). The unit cell is shown as a black box and contains two P3HT chains, but for the sake of clarity, only one is shown here.
Figure 3. Illustration of how the carbon atoms of P3HT can be split into at most six atom types. Sulfur and hydrogen atoms are white. Carbon atoms with up to six different chemical environments can be distinguished: sp2-hybridized carbon atoms bonded to two other carbon atoms and one sulfur (orange), bonded to two carbon atoms and one hydrogen (dark red), and bonded to three carbon atoms (light red); furthermore, sp3-hybridized carbon atoms bonded to the sp2-hybridized backbone (medium bright blue), sp3-hybridized carbons in the middle of the sidechain (light blue), and at the end of the side chain (dark blue). The unit cell is shown as a black box and contains two P3HT chains, but for the sake of clarity, only one is shown here.
Molecules 29 03724 g003
Figure 4. Influence of the training data and the chosen level on the accuracy of MTPs. Accuracy is measured by the root mean square deviation (RMSD) between MTP- and DFT-calculated phonon frequencies sampled in the whole Brillouin zone up to 12.5 THz (RMSDphonon) and by the RMSD of forces on atoms calculated for validation structures generated in a 300 K active learning MD run (RMSDMD). Panels (a,b) show RMSDphonon and RMSDMD for different training data sets for PE (blue) and PT (orange) for the MTP level set to 22. The values on the x-axis denote the maximum temperature in the generation of the respective reference data during active learning runs (for details see main text). Panel (c) shows RMSDphonon for PE for different MTP levels (generated with 15–100 K training data), and panel (d) RMSDMD for PE for different MTP levels (again generated with 15–100 K training data). For each data set/level, five MTPs are parametrized and the associated RMSD values are shown as small, filled symbols. Datapoints for MTPs with median RMSDphonon/RMSDMD are connected by solid lines, while data for the “best” MTPs are denoted by large, dark, open symbols and are connected by dashed lines. The numerical values for these plots are provided in Tables S3–S5 of the Supplementary Materials.
Figure 4. Influence of the training data and the chosen level on the accuracy of MTPs. Accuracy is measured by the root mean square deviation (RMSD) between MTP- and DFT-calculated phonon frequencies sampled in the whole Brillouin zone up to 12.5 THz (RMSDphonon) and by the RMSD of forces on atoms calculated for validation structures generated in a 300 K active learning MD run (RMSDMD). Panels (a,b) show RMSDphonon and RMSDMD for different training data sets for PE (blue) and PT (orange) for the MTP level set to 22. The values on the x-axis denote the maximum temperature in the generation of the respective reference data during active learning runs (for details see main text). Panel (c) shows RMSDphonon for PE for different MTP levels (generated with 15–100 K training data), and panel (d) RMSDMD for PE for different MTP levels (again generated with 15–100 K training data). For each data set/level, five MTPs are parametrized and the associated RMSD values are shown as small, filled symbols. Datapoints for MTPs with median RMSDphonon/RMSDMD are connected by solid lines, while data for the “best” MTPs are denoted by large, dark, open symbols and are connected by dashed lines. The numerical values for these plots are provided in Tables S3–S5 of the Supplementary Materials.
Molecules 29 03724 g004
Figure 5. Influence of increasing the number of atom types on MTP accuracy for P3HT. Panel (a) shows RMSDphonon, panel (b) RMSDMD for treating all carbon atoms in P3HT the same (one atom type), when distinguishing between sp3-hybridized carbons and sp2-hybridized ones, bonded either only to other carbons or to one sulfur and two carbons (three atom types), and when distinguishing between all chemically inequivalent carbons (six atom types), as shown in Figure 3. The level of the MTP is set to 22 and the training data are taken from an active learning run between 15 K and 500 K. For each number of atom types, five MTPs are parametrized and the associated MTPs are shown as small, filled symbols. MTPs with median RMSDphonon/RMSDMD are connected by solid lines, while data for the “best” MTPs are denoted by large, dark, open symbols and are connected by dashed lines. Numerical values associated with these plots are given in Table S9 of the Supplementary Materials.
Figure 5. Influence of increasing the number of atom types on MTP accuracy for P3HT. Panel (a) shows RMSDphonon, panel (b) RMSDMD for treating all carbon atoms in P3HT the same (one atom type), when distinguishing between sp3-hybridized carbons and sp2-hybridized ones, bonded either only to other carbons or to one sulfur and two carbons (three atom types), and when distinguishing between all chemically inequivalent carbons (six atom types), as shown in Figure 3. The level of the MTP is set to 22 and the training data are taken from an active learning run between 15 K and 500 K. For each number of atom types, five MTPs are parametrized and the associated MTPs are shown as small, filled symbols. MTPs with median RMSDphonon/RMSDMD are connected by solid lines, while data for the “best” MTPs are denoted by large, dark, open symbols and are connected by dashed lines. Numerical values associated with these plots are given in Table S9 of the Supplementary Materials.
Molecules 29 03724 g005
Figure 6. Independent elements of the elastic tensor calculated with MTPMD and MTPphonon compared to the results obtained with DFT. All simulations are performed applying the “clamped ion method” as explained in the main text. The investigated materials are (a) PE and (b) PT. The black line with a slope of 1 indicates perfect agreement between MTP and DFT. The displayed data are obtained with the “best” MTPs as defined at the end of Section 2.2. The large open symbols represent Cxx, Cyy, and Czz, while the small symbols represent the other independent components of the elastic tensor. An analogous plot containing the mean values over elastic constants calculated with all parametrized (typically five) MTPs, and corresponding error bars can be found in Figure S6.
Figure 6. Independent elements of the elastic tensor calculated with MTPMD and MTPphonon compared to the results obtained with DFT. All simulations are performed applying the “clamped ion method” as explained in the main text. The investigated materials are (a) PE and (b) PT. The black line with a slope of 1 indicates perfect agreement between MTP and DFT. The displayed data are obtained with the “best” MTPs as defined at the end of Section 2.2. The large open symbols represent Cxx, Cyy, and Czz, while the small symbols represent the other independent components of the elastic tensor. An analogous plot containing the mean values over elastic constants calculated with all parametrized (typically five) MTPs, and corresponding error bars can be found in Figure S6.
Molecules 29 03724 g006
Figure 7. Panel (ac) phonon band structures calculated with DFT (solid blue lines) and the “best” MTPphonon (dashed orange lines) in the low-frequency region. Materials are (a) PE, (b) PT, and (c) P3HT. The frequency range for P3HT is half that for PE and PT, such that the many bands are resolved more clearly. A corresponding plot for the extended frequency range is given in Figure S9. Panel (d) compares the phonon band structures of PE calculated with DFT and with the AIREBO potential.
Figure 7. Panel (ac) phonon band structures calculated with DFT (solid blue lines) and the “best” MTPphonon (dashed orange lines) in the low-frequency region. Materials are (a) PE, (b) PT, and (c) P3HT. The frequency range for P3HT is half that for PE and PT, such that the many bands are resolved more clearly. A corresponding plot for the extended frequency range is given in Figure S9. Panel (d) compares the phonon band structures of PE calculated with DFT and with the AIREBO potential.
Molecules 29 03724 g007
Figure 8. Phonon bands along the chain direction in PE for which the lifetimes of the phonons are encoded by the size of the symbols and by their colors. They have been calculated with (a) DFT and with (b) the “best” MTPphonon. While the color scale is logarithmic, the areas of the spherical symbols are scaled linearly with the phonon lifetimes. Data points depicted as smaller circles are plotted in front, such that they are not hidden by larger circles. We note that a similarly good agreement is obtained for the other four variants of MTPphonon that were trained.
Figure 8. Phonon bands along the chain direction in PE for which the lifetimes of the phonons are encoded by the size of the symbols and by their colors. They have been calculated with (a) DFT and with (b) the “best” MTPphonon. While the color scale is logarithmic, the areas of the spherical symbols are scaled linearly with the phonon lifetimes. Data points depicted as smaller circles are plotted in front, such that they are not hidden by larger circles. We note that a similarly good agreement is obtained for the other four variants of MTPphonon that were trained.
Molecules 29 03724 g008
Figure 9. Contributions to the thermal conductivity of PE along the chain (κzz) in the RTA for different frequency ranges. Data calculated with DFT are compared to those obtained with the “best” MTPphonon.
Figure 9. Contributions to the thermal conductivity of PE along the chain (κzz) in the RTA for different frequency ranges. Data calculated with DFT are compared to those obtained with the “best” MTPphonon.
Molecules 29 03724 g009
Figure 10. Thermal expansion of PE along its lattice vectors (a) a, (b) b, and (c) c. The purple line labelled ’MTPMD, “best”’ is the MTPMD with lowest RMSDMD. The blue line labelled “MTPMD, mean ± 2σ” is the mean from five molecular dynamics runs each with a different MTPMD; whereby the error bars are two times the standard deviation (95% confidence interval). Experimental data from Davis et al. [78], Shen et al. [55], and Takahashi [53] are shown. Davis measured two samples, which they called “32 days” and “2 h” [78]. These names refer to how long the samples remained in constant temperature baths during their production. Shen et al. measured PE single crystals.
Figure 10. Thermal expansion of PE along its lattice vectors (a) a, (b) b, and (c) c. The purple line labelled ’MTPMD, “best”’ is the MTPMD with lowest RMSDMD. The blue line labelled “MTPMD, mean ± 2σ” is the mean from five molecular dynamics runs each with a different MTPMD; whereby the error bars are two times the standard deviation (95% confidence interval). Experimental data from Davis et al. [78], Shen et al. [55], and Takahashi [53] are shown. Davis measured two samples, which they called “32 days” and “2 h” [78]. These names refer to how long the samples remained in constant temperature baths during their production. Shen et al. measured PE single crystals.
Molecules 29 03724 g010
Figure 11. Histograms of the deviations to DFT reference data for equilibrium energies, ΔE, forces on individual atoms, ΔF, and stresses, Δσ. The compared quantities were obtained with the “best” MTPMD (dark blue shading) and with the AIREBO potential (light blue shading). Panels (a,c,e) show the data for the “best” MTPMD. Panels (b,d,f) show the same data for MTPMD (albeit on a different scale) and compare them to data obtained with the AIREBEO potential. The number of counts was scaled by a factor of five for the AIREBO potential to improve visibility. Gray-shaded areas represent ranges where errors are within 5% and 10% of the absolute DFT values for forces and stresses, and within 1% for energy differences.
Figure 11. Histograms of the deviations to DFT reference data for equilibrium energies, ΔE, forces on individual atoms, ΔF, and stresses, Δσ. The compared quantities were obtained with the “best” MTPMD (dark blue shading) and with the AIREBO potential (light blue shading). Panels (a,c,e) show the data for the “best” MTPMD. Panels (b,d,f) show the same data for MTPMD (albeit on a different scale) and compare them to data obtained with the AIREBEO potential. The number of counts was scaled by a factor of five for the AIREBO potential to improve visibility. Gray-shaded areas represent ranges where errors are within 5% and 10% of the absolute DFT values for forces and stresses, and within 1% for energy differences.
Molecules 29 03724 g011
Table 1. Comparison of the lengths of the lattice vector a, b, and c and the monoclinic angle (between the a- and b-vectors), α, calculated with MTPs and with DFT by minimizing the energy, as detailed in Section 3. For the two types of MTPs (MTPphonon and MTPMD), the results for the “best” MTP (with the lowest RMSDphonon/RMSDMD values) are given. Additionally, the arithmetic mean values for the five parametrized potentials for each system are reported. The uncertainties of the mean MTP values only refer to the statistical errors caused by the MTP parametrization, as they are the standard deviation for the data from the five MTPs parametrized with different seeds. For MTPphonon in the case of P3HT, no uncertainties are reported, as there only two MTPs are parametrized due to the associated high computational costs associated with generating level 26 potentials for such a complex material. Experimental data are taken from the literature. Temperatures at which the experiments are performed are given in parentheses (with RT meaning room temperature) [41,52,53,54,55]. For the specific conformation of P3HT studied here (see Section “Structures of Polyethylene, Polythiophene and Poly(3-Hexylthiophene-2,5-Diyl) (P3HT)”), no experimental data are found.
Table 1. Comparison of the lengths of the lattice vector a, b, and c and the monoclinic angle (between the a- and b-vectors), α, calculated with MTPs and with DFT by minimizing the energy, as detailed in Section 3. For the two types of MTPs (MTPphonon and MTPMD), the results for the “best” MTP (with the lowest RMSDphonon/RMSDMD values) are given. Additionally, the arithmetic mean values for the five parametrized potentials for each system are reported. The uncertainties of the mean MTP values only refer to the statistical errors caused by the MTP parametrization, as they are the standard deviation for the data from the five MTPs parametrized with different seeds. For MTPphonon in the case of P3HT, no uncertainties are reported, as there only two MTPs are parametrized due to the associated high computational costs associated with generating level 26 potentials for such a complex material. Experimental data are taken from the literature. Temperatures at which the experiments are performed are given in parentheses (with RT meaning room temperature) [41,52,53,54,55]. For the specific conformation of P3HT studied here (see Section “Structures of Polyethylene, Polythiophene and Poly(3-Hexylthiophene-2,5-Diyl) (P3HT)”), no experimental data are found.
a [Å]b [Å]c [Å]α [°]
Polyethylene
DFT7.0744.8532.55490
MTPphonon, “best”7.0624.8472.554390
MTPphonon, mean7.060 ± 0.0034.849 ± 0.0022.5543 ± 0.000190
MTPMD, “best”7.0504.8532.554190
MTPMD, mean7.052 ± 0.0134.849 ± 0.0152.5542 ± 0.000290
Experiment (4 K) [52]7.1214.8512.54890
Experiment (10 K) [53]7.1204.842-90
Experiment (10 K) [55]7.164.862.53490
Polythiophene
DFT7.5305.5427.78590
MTPphonon, “best”7.4675.5087.78290
MTPphonon, mean7.472 ± 0.0035.506 ± 0.0027.782 ± 0.00190
MTPMD, “best”7.5005.4887.78290
MTPMD, mean7.449 ± 0.0435.524 ± 0.0277.781 ± 0.00190
Experiment 1 [54]7.805.558.0390
Experiment (RT) [41]7.795.53--
P3HT
DFT7.57514.7317.81688.75
MTPphonon, “best”7.56114.7057.81588.81
MTPphonon, mean7.56414.7037.81688.74
MTPMD, “best”7.57314.6577.81688.7
MTPMD, mean7.569 ± 0.00714.675 ± 0.0217.816 ± 0.00189.13
±0.46
1 Temperature not specified.
Table 3. RMSDphonon values for the “best” and median MTPs obtained when comparing DFT- and MTP-calculated phonon band structures for a homogeneous sampling of reciprocal space. The values are given in THz and in parentheses also in cm−1. Note that the values reported here are obtained by fixing the unit cell to the DFT-optimized parameters to avoid simultaneously assessing the ability of the MTPs reproducing unit cell parameters and phonon properties. RMSDphonon values for MTP optimized unit cells are contained in Table S15. Due to the high computational cost, only two MTPphonon are parametrized for P3HT. Therefore, there is no median MTP in that case.
Table 3. RMSDphonon values for the “best” and median MTPs obtained when comparing DFT- and MTP-calculated phonon band structures for a homogeneous sampling of reciprocal space. The values are given in THz and in parentheses also in cm−1. Note that the values reported here are obtained by fixing the unit cell to the DFT-optimized parameters to avoid simultaneously assessing the ability of the MTPs reproducing unit cell parameters and phonon properties. RMSDphonon values for MTP optimized unit cells are contained in Table S15. Due to the high computational cost, only two MTPphonon are parametrized for P3HT. Therefore, there is no median MTP in that case.
PE [THz and
(cm−1)]
PT [THz and
(cm−1)]
P3HT [THz and
(cm−1)]
MTPphonon, “best”0.043 (1.43)0.029 (0.97)0.036 (1.20)
MTPphonon, median0.053 (1.77)0.032 (1.07)-
MTPMD, “best”0.152 (5.07)0.080 (2.67)0.059 (1.97)
MTPMD, median0.163 (5.44)0.086 (2.87)0.089 (2.97)
Table 4. Thermal conductivity of PE at 300 K, calculated with the BTE employing DFT and MTPs. RTA refers to the use of the relaxation time approximation. “Full BTE” refers to BTE calculations disregarding the RTA. For the result of the mean MTP, the arithmetic mean is taken over the thermal conductivities calculated with all five MTPs. The uncertainty is the standard deviation calculated for the thermal conductivities obtained with these five MTPs.
Table 4. Thermal conductivity of PE at 300 K, calculated with the BTE employing DFT and MTPs. RTA refers to the use of the relaxation time approximation. “Full BTE” refers to BTE calculations disregarding the RTA. For the result of the mean MTP, the arithmetic mean is taken over the thermal conductivities calculated with all five MTPs. The uncertainty is the standard deviation calculated for the thermal conductivities obtained with these five MTPs.
κxx [Wm−1K−1]κyy [Wm−1K−1]κzz [Wm−1K−1]
RTA
DFT0.540.46306
MTPphonon, “best”0.710.54307
MTPphonon, mean0.67 ± 0.090.53 ± 0.08292 ± 33
Full BTE
DFT0.520.47398
MTPphonon, “best”0.710.58408
MTPphonon, mean0.66 ± 0.090.58 ± 0.09388 ± 40
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Reicht, L.; Legenstein, L.; Wieser, S.; Zojer, E. Designing Accurate Moment Tensor Potentials for Phonon-Related Properties of Crystalline Polymers. Molecules 2024, 29, 3724. https://doi.org/10.3390/molecules29163724

AMA Style

Reicht L, Legenstein L, Wieser S, Zojer E. Designing Accurate Moment Tensor Potentials for Phonon-Related Properties of Crystalline Polymers. Molecules. 2024; 29(16):3724. https://doi.org/10.3390/molecules29163724

Chicago/Turabian Style

Reicht, Lukas, Lukas Legenstein, Sandro Wieser, and Egbert Zojer. 2024. "Designing Accurate Moment Tensor Potentials for Phonon-Related Properties of Crystalline Polymers" Molecules 29, no. 16: 3724. https://doi.org/10.3390/molecules29163724

APA Style

Reicht, L., Legenstein, L., Wieser, S., & Zojer, E. (2024). Designing Accurate Moment Tensor Potentials for Phonon-Related Properties of Crystalline Polymers. Molecules, 29(16), 3724. https://doi.org/10.3390/molecules29163724

Article Metrics

Back to TopTop