Tutorial: AI-assisted exploration and active design of polymers with high intrinsic thermal conductivity

Polymers have attracted extensive attention in various fields such as energy,^1,2 environment,^3,4 electronics,⁵ biologies,⁶ medicine,⁷ and engineering⁸ thanks to their light weight, low cost, excellent mechanical ductility, superior biocompatibility, and good chemical and thermal stability.^9–11 However, intrinsic bulk polymers are thermally insulating and have a low thermal conductivity (TC) in a narrow range of 0.1–0.3 W m⁻¹ K⁻¹,^12–14 which restricts the heat dissipation of organic equipment and severely obstructs the process of miniaturization and integration of flexible electronic and optoelectronic devices.^15–22 Achieving high TC in polymers for industrial applications is an urgent demand, and some progress has been realized recently.^23,24 Different fabrication techniques such as micromechanical stretching,^25–28 electrostatic spinning,^29–31 and nanoscale templating^32,33 have been employed in the exploitation of polymers, which effectively improve the ordering and crystallinity of the polymer chains and, thus, exhibit high TC. Moreover, the construction of thermal networks by engineering interchain interactions such as hydrogen bonding,^34, $π−π$ stacking,³⁵ and side chain modifying³⁶ in polymer blends and copolymers has also been demonstrated to be beneficial for enhanced TC. Taking polyethylene (PE) as an example, the TC of PE films²⁷ and nanofibers²⁶ by mechanical stretching was found to be as high as ∼62 and ∼104 W m⁻¹ K⁻¹, respectively, over two or three orders-of-magnitude greater than that of typical polymers.

Despite the engineered polymers that can be produced experimentally to achieve increased TC, this requires a strong chemical background from investigators and is limited by process and characterization instruments. Further, the applicability of different techniques is restricted, e.g., micromechanical stretching is unsuitable for brittle polymers.³⁷ Accompanied by the evolution of high-performance computers and the revolution of multiscale simulation methodology, in silico experiments are playing an important role in the study of thermal transport in polymers.^38–42 Computational approaches including first-principles calculations^43,44 and molecular dynamics (MD) simulations⁴⁵ have led the way in revealing the effect of polymer nanostructures on TC. The first-principles calculation to TC is based on the computation of interatomic force constants via density-functional theory (DFT). On this basis, all relevant phonon properties can be calculated using lattice dynamics and the Boltzmann transport equation.⁴⁴ This method has been successfully applied to molecular crystals such as PE,⁴³ polyvinylidene fluoride,⁴³ and polythiophene,⁴⁴ but is challenging to apply in amorphous systems, owing to quantum nuclear motion and their complex primitive cells.

MD simulations employ classical force fields combined with Newtonian mechanics and statistical physics to derive the macroscopic properties of systems.⁴⁶ MD simulations can handle polymer systems containing tens of thousands of atoms, which are widely used in polymer thermal transport studies, not only to predict the TC of hierarchical structures, but also to probe the linkages between microstructures and TCs. Using MD simulations, an individual PE chain exhibits a very high TC of 350 W m⁻¹ K⁻¹, even divergent in some cases.^47–49 These results are encouraging and inspire researchers to make further efforts to develop polymers with high TC. Moreover, properties of polymers, such as molecular weight,^50,51 chain length,^52–56 side chains,^36,57–59 and chain conformation,^51,60–63 intra-chain effects, such as bonds,^64–66 angles,^67,68 and dihedrals,^69–71 inter-chain behaviors, such as molecular cross-linking,^72–77 hydrogen bonding networks,^78–81 and $π−π$ stacking,^58,82–84 on TC were extracted in some separate MD simulations. However, current investigations on the thermal transport mechanisms of polymers have mainly focused on common polymers such as PE, polytetrafluoroethylene (PTFE), polyvinylidene fluoride, and other conjugated structures.⁸⁵ 

For a long time, researchers have been working on exploring quantitative structure–activity relationships (QSAR) from chemical data, which, in turn, has enabled the rational design of innovative materials, including polymers and inorganic and biological components.^86,87 QSAR models reveal empirical, linear, or non-linear relationships between descriptors extracted from chemical structures and computational/experimental properties or activities.⁸⁸ As modern research methods have facilitated the proliferation in the amount of chemical data, the data-driven research paradigm is of critical importance for QSAR modeling.^89–91 In terms of polymers, the chemical space is enormous,⁹² corresponding to potential candidates of small organic molecules as many as 10⁶⁰ whereas the known organic compounds are more than 10⁸ recorded in the PubChem database.^93–95 Moreover, virtual chemical reactions⁹⁶ or generative algorithms⁹⁷ for small molecules can create a nearly infinite chemical space. Polymer informatics is a data-centric technology equipped with artificial intelligence and machine learning (ML) as powerful engines to accelerate the optimization of organic materials and the development of novel macromolecules.^98–101 Polymer informatics has been effective in facilitating polymer innovation, and has achieved a series of successful applications, involving optical,^102–104 electrical,^105–108 thermal,^109–112 mechanical,^113,114 and other properties.^115–118 Over the past five years, several efforts have been made to apply ML to the exploitation of polymers with high TC (TC > 0.30 W m⁻¹ K⁻¹), with essential contributions for expanding the potential candidates and revealing underlying physical mechanisms.^{85,94,119–124}

