Molecular Dynamics Investigations of Ozone on an Ab Initio Potential

ARTICLE pubs.acs.org/JPCA

Molecular Dynamics Investigations of Ozone on an Ab Initio Potential Energy Surface with the Utilization of Pattern-Recognition Neural Network for Accurate Determination of Product Formation Hung M. Le,* Thach S. Dinh, and Hieu V. Le Faculty of Materials Science, College of Science, Vietnam National University, Ho Chi Minh City, Vietnam 750000

bS Supporting Information ABSTRACT: The singlet
triplet transformation and molecular dissociation of ozone (O3) gas is investigated by performing quasi-classical molecular dynamics (MD) simulations on an ab initio potential energy surface (PES) with visible and near-infrared excitations. MP4(SDQ) level of theory with the 6-311g(2d,2p) basis set is executed for three different electronic spin states (singlet, triplet, and quintet). In order to simplify the potential energy function, an approximation is adopted by ignoring the spin
orbit coupling and allowing the molecule to switch favorably and instantaneously to the spin state that is more energetically stable (lowest in energy among the three spin states). This assumption has previously been utilized to study the SiO2 system as reported by Agrawal et al. (J. Chem. Phys. 2006, 124 (13), 134306). The use of such assumption in this study probably makes the upper limits of computed rate coefficients the true rate coefficients. The global PES for ozone is constructed by fitting 5906 ab initio data points using a 60-neuron two-layer feed-forward neural network. The meanabsolute error and root-mean-squared error of this fit are 0.0446 eV (1.03 kcal/mol) and 0.0756 eV (1.74 kcal/mol), respectively, which reveal very good fitting accuracy. The parameter coefficients of the global PES are reported in this paper. In order to identify the spin state with high confidence, we propose the use of a pattern-recognition neural network, which is trained to predict the spin state of a given configuration (with a prediction accuracy being 95.6% on a set of testing data points). To enhance the prediction effectiveness, a buffer series of five points are validated to confirm the spin state during the MD process to gain better confidence. Quasi-classical MD simulations from 1.2 to 2.4 eV of total internal energy (including zero-point energy) result in rate coefficients of singlet
triplet transformation in the range of 0.027 ps
1 to 1.21 ps
1. Also, we find very low dissociation probability up to 2.4 eV of internal energy during the investigating period (5 ps), which suggests that dissociation does not occur directly from the singlet ground-state, but it involves the excited triplet-state as an intermediate step and requires more reaction time to occur.

I. INTRODUCTION During the past recent years, noticeable attention has been given to the chemistry of ozone gas (O3) due to its important role in environmental issues. The molecular absorption of ozone varies from a wide range of bands that consist of the near-infrared Wulf band up to the extension of visible and ultraviolet bands.1 It is of importance to notice the absorption of those ultraviolet bands of ozone that provide an indispensible vitality to life on Earth. The destruction of ozone nowadays due to air pollution thus leads to a potential serious hazard of living species in the regions near the Earth's surface, and therefore has become a major concern for many environmentalists and research scientists. In its most stable ground state, the electronic spin of ozone adopts a singlet configuration and a C2v structural geometry as reported by Tanaka and Morino.2 Configurations of the higherenergetic spin state, which is triplet, have also been the focus of many theoretical and experimental scientists during the past few decades. There are three low-lying triplet states observed and reported by Braunstein and Pack in an early experimental study,3 which are 3B2, 3A2, and 3B1. In that study, a photodissociation r 2011 American Chemical Society

spectrum was applied to the triatomic molecule in the 10 000
22 000 cm
1 band, namely, the Wulf
Chappuis band, and the lowest excited triplet states were observed. In a subsequent study, a computational effort was reported that discussed the spin
orbit interaction of the above low-lying triplet states and the O2+O dissociation threshold.4 A noticeable effort in the early era of the approach to investigate the available triplet states is the contribution to develop a potential energy surface (PES) for the ground-state ozone and its eight low-lying triplet states.5 Tsuneda et al.6 reported the energy levels of various triplet structures using the expensive complete-active-space calculations and multireference Møller
Plesset perturbation. More importantly, the theoretical calculations predicted that no triplet states were found to be located below the molecular dissociation threshold of ozone to form O2 and radical O. Another wellestablished study of the low-lying triplet state of ozone has been reported by Borowski et al.7 Received: July 11, 2011 Revised: August 24, 2011 Published: September 05, 2011 10862

