DFTB and NN methods

A detailed description of the DFTB formalism is given in Refs. Reference Elstner, Porezag, Jungnickel, Elsner, Haugk, Frauenheim, Suhai and Seifert4 and Reference Gaus, Cui and Elstner21. For the sake of brevity, we only introduce here essential traits of the DFTB family of approximations. The members of the hierarchy of DFTB methods are related to the order up to which the exact Kohn–Sham DFT energy of a molecular system E[ρ] is expanded in a Taylor series around a reference electron density ρ ⁰ with respect to a difference electron density δρ as

(1)$$\eqalign{E[ \rho ] = & E^0[ {\rho^0} ] + E^1[ {\rho^0,\delta \rho} ] + E^2[ {\rho^0,{\lpar {\delta \rho} \rpar }^2} ] \cr & + E^3[\rho ^0,\lpar {\delta \rho )^3} ] + \ldots,}$$

where ρ = ρ ⁰ + δρ.^{[Reference Foulkes and Haydock1]} The first-order DFTB (DFTB or DFTB1) method^{[Reference Porezag, Frauenheim, Köhler, Seifert and Kaschner2]} includes only the first two terms in Eq. (6), the second-order DFTB [self-consistent charge DFTB (SCC-DFTB) or DFTB2] method^{[Reference Elstner, Porezag, Jungnickel, Elsner, Haugk, Frauenheim, Suhai and Seifert4]} includes the first three terms, and the third-order DFTB (DFTB3) method^{[Reference Gaus, Cui and Elstner21]} includes the first four terms.

In our formulation of the DFTB-NN methodology, we follow the work of Shen et al.^{[Reference Shen, Wu and Yang19,Reference Shen and Yang20]} and utilize the difference between a given Kohn–Sham DFT energy and the ordinary DFTB energy following the Δ-machine learning approach.^{[Reference Ramakrishnan, Dral, Rupp and von Lilienfeld17]} Note that we could, in principle, select any high-level method to define this difference. However, following the true spirit of DFTB, it is natural to utilize the NN to estimate the remainder of higher-order terms of the Taylor expansion in Eq. (6), as well as the effects of minimal basis set, two-center approximations, etc., as they can be thought of as introducing systematic errors that an NN should be easily able to spot as “patterns” in the data points it is presented. Consequently, the DFTB-NN energy prediction is simply the sum of the DFTB energy and the NN energy correction. For a given chemical system with a geometry “0” chosen as reference geometry, the DFT-relative energies ΔE _DFT(X) for arbitrary geometries X are calculated as follows:

(2)$${\rm \Delta }E_{{\rm DFT}}\left( X \right) = E_{{\rm DFT}}\left( X \right)-E_{{\rm DFT}}\left( 0 \right),$$

where E _DFT(X) and E _DFT(0) are DFT energies of geometry X and the reference geometry 0, respectively. Correspondingly, the DFTB relative energy ΔE _DFTB(X) for a geometry X is calculated as follows:

(3)$${\rm \Delta }E_{{\rm DFTB}}\left( X \right) = E_{{\rm DFTB}}\left( X \right)-E_{{\rm DFTB}}\left( 0 \right),$$

where E _DFTB(X) and E _DFTB(0) are DFTB energies of geometry X and the reference geometry “0”. Both the geometry X and the reference geometry are the same for both DFT and DFTB calculations. The NN correction target energy for a geometry X is defined as the difference between DFT- and DFTB-relative energies:

(4)$$E_{{\rm NN}}\left( X \right) = {\rm \Delta }E_{{\rm DFT}}\left( X \right)-{\rm \Delta }E_{{\rm DFTB}}\left( X \right).$$

The DFTB-NN energy for geometry X is then simply calculated as follows:

(5)$$\eqalign{E_{{\rm DFTB - NN}}\left( X \right) = &E_{{\rm DFTB}}\left( 0 \right) + {\rm \Delta }E_{{\rm DFTB}}\left( X \right) \cr & + E_{{\rm NN}}\left( X \right).} $$

Following Behler–Parrinello,^{[Reference Behler and Parrinello11]} this “NN energy” E _NN for geometry X is a sum of atomic energy contributions