In this Tutorial, we introduce some development paradigms combining high-throughput MD simulations and ML for high TC polymers in the hope of inspiring new researchers who are interested in becoming involved in this field. We start with a description of three core components in polymer informatics of polymer datasets, TC simulation methods, and polymer representations in Sec. II. Following this, interpretable regression models are constructed for mapping polymer microstructures to TCs in Sec. III. Next, active optimization algorithms are utilized for the design of polymers with high thermal conductivity in Sec. IV, including single- and multi-objective cases. Our conclusions and outlooks for this area are provided in Sec. V.

The principle of polymer informatics is to establish patterns from a sufficient amount of existing or generated polymer data, thus facilitating the design/discovery of new functional polymers with improved target properties.⁹²  Figure 1 illustrates a mainstream informatics framework for the development of polymers with high TC, consisting of four elements: (1) polymer data sets; (2) polymer modeling and TC calculations; (3) feature engineering; and (4) informatics algorithms. In the following, we explain the implementation of the design framework in these four aspects.

Well-organized and clean data are a fundamental prerequisite for the accomplishment of high-fidelity informatics algorithms. To promote the openness and sharing of scientific data, the findable, accessible, interoperable, and reusable (FAIR) principle has been proposed for data storage and management,^125–127 which also greatly contributed to the progression of material genetic engineering.¹²⁸ Inorganic materials have a more extensive and accessible database than polymeric materials, such as Materials Project,¹²⁹ Atomly,¹³⁰ ICSD,¹³¹ and AFLOW.¹³² These inorganic databases usually follow the FAIR guidelines and provide application programming interfaces (APIs) for access and download. However, the development process of polymer databases is relatively slow and is mostly limited to access. This is attributed to the fact that large polymer primitive cells make experimentation/computation difficult and costly,¹³³ and empirical nomenclature¹³⁴ and diverse forms of expression (strings and graphs) of polymers prevent text-mining techniques from obtaining property data from the published scientific literature.¹³⁵ 

A few representative databases include PoLyInfo,¹³⁶ Khazana,¹³⁷ Polymers: A Property Database,¹³⁸ Polymer Property Predictor and Database,¹³⁹ CAMPUS,¹⁴⁰ and PI1M,⁹⁷ which are listed in Table I. PoLyInfo¹³⁶ is one of the largest polymer experimental databases, with over 18 000 homopolymers and 7000 copolymers from about 20 000 scientific literature, including hundreds of properties such as refractive index, dielectric constant, and glass transition temperature. Nevertheless, the abundance of different properties is quite variable, for example, more than 8000 homopolymers contain recorded glass transition temperature, while only 84 homopolymers expose TC (accessed 14 December, 2023). Moreover, PoLyInfo is contributed by the National Institute for Materials Science of Japan and is prohibited for downloading large amounts of data. To address this concern, Hayashi et al.¹⁴¹ selected 1070 amorphous polymers from the PoLyInfo database and calculated the associated 15 properties, including TC, using all-atom classical MD simulations. Ma and Luo⁹⁷ trained a recurrent neural network based on ∼12 000 polymers collected from PoLyInfo, and then generated ∼1 × 10⁶ polymers to form a benchmark database of PI1M, which covers a similar chemical space as the training data sets.

Databases of computational properties of polymers are rarer than experimental databases, of which Khazana¹³⁷ is a typical one. The Khazana database was supported by the Ramprasad's group, including DFT computed refractive index, dielectric constant, and bandgap, and was used as training data sets to build an open informatics platform of Polymer Genome.¹⁴² In Table I, we also describe three databases containing small molecule data sets, namely, PubChem,⁹³ eMolecules,¹⁴³ and Material Project.¹⁴⁴ With a grasp of the laws of chemical reactions, virtual reactions can convert small molecules into polymers through ML^145,146 or specific grammar rules.^96,116,147 Yang et al.¹¹⁶ established a polymer data set containing more than 8 × 10⁶ hypothetical polyimides formed by known dianhydride and diamine/diisocyanate pairs from PubChem.⁹³ Kim et al.⁹⁶ released a vast polymer database known as the Open Macromolecular Genome, which began with approximately 24 × 10⁶ potential reactant molecules of eMolecules¹⁴³ and then formed synthesizable polymer chemistries compatible with 17 polymerization reactions that cover a variety of step growth, chain growth, ring opening, and metathesis polymerization reactions.

Polymer simulation includes modeling, force field assignment, equilibrium simulation, and TC calculations. While some discrete software such as RDKit,¹⁴⁸ Materials Studio,¹⁴⁹ CHARMM-GUI,¹⁵⁰ Packmol,¹⁵¹ Moltemplate,¹⁵² Pysimm,¹⁵³ Polymer Structure Predictor (PSP),¹⁵⁴ Enhanced Monte Carlo (EMC),¹⁵⁵ and Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)¹⁵⁶ enable the above processes to be realized collaboratively, open-source toolkits that facilitate the building of the entire workflow are of great significance in generating the training data required for polymer informatics. The Polymer Molecular Dynamics (PMD) package¹⁵⁷ integrates with PSP, EMC, and LAMMPS software to realize polymer modeling and high-throughput MD simulations for various properties, including glass transition temperature, viscosity, TC, and so on. RadonPy¹⁴¹ is a robust open-source Python library that fully automates the calculation of various properties of polymers using all-atom classical MD simulations. The entire simulation process, including modeling and equilibrium and non-equilibrium MD simulations, can be automated by taking only simplified molecular-input line-entry system (SMILES)¹⁵⁸ string of the polymer repeating unit as input. Furthermore, a computational database has been released containing more than 1000 unique amorphous polymers with various thermophysical properties calculated using RadonPy. Given the powerful and convenient capability of RadonPy, the polymer simulation and ML training data for this Tutorial relies primarily on RadonPy and its derived computational data set.