dx.doi.org/10.1021/jp206531s | J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A The ozone singlet
triplet transition has been the subject of experimental and theoretical studies8 for decades. In a theoretical approach reported by Minaev and Agren,9 the energy barriers for singlet
triplet transitions were computed based on multiconfiguration wave functions, and explained the absence of 3B2 state absorption in theory. In a subsequent study, Minaev and Khomenko10 conducted a multiconfiguration self-consistent field method on complete-active-space calculations to reinvestigate ozone singlet
triplet transitions with the utilization of quadratic response functions. It is also implied that the triplet state plays an important role as an intermediate step in the 3 dissociation of singlet-state O3 to O2(3∑
g ) and O( P) (and vice 3
3 versa, the recombination of O2( ∑g ) and O( P)). Grebenshchikov et al.11 provided an excellent study on ozone photodissociation in which they conducted both quantum and classical dynamics, and they focused on the nonadiabatic coupling of singlet and triplet states. However, in that study, the construction of the ab initio PES and its fitting quality were not clearly indicated. In many chemical reactions (in this case, a simple transformation from singlet to triplet), there exists a quantum electronic switch among two or multiple spin states, which leads to a change in the electronic configuration structure of the system. Consideration of such nonadiabatic processes has motivated many theoretical chemists to design a method that treats the classical systems using a quantum mechanical assumption. Consequently, the surface hopping model developed by Tully and co-workers12 since the early 1970s has been employed to study many molecular systems. The molecular dynamics (MD) investigation of ozone presents a particularly difficult challenge because of the complexity of three spin states of ozone. In reality, a “spin-orbit coupling” model is a more realistic scheme to deal with the problem of this type. However, for simplicity of our calculations, we adopt an assumption that was previously utilized by Agrawal and co-workers13 to treat the SiO2 system. It is very important to imply at this point that we treat the spin switch as an instantaneous process that favors the lowest energy level, which is in the same manner that Agrawal et al.13 handled the SiO2 case. In most cases, including our study of ozone singlet
triplet transformation and molecular dissociation, the utilization of this assumption would make the computed rate coefficient upper limits the true rate coefficients. The construction of such PESs that fully describes all available geometric configurations and favors the lowest energy level among three spins is required prior to MD simulations. The construction procedure of such a potential function will be discussed in a later section. There are numerous efforts of PES development for ozone that involve different levels of theory and accuracy. The oldfashion empirical PES’s and classical MD simulations have been utilized to investigate the transition (and dissociation) process of ozone with limited accuracy.14 Xie et al.15 performed several different types of calculations using various levels of computationally expensive ab initio theories, and reported two different PES fits for the near-equilibrium ground-state ozone and lowlying singlet
triplet transition. Nowadays, the methods for PES fitting have been developed in advance along with the growth of machine learning algorithms. The interpolating moving least-squares (IMLS) method16 has been developed for years and applied to a wide variety of molecular systems. In this method, the fitting accuracy depends on the parameters in use and the size of the “cut-off radius”, which determines the size of the database in consideration.

ARTICLE

Another method has been developed by Collins and co-workers, namely, the Shepard interpolation method.17 This method is designed to interpolate the potential energy of a given configuration from a series of data points surrounding it as a combination of Taylor’s expansions, which involve the zeroth-, first(gradient), and second-order (Hessian) terms, and thereby provide excellent accuracy during gradient analysis for differential equation integration. The Shepard interpolation method has a noticeable advantage in that the number of data points to construct a PES is much smaller than that required by other methods; however, the Hessian matrix and gradients are required to be computed at every configuration in the database. The last method in our discussion is the neural network method,18 which employs a machine learning algorithm to perform generalized fits for potential energy based on geometric configuration parameters. This method has been proven to be robust and powerful in numerical analysis problems because parameters are automatically optimized by the training algorithm using some strict convergence criteria, and the mathematical formation of the neural network is of ease to interpret. Recently, the neural network method has been successfully utilized to study many molecular systems, whose examples include vinyl bromide (CH3CH2Br),19 nitrous acid (HONO),20 hydrogen peroxide (HOOH),21 and BeH3.22 Pukrittayakamee et al.23 have made the neural network more powerful in gradient fitting by using both the potential energy and its first derivatives as fitting objectives. This technique was applied to the H + HBr illustration system and resulted in excellent fitting accuracy. In this research, we aim to execute two tasks: (1) The construction of a neural network PES for the ozone system that efficiently describes the three multiplicity states (the singlet ground-state, the triplet excitation-state, and the quintet dissociation-state) as shown in the following chemical reactions: O3 ðX1 A 1 Þ f O3 ð3 A 2 Þ

ð1Þ

O3 ðX 1 A 1 Þ f O3 ð3 B1 Þ

ð2Þ

O3 ðX1 A 1 Þ f O2 ð

3

∑ g
Þ þ Oð3PÞ

ð3Þ

Following the first task, our second objective (2) is subsequently executed when the singlet
triplet transition as well as molecular dissociation at various levels of excitation energy introduced into three vibrational modes are investigated. This second task is performed using the quasi-classical MD simulation described by Raff and Thompson.24

II. COMPUTATIONAL DETAILS In this work, we utilize fourth-order Møller
Plesset perturbation theory25 calculations with single, double, and quadruple excitations26 (MP4(SDQ)) on the O3 system using the 6-311g(2d,2p)27 basis set. This feature is provided in the Gaussian 03 suite of program,28 and the job is executed on our workstation facility with an Intel i7-990X processor. In this three-body system, we define a set of three internal coordinates, R1, R2, and θ, as the two bond lengths and the angle between the two bonds, respectively. Since the aim of this research project is to construct a generalized PES that describes both the singlet
triplet transformation and molecular dissociation, we perform a grid scan of those three parameters in reasonably and sufficiently wide ranges. As mentioned above, the equilibrium 10863

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A

ARTICLE

Figure 2. A two-layer feed-forward neural network model. Figure 1. A singlet
quintet hopping illustration. In this example, the
, and we O
O bond is calculated at various distances from 1.3 to 2.3 Å see clearly that the O
O bond is broken before 1.7 Å
. Note that this plot does not represent the reactive potential barrier of ozone dissociation.

singlet state of O3 adopts a C2v geometric configuration, and therefore has two equivalent O
O bonds. In our calculations, R1 is
to 1.657 Å
, R2 is varied in the range of varied in the range of 1.057 Å
to 2.157 Å
, while θ is varied in the range of 73.7° to 157.7°. 1.057 Å It should be noted that these ranges sufficiently cover the required configuration hyperspace, and satisfy the required numerical criteria to construct a full-scale PES for the reactions of interest. For each configuration point that we generate in the defined ranges, three different calculations are to be performed at three different spin states, which include the singlet, triplet, and quintet states. During the dynamics process, with a certain level of excitation energy, it is possible for the molecule to “cross” to another spin state that is distinct from its current state. The generalized energetic scheme that handles the spin switching is considered with simplicity in our case, i.e., the lowest MP4 energy and its corresponding multiplicity will be chosen among three states, and the information (of energy and spin) will be stored in the database. In our approach, an important assumption has been made, i.e., that spin switching is very rapid, such that the current spin state changes instantaneously if the energy level of another spin is favorably lower. This statement is made to be a primary phenomenon to construct the PES for ozone in our problem, and it was also previously employed to treat the SiO2 system in a theoretical study by Agrawal et al.13 An illustration of this process is presented in Figure 1 with an illustration of the singlet-quintet switching when O3 undergoes a dissociation process. In this study, it is simple to sample the data in hyperspace by performing a grid scan in three-dimensional hyperspace. In particular, we believe it is beneficial to discuss some methods that have been employed in the sampling procedure in hyperspace. In several previous studies that describe more complicated systems (four-body and six-body systems), it is very important to develop a sampling scheme to select the geometric configurations wisely and effectively. Over the past recent years, at least three methods have been proposed and successfully demonstrated in PES development. The first traditional method in our discussion here is used to construct the PES of the vinyl bromide system,19 and the MD technique is employed to collect data points in the multidimensional hyperspace. In the very early stage