(6)$$E_{{\rm NN}}\lpar X \rpar = \mathop \sum \limits_{e = 1}^M \mathop \sum \limits_{\,j = 1}^{N_{{\rm atom}}^e} E_j^e \lpar X \rpar,$$

where $E_j^e$ is the atomic energy of atom j of chemical element e, M is the number of chemical elements, and $N_{{\rm atom}}^e$ is the number of atoms belonging to element e the given system.

For each chemical element present in the molecular system, an atomistic feed-forward BPNN is built to predict the atomic energy of the atomic energy contribution $E_j^e$. For a high-dimensional NN potential, its output must be kept invariant with respect to permutation of the atomic index and to the rotation and translation of the input geometry to accurately represent the potential energy surface. Following Nguyen et al.’s approach,^{[Reference Nguyen, Székely, Imbalzano, Behler, Csányi, Ceriotti, Götz and Paesani18]} we select two types of permutational, rotational, and translational invariant atomic number differentiated symmetry function vectors as the input of the NN. For each atom, a set of radial and angular symmetry function vectors are calculated from the neighboring atomic coordinates to describe its immediate chemical environment. The radial symmetry function vector $G_{i,m}^{{\rm Rad}}$ for atom i is a set of element-based contributions $G_{i,m,Z}^{{\rm Rad}}$ defined as follows:

(7)$$G_{i,m,Z}^{{\rm Rad}} = \mathop \sum \limits_{\,j\ne i}^{N_{{\rm atom}}^Z} e^{-\eta _m{(R_{ij}-R_{s,m})}^2},$$

where j is the index of atoms with the same atomic number Z, R _ij is the distance between atom i and j, and the η _m and R _s,m are a set of arbitrary hyperparameters used in the mth symmetry function vector.^{[Reference Behler13]} The angular symmetry function vector $G_{i,m}^{{\rm Ang}}$ for atom i is a set of element-based contribution $G_{i,m,Z_1,Z_2}^{{\rm Ang}}$ calculated as follows:

(8)$$G_{i,m,Z_1,Z_2}^{{\rm Ang}} = \mathop \sum \limits_{\,j\ne i}^{N_{{\rm atom}}^{Z_1}} \mathop \sum \limits_{k\ne i,j}^{N_{{\rm atom}}^{Z_2}} \displaystyle{{2^{1-\zeta _m}{\lpar {1 + \lambda_m{\rm cos}\theta_{ijk}} \rpar }^{\zeta _m}} \over {e^{\eta _m{(R_{ij}^2 + R_{ik}^2 + R_{\,jk}^2 )}^2}}},$$

where j and k correspond to the index of atoms with the same atomic number Z ₁and Z ₂, R _ik is the distance between atom i and k, R _jk is the distance between atom j and k, θ _ijk is the angle formed by atoms i, j, and k, and ζ _m and λ _m are hyperparameters used in the mth symmetry function vector.^{[Reference Behler13]} Similar to Smith et al.’s atomic environment vector,^{[Reference Smith, Isayev and Roitberg15]} only contributions from neighboring atoms of the same chemical element are summed into the same symmetry function vector in order to differentiate the atom type. In this proof-of-concept work, we do not introduce any cutoff for symmetry functions. The implementation of cutoff functions will be exploited in the future work.

A set of sub-NNs are built for all the elements present in the molecular system. The elemental sub-NNs use the symmetry function for the input and output of the atomic energy $E_j^e$ in Eq. (6). Each feed-forward sub-NN consists of three layers with tanh as an activation function and a linear layer to sum up all the node outputs into one atomic energy $E_j^e$. The output of node n in the mth layer is given by the following equation:

(9)$$y_n^m = f_n^m \left( {b_n^m + \mathop \sum \limits_{\,p = 1}^{N_{m-1}} a_{\,p,n}^{m-1,m} \,\ast\ y_n^{m-1}} \right),$$

where $y_n^m$ is the node output, $f_n^m$ is the activation function, $a_{q,n}^{m-1,m}$ is the weight-connecting nodes between neighboring layers, and $y_n^{m-1}$ is the output of the node in the previous layer. A detailed review of the BPNNs can be found in Refs. Reference Behler and Parrinello11 and Reference Behler13.