For obtaining a simulation cell, the single polymer chain was duplicated into ten copies by translational and rotational operations to prevent overlap with each other and was placed in a large box with a density of ∼0.05 g/cm³. The packing simulation was performed to increase the density of the amorphous systems, to be adjusted to a suitable level. An NVT (constant number of atoms, volume, and temperature) simulation with a Nosé−Hoover thermostat was applied to the system in three sequential stages at a temperature of 300 K, from 300 to 700 K, and held at 700 K, under periodic boundary conditions (PBC) and a time step of 1 fs. Each of these stages took 1 ns, and all the bonds and angles were constrained by the SHAKE algorithm, resulting in a packaged cell with a density of around 0.80 g/cm³.

Equilibrium simulation was executed for the structural relaxation of amorphous polymers, which follows the 21-step compression/relaxation scheme.¹⁶⁶ During the simulation, by combining NVT and NPT (constant number of atoms, pressure, and temperature) simulations with a Noose–Hoover thermostat, the temperature rise to 600 K and fall to 300 K was repeated for about 1.5 ns while the system was compressed to 50 000 atm and then depressurized to 1 atm. In RadonPy, the amorphous system was considered to be in equilibrium when it satisfies the following conditions: the fluctuations in the total, kinetic, bonding, bond angle, dihedral, van der Waals, and long-range Coulomb energies with relative standard deviations (RSDs) of less than 0.05%, 0.05%, 0.1%, 0.1%, 0.2%, 0.2%, and 0.1%, respectively. At the same time, the RSDs for the fluctuations in density and the radius of gyration were less than 0.1% and 1%, respectively. NPT simulation was run at 300 K and 1 atm with a time step of 1 fs. The simulated system was checked for equilibrium states every 50 ns until the equilibrium requirements were realized.¹⁴¹ 

Polymer descriptors translate structural and chemical information about polymers into a machine-readable form for ML model training.¹⁷⁴ Successful descriptors are required to uniquely, completely, and minimally express polymer information and are a prerequisite and important condition for guaranteeing high accuracy of ML models. To date, many polymer representations have emerged to be exploited for polymer informatics, which can be categorized into three categories:¹⁷⁵ string-based descriptors, graph-based descriptors, and physical-based descriptors, as illustrated in Fig. 3.

String-based descriptors are efficient and convenient line notations, such as SMILES^158,176 and self-referencing embedded strings (SELFIES).^177,178 SMILES is a popular polymer representation that uniquely encodes atoms, bonds, rings, and branches in polymer monomers by ASCII string. SMILES allows a relatively uniform expression of polymers, does not depend on a strong chemical background, and has good readability for both humans and computers. Thus, SMILES is widely used in data storage,¹⁴² modeling,^{141,153,154,179} and ML^174,180 of polymers and has become a standard tool in computational chemistry. However, along with the progress of polymer informatics, SMILES has also exposed some limitations. For instance, in some inverse design tasks based on evolutionary or deep learning algorithms, a lot of the generated SMILES cannot correspond to valid polymers.¹⁸¹ An improved version of SELFIES⁹⁶ was proposed to address this issue, which can represent every polymer and be directly applied in arbitrary ML models without the adaptation of the models. Moreover, SMILES cannot usually be directly fed into the regression models and need to be further transformed through one-hot encoding¹¹⁰ or chemical language models.^182–184 

Graph-based descriptors are organic chemical representations based on topological information, involving substructure statistics, interatomic connections, and relative positional relationships. Molecular access system (MACCS) keys¹⁸⁵ are one of the most commonly used structural keys and have been integrated into some open-source cheminformatics software, including RDKit,¹⁴⁸ CDK,¹⁸⁶ and OpenBabel.¹⁸⁷ The generated MACCS keys contain 167 bits, the first bit determines whether the molecule has a predefined feature (if exists is set to 1, else is set to 0), and the last 166 bits correspond to 166 substructures. Morgan fingerprints, also known as extended connectivity fingerprints, are adapted from the Morgan algorithm and are one of the most popular fingerprints in chemical information.^188–190 The generation of Morgan fingerprints of a polymer monomer requires three steps.¹⁹¹ (1) Initialization: initializes each atom to be encoded as a unique integer identifier. (2) Iteration: in each iteration, the identifier of each atom is updated to a combination of its own and its neighbors' identifiers. The emerging identifiers are hashed to yield a fixed-length bit vector. Once all atoms have been given new identifiers, the old ones are replaced and the new identifiers are (3) post-processing: after a prespecified number of iterations, duplicate atom identifiers are removed, and the Morgan fingerprints are formed by the retained identifiers. Mogan fingerprints for polymer representation have the advantages of high efficiency, convenience, and absence of a pre-training process, but their feature dimensions are very high and sparse, and possibly introduce bit collisions caused by the hashing process. Mol2vec¹⁹² is an unsupervised method inspired by natural language processing techniques, that considers compound substructures derived from the Morgan algorithm as words and compounds as sentences. Ma et al.^120,193 trained the Mol2vec model on more than 1 × 10⁶ monomer structures from a combination of PoLyInfo and PI1M databases and obtained a pre-trained model called polymer embedding. Polymer embedding is a 300-dimensional continuous-valued vector and has been successfully used for the predictions of polymer density, melting temperature, glass transition temperature, and TC. Morgan fingerprint with frequency (MFF)¹¹⁶ is an expansion of Morgan fingerprints, which captures the frequency of chemical moieties (substructures) present in monomers. Across the entire explored chemical space, substructures (identifiers) with a frequency of occurrences larger than a predefined threshold are preserved to compose descriptor vectors. MFF has much lower feature dimensions compared to Morgan fingerprints, which can effectively suppress the overfitting of ML models and is widely utilized in the prediction of various properties of polymers.^114,115,194

