Nonadiabatic Derivative Couplings Calculated Using Information of Potential Energy Surfaces without Wavefunctions: Ab Initio and Machine Learning Implementations

In this work, we implemented an approximate algorithm for calculating nonadiabatic coupling matrix elements (NACMEs) of a polyatomic system with ab initio methods and machine learning (ML) models. Utilizing this algorithm, one can calculate NACMEs using only the information of potential energy surfaces (PESs), i.e., energies, and gradients as well as Hessian matrix elements. We used a realistic system, namely CH2NH, to compare NACMEs calculated by this approximate PES-based algorithm and the accurate wavefunction-based algorithm. Our results show that this approximate PES-based algorithm can give very accurate results comparable to the wavefunction-based algorithm except at energetically degenerate points, i.e., conical intersections. We also tested a machine learning (ML)-trained model with this approximate PES-based algorithm, which also supplied similarly accurate NACMEs but more efficiently. The advantage of this PES-based algorithm is its significant potential to combine with electronic structure methods that do not implement wavefunction-based algorithms, low-scaling energy-based fragment methods, etc., and in particular efficient ML models, to compute NACMEs. The present work could encourage further research on nonadiabatic processes of large systems simulated by ab initio nonadiabatic dynamics simulation methods in which NACMEs are always required.


Introduction
The Born-Oppenheimer (B-O) approximation is widely used in the research of computational and theoretical chemistry when potential energy surfaces (PESs) are far away from each other. However, the B-O approximation will break down when two or more PESs come close or even cross. In such a situation, nonadiabatic effects are non-negligible and nonadiabatic couplings between different PESs should be considered seriously. Nonadiabatic coupling matrix elements (NACMEs) are important physical quantities that are used for measuring coupling strengths between different adiabatic electronic states. Additionally, NACMEs play important roles in calculating non-radiative rates and nonadiabatic molecular dynamic (NAMD) simulations [1][2][3][4][5][6][7][8][9][10]. However, the computational cost of NACMEs is generally pretty expensive using ab initio methods. Moreover, not all electronic structure programs and methods have implemented wavefunction-based algorithms for calculating NACMEs yet. From this point of view, it is desirable to develop a simple and generalized method for obtaining NACMEs which can be easily adopted by most excited-state electronic structure methods. To achieve this target, Köppel et al. proposed a strategy for constructing diabatic states from PESs of adiabatic states, which can compute nonadiabatic couplings without wavefunction information [11]. Lasorne et al. re-investigated based algorithms for calculating NACMEs yet. From this point of view, it is desirable to develop a simple and generalized method for obtaining NACMEs which can be easily adopted by most excited-state electronic structure methods. To achieve this target, Köppel et al. proposed a strategy for constructing diabatic states from PESs of adiabatic states, which can compute nonadiabatic couplings without wavefunction information [11]. Lasorne et al. re-investigated and compared wavefunction-and PES-based NACMEs and suggested that PES-based NACMEs are accurate in regions around conical intersections [12]. Richardson reported a new machine learning (ML) approach for eliminating the issue with the double-valued nature of NACMEs by a set of auxiliary single-valued functions [13]. Recently, Baeck and An developed a practical approximation for the calculation of NACMEs and successfully applied it in one-dimensional model systems [14,15]. This method is easy to implement and has the potential to combine with electronic structure methods that can provide PES information. Very recently, Gastegger and coworkers used this approximate PES-based algorithm in their developed SchNarc machine learning approach to calculate NACMEs for NAMD simulations [16]. However, the accuracy of this approximate method to calculate NACMEs for polyatomic molecular systems still needs comprehensive benchmarks and in-depth analysis.
Motivated by these facts, we have in this work implemented an approximate PESbased method for NACMEs using a realistic system CH2NH ( Figure 1) and proved that this algorithm can give accurate NACMEs if PES information (i.e., energies, gradients and Hessian matrices) can be accurately provided by either ab initio methods or ML-trained models. Our work indicates that this approximate method performs very well even near PES regions with quasi-degenerate energies. However, this approximate algorithm will fail and become divergent if relevant PESs are rigorously degenerate in energy, i.e., truly conical intersections (CIs). Nevertheless, this PES-based algorithm remains useful and beneficial for combining electronic structure methods, energy-based fragment methods, and ML models, which can provide PES information but with either difficult-or impossible-to-supply wavefunctions, and has potential applications in nonadiabatic dynamics simulations, etc.  Figure 2 shows the wavefunction-and PES-based NACMEs with structures near the optimized CI structure of CH2NH. It can be found that the PES-based NACMEs at the optimized CI structure (energy gap less than 0.17 kcal/mol) not only have large values but also are far away from the wavefunction-based NACMEs (red points in Figure 2). This is originated from the fact that NACMEs at truly degenerate CI points become divergent because their values are inversely proportional to related energy gaps. Note that at the optimized CI structure of CH2NH, the energy gap is only 0.17 kcal/mol, which makes NACMEs a little divergent. On the other hand, the PES-based NACMEs at the structures  Figure 2 shows the wavefunction-and PES-based NACMEs with structures near the optimized CI structure of CH 2 NH. It can be found that the PES-based NACMEs at the optimized CI structure (energy gap less than 0.17 kcal/mol) not only have large values but also are far away from the wavefunction-based NACMEs (red points in Figure 2). This is originated from the fact that NACMEs at truly degenerate CI points become divergent because their values are inversely proportional to related energy gaps. Note that at the optimized CI structure of CH 2 NH, the energy gap is only 0.17 kcal/mol, which makes NACMEs a little divergent. On the other hand, the PES-based NACMEs at the structures far away from the CI structure (energy gap larger than 0.17 kcal/mol) are very close to the wavefunction-based NACMEs (see blue points and insert panel of Figure 2). The small deviations may be caused by numerically calculated Hessian matrix elements. In total, the PES-based NACMEs become much better compared with the wavefunction-based NACMEs when the energy gap is increased. far away from the CI structure (energy gap larger than 0.17 kcal/mol) are very close to the wavefunction-based NACMEs (see blue points and insert panel of Figure 2). The small deviations may be caused by numerically calculated Hessian matrix elements. In total, the PES-based NACMEs become much better compared with the wavefunction-based NACMEs when the energy gap is increased. To give insight into the accuracy of the PES-based NACMEs, we separated the structures near the optimized CI structure into several groups according to their energy gaps Δ (i.e., 0.17-1 kcal/mol, 1-3 kcal/mol, 3-5 kcal/mol, 5-10 kcal/mol, 10-15 kcal/mol and >15 kcal/mol). Then, the average norms (see definition in Section 3.3) of the wavefunctionand PES-based NACMEs in each group were calculated (see Table 1). Firstly, the average norm gets smaller when the energy gap Δ is increased, which is consistent with the discussion above(see supplementary materials). In addition, as mentioned above, the approximate PES-based algorithm fails to give good enough norm of NACMEs at the optimized CI (363.4 Bohr −1 vs. 240.9 Bohr −1 ). As to the other groups, the deviations of the average norms between the wavefunction-and PES-based NACMEs becomes much smaller with the increasing energy gap. For instance, the average norms of the wavefunction-and PES-based NACMEs were estimated to be 78.2 and 80.4 Bohr −1 when the energy gap between S0 and S1 was within 0.17-1 kcal/mol. The relative deviations for the other groups, i.e., 0.17-1 kcal/mol, 1-3 kcal/mol, 3-5 kcal/mol, 5-10 kcal/mol, 10-15 kcal/mol and > 15 kcal/mol, were calculated to be 2.8%, 3.6%, 3.8%, 4.0%, 3.6%, and 3.4%, respectively (see the last column in Table 1).  To give insight into the accuracy of the PES-based NACMEs, we separated the structures near the optimized CI structure into several groups according to their energy gaps ∆E S 0 S 1 (i.e., 0.17-1 kcal/mol, 1-3 kcal/mol, 3-5 kcal/mol, 5-10 kcal/mol, 10-15 kcal/mol and >15 kcal/mol). Then, the average norms (see definition in Section 3.3) of the wavefunction-and PES-based NACMEs in each group were calculated (see Table 1). Firstly, the average norm gets smaller when the energy gap ∆E S 0 S 1 is increased, which is consistent with the discussion above(see supplementary materials). In addition, as mentioned above, the approximate PES-based algorithm fails to give good enough norm of NACMEs at the optimized CI (363.4 Bohr −1 vs. 240.9 Bohr −1 ). As to the other groups, the deviations of the average norms between the wavefunction-and PES-based NACMEs becomes much smaller with the increasing energy gap. For instance, the average norms of the wavefunction-and PES-based NACMEs were estimated to be 78.2 and 80.4 Bohr −1 when the energy gap between S 0 and S 1 was within 0.17-1 kcal/mol. The relative deviations for the other groups, i.e., 0.17-1 kcal/mol, 1-3 kcal/mol, 3-5 kcal/mol, 5-10 kcal/mol, 10-15 kcal/mol and > 15 kcal/mol, were calculated to be 2.8%, 3.6%, 3.8%, 4.0%, 3.6%, and 3.4%, respectively (see the last column in Table 1). In order to give more in-depth comparison between the wavefunction-and PES-based NACMEs, the scattering plots for each group with the different ∆E S 0 S 1 values (including all x/y/z components of each atom) are given in Figure 3. It is obvious that: (1) the PES-based NACMEs at the optimized CI have large deviations (red points in Figure 3a); and (2) the PESbased NACMEs are close to the wavefunction-based NACMEs for the structures far away from the optimized CI structure (see Figure 3a-f). In one sentence, these data demonstrate that the PES-based algorithm can give as accurate NACMEs as the wavefunction-based algorithm. Finally, one should also note that NACMEs calculated by both the wavefunctionand PES-based algorithms will become divergent at truly degenerate CIs. So, one can see a large deviation between both the wavefunction-and PES-based NACMEs at the optimized CI structure (energy gap less than 0.17 kcal/mol; see above). In order to give more in-depth comparison between the wavefunction-and PESbased NACMEs, the scattering plots for each group with the different Δ values (including all / / components of each atom) are given in Figure 3. It is obvious that: (1) the PES-based NACMEs at the optimized CI have large deviations (red points in Figure  3a); and (2) the PES-based NACMEs are close to the wavefunction-based NACMEs for the structures far away from the optimized CI structure (see Figure 3a-f). In one sentence, these data demonstrate that the PES-based algorithm can give as accurate NACMEs as the wavefunction-based algorithm. Finally, one should also note that NACMEs calculated by both the wavefunction-and PES-based algorithms will become divergent at truly degenerate CIs. So, one can see a large deviation between both the wavefunction-and PES-based NACMEs at the optimized CI structure (energy gap less than 0.17 kcal/mol; see above).