of this approach, the first database for a complex molecular system is initialized by performing MD on an empirical PES, then a temporary neural-network PES is constructed. Subsequently, the data are selectively generated and added to the database until some convergence criteria are satisfied. Such procedure has been termed “novelty sampling”.29 The second method is developed by Agrawal et al. and does not require an initial empirical PES to be preconstructed nor available in the literature, but the database is initialized by performing direct dynamics.30 The last method discussed in this paper is the “gradient sampling” method developed by Le et al.,21 which allows configurations points to be statistically selected with a guarantee of accurate gradient fitting. The advantage when using these methods is clear, i.e., the configuration space is well-covered. However, in a simple case such as ozone with only three parameters, the utilization of such method is unnecessary, and we can still achieve the appropriate accuracy by executing the sampling using a simple grid scan.

III. CONSTRUCTION OF THE POTENTIAL ENERGY SURFACE After ab initio calculations of the generated configurations are finished, we obtained a set of 5906 data points with three energy levels (corresponding to singlet, triplet, and quintet spin states). Among the three energies, we choose the lowest energy values and its corresponding multiplicity and store them in the database in accordance with the structural data (two bond distances and the bending angle). We then develop a neural network function that fits the energies to the inputs to produce a global PES, and a pattern-recognition neural network for the multiplicity identifying purposes. The development of these two neural networks is described in detail in this session. 1. The Construction of a Two-Layer Neural Network as a Global Potential Function. At this point, we have the potential

energies for all configurations in hand, and a mathematical function is required to generalize all these points. Among the well-established methods, the neural network method is chosen due to its robustness, stability, and affordable computational expense. In this work, a twolayer feed-forward neural network18 (Figure 2) incorporated in MATLAB31 is employed to fit the database. In a two-layer neural network, a mathematical function is used to comprise the data at the output point of each layer. It is shown by Hornik et al.32 that a multilayer neural network that utilizes a sigmoid transfer function (which is the logistic or hyperbolic 10864

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A

ARTICLE

Table 1. Weight, Bias Value of the Neural Network Hidden Layer and Output Layera w1

b1

w2


3.586180


0.507474

1.556666

4.365270


0.889676

0.013313 0.142110


2.980783
0.530101

0.461965 1.596286

2.589789
2.118280


6.979013
2.930647

4.178026

0.914849


0.749920


2.139668

0.540719

0.036177


3.783701

0.589491

3.231187

3.702781


1.471031


1.477788

0.290911

0.261044


4.694358


0.952020

1.917824

1.984804


1.926846

3.899667
3.494932

3.503966


0.662238


1.755436


1.399616


1.182778


0.038761

6.292202

6.722517

0.231634

5.761412 2.789320


16.765492
1.069637


1.999299 1.606116


14.983415
4.975135

0.171531 1.030260 1.712624

4.110416


8.620036


1.537525


8.515514

1.697334

2.121878


0.250906

1.165366

1.309999

3.813647


8.444381


1.420702


8.171693


2.003579

8.286917

0.075225

2.883465


1.715616


0.068968

4.597119

4.997342


8.010269

5.183884


0.038436

1.481080


0.853177

1.424392


1.285252


3.628971

5.837761 0.224486


3.130490
3.465725

0.494304 1.155868

6.820777
1.710902

1.853528
1.647621

3.397780


2.161648

0.329065


1.943898


0.618744

6.089249


2.956209

4.407074


2.528733


0.122228

6.349488

1.243148


2.537307


3.022796


0.158568


0.917570

0.572984


1.061158

0.817996


7.289877
1.987767

5.669861


3.123531

0.540623

6.679986

27.293490


1.446717

0.455340

13.719343

0.094797

1.036277
0.930713


3.173065 7.478319


0.706004
3.431481

0.131221
1.241259


0.803490 0.399432


1.008848


0.860771

0.170494

0.098127

14.139998


0.349061

4.290373


1.656082

0.491423

1.261193


2.023635


2.010816


5.908378


0.210538


2.330207

3.271504


0.622536


1.723671


1.245427

3.706852


2.091480


2.065109


5.867482


0.200003

2.417960

22.343306

2.662149


0.842448

12.849888

4.978558


0.854830
0.644351

1.261765
0.291914

2.934264 2.057916


0.413686
0.274056

0.514588
0.818341
17.347112


1.162524

0.325412

0.042953


0.587039

1.917997


0.502696

0.125991

0.825124


5.960402

5.065570

5.082629

9.907110

3.781791


0.034063

0.165371


3.851599

1.861966

1.521588


0.655345

0.453965


1.584599


1.741619


1.132183


0.712761


17.997331


1.395359

0.396536


10.097165

1.080523

0.161556
0.572108


3.404615
0.312341

1.003647 2.353646


3.058539 1.003286


1.375830
2.194331
0.023259

10.337297


4.813570


2.338565

10.059357


2.978025


87.367281

0.952322


51.968073

7.563908

3.269906

0.892203


0.142491

2.376933


0.922522
0.129576


6.347758

6.136418


0.070433


5.165272

2.332873


4.776185


5.058059

2.711062

0.238507


0.906192

1.718923

2.640756


1.969971

5.167740


2.630881

0.647547


1.855399


1.958756


0.394847


2.762073


81.626516

0.894195


48.645140


7.556081


0.328277


0.173489

2.270328

0.990881

2.320027

1.020820

1.019707


3.131851

3.570778

6.973449

Table 1. Continued w1

b1

w2


2.245321
4.974509

3.174784
8.564097


4.902874
0.060144

1.269038
1.150876

1.202540 15.299107


0.971378


1.014220

4.094778


4.277514

3.513516

24.156593

3.016682


0.969102

13.926658


3.964885
8.200755


0.979408

1.894073

2.436964


2.049691

0.001332


2.460593


0.220389


3.457202

8.093127

2.615988

0.405050

0.009801

3.317430


8.506810

a 2 b = 1.459784, r1 and r2 are in the range of [1.056622; 2.156622] (Å
), and r3 (cos of the bending angle) is in the range of [
0.925316; 0.280398] (unitless). The range of the potential energy is [0; 2.499920] (eV).

tangent function) in the first layer and a linear function in the second (output) layer provides a universal approximation to many analytic functions. The choice of such functions is advantageous as it provides smooth and easily differentiable potential energy functions. This statement is validated through a series of studies in which several potential energy problems have been empirically demonstrated,13,19
23,30,33 and the results are mostly revealed with relatively low fitting errors (in terms of both potential energies and gradients with respect to the input parameters). The input vector in this problem comprises three parameters, which includes two bond lengths and the cosine (cos) of the bending angle between two bonds. For simplicity, all of these parameters are denoted as ri. It is particularly advised that those input values need to be scaled in the range of [
1;1] for the unitless guarantee to achieve better fitting accuracy. The scaling formula for the input values is as follows: rscaledi ¼