Exploring the collection of physically independent descriptors for characterizing molecular structures is important for qualitative structure–property relationship building and provides a more intuitive guide to molecular performance assessment.¹⁹⁵ Benefiting from the advances in the feature engineering of drug-like molecules, some chemoinformatics software^196–198 is available for the automatic calculation of molecular descriptors. For example, Mordred¹⁹⁸ is a mainstream descriptor-calculation software that can calculate more than 1800 descriptors, including constructional descriptors, 2D topological descriptors, and 3D geometric descriptors.¹⁷⁴ However, polymer monomers are different from small molecules due to the presence of connection site information, which prevents the creation of 3D structures and leads to the inability to obtain some geometric descriptors. To compensate for the lack of 3D information, we incorporated molecular force field (FF) parameters as added descriptors and designed automated physical feature engineering for polymer descriptor optimization in our previous work.^85,124 The initial descriptors are a joint collection, which are calculated by the Mordred software and extracted from the parameters of the polymer structural data file after the FF assignment. Afterward, feature down-selection technology is employed to acquire the optimized descriptors, which consists of three stages.

1. Evaluate the numerical fluctuations of each physical descriptor itself using a variance indicator and remove features with low variance, since these features have less impact on target properties. This is beneficial in reducing the complexity and improving the performance of ML models. For a specific physical descriptor, the variance can be denoted as

   $S 2 = 1 n u m − 1 ∑ i = 1 n u m ⁡ ( x i − x ¯ ) 2 ,$

   (12)

   where $n$ is the number of candidates, $x i$ is the value of the feature for each candidate, and $x ¯$ is the average of the feature values.
2. Filtering and removal of features with poor correlation to target property by four correlation coefficients of the four metrics of Pearson, Spearman, distance, and maximum information coefficients (MIC). For two variables x and y, the Pearson coefficient can be solved by

   $P e a x , y = ∑ i = 1 n u m ⁡ ( x i − μ x ) ( y i − μ y ) ( n u m − 1 ) S x S y ,$

   (13)

   $μ k = 1 n u m ∑ i = 1 n u m ⁡ k ,$

   (14)

   $S k = 1 n u m − 1 ∑ i = 1 n u m ⁡ ( k i − k ¯ ) 2 ,$

   (15)

   where $μ x$ and $μ y$ are the mean values of the two variables x and y, respectively. $S x$ and $S y$ are the standard deviations of the two variables x and y, respectively. k can indicate any of the variables x and y.
   The Spearman coefficient is defined as

   $r s = 1 − 6 ∑ i = 1 n u m ⁡ d i 2 n u m ( n u m 2 − 1 ) ,$

   (16)

   where $d i$ is the difference in the ranks of the corresponding variables, i.e., the difference in the positions (ranks) of pairs of variables after the two variables have been individually ranked. In addition, the distance correlation coefficient¹⁹⁹ and MIC²⁰⁰ are used to measure the correlation of possible nonlinear variables.

3. Feature optimization using ML models to better match model performance. Recursive feature elimination (RFE)²⁰¹ is a widely employed technique for feature selection, which selects features by removing features with smaller absolute weights while ensuring accuracy in repetitive training of the ML model.²⁰² 

Establishing the mapping from microstructures to TCs using interpretable ML contributes to the understanding of the intrinsic thermal transport mechanism and guides the design of novel promising structures. Shapley additive explanations (SHAP)²⁰³ is a powerful tool for explaining ML models to alleviate black-box challenges, which links the optimal credit allocation of the model's input features to local explanations. Besides, the symbolic regression (SR)^204–207 technique can construct analytical models of key physical parameters for TC prediction and intuitively assists in capturing the underlying physical correlations. Some open-source frameworks for SR, such as gplearn,²⁰⁵ PySR,²⁰⁸ and SISSO,²⁰⁹ facilitate the discovery of optimal combinations between features and arithmetic symbols, enabling the creation of interpretable explicit mathematical models. Moreover, these tools are user-friendly and provide helpful step-by-step guidance.^210–212 Here, we presented three case studies of interpretable TC prediction models constructed using the SHAP or SR approaches, respectively.