PES-Based NACMEs with ML Models
Although the PES-based NACMEs give good enough results compared with the wavefunction-based NACMEs, this PES-based algorithm remains expensive because it needs full Hessian matrix elements, which still cannot be analytically calculated for many electronic structure methods and packages nowadays, to our best knowledge. Fortunately, this shortcoming can be avoided by using trained ML models, which can calculate Hessian matrix elements analytically and efficiently, as seen in our recent papers [17]. Thus, we next used a machine learning technique, i.e., the embedding atom neural network (EANN) method [18][19][20], as an alternative to calculate the relevant information of PESs, i.e., energies, forces, Hessian matrix elements, etc.
It should be noted that the ML models trained by the EANN method can provide analytical Hessian matrix elements while this ML method does not need ab initio Hessian matrix elements as input in the training stage [17,19]. The trained ML model reproduces the PESs near the optimized CI structure, calculated with the SA-CASSCF method (Figure 4), very well, which demonstrates that the trained ML model is accurate enough and can be used for the PES-based algorithm to calculate NACMEs. In fact, the previous works have proved that the ML model can give as accurate energies, forces, and Hessian matrix elements as the ab initio method as long as a sufficient quantity and quality data are provided for the ML training [16,17,21].
Although the PES-based NACMEs give good enough results compared with the wavefunction-based NACMEs, this PES-based algorithm remains expensive because it needs full Hessian matrix elements, which still cannot be analytically calculated for many electronic structure methods and packages nowadays, to our best knowledge. Fortunately, this shortcoming can be avoided by using trained ML models, which can calculate Hessian matrix elements analytically and efficiently, as seen in our recent papers [17]. Thus, we next used a machine learning technique, i.e., the embedding atom neural network (EANN) method [18][19][20], as an alternative to calculate the relevant information of PESs, i.e., energies, forces, Hessian matrix elements, etc.
It should be noted that the ML models trained by the EANN method can provide analytical Hessian matrix elements while this ML method does not need ab initio Hessian matrix elements as input in the training stage [17,19]. The trained ML model reproduces the PESs near the optimized CI structure, calculated with the SA-CASSCF method ( Figure  4), very well, which demonstrates that the trained ML model is accurate enough and can be used for the PES-based algorithm to calculate NACMEs. In fact, the previous works have proved that the ML model can give as accurate energies, forces, and Hessian matrix elements as the ab initio method as long as a sufficient quantity and quality data are provided for the ML training [16,17,21]. Similarly, we first compare the average norms of the PES-based NACMEs calculated by the ab initio (SA-CASSCF) method and the trained ML model ( Table 2). In total, the deviation between the ab initio and ML-calculated NACMEs (PES-based algorithm) is very small when the energy gap is greater than 1 kcal/mol except at the optimized CI structure with an energy gap of 0.17 kcal/mol. In detail, the deviation of the average norms at the optimized CI structure is as large as 73.9% (SA-CASSCF vs. ML: 240.9 Bohr −1 vs. 62.9 Bohr −1 ). Moreover, the PES-based NACMEs with the ML model do not perform very well when the energy gap is in the range of 0.17-1 kcal/mol, which gives a deviation of 20.2%. Nevertheless, when the energy gap is larger than 1 kcal/mol, the deviations of the average norms become much smaller. Specifically, the deviations of average norms related to the energy gaps of 0.17-1 kcal/mol, 1-3 kcal/mol, 3-5 kcal/mol, 5-10 kcal/mol, 10-15 kcal/mol and >15 kcal/mol are estimated to be 3.6%, 9.1%, 0.6%, 2.3% and 2.7%, respectively. In terms of the above results, one can see that the PES-based NACMEs calculated by the ML model can also give accurate results, except at certain structures with extremely small energy gaps (i.e., around 1 kcal/mol). Similarly, we first compare the average norms of the PES-based NACMEs calculated by the ab initio (SA-CASSCF) method and the trained ML model ( Table 2). In total, the deviation between the ab initio and ML-calculated NACMEs (PES-based algorithm) is very small when the energy gap is greater than 1 kcal/mol except at the optimized CI structure with an energy gap of 0.17 kcal/mol. In detail, the deviation of the average norms at the optimized CI structure is as large as 73.9% (SA-CASSCF vs. ML: 240.9 Bohr −1 vs. 62.9 Bohr −1 ). Moreover, the PES-based NACMEs with the ML model do not perform very well when the energy gap is in the range of 0.17-1 kcal/mol, which gives a deviation of 20.2%. Nevertheless, when the energy gap is larger than 1 kcal/mol, the deviations of the average norms become much smaller. Specifically, the deviations of average norms related to the energy gaps of 0.17-1 kcal/mol, 1-3 kcal/mol, 3-5 kcal/mol, 5-10 kcal/mol, 10-15 kcal/mol and >15 kcal/mol are estimated to be 3.6%, 9.1%, 0.6%, 2.3% and 2.7%, respectively. In terms of the above results, one can see that the PES-based NACMEs calculated by the ML model can also give accurate results, except at certain structures with extremely small energy gaps (i.e., around 1 kcal/mol).
In addition, Figure 5 compares the SA-CASSCF-and ML-calculated NACMEs (PESbased algorithm). Figure 5a shows that the ML-calculated NACMEs are to some extent far from the SA-CASSCF-calculated NACMEs at the optimized CI structure. This can be understood very well when considering that PESs are discontinuous at these energetically degenerate points (i.e., singularities). In addition, the accurate fitting of PESs at CI points is also extremely difficult for ML techniques due to this discontinuous character. Although the ML-calculated NACMEs are not so accurate at the optimized CI structure, they can give fairly good results for the structures where the energy gap is larger than 3 kcal/mol (Figure 5c-f). In summary, ML-calculated NACMEs provide reasonably accurate results except at truly degenerate points.  Finally, the computational cost of both the wavefunction-and PES-based NACMEs is discussed. As we have mentioned, the calculation of Hessian matrix elements at ab initio level is very expensive, which makes the computational cost of PES-based NACMEs high using ab initio methods. However, the combination of efficient ML models with the PESbased algorithm brings advantages because the calculation of Hessian matrix elements is cheap using ML models. Finally, the computational cost of both the wavefunction-and PES-based NACMEs is discussed. As we have mentioned, the calculation of Hessian matrix elements at ab initio level is very expensive, which makes the computational cost of PES-based NACMEs high using ab initio methods. However, the combination of efficient ML models with the PES-based algorithm brings advantages because the calculation of Hessian matrix elements is cheap using ML models.