2ðri
rmini Þ
1 rmaxi
rmini

ð4Þ

where rmax and rmini are the maximum and minimum of the ith input parameter in the database, respectively. Similarly, the output potential energy is also scaled using the same method: Vscaled ¼

2ðV
Vmin Þ
1 Vmax
Vmin

ð5Þ

Vmax and Vmin in the above equation are the maximum and minimum potential energies, respectively. A 60-neuron feed-forward neural network is employed to fit the data in our study. The output of the first layer is produced using the scaled input parameters as demonstrated in eq 6: ni ¼ tanhð

3

w1j, i rscaled ∑ j¼1

j

þ b1i Þ for i ¼ 1, :::, 60

ð6Þ

In the above equation, we denote w1 as the “weight” matrix of size (60  3), and b1, a column vector of size (60  1), as the bias vector. In the first hidden layer, the hyperbolic tangent function (tanh) is used as the transfer function. The outputs of the first layers are then introduced into the second layer to produce the final output, which resembles the scaled potential energy as shown in the following equation: Vscaledo ut ¼

60

∑ w2i ni þ b i¼1

ð7Þ

In this equation, w2 is a row vector of size (1  60) representing the weight values of the second layer, and b is the bias value of the 10865

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A

ARTICLE

Table 2. Bond, Angles, and Wavenumbers of Singlet and Triplet States R 1

MP4(SDQ) singlet (X A1) triplet (3A2) triplet (3B1) MRCI+Q

exp a

Figure 3. (a) Plot of MP4(SDQ) energies versus neural network predicted energies. The majority of data are fitted with low statistical errors, while a small number of data points are not well fitted. Overall, the fitting mean-absolute error is still impressive (as low as 0.0446 eV, or 1.21 kcal/mol). (b) Distribution plot of absolute fitting errors. The majority of errors reside in the low fitting error region, which results in a good mean-absolute error.

second layer. The transfer function in the second layer is just simply a linear function. All required parameters for this neural network are listed in Table 1. When a training job is performed, all data is randomly divided into three subsets: a training set with approximately 80% of the data, a testing set with approximately 10% of the data, and the remaining data is used as the validation set. In the training process, all parameters are optimized to reduce the mean squared error of the training set using the Levenberg
Marquardt algorithm.18 Technically, a training iteration is named as an “epoch”, and the training process may take a maximum of 1000 epochs if the convergence criteria have not been satisfied. In the neural network construction, the phenomenon that a larger number of hidden neurons results in lower statistical is not necessarily true when “over-fitting”34 becomes a major concern if an excessive number of neurons is used. It is advisible to estimate a sufficient number of hidden neurons; however, there are in fact no criteria that have provided a complete guarantee with total satisfaction. Indeed, a very efficient and well-established

θ

ν1

1.257 117.7 756.2

ν2

ν3

1281.6

1340.6

1.314 97.6 1.285 129.6

singlet (X1A1)15a 1.274 116.9 718.8 1105.28 1125.16 triplet (3A2)15b

1.332

triplet (3B1)15b

1.315 123.8

singlet (X1A1)a

1.272 116.8

triplet (3A2)8c

1.345

99.1

98.9

J. C. Depannemaecher and J. Bellet, J. Mol. Spectrosc. 1977, 66, 106.