The training of the regression model started with 1051 polymers sourced from a computational database, and all candidates were labeled TCs by NEMD simulations in RadonPy.¹⁴¹ We then calculated 325 initial descriptors, of which 294 are Mordred-based descriptors and 31 are MD-based descriptors. The down-selection process is displayed in Fig. 4(a). The threshold for variance assessment was 0.10 and a total of 202 descriptors were reserved. A weight assignment mechanism was developed for filtering descriptors, each metric was assigned a factor of 0.25, and descriptors with a cumulative weighting factor of 1 were retained. That is, the descriptor is valid only if all of the four correlation coefficients reach the corresponding thresholds. We filtered out 53 descriptors with a cumulative factor of 1 using the thresholds of 0.050, 0.050, 0.213, and 0.186 for Pearson, Spearman, Distance, and maximum information coefficients. Ultimately, the random forest (RF) model combined with RFE was employed in the Scikit-learn²¹³ for descriptor optimization, where 25 descriptors were determined based on the evaluation of the model's accuracy, as listed in Table II.

The 1051 polymers were represented by optimized physical descriptors and randomly split by training/test set as 80%/20%. We constructed an RF model using those optimized descriptors, where the hyperparameters for RF were optimized with Bayesian optimization with R² as a target.²¹⁴ The Gaussian regression process and acquisition function with ten random pairs of parameters were selected for initial training, and the ideal parameters for each ML model were determined after 100 optimization iterations.^215, Figure 4(b) compares the TCs predicted by RF with those calculated by NEMD simulations, with the training and test R² of 0.88 and 0.78, respectively. Moreover, we conducted 20 evaluations for RF models with optimized descriptors. During each evaluation, the training and test data were randomly sampled from a total of 1051 benchmark data at a ratio of 80/20%. The test R² of 20 RFs is 0.72 ± 0.02 (mean value and one standard deviation). The prediction error of the RF model in the high TC region (TC > 0.40 W m⁻¹ K⁻¹) is relatively large, since the percentage of these candidates among 1051 polymers is small (∼2.28%).¹²⁴ In addition, we verified the robust performance of the optimization descriptors and on other ML models such as multilayer perceptron and kernel ridge regression, and confirmed their superiority over graph descriptors of MACCS keys and Morgan fingerprints in our previous work.¹²⁴ 

The trained RF model was explained by the SHAP approach, and the importance of the features was analyzed based on the output SHAP values. Figure 4(c) exhibits the eight most important physical descriptors related to properties such as atomic number, atomic mass, bond connection, bond strength, sigma electrons, and monomer length. Combined with the distribution of SHAP values for the first eight physical descriptors of training candidates in Fig. 4(d), the promotion/inhibition of those key descriptors for TC can be recognized in general. The most significant descriptor is BCUTZ-1 h, which is the first highest eigenvalue of the Burden matrix weighted by atomic number and is associated with atomic and bonding properties.²¹⁶ BCUTZ-1 h was observed to have a strong positive correlation with the maximum atomic mass (Mass_max) in the monomer.¹²⁴ Typically, the presence of large masses of atoms such as chlorine, bromine, and iodine in the system suppresses lattice vibrations, resulting in small phonon group velocities and low TC. We also analyzed two MD-inspired descriptors of MW_ratio and K_bond_ave in detail, as shown in Figs. 4(e) and 4(f). The MW_ratio is the ratio of the molecular weight of the backbone to the total molecular weight of the repeating unit, and the K_bond_ave indicates the average of different bond force constants. Both descriptors have a positive effect on TC, as polymers with fewer side chains and stronger chain stiffness are favorable for efficient thermal transport.

The utilization of substructures of polymers as descriptors is more intuitive than physical descriptors in revealing structure–TC relationships and facilitates the design of new structures. We give a guide to building a deep neural network (DNN) model using 1144 polymer data with MFF. Among the 1144 polymers, the TCs of most of the structures were collected in a computational database,¹⁴¹ and the rest were calculated in Radonpy using the same setup.¹²⁴ Those polymers were expressed in the form of the SMILES and converted to MFFs. The MFF has 194 dimensions, mapped to the counts of the 194 most frequent substructures in the entire training data set. For DNN model training, the hyperparameters were optimized by Keras Tuner²¹⁷ Toolkit with Adam optimizer and mean squared error loss in TensorFlow,²¹⁸ using the 1144 polymers with a training/testing ratio of 80%/20%. After optimization of hyperparameters, the trained DNN model has four hidden layers with 416, 256, 244, and 256 nodes, respectively; ReLU activation; and dropout of 0.5.

The pairs of DNN-predicted and MD-calculated TCs are plotted in Fig. 5(a), with good consistency and a test root mean square error (RMSE) of 0.040 W m⁻¹ K⁻¹ and test R² of 0.79. We performed an additional fivefold cross-validation (CV) to evaluate the accuracy of the DNN model. In the fivefold CV, the test R² of the DNN models is 0.72 ± 0.03 (mean value and one standard deviation), reflecting that the trained DNN model is reliable.²¹⁹ Further, the trained DNN model was interpreted by SHAP, and the role of the most important 16 substructures on TC is shown in Fig. 5(b). When the descriptor values follow the same trend as the SHAP values, it suggests that the substructure is contributing to the realization of high TC. Thus, eight substructures were found to have a positive effect on TC, while six substructures suppressed TC, which are marked in Fig. 5(c) in blue and red text, respectively. The insights that can be extracted from these substructures coincide with the conclusions previously gained from the RF model trained with physical descriptors in Sec. III A, namely, that conjugated, linear side-chain-free polymers are beneficial for maintaining large chain stiffness and high TC. Additionally, when the polymer system contains heavy atoms such as fluorine (F), it hinders the efficient heat flux transport, preventing the achievement of high TC.