The Approximate PES-Based Algorithm for NACMEs
Here, we briefly introduce the approximate algorithm for NACMEs proposed by Baeck and An [14,15]. Diabatic states could be transformed from adiabatic states by a unitary transformation with an adiabatic-to-diabatic (ADT) mixing angle θ [22]: where Ψ i Ψ j represent diabatic states while Φ i Φ j represent adiabatic states, respectively. The index i and j represent the indexes of different electronic states. Meanwhile, the relationship between diabatic potentials V i V j and adiabatic potentials E i E j can be written as follows [11]: Additionally, θ can be expressed by in which Q 0 is a reference coordinate and d ij (Q) = Φ a j ∂/∂Q |Φ a i is a nonadiabatic coupling term. On the other hand, in 1981, Werner and Meyer proposed a Lorentz function d Lo ij (Q) to approximate d ij (Q) as [23] d Lo ij (Q) = in which Q are nuclear coordinates of the current structure while Q c are coordinates of the crossing point. When Q = Q c , one can obtain α = 2 · d Lo ij (Q). However, this approximate method was not developed further until the work by Baeck and An in 2017, [15] who showed that the approximate NACMEs d appr ij can be obtained by in which ∆E ij is the energy gap between electronic states i and j, R represents the coordinates of the atoms. In addition, a more practical formula was adopted by Lasorne et al. and Westermayr et al. [12,16], which is shown below: in which ⊗ means a tensor product. The final approximate PES-based NACMEs are calculated via in which v ij and λ ij are eigenvectors and eigenvalues of the right side of Equation (6) calculated by the singular value decomposition (SVD). Additionally, the signs (i.e., positive or negative) of the PES-based NACMEs are determined by Equation (7) directly.

