Three-Dimensional Diabatic Potential Energy Surfaces of Thiophenol with Neural Networks†

Three-dimensional (3D) diabatic potential energy surfaces (PESs) of thiophenol involving the S0, and coupled 1ππ∗ and 1πσ∗ states were constructed by a neural network approach. Specifically, the diabatization of the PESs for the 1ππ∗ and 1πσ∗ states was achieved by the fitting approach with neural networks, which was merely based on adiabatic energies but with the correct symmetry constraint on the off-diagonal term in the diabatic potential energy matrix. The root mean square errors (RMSEs) of the neural network fitting for all three states were found to be quite small (<4 meV), which suggests the high accuracy of the neural network method. The computed low-lying energy levels of the S0 state and lifetime of the 0 state of S1 on the neural network PESs are found to be in good agreement with those from the earlier diabatic PESs, which validates the accuracy and reliability of the PESs fitted by the neural network approach.

incidently, in which nonadiabatic couplings have to be considered and a single PES is not adequate anymore in the description of nuclear dynamics. As a result, the multi-state coupled PESs are necessarily demanded in the accurate simulation of the nonadiabatic dynamics processes.
In particular, the nonadiabatic effects are introduced in the simulation by means of the multiple adiabatic PESs as well as the nonadiabatic couplings, or the diabatic coupled PESs alternatively [2]. Although two treatments in different (adiabatic and diabatic) representations are identical numerically [3,4], the diabatic approach based on the diabatic PESs is preferable due to the fact that the diabatic potential energy matrix elements are all smooth in the whole configuration space of the nuclei and no divergence problems exist in the calculations [5]. Nevertheless, the rigorous diabatic PESs are built up from both the adiabatic energies and nona-diabatic couplings [6,7], which is usually expensive and challenging. And the adiabatic and diabatic representations are connected by a unitary transformation [2].
Thiophenol has received increasing attention in the past several decades owing to its important role of the chromophores in the biological molecules [39−42]. The π * ←π transition in thiophenol is responsible for the strong ultraviolet (UV) absorption [41], and the photochemistry of thiophenol is dominated by the 1 πσ *mediated photodissociation due to the repulsive nature of the πσ * state [41,42]. Since there are 33 degrees of freedom (DOFs) in thiophenol, it is very challenging to construct the full dimensional PESs [43], especially using the high-level ab initio method. There have been two reduced diabatic PESs for thiophenol [44,45], including the two-dimensional (2D) [44] and threedimensional (3D) [45,46] DOFs, which can be used in the photodissociation dynamics calculations. In previous 3D diabatic PESs of thiophenol [45,46], the diabatization of PESs was achieved using the linear-vibronic coupling method [47,48], in which the off-diagonal terms of the diabatic PESs were determined solely by the adiabatic electronic energies. And the PESs were finally fitted by the spline interpolation method. Based on the constructed 3D diabatic PESs, the calculated lifetimes of the S 1 vibronic states of thiophenol were found to be in good agreement with experiment [46], which validated the reliability of the diabatic PESs. In this work, we constructed a new 3D diabatic PESs of the ππ * and πσ * states of thiophenol using the NN method based on the ab initio energy points and some more sampled points on previous 3D diabatic PESs [45]. We focus here on the strongly coupled conical intersection region between the 1 ππ * and 1 πσ * states, which locates nearby the Franck-Condon region and is responsible for the lifetimes of the S 1 vibronic states. Thus, only the coupling between the 1 ππ * and 1 πσ * states was considered. The diabatization was achieved by the NN fitting with the off-diagonal term constrained by the correct symmetry of the system.

A. Neural networks
Since the diabatic potential energy matrix elements are smooth and continuous as a function of nuclear coordinates, in this work we used the NN functions to represent each element of diabatic potential energy matrix.
A multilayer forward-feed NN is built up of an input layer, one or more hidden layers, and an output layer. For a given layer (i+1), the form of the kth neuron can be written as [21,49] y i+1 in which N i is the number of neurons in the ith layer, w i+1 k,j is the weights connecting the jth neuron in the ith layer and the kth neuron in the (i+1)th layer, and b i+1 k is the biases of the kth neuron in the (i+1)th layer. f i+1 is the transfer function for the (i+1)th layer. It can be readily seen that the output of one layer is the input of the next layer. In our calculations, there are two layers in the hidden layer. The input vector includes three variables of the coordinates: where R, θ, and φ (inÅ, degree, degree) denote S−H stretch, C−S−H bend, and C−C−S−H torsion coordinates of thiophenol.

