J. Phys. Chem. A 2010, 114, 45–53
45
Molecular Dynamics Investigation of the Bimolecular Reaction BeH + H2 f BeH2 + H on an ab Initio Potential-Energy Surface Obtained Using Neural Network Methods with Both Potential and Gradient Accuracy Determination Hung M. Le and Lionel M. Raff* Chemistry Department, Oklahoma State UniVersity, Stillwater, Oklahoma 74078 ReceiVed: August 4, 2009; ReVised Manuscript ReceiVed: September 25, 2009
The classical reaction dynamics of a four-body, bimolecular reaction on a neural network (NN) potentialenergy surface (PES) fitted to a database obtained solely from ab initio MP2/6-311G(d,p) calculations are reported. The present work represents the first reported application of ab initio NN methods to a four-body, bimolecular, gas-phase reaction where bond extensions reach 8.1 Å for the BeH + H2 f BeH2 + H reaction. A modified, iterative novelty sampling method is used to select data points based on classical trajectories computed on temporary NN surfaces. After seven iterations, the sampling process is found to converge after selecting 9604 configurations. Incorporation of symmetry increases this to 19 208 BeH3 configurations. The analytic PES for the system is obtained from the ensemble average of a five-member (6-60-1) NN committee. The mean absolute error (MAE) for the committee is 0.0046 eV (0.44 kJ mol-1). The total energy range of the BeH3 database is 147.0 kJ mol-1. Therefore, this MAE represents a percent energy error of 0.30%. Since it is the gradient of the PES that constitutes the most important quantity in molecular dynamics simulations, the paper also reports mean absolute error for the gradient. This result is 0.026 eV Å-1 (2.51 kJ mol-1 Å-1). Since the gradient magnitudes span a range of 15.32 eV Å-1 over the configuration space tested, this mean absolute gradient error represents a percent error of 0.17%. The mean percent absolute relative gradient error is 4.67%. The classically computed reaction cross sections generally increase with total energy. They vary from 0.007 to 0.030 Å2 when H2 is at ground state, and from 0.05 to 0.10 Å2 when H2 is in the first excited state. Trajectory integration is very fast using the five-member NN PES. The average trajectory integration time is 1.07 s on a CPU with a clock speed of 2.4 GHz. Zero angular momentum collisions are also investigated and compared with previously reported quantum dynamics on the same system. The quantum reaction probabilities exhibit pronounced resonance effects that are absent in the classical calculations. The magnitudes of quantum and classical results are in fair accord with the classical results being about 30-40% higher due to the lack of quantum restrictions on the zero-point vibrational energy. I. Introduction 1
Neural network (NN) fitting has been shown to a powerful and robust method for the execution of many studies associated with chemical reaction dynamics. These include the fitting of potential-energy databases obtained from ab initio electronic structure calculations,2-23 parameter fitting of empirical surfaces,24 statistical error reduction in molecular dynamics (MD) investigations,25 acceleration of genetic algorithm (GA) studies,26 and most recently the prediction of high-level electronic structure ab initio energies from lower-level, Hartree-Fock energies.27 The research reported in this paper focuses on the first of these applications, the fitting of ab initio databases to obtain an analytic expression for the potential-energy surface (PES) of the system. If the PES is to be used to compute the vibrationalenergy levels of the system, Manzhos and Carrington14 have shown that it is possible to obtain accuracies on the order of 1-2 cm-1 for three- and four-atom molecules by using a nested NN approach in which one NN is employed to obtain an approximate fit to the database and the second NN is then used to fit the differences between the target ab initio energies and the predictions of the first NN. When chemical reactions are involved for three- and fouratom molecules, surface fitting becomes more difficult and the accuracy decreases accordingly. For example, the unimolecular
dissociation of H2O2 by O-O bond rupture was recently examined by Le et al.20 using a PES fitted to an ab initio database obtained using a new gradient sampling method. With only a single bond rupture reaction involved, the mean absolute fitting error (MAE) of the NN to the database was 0.58 kJ mol-1. Nitrous acid (HONO), which was studied by Le and Raff4 using NN methods, has two reaction channels: one dissociation channel (the breaking of N-O bond) and cis-trans isomerization. With two open reaction channels, fitting is more difficult. As a result, the MAE was 1.07 kJ mol-1. Later, a NN PES for HONO was developed by Agrawal et al.5 using a different technique for data sampling. The results of this investigation were consistent with those reported by Le and Raff.4 Pukrittayakamee et al.6,7 have recently made NN methods more powerful and accurate by requiring the training process to simultaneously fit both the potential energy and its gradients with respect to the input variables. Application of this method to inelastic scattering, exchange, and abstraction reactions occurring in the H + HBr system has shown excellent fitting accuracy with a median root-mean-square (rms) fitting error of 0.015 kJ mol-1 or 1.2 cm-1. The derivatives were also fitted with excellent accuracy with a median rms error of 0.162 kJ mol-1 Å-1. This NN method is denoted as the combined function derivative approximation (CFDA).
10.1021/jp907507z 2010 American Chemical Society Published on Web 10/23/2009
46
J. Phys. Chem. A, Vol. 114, No. 1, 2010
Figure 1. Definition of atomic distance identities of BeH3 system. Atom number 1 is Be.
Le and Raff The dynamics of the reaction were investigated using quasiclassical MD. In an extended work,28b the PES was further improved by employing a database containing 1300 points using the same technique. The electronic structure calculations in these studies were executed using second-order Møller-Plesset perturbation theory (MP2)31-35 with a 6-311G(d,p) basis set. It was reported that this method converges well for this BeH3 system;28a,b therefore, in this work, we also employ the MP2/ 6-311G(d,p) level of theory incorporated in the Gaussian suite of programs36 for electronic structure calculations. Collins and Zhang28b reported quantum37-39 reactive scattering calculations of the reaction probabilities for H-atom abstraction reactions of BeH with H2 as a function of translational energy for collisions with zero angular momentum (J ) 0). II. Neural Network Fitting
The most complex and difficult system yet studied by any fitting technique is the unimolecular dissociation of vinyl bromide (VB) into six, simultaneously open reaction channels that include four two-center bond-rupture processes and two three-center dissociations leading to HBr and H2. This 12dimensional system was originally investigated by Doughan et al.2 and by Malshe et al.3 Because of its high dimensionality and the large configuration-space volume of the system, NN fitting is much more difficult. Malshe et al.3 originally obtained an MAE with a (15-140-1) NN of 6.28 kJ mol-1. This has recently been improved27 to 6.16 kJ mol-1. Manzhos and Carrington17 used a high-dimensional model representation (HDMR) with NNs as six-dimensional component functions to obtain a fit to the VB database whose error was very similar to that reported by Malshe et al.3 even though the problem was treated with six-dimensional NNs rather than the 15-dimensional NN used by Malshe et al.3 When Malshe et al.23 used a manybody expansion formalism truncated after the four-body terms with NNs as the many-body expansion functions, the MAE was about 15% higher than that achieved by Manzhos and Carrington.17 In this paper, we report the first application of NN fitting of an ab initio database for a four-atom bimolecular reaction. The test system is the bimolecular reaction of BeH with H2 to yield BeH2 and H atoms: BeH + H2 f BeH2 + H (Figure 1). The configuration space of this system is expected to be much larger than the PESs that describe unimolecular dissociations. Therefore, a NN with more hidden neurons is required to fit the BeH3 system. The PES property of greatest importance in an MD study is the surface gradient, not the potential itself. This is the case because the MD trajectories do not depend upon the potential. They are determined solely by the intermolecular forces which are given by the surface gradients at all integration points. In spite of this fact, very few of the previously reported papers dealing with PES fitting by any method have reported data on the accuracy of the gradients obtained from their fitted PESs. Some exceptions include the H2Br system studied by Pukrittayakamee et al.,6,7 the BeH3 and CH3+ systems investigated by Bettens and Collins,28c and the OH + H2 f H2O + H reaction examined by Thompson and Collins.29 The results of an investigation of the gradient accuracy are reported in this paper. Theoretical studies of the adiabatic chemical reaction BeH + H2 f BeH2 + H, where the multiplicities of the reactant and product systems are both doublet, have been reported by Collins and co-workers.28 In their first investigation,28a an ab initio PES was obtained using a modified Shepard interpolation (MSI) method30 with a database comprising 438 BeH3 configurations.
Ab initio PESs for several systems have been obtained using NN methods to fit the databases.2-23 This technique has been proven to be powerful and robust since the fit is done only once, and retrieving stored data points is not required during trajectory calculations. We employ the NN toolbox in Matlab to perform the fit. Full features and properties of this toolbox are discussed by Hagan, Demuth, and Beale.1 A two-layer, feed-forward neural network is used to fit the ab initio database. For the BeH3 system, we choose six atomic distances as shown in Figure 1 as inputs to the NN. The output layer yields the ab initio potential energy. Before fitting, the distances, ri, and the output potential energy V are scaled using
rscaledi )
2(ri - ri_min) -1 ri_max - ri_min
(1)
Vscaled )
2(V - Vmin) -1 Vmax - Vmin
(2)
where ri_max, ri_min and Vmax, Vmin are the maximum and minimum values of the atomic distances and potential energy, respectively. Equations 1 and 2 scale all input and output data to the range [-1, 1]. The reduced forms of eqs 1 and 2 ensure that the input and output variables are unitwise consistent. With m neurons in the first layer, the output of this layer is 6
ni )
∑ w1j,irscaled j + b1i
for i ) 1, ..., m
(3)
j)1
where w1 is the “weight” matrix of the first layer of size m × 6 and b1 is an m × 1 bias vector for the first layer. The ni are then converted to the input vector for the second layer by using a transfer function f1, for which we employ a hyperbolic tangent function. Finally, the output of the two-layer network is given by m
Vout )
∑ w2i f 1(ni) + b
(4)
i)1
where w2 is a 1 × m weight matrix and b is a 1 × 1 bias vector. The transfer function for the output layer is taken to be a linear function. Hornik et al.40 have shown feed-forward NNs that use
MD Investigation of BeH + H2 f BeH2 + H
J. Phys. Chem. A, Vol. 114, No. 1, 2010 47
TABLE 1: Average Absolute Error of the Networks network
average absolute error (eV)
1 2 3 4 5 NN committee
0.0046 0.0051 0.0046 0.0049 0.0047 0.0046
a sigmoid transfer function in the first layer and a linear function in the second layer are universal approximators for analytic functions. The number of neurons in the hidden layer is selected using an empirical investigation. In this procedure, an initial value of m is selected to train the network. This initial value is gradually increased until a point is reached at which further increases in m produce only marginal improvement in the rms error for an independent testing set. For the complex, six-body vinyl bromide (CH2CHBr) molecule,2,3 the number of neurons used by Malshe et al.3 for the final fitting was 140. Less-complicated molecules such as HONO4,5 and SiO222 require approximately 40 neurons to fit the data. In this study, we employ a NN with 60 neurons in the hidden layer. Once an appropriate number for hidden neurons is chosen, the fitting process is done only once, and all parameters are saved for MD investigation. The final database is divided into three sets: training, testing, and validation. Empirically, the training set contains about 80% of the data, while the other two sets contain about 10%. The iterative Levenberg-Marquardt algorithm is used to train the NN so as to minimize the mean square error of the training set. In the present work, the maximum number of iterations permitted was 1000. During the training process, an important concern is overfitting.41 To prevent overfitting, we employ “early-stopping”,1 to terminate the training process. In this method, the rms fitting error for the validation set is monitored as the training proceeds. If overfitting begins, it will be seen by a decreasing fitting error on the training set but an increasing error on the validation set. When such an increase in the validation set error occurs in six consecutive training iterations, training is terminated, and the weights and biases in existence when the validation set error first began to increase are taken to be those for the trained NN. The actual potential surface of BeH3 is a committee of five, (6-60-1) feed-forward NNs. After fitting each different NN individually by selecting different sets of training data, the NN committee result is calculated by taking the average of the energy and gradients given by the five networks. Recently, NN committees have been employed in HONO4 and HOOH20 studies. The reported fitting error of the NN committee is usually significantly less than the errors of the individual networks because of the cancellation of nearly random errors. In this study, the MAEs for each network vary from 0.0046 to 0.0051 eV. The average absolute error of the NN committee average is 0.0046 eV, which is slightly improved in comparison with the errors of individual networks. The fitting errors for all members of the committee as well as that for the committee are given in Table 1. The use of NN committee in this study does not reduce the error as much as has been observed for other systems.4,20 The absolute error distribution for all data is shown in Figure 2.
Figure 2. Absolute error distribution of all data using the five-member NN committee. The number of small errors ( 0.13 (c) avgnew > 0.3
or or
(8)
In the second iteration, 5181 new configurations that satisfy the selection criteria given in condition (8) are identified using sampling trajectories. Subsequently, ab initio calculations are executed for these configurations. A new configuration is added to the database if the difference between its NN energy (predicted by the temporary surface) and the ab initio energy is greater than the MAE of the temporary fit for all previous data points. During the second iteration, these methods result in the addition of 3178 new data points to the overall database. For the next iteration, a temporary NN is now fitted to the current database that now comprises 4478 BeH3 configurations. The process is repeated iteratively until the number of new configurations added to the database becomes very small. During the seventh iteration for the present study, 5027 new configurations were identified by the trajectory sampling. After making comparisons between the NN-predicted energies and ab initio energies, 484 configurations were added to the database. At this point, the configuration sampling process was terminated. The difficulty associated with identifying new configurations to add to the database suggests that we have reached convergence. The detailed numbers of selected data points during this sampling process are reported in Table 2. After seven iterations, the BeH3 database contains 9642 configurations. This number is reduced to 9604 by removing 38 points that are associated with energies too high to be significant in the reaction dynamics of the BeH + H2 system. By switching r4 and r5, r2 and r6 for symmetry considerations, the final database contains 19 208 configurations. As previously mentioned, the BeH3 database is randomly partitioned between training, validation, and testing sets five times to generate a five-member NN committee. The average energy and gradients predicted by this committee are taken to be those for the BeH3 system. The MAE results for the potential for each member of the committee are given in Table 1. The committee MAE is 0.0046 eV or 0.44 kJ mol-1. Since the potential energies of the points in the database span an energy range of 147.0 kJ mol-1, this MAE represents an average percent energy error of 0.30%. The gradients calculated by the NN committee are very important in MD investigations. Although they are not directly fitted by the NN, we expect them to be well predicted. A set of 965 points is randomly chosen as a gradient testing set. The energy gradients with respect to 12 Cartesian coordinates are computed and compared to the MP2 gradients. The average
Figure 4. Distribution plots of absolute gradient errors.
absolute error is 0.026 eV Å-1. Since the magnitudes of the gradients at the testing points span a range of 15.32 eV Å-1, this mean absolute error represents a percent gradient error of 0.17%. To facilitate comparison with the results reported by Bettens and Collins,28c we have also computed the percent absolute relative gradient error, E∇V, and its mean. As defined by Bettens and Collins,28c E∇V is given by
E∇V ) 2
|
∇V - ∇VNN
|
× 100
(9)
∇V + ∇VNN
| | | | where ∇V and ∇VNN are the MP2 and NN predicted gradients, respectively. The mean percent absolute relative gradient error, 〈E∇V〉, is obtained by averaging E∇V over the 965-point gradient testing set. This result is 〈E∇V〉 ) 4.67%. The value of 〈E∇V〉 reported by Bettens and Collins28c depends upon the parameters chosen for the MSI method. Their best result was 5.58%. Figure 4 shows the distribution of absolute gradient errors. The fact that a preponderance of the test points exhibit very small absolute errors shows that not only is the potential accurately fitted, but the gradients are also predicted with very good accuracy. IV. Classical Dynamics Investigation of BeH + H2 f BeH2 + H 1. Initializing the BeH + H2 Reaction. Prior to the integration of BeH + H2 trajectories, the initial states of the reactants are sampled from appropriate distributions. In order to do this effectively, we have constructed two simple ab initio NN PESs for BeH and H2. BeH has a doublet multiplicity in its ground state. Ab initio calculations were executed at the MP2 level of theory using the 6-311G(d,p) basis set with various atomic distances between Be and H from 0.9 to 3.15 Å. A feed-forward NN fit of BeH potential was obtained using four neurons in the hidden layer. This gives a BeH potential having the form
50
J. Phys. Chem. A, Vol. 114, No. 1, 2010
{∑
[ (
4
VBeH(r) )
BeH wBeH tanh w1i 2 2
i)1
BeH r - rmin BeH BeH rmax - rmin
Le and Raff
) ] }
TABLE 4: Fitting Parameters of H2 Potential Surfacea
- 1 + bBeH + 1 bBeH +1 2
BeH Vmax 2
(10)
BeH BeH BeH where rmax ) 3.15 Å and rmin ) 0.9 Å, and Vmax ) 2.35 eV. BeH BeH BeH BeH The values of w1 , w2 , b1 , and b2 are reported in Table 3. We obtain an excellent absolute average error of 0.0008 eV. The plot of the fit is shown in Figure 5. The BeH molecule is initially aligned along the Z-axis and placed in its equilibrium configuration with r ) 1.34 Å. A ground state vibrational energy of 0.127 eV (12.3 kJ mol-1) is introduced as kinetic energy, and the trajectory is integrated for a random period of time using the fourth-order Runge-Kutta integration method45 and eq 10 in order to average over the BeH vibrational phase angle. Subsequently, the molecular Cartesian coordinates and momenta are rotated about the y-axis by an angle chosen randomly in the range [0, π] and then rotated about the z-axis by a randomly chosen angle over the range [0, 2π]. Molecular hydrogen, H2, is singlet in its ground electronic state. We also execute ab initio calculations for H2 using the same level of theory at various atomic distances from 0.45 to 1.54 Å. The potential of H2 has the form
VHH(r) )
{
4
∑w i)1
HH 2
[ (
HH tanh w1i 2
HH r - rmin HH rmax
HH - rmin
) ] }
bHH 1
wHH 2
bHH 2
–1.04 –1.65 –2.63 –2.40
0.920 0.146 –0.845 –3.631
–1.181 –0.508 –0.130 15.068
15.33
All parameters are unitless.
- 1 + bHH + 1 bHH 2 + 1
HH Vmax 2
(11)
HH HH HH where rmax ) 1.54 Å and rmin ) 0.45 Å, and Vmax ) 3.45 eV. HH HH HH HH The values of w1 , w2 , b1 , and b2 are shown in Table 4.
TABLE 3: Fitting Parameters of BeH Potential Surfacea
a
a
wHH 1
wBeH 1
bBeH 1
wBeH 2
bBeH 2
6.25 2.62 4.14 3.51
–5.64 0.30 2.09 4.42
–0.006 0.817 0.186 –8.584
8.18
All parameters are unitless.
Figure 5. Potential energy curve of BeH. The solid curve represents the fit by eq 10, and the circles (O) represent the computed ab initio potential energies of BeH.
Figure 6. Potential energy curve of H2. The solid curve represents the fit by eq 11, and the circles (O) represent the computed ab initio points of H2 potential.
The MAE is 0.0004 eV. Figure 6 compares the fitted potential curve with the ab initio energies. The sampling process of H2 is done similarly to the sampling process of BeH with a zero-point energy of 0.261 eV (25.2 kJ mol-1) being used for H2 trajectories. After averaging over vibrational phase angle and orientation, the two molecular centers of mass are separated by a distance of 5.82 Å, and a random impact parameter, b, is chosen using b ) ξ1/2bmax, where bmax is the maximum value of 0.265 Å. The desired amount of relative translational energy is then inserted into the system with no center-of-mass motion. 2. Classical Dynamics Investigation. The results reported by Collins and Zhang28b indicate that the reaction probability derived from quantum dynamics calculations37-39 has a maximum value of 0.11 for translational energies in the range 40-77 kJ mol-1 when the total angular momentum of the system is zero. Classically, this corresponds to an initial state with zero rotational energy for both BeH and H2 and b ) 0. The potential barrier has been investigated by Collins and Bettens.28a At the MP2/6-311G(d,p) level of theory, the potential barrier is reported to be around 51 kJ mol-1 (0.529 eV). In our first investigation, we use classical dynamics trajectories at zero impact parameter to obtain the classical analogue of the J ) 0 quantum result. The translational energy is varied from 0.415 to 0.829 eV. The Hamiltonian motion equations are integrated using the fourth-order Runge-Kutta method which yields energy conservation to five significant digits. A trajectory is terminated by one of two criteria: (1) The distance between BeH and H-H centers of mass continues to increases and reaches 5.83 Å. In this case, no chemical reactions are found; the two molecules scatter inelastically.
MD Investigation of BeH + H2 f BeH2 + H
J. Phys. Chem. A, Vol. 114, No. 1, 2010 51
Figure 8. Plot of reaction cross sections vs various relative translational energies when H2 is in its ground or first excited vibrational state. The reaction is investigated when BeH and H2 are both in their ground rotational states (j ) 0) and BeH is in its vibrational ground state. The maximum impact parameter in this investigation is 0.265 Å.
Figure 7. Plot of reaction probability vs various relative translational energy when both gases are at their ground vibrational states (V ) 0). The reaction is investigated when BeH and H2 both have no rotations (j ) 0), and the total angular momentum is 0 (J ) 0). The quantum dynamics results are reported by Collins and Zhang.28b The error bars on the classical calculations represent one sigma limit of statistical uncertainty using 5000 trajectories for each relative translational energy.
(2) The distance between Be and one H from the H2 molecule is less than 1.5 Å, and the distance between two hydrogen atoms in H2 is greater than 2.5 Å. In this case, the formation of BeH2 is found. For each translational energy level, 5000 trajectories are integrated to obtain average reaction probabilities. The reaction probabilities at relative translational energies from 0.41 to 0.83 eV at 0.01 eV intervals have been computed for the J ) 0 case. The results are shown in Figure 7. The error bars represent one sigma limit of statistical uncertainty. Figure 7 also shows the quantum mechanical results reported by Collins and Zhang.28b The agreement between classical and quantum results is fair. The large resonance effects seen in the quantum results are, of course, absent in the classical MD simulations. In addition, the availability of the zero-point vibrational energy in classical calculations to assist in traversing the reaction barrier is seen to artificially enhance the reaction probabilities by about 30-40% over those obtained in the quantum studies.28b The fact that this enhancement is due to the lack of quantum mechanical restrictions on the zero-point energy is shown clearly by the classical MD calculations carried out in the absence of zero-point energy. These results are shown in Figure 7 by the dashed-dotted curve. When zero-point energy is removed, the classical reaction probabilities fall significantly below the quantum results. In quantum calculations, some of the reactant zero-point energy is available to the system to assist in traversing the reaction barrier and thereby enhancing the reaction prob-
ability. This fraction is roughly equal to the difference in vibrational zero-point energy present in the reactants and that present in the transition-state complex. As a result, the quantum probabilities are intermediate between classical results with no zero-point energy and those with full zero-point energy. For the same reasons, the reaction threshold for classical calculations with full zero-point energy falls below that for the quantum studies while the classical threshold lies above the quantum threshold when all zero-point energy is omitted. Figure 8 and Table 5 report the computed reaction cross sections for hydrogen initially in its vibrational ground state, V ) 0, and its first excited vibrational state, V ) 1, with full zeropoint energy and a range of total angular momentum states obtained by random selection of the impact parameter from zero to a maximum value of 0.265 Å. As expected with a late transition state located in the product valley, dynamics effects play an important role in the reaction with the result that reactant H2 vibrational energy is much more effective in promoting reaction than is relative translational energy.46,47 Quantitatively, it is found that when H2 is in its ground state, the cross sections vary from 0.007 to 0.030 Å2. On the other hand, when H2 is promoted to the first excited state, the cross sections vary from 0.05 to 0.10 Å2. V. Summary and Conclusions The classical reaction dynamics of a four-body, bimolecular reaction on a NN PES fitted to a database obtained solely from ab initio MP2/6-311G(d,p) calculations are reported. Previous applications of NN methods to four-body and larger systems have been to unimolecular reactions where the atomic distances are more restricted. The present work represents the first reported application of ab initio NN methods to a four-body, bimolecular, gas-phase reaction where bond extensions reach 8.1 Å for the BeH + H2 f BeH2 + H reaction. A modified, iterative novelty sampling19 method is used in this study to select data points based on classical trajectories computed on temporary NN surfaces. After seven iterations, the sampling process is found to converge after selecting 9642 configurations. Thirty-eight of these configurations with high potential energies are eliminated as they are unimportant in the MD calculations. This database is then duplicated to incorporate symmetry considerations to yield a final database of 19 208 BeH3 configurations. The analytic PES for the system is obtained from the ensemble average of a five-member NN committee in which
52
J. Phys. Chem. A, Vol. 114, No. 1, 2010
Le and Raff
TABLE 5: Averaged Reaction Probabilities and Cross Sections for BeH + H2f BeH2 + H as a Function of Relative Translational Energy and H2 Vibrational State for Initially Nonrotating BeH and H2 Reactantsa reaction probability cross section translational energy (eV) V ) 0, j ) 0 V ) 1, j ) 0 V ) 0, j ) 0 V ) 1, j ) 0 0.41 0.42 0.44 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.60 0.61 0.62 0.63 0.64 0.65 0.66 0.67 0.68 0.69 0.70 0.72 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83
0.032 0.040 0.047 0.051 0.050 0.053 0.065 0.072 0.080 0.077 0.084 0.079 0.086 0.087 0.086 0.097 0.088 0.095 0.093 0.095 0.098 0.104 0.103 0.102 0.103 0.102 0.104 0.106 0.112 0.109 0.120 0.115 0.120 0.115 0.121 0.116 0.125 0.120 0.129 0.131 0.136
0.007 0.009 0.010 0.011 0.011 0.012 0.014 0.016 0.018 0.017 0.018 0.017 0.019 0.019 0.019 0.021 0.019 0.021 0.020 0.021 0.022 0.023 0.023 0.022 0.023 0.022 0.023 0.023 0.025 0.024 0.026 0.025 0.026 0.025 0.027 0.026 0.027 0.026 0.028 0.029 0.030
0.236 0.229 0.242 0.229 0.243 0.242 0.239 0.239 0.233 0.256 0.259 0.264 0.261 0.276 0.275 0.290 0.293 0.298 0.312 0.315 0.333 0.334 0.336 0.348 0.343 0.361 0.368 0.366 0.368 0.387 0.381 0.389 0.399 0.404 0.414 0.421 0.441 0.432 0.436 0.456 0.470
0.052 0.050 0.053 0.050 0.053 0.053 0.053 0.053 0.051 0.056 0.057 0.058 0.057 0.061 0.060 0.064 0.064 0.066 0.069 0.069 0.073 0.074 0.074 0.077 0.075 0.079 0.081 0.081 0.081 0.085 0.084 0.085 0.088 0.089 0.091 0.093 0.097 0.095 0.096 0.100 0.103
a The (V, j) quantum numbers given in the column headings refer to the initial quantum states for H2 in the collisions. BeH is initially in its ground vibrational state in all calculations.
each (6-60-1) NN is fitted to the ab initio database. The NNs are obtained by randomly partitioning the database into training (80%), validation (10%), and testing (10%). The MAEs for each network in the committee range from 0.0046 to 0.0051 eV. The committee average is 0.0046 eV (0.44 kJ mol-1). This testing set MAE is very close to the testing error reported by Collins and Zhang,28b which is 0.41 kJ mol-1. The total energy range spanned by the present database is 147.0 kJ mol-1, so our MAE represents a percent error of 0.30%. The energy range employed by Collins and Bettens28c in their study of the BeH3 system was 105.6 kJ mol-1. Their reported mean percent energy error was, therefore, 0.39%. Since it is the gradient of the PES that constitutes the most important quantity in MD simulations, reporting the average fitting error for the gradient should be standard practice in all surface fitting papers. Comparison of the magnitudes of the 12 Cartesian derivatives obtained from the NN committee with corresponding MP2 results at 965 randomly selected testing ponts yields a mean absolute gradient error of 0.026 eV Å-1 (2.51 kJ mol-1 Å-1). This is 0.17% of the total range of derivative magnitudes present at the testing points. The mean
percent absolute relative gradient error evaluated at the same testing points is 4.67%. The MSI result28c is 5.58%. The investigation of the reaction dynamics is executed using classical trajectories. With a maximum impact parameter of 0.265 Å, reaction cross sections at relative translational energies in the range 0.415-0.829 eV are computed by using 5000 trajectories at each energy to average over vibrational phase angle, impact parameter, and molecular orientation. The reaction cross sections generally increase with total energy. They vary from 0.007 to 0.030 Å2 when H2 is at ground state, and from 0.05 to 0.10 Å2 when H2 is in the first excited state. Trajectory integration is very fast using the five-member NN PES. The average trajectory integration time is 1.07 s on a CPU with a clock speed of 2.4 GHz. If a single NN were employed, this time would decrease to about 0.21 s. Zero angular momentum collisions are also investigated and compared with previously reported quantum dynamics on the same system where an MSI method was employed to fit the MP2/6-311G(d,p) database.28b The quantum reaction probabilities exhibit pronounced resonance effects that are absent in the classical calculations. The magnitudes of quantum and classical results are in fair accord with the classical results being about 30-40% higher due to the lack of quantum restrictions on the zero-point vibrational energy. Acknowledgment. We sincerely thank Professor M. A. Collins from the National University of Australia for graciously providing his database of 1300 BeH3 configurations and associated energies. We also thank Professor Collins for providing other helpful advice and cogent comments related to this research. References and Notes (1) Hagan, M. T.; Demuth, H. B.; Beale, M. Neural Network Design; Colorado University Bookstore: Boulder, 1996. (2) Doughan, D. I.; Raff, L. M.; Rockley, M. G.; Hagan, M. T.; Agrawal, P. M.; Komanduri, R. J. Chem. Phys. 2006, 124, 054321. Doughan, D. I.; Raff, L. M.; Rockley, M. G.; Hagan, M. T.; Agrawal, P. M.; Komanduri, R. J. Chem. Phys. 2006, 125, 079901. (3) Malshe, M.; Raff, L. M.; Rockley, M. G.; Agrawal, P. M.; Komanduri, R. J. Chem. Phys. 2007, 127, 134105. (4) Le, H. M.; Raff, L. M. J. Chem. Phys. 2008, 128, 194310. (5) Agrawal, P. M.; Malshe, M.; Narulkar, R.; Raff, L. M.; Hagan, M. T.; Bukkapatnum, S.; Komanduri, R. J. Phys. Chem. A 2009, 113, 869. (6) Pukrittayakamee, A.; Malshe, M.; Hagan, M. T.; Raff, L. M.; Narulkar, R.; Bukkapatnum, S.; Komanduri, R. J. Chem. Phys. 2009, 130, 134101. (7) Pukrittayakamee, A.; Hagan, M. T.; Raff, L. M.; Bukkapatnam, S.; Komanduri, R. Intell. Eng. Syst. Artif. Neural Networks 2007, 17 (ANNIE2007). (8) Blank, T. B.; Brown, S. D. Anal. Chem. Acta 1993, 277, 273. (9) Blank, T. B.; Brown, S. D.; Calhoun, S. W.; Doren, D. J. J. Chem. Phys. 1995, 103, 4129. (10) Hobday, S.; Smith, R.; BelBruno, J. Nucl. Instrum. Methods Phys. Res., Sect. B 1999, 153, 247. (11) Tafeit, E.; Estelberger, W.; Horejsi, R.; Moeller, R.; Oettl, K.; Vrecko, K.; Reibnegger, G. J. Mol. Graphics 1996, 14, 12. (12) Gassner, H.; Probst, M.; Lauenstein, A.; Hermansson, K. J. Phys. Chem. A 1998, 102, 4596. (13) Brown, D. R.; Gibbs, M. N.; Clary, D. C. J. Chem. Phys. 1996, 105, 7597. (14) Manzhos, S.; Wang, X.; Dawes, R.; Carrington, T., Jr. J. Phys. Chem. A 2006, 110, 5295. (15) Manzhos, S.; Carrington, T. J., Jr. J. Chem. Phys. 2006, 125, 194105. (16) Manzhos, S.; Carrington, T. J., Jr. J. Chem. Phys. 2006, 125, 084109. (17) Manzhos, S.; Carrington, T. J., Jr. J. Chem. Phys. 2008, 129, 224104. (18) Lorenz, S.; Gross, A.; Scheffier, M. Chem. Phys. Lett. 2004, 395, 210. (19) Raff, L. M.; Malshe, M.; Hagan, M. T.; Doughan, D. I.; Rockley, M. G.; Komanduri, R. J. Chem. Phys. 2005, 122, 084104.
MD Investigation of BeH + H2 f BeH2 + H (20) Le, H. M.; Huynh, S.; Raff, L. M. J. Chem. Phys. 2009, 131, 014107. (21) Prudente, F. V.; Acioli, P. H.; Neto, J. J. S. J. Chem. Phys. 1998, 109, 8801. (22) Agrawal, P. M.; Raff, L. M.; Hagan, M. T.; Komanduri, R. J. Chem. Phys. 2006, 124, 134306. (23) Malshe, M.; Narulkar, R.; Raff, L. M.; Hagan, M. T.; Bukkapatnam, S.; Agrawal, P. M.; Komanduri, R. J. Chem. Phys. 2009, 130, 184102. (24) Malshe, M.; Narulkar, R.; Raff, L. M.; Hagan, M. T.; Bukkapatnam, S.; Komanduri, R. J. Chem. Phys. 2008, 129, 044111. (25) Agrawal, P. M.; Samadh, A. N. A.; Raff, L. M.; Hagan, M. T.; Bukkapatnam, S.; Komanduri, R. J. Chem. Phys. 2005, 123, 1. (26) Bukkapatanam, S.; Malshe, M.; Agrawal, P. M.; Raff, L. M.; Komanduri, R. Phys. ReV. B 2006, 74, 224102. (27) Malshe, M.; Pukrittayakamee, A.; Hagan, M. T.; Raff, L. M.; Bukkapatnum, S.; Komanduri, R. J. Chem. Phys. 2009, 131, 124127. (28) (a) Collins, M. A.; Bettens, R. P. A. Phys. Chem. Chem. Phys. 1999, 1, 939. (b) Collins, M. A.; Zhang, D. H. J. Chem. Phys. 1999, 111, 9924. (c) Bettens, R. P. A.; Collins, M. A. J. Chem. Phys. 1999, 111, 816. (29) Thompson, K. C.; Collins, M. A. Faraday Trans. 1997, 93, 871. (30) Ischtwan, J.; Collins, M. A. J. Chem. Phys. 1994, 100, 8080. (31) Head-Gordon, M.; Pople, J. A.; Frisch, M. J. Chem. Phys. Lett. 1988, 153, 503. (32) Frisch, M. J.; Head-Gordon, M.; Pople, J. A. Chem. Phys. Lett. 1990, 166, 275.
J. Phys. Chem. A, Vol. 114, No. 1, 2010 53 (33) Frisch, M. J.; Head-Gordon, M.; Pople, J. A. Chem. Phys. Lett. 1990, 166, 281. (34) Head-Gordon, M.; Head-Gordon, T. Chem. Phys. Lett. 1994, 220, 122. (35) Saebo, S.; Almlof, J. Chem. Phys. Lett. 1989, 154, 83. (36) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; et al. Gaussian 98, Revision A.1; Gaussian, Inc.: Pittsburgh, PA, 2001. (37) Zhang, D. H.; Zhang, J. Z. H. J. Chem. Phys. 1994, 101, 1146. (38) Zhang, D. H.; Light, J. C.; Lee, S. Y. J. Chem. Phys. 1998, 109, 79. (39) Zhang, D. H.; Lee, S. Y. J. Chem. Phys. 1998, 110, 4435. (40) Hornik, M.; Stinchcombe, M.; White, H. Neural Networks 1989, 2, 359. (41) Lawrence, S.; Giles, C. L. Int. Jt. Conf. Neural Networks 2000, 1, 1114. (42) Collins, M. A. Theor. Chem. Acc. 2002, 108, 313. (43) Rahaman, A.; Raff, L. M. J. Phys. Chem. A 2001, 105, 2156. (44) Guan, Y.; Thompson, D. L. Chem. Phys. 1989, 139, 147. (45) Raff, L. M.; Thompson, D. L. Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC: Boca Raton, FL, 1985; Vol. 3, p 1. (46) Raff, L. M.; Porter, R. N.; Thompson, D. L.; Sims, L. B. J. Am. Chem. Soc. 1970, 92, 3208. (47) Raff, L. M.; Sims, L. B.; Thompson, D. L.; Porter, R. N. J. Chem. Phys. 1970, 53, 1606.
JP907507Z
