Exploration of elastic moduli of molecular crystals via database screening by pretrained neural network potential

This paper presents the database screening of elastic moduli of molecular crystals using a pretrained neural network potential. The calculated elastic moduli were sufficiently consistent with experiment values, better than Hartree-Fock calculations. The database screening suggested crystals with large and small moduli. The elastic constant tensor, commonly referred to as the stiffness tensor, serves for characterizing the elastic behavior of materials. Once the components of the stiffness tensor are ascertained, mechanical properties including the bulk modulus, Young's modulus along various crystallographic axes, the Poisson’s ratio, and material anisotropy can be computed. The elastic modulus, a measure of a material's resistance to elastic deformation, is a crucial parameter for assessing its suitability for specific applications. For instance, engineered structures subjected to high mechanical stress necessitate materials with substantial elastic moduli, typically ranging from tens to hundreds of GPa. 1 Conversely, materials intended for biocompatible applications are often polymeric in nature, exhibiting elastic moduli on the order of several to hundreds of MPa. 2 Molecular crystals generally have medium Young's modulus between hard and soft materials. 3-5 Molecular crystals have not been the mainstream of industrial materials, but they are expected to be used in optoelectronics and actuators, 6-8 and so it is significant to clarify their elastic properties. Although there have been many reports in which Young's modulus in a specific direction has been measured by nanoindentation or bending tests, 9-11 there have been a limited number of literatures in which the elastic constant tensor has been experimentally determined. Spackman et al. reported that only 80 crystals have been determined their elastic constant tensors since 1950. 12 This is due to the difficulty of obtaining large-size crystals suitable for measurement and the need to determine many tensor components (21 components at maximum) due to their low symmetry compared to inorganic crystals. If the elastic constant

