# Gaussian Process Model for Collision Dynamics of Complex Molecules

###### Abstract

We show that a Gaussian Process model can be combined with a small number (of order 100) of scattering calculations to provide a multi-dimensional dependence of scattering observables on the experimentally controllable parameters (such as the collision energy or temperature) as well as the potential energy surface (PES) parameters. For the case of Ar - CH collisions, we show that 200 classical trajectory calculations are sufficient to provide a 10-dimensional hypersurface, giving the dependence of the collision lifetimes on the collision energy, internal temperature and 8 PES parameters. This can be used for solving the inverse scattering problem, the efficient calculation of thermally averaged observables, for reducing the error of the molecular dynamics calculations by averaging over the PES variations, and the analysis of the sensitivity of the observables to individual parameters determining the PES. Trained by a combination of classical and quantum calculations, the model provides an accurate description of the quantum scattering cross sections, even near scattering resonances.

###### pacs:

The reliable scattering calculations of dynamical properties of molecules are required in almost any research field related to molecular physics. In particular, the experiments on collisional cooling of molecules to cold and ultracold temperatures cold-molecules (), chemical reaction dynamics resonances-1 (), the development of new pressure standards kirk-1 (); kirk-2 (), astrophysics and astrochemistry astrophysics () rely on accurate calculations of molecular collision cross sections. Currently, there are two major problems with the ab initio calculations of molecular dynamics observables. The first problem is the inaccuracy of the potential energy surfaces (PES). Unfortunately, even the most sophisticated quantum chemistry calculations produce the PES with uncertainties that lead to significant (and often unknown) errors in the dynamical calculations. This sensitivity to PES inaccuracies is especially detrimental for low temperature applications (cold molecules, ultracold chemistry, astrophysics and pressure standards) hutson-pes (); hutson-pes2 (); nh-nh (). The second problem is related to the numerical complexity of the quantum dynamics calculations ct (); mott-massey (). For complex molecules with many degrees of freedom, accurate dynamical calculations are extremely time-consuming and it is often impossible to compute enough results for accurate averaging over the collision or internal energies of the colliding partners.

In the present work we propose a solution to these two problems. In order to account for the PES uncertainties, the dynamical results can be averaged over variations of the PES. If the computed observables are averaged over variations of each individual PES parameter, producing an expectation interval of the observables, the ab initio dynamical calculations can have fully predictive power (with error bars). However, the outcome of a molecular collision is generally a complicated function of many (ten or more) PES parameters. It is impossible to obtain the dependence of the collision observables on the individual PES parameters by the direct scattering calculations. We show that such a dependence can be obtained by a combination of a small number (on the order of 100) of scattering calculations with a Gaussian Process (GP) model sacks1989 (); rasmussen2004gaussian (). We show that the same model can be used to obtain the accurate dependence of the scattering observables on the collision or internal energies of the molecules, with a small number of scattering calculations. The result is an accurate global dependence of the scattering observables on the collision energy, internal energy and every individual parameter of the PES surface. This global dependence can be used to average the computed observables over variations of the individual PES parameters, as well as over the collision and internal energies in order to produce thermally average observables. It can also be used to analyze the influence of the individual PES parameters on the scattering outcome. This makes the model proposed here a unique tool for the analysis of the effects of the PES topology on the molecular scattering dynamics.

Widely used in engineering technologies gp-applications-1 (); gp-applications-2 (), the GP model can be viewed as a technique for interpolation in a multi-dimensional space. We choose the GP model because it is an efficient non-parametric method. There is no need to fit data by analytical functions so the model is expected to work for any distribution of scattering observables and to become more accurate when trained by more computed observables. Given the scattering observables computed at a small number of randomly chosen points in the multi-dimensional parameter space, the GP model learns from correlations between the values of these scattering observables to produce a smooth dependence on all the underlying parameters. As an illustrative example, we consider the scattering of benzene molecules CH by rare gas (Rg) atoms He - Xe. The PES surface for CH - Rg interactions is characterized by 8 parameters. We consider two scattering observables buffer-gas-cooling (); bgc-2 (); bgc-3 (): the collision lifetimes and the scattering cross sections. We address the following questions: how many scattering calculations are sufficient to train a GP model to produce an accurate global dependence on all the underlying parameters? Can the GP model be used to make predictions of the scattering observables for one collision system based on the known properties of another collision system? Can the GP model be used to characterize the scattering observables near quantum resonances?