A four-layer NN is constructed for each element from the open-source machine learning package Keras^[22] using the TensorFlow backend.^{[Reference Abadi, Barham, Chen, Chen, Davis, Dean, Devin, Ghemawat, Irving, Isard, Kudlur, Levenberg, Monga, Moore, Murray, Steiner, Tucker, Vasudevan, Warden, Wicke, Yu and Zheng23]} The tanh function is selected as the activation function for the first three layers and end with a linear layer to sum up into one atomic energy. Each layer consists of 34 nodes and the Nguyen–Widrow weight initialization^{[Reference Nguyen and Widrow24]} procedure was applied to initialize the weights for each elemental sub-NN.

In the DFTB-NN scheme, the force on atom i with respect to changes in a component of its Cartesian coordinate α is the sum of the ordinary DFTB force and the analytic NN gradient correction, obtained via a Jacobian matrix ∂G _j,m/∂α as

(10)$$F_{{\rm DFTB - NN}}^{i,\alpha} = F_{{\rm DFTB}}^{i,\alpha} -\mathop \sum \limits_{\,j = 1}^{N_{{\rm atom}}} \mathop \sum \limits_{m = 1}^{N_{{\rm sym}}} \displaystyle{{\partial E_j^e} \over {\partial G_{\,j,m}}}\displaystyle{{\partial G_{\,j,m}} \over {\partial \alpha}}. $$

The symbolic differentiation in TensorFlow was used to calculate the gradient of the NN energy with respect to each symmetry function vector. We employed the Adam optimizer implemented in the Keras python package to train the DFTB-NN model, minimizing the loss function f _loss, which is defined as follows:

(11)$$f_{{\rm loss}} = \sqrt {\displaystyle{1 \over N}\mathop \sum \limits_{X = 1}^L {\lpar {{\rm \Delta} E_{{\rm DFTB}}\lpar X \rpar + E_{{\rm NN}}\lpar X \rpar -{\rm \Delta} E_{{\rm DFT}}\lpar X \rpar } \rpar }^2},$$

where L is the number of geometries in the training set, and ΔE _DFTB(X) and ΔE _NN(X) are the DFTB-relative energy and NN energy corrections for geometry X.

Computational details

As an illustrative example, the accuracy of the implemented DFTB-NN energy and the gradient was evaluated for the glycine molecule in both static and dynamic calculations. To simplify the discussion, we will only present data based on the DFTB2 method, using the MIO Slater–Koster parameter set.^{[Reference Elstner, Porezag, Jungnickel, Elsner, Haugk, Frauenheim, Suhai and Seifert4]} However, our implementation is also applicable to other members of the DFTB hierarchy of methods. As target high-level method, we selected the B3LYP functional^{[Reference Becke25]} with Pople's 6-31G** basis set. For DFTB- and DFTB-NN MD simulations, we used the Velocity Verlet algorithm implemented in the Atomic Simulation Environment (ASE) library.^{[Reference Hjorth Larsen, Jørgen Mortensen, Blomqvist, Castelli, Christensen, Dułak, Friis, Groves, Hammer, Hargus, Hermes, Jennings, Bjerre Jensen, Kermode, Kitchin, Leonhard Kolsbjerg, Kubal, Kaasbjerg, Lysgaard, Bergmann Maronsson, Maxson, Olsen, Pastewka, Peterson, Rostgaard, Schiøtz, Schütt, Strange, Thygesen, Vegge, Vilhelmsen, Walter, Zeng and Jacobsen26]} All DFT calculations were computed using the NWChem program.^{[Reference Valiev, Bylaska, Govind, Kowalski, Straatsma, Van Dam, Wang, Nieplocha, Apra, Windus and de Jong27]}