technique has been utilized to prevent overfitting even when an excessive number of neuron is used. This technique is known as “early-stopping”, which terminates the training process if the mean squared error of the validation set increases in six consecutive epochs (this number is chosen empirically, and may be changed for various purposes). The neural network fit has been executed with a relatively low root-mean-squared error with a total database of 5906 points (including the three sets that we mention above). In some regions, the statistical errors are high, but the overall fitting error is still in good agreement with the MP4 energies. A plot of neuralnetwork-predicted outputs versus MP4 energies is shown in Figure 3a. The mean-absolute error and root-mean-squared error for the fit of all data points in this study are calculated as 0.0446 eV (1.03 kcal/mol) and 0.0756 eV (1.74 kcal/mol), respectively. In Figure 3b, we prepare a plot of absolute error distribution, which clearly shows a majority domination of small absolute errors. It can be seen that there are a few data points that are not wellfitted, and this leads to the domination in root-mean-squared error over the mean-absolute error, and makes the root-mean-squared error almost 1.7 times larger than the mean-absolute error. When making a judgmental analysis for this fit, we have to keep in mind that the electronic configuration of the ozone system is too complicated with the consideration of two excited multiplicity states (triplet and quintet) as well as the singlet ground state in one single potential energy function. Nevertheless, those two reported errors in our study represent an excellent statistical fit for a PES construction that can be compared favorably to the reported fitting errors in the literature.17b,20
22,33a 2. Identifying the Spin State Using a Pattern-Recognition Neural Network. The identification of a configuration during the MD process is very important as we will discuss later in more detail. Therefore, the determination of the ozone spin state for a given geometric configuration becomes a challenging task for us. As shown in Table 2, the structure of singlet ozone (X1A1) and the two triplet structures (3A2 and 3B1) are difficult to distinguish, especially during the MD process because of the small distinction among them. Consequently, a computational pattern recognition technique must be applied to deal with the identification problem. The architecture of the pattern-recognition neural network is of similar type as the structure of the previous neural network. The only distinction between the two neural networks is the transfer function. In most pattern classifications, a threshold function is used at the last output layer to set the final output at one of the multiple defined values. An illustrative example of this 10866

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A

ARTICLE

Figure 4. Three multiplicity states are classified according to the geometric configuration. In some regions, the prediction line is blurred and does not indicate the state clearly; as a result, 95.6% of points are classified correctly. We prevent this source of error by using a buffer series of five points in our multiplicity examination.

application in our study is presented in Figure 4. Overall, when we produce a pattern recognition fit, the spin of 95.6% of a testing data set is predicted correctly.

IV. MD SIMULATIONS In this study, the dynamics of singlet
triplet transformation is the most important process. As mentioned elsewhere in the literature,15 this transformation may also involve the ozone dissociation from the original singlet state. Energetically, the energy levels of two triplet states in this study (3A2 and 3B1) locate above the ground-state dissociation. Theoretically, the potential energy barrier of the ground-state dissociation remains unknown because of the complicated scheme of surface hopping, but it is predicted that the ozone system (originated from the singlet ground state) has to undergo one of the triplet states before reaching a final dissociation (quintet). Therefore, it is very important to detect the existence of the triplet states precisely after the excitation energy is introduced into the system. The quasi-classical MD is executed using the fourth-order Runge
Kutta method with a fixed step size of 1.018  10
16ps. In total, there are 18 differential equations in the system to be solved simultaneously in each step of trajectory calculations. The requirement of total energy conservation is satisfied as the numerical energy is conserved to the sixth decimal place. During the integration, the system multiplicity is monitored at every step by the pattern-recognition neural network with a buffer series of five latest configurations. In the buffer examination, we record the spin states of the five latest structures from the investigated trajectory, and if all of them turn to triplet or quintet, we conclude the reaction to be a singlet
triplet transformation or a molecular dissociation, respectively. This technique allows us to reduce minor prediction errors by the pattern-recognition neural network, and thus enhance the accuracy of reaction time. In the first stage, we need to initialize a random structure with the structural and kinetic variables (Cartesian coordinates and momenta) being statistically randomized in hyperspace. Initially, the optimized equilibrium structure is assigned to the molecule, and the zero-point energy of 0.182 eV is introduced into the system by distributing the corresponding vibrational energy into each mode using the projection method.35 The system is then allowed to vibrate for a randomized period of time (less than 3 times the longest vibrational period or 0.15 ps), and the total

Figure 5. First-order decay plot of O3 singlet
triplet transformation at the internal energy of 2.0 eV.

angular momentum should vanish during the entire integration. After executing the integration, a randomized configuration with randomly assigned momenta (with correlation to zero vibrational state) is well-prepared for MD investigations. The excitation energy is introduced to each vibrational mode also using the projection method. From our MP4(SDQ) calculations, the two triplet states are located at 1.15 eV (3A2) and 1.39 eV (3B1) above the energy of the equilibrium singlet groundstate. We carry out several trajectory calculations with various excitations of three vibrational modes (the total energy is up to 0.687 eV) in order to validate the low-energetic dynamics of ozone. The sets of vibrational quantum numbers for those testing excitations include (0,0,0) (ground state), (1,0,0), (0,1,0), (0,0,1), (1,1,0), (1,0,1), (0,1,1), (1,1,1), (2,1,1), (1,2,1), and (1,1,2). We observe in these cases that no chemical reactions occur according to the energetic scheme provided by the PES. In each excitation, 500 trajectories are integrated, the elapse time for each trajectory is 5 ps, and the resulting reaction periods are recorded. It is observed in this first test that no singlet
triplet transformation or dissociation occurs due to the insufficient excitation energy to activate barrier overcoming. Recall that the zero-point energy for ozone is 0.182 eV, and the maximum internal energy of those excitation test cases (0.687 eV) is far above the vibrational ground state energy. In the later investigation, when a larger amount of excitation energy is introduced into the vibrational modes, we assume that the system behaves in a classical manner, i.e. the excitation energy is spread equally among three modes. In this investigation, we perform MD simulation at different total energy levels from 1.2 to 2.4 eV (including zero-point energy). The process will be terminated if a chemical reaction occurs (which means the multiplicity has changed to triplet or quintet) or 5 ps elapses. To confirm the spin change with high accuracy, we monitor the multiplicity of the systems using the pattern-recognition neural network with a buffer of 5 integration steps as discussed above. For each internal energy level, 500 trajectories are investigated. Kinetically, we observe that the singlet
triplet transformation is a first-order reaction when first-order decay plots for all investigated cases are made based on our results, and the firstorder decay plots are all linear. In Figure 5, a first-order decay plot of the singlet
triplet transformation is shown with excellent 10867

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A

ARTICLE