The elastic constant tensor, commonly referred to as the stiffness tensor, serves for characterizing the elastic behavior of materials.Once the components of the stiffness tensor are ascertained, mechanical properties including the bulk modulus, Young's modulus along various crystallographic axes, the Poisson's ratio, and material anisotropy can be computed.The elastic modulus, a measure of a material's resistance to elastic deformation, is a crucial parameter for assessing its suitability for specific applications.4][5] Molecular crystals have not been the mainstream of industrial materials, but they are expected to be used in optoelectronics and actuators, [6][7][8] and so it is significant to clarify their elastic properties.Although there have been many reports in which Young's modulus in a specific direction has been measured by nanoindentation or bending tests, [9][10][11] there have been a limited number of literatures in which the elastic constant tensor has been experimentally determined.Spackman et al. reported that only 80 crystals have been determined their elastic constant tensors since 1950. 12his is due to the difficulty of obtaining large-size crystals suitable for measurement and the need to determine many tensor components (21 components at maximum) due to their low symmetry compared to inorganic crystals.If the elastic constant tensors of molecular crystals can be calculated, it may be possible to develop molecular crystals with desired elastic properties.The material screening should be effective approach to develop functional molecular crystals.
4][15][16] DFT-D calculations offer a notable advantage in calculation accuracy, of course depending on the selected basis and method.However, they are computationally intensive, making it challenging to screen large datasets efficiently.Recognizing this limitation, Spackman et al. employed the corrected small basis set Hartree-Fock method (S-HF-3c) for elasticity calculations, 12 which has reduced computational cost in comparison to DFT-D.Their study provided a comprehensive comparison of experimentally determined elastic constant tensors with those obtained from S-HF-3c calculations, revealing that the elastic moduli determined by S-HF-3c were broadly close to experimental values. 12The importance of their work lies not only in highlighting the efficacy of S-HF-3c but also in establishing a publicly accessible benchmark dataset (elasticity dataset) comprising both experimental and calculated elastic tensors.
8][19] Neural network potentials (NNPs) can approximate the relationship between material structures and potential energies with an accuracy close to that of DFT-D, provided sufficient training data is available.][25][26][27] These NNPs have been used for screening and molecular dynamics simulation, mostly for inorganic crystals and metalorganic frameworks.However, the use of NNP for predicting the properties of molecular crystals is limited, and its effectiveness remains unclear.If NNP proves to be effective, it can be utilized for the screening of molecular crystals.
In this research, the prediction accuracy of NNP for the elastic constant tensor of molecular crystals was assessed using the elasticity benchmark dataset. 12The pretrained NNPs available from the Atomic Simulation Environment (ASE) calculator, namely PreFerred Potential (PFP) based on TeaNet and the Crystal Hamiltonian Graph Neural Network (CHGNet), were employed. 24,27Their predictive capabilities were first evaluated using metrics associated with structure relaxation and subsequently against experimental elastic moduli.A comparison of NNP results with experimental results and theoretical calculations (DFT-D and S-HF-3c) revealed that the NNP's prediction accuracy was more akin to DFT-D than S-HF-3c.Moreover, this study inferred various moduli from a large dataset, highlighting molecular crystals with both larger and smaller moduli owing to low computation cost of NNP.
The benchmark dataset comprises 123 experimental elastic tensors corresponding to 80 distinct compounds.S-HF-3c calculations have been carried out for 44 compounds and DFT-D results for 11 compounds have been summarized, in the literature. 12Comprehensive details can be found in the reference. 12For comparative purposes, NNP calculations were executed on the same dataset with S-HF-3c.When multiple experimental results were reported for a single compound, the average was considered as true value in the current study.
To remove residual strain prior to the calculation of elastic constant tensors, the optimization of atomic coordinates and cell parameters was conducted by pretrained PFP and CHGNet models using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method, 28 with a maximum of 2000 iterations and a residual force threshold at 0.02 eV/Å.The optimized crystal structure was assessed by the relative deviation ΔV of the cell volume from the experimental volume.
The experimental-optimized scatter plot of cell volume indicates that the PFP optimization outperformed others (Figure 1).CHGNet generally overestimated the volume.S-HF-3c tended to underestimate it as pointed out in previous literature. 12All the optimized cell parameters using PFP and CHGNet can be found in Figure S1 and S2.The mean absolute error (MAE) for ΔVPFP was 2.15% (Table 1), which is lower than the 6.5% error of ΔVS-HF-3c (Table S1).These findings validate that PFP calculates potential energy with greater precision compared to S-HF-3c.
In contrast, the MAE for ΔVCHGNet was 19.88% (Table 1).This considerable discrepancy is likely attributable to the training dataset.The pretrained CHGNet has been trained using crystal structures in Materials Project, which predominantly features inorganic crystals. 27Conversely, PFP has been trained using a broader dataset, encompassing both Materials Project and organic dataset. 24 further assess the NNPs, the X23 dataset-commonly employed to evaluate theoretical calculations of molecular crystals-was utilized.This dataset is comprised of crystals characterized by typical intermolecular interactions: van der Waals forces and hydrogen bonds.The MAEs of ΔVPFP and ΔVCHGNet for the X23 dataset were 4.27% and 25.19%, respectively (Table 1 and Table S2).The literature reported that the MAE of ΔV was 5.2±2.3% for S-HF-3c and 2.7±3.2% for DFT-D (PBEh-3c). 12This result also suggested the generalization ability of PFP.Moreover, structure optimizations were conducted on molecular crystals whose papers have been published in 2023 (Y2023 dataset) to guarantee that inference data is not used for training.The resultant MAE for ΔVPFP was 1.88%, closely mirroring the MAEs of ΔVPFP in the elasticity and X23 datasets (Table 1 and Table S3).Consequently, it has been determined that the PFP is suited for molecular crystals as pretrained neural network potential.Here, it is posited that the worse performance of CHGNet should not be inherent to its model architecture but is rather dependent on the training data.It is anticipated that CHGNet could yield better results if provided with enough training data of molecular crystals.The elastic constant tensor is determined by the secondrank derivative of the potential energy, and yields values for the bulk, Young's, and shear moduli, as well as the Poisson's ratio and anisotropy (see Supporting Methods).Elastic moduli can be compared using average properties, and multiple averaging schemes are available such as Voigt and Reuss schemes.The Hill averaging employs the arithmetic mean of the Voigt and Reuss averages.In addition, the arithmetic mean of the Reuss and Hill

