Summary, Conclusions, and Future Trends

Since the introduction of classical and semiclassical molecular dynamics (MD) methods in the 1960s and Gaussian procedures to conduct electronic structure calculations in the 1970s, a principal objective of theoretical chemistry has been to combine the two methods so that MD and quantum mechanical studies can be conducted on ab initio potential surfaces. Although numerous procedures have been attempted, the goal of first principles, ab initio dynamics calculations has proven to be elusive when the system contains five or more atoms moving in unrestricted three-dimensional space. For many years, the conventional wisdom has been that ab initio MD calculations for complex systems containing five or more atoms with several open reaction channels are presently beyond our computational capabilities. The rationale for this view are (a) the inherent difficulty of high level ab initio quantum calculations on complex systems that may take numerous, large-scale computations impossible, (b) the large dimensionality of the configuration space for such systems that makes it necessary to examine prohibitively large numbers of nuclear configurations, and (c) the extreme difficulty associated with obtaining sufficiently converged results to permit accurate interpolation of numerical data obtained from electronic structure calculations when the dimensionality of the system is nine or greater. Neural networks (NN) derive their name from the fact that their interlocking structure superficially resembles the neural network of a human brain and from the fact that NNs can sense the underlying correlations that exist in a database and properly map them in a manner analogous to the way a human brain can execute pattern recognition. Artificial neurons were first proposed in 1943 by Warren McCulloch, a neurophysiologist, and Walter Pitts, an MIT logician. NNs have been employed by engineers for decades to assist in the solution of a multitude of problems. Nevertheless, the power of NNs to assist in the solution of numerous problems that occur in chemical reaction dynamics is just now being realized by the chemistry community.

Empirical Potential-Energy Surfaces Fitting Using Feed forward Neural Networks

When the system of interest becomes too complex to permit the use of ab initio methods to obtain the system potential-energy surfaces (PES), empirical potential surfaces are frequently employed to represent the force fields present in the system under investigation. In most cases, the functional forms present in these potentials are selected on the basis of chemical and physical intuitions. The parameters of the surface are frequently adjusted to fit a very small set of experimental data that comprise bond energies, equilibrium bond distances and angles, fundamental vibrational frequencies, and perhaps measured barrier heights to reactions of interest. Such potentials generally yield only qualitative or semiquantitative descriptions of the system dynamics. Several research groups have significantly improved the accuracy of the values of the experimental properties computed using empirical potential surfaces by fitting the chosen functional form for the potential to the force fields obtained from trajectories using ab initio Car-Parrinello molecular dynamics simulations. The fitting to the force fields is usually done using a least-squares fitting approach. This method has been employed by Izvekov et al. to obtain effective non-polarizable three-site force fields for liquid water. Carré et al. have employed such a procedure to obtain a new pair potential for silica. In their investigation, the vector of potential parameters was fitted using an iterative Levenberg-Marquardt algorithm. Tangney and Scandolo have also developed an interatomic force field for liquid SiO2 in which the parameters were fitted to the forces, stresses, and energies obtained from ab initio calculations. Ercolessi and Adams have used a quasi-Newtonian procedure to fit an empirical potential for aluminum to data obtained from first-principals computations. Empirical potentials can be improved by making the parameters parameterized functions of the coordinates defining the instantaneous positions of the atoms of the system. This approach has been successfully employed by numerous investigators The difficulty with this procedure is that the number of parameters that must be adjusted increases rapidly. Appropriate fitting of these parameters requires a much more extensive database. Finally, the actual fitting process can often be tedious, difficult, and time-consuming.

Configuration Space Sampling Methods

In order to achieve the maximum accuracy in characterizing the PES and the associated force fields for an MD investigation, careful preparation of the database is an essential step in the process. The points that must be addressed include the following: 1. The total volume of configuration space is extremely large, and its size increases as the internal energy of the system rises. For example, consider a four-atom system. For this system, at least six internal coordinates must be specified to determine the spatial configuration of the molecular system. At a given internal energy, each of these six coordinates can span a continuous range of values from some minimum to some maximum. If each variable range is divided into 100 equal increments and the potential energy of the system computed by some ab initio method for all possible configurations of the system, a total of 1006or 1012 electronic structure calculations would need to be executed. This is clearly beyond the computational capabilities of any computational system currently in existence. Grid sampling methods can and have been used effectively for three atom systems. However, for more complex systems, it is essential that procedures be developed that permit the regions of configuration space that are important in the reaction dynamics to be identified. 2. Sampling methods usually should be optimized to produce a reasonably uniform density of data points in those regions of configuration space that are important in the dynamics. If this is not done and there are regions of very high point density and others with low point density, no fitting technique will function well. The parameters of the method will adjust themselves to fit regions of high density preferentially over those with low density even when the low-density regions may be more important in the dynamics. An exception to the need to have an approximately uniform density of points in the database occurs in regions where the potential gradient is large. In such regions, the density of points in the database will need to be larger than in regions in which the gradient is small.