Table 3. Singlet
Triplet Transformation Rate Coefficients and Probability, and Dissociation Probability at Various Total Internal Energy Levels total energy (including zero-point energy) (eV)

singlet
triplet standard singlet
triplet dissociation rates (ps
1) deviation

probability

probability

1.2

0.027

0.002

0.104

0.000

1.4 1.6

0.041 0.049

0.001 0.002

0.148 0.182

0.000 0.002

1.8

0.060

0.002

0.222

0.006

2

0.078

0.003

0.294

0.006

2.2

0.083

0.002

0.322

0.010

2.4

0.121

0.007

0.432

0.006

statistical errors. All first-order rate coefficients corresponding to the assigned total energy levels are shown in Table 3. The statistical errors for the calculated rate coefficients vary from 2% to 6% when the total energy is in the range of [1.2 eV, 2.4 eV]. The highest percent error is observed when we investigate the total energy of 2.4 eV (as high as 7%). However, we believe that the percent error for rate constants at low energy levels should be high, as the occurrence of chemical reactions at this stage is rare because of insufficient activation energy, which causes high statistical errors. From the trajectories, it is also possible to monitor the dissociation process (as the multiplicity changes to quintet). However, we observe that the product of dissociation is very limited. As also listed in Table 3, the conversion to triplet dominates largely when it is compared to quintet conversion (dissociation evidence). Therefore, it can be concluded from our study that direct decomposition of ozone in low-excitation energy is very unfavorable, and this statement provides a valid evidence to strongly supports the suggestion made by Xie et al.15b that we have previously discussed in this paper.

V. SUMMARY AND CONCLUSIONS In this study, the singlet
triplet transition and molecular dissociation are investigated using the quasi-classical MD method on a neural network PES constructed by our group. When the PES is developed, a surface hopping consideration is made, and we simplify the problem by assuming that the electronic multiplicity automatically switches to the more energetically favored spin when the molecular configuration reaches the “sensitive region” in the three-dimensional hyperspace (in which the spin hopping probability is assuming to be nearly unity), and the switching process is extremely rapid. This assumption was also adopted in another previous study reported by Agrawal et al.13 to treat the SiO2 systems. The use of this assumption in our study might lead to the computed rate upper limits to reach the true, but unknown reaction rate coefficients. The PES is constructed by fitting a set of 5609 data points, which is primarily collected by performing a grid scan, and the MP4(SDQ) level of theory is chosen to perform ab initio energy calculations. For each generated data point, three different calculations are performed for the three considered spin states (singlet, triplet, and quintet), and the lowest energy among the three values is then selected and stored in the database. Although there are several sampling methods that have been developed and guarantee the effectiveness in data sampling in

multidimensional configuration hyperspace,19,21,30 we realize that such techniques are not necessarily required in this case of study, and a simple grid search would satisfactorily provide the required statistical accuracy to the PES. As previously mentioned, 5906 data points are generated with a symmetry consideration of the two equivalent bond distances. A 60-neuron neural network fit is then performed with 80% percent of the data points (the remaining 20% serve as the testing and validation sets). The mean-absolute error and root-mean-squared errors of this fit are 0.0446 eV (1.03 kcal/mol) and 0.0756 eV (1.74 kcal/mol), respectively. This fitting result reveals a good accuracy in terms of energy to perform classical MD simulations. In this work, we employ pattern recognition in the neural network feature to identify the multiplicity states for the first time, and, experimentally, the technique is shown to be powerful and accurate in reaction time determination. The pattern recognition neural network shares a similar construction with the numerical fitting application, with the threshold transfer function in the last layer being the only distinction between the two. This application is successfully utilized in our study to identify the electronic state of the system, thereby making the investigation scheme more accurate in terms of molecular reaction time, and yields reaction rate coefficients with better-accuracy. Several molecular tests have been executed with zero-point energy and low-excitation energy (with the highest total energy being 0.687 eV), and we observe no chemical reactions occurring in 5 ps. The primary investigation is executed when we investigate several internal energy levels from 1.2 to 2.4 eV, and the reaction rate constants are calculated in the 5 ps time range. The rate coefficients are reported to be in the range of 0.027 to 0.121 ps
1 with very good statistical errors. The highest error is associated with the total energy level of 2.4 eV. The dissociation probability is also monitored in this process, and a very limited number is recorded (less than 1%) below the internal energy of 2.4 eV. As a result, we believe that the singlet
triplet transition is an intermediate step for molecular dissociation when the excitation energy is low (the excitation photon is in the visible and near-infrared regions), which has been previously suggested by Xie et al.15a in a theoretical study that the direct dissociation from singlet ground state is very unfavorable. The direct dissociation of ozone from its singlet state is suggested to occur at high excitation energy and requires a much longer period rather than the period of 5 ps we use in this study.

’ ASSOCIATED CONTENT

bS1

Supporting Information. The weight and bias matrices (w1, b , w2, b2) with more significant digits. The minimum and maximum scaling values (as used in eqs 4 and 5) are stored in the file minmax.dat
for the first two bond distances and degrees for the last in units of Å variable representing the O
O
O bending angle. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*Electronic mail: [email protected]; phone: 84 838350831.

’ ACKNOWLEDGMENT The authors thank the Faculty of Materials Science and the Department of Theoretical Physics, Faculty of Physics, College 10868

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A of Science for their computing facility supports. This work is supported by a research grant funded by the Vietnam National University in Ho Chi Minh City, Vietnam.