SR is another interpretable ML method for discovering specific mathematical expressions to match the fit of a data set.²⁰⁶ SR simultaneously searches for a set of parameters and the optimal mathematical formula of a function.²²⁰ Reliable training data are critical for SR without requiring massive amounts of data. We calculated the TC of 104 promising amorphous polymers recommended by the RF model¹²⁴ and found that their $R g$ extracted from equilibrium systems have a strong linear positive correlation with the TCs, as shown in Fig. 6(a). Besides the radius of gyration, we additionally acquired ten parameters from the equilibrium amorphous systems for SR, as listed in Table III. The SR was implemented in gplearn,²⁰⁵ and the hyperparameters settings are listed in Table IV. gplearn is a proven tool based on genetic programming that provides a scheme for optimizing mathematical expressions using genetic algorithms. Thus, the main input parameters of the genetic algorithm in gplearn include the optimization generations, the number of formulas produced in each generation, the crossover probability, and the probability of mutation at each node. The formulas were selected based on the criterion of simultaneously having low complexity and high fitting accuracy. The complexity was defined as the length of the formula, with each constant, variable, or operation symbol being represented as a unit of length. Figure 6(b) shows statistics of the 3364 mathematical formulae with complexity within 15 and R² over 0.70 by density plot. The six formulas F1–F6 at the Pareto front are marked by stars, and their corresponding analytic functions are listed in Table V. The Pareto front is the total set of all feasible and Pareto optimal solutions and is often considered the optimal trade-off between various objectives.^221,222 Overall, formulas with large complexity are more likely to yield large prediction accuracies. Of the six formulas, the smallest complexity is only 5, and the highest precision is 0.876, and their predicted TCs are both in good agreement with those calculated by NEMD, as depicted in Figs. 6(c) and 6(d). All six formulas capture the positive correlation between $R g$ and TC, and three of them reflect the fact that it is unfavorable to have large masses of atoms in the systems for TC. $R g$ was calculated to express the spatial extent of the molecular chains. When the molecular chains in the amorphous system have a high $R g$⁠, it is beneficial to maintain large rigid chain segments and enhance the heat transport along the chain backbone through intra-chain bonding interactions, thus increasing the TC.^60,124 In addition, it is worth mentioning that since some of the variables utilized for SR, such as $R g$⁠, were extracted from the equilibrated amorphous systems and are closely related to the TC, resulting in higher prediction accuracy of the obtained analytic mathematical models than that of the RF model in Fig. 4(b).

Active learning is oriented to the design of new polymers driven by target properties, which breaks through the limitations of regression tasks restricted to a fixed exploration chemical space. The inverse design of polymers with high TC can be achieved by some lightweight and smart optimization algorithms, such as genetic algorithm, Bayesian optimization, and quantum annealing. In this section, we not only present some cases of polymer design with a single target of high TC, but also additionally consider the synthesizability of polymers in multi-objective optimization trials.

Inspired by the knowledge of the interpretable DNN model outcomes in Sec. III B and the seen high TC polymers, we constructed a library containing 32 polymer motifs, as listed in Table VI. To ensure the uniqueness of the identification for each fragment, these motifs were binary coded from [00000] to [11111]. Theoretically, we can construct an infinite number of multi-block polymers by controlling the number and order of fragments. However, in consideration of the computational cost and hardware capabilities, we built a benchmark data set of triblock polymers. Figure 7(a) depicts the process of producing a triblock polymer, which was characterized by a binary sequence of 15 bits in length. The 15-bit binary sequence was divided equally into three equal parts, each one corresponding to a fragment. The composition of the polymer fragments is directionless, i.e., two polymers consisting of the fragments [0, 2, 29] and [29, 2, 0] are equivalent. The entire benchmark data set contains 16 896 triblock polymers, which were classified into 13 categories referring to the same classification method as PoLyInfo, such as polyolefins, polyethers, and polyimides [in Fig. 7(c)]. The TC of emerging polymers was estimated by a high-fidelity DNN model trained in Sec. III B and ranged from 0.16 to 1.03 W m⁻¹ K⁻¹, with 42.6% of the TC greater than 0.40 W m⁻¹ K⁻¹ [in Fig. 7(d)]. Moreover, the synthesizability of these polymers was evaluated by the synthetic accessibility (SA) score.²²³ The SA score, which ranges from 1 (easy) to 10 (hard), is calculated by considering both fragment contribution and complexity penalties. This score is used to evaluate the synthesizability of molecules or polymer repeating units. The SA scores for polymers are within a range of 2.28–6.21, with 6.3% of them having scores less than 3.0. Out of all 16 896 generated polymers, only 4.5% meet the predefined criteria for both TC and SA, which are known as ideal polymers.

The genetic algorithm (GA) is a heuristic search algorithm that is inspired by the biological evolutionary process and based on the mechanics of natural genetics and selection.^224–226 The operation of the GA algorithm involves various stages of population initialization, fitness evaluation, selection, crossover, and mutation. The core operators of crossover and mutation are illustrated in Fig. 8(a). Once the parents are selected according to the fitness function, the crossover operator combines parents into one or several offspring. After that, the mutation will execute with a predefined probability to increase the diversity in the population.