B. Diabatization with neural networks
In this work, the diabatic PESs of the coupled excited states 1 ππ * and 1 πσ * were fitted by the NN method, and the single PES of the ground state 1 ππ was fit-ted individually. The reliability of the new NN diabatic PESs was examined by comparing the calculated lifetime of the S 1 vibronic state with that on previous diabatic PESs. The topography of the 1 ππ * and 1 πσ * PESs in Franck-Condon region is responsible for the lifetimes of the S 1 vibronic states, since the S 1 /S 2 conical intersection on the PESs of thiophenol locates at a short distance (R=2.8 bohr). The rigorous construction of the diabatic PESs is based on the adiabatic energies and nonadiabatic couplings, which is very challenging for thiophenol due to quite many DOFs. Following earlier work [45,46], here we build up the diabatic PESs merely based on the adiabatic energies with NN fitting method.
The 2×2 diabatic potential energy matrix for the 1 ππ * and 1 πσ * states can be written as Each matrix element is represented by a NN function in the fitting procedure. And the symmetry of each element has to be correctly considered. At the planar geometry (φ=0), the 1 ππ * (V 11 ) and 1 πσ * (V 22 ) states are of A ′ and A ′′ symmetries, respectively. But the symmetries of two states exchange at φ=90 • . As a result, the coupling term V 12 carries the A ′ ×A ′′ =A ′′ symmetry. To correctly constrain the symmetry, the NN function for the coupling term V 12 was multiplied by a factor sin2φ. Also the PESs are assumed to be symmetric with φ=90 • in the range of [0 • , 180 • ] [46]. In the fitting procedure, the weights and biases of the NN functions are determined by minimizing the residual sum of squares: in which N dat is the number of the energy points in the fitting. E n,1 and E n,2 are the adiabatic energies of S 1 and S 2 states, which can be easily obtained by diagonalizing the diabatic potential energy matrix in Eq.(3). The Levenberg-Marquardt algorithm [50] was used to iteratively optimize the weights and biases of the NN functions and the Jacobian matrix was calculated numerically by the central-difference algorithm with the displacements of 10 −5 .