The Embedding Atom Neural Network (EANN) Method
Nowadays, machine learning techniques are widely used in the research fields of physics, chemistry, biology, and materials [24][25][26][27][28][29][30][31][32][33][34]. In particular, using ML methods for improving computational efficiency is one of the most popular fields [21,[35][36][37][38][39][40][41][42][43]. In this work, the EANN method proposed by Jiang et al. was adopted for training potential energy surfaces of polyatomic molecules because of its available Hessian matrices [20]. In the EANN method, the total energy of a system is written as the sum of all atomic energies E i : (8) in which N atom is the total number of atoms; NN i is the atomic neural network and ρ i is the function of the embedding electron density of the atom i, which is given by the superposition of electron densities of the surrounding atoms. A set of Gaussian-type orbitals centered at each atom is used for calculating the embedding electron density associated with each atom: Ψ α,r s l x l y l z = x l x y l y z l z exp −α|r − r s | 2 (9) in which x, y and z are coordinate components of an electron coordinate vector r, while |r| is the norm of r. α and r s are parameters for determining the radial distribution of the atomic orbitals Ψ, and l x , l y , l z are indicating the orbital angular momentum components with the relationship L = l x + l y + l z . With the definition of these Gaussian-type orbitals, the embedding electron density of the atom i can be calculated as the square of all the linearly combined atomic orbitals from the neighboring atoms. Then, the individual electron density from Gaussian-type orbitals is written as in which r ij is the distance between the atoms i and j, N neighbor is the number of the neighboring atoms of the atom i within the sphere with the radius equaling the cutoff radius r c , and c j is the expansion coefficient of the Gaussian-type orbitals. Finally, the energy for each atom can be obtained from the atomic neural networks by taking ρ as its input vectors. Moreover, in order to decay the interactions between the central atom with neighbor atoms to zero smoothly, a cutoff function f c (r) is used in the EANN method: (11) in which r c is the cutoff radius. Note that the same elements share the same atomic neural network. After obtaining the input vectors of the EANN method, the atomic neural network models are trained via minimizing the loss function δ(w, b) by the extreme learning machine (ELM) and Levenberg-Marquardt (LM) algorithm (ELM-LM) developed by Jiang et al. [44]. In addition, to improve the efficiency of the ML model training processes and reduce the size of the training data set, the deviations between ML-predicted and SA-CASSCF-calculated atomic forces are added into δ (w, b). Thus, the loss function δ(w, b) is expressed as follows: in which w and b are the weight and bias parameters of atomic neural network models, and N data is the size of the training data set. η is a parameter to control the importance of deviations of atomic forces during the training process. E EANN a and E REF a are energies of the ath configuration in the training data set calculated by EANN models and the ab initio method, while F EANN a,b,c and F REF a,b,c are the forces of the bth atom in c direction for the ath configuration, evaluated by EANN models and the ab initio method, respectively. The accuracy and efficiency of the EANN method for both isolated and periodic systems have been demonstrated in previous works [18][19][20].

Definition of the Average Norm of Matrices
In this work, the average norms of NAC matrices are calculated. The norm of a matrix |A i | is calculated as follows: (13) in which A i,nm is the matrix element of A i . The average norm is calculated by dividing the sum of the norms of A i (I = 1 − n)

Computational Details
In this work, a CH 2 NH molecule is used for demonstrating the accuracy of the approximate PES-based algorithm of NACMEs. To make direct comparisons, the wavefunctionbased NACMEs of CH 2 NH are also calculated at the same SA-CASSCF level. The approximate PES-based NACMEs are calculated using the PES information (energies, gradients, and Hessian matrices) provided by the SA-CASSCF method or the trained ML model, respectively. In this work, we only consider NACMEs between S 0 and S 1 of CH 2 NH. An active space including two electrons and two orbitals, as well as 6-31G* basis sets, are adopted for all the SA-CASSCF calculations. All the SA-CASSCF calculations are performed with OpenMolcas-v18.09 [45,46]. The S 1 /S 0 conical intersection structure of CH 2 NH is optimized by Gaussian16 [47].
In ML model training processes using the EANN method, the atomic neural networks with two hidden layers of 40 and 50 nodes are adopted as the atomic neural network models for C, H, and N elements. The cutoff radius r c is set to 9 , which is large enough to take all other atoms as the neighbor atoms of the center atom. The maximum orbital angular momentum L is set to 2, which is suitable for training models of PESs of CH 2 NH.
On the other hand, β = 0.2 and ∆r s = 1 Å leads to α = β (∆r s ) 2 = 0.2 Å −2 , which are used for determining the Gaussian-type-orbital input vectors in Equation (10) for all atomic neural networks in this work. In addition, η is set to 1 during the model training processes for balancing the importance of deviations of energies and atomic forces in the present work. More details of data collection and training process for the ML models can be found in our previous works [17,21,48]. The purpose of this work is to prove the accuracy and efficiency of the approximate PES-based method for calculating NACMEs. CH 2 NH is chosen as a test system because its excited-state properties have been reported and studied, and it has only 5 atoms [48][49][50]. Thus, the computational cost is affordable for calculating Hessian matrix elements numerically at the SA-CASSCF level.

Conclusions
In this work, we have implemented an approximate algorithm for calculating NACMEs of polyatomic systems based only on PES information, i.e., energies, gradients and Hessian matrix elements. The advantage of this algorithm is that it only demands the information of PESs without any wavefunction information, which is very suitable for machine learning techniques or low-scaling energy-based fragment methods. The approximate algorithm is used for calculating NACMEs with the ab initio SA-CASSCF method and the ML model. The results show that the PES-based algorithm for NACMEs, whether using the ab initio method or the ML model, performs very well except at the CI structure. The present work demonstrates the accuracy of the approximate PES-based algorithm for NACMEs of polyatomic systems. It encourages the combination of this approximate PES-based algorithm with electronic structure methods that do not code wavefunction-based algorithms for NACMEs. Most importantly, this PES-based algorithm provides a good opportunity to calculate NACMEs for energy-based fragment methods and ML models, in which the computational cost of Hessian matrix elements is significantly reduced, as demonstrated in our recent work, and in particular for ML models [17]. Finally, this PES-based algorithm for NACMEs in combination with either low-scaling energy-based fragment methods or efficient ML models also brings an economic but accurate nonadiabatic dynamics simulation method to investigate the nonadiabatic process of large systems.