Accurate prediction of structures in chemical compound space (CCS) is crucial for the design and synthesis of novel materials [1]. Although theoretical structure prediction methods have been widely adopted in material research [2,3], the computational cost associated with large-scale structural systems remains a challenge. A CALYPSO software based on swarm intelligence has been developed to accelerate structure prediction [4,5]. The existing first-principles-based (DFT) methods are time-consuming and economically expensive when applied to predict large-scale structures. Consequently, there is an urgent need to accelerate the process of material structure prediction [6].

Typically, structure prediction methods involve several steps [6]: (1) Randomly generating an initial population of structures under symmetry constraints; (2) Assessing the similarity among generated structures and eliminating duplicates; (3) Performing local optimization and energy calculations using the DFT on generated structures; (4) Constructing the next generation of structures using swarm intelligence algorithms, where low-energy structures evolve into new structures while high-energy structures are replaced by randomly generated ones. However, for large-scale structures containing hundreds or thousands of atoms, the computational time required for energy calculations is significant. Furthermore, performing energy calculations and optimizations for multiple randomly generated structures further compounds the time cost. This limitation significantly impedes the prediction of large-scale material structures in most structure prediction software. Deep learning (DL) methods have gained substantial success in computer vision, pattern recognition, and other fields [7,8]. However, their application in predicting and generating new structures in computational chemistry remains limited [9].

In conclusion, the acceleration of material structure prediction is a critical problem that needs to be addressed. While traditional methods face challenges in dealing with large-scale structures, DL methods offer a potential solution. By leveraging DL techniques, it may be possible to overcome the time constraints associated with energy calculations and optimize the prediction process. This application of DL in computational chemistry holds promise for advancing the field of structure prediction and generating new structures.

In recent years, the rapid development of artificial intelligence has made ML a new method that researchers are eager to explore. ML has become popular in various application scenarios that require prediction, classification, and decision-making [4,10]. With the availability of large-scale quantum mechanical data [11–13], researchers have established ML models and used them to predict material properties such as formation energies, defect energies, elasticity, and other mechanical properties [4]. For instance, Hansen et al. [14] employed a linear regression (LR) algorithm to learn the relationship between structural information and cluster energy and predicted the energy of new clusters based on atomic Cartesian coordinates and nuclear charge. Meanwhile, neural network (NN) methods [15,16] have been leveraged to accelerate the construction of potential energy surfaces. Rupp et al. [17] introduced a machine learning model that predicts the AEs of different organic molecules. Gupta et al. [18,19] utilized the DL method to establish the correlation between molecular topological features and accurate AE predictions, achieving an impressive MAE.

In addition, some DL methods were used to replace the DFT calculation process to reduce the computational complexity of the system and accelerate structure prediction in an end-to-end way [20]. Due to the strong feature extraction and feature learning abilities of the DL methods, we have introduced a new NPN [21] model that predicts the AE of different molecules in an end-to-end manner based solely on nuclear charge and atomic position. To address this task, an end-to-end energy prediction model was constructed, and its schematic diagram is shown in Fig. 1.

Let Φ be an AE prediction model, and the model Φ_NPN utilizes a neural network architecture for probabilistic modeling. The model leverages the properties of exponential family distributions to transform the neural network output into the natural parameters of the exponential family distribution. This allows NPN to conveniently parameterize the probabilistic model and maximize the log-likelihood function using gradient-based optimization algorithms during training. NPN transforms the conditional probability p(Y/X) between the input variable and the output variable into the form of an exponential family distribution. According to all the molecular structure information, we extract the coordinate information of all atoms of each molecule as the input of Φ_NPN, which is represented by X. The AE of each molecule is represented by Y. Here we have:

(1)

where X is the coordinate information of all atoms of each molecule, and θ_NPN is the parameters of the NPN model. NPN consists of both linear and nonlinear transformation layers. Here, we first introduce the linear form of NPN. For simplicity, let's assume a Gaussian distribution ξ = (q, d)^T with two natural parameters. We decompose the distribution over a weight matrix p(W^l|W_q^l, W_d^l) = Π_i,jp(W_ij^l|W_q,ij^l, W_d,ij^l), where (W_q,ij^l, W_d,ij^l) is the corresponding natural parameter. We assume similar decomposed distributions for b^l, z^l, and x^l in NPN, which are all exponential family distributions. For computational convenience, we wish to approximate z^l using another exponential family distribution by matching its mean and variance. To compute (W_m^l, W_s^l) = f(W_q^l, W_d^l) and (b_m^l, b_s^l) = f(b_q^l, b_d^l), we can obtain z_q^l and z_d^l from the mean z_m^l and variance z_s^l of z^l as follows:

(2)

(3)

(4)

where, the symbol ⊙ represents the product of the elements in a distribution, while the bijective function f(., .) maps its natural parameters to the distribution's mean and variance. Similarly, symbol f⁻¹(., .) denotes the inverse transformation. The means and variances of W^l and b^l, which are obtained from the natural parameters, are represented by W_m^l, W_s^l, b_m^l and b_s^l, respectively. The values of z_m^l and z_s^l that are computed can be used to recover z_q^l and z_d^l, thereby streamlining the feed-forward calculation of the nonlinear transformation.

Once we have the distribution over z^l, defined by natural parameters z_q^l and z_d^l, that has been linearly transformed, we then apply an element-wise, nonlinear transformation v(·) to it. The resulting distribution of activations is denoted as p_x(x^l) = p_z(v⁻¹(x^l))|v^−1′(x^l)|, where p_z is the factorized distribution over z^l defined by (z_q^l, z_d^l). Even if p_x(x^l) is not an exponential-family distribution, we can still approximate it using p(x^l|x_q^l, x_d^l) by matching the first two moments. Once we obtain the mean x_m^l and variance x_s^l of p_x(x^l), we can compute the corresponding natural parameters using f⁻¹(., .). The feed-forward computation can be expressed as:

(5)

(6)

To model distributions over weights and neurons, the natural parameters must be learned. The NPN is designed to take a vector random distribution as input, such as a multivariate Gaussian distribution. It then multiplies this input by a matrix random distribution and applies a nonlinear transformation before outputting another distribution. Since all three distributions in this process can be characterized by their natural parameters, learning and prediction of the network can be performed in the space of natural parameters. During back-propagation, for distributions characterized by two natural parameters, the gradient is composed of two terms. For instance, in the equation ∂E/∂z_q = (∂E/∂x_m)⊙(∂x_m/∂z_c) + (∂E/∂x_s)⊙(∂x_s/∂z_q), where E represents the error term of the network. Naturally, nonlinear NPN layers can be stacked to form deep NPN, as shown in Algorithm 1. NPN does not need expensive reasoning algorithms, such as variational reasoning or Markov chain Monte Carlo (MCMC). Moreover, in terms of flexibility, it can choose different types of exponential family distributions for weights and neurons.

We know that χ can be converted using a previous study [22] with Coulomb matrices (CM). The CM are matrices that contain information about the atomic charges and atomic coordinates. The CM can be calculated using the following formula:

(7)

among them, 1≤i, j≤23, Z_i is the nuclear charge of atom i, and R_i is its Cartesian coordinate.

Specifically, QM7 is a dataset consisting of 7165 molecules taken from the larger GDB-13 database, which comprises nearly one billion synthetic organic molecules. QM7 has information about the molecule that consists of the atoms H, C, N, O, and S. The maximum number of atoms in a molecule is 23. The molecules in QM7 exhibit a diverse range of chemical structures, such as double and triple bonds, cycles, carboxylic acids, cyanides, amides, alcohols, and epoxides [23]. The CM representation is used to describe the molecular structures, and this approach is invariant to translations and rotations of the molecule. In addition to the CM, the dataset includes AEs ranging from −800 kcal/mol to −2000 kcal/mol [14].

The main focus of this article is on the χ and components in the QM7 dataset. Due to the varying number of atoms in each molecule, the sizes of their true CMs are also different. To make all CMs with the same size of 23 × 23, we extended the size of the original CM and filled the expanded space in the CM with 0.

