Letter pubs.acs.org/JPCL
High-Dimensional Atomistic Neural Network Potentials for Molecule−Surface Interactions: HCl Scattering from Au(111) Brian Kolb,†,‡ Xuan Luo,§ Xueyao Zhou,§ Bin Jiang,*,§ and Hua Guo*,† †
Department of Chemistry and Chemical Biology, University of New Mexico, Albuquerque, New Mexico 87131, United States Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, United States § Department of Chemical Physics, University of Science and Technology of China, Hefei, Anhui 230026, China ‡
ABSTRACT: Ab initio molecular dynamics (AIMD) simulations of molecule−surface scattering allow first-principles characterization of the dynamics. However, the large number of density functional theory calculations along the trajectories is very costly, limiting simulations of long-time events and giving rise to poor statistics. To avoid this computational bottleneck, we report here the development of a high-dimensional molecule−surface interaction potential energy surface (PES) with movable surface atoms, using a machine learning approach. With 60 degrees of freedom, this PES allows energy transfer between the energetic impinging molecule and thermal surface atoms. Classical trajectory calculations for the scattering of DCl from Au(111) on this PES are found to agree well with AIMD simulations, with ∼105-fold acceleration. Scattering of HCl from Au(111) is further investigated and compared with available experimental results.
C
DFT calculations. Thus, it is advantageous to construct analytical potential energy surfaces (PESs) using DFT data to achieve higher efficiency for MD calculations. Indeed, DFT-based PESs for molecule−surface interactions have been developed by many authors5,6,20,21 but have so far been mostly limited to the rigid surface approximation, which disallows energy transfer between the impinging molecule and surface phonons. Progress has been made to include the motion of surface atoms. For example, a reaction force field approach was advanced for the dissociative chemisorption (DC) of CH4 on metals.22 More recently, an impressive high-dimensional atomistic PES was developed by Natarajan and Behler to describe the liquid−solid interface between water and Cu surfaces.23 Machine learning, which is an artificial intelligence approach that extracts patterns from large data sets, has been used in many fields in physics and chemistry.24−29 Its applications to represent high-dimensional PESs with large sets of ab initio points have recently witnessed rapid growth.30−33 In this Letter, we report, as a proof of concept, the first high-dimensional accurate first-principles PES for molecular scattering from a metal surface with movable surface atoms, using the Behler− Parrinello atomistic neural network (AtNN) approach.24,34 Such an approach differs fundamentally from force-field-based PESs22 as NNs assume no particular functional form in representing the PES and are thus extremely flexible in representing DFT points. Compared to the existing AtNN molecule−surface
ollisional dissociation of molecules on metal surfaces represents the initial and often rate-limiting step in many industrially important heterogeneous catalytic processes. Hence, it is not surprising that many recent surface dynamics studies have focused on a fundamental understanding of such processes.1−4 The collisional process of a light molecule such as H2 on surfaces is impulsive and can often be treated within a rigid surface approximation as the energy transfer to surface phonons is negligible.5,6 However, this is not the case for heavier molecules, such as H2O, CH4, CO2, and HCl, because they readily transfer energy to surface atoms upon collision.7,8 The amount of energy transfer can be estimated by the Baule model based on collision of two hard spheres,9 namely, ΔE = 4 μ/(1 + μ)2Ei, where Ei is the incidence energy and μ is the ratio between the effective masses of the molecule and surface atom. A quantitative characterization of the scattering, dissociation, and energy transfer can now be achieved by ab initio molecular dynamics (AIMD),10 in which the classical nuclear motion is driven by quantum mechanical forces computed on-the-fly using density functional theory (DFT). Indeed, there has recently been a growing interest in AIMD simulations of molecular scattering and dissociation on metal surfaces.11−19 AIMD requires a large number of DFT calculations, making it computationally expensive. Typically, only relatively short (∼1 ps) dynamical events can be simulated, and even for such processes, only a small (∼103) number of AIMD trajectories is affordable, making a good statistical sampling difficult. Yet, these AIMD trajectories traverse the same coordinate space for scattering or dissociation, thus involving duplicating © XXXX American Chemical Society
Received: December 20, 2016 Accepted: January 19, 2017 Published: January 19, 2017 666
DOI: 10.1021/acs.jpclett.6b02994 J. Phys. Chem. Lett. 2017, 8, 666−672
Letter
The Journal of Physical Chemistry Letters
with limited sampling), the configuration space was first partitioned based on the distance of the HCl from the surface (ZHCl) into three regions: asymptotic (ZHCl > 4.25 Å), weakly interacting (2.5 Å < ZHCl < 4.25 Å), and strongly interacting (ZHCl < 2.5 Å). For each of these regions containing N points, a subset of a given size (n) was selected using a Monte Carlo scheme. The algorithm proceeds by first choosing n points randomly from the region. The fictitious “energy” of this subset is defined as
PESs,23 our PES has to cover a large coordinate space far away from equilibrium because of the high energy of the impinging molecule. The scattering of energetic HCl from Au(111) has been extensively investigated using a molecular beam approach in high vacuum,35−40 which revealed large energy loss of the scattered HCl and vibrational inelasticity. AIMD simulations of the scattering dynamics have confirmed significant energy transfer, though quantitative differences exist.18,19 This system thus serves as an ideal benchmark for developing DFT-based PESs. The AtNN PES for HCl/DCl scattering from Au(111) is developed using the Behler−Parrinello approach, which expresses the total energy as a sum of atomic contributions24
E=
∑ Ei i
n
ϵ=
∑ i,j
1 |rij|2
(2)
The idea is to find the subset of size n that minimizes this “energy”. At random, a point in the subset is exchanged with a point outside of the subset (but still in the relevant region), and the change in energy Δϵ is computed. Note that the computation of Δϵ requires only 2n distance computations. The swap is kept with probability according to the Metropolis criterion
(1)
The atomic contributions (Ei) are approximated by NNs that include information of the environment. The AtNN approach has two distinct advantages.41 First, it always satisfies the permutation symmetry as the contribution to the total energy of each atom in its chemical environment is considered separately. This also guarantees smooth connections at the boundaries of two symmetry-equivalent regions. Second, the complexity of the NNs increases with the number of atomic types, not the number of coordinates. Consequently, it is well suited for large systems, particularly those with few atomic types. Indeed, its applicability has been demonstrated for many extended and molecular systems.41,42 Very recently, this approach has also been extended to reactive systems involving bond breaking and formation.43,44 This work represents the first attempt to extend AtNNs to describe molecular scattering from a solid surface. The generation of the DFT energies used here has been described previously for DCl scattering from Au(111).18 Briefly, non-spin-polarized DFT calculations using the RPBE functional45 were performed within the Vienna Ab Initio Simulation Package (VASP)46,47 using a plane-wave cutoff of 350 eV, and the core electrons were treated using the projector augmented wave (PAW) approach.48 The Brillion zone was sampled with a Monkhorst−Pack 4 × 4 × 1 k-point mesh,49 and the first-order Methfessel−Paxton smearing50 was used with σ = 0.2 eV. The surface model consisted of four layers of Au atoms (each 3 × 3), with the bottom two layers fixed at bulk geometry and the top two fully movable. Structures were generated via AIMD using a 1 fs time step. Initial conditions were chosen by placing the DCl molecule 5 Å above the surface with a randomly selected angular orientation. Initial velocities for D, Cl, and the moving Au atoms were selected randomly from a Maxwell−Boltzmann distribution at 300 K, with a center-of-mass velocity in the Z direction replaced by that corresponding to the predetermined incidence energy. One thousand trajectories were run at each incidence energy (1.25, 2.0, and 2.5 eV) with propagation lengths of 400, 350, and 300 fs for the three energies, respectively. The DFT data so generated consisted of around 900 000 structures, just over 800 000 of which correspond to nonreacted geometries. (The dissociation channel will be considered in the future.) While a fit could have been performed with the entire set, it was desirable to select a subset that was representative of the whole. To this end, a Monte Carlo method was devised to cull the points. To select a subset of the data that maximizes the coverage of the relevant region of configuration space (e.g., minimize sampling clumps or regions
⎧ 1 Δϵ ≤ 0 ⎪ P(Δϵ) = ⎨ ⎛ Δϵ ⎞ ⎟ Δϵ > 0 ⎪ exp⎜⎝ − T ⎠ ⎩
(3)
The fictitious “temperature” T is a parameter in units of 1/Å2 that starts out large, allowing many unfavorable swaps to take place, but is slowly reduced as the simulation progresses, increasingly requiring maintained swaps to be energetically favorable. In the present case, the temperature was lowered from an initial value of 0.1 to a final value of 0.0008 by a factor of 0.95 at each step. Twenty-five thousand structure swaps were attempted at each temperature. A total of 115 000 points were selected by this method, 25 000 from the asymptotic region and 45 000 each from the strongly and weakly interacting regions. These points were used to fit the AtNN PES, which was carried out using the opensource code PROPhet.51 In this approach, a unique NN is created for each atom type in the system and the chemical environment of each atom is described by a vector of two- and three-body mapping functions (called atomic-centered symmetry functions in Behler’s work34), each averaging over interactions with atoms within a chosen cutoff radius (Rc).24,34 Specifically, two-body mapping functions are defined as Gtwo‐body = exp( −α(R ij − R s)2 )fc (R ij)
(4)
while three-body functions take the form Gthree‐body =
1 (1 + λ cos(θ ))ξ exp(−α(R ij 2 + R ik 2))fc (R ij)fc (R jk) 2ξ − 1
(5)
where the cutoff function ⎧ ⎡ ⎛ ⎞ ⎤ 1 πR ⎪ ⎪ ⎢cos⎜ ⎟ + 1⎥ R < R c ⎥⎦ fc (R ) = ⎨ 2 ⎢⎣ ⎝ R c ⎠ ⎪ ⎪ 0 R ≥ Rc ⎩
(6)
ensures that the mapping functions smoothly go to 0 at the chosen cutoff distance and the parameters α, Rs, ξ, λ, and Rc define the specific mapping functions. These mapping functions are summed over all atoms in the chemical environment of the central atom to produce the input to the NN. 667
DOI: 10.1021/acs.jpclett.6b02994 J. Phys. Chem. Lett. 2017, 8, 666−672
Letter
The Journal of Physical Chemistry Letters
the 3 × 3 unit cell, but the scattering at 2.5 eV was not included because it is too close to the energy limit of the DFT data set used in the fitting. A 0.5 fs time step was used in the velocity Verlet algorithm, which converged the total energy within 0.5 meV on average for the majority of trajectories, better than that in AIMD simulations. As shown in Figure 2, there is significant
The PES for HCl on Au(111) was constructed from NNs containing 2 hidden layers of 55 nodes each. Other NN structures have been tried but found not to produce dramatically different outcomes. The cutoff radius (Rc) was taken to be 4.25 Å. A total of 12 two-body functions were constructed with α = 19.2811 Å−2 and Rs values on a uniform grid ranging from 0 to 4.25 Å. A total of six three-body functions were chosen with α = 0.2549 Å−2, λ = ±1, and ξ = 1, 8.5, and 16. Training was conducted with the resilient backpropagation (Rprop) algorithm52−54 in the early stages and switched to the limited-memory Broyden−Fletcher− Goldfarb−Shanno (BFGS) algorithm55 as the minimum was neared. The root-mean-square (RMS) error for fitting the training set was used as the penalty function, and the final value was 15 meV. The PES was validated by calculating the energies of the ∼700 000 structures not used in the fitting and comparing them to the corresponding DFT values, resulting in a RMS error of 11 meV. Note that this error is actually lower than the fitting error, demonstrating that the subset chosen for fitting is indeed representative of the whole. In Figure 1, several cuts of the PES are shown for HCl approaching the surface at different sites, with the surface atoms frozen at equilibrium geometries and HCl being perpendicular to the surface with H down. Not surprisingly, the top site is most repulsive, while the bridge and hollow sites are more accessible. In Figure 1d, a cut of the PES with HCl (rHCl = 1.288 Å) at ZHCl = 2.5 Å is displayed in X and Y coordinates. To validate the fit, we have first redone MD calculations of the DCl scattering using the AtNN PES in exactly the same conditions as the AIMD calculations. Using an adapted VENUS package,56 classical trajectory (CT) calculations were carried out on the analytical PES. One thousand trajectories were run at each incidence energy (1.25 and 2.0 eV), randomly covering
Figure 2. Kinetic energy distributions of the impinging (red) and scattered (blue) DCl obtained from AIMD18 and CT calculations at two collision energies, with Ts = 300 K.
energy loss in the scattered DCl from the surface, and the final kinetic energy distributions agree very well with AIMD results, demonstrating the accuracy of the AtNN PES. As discussed in our previous publication,18 the fraction of kinetic energy loss is 33% and 40% for incidence energies of 1.25 and 2.0 eV, respectively. The acceleration over the AIMD calculations is on the order of 105.
Figure 1. HCl−Au(111) PES as a function of the distance of the molecular center above the surface (Z) and HCl bond length (r) on the top (a), bridge (b), and hollow (c) sites and as a function of the lateral position of the molecular center (X and Y) at Z = 2.5 Å and r = re = 1.288 Å (d). The surface is fixed at equilibrium, and HCl is perpendicular to the surface with H down. 668
DOI: 10.1021/acs.jpclett.6b02994 J. Phys. Chem. Lett. 2017, 8, 666−672
Letter
The Journal of Physical Chemistry Letters
confirm the trend observed in experiment: Increasing the incidence energy leads to a broader kinetic distribution in which the peak shifts to a higher energy. Quantitatively, the calculated mean energy loss, ranging from ∼23% of incidence energy at 0.28 eV to ∼37% at 1.27 eV, is consistent with the mean energy loss of 30% at 1.4 eV reported in a recent AIMD study of HCl scattering from Au(111) using the PBE functional.19 However, these theoretical predictions are much smaller than the measured ones, ∼55% regardless of the incidence energy.39 This becomes clearer in Figure 3h, where the theoretical and experimental mean kinetic energies of the scattered HCl are compared as a function of incidence energy, along with the Baule limit (dotted line). The calculated mean kinetic energy of scattered molecules is too high by 0.1−0.2 eV depending on the incidence energy, compared to experimental counterparts, which are close to the Baule limit. To further understand the effect of surface temperature on the energy transfer process, we calculated the kinetic energy distribution of the HCl(ν = 0 → 0) scattering with the 1.27 eV incidence energy at several Ts values using the same PES. Figure 3g shows that mean energy loss depends very weakly on Ts, but the distribution becomes broader with increasing Ts, which is understandable as higher Ts samples more nonequilibrium configurations, resulting in larger fluctuations of the energy exchange with the molecule. These results suggest that the surface temperature has a relatively minor effect on the scattering, consistent with recent AIMD findings.19 In the upper panel of Figure 4, the angular distribution of the scattered HCl is displayed for several incident energies.
The availability of the PES allows us to explore the scattering of HCl from Au(111), aimed at a more in-depth understanding of the state-to-state scattering experiments of Wodtke and co-workers.39 To this end, quasi-classical trajectory (QCT) calculations were performed on the same PES. Differing from CT, the vibrational energy of HCl was determined by semiclassical quantization of the classical action,57 using the asymptotic diatomic potential of the PES. Trajectories were generated from the ground rovibrational state of HCl located 6 Å above the surface with its orientation and position randomly sampled. The surface configurations were sampled as follows.13 The equilibrium coordinates and momenta of surface atoms were initiated with the Maxwell−Boltzmann distribution at a given surface temperature (Ts) and equilibrated for 1.5 ps with the velocities of surface atoms scaled to the desired Ts every 150 fs. This was followed by a 0.5 ps run, which verified that the ensemble temperature was close to the desired one. Subsequent surface configurations were then chosen for the collision simulation, in which the molecule was allowed to impinge on the surface with an incidence angle identical to that in the experimental setup (θ = 3°).39 The time step for the surface equilibration was 0.5 fs, while 0.1 fs was used for the scattering process to converge to a total energy of