Neural Network Methods for Data Analysis and Statistical Error Reduction

The use of neural networks (NNs) to predict an outcome or the output results as a function of a set of input parameters has been gaining wider acceptance with the advance in computer technology as well as with an increased awareness of the potential of NNs. A neural network is first trained to learn the underlying functional relationship between the output and the input parameters by providing it with a large number of data points, where each data point corresponds to a set of output and input parameters. Sumpter and Noid demonstrated the use of NNs to map the vibrational motion derived from the vibrational spectra onto a PES with relatively high accuracy. In another application, Sumpter et al. trained an NN to learn the relation between the phase-space points along a trajectory and the mode energies for stretching, torsion, and bending vibrations of H2O2. Likewise, Nami et al. demonstrated the use of NNs to determine the TiO2 deposition rates in a chemical vapor deposition (CVD) process from the knowledge of a range of deposition conditions. In view of the success achieved in obtaining interpolated values of the PESs for multi-atomic systems using an NN trained by the ab initio energy values for a large number of configurations, it is reasonable to ask whether we can successfully compute the results of an MD trajectory for a chemical reaction using an NN trained by the data obtained by previous MD simulations. If this can be done successfully, it becomes possible to execute a small number of trajectories, M, and then utilize the results of these trajectories as a database to train an NN to predict the final results of a very large number of trajectories N, where N >> M, that can be used to increase the statistical accuracy of the MD calculations and to further explore the dependence of the trajectory results upon a wide variety of variables without actually having to perform any further numerical integrations. In effect, the NN replaces the computationally laborious numerical integrations.

Genetic Algorithm (GA) and Internal Energy Transfer Calculations Using Neural Network (NN) Methods

Genetic algorithms (GA), like NNs, can be used to fit highly nonlinear functional forms, such as empirical interatomic potentials from a large ensemble of data. Briefly, a genetic algorithm uses a stochastic global search method that mimics the process of natural biological evolution. GAs operate on a population of potential solutions applying the principle of survival of the fittest to generate progressively better approximations to a solution. A new set of approximations is generated in each iteration (also known as generation) of a GA through the process of selecting individuals from the solution space according to their fitness levels, and breeding them together using operators borrowed from natural genetics. This process leads to the evolution of populations of individuals that have a higher probability of being “fitter,” i.e., better approximations of the specified potential values, than the individuals they were created from, just as in natural adaptation. The most time-consuming part in implementing a GA is often the evaluation of the objective or the fitness function. The objective function O[P] is expressed as sum squared error computed over a given large ensemble of data. Consequently, the time required for evaluating the objective function becomes an important factor. Since a GA is well suited for implementing on parallel computers, the time required for evaluating the objective function can be reduced significantly by parallel processing. A better approach would be to map out the objective function using several possible solutions concurrently or beforehand to improve computational efficiency of the GA prior to its execution, and using this information to implement the GA. This will obviate the need for cumbersome direct evaluation of the objective function. Neural networks may be best suited to map the functional relationship between the objective function and the various parameters of the specific functional form. This study presents an approach that combines the universal function approximation capability of multilayer neural networks to accelerate a GA for fitting atomic system potentials. The approach involves evaluating the objective function, which for the present application is the mean squared error (MSE) between the computed and model-estimated potential, and training a multilayer neural network with decision variables as input and the objective function as output.

Applications of Neural Network Fitting of Potential-Energy Surfaces

In this chapter, several examples of NN fitting of databases obtained using either ab initio electronic structure methods or an empirical potential will be discussed. The objective of this presentation is not to provide a complete and comprehensive review of the field nor is it to acquaint the reader with the details of the reaction dynamics of the particular systems employed as examples. It is rather to provide a clear picture of the power and limitations of NN methods for the investigation of reaction dynamics. We begin with a brief overview of the literature in the field. Neural networks provide a powerful method to effect the fitting of an ensemble of potential energy points in a database. In 1993, Blank et al. employed an NN to fit data derived from an empirical potential model for CO chemisorbed on a Ni(111) surface. Two years later, these same investigators also examined the interaction potential of H2 on a Si(100)-2 × 1 surface using a data set comprising 750 energies computed using local density functional theory. To the best of our knowledge, these were the first two examples in which NNs were employed to provide the PES for a dynamics study. Hobday et al. have investigated the energies of C-H systems by using a Tersoff potential form in which the three-body term is replaced by an NN comprising five input nodes, one hidden layer with six nodes, and an output layer. In this work, the five input elements are computed by consideration of the bond type, i.e., C-C or C-H, the three-body bond angle θ, which is input to the NN in the form (1 + cos θ)2, the connectivity of the local environment, and the second neighbor information. The method was applied to carbon clusters and a wide variety of alkanes, alkenes, alkynes, aromatics, and radicals. Comparison of the atomization energies obtained using the NN potential surfaces with experimental values showed the errors for 12 alkanes, 13 alkenes, 4 alkynes, 7 aromatics, and 12 radicals to lie in the ranges zero to 0.3 eV (alkanes), 0.1 to 1.5 eV (alkenes), 0 to 0.5 eV (alkynes), zero to 1.0 eV (aromatics), and zero to 2.8 eV (radicals).