We consider a scattering observable as a function of parameters described by vector . The components of the vector can be the collision energy, the internal energy and/or the parameters representing the PES. We assume that is known from a classical or quantum dynamics computation at a small number of values. Our first goal is to construct an efficient model that, given a finite set of , produces a global dependence of the scattering observable on . If the observable is known from a measurement or a rigorous quantum calculation as a function of some parameters – e.g., the collision energy – we show that the model can be adjusted to produce the global dependence of on that reproduces the accurate data, even if the dynamical calculation method is inaccurate.

We assume that the scattering observable of interest at any is a realization of a Gaussian process , characterized by a mean function , constant variance and correlation function . For any fixed , is a value of a function randomly drawn from a family of functions Gaussian-distributed around . Consequently, the multiple outputs and at and jointly follow a multivariate normal distribution defined by , , and adler1981geometry (); cramer2013stationary (). We assume the following form for the correlation function mitchell1990existence (); cressie1993statistics (); stein1999interpolation (); abt1999estimating ():

(1) |

and write

(2) |

where is a vector of regression functions regressors (), is a vector of unknown coefficients, and is a Gaussian random function with zero mean. The problem is thus reduced to finding , and .

We spread input vectors evenly throughout a region of interest and compute the desired observable at each with a classical or quantum dynamics method. The outputs of a GP at these points follow a multivariate normal distribution with the mean vector and the covariance matrix . Here, is a design matrix with th row filled with the regressors at site , and is a matrix with the elements .

Given , the maximum likelihood estimators (MLE) of and have closed-form solutions sacks1989 ():

(3) |

(4) |

To find the MLE of , we fix and maximize the log-likelihood function

(5) |

numerically by an iterative computation of the determinant and the matrix inverse .

The goal is to make a prediction of the scattering observable at an arbitrary . Because the values at and the outputs at training sites are jointly distributed, the conditional distribution of possible values given the values is a normal distribution with the conditional mean and variance

(6) | |||||

(7) |

where is specified by the now known correlation function . Eq. (6) provides the GP model prediction for the value of the scattering observable at

To illustrate the applicability and accuracy of the GP model, we first compute the collision lifetimes of benzene molecules with Rg atoms zhiying (); cui2012collisional (); zrr (). We use the classical trajectory (CT) method described in Ref. zhiying (). As shown in Ref. potentials (), the CH - Rg PES can be expressed as a sum over terms describing the interaction of Rg with the C-C and C-H bond fragments, characterized by parameters. We first fix the PES parameters to describe the CH - Ar system and focus on the dependence of the lifetimes on two parameters: the collision energy and the rotational temperature . Figure 1 shows the results of the CT calculations illustrating that the collision lifetime exhibits an inverse correlation with , while no apparent correlation with . Figure 1 (c) shows the global surface of the lifetime as a function of and obtained from the GP model with and set to 1.95. To quantify the prediction accuracy of the GP model, we calculate the errors and , where are the computed values and are the GP model predictions. For the model with only 20 scattering calculations used as training points, ps and . With the number of the scattering calculations increased to 50, the errors decrease to ps and .

The scattering calculations presented in the upper panels of Figure 1 cannot be interpreted to assume any simple functional form. In addition, the vastly different gradients of the and dependence may make the conclusions based on calculations at fixed values of one of the parameters misleading. In contrast, the surface plot in Figure 1(c) clearly illustrates that the collision lifetimes decrease monotonically with both and . The effect of the rotational temperature is much weaker especially when cm and there is no strong two-way interaction between and . The GP model surface can be used to evaluate thermally averaged collision lifetimes by integrating the -dependence at given .

The GP model can be extended to multiple collision systems for the predictions of the collision properties of a specific collision system based on the known collision properties of another system. To illustrate this, we consider the lifetimes of the long-lived complexes formed by benzene in collisions with Rg atoms He – Xe. As the collision system is changed from CH - He to CH - Xe, there are two varying factors that determine the change of the collision dynamics: the reduced mass and the PES.

As before, we use the GP model , with now representing the atomic mass and the interaction strength at the global minimum of the atom - molecule PES obtained by scaling the Ar - CH PES. We fix K and cm, and compute the collision lifetimes at 40 randomly chosen points in the interval of and g/molg/mol, which covers all of the Rg – CH systems. These 40 calculation points are then used to train the GP model to produce the surface plot shown in Figure 1 (d). The error of the surface is %. The plot reveals that increasing both and enhances the collision lifetimes and that the reduced-mass dependence of the collision lifetimes is very weak compared to the dependence on the interaction strength.