Please do not adjust margins
Please do not adjust margins averages has been mentioned as a better approximation in the literature, 12 and used in this study.The elastic properties were predicted and subsequently compared with experimental values.The error in predicting elastic properties using PFP was found to be lower than with S-HF-3c, and comparable to the error associated with DFT-D (Table 2).For instance, the MAE of ERH calculated by PFP was 3.55 GPa, which is less than that of S-HF-3c (7.53 GPa) and even DFT-D (3.87 GPa).The observed-predicted plot of ERH indicated superior predictions using PFP (Figure 2a), with errors being distributed in an almost random fashion (Figure 2b).There is no evident bias in PFP unlike S-HF-3c, which remains overestimation bias (Figure 2b).Such a trend was consistent across different averaging approaches for E and G, surpassing the MAE of the mean model, which assumes that there is no relationship between the input structure and output.Thereby it indicates a successful prediction of E and G by PFP, which should capture of the underlying relationship between crystal structure and potential energy.Observed-predicted plots of other elastic properties were summarized in Figure S3-5.
The prediction of the bulk modulus K using PFP did not surpass the mean model but was still superior to S-HF-3c (Table 2).The observed-predicted plot demonstrates that the predicted values from PFP align closely with the reference line, though there seems to be a slight intercept, possibly attributed to some inherent bias (Figure S3).As the Poisson's ratio and anisotropy have relatively smaller values than elastic moduli, assessing prediction accuracy becomes challenging.However, the prediction errors from PFP align closely with other computational methods (Table 2).Since PFP afforded better prediction ability than S-HF-3c, the author screened 5000 crystal structures, which were downloaded from the Cambridge Structural Database (CSD), to identify molecular crystals with large and small Young's modulus.The deviation of cell volumes of the optimized structures from experiments was sufficiently small, 1.88% (Figure S6).The predicted ERH is distributed with a mean of 9.60 GPa, a mediun of 9.17 GPa, a maximum of 92.5 GPa, and a minimum of 1.88 GPa, respectively (Figure 3a).Other elastic properties are also summarized in Figure S7.Some crystal structures were visualized to find structural features from the screening.The largest ERH was found in WAWNAL (CSD refcode), which formed hydrogen-bonding chains and p-p interactions along the b-axis of 2-fold screw axis (Figure 3b).This result is agreed with the knowledge that hydrogen bonds promote the larger Young's modulus. 4The second crystal (IVENAZ) also forms hydrogen-bonded chain of C=O…H along the b-axis (Figure 3b).The third crystal (KUCWEN) is formed by ionic bonds between C2N 3-and Sr 2+ , which are stronger interactions than most intermolecular interactions of molecular crystals.Other high-ranking crystals also have similar interactions in common.The largest modulus E of molecular crystals is reported to be 46.8GPa, 2 and thus the suggested crystals may have larger modulus than known if measurement is conducted.
The smallest predictions are located at a left tail of the distribution centered around 10 GPa, resembling a normal distribution (Figure 3a).Molecular crystals often have weak intermolecular interactions, leading to a greater number of examples than larger Young's modulus.The smallest crystal (ECUSAZ) is formed by van der Waals interactions without evident stronger interactions (Figure 3c).The other smallest predictions are also found in the crystals without stronger interactions such as hydrogen bonds (Figure 3c).The smallest modulus of molecular crystal is known to be 0.27 GPa, 2 and no crystals smaller than the smallest known modulus could be

Figure 1 .
Figure1.Experimental-optimized plot of cell volumes.The data of the largest volume was not calculated using CHGNet due to computer memory limit.The dashed line is the reference line when predictions perfectly match with experimental data.

Figure 2 .
Figure 2. Comparison of experimental and calculated Young's modulus.(a) The observedpredicted plot of ERH.(b) Error plot of the predictions.Error was defined by calculated value minus experimental value.The dashed lines are the reference lines when predictions perfectly match with experimental data.

Table 2 .
MAE of predictions compared with experimental values.