For NPN, GCN, CNN, and LSTM models, we use the rectified linear unit (ReLU) activation function in the hidden layer. We used the ADAM optimizer [24] with a learning rate of 0.0001 and a decay rate of 0.00001, and based on our experience and knowledge, we have selected different combinations of network hyper-parameters for model training and fine-tuning. We set the learning rate to [0.01 0.001 0.0001], the number of training rounds [100 200 300], the batch size [32 64 128], and the hidden layer nodes [100 200 300]. Then, we combined the values of different network hyper-parameters. The experimental results showed that when the learning rate was 0.0001, the batch size was 64, the hidden layer nodes were 200, and the number of training rounds was 200, the model performance was optimal. The proposed models have been implemented in the Pytorch library. As for the hardware system configuration, the processor is an Intel (R) Xeon (R) Silver 4110 CPU@2.10 GHz, the RAM has a capacity of 32.0 GB, and the graphics card is an NVIDIA RTX 2080 Ti [25]. We randomly divided the dataset into a training dataset and a test dataset in the ratio of 8:2.

We used the proposed method to train the AE prediction models on training data and validated the performance of the models using testing data, and the losses of the proposed models during the training are shown in Fig. 2.

By training on a large set of different molecules in QM7, we use MAE and root mean square error (RMSE) as evaluation indicators to evaluate the prediction performance of various algorithms under different molecule presentation methods. The comparison results between our method and other methods [22] are shown in Table 1. Our results indicate that the proposed method exhibits a lower prediction error, indicating excellent prediction performance. Based on this, we trained and tested different DL models we designed on the QM7 dataset [26], and obtained some information that reflects the performance of the models, such as model parameters, testing time, and testing errors. Experiments on real-world datasets [26] show that NPN can achieve state-of-the-art performance on regression tasks. The performance of different methods are compared in Table 2.

In previous experiments [14], researchers often merged feature values/vectors and feature centers with flattened CM to form a 7165-sample dataset with N-dimensional features. However, the integration of new molecular property information has brought about a new problem: the concatenation of old and new features may lead to unwanted "heterogeneity" in the feature vectors, and including more features may actually increase noise in the dataset [30]. It should be noted that some traditional ML techniques may not be able to identify meaningful patterns from these newly added features, so merging new features may result in poorer results [31].

We present the results analysis of AEs prediction for four different molecules and their isomers, it can be shown in Fig. 3. Based on the structural information, our proposed method achieved good AE prediction results by using the extended CM method to characterize molecules. However, we also found that the model has a significant prediction error for molecules [32] with fewer atoms, which may be due to the interference of additional information supplemented in the extended CM to the model's learning of structural information features. For some molecules with compact structures relative to sparse spatial structures, the prediction errors were relatively larger. We analyzed that this may be due to the high density distribution of a large number of identical atoms in the molecule, which led to the model's unsatisfactory extraction of the CM features.

To validate the robustness and effectiveness of our proposed method, we randomly divided the BC2P data from Fu et al. [33] into a train dataset and a test dataset with an 8:2 ratio. We use our proposed method to train a model on train data, and the trained model accurately predicts the corresponding energy [33]. The training loss curve can be seen in Fig. 4. Moreover, we tested the model on test data, and the results showed that our model had an average prediction time of 0.391 ms for the energy of each molecular conformation and an average test loss of 0.236 eV/atom. The results of performance evaluation of NPN model on BC2P is shown in Table 2.

The study aims to propose and validate a method for predicting atomic energy. Initially, we introduce the extended CM representation of molecular structure along with its inherent characteristics. Subsequently, we compare the performance of various ML and DL models for predicting AEs. These models not only enable fast AE prediction but also expedite the structure prediction process. Furthermore, to overcome the limitations of our current method, we plan to incorporate the three-dimensional spatial pattern of material structure information in future research. We intend to employ the graph convolutional network (GCN) method to mathematically represent atoms in chemical molecules. In addition to model development, we informally discuss the results of numerical experiments conducted in this study. Through extensive comparative experiments and analyses, we gain valuable insights into the performance and potential applications of the proposed method. Some potential applications include catalyst design, materials discovery, optimization of energy storage and conversion, as well as material performance prediction. In the future, we envision conducting further in-depth research on material performance prediction and material discovery. By leveraging the power of machine learning and data mining techniques on large datasets, we can uncover hidden correlations within complex data structures. This will allow us to predict new material structures and their corresponding properties. Ultimately, this research contributes to the advancement of materials science and opens up new possibilities for designing innovative materials with tailored properties.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment

This work was supported by the Nature Science Foundation of China (Nos. 61671362 and 62071366).