’ REFERENCES (1) (a) Flaud, J. M.; Bacis, R. The ozone molecule: Infrared and microwave spectroscopy. Spectrochim. Acta, Part A: Mol. Biomol. Spectrosc. 1998, 54 (1), 3–16. (b) Bacis, R.; Bouvier, A. J.; Flaud, J. M. The ozone molecule: electronic spectroscopy. Spectrochim. Acta, Part A: Mol. Biomol. Spectrosc. 1998, 54 (1), 17–34. (2) Tanaka, T.; Morino, Y. Coriolis interaction and anharmonic potential function of ozone from the microwave spectra in the excited vibrational states. J. Mol. Spectrosc. 1970, 33 (3), 538–551. (3) Braunstein, M.; Pack, R. T. Simple theory of diffuse structure in continuous ultraviolet spectra of polyatomic molecules. III. Application to the Wulf
Chappuis band system of ozone. J. Chem. Phys. 1992, 96 (9), 6378–6388. (4) Braunstein, M.; Martin, R. L.; Hay, P. J. Investigation of the role of triplet states in the Wulf bands of ozone. J. Chem. Phys. 1995, 102 (9), 3662–3666. (5) Banichevich, A.; Peyerimhoff, S. D.; Grein, F. Potential energy surfaces of ozone in its ground state and in the lowest-lying eight excited states. Chem. Phys. 1993, 178 (1
3), 155–188. (6) Tsuneda, T.; Nakano, H.; Hirao, K. Study of low-lying electronic states of ozone by multireference Møller
Plesset perturbation method. J. Chem. Phys. 1995, 103 (15), 6520–6528. (7) Borowski, P.; F€ulscher, M.; Malmqvist, P.-Å.; Roos, B. O. A theoretical study of the low-lying excited states of ozone. Chem. Phys. Lett. 1995, 237 (3
4), 195–203. (8) (a) Bouvier, A. J.; Veyret, V.; Russier, I.; Inard, D.; Churassy, S.; Bacis, R.; Brion, J.; Malicet, J.; Judge, R. H. A comparative rotational analysis of the 000 bands of the 3A2rX1A1 Wulf transition for the isotopomers 16O3 and 18O3 of ozone by high resolution Fourier transform spectrometry. Spectrochim. Acta, Part A: Mol. Biomol. Spectrosc. 1999, 55 (14), 2811–2821. (b) Bouvier, A. J.; Wannous, G.; Churassy, S.; Bacis, R.; Brion, J.; Malicet, J.; Judge, R. H. Spectroscopy and predissociation of the 3A2 electronic state of ozone 16O3 and 18O3 by high resolution Fourier transform spectrometry. Spectrochim. Acta, Part A: Mol. Biomol. Spectrosc. 2001, 57 (3), 561–579. (c) Bouvier, A. J.; Inard, D.; Veyret, V.; Bussery, B.; Bacis, R.; Churassy, S.; Brion, J.; Malicet, J.; Judge, R. H. Contribution to the analysis of the 3 A2rX1A1 “Wulf” transition of ozone by high-resolution Fourier transform spectrometry. J. Mol. Spectrosc. 1998, 190 (2), 189–197. (9) Minaev, B.; Ågren, H. The interpretation of the Wulf absorption band of ozone. Chem. Phys. Lett. 1994, 217 (5
6), 531–538. (10) Minaev, B.; Khomenko, E. Study of singlet
triplet transitions in the ozone molecule using the multiconfigurational self-consistent field theory. High Energy Chem. 2006, 40 (4), 230–233. (11) Grebenshchikov, S. Y.; Qu, Z. W.; Zhu, H.; Schinke, R. New theoretical investigations of the photodissociation of ozone in the Hartley, Huggins, Chappuis, and Wulf bands. Phys. Chem. Chem. Phys. 2007, 9 (17), 2044–2064. (12) (a) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990, 93 (2), 1061–1071. (b) Kohen, D.; Stillinger, F. H.; Tully, J. C. Model studies of nonadiabatic dynamics. J. Chem. Phys. 1998, 109 (12), 4713–4725. (13) Agrawal, P. M.; Raff, L. M.; Hagan, M. T.; Komanduri, R. Molecular dynamics investigations of the dissociation of SiO2 on an ab initio potential energy surface obtained using neural network methods. J. Chem. Phys. 2006, 124 (13), 134306. (14) (a) Hay, P. J.; Pack, R. T.; Walker, R. B.; Heller, E. J. Photodissociation of ozone in the Hartley band. Exploratory potential energy surfaces and molecular dynamics. J. Phys. Chem. 1982, 86 (6), 862–865. (b) Xantheas, S. S.; Atchity, G. J.; Elbert, S. T.; Ruedenberg, K. Potential energy surfaces of ozone. I. J. Chem. Phys. 1991, 94 (12), 8054–8069. (c) Atchity, G. J.; Ruedenberg, K. Global potential energy surfaces for the lowest two 1A0 states of ozone. Theor. Chem. Acc.: Theory, Comput., Model. (Theor. Chim. Acta) 1997, 96 (3), 176–194.

ARTICLE