To construct a DFTB-NN model for the glycine molecule, a total of 276 symmetry function vectors are generated for each atom. Following Nguyen et al.’s approach,^{[Reference Nguyen, Székely, Imbalzano, Behler, Csányi, Ceriotti, Götz and Paesani18]} the radial component uses 24 vectors for each element calculated by Eq. (7) with the R _s,m evenly spread from 0.5 to 5 Å and $1/\sqrt {2\eta _m} = R_{s,m}$. The angular symmetry functions calculated from Eq. (8) adopt the parameters η _m = [0.001, 0.01, 0.05] Å⁻¹, λ _m = ±1, and ζ _m = [1, 4, 6]. To prevent overfitting, the early stopping method was employed and two dropout layers were applied to the second and third hidden layers of the elemental sub-NN. During the training process, the NN was continuously trained with a descending learning rate of [5 × 10⁻⁵, 1.5 × 10⁻⁵] with corresponding epochs of [4000, 2000] and a decay rate of 1 × 10⁻⁵.

Accuracy of DFTB-NN-relative energies

Our choice for using glycine as a test case of the DFTB-NN approach was motivated by the fact that DFTB2 (and likewise DFTB3) tends to underestimate the strength of the O–H · · · N hydrogen bond relative to that of the O–H · · · O and N–H · · · O hydrogen bond, as can be seen in Supplementary Fig. S1.

A training and two test sets were generated focusing on variations in the dihedral angles of the glycine molecule, followed by a short DFTB2 MD simulation. For the training set, 50 initial conformations were randomly generated using a modified version of the Kick³ program^{[Reference Addicoat, Fukuoka, Page and Irle28]} and followed by 100 steps of canonical ensemble (NVT) MD simulations at 300 K with a step size of 1 fs. In total, 5000 structures and DFT energies were obtained for the training set. For the first test set named Glycine300K, 100 initial conformations were generated and followed by 1000 steps NVT MD simulations at 300 K. In total, 100,000 structures were obtained for the Glycine300K test set. The second test set named Glycine400K was generated in the same way as the first test set, but the NVT MD simulations were carried out at 400 K in order to generate geometries with higher variations than the training set. In total, we generated 200,000 DFT reference points for the test set. The size of the training set of glycine conformations was kept deliberately much smaller than that of the test set in order to avoid overfitting of the NN.

Figure 1 shows the comparison of DFTB-, DFTB-NN-, and DFT-relative energies ΔE for training, Glycine300K, and Glycine400K test sets, respectively. Overall, the DFTB-NN improves the accuracy of the DFTB-relative energy prediction. The DFTB method achieves a root mean square error (RMSE) of 2.39 and 2.54 kcal/mol for the Glycine300K and Glycine400K test sets, respectively. The DFTB-NN approach decreases the RMSE values significantly to 0.59 and 0.71 kcal/mol for the test sets.

Figure 1.

Comparison of DFTB and DFTB-NN ΔE to DFT ΔE for (a) the training set, (b) the Glycine300K test set, and (c) the Glycine400K test set.

Even though the DFTB-NN was trained using MD snapshots generated at 300 K with B3LYP ΔE in a range of [4.67, 35.80] kcal/mol, its RMSE for the Glycine400K test set, MD snapshots generated at 400 K with ΔE in the range of [2.89, 39.92] kcal/mol, is only 0.12 kcal/mol higher than that for the Glycine300K test set. The small difference in RMSE for the two test sets indicates that the trained DFTB-NN model is applicable within this temperature span. The NN is thus shown to be capable of learning from the relatively small training set and accurately predicts the high-dimensional potential energy surface.

Accuracy of the DFTB-NN gradient

In order to validate the implemented DFTB-NN analytic gradient, the difference between the analytic and numerical gradients was evaluated for a particular geometry of the glycine molecule. The maximum difference between the analytic and the numerical gradients was found to be only 6.0 × 10⁻⁶ hartree/bohr. Furthermore, the accuracy of the DFTB-NN gradient was evaluated by examining the energy conservation of DFTB-NN and traditional DFTB MD simulations in the microcanonical ensemble (NVE). NVE MD simulations were performed for 10 ps with the time step of 0.5 fs using velocities initiated at 300 K. As shown in Fig. 2, the DFTB-NN was able to keep the variation in total energy (kinetic MD energy plus total electronic + nuclear energy as potential energy) during the MD simulation lower than ± 0.03 kcal/mol, which is at the same level as DFTB. Ideally, the total energy needs to be perfectly conserved in such NVE simulations, and the deviation of actual total energy from 0 is a good measure of the accuracy of the MD simulation. Figure 2 also shows that the mean of energy fluctuation is very close to zero and of the same magnitude for both DFTB and DFTB-NN, indicating the good accuracy of the DFTB-NN force calculation.