Fitting Potential Energy Hypersurfaces

Molecular dynamics (MD) and Monte Carlo (MC) simulations are the two most powerful methods for the investigation of dynamic behavior of atomic and molecular motions of complex systems. To date, such studies have been used to investigate chemical reaction mechanisms, energy transfer pathways, reaction rates, and product yields in a wide array of polyatomic systems. In addition, MD/MC methods have been successfully applied for the investigation of gas-surface reactions, diffusion on surfaces and in the bulk, membrane transport, and synthesis of diamond using chemical vapor deposition (CVD) techniques. The structure of vapor deposited rare gas matrices has been studied using trajectories procedures. If the chemical reaction of interest contains three atoms or fewer, various types of quantum and semiclassical calculations can be brought to bear on the problem. These methods include wave packet studies, close-coupling calculations at various levels of accuracy, and S-matrix theory. Several excellent review articles have been published describing the principal techniques and problems involved in conducting MD studies; the reader may wish to consult these as background material for this discussion. With the advent of relatively inexpensive, powerful personal computers, MD/MC simulations have become routine. Once the potential-energy hypersurface for the system has been obtained, the computations are straightforward, though time-consuming. In the majority of cases, the computational time required is on the order of hours to a few days. However, the accuracy of these simulations depends critically on the accuracy of the potential hypersurface used. The major problem associated with MD/MC investigations is the development of a potential-energy hypersurface whose topographical features are sufficiently close to those of the true, but unknown, surface that the results of the calculations are experimentally meaningful. Once the potential surface is chosen or computed, all the results from any quantum mechanical, semiclassical, or classical scattering or equilibrium calculation are determined. The only purpose of the MD calculations is to ascertain what these results are. Therefore, the most critical part of any MD/MC study is the development of the potential-energy hypersurface and the associated force field.

Potential-Energy Surfaces Using Expansion Methods and Neural Networks

Expansion methods have been employed for some time to represent the potentialenergy surface for molecular systems. The basic concept involved with any expansion method is to write the PES expression for an N-atom system that requires the specification of 3N-6 internal coordinates as a sum of terms each of which involves fewer than N atoms and/or fewer than (3N-6) coordinate variables. Two approaches to the implementation of this concept have been suggested. In the first approach, the focus of attention is the number of internal coordinates upon which each term in the expansion depends.

Other Applications of Neural Networks to Quantum Mechanical Problems

The usual method for solving the vibrational Schrödinger equation to obtain molecular vibrational spectra and the associated wave functions generally involves the expansion of the vibrational wave function, ψk(y), in terms of a linear combination of a set of basis functions.

Feed Forward Neural Networks

In this section, we want to give a brief introduction to neural networks (NNs). It is written for readers who are not familiar with neural networks but are curious about how they can be applied to practical problems in chemical reaction dynamics. The field of neural networks covers a very broad area. It is not possible to discuss all types of neural networks. Instead, we will concentrate on the most common neural network architecture, namely, the multilayer perceptron (MLP). We will describe the basics of this architecture, discuss its capabilities, and show how it has been used on several different chemical reaction dynamics problems (for introductions to other types of networks, the reader is referred to References 105-107). For the purposes of this document, we will look at neural networks as function approximators. As shown in Figure 3-1, we have some unknown function that we wish to approximate. We want to adjust the parameters of the network so that it will produce the same response as the unknown function, if the same input is applied to both systems. For our applications, the unknown function may correspond to the relationship between the atomic structure variables and the resulting potential energy and forces. The multilayer perceptron neural network is built up of simple components. We will begin with a single-input neuron, which we will then extend to multiple inputs. We will next stack these neurons together to produce layers. Finally, we will cascade the layers together to form the network. A single-input neuron is shown in Figure 3-2. The scalar input p is multiplied by the scalar weight w to form wp, one of the terms that is sent to the summer. The other input, 1, is multiplied by a bias b and then passed to the summer. The summer output n, often referred to as the net input, goes into a transfer function f, which produces the scalar neuron output a.