In our case, the GA was realized in the pymoo²²⁷ package with simulated binary crossover and polynomial mutation.²²⁸  Figure 8(b) depicts the convergence curve of the GA in a single optimization run, using the setups of ten initial random structures, 200 iterations × 10 candidates per batch, and both crossover and mutation probabilities of 100%. Thanks to the inheritance of excellent fragments from the parents, GA can efficiently explore the optimized polymers by only simulating a few candidates (64 non-repeating structures). To probe the effect of the initial structures on the convergence capacity, 20 GAs with different initial structures were executed, and the results are exposed in Fig. 8(c). The TC of the best structures for 20 runs ranges from 0.72 to 1.03 W m⁻¹ K⁻¹, within the top 0.3% of all possible candidates. A total of 6 rounds out of 20 optimizations yielded polymers with the global optimal TC, which indicates that the initial population has a significant influence on the optimization performance of the GA.

Bayesian optimization (BO) is a method for finding the optimal solution to black-box functions by using sequential design strategies that rely on the probabilistic surrogate models and acquisition functions.^229,230 Previously, we have released a tutorial on the design of thermal functional materials by coupling thermal transport calculations and BO.²³¹ Moreover, additional instruction is available to assist in understanding the core components of BO more intuitively.²³² The general process of BO can be described as²²⁵ (1) training a surrogate probabilistic model using several structures with the labeled property; (2) fitting the function using the surrogate model and recommending new structures by the acquisition function; (3) observing the property of the emerging structures and adding them to the training data set; and (4) repeating the above process until the predefined iteration times are reached.

As shown in Fig. 9(a), the single-task Bayesian optimization with a Gaussian process regression model and the Monte Carlo acquisition function of qEI was implemented in BoTorch.²³³  Figures 9(b) and 9(c) count the convergence curves of TCs for single and 20 parallel BO runs, using the hyperparameters set of initial random structures, 200 iterations × 10 candidates per batch, respectively. Along with optimization iterations, BO does search for optimized structures with enhanced TCs. However, many candidates with low TC are invalidly selected due to the uncertainty of Monte Carlo sampling. The best candidates in 20 runs have TC changed from 0.77 to 1.03 W m⁻¹ K⁻¹, of which 14 runs reach the global optimal structure.

Figure 11(a) compares the optimization performance of the three algorithms based on the averaged TCs from 20 separate optimization runs at different random states, and the shadows correspond to a standard deviation. Since BO has a robust Gaussian kernel and acquisition function that comprehensively evaluates the TC and uncertainty of the candidates, it has the strongest global optimization capability and the best overall performance. GA is able to inherit some fragments that have a positive effect on TC; thus, it enables a rapid increase in the TC of the optimized structures in the early optimization stages, but its ability is confined to the initial populations. Limited to the accuracy of the FM-trained surrogate model and the hardware of Ising machines, the simulated QA performs worse. However, affected by the stochasticity and uncertainty of the Monte Carlo sampler, BO simulated significantly more structures (after de-weighting) in a single optimization run compared to GA and QA [Fig. 11(b)].

Two state-of-the-art multi-objective optimization algorithms of unified non-dominated sorting genetic algorithm III (U-NSGA-III) and q-noisy expected hypervolume improvement (qNEHVI) were employed for the design of triblock polymers with both high TC and synthetic possibility. As depicted in Fig. 12(a), the U-NSGA-III is a multi-objective evolutionary algorithm (MOEA) that is an updated version of NSGA-III.^240,241 It generalizes different dimensional objective problems by increasing the selection pressure by a scalar selection operator. U-NSGA-III was implemented in the pymoo²²⁷ and kept all hyperparameters with default values. qNEHVI is a multi-objective Bayesian optimization (MOBO) that extends the acquisition function of expected improvement to hypervolume (HV) as an objective. It evaluates samples collected by the Quasi-Monte Carlo (QMC) sampler from the model posterior and identifies the candidate with the largest objective value. qNEHVI was operated in BoTorch,²³³ with the base and raw sampling set at 256 and 128, respectively, to enhance computational efficiency.

Figures 12(b) and 12(c) illustrate the optimization trajectories for a single run of MOEA and MOBO with ten random initial structures and 200 iterations × 10 candidates per batch. Nine gray stars indicate the sites of global optimal polymers, while the polymer dots are color-coded according to the generations. The distribution of searched non-duplicated polymer structures in a MOBO run is much denser than those in a MOEA run. qNEHVI integrates HV into the expected improvement acquisition function as an objective to evaluate the randomized QMC samples sourced from the model posterior, generating non-duplicated candidates in almost every generation. This enables the models to break out of local optimal solutions and further increases HV. In contrast, the optimization strategy of U-NSGA-III is inspired by the behavior of genes in organisms that crossover and mutate during evolution. The optimal polymers are designed by randomly selecting parents for matching and introducing a tournament operator. However, the performance of U-NSGA-III is affected by the initial polymer structures, as the optimization process primarily accumulates previous polymer units with positive contributions, making it easy to become trapped in locally optimal solutions.

To obtain statistical outcomes, we performed 20 runs of the MOEA and MOBO algorithms with different initial candidates, respectively. The HV convergence curves are displayed in Figs. 13(a) and 13(b). The HVs of U-NSGA-III can rapidly rise to a certain level within 20 generations, but it is difficult to increase again in subsequent generations. However, three qNEHVI runs identified nine global optimal polymers within 200 generations, and almost all of the HVs experienced a secondary boost after reaching a certain level for the first time. This enhancement difference depends on the stochastic nature of QMC sampling.²²¹ All the HVs of optimization algorithms reach a referred value calculated by the five ideal global optimal Pareto polymers (TC > 0.4 W m⁻¹ K⁻¹ and SA < 3.0) and the referred point ([0,−10] for TC and SA), although the mean HV of MOBO is greater than that of MOEA [see Figs. 13(c) and 13(d)].