Figure 2.

DFTB and DFTB-NN total energy of the glycine molecule fluctuate as a function of NVE simulation time.

To demonstrate the computational performance of the implemented DFTB-NN, we compared the wall-clock time required to carry out a single-point energy + gradient calculation for the glycine molecule. On a single Xeon 3.00 GHz (E3-1505M) CPU core, the calculation took 0.256 s with DFTB2 and 0.302 s with DFTB2-NN. The wall-clock time was averaged over 1000 steps of a MD trajectory.

Potential energy profile of the N–C–O–H dihedral angle

DFTB neglects three-center, four-center, and higher-order many-body contributions as well as multipole interactions, which particularly affect the description of the chemical environment related to dihedral angles and (intramolecular) hydrogen bonds.^{[Reference Lee, Schnupf, Sumpter and Irle8]} To examine the accuracy of our trained DFTB-NN model, the DFTB-NN potential energy profile for a 360° rotation of the N–C–O–H dihedral angle of the glycine molecule was compared to the corresponding PES computed by B3LYP. At first, the glycine geometry was optimized using the B3LYP/6-31G** method, following a partially relaxed potential energy surface scan at the same level of theory for the N–C–O–H dihedral angle of the molecule. The corresponding energies for each geometry along the partially relaxed energy scan were calculated at the DFTB2 and DFTB2-NN levels of theory. The only other frozen parameter during the relaxed scan was the N–C–C–O dihedral angle, which we fixed at −7.0 in order to avoid sudden jumps between different conformations along the energy profile One sudden geometrical change remains for which the NH₂ group suddenly changes its orientation, see Fig. 3 (small insets at 282.3). The geometries of this energy scan are given in Supplementary Material.

Figure 3.

Potential energy profile for DFTB2 and DFTB2-NN in comparison to B3LYP/6-31G** as a function of the N–C–O–H dihedral angle in glycine.

Figure 3 illustrates that DFTB2 incorrectly predicts the global minimum at the position of the N–H · · · N intramolecular hydrogen bond, whereas DFT predicts, in agreement with other higher levels of theory, that the true global minimum should be located at the position of the O–H · · · N bond. As a consequence, DFTB2 predicts incorrectly that the cis-isomer of glycine with N–C–O–H dihedral angle of 0° is less stable than the trans-isomer of 180°. DFT, on the other hand, predicts that the cis-isomer is lower in energy. Clearly, DFTB2 fails to yield accurate prediction, because it uses Coulombic interactions between atomic point charges that cannot accurately model the hydrogen atom interactions.^{[Reference Christensen, Kubar, Cui and Elstner7]} As a remedy, the DFTB-NN is capable of reproducing the DFT PES with a maximum absolute error of 1.7 kcal/mol (9.4 kcal/mol in the case of DFTB) and predicts that the cis-isomer is more stable. It is apparent that DFTB-NN corrects faithfully flaws in the DFTB hydrogen bonding energies and approaches the accuracy of DFT without modifying the DFTB formalism itself, as for instance by the use of hydrogen corrections.^{[Reference Rezáč, Riley and Hobza29,Reference Rezáč30]}

Energy landscape of glycine

To further evaluate the accuracy of DFTB-NN for the description of intramolecular hydrogen bonds and the dihedral angle, we examined the potential energy landscape of the C–C–O–H and N–C–C–O dihedral angles using all geometries of the Glycine300K test set. Figure 4 shows the relative energy of DFTB, DFTB-NN, and B3LYP varying with the two dihedral angles. The scale of ΔE is mapped according to the color scale shown in the plot.

Figure 4.

DFTB-, DFT-, and DFTB-NN-relative energies as a function of C–C–O–H and N–C–C–O dihedral angles of glycine.

DFTB systematically underestimates the potential energy when compared with DFT. Overall, the DFTB-NN energy landscape matches the B3LYP energy landscape precisely. Remaining deviations exist for O–C–C–H dihedral angles in the range of [0°, 60°] and [90°, 120°], although the magnitude of these deviations is relatively small.