III. RESULTS AND DISCUSSION
In our calculations, 3952 adiabatic energy points were chosen from previous 3D diabatic PESs [46] to fit the ground state and the diabatic matrix elements of the 1 ππ * and 1 πσ * states of thiophenol. 16  In the NN fitting procedure, the transfer functions in the first and second hidden layers are the hyperbolic tangent function (f (x)=tanh(x)) and is a linear function in the output layer. 20 neurons in each of the two hidden layers were used in the NN fitting of the single S 0 PES and the diabatic PESs of two coupled 1 ππ * and 1 πσ * states. FIG. 1 shows the root mean square error (RMSE) as a function of the iterative steps in the fitting for the S 0 , S 1 , and S 2 states. It can be readily seen that the RMSEs quickly decrease below 10 meV for the two sets of the NN fitting. The convergence for a single PES fitting of S 0 is achieved within 800 iterations and the RMSE is 2.1 meV. For the NN fitting of the diabatic PESs of 1 ππ * and 1 πσ * states, it converges quickly with 600 interactions and the RMSEs are 2.9 and 3.8 meV for the S 1 and S 2 states, respectively. FIG. 2 displays the fitting error distributions for the S 0 , S 1 , and S 2 states. The distributions of the errors are quite balanced in the energy ranges (∼5 eV). The largest fitting errors for the S 0 , S 1 , and S 2 states are found to be 13.9, 25.9, and 27.4 meV, respectively. Since the range of the energy points in the NN fitting is not large, the weight factors for the energies in the NN fittings are not included.
To examine the accuracy of the NN PES of the S 0 state, we calculated the vibrational energy levels of the excited states for the three vibrational modes (C−C−S−H torsion, C−S−H bend, and S−H stretch). The Lanczos method was used to obtain the 3D vibrational energy levels and the grids in the calculations were the same as those described in Ref. [45]. As shown in Table I the C−S−H angle fixed at the equilibrium (θ=96.2 • ) of S 0 state. As expected, both the NN PESs of the single S 0 and the coupled S 1 and S 2 states are fitted very well. The S 1 /S 2 conical intersection at φ=0 • is located at R=2.795 bohr on the diabatic PESs, which is in excel-  1 , v 2 , and v 3 denote the quantum numbers of the S−H stretching, C−S−H bending, and C−C−S−H torsional vibrational modes, respectively. lence with the location (R=2.786 bohr) of the conical intersection on previous diabatic PESs [46]. FIGs. 4 and 5 display the 3D plot of the adiabatic PESs of the S 1 and S 2 states as a function of R and φ with θ=96.2 • and of θ and φ with R=2.80 bohr, respectively, fitted by the NN method, and the spline PESs from Ref. [46] are also shown for comparison. The locations and shapes of the conical intersections are clearly shown in FIGs. 4 and 5. It is clear that the NN PESs reproduce the topographies of the spline ones especially for the conical intersection regions quite well, which validates the high accuracy of the NN approach.
In FIG. 6, the 1D NN diabatic potential energy curves along R for four different φ values (0 • , 30 • , 60 • , and 90 • ) of the 1 ππ * and 1 πσ * states (θ=96.2 • ) are compared to the spline ones. It can be seen that the 1 ππ * state is very close for the NN and spline PESs. For the 1 πσ * state, two sets of the PESs are similar with some differences for longer R>3.0 bohr. The differences should be caused by two aspects. One is different diabatization methods utilized for two sets of PESs. Furthermore, the coupling between the 1 ππ and 1 πσ * states at longer R was considered in the spline PESs [46], while it was not included in our NN PESs. FIG. 7 displays the 3D plot of the diabatic PESs of the 1 ππ * and 1 πσ * states (θ=96.2 • ) along R and φ fitted by the NN method. It shows that the NN and spline diabatic PESs are quite similar as discussed above. FIG. 4 3D plot of the adiabatic PESs of the S1 and S2 states with θ=96.2 • fitted by the NN method and the spline PESs in Ref. [46]. 1 ππ * and 1 πσ * states are well represented by the NN approach. On the other hand, the diabatization of previous spline PESs was achieved by the linear-vibronic coupling method [46], while the NN diabatic PESs in this work were done by the fitting method with the correct symmetry constraint on the off-diagonal term. Good agreement between two sets of the PESs indicates the reasonable treatment of the nonadiabatic coupling of the conical intersection region of the 1 ππ * and 1 πσ * states by two diabatization approaches.

IV. CONCLUSION
To summarize, in this work we report the NN method to calculate the single S 0 PES and the diabatic PESs of two coupled 1 ππ * and 1 πσ * states of thiophenol. The potential energy points were chosen from previous spline PESs [46] in the fitting, and no ab initio calculations were done in this work. As expected, the NN method behaves quite well on the fitting of the single PES of S 0 and the RMSE is 2.1 meV. The calculated low-lying energy levels of the S 0 state on the new NN PES are found to close to those on previous PES, which suggests the high accuracy of the NN method. Specifically, for the the NN fitting of the diabatic PESs of 1 ππ * and 1 πσ * states, the adiabatic energies were merely in- cluded. The nonadiabatic coupling was not included but with the correct symmetry constraint on the offdiagonal term. The RMSEs are 2.9 and 3.8 meV for the S 1 and S 2 states, respectively. The energy and lifetime of the ZPE state of the S 1 state were calculated on the new NN diabatic PESs, which are in good agreement with previous spline diabatic PESs, validating the reasonability of the diabatic PESs constructed by the NN diabatization method.

V. ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China (No.22073073). Changjian Xie thanks the Startup Foundation of Northwest University. The Double First-class University Construction Project of Northwest University is acknowledged.