Figures 14(a) and 14(b) present the number of explored de-duplicate polymers in 20 MOEA and MOBO runs, respectively. The number of unique polymers generated per MOEA run is significantly lower than that of MOBO, with a mean value of approximately 77, which is less than 5.0% of the average value for MOBO. To address this issue, an effective approach is to design high TC polymers through multiple parallel MOEAs with different random states, thereby reducing the influence of initial structures. Furthermore, we calculated the TC of 20 MOEA-designed polymers (red dots) using NEMD in Fig. 14(c), which indeed enhances the Pareto front (marked by stars) formed with 1144 raw polymers (blue dots). The MD validation of polymers with top 3 TCs is displayed in Fig. 14(d). These polymers are all combinations of fragments of the benzene and 3,5-dihydroimidazo[4,5-f]benzimidazole, with conjugated structures and large backbone stiffness. Consequently, intra-chain interactions of bonds, angles, and dihedrals dominate the contribution to TC.

Over the past few years, data-driven informatics algorithms have contributed to a revolution in the materials development paradigm, greatly facilitating the design of polymers and enhancing our understanding of their underlying mechanisms. In this Tutorial, we discuss the basic principles and implementation of ML for the exploitation of high thermal conductivity polymers, covering polymer datasets, polymer modeling and TC calculation, feature engineering, as well as informatics algorithms. We begin by describing the construction of interpretable regression models via physical or graph descriptors, and reveal the mapping between polymer microstructures and TCs. Based on the trained surrogate prediction model and the knowledge derived from the ML, we create a library containing 32 motifs and employ lightweight active learning algorithms to design sequence-ordered triblock polymers with high TCs. We not only focus on designing polymers with a single optimization target of high TC using GA, BO and simulated QA, but also consider the synthetic feasibility of polymers in multi-objective optimization trials that are realized by two state-of-the-art algorithms of U-NSGA-III and qNEHVI, respectively.

Although ML has facilitated the development of macromolecules with high TC, there is still a large gap in satisfying the various demands of realistic engineering applications, which also provides great opportunities for future investigations.

1. Sufficient high-quality polymer data are a fundamental prerequisite for polymer informatics. Accessible polymer databases are rare compared to inorganic databases, and the recorded data are rather sparse with strict acquisition rules. The development and preservation of publicly accessible polymer databases that adhere to the FAIR principles and encompass a wide range of properties that necessitate collaboration and consensus among chemical researchers. In addition, providing APIs for automated batch downloading of data is favored by polymer informatics.

2. Several open-source software^141,157 enable automated modeling and TC calculations of polymers through classical MD simulation, which promotes the development of polymer informatics. It is crucial to ensure the reliability of the obtained TCs of polymers. Therefore, efforts are being made at the computational level to enhance the accuracy of force fields for MD or to develop first-principles computational methods to be efficiently and economically applicable to the simulation of macromolecular systems.

3. Most current work on ML in polymer science focuses on the computational TC given the convenience and consistency of the research. Subsequently, with the success of applying ML to automated chemistry experiments,^242–244 we look forward to the emergence of automated platforms that integrate polymer literature mining, polymer synthesis, TC measurement, data storage and analysis, as well as novel structure generation and TC evaluation. This will enable the expansion of reliable polymer experimental data and the identification of promising polymers with high TCs.

4. State-of-the-art informatics algorithms are always sought after by the polymer community. Deep learning algorithms such as transfer learning,⁹⁴ recurrent neural networks, and reinforcement learning¹²¹ have been successfully applied to the exploitation of polymers with high TC. On the one hand, more intelligent molecular generation algorithms are required for the design of polymers with high TC; on the other hand, efforts are made to explore the application of large language models and multitask learning to the design of multifunctional polymers with enhanced TC.

Last but not the least, with the rapid advancement of artificial intelligence and automated experiments, we foresee that ML will become a powerful driving force in accelerating the design of advanced polymers to meet the immense demand in various fields. The attractive characteristics of polymers for ML extend beyond TC to encompass other properties such as optical, electrical, and mechanical properties.

This work was supported by Shanghai Key Fundamental Research Grant (No. 21JC1403300), Shanghai Pujiang Program (No. 20PJ1407500), the National Natural Science Foundation of China (NSFC) (Nos. 92366203 and 52006134.), and the SJTU Global Strategic Partnership Fund (2022 SJTU-Warwick).

The authors have no conflicts to disclose.

Xiang Huang: Investigation (lead); Methodology (lead); Writing – original draft (lead). Shenghong Ju: Conceptualization (lead); Funding acquisition (lead); Investigation (equal); Resources (lead); Supervision (lead); Writing – review & editing (lead).

Implementation of polymer physical feature engineering and the creation of interpretable machine learning models are available at https://github.com/SJTU-MI/APFEforPI, and the polymer inverse design cases are available at https://github.com/SJTU-MI/Inverse_Design_of_Polymers. More details can be found in our previous publications^124,219 or on reasonable request from the corresponding author.