The GP model can be exploited to explore the role of the individual PES parameters on the observables. To illustrate this, we now consider that contains parameters giving the analytical form of the Rg - CH PES potentials (), in addition to and . We calculate the lifetimes at 200 randomly selected points in this parameter space and use these points to train the GP model. Figure 2 compares the predicted values with the calculated values for another set of 70 randomly selected points. The plot corresponds to the model error .

The 10-parameter GP model contains a wealth of information on the dependence . For example, one can perform a sensitivity analysis by using the functional analysis of variance decomposition pujol2008sensitivity (); saltelli2008global (); roustant2012dicekriging () to determine, which of the PES parameters have the strongest impact on the observable (right inset of Figure 2). Of the 8 PES parameters, the location of the potential well due to the interactions of Rg with the C-C bonds for the parallel approach potentials () is the most important factor determining the collision lifetime. The model can also use be used to compute the uncertainties due to global variation of the PES. Figure 2 (left inset) shows the interval of the lifetimes obtained by the simultaneous variation of all PES parameters.

We now consider the applicability of the GP model to quantum scattering calculations. The quantum results are often affected by resonances resonances-1 (); resonances-2 (), leading to wild variations of the scattering observables in a small range of the underlying parameters. If applied directly to such the case, the GP model is unstable because steep variation of the correlations leads to singularities in hankin (). This is illustrated in Figure 3, showing the GP model predictions trained directly by 60 quantum calculations of cross sections for rotationally inelastic He - CH scattering, randomly chosen at between 1 and 10 cm. The instability of the GP model arises from the wild variations of the scattering cross sections near resonances. We repeated these calculations for the elastic and state-resolved rotationally inelastic cross sections shown in Figure 4 (a-c) of Ref. zrr (). In each case, we found that the wild variation of the quantum results leads to unstable GP model predictions.

However, the GP model can be extended to model the time-consuming quantum scattering calculations with the help of efficient classical dynamics calculations. To do this, we introduce a more complex GP as kennedy2001bayesian ()

(8) |

where and are independent Gaussian random functions, with characterizing the difference between the CT and QM calculations and effectively describing the inaccuracy of the classical trajectory method. The calculations are performed in two steps. First, the CT calculations are used to train the GP model . In the second step, the QM and CT calculations are used together to train the model in Eq. (8), using the parameters of and treating and as variable parameters. This fixes the models and as well as and .

The accuracy of this combined quantum - classical model is illustrated in Figure 3, showing that the model provides an accurate energy dependence of the cross sections, even near scattering resonances. The CT calculations in a two-function model (8) stabilize the model, removing the errors arising from the resonant variation of the quantum results. We applied the two-step model (8) to the calculations for the elastic and state-resolved rotationally inelastic cross sections shown in Figure 4 (a-c) of Ref. zrr () and found a similar improvement in each case.

In summary, we have shown that a Gaussian Process model combined with a small number of scattering calculations can be used to obtain an accurate multi-dimensional dependence of the scattering observables on the experimentally controllable parameters and the PES parameters. Specifically, we showed that the GP model trained only by 20 CT calculations produces a dependence of the CH - Ar collision lifetimes on the collision energy and the rotational temperature of benzene, with the normalized error . Trained by 200 calculations, the GP model produces a 10-dimensional dependence of the collision lifetimes on the collision energy, the rotational temperature and 8 individual PES prameters, with the error %. We have introduced a hybrid GP model that can be trained by a combination of classical and quantum dynamics calculations in order to model the quantum results. We showed that this model works even in the vicnity of quantum scattering resonances, where the direct fit of the quantum results by means of a GP model is unstable. The models described here are expected to find a wide range of applications, from fitting the interaction potentials by solving the inverse scattering problem, to analyzing the dependence of scattering observations on external parameters, to calibrating the accuracy of the scattering calculation methods. For example, the inverse scattering problem can be approached with the help of Eq. (8), where is parametrized by unknown PES parameters and models the experimental data. The best estimates of the unknown PES parameters can then be found by a Markov-chain Monte Carlo method MCMC (), in a procedure similar to one recently applied in Ref. hep ().

###### Acknowledgements.

We thank Dr. Zhiying Li for allowing us to use her codes for the classical and quantum dynamics calculations. This work is supported by NSERC of Canada.## References

- (1) R. V. Krems, W. C. Stwalley, and B. Friedrich (eds.), “Cold molecules: Theory, Experiment, Applications”, CRC Press (2009).
- (2) K. Liu, R. T. Skodje and D. E. Manolopoulos, Phys. Chem. Comm. 5, 27 (2002).
- (3) J. Van Dongen, C. Zhu, D. Clement, G. Dufour, J. L. Booth, and K. W. Madison, Phys. Rev. A 84, 022708 (2011).
- (4) J. L. Booth, D. E. Fagnan, B. Klappauf, K. W. Madison, and J. Wang, “Method and device for accurately measuring the incident flux of ambient particles in a high or ultra-high vacuum environment”, Google Patents, 2011.
- (5) D. Flower, “Molecular Collisions in the Interstellar Medium”, 2nd edition, Cambridge University Press, Cambrdige (2011).
- (6) M. D. Frye and J. M. Hutson, Phys. Rev. A 89, 052705 (2014).
- (7) M. L. González-Martínez and J. M. Hutson, Phys. Rev. Lett. 111, 203004 (2013).
- (8) Y. V. Suleimanov, T. V. Tscherbul, and R. V. Krems, J. Chem. Phys. 137, 024103 (2012).
- (9) R. B. Bernstein (ed.), “Atom-Molecule Collision Theory”, Plenum, New York (1979).
- (10) N. F. Mott and H. S. W. Massey, “The theory of atomic collisions”, Oxford University Press (1965).
- (11) J. Sacks, W. J. Welch, T. J. Mitchell, and H. P. Wynn, Stat. Sci. 4, 409 (1989).
- (12) C. E. Rasmussen, “Gaussian processes in machine learning”, In Advanced lectures on machine learning, Springer (2004).
- (13) D. Higdon, M. Kennedy, J. C. Cavendish, J. A Cafeo, and R. D Ryne, SIAM J. Sci. Comput. 26, 448 (2004).
- (14) D. Higdon, J. Gattiker, B. Williams, and M. Rightley, J. Am. Stat. Assoc. 103, 482 (2008).
- (15) D. Patterson, E. Tsikata, and J. M. Doyle, Phys. Chem. Chem. Phys. 12, 9736 (2010).
- (16) D. Patterson and J. M. Doyle, Phys. Chem. Chem. Phys., in press (2015).
- (17) D. Patterson, M. Schnell, and J. M. Doyle, Nature 497, 475 (2013).
- (18) R. J. Adler, “The geometry of random fields”, SIAM (2010).
- (19) H. Cramér and M. R. Leadbetter, “Stationary and related stochastic processes: Sample function properties and their applications”, Courier Corporation (2013).
- (20) T. Mitchell, M. Morris, and D. Ylvisaker, Stoch. Proc. Appl. 35, 109 (1990).
- (21) N. A. Cressie, “Statistics for spatial data”, Wiley, New York (1993).
- (22) M. L. Stein, “Interpolation of spatial data: some theory for kriging”, Springer (1999).
- (23) M. Abt, Scand. J. Stat. 26, 563 (1999).
- (24) The regression functions can be chosen to follow a specific parametric dependence of the calculated data in order to make the GP model more efficient. If no such dependence is known or can be determined, the regression functions can simply be chosen as and , reducing the first term in Eq. (2) to a constant .
- (25) Z. Li and E. J. Heller, J. Chem. Phys. 136, 054306 (2012).
- (26) J. Cui, Z. Li and R. V. Krems, J. Chem. Phys. 141, 164315 (2014).
- (27) Z. Li, R. V. Krems and E. J. Heller, J. Chem. Phys. 141, 104317 (2014).
- (28) F. Pirani, M. Albertí, A Castro, M. Moix Teixidor, and D Cappelletti, Chem. Phys. Lett. 394, 37 (2004).
- (29) G. Pujol, “Sensitivity: Sensitivity analysis. R package version 1.4–0 (2008).
- (30) A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, and S. Tarantola, “Global sensitivity analysis: the primer”, Wiley (2008).
- (31) O. Roustant, D. Ginsbourger, and Y. Deville, J. Stat. Softw. 51, 54 (2012).
- (32) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
- (33) R. K. S. Hankin, J. Stat. Softw. 14, 16 (2005).
- (34) M. Kennedy and A. O’Hagan, J. Roy. Statist. Soc. Ser. B 63, 425 (2001).
- (35) W. K. Hastings, Biometrika 57, 97 (1970).
- (36) J. Novak, K. Novak, S. Pratt, J. Vredevoogd, C. E. Coleman-Smith, and R. L. Wolpert, Phys. Rev. C 89, 034917 (2014).