(15) (a) Xie, D.; Guo, H.; Peterson, K. A. Accurate ab initio nearequilibrium potential energy and dipole moment functions of the ground electronic state of ozone. J. Chem. Phys. 2000, 112 (19), 8378–8386. (b) Xie, D.; Guo, H.; Peterson, K. A. Ab initio characterization of low-lying triplet state potential-energy surfaces and vibrational frequencies in the Wulf band of ozone. J. Chem. Phys. 2001, 115 (22), 10404–10408. (16) (a) Guo, Y.; Kawano, A.; Thompson, D. L.; Wagner, A. F.; Minkoff, M. Interpolating moving least-squares methods for fitting potential energy surfaces: Applications to classical dynamics calculations. J. Chem. Phys. 2004, 121 (11), 5091–5097. (b) Maisuradze, G. G.; Kawano, A.; Thompson, D. L.; Wagner, A. F.; Minkoff, M. Interpolating moving least-squares methods for fitting potential energy surfaces: Analysis of an application to a six-dimensional system. J. Chem. Phys. 2004, 121 (21), 10329–10338. (17) (a) A. Collins, M.; P. A. Bettens, R. Potential energy surface for the reactions BeH2+HBeH+H2. Phys. Chem. Chem. Phys. 1999, 1 (6), 939–945. (b) Collins, M. A.; Zhang, D. H. Application of interpolated potential energy surfaces to quantum reactive scattering. J. Chem. Phys. 1999, 111 (22), 9924–9931. (c) Bettens, R. P. A.; Collins, M. A. Learning to interpolate molecular potential energy surfaces with confidence: A Bayesian approach. J. Chem. Phys. 1999, 111 (3), 816–826. (18) Hagan, M. T.; Demuth, H. B.; Beale, M. Neural Network Design; Martin Hagan: Boulder, CO, 1996. (19) Malshe, M.; Raff, L. M.; Rockley, M. G.; Hagan, M.; Agrawal, P. M.; Komanduri, R. Theoretical investigation of the dissociation dynamics of vibrationally excited vinyl bromide on an ab initio potential-energy surface obtained using modified novelty sampling and feedforward neural networks. II. Numerical application of the method. J. Chem. Phys. 2007, 127 (13), 134105. (20) Le, H. M.; Raff, L. M. Cis f trans, trans f cis isomerizations and N
O bond dissociation of nitrous acid (HONO) on an ab initio potential surface obtained by novelty sampling and feed-forward neural network fitting. J. Chem. Phys. 2008, 128 (19), 194310. (21) Le, H. M.; Huynh, S.; Raff, L. M. Molecular dissociation of hydrogen peroxide (HOOH) on a neural network ab initio potential surface with a new configuration sampling method involving gradient fitting. J. Chem. Phys. 2009, 131 (1), 014107. (22) Le, H. M.; Raff, L. M. Molecular dynamics investigation of the bimolecular reaction BeH + H2 f BeH2 + H on an ab initio potentialenergy surface obtained using neural network methods with both potential and gradient accuracy determination. J. Phys. Chem. A 2009, 114 (1), 45–53. (23) Pukrittayakamee, A.; Malshe, M.; Hagan, M.; Raff, L. M.; Narulkar, R.; Bukkapatnum, S.; Komanduri, R. Simultaneous fitting of a potential-energy surface and its corresponding force fields using feedforward neural networks. J. Chem. Phys. 2009, 130 (13), 134101. (24) Raff, L. M.; Thompson, D. L. Theory of Chemical Reaction Dynamics; CRC Press: Boca Raton, FL, 1985; Vol. 3. (25) (a) Møller, C.; Plesset, M. S. Note on an approximation treatment for many-electron systems. Phys. Rev. 1934, 46 (7), 618. (b) Krishnan, R.; Pople, J. A. Approximate fourth-order perturbation theory of the electron correlation energy. Int. J. Quantum Chem. 1978, 14 (1), 91–100. (26) (a) Trucks, G. W.; Watts, J. D.; Salter, E. A.; Bartlett, R. J. Analytical MBPT(4) gradients. Chem. Phys. Lett. 1988, 153 (6), 490–495. (b) Trucks, G. W.; Salter, E. A.; Sosa, C.; Bartlett, R. J. Theory and implementation of the MBPT density matrix. An application to oneelectron properties. Chem. Phys. Lett. 1988, 147 (4), 359–366. (27) (a) McLean, A. D.; Chandler, G. S. Contracted Gaussian basis sets for molecular calculations. I. Second row atoms, Z = 11
18. J. Chem. Phys. 1980, 72 (10), 5639–5648. (b) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions. J. Chem. Phys. 1980, 72 (1), 650–654. (28) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, 10869

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870

The Journal of Physical Chemistry A

ARTICLE

J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Laham, A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision E.01; Gaussian, Inc.: Wallingford, CT, 2003. (29) Raff, L. M.; Malshe, M.; Hagan, M.; Doughan, D. I.; Rockley, M. G.; Komanduri, R. Ab initio potential-energy surfaces for complex, multichannel systems using modified novelty sampling and feedforward neural networks. J. Chem. Phys. 2005, 122 (8), 084104. (30) Agrawal, P. M.; Malshe, M.; Narulkar, R.; Raff, L. M.; Hagan, M.; Bukkapatnum, S.; Komanduri, R.; Self-Starting, A Method for obtaining analytic potential-energy surfaces from ab initio electronic structure calculations. J. Phys. Chem. A 2009, 113 (5), 869–877. (31) MATLAB; MathWorks: Natick, MA, 2011. (32) Hornik, K.; Stinchcombe, M.; White, H. Multilayer feedforward networks are universal approximators. Neural Networks 1989, 2 (5), 359–366. (33) (a) Malshe, M.; Narulkar, R.; Raff, L. M.; Hagan, M.; Bukkapatnam, S.; Agrawal, P. M.; Komanduri, R. Development of generalized potentialenergy surfaces using many-body expansions, neural networks, and moiety energy approximations. J. Chem. Phys. 2009, 130 (18), 184102. (b) Manzhos, S.; Tucker Carrington, J. Using neural networks to represent potential surfaces as sums of products. J. Chem. Phys. 2006, 125 (19), 194105. (c) Agrawal, P. M.; Samadh, A. N. A.; Raff, L. M.; Hagan, M. T.; Bukkapatnam, S. T.; Komanduri, R. Prediction of molecular-dynamics simulation results using feedforward neural networks: Reaction of a C2 dimer with an activated diamond (100) surface. J. Chem. Phys. 2005, 123 (22), 224711. (34) Caruana, R.; Lawrence, S.; Giles, C. L. Overfitting in neural nets: Backpropagation, conjugate gradient, and early stopping. In Proceedings of the Neural Information Processing Systems Conference; MIT Press: Cambridge, MA, 2000; pp 402
408. (35) Raff, L. M. Projection methods for obtaining intramolecular energy transfer rates from classical trajectory results: Application to 1, 2-difluoroethane. J. Chem. Phys. 1988, 89 (9), 5680–5691.

10870

dx.doi.org/10.1021/jp206531s |J. Phys. Chem. A 2011, 115, 10862–10870