Yuan Qi, Hong Ren, Hong Li, Ding-lin Zhang, Hong-qiang Cui, Jun-ben Weng, Guo-hui Li, Gui-yan Wang, Yan Li. Interaction Energy Prediction of Organic Molecules using Deep Tensor Neural Network[J]. Chinese Journal of Chemical Physics , 2021, 34(1): 112-124. DOI: 10.1063/1674-0068/cjcp2009163
Citation:
Yuan Qi, Hong Ren, Hong Li, Ding-lin Zhang, Hong-qiang Cui, Jun-ben Weng, Guo-hui Li, Gui-yan Wang, Yan Li. Interaction Energy Prediction of Organic Molecules using Deep Tensor Neural Network[J]. Chinese Journal of Chemical Physics , 2021, 34(1): 112-124. DOI: 10.1063/1674-0068/cjcp2009163
Yuan Qi, Hong Ren, Hong Li, Ding-lin Zhang, Hong-qiang Cui, Jun-ben Weng, Guo-hui Li, Gui-yan Wang, Yan Li. Interaction Energy Prediction of Organic Molecules using Deep Tensor Neural Network[J]. Chinese Journal of Chemical Physics , 2021, 34(1): 112-124. DOI: 10.1063/1674-0068/cjcp2009163
Citation:
Yuan Qi, Hong Ren, Hong Li, Ding-lin Zhang, Hong-qiang Cui, Jun-ben Weng, Guo-hui Li, Gui-yan Wang, Yan Li. Interaction Energy Prediction of Organic Molecules using Deep Tensor Neural Network[J]. Chinese Journal of Chemical Physics , 2021, 34(1): 112-124. DOI: 10.1063/1674-0068/cjcp2009163
The interaction energy of two molecules system plays a critical role in analyzing the interacting effect in molecular dynamic simulation. Since the limitation of quantum mechanics calculating resources, the interaction energy based on quantum mechanics can not be merged into molecular dynamic simulation for a long time scale. A deep learning framework, deep tensor neural network, is applied to predict the interaction energy of three organic related systems within the quantum mechanics level of accuracy. The geometric structure and atomic types of molecular conformation, as the data descriptors, are applied as the network inputs to predict the interaction energy in the system. The neural network is trained with the hierarchically generated conformations data set. The complex tensor hidden layers are simplified and trained in the optimization process. The predicted results of different molecular systems indicate that deep tensor neural network is capable to predict the interaction energy with 1 kcal/mol of the mean absolute error in a relatively short time. The prediction highly improves the efficiency of interaction energy calculation. The whole proposed framework provides new insights to introducing deep learning technology into the interaction energy calculation.
Molecular dynamics (MD) simulations can provide detailed information concerning individual particle motions as a function of time.Moreover, it can be applied to address the specific questions about the properties of physics, chemistry, biology, and material, especially the organic related systems.
Generally, the MD technology is more easily implemented than the experiment tools on most of the practical organic systems.Therefore, MD has become an essential tool for understanding the physical basis of the structure and function of organic molecules and its related systems [1-3].In MD simulation, the interaction between molecules depicts the motions of atoms and electron based on the Schrdinger equations, and it also depicts the influence between molecules.Therefore, the interaction energy is one of the most fundamental variable to simulate the atomic movements iteratively.
The classical methods of studying the interaction energy are mainly based on quantum mechanics (QM) and force-fields (FF).In one aspect, the QM, as a classical theoretical calculating framework, is hard to be entirely implemented on the whole MD simulation of large organic systems for its limitations of sheer calculating time and source.In the other aspect, the force field represents the molecular features with simplicity parameters. The springs for bond length and angles, the periodic functions for bond rotations and Lennard-Jones potentials, and the Coulomb's law for van der Waals and electrostatic interactions, respectively are also mature theoretical frameworks. When an organic system is built, the forces acting on each atom are obtained by a specific force field, where potential energy is deduced from the molecular structure.
The calculation of classical potential based on the force field is extremely fast, which enables the MD simulation of the organic system with 109 atoms most.However, the classical force field suffers from the lack of transferability and can yield the accurate results only close to the conformations they have been fitted to. Therefore, the parameterization process is empirical, and the resulting potentials may not transfer to the dynamic process of new organic systems.One feasible solution to the accuracy problem is implemented by ab initio molecular dynamics (AIMD) simulations, where the quantum-mechanical energy and forces are computed for each atom in the organic system at each iteration.In recent years, the QM and FF are combined as QM/MM to simulate the specific area with QM precision, and the left system is simulated via MM based on FF [4, 5].In FF improvements, there are already many mechanisms introduced into the MD simulations process. The polarizable force field, such as the Drude polarizable force field [6], AMOEBA [7-9], depicts the polarization more accurately and have already been applied to the organic related systems. Furthermore, the fluctating charge model [10, 11] and other charge related models [12, 13] are introduced to illustrate the polarized effect in certain process accurately.To reduce computational time, the CG model is applied [14].In the hardware aspect, the Anton machine [15, 16] is a special-purpose system for molecular dynamics simulations of proteins and other biological macromolecules.The widely accepted scheme is to use the graphics processing unit (GPU) into MD simulation [17-24] and parallelly calculate the biomolecular and organic systems in different GPU accerlerating cards.In QM improvements, many different models are proposed to fit various kinds of systems with high precision and reduce calculating source and time relatively.Moreover, the proposed calculating method and basis sets are already applied to drug discovery [25] and analyzing the conformation changing with electron transfer [26]. Though there exist many successful applications, it is still urgent to develop a theoretical method which predicts the interaction energy large organic system in QM precision.
In the last two decades, the computer hardware, especially the GPU, undergoes a high development, along with the machine learning technology, which has boomed rapidly.The deep learning (DL), as an important part of machine learning technology, with its ability to draw the abstract information from big data such as images, audios, and videos metals, has been applied into countless fields [27-32].
In the last few years, the machine learning including deep learning, has also been used into the computational chemistry researches.Machine learning, in general, and deep learning, in particular, are ideally suitable for representing quantum-mechanical, enabling the researchers to model nonlinear potential-energy surfaces and enhancing the exploration of chemical compound space.The machine learning methods have been implied as an effective tool to represent the physicochemical properties of molecular systems recently.Most of these studies have achieved remarkable success.These studies have demonstrated the potential of using machine learning methods and particularly deep learning models to represent the potential energy surface (PES).Montavon et al. applied a deep multi-task artificial neural network to predict multiple electronic ground and excited state properties of organic molecules with ab initio calculation accuracy [33].Behler analyzed the changing process of smaller molecular via high dimensional potential and achieved great achievements [34-38].The conformation structure calculation and machine learning technology have been combined aiming toward the target of accelerated analyzing of molecular properties.Rupp et al. applied machine learning technology to predict atomic energy [39].Schnet et al. [40] used a continuous-filter convolutional neural network to predict the interaction of the molecular system.Chmiela et al. developed spatial and temporal physical symmetries with a gradient-domain machine learning (sGDML) model to build a CCSD(T) level force field [41].Based on locality and active learning, Gubaev et al. predicted molecular properties [42].Inspired by the many-body expansion, Lubbers et al. introduced the hierarchically interacting particle neural network (HIP-NN) to predict the molecular properties from datasets of quantum calculations [43].Schütt et al. described an interpretation technique for atomistic neural networks on the example of Behler-Parrinello networks as well as the Schnet. Both models obtain predictions of chemical properties by aggregating atom-wise contributions [44].Zhang et al. proposed a framework, named DeepMD-kit to perform molecular dynamics based on deep learning potential energy and force field [45, 46].In 2018, Zhang et al. modified the DeepMD framework and proposed a new smooth version (DeepPot_SE) to predict the energy of the molecular system and force of atoms.Kearnes et al. proposed to encode the molecular as a molecular graph to fed a convolution network [47].Gastegger et al. applied weighted atom-centered symmetry function (wACSFs) as descriptors of the geometric structure of a chemical system [48], and predicted the chemical properties and potential energies via the high-dimensional neural network in machine learning.
The deep tensor neural network (DTNN) was first proposed to speech recognition [49].In 2017, Kristof and Farhad et al. simplified the connect of tensor layers of DTNN and applied this deep neural network model into the energy prediction of small molecules in GDB9 [50], a database contains more than 2×105 molecules with their physical and chemical properties.In 2018, we combined the DTNN and MFCC mechanisms to predict the energy of amino acid fragments of protein.And the prediction results showed a good correlation between the DTNN and QM calculation in B3LYP/6-31+G∗ precision [51].
In this work, to solve the accuracy and molecular size dilemma and furthermore to provide a possible way to enable the AIMD simulation into analyzing the organic related system, we extend our work and introduce it into the interaction energy calculation process.The DTNN network is chosen as the main frame work as the deep learning technology.Here, we used deep learning into a little more complex molecular system with two molecules, including at least one organic molecule.In this introduced DTNN framework, the output is no longer the energy of a single molecule, but being replaced by the interaction energy between two molecules.
The paper is organized as follows: In the Methods section, the three kinds of organic related systems with six testing cases are generated as the testing systems.The interaction energy between two molecules via QM is introduced.The architecture of DTNN is depicted in detail.And all the calculating details, which are for readers to reproduce our work, are also provided. In the Results section, the training, validating, and testing results of the interaction energy prediction via DTNN framework are staticly analyzed and summarized.In the Conclusion section, the conclusion of this proposed deep learning framework is given; and some aspects of possible improvements are proposed.
Ⅱ.
METHODS
The predictive power of a DL network can only be assessed by the precondition of building an enough reliable and convergent data set.In this session, the data generation process is introduced. The input data of the DTNN, including the two molecular conformations, are generated by different hierarchies distribution.The output data, the interaction energy between two molecules is calculated by QM.The illustration of the generation process is shown in FIG. 1.
Figure
1.
The schematic plot of two molecules system generation. (a) The interaction energy decays with the minimum distance of two molecules increasing. The scatter points are the interaction energies of two molecules in different ranges. (b) Different hierarchies generations of two molecules systems with increasing Radcutoff. (c) The distribution of Dismin of two molecules in the system of the whole data set with MaxRadcutoff.
In most organic related systems, the typical interaction is the multi-influence between two molecules.So the interaction energy of such two-organic molecule system is applied as the main study case.First, two-simple organic molecule system is introduced, and it is named as Org-org C4.Here, succinonitrile (C4N2H4) and propanenitrile (C3NH5) are chosen.Then, a more complex two organic molecular system is introduced named as Org-org C6. The 1, 1-diethoxyethane (C6O2H14) and C6NO4H13 are chosen.It contains molecule with the benzene ring.
In the organic MD simulation, the organic molecule is placed in a box full of solvating molecules.To study the interaction between the organic molecule and water molecule is an essential process to invest the organic-water system, named as Org-wat.To neutralize the pH value of the MD systems, the electronical ions are introduced into the solvating environment.We choose Na+ and K+ as the examples of positive ion andthe Cl− as the example of negative ion.The three cases are studied to present the interaction energy of organic molecule and each ion.All three of them are named as Org-ion in this work.The three testing systems related organic molecules, which are applied to valid predicting results of the proposed DL framework, are shown in FIG. 2(a).More details of the six testing cases are shown in supplementary materials.
Figure
2.
The schematic plot of the prediction interaction energy of organic related systems and workflow of the DTNN. (a) Three main organic related testing systems, Org-ion, Org-org, and Org-wat. They are all two molecules systems. (b) The inputs of DTNN are the atomic type and pair-wise distance matrix. (c) The output of DTNN is the interaction energy, which is calculated via QM. (d) The main network structure of DTNN.
The predicted results of supervised learning in deep learning mainly depend on two aspects.The first one is whether the network model fits the prediction problem.The second one is the good quality of the data set.The DTNN framework, which has been proven by previous works, is capable of predicting the physical and chemical properties in the molecular system.In this work, we focus on constructing a reliable data set for training the DTNN model.
In the data set, each pair of molA and molB is one input of the DTNN, and the interaction energy between these two molecules is its relative output of DTNN.The most reliable data set covers nearly all the possible molecular pairs and their relative positions.This kind of sampling can not be fulfilled in practice.By balancing both the calculating consumes and prediction precision, we applied the hierarchy generation to build the data set of three main testing systems related to the organic molecules.
The interaction energy between two molecules generally decays with the increasing of the minimum distance Dismin of the two molecules.
The minimum distance between two molecules is defined as the minimum of the set of distances between atomA and atomB, which is shown in Eq.(1).The atomA is an arbitrary atom of molA, which is one molecular in the testing system.In the same pattern, atomB is an arbitrary atom of molB, the other molecule is the same testing system.When the Dismin exceeds the threshold MaxRadcutoff, the interaction energy between two molecules decreases to 0.It is supposed that there exists no interaction between the two molecules in the system.After the position of molA is fixed, the molB is placed further.The MaxRadcutoff is decided via the changing tendency of interaction energy.In the range of MaxRadcutoff, different numbers of molecular system are generated with increasing Radcutoffi, which is shown in FIG. 1(b). The density of molB is higher within a small Radcutoffi.Within a large Radcutoffi, more positions of molB are easy to be generated.
In the hierarchy generation, even two molecules are in the same Dismin distance, they are still different relative positions samples and relative interaction energies in the data set.It is shown in FIG. 1(a).Therefore, the diversity of the data set increases.
After all two molecules systems are generated, the distribution of Dismin is shown in FIG. 1(c).It is easy to observe that the changing tendency of the generated data set fits the physical truth.
To the Org-org C4, four hierarchies of data are generated. And the MaxRadcutoff is set as 15 Å.There are totally 10×105 pairs of conformations in data set.To the Org-org C6, 5 hierarchies are applied.The minimum Radcutoff is extend to 6 Å. 2.5×104 pairs of conformations are generated.All the parameters used to generate the Org-org C4 and C6 systems are listed in Table Ⅰ.The parameters for data set generation of the other four testing systems are listed in supplementary materials.
Table
Ⅰ.
The parameters in the hierarchy generation of two organic molecules systems.
Since the electronic properties of molecules are determined by quantum mechanics, the chemical analyzing process has to be based on the fundamental quantum principle.In the testing system, each pair of molecules relates one interaction energy.Furthermore, this interaction energy characterizes the mutual influence between the two molecules in the system.An approach to calculating it is to calculate the difference between the energies of isolated molecules and their assembly [52].
In each testing case, the interaction energy ΔE(molA, molB)int between molA and molB can be calculated as Eq.(2).
ΔE(molA,molB)int=E(molA,molB)−(E(molA)+E(molB))
(2)
where E(molA) and E(molB) are the energies of the isolated molecules, and E(molA,molB) is the energy of their interacting assembly.
The deep tensor neural network is a deep learning framework.It was proposed for audio recognition first in 2017.In the hidden layers, the tensor production is simplified to improve the training efficiency in the same prediction precision.In 2017, the DTNN was introduced to predict the physical and chemical properties of molecules in GDB9.
Here, we introduce the details of this fundamental and how to apply it into the prediction of the interaction energy of two molecules.
C
DTNN
The main thought of the original DTNN is to decompose the energy of a whole molecular system into the sum of energy of individual "atomic" contributions in the environment bounded by a cutoff R.And this mechanism is introduced into the interaction energy calculating process of organic related systems.The DTNN applied in this model directly maps the organic relating two molecules system to the interaction energy, as shown in Eq.(3).
{RI,RII,ZI,ZII}→DTNNΔE
(3)
In which RI and RII denote the conformation relative position information (geometric information) of molA and molB in the system.The ZI and ZII denote the atomic type information of molA and molB. The distance matrix of the two molecules in each system and their atom types are the inputs of the DTNN framework.The interaction energy replaces the molecular energy and becomes the output of the DTNN applied in this work.
1
Network structure
In this subsection, the flowchart of calculating the interaction energy of organic molecule related cases is introduced, and it is shown in FIG. 2.
To three main test organic related systems, the organic molecule, the water, the ions, the coordinates of each atom in the system are provided to generate the distance matrix as a part of the input of the DTNN.
In this proposed calculating framework, the deep learning modular, DTNN, is the main part to predict the interaction energy.Here, we illustrate this part in detail and list the difference between the original DTNN and the framework proposed here.DTNN was first proposed by Schütt et al. [53].The purpose of this deep learning framework is to predict the energy of small molecular (20 is the highest atom number in GDB9 database [50]), and the physicochemical properties, such as the potential of aromatic rings.DTNN has also shown good predicting ability.We extend the DTNN to organic related systems and apply different DTNNs to its relative testing systems.
2
The inputs of DTNN
One of the essential aspects of creating a general deep learning model is the choice of appropriate data descriptors, which reflects the various properties of the two molecules system.As in most cases, the deep learning method efficiently handles the high-throughout data set with related properties, such as the huge amount of image files.In the DTNN module, the two molecules in testing systems are also converted into image format with serval layers of gray degree data, but not the typical color image with RGB (red, green, blue) three layers.The basis number (layers number) is decided by the number of the atoms in the testing systems.In the converting process, the two molecules are converted into the pair-wise inter-atom distance matrix D, which is shown in Eq.(4) and in FIG. 2(b).
D=[D11D12⋯D1nD21D22⋯D2n⋯⋯⋯⋯Dn1Dn2⋯Dnn]n×n
(4)
in which, n is the number of atom in two molecules system.
Dij is the pair-wise distance of atomi and atomj rather than the atomic positions.The map from the organic molecular system to pair-wise distance matrix is shown as Eq.(5).
Dij=distance(atomi,atomj)
(5)
Under rigid transformation, the matrix construction satisfies the rotational and translational invariance, which is convenient for both the MD simulation and the network training.As shown in FIG. 2(b), the distance between two atoms, which are in the same molecule, is usually smaller than the distance of two atoms belonging to two different molecules.But there are still cells in the mesh plot with a smaller distance, and the two atoms are in the different two molecules.These kinds of two atoms, which are in a short distance, commonly make a big contribution to the total interaction energy between the two molecules.
Only the geometric information of two molecules are not sufficient to distinguish the two kinds of molecules system.We encode the two molecules in the system by their atomic coordinates and atomic types.We also feed the atom types, which is a vector of nuclear charges, as the input of DTNN.It is shown as the atom type vector in FIG. 2(b, d).
Z=[Z1,Z2,⋯,Zn]
(6)
3
Gaussian feature expansion
Each element Dij∈D is spread across multi-dimensions of a uniform grid of Gaussians basis. This process is shown in Eq.(7).And the illustration of Gaussian feautre expansion is shown in supplementary materials.
ˆdij={exp[−Dij−(μmin+kΔμ)22σ2]}
(7)
in which, Δμ is the gap between two Gaussian basises with width σ.μmin is the center of minimum Gaussian basis.μmax is the center of maximum Gaussian basis.These two values bound the Gaussian basis range.
In this process, a conformation is converted to an image with multi-layers.And the pixel value in each layer is decided by the pair-wise distance of two atoms in conformation.
After the decomposition, the atom is characterized by atomic descriptors with B coefficients.
c(0)i=CZi∈RB
(8)
And B is the number of basis functions.The initial value of Cz is set by a random coefficient vector according to Eq.(9), and it is optimized during the neural network training process.
Cz≈N(0,1/√B)
(9)
4
T interaction passes
Each coefficient vector c(t)i, corresponding to atom i after t passes, is corrected by the interactions with other atoms of the two molecules in the testing system.
c(t+1)t=c(t)i+∑j≠ivij
(10)
Non-linear coupling between the atomic vector features and the interatomic distances is achieved by a tensor layer, such that the coefficient k of the refinement is given by
υijk=tanh[c(t)jVkˆdij+(Wcc(t)j)k+(Wdˆdij)k+bk]
(11)
where bk is the bias of feature k, Wc and Wd are the weights of atom representation and distance, respectively.Eq.(11) is simplified into Eq.(12) according to low-rank tensor factorization.
vij is defined as
vij=tanh{Wfc[(Wcfcj+bf1)∘(Wdfˆdij+bf2)]}
(12)
where '∘' is the element-wise multiplication.And Wcf, bf1, Wdf, bf2 and Wfc are the weight matrices and corresponding biases of atom representations, distances and resulting factors, respectively [54].The simplified weight tensors in the DTNN module are illustrated as eight cubics in FIG. 2(d).
5
Sum of interaction energy
The optimized vij is activated by an activation function in DTNN.This activation function is a hybrid tangent function.
The interaction energy of two molecules is the sum of energy contribution of each atom ΔEi.
oi=tanh(Wout1c(T)i+bout1)
(13)
The final stage is nl fully-connected layers, which is a commonly used layers in DL network, nl=2.The hidden layer oi possesses nhid neurons.
ΔˆEi=Wout2oi+bout2
(14)
The ΔˆEi is shifted to the mean ΔEμ and scaled by the s.d. ΔEσ of the interaction energy per atom estimated on the training set.
ΔEi=ΔEσΔˆEi+ΔEμ
(15)
The interaction energy of two molecules is obtained by Eq.(16).
E=∑iΔEi+∑jΔEj
(16)
In which, i∈(1,nmolA), and j∈(1,nmolB).nmolA is the number of atoms in molA.nmolB is the number of atoms in molB.
6
The output of DTNN
The interaction energies are the output of DTNN. They are calculated via high efficient QM soft-package Gaussian 09 [55].It should be mentioned that the deep learning methods can, in theory, be applied to any level of approximation of the interaction energy.The precondition is that the training data set with high-level QM calculation is convergent enough and the network structure is complex enough to descript the map from molecular conformations to interaction energy.But in practice, limited by the calculating source and calculating time-consuming, it is reasonable to balance the accuracy and QM calculating speed.Most of the organic data sets, such as GDB9, GDB13 [56], and GDB17 [57], are calculated with the B3LYP as the method and (6-311+G∗) as basis set.Here, we choose the same level as it is easy to merge the data set of hierarchy which can generate into the commonly used data set in future work.
7
Loss function
The loss function is defined as root mean square deviation (RMSD) of the interaction energy between the DTNN predicted and QM calculated, which is shown in Eq.(17).
lossfunction=√∑i(EDTNN−EQM)2ntrain
(17)
in which, ntrain is the number of training data.
8
K-fold
To eliminate the influence from the random choosing of data, and verify the robustness of the DTNN module, the K-fold mechanical is introduced.The data set is randomly split into K-folds.With this kind of segmentation, the training data set, validation set are trained and validated to prevent the early stopping and over-fitting.The remaining data are used as the testing data set.In all testing cases of three type testing systems, k=4.Two 1/4 data are merged as the training set.One 1/4 data are used as the validating data set.The left 1/4 data are for testing the model.
After the prediction, the interaction energy obtained from DTNN is compared with its related value calculated from QM.And the static analysis is implemented to evaluate the training of the neural network.
D
Simulation environment
In order that the two molecules system is not disrupted by the large van der Waals repulsive interactions, atoms from each molecule must keep at least a certain pairwise distance.Then, the two molecules in all six testing systems are generated via the open-source soft PackMol [58, 59], which is a package for building initial configurations from MD.
To store and transfer the molecular information, including the molecular geometric, physical, and chemical properties, the python-based atom simulation environment (ASE) [60, 61] is applied.All the data are stored in a binary database file, which is easy to be trained and validated.
The DTNN is trained and validated in the deep learning framework TensorFlow 1.3 [62].And all the training processes are implemented on Tesla K80 GPU with 8GB ram.
All of the DTNN models are trained with stochastic gradient descent using the ADAM optimizer [63] with mini-batches of 100 examples.The learning rate exponentially decays with a ratio 0.98 every 105 steps.The bias parameters bout1, bout2, and tensor parameters wout1 and wout2 are randomly initialized.All other weight matrices are initialized with uniform distributions.All the parameters to implement the training process are shown in supplementary materials.
Ⅲ.
RESULTS
The methology details of data set generation and DTNN architecture are given in the Methods section.In this section, the interaction energies, which are predicted via DTNN, are compared with the ground truth energies from QM calculation.Furthermore, the static analysis is implemented to evaluate the training of the neural network.Three organic related systems, including six testing cases, are applied.They are Org-org C4, Org-org C6, Org-wat, Org-Na, Org-K, and Org-Cl.The details are illustrated in Methods section and supplementary materials.The molecules in each testing cases are hierarchy generated, which is introduced in Methods section.We measured the performance of the DTNN prediction using K-fold cross-validation with 50% training, 25% validating, and 25% testing separation of the data set to observe the best hidden layer weights and prediction result.
We evaluate the predictions both from the aspects of regression and static analysis.Four static parameters of regression and error distribution (paraa, parab, errmean, errstd) are utilized to judge the testing results of DTNN.
In each panel of FIG. 3, the regression data are randomly chosen from one fold of the K-fold of data sets.Only the testing regression results are shown in FIG. 3.All four fold regressions results of six testing cases, including training, validating, and testing, are shown in supplementary materials.
Figure
3.
The regression results of QM calculation and DTNN prediction of six testing cases. (a) The comparison result of two simple organic molecules C4. (b) The comparison result of two complex organic molecules C6. (c) The comparison result of organic molecule and water. (d) The comparison result of organic molecule and Na+ ion. (e) The comparison result of organic molecule and K+ ion. (f) The comparison result of organic molecule and Cl− ion.
The regression results of the testing system of Org-org are shown in FIG. 3(a, b).The regression result of Org-wat is shown in FIG. 3(c).The regression results of Org-ion related systems, are shown in FIG. 3(d, e, f).The horizontal coordinate of each panel is the interaction energies calculated from QM (engQM).The vertical coordinate is the interaction energies predicted by the DTNN (engDTNN).
First, the predicted result is judged by a regression process, which is shown in Eq.(18).
Fit=paraa×QM+parab
(18)
The paraa is the slope of the fitted regression line.The parab is the intercept of the fitted regression line.Both the paraa and parab are fitted by the least-square method simultaneously.The fitted line is shown in red color in each panel of FIG. 3.Both conditions, paraa is close to 1 and parab is close to 0, that mean the DTNN predicted result has a good correlation to QM values.
The scatter regression points are shown as green circles in each panel of FIG. 3.As the predicted results exhibit good correlations between QM and DTNN predicted energies.The scatter points are near the line Fit=1×QM+0, which is shown as the black line in each panel.
Then, the error of DTNN prediction is analyzed statically.The absolute prediction error is defined as the interaction energy difference between QM and DTNN prediction, which is shown in Eq.(19).The errmean is the mean value of the absolute prediction error.The errstd is the standard deviation (std) value of the absolute prediction error.
Error=engQM−engDTNN
(19)
Because the Org-org C4 and Org-org C6 contain more atoms than the other four testing cases.The mean value of absolute errors is larger relatively.But the errmean of these two testing cases are still in the range of 0.4 kcal/mol.To Org-wat testing cases, the absolute errors made a gaussian distribution.The errmean is less than 0.1 kcal/mol.With a proper Rcutoff, the DTNN describes the relationship between atoms in organic and the ion accurately.With a wider changing range of interaction energy, the DTNN predict it within a high precision.In Org-ion cases, the mean value of absolute error reaches 0.05 kcal/mol.
The predicted results of testing data of all six testing cases in K-fold cross validations are statically analyzed and shown as boxplots in each panel of FIG. 4. The paraa and parab are shown in FIG. 4(a, b).Both of the two error static parameters errmean and errstd are shown in FIG. 4(c, d).
Figure
4.
Four static parameters to evaluate the DTNN prediction results. (a) The slope of the fitting line with QM and DTNN. (b) The intercept of the fitting line with QM and DTNN. (c) The mean value of the absolute value of absolute errors. (d) The std value of the absolute value of absolute errors.
As shown in FIG. 4 (a), nearly all the paraa in each testing case are in the range of [0.995 1.005].Only in one fold of Org-org C6, the paraa is near 0.985, which is still a high correlation coefficient.The parameters parab in FIG. 4(b) show the same conclusion.Only in Org-org C6, this parameter has a big fluctuation.
Ⅳ.
DISCUSSION
We have provided a general protocol, based on deep learning model (DTNN) for predicting the interaction energy of two molecules in organic related systems.We provided brief instructions on each step of the framework, including the data generation and preparation, the model architecture and building, the neural networks training and testing.And, the accuracy and efficiency of the proposed framework are illustrated by three organic related systems.The results show that the proposed framework is accurate, computationally efficient and scalable.
As a newly proposed framework, the whole framework of applying DTNN to predict the interaction energy is only an exploration to study the interaction prediction of two molecules.
The advantages of this framework are the following, (i) The high calculating speed: the main calculating time is spent on the DTNN module, which can nearly be ignored compared with the QM calculation.(ii) The reduction in computational cost, including the calculating time and computer power: after the DTNN is trained, there is no extra cost of predicting.(iii) The high precision of interaction energy prediction: because the DTNN module is trained based on the QM precision, the prediction precision is the same as the QM result.
Although the DTNN model does not outperform some other DL models, it still represents a new paradigm in predicting the interaction energy, with exciting opportunities for future improvements.And some considerable scopes still exist for their improvements.(1) More generalization: only three testing systems, which represent three main organic related environments in MD, are applied to validate and got acceptable results.(2) One network: each testing case requires its own network, which has its specific weights in hidden layers, to predict the interaction energy.In the future, all these prediction works can be merged into one network with only one group of weights of hidden layers.(3) High precision: the QM calculation is based on B3LYP/6-31+G∗. It is a popular DFT in numbers of MD processes. This precision of DFT is chosen because of the consideration of the calculating source. And most deep learning models related to the molecular properties are based on this precision.If the calculation resources are sufficient enough, a high precision, such as CCSD can be applied into this framework.
Supplementary materials: The parameters to re-implement the DTNN training process and the prediction results of six different testing cases, and serval illustrations of definition in the DTNN are given.
Ⅴ.
ACKNOWLEDGMENTS
The authors would like to thank Xiang-Da Peng at The University of Chicago for his fruitful discussion and encouragement for this work. This work was supported by the National Natural Science Foundation of China (No.21933010 to Guo-hui Li).
Author contributions statement
Yuan Qi, Hong Ren, and Yan Li for conceptualization, Hong Li and Ding-lin Zhang for methodology, Yuan Qi and Yan Li for software, Yuan Qi and Gui-yan Wang for validation, Yan Li and Gui-yan Wang for formal analysis, Yuan Qi for investigation, Yuan Qi, Jun-ben Weng, and Yan Li for data curation, Yan Li and Guo-hui Li for writing-original draft preparation, Yan Li, Yuan Qi, and Guo-huiLi for writing-review and editing, Yuan Qi, Yan Li, and Guo-hui Li for visualization, Guo-hui Li for project administration. All authors reviewed the manuscript.
SI of Interaction energy prediction of organic molecules using Deep Tensor Neural Network
ABSTRACT
The support informations of DTNN architecture and prediction results of six testing cases in K-fold are listed here.
Figures
0.1 Three testing systems
Three typical testing systems represent the three main MD environments of organic molecules.
Figure
1.
Three testing systems and six testing cases of organic molecules.
In the DTNN modular, the value of loss function depicts the training process with optimizing parameters in deep learning modular. The loss function is decreasing with the iteration number increase and converges into a fluctuating stage in the later training process.
In the beginning stage of optimization, the optimization convergent quickly. The loss value decreases dramatically from 3 × 104 to 1 × 104. The loss value of the first 1 × 104 iterations is not shown in figure A2 for its large change magnitudes. In the middle stage of optimizations, the optimization searches slowly of the weighs space of hidden layers. And the rate of the changing becomes slow, decreasing loss value slow. In the late stage optimization, the weights of the hidden layers optimization process are trapped into some local minimum in the weight spaces. The loss value fluctuates into a convergence stage, and gets to an acceptable optimization. This optimization result may not be (in most cases) a global optimization. But it is minimization enough for practical usage.
Figure
3.
The loss function changes with iteration. The zoom-in panel shows the iteration 4×104 to 8×104. In this stage, the training process gets into a convergent feasible region. And no better weights of hidden layers are searched.
In each fold of training process, there are two columns of sub-figures. The first column contains three vertical panels including the regression plots of training, validating, and testing results. The horizontal coordinate is the interaction energies calculated by QM. The vertical coordinate is the interaction energies predicted by the DTNN framework. The second columns contain three vertical panels including the error distribution bar plots of training, validatin, and testing results. The vertical coordinate is the distribution of error values.
The scatter plot of regression and bar plot of error distribution of Org-org C4 testing system are shown in each panel of Figure 6.
Figure
6.
The regression and error distribution of Org-org C4 of K-fold.
The static parameters of the regression and error analysis of K-folder in Org-org C4 system. The 4 panels are paraa, parab, errmean and errstd. In each panel, the three boxplots show the parameter of the training, validating and testing processes.
The box plot of static parameters of Org-org C4 testing system are shown in each panel of Figure 7.
Figure
7.
The boxplot of static parameters of Org-org C4 of K-fold.
The parameters in generation of the organic related data sets are shown in Table AI. Since the water and ion are easily generated and they are uniformly distributed in the environment of MD, only one R is applied.
Table
1.
The parameters in the generation of organic related systems.
H. Li, J. Chowdhary, L. Huang, X. B. He, A. D. MacKerell, and B. Roux, J. Chem. Theory Comput. 13, 4535 (2017). doi: 10.1021/acs.jctc.7b00262
[7]
J. W. Ponder, C. J. Wu, P. Y. Ren, V. S. Pande, J. D. Chodera, M. J. Schnieders, I. Haque, D. L. Mobley, D. S. Lambrecht, R. A. DiStasio Jr., M. Head-Gordon, G. N. I. Clark, M. E. Johnson, and T. Head-Gordon, J. Phys. Chem. B 114, 2549 (2010). doi: 10.1021/jp910674d
[8]
J. C. Wu, G. Chattree, and P. Y. Ren, Theor. Chem. Acc. 131, 1138 (2012). doi: 10.1007/s00214-012-1138-6
[9]
C. S. Zhang, C. Lu, Z. F. Jing, C. J. Wu, J. P. Piquemal, J. W. Ponder, and P. Y. Ren, J. Chem. Theory Comput. 14, 2084 (2018). doi: 10.1021/acs.jctc.7b01169
[10]
J. L. Gao, D. G. Truhlar, Y. J. Wang, M. J. M. Mazack, P. Löffler, M. R. Provorse, and P. Rehak, Acc. Chem. Res. 47, 2837 (2014). doi: 10.1021/ar5002186
M. J. M. Mazack, A. Cembran, and J. L. Gao, J. Chem. Theory Comput. 6, 3601 (2010). doi: 10.1021/ct100426m
[15]
D. E. Shaw, M. M. Deneroff, R. O. Dror, J. S. Kuskin, R. H. Larson, J. K. Salmon, C. Young, B. Batson, K. J. Bowers, J. C. Chao, M. P. Eastwood, J. Gagliardo, J. P. Grossman, C. R. Ho, D. J. Ierardi, I. Kolossváry, J. L. Klepeis, T. Layman, C. McLeavey, M. A. Moraes, R. Mueller, E. C. Priest, Y. B. Shan, J. Spengler, M. Theobald, B. Towles, and S. C. Wang, Commun. ACM 51, 91 (2008).
[16]
D. E. Shaw, J. P. Grossman, J. A. Bank, B. Batson, J. A. Butts, J. C. Chao, M. M. Deneroff, R. O. Dror, A. Even, C. H. Fenton, A. Forte, J. Gagliardo, G. Gill, B. Greskamp, C. R. Ho, D. J. Ierardi, L. Iserovich, J. S. Kuskin, R. H. Larson, T. Layman, L. S. Lee, A. K. Lerer, C. Li, D. Killebrew, K. M. Mackenzie, S. Y. H. Mok, M. A. Moraes, R. Mueller, L. J. Nociolo, J. L. Peticolas, T. Quan, D. Ramot, J. K. Salmon, D. P. Scarpazza, U. Ben Schafer, N. Siddique, C. W. Snyder, J. Spengler, P. T. P. Tang, M. Theobald, H. Toma, B. Towles, B. Vitale, S. C. Wang, and C. Young, Proceedings of 2014 International Conference for High Performance Computing, Networking, Storage and Analysis, New Orleans, LA, USA, (2014).
[17]
J. C. Phillips, J. E. Stone, and K. Schulten, Proceedings of 2008 ACM/IEEE Conference on Supercomputing, Austin, TX, USA, (2008).
[18]
N. Goga, S. Marrink, R. Cioromela, and F. Moldoveanu, Proceedings of the IEEE 12th International Conference on Bioinformatics & Bioengineering, Larnaca, Cyprus, (2012).
[19]
C. Kutzner, S. Páll, M. Fechner, A. Esztermann, B. L. de Groot, and H. Grubmüller, J. Comput. Chem. 36, 1990 (2015). doi: 10.1002/jcc.24030
[20]
J. E. Stone, A. P. Hynninen, J. C. Phillips, and K. Schulten, ISC High Performance 2016 International Workshops on High Performance Computing, Frankfurt, Germany, 188 (2016).
[21]
Y. González, P. Ezzatti, and M. Paulino, Proceedings of the 4th International Conference on Bioinformatics and Biomedical Engineering, Granada, Spain, 639 (2016).
S. Purawat, P. U. Ieong, R. D. Malmstrom, G. J. Chan, A. K. Yeung, R. C. Walker, I. Altintas, and R. E. Amaro, Biophys. J. 112, 2469 (2017). doi: 10.1016/j.bpj.2017.04.055
[24]
D. J. Mermelstein, C. Lin, G. Nelson, R. Kretsch, J. A. McCammon, and R. C. Walker, J. Comput. Chem. 39, 1354 (2018). doi: 10.1002/jcc.25187
[25]
A. Crespo, A. Rodriguez-Granillo, and V. T. Lim, Curr. Top. Med. Chem. 17, 2663 (2017).
[26]
S. M. Bachrach, Computational Organic Chemistry (2nd Edn. ), S. M. Bachrach Ed., Hoboken: John Wiley & Sons, Inc., 1 (2014).
[27]
Y. LeCun, Y. Bengio, and G. Hinton, Nature 521, 436 (2015). doi: 10.1038/nature14539
[28]
M. Wainberg, D. Merico, A. Delong, and B. J. Frey, Nat. Biotechnol. 36, 829 (2018). doi: 10.1038/nbt.4233
[29]
A. Gupta, P. J. Harrison, H. Wieslander, N. Pielawski, K. Kartasalo, G. Partel, L. Solorzano, A. Suveer, A. H. Klemm, O. Spjuth, I. M. Sintorn, and C. Wahlby, Cytom. Part A 95, 366 (2019). doi: 10.1002/cyto.a.23701
[30]
M. Z. Alom, T. M. Taha, C. Yakopcic, S. Westberg, P. Sidike, M. S. Nasrin, M. Hasan, B. C. Van Essen, A. A. S. Awwal, and V. K. Asari, Electronics 8, 292 (2019). doi: 10.3390/electronics8030292
G. Montavon, M. Rupp, V. Gobre, A. Vazquez-Mayagoitia, K. Hansen, A. Tkatchenko, K. R. Müller, and O A. von Lilienfeld, New J. Phys. 15, 095003 (2013). doi: 10.1088/1367-2630/15/9/095003
[34]
J. Behler, S. Lorenz, and K. Reuter, J. Chem. Phys. 127, 014705 (2007).
M. Rupp, A. Tkatchenko, K. R. Müller, and O. A. von Lilienfeld, Phys. Rev. Lett. 108, 058301 (2012). doi: 10.1103/PhysRevLett.108.058301
[40]
K. T. Schütt, P. J. Kindermans, H. E. Sauceda, S. Chmiela, A. Tkatchenko, and K. R. Müller, Advances in Neural Information Processing Systems 30, Long Beach, CA, USA, (2017).
[41]
S. Chmiela, H. E. Sauceda, I. Poltavsky, K. R. Müller, and A. Tkatchenko, Comput. Phys. Commun. 240, 38 (2019). doi: 10.1016/j.cpc.2019.02.007
[42]
K. Gubaev, E. V. Podryabinkin, and A. V. Shapeev, J. Chem. Phys. 148, 241727 (2018).
[43]
N. Lubbers, J. S. Smith, and K. Barros, J. Chem. Phys. 148, 241715 (2018).
[44]
K. T. Schütt, M. Gastegger, A. Tkatchenko, and K. R. Müller, Explainable AI: Interpreting, Explaining and Visualizing Deep Learning, W. Samek, G. Montavon, A. Vedaldi, L. K. Hansen, K. R. Mülle Eds., Cham, Switzerland: Springer International Publishing, 311 (2019).
[45]
L. F. Zhang, J. Q. Han, H. Wang, R. Car, and E. Weinan, Phys. Rev. Lett. 120, 143001 (2018). doi: 10.1103/PhysRevLett.120.143001
[46]
H. Wang, L. F. Zhang, J. Q. Han, and E. Weinan, Comput. Phys. Commun. 228, 178 (2018). doi: 10.1016/j.cpc.2018.03.016
[47]
S. Kearnes, K. McCloskey, M. Berndl, V. Pande, and P. Riley, J. Comput. Aided Mol. Des. 30, 595 (2016). doi: 10.1007/s10822-016-9938-8
[48]
D. Yu, L. Deng, and F. Seide, IEEE Trans. Audio, Speech, Language Process. 21, 388 (2013). doi: 10.1109/TASL.2012.2227738
[49]
M. Gastegger, L. Schwiedrzik, M. Bittermann, F. Berzsenyi, and P. Marquetand, J. Chem. Phys. 148, 241709 (2018).
[50]
R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, Sci. Data 1, 140022 (2014). doi: 10.1038/sdata.2014.22
[51]
Y. Li, G. H. Li, H. Y. Min, Z. B. Dong, T. Yuan, X. Q. Li, and P. J. Xu, Commun. Inf. Sys. 18, 229 (2018).
K. T. Schütt, F. Arbabzadah, S. Chmiela, K. R. Müller, and A. Tkatchenko, Nat. Commun. 8, 13890 (2017). doi: 10.1038/ncomms13890
[54]
G. W. Taylor and G. E. Hinton, Proceedings of the 26th Annual International Conference on Machine Learning, Montreal, Quebec, Canada, 1025 (2009).
[55]
M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, R. W. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo-Merchan, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, Gaussian 09, Revision A02, Wallingford CT: Gaussian, Inc., (2009).
[56]
L. C. Blum and J. L. Reymond, J. Am. Chem. Soc. 131, 8732 (2009). doi: 10.1021/ja902302h
[57]
L. Ruddigkeit, R. van Deursen, L. C. Blum, and J. L. Reymond, J. Chem. Inf. Model. 52, 2864 (2012). doi: 10.1021/ci300415d
[58]
J. M. Martínez and L. Martínez, J. Comput. Chem. 24, 819 (2003). doi: 10.1002/jcc.10216
[59]
L. Martínez, R. Andrade, E. G. Birgin, and J. M. Martínez, J. Comput. Chem. 30, 2157 (2009).
[60]
S. R. Bahn and K. W. Jacobsen, Comput. Sci. Eng. 4, 56 (2002).
[61]
A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. H. Zeng, and K. W Jacobsen, J. Phys. : Condens. Matter 29, 273002 (2017). doi: 10.1088/1361-648X/aa680e
[62]
M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. F. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Q. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Q. Zheng, arXiv: 1603.04467 (2016).
Xu, P., Mou, X., Guo, Q. et al. Coarse-grained molecular dynamics study based on TorchMD†. Chinese Journal of Chemical Physics, 2021, 34(6): 957-969.
DOI:10.1063/1674-0068/cjcp2110218
DownLoad:
