Gaussian Process Model for Extrapolation of Scattering Observables for Complex Molecules: from Benzene to Benzonitrile
We consider a problem of extrapolating the collision properties of a large polyatomic molecule A-H to make predictions of the dynamical properties for another molecule related to A-H by the substitution of the H atom with a small molecular group X, without explicitly computing the potential energy surface for A-X. We assume that the effect of the H X substitution is embodied in a multidimensional function with unknown parameters characterizing the change of the potential energy surface. We propose to apply the Gaussian Process model to determine the dependence of the dynamical observables on the unknown parameters. This can be used to produce an interval of the observable values that corresponds to physical variations of the potential parameters. We show that the Gaussian Process model combined with classical trajectory calculations can be used to obtain the dependence of the cross sections for collisions of CHCN with He on the unknown parameters describing the interaction of the He atom with the CN fragment of the molecule. The unknown parameters are then varied within physically reasonable ranges to produce a prediction uncertainty of the cross sections. The results are normalized to the cross sections for He - CH collisions obtained from quantum scattering calculations in order to provide a prediction interval of the thermally averaged cross sections for collisions of CHCN with He.
This work is motivated by recent measurements of the thermalization dynamics of benzonitrile CHCN in the cold gas of He (1). The interpretation of the experiments requires the knowledge of the thermally averaged cross sections for collisions of CHCN with He. However, the potential energy surface for He – CHCN interactions is, at present unknown and the accurate quantum dynamical calculation of the cross sections for He – CHCN collisions is a formidable task. On the other hand, the interaction of a closely-related molecule CH with He is well characterized and accurate quantum dynamical calculations of the cross sections for CH – He collisions have been previously reported (2); (3). This raises the question: given the cross sections for CH – He scattering is it possible to make predictions of the cross sections for CHCN – He collisions?
The experimentally measurable cross sections for elastic or inelastic collisions of molecules in the gas phase can be computed by means of quantum scattering or classical trajectory calculations (4). Any such calculation uses as input the interaction potential energy surface (PES) that governs the molecular dynamics. In general, the PES must be computed for every collision system before the dynamical calculations. The collision dynamics is typically sensitive to details of the PES and the collision properties of different molecules must be obtained from independent scattering calculations with the corresponding PESs. However, for large polyatomic molecules, the substitution of one atom or molecular group with another molecular group may alter only a small part of the global PES. In this case, instead of calculating the PES and the dynamical properties for the molecule with the substituted molecular group, one may consider to determine the effects of the substitution on the collision properties. This would allow one to make predictions about the collision properties of a specific molecule based on the known collision observables for another molecule.
Extrapolating the collision properties between different molecular systems is important for multiple applications. First of all, it can be used to significantly increase the range of molecules amenable to rigorous scattering theory analysis. Consider, for example, the collision properties of benzene CH and benzene substitutes CH-X, where X is a halogen atom or a molecular group. The symmetry of benzene reduces the numerical complexity of the quantum scattering calculations to a great extent (2). The absence of this symmetry in CH-X makes the scattering calculations much more computationally difficult, often impossible. It is, therefore, desirable to develop an approach for predicting the scattering observables for benzene substitutes, given the computed scattering observables for CH. Second, it is often necessary to compare the collision properties of different molecules for calibrating the experimental observations and/or the design of new experiments (5); (6); (7). However, the direct measurements for some molecular species may often be difficult or impossible.
The extrapolation of the collision observables from one polyatomic molecule to another is a difficult task. Consider, for example, the interaction of CH and CH-X with a He atom. The substitution of H by X distorts the highly symmetric PES for the He - molecule interaction. If X is a halogen atom, it may be possible to evaluate the effect of the substitution by computing the scattering cross sections for collisions of CHX with He on a distorted benzene - He PES as functions of the distortion within a physically reasonable range of distortions. However, the problem becomes much more complex if X is a molecular group, such as CH or CN. The introduction of the molecular group X changes the number of the internal degrees of freedom and the distortion of the PES becomes a function of many parameters, determined by the relative positions of the atoms in the molecular group X. For example, if X is CN, the distortion of the PES is a function of 6 parameters. It is impossible to analyze the effect of the H CN substitution on the scattering observables by direct scattering calculations on a grid of these six parameters.
To overcome this problem, we propose to apply the Gaussian Process (GP) model (8); (9); (10), used for machine learning applications in engineering technologies (11); (12); (13). Given a set of data depending simultaneously on multiple parameters, the GP model determines the correlations between the data points in the entire data sample and provides a highly efficient, non-parametric interpolation between the data points in the multi-dimensional space. The Gaussian Process model is known to yield good prediction accuracy, when trained by 10 to 100 data points per dimension (14). An accurate six-dimensional hypersurface of scattering observables can thus be obtained with only 60 to 600 dynamical calculations at randomly chosen values of the 6 PES parameters. Here, we show that this can be exploited to provide a prediction interval of the scattering observables for a molecule A-X, given the scattering observable of the molecule A-H, where A and X are some molecular groups. To provide predictions useful for the experiments in Ref. (1), we use the GP model for the extrapolation of the scattering cross sections for benzene - He collisions to determine the prediction range of the scattering cross sections for benzonitrile - He collisions.
Ii Gaussian Process Model of Scattering Cross sections
This work extends our previous contribution (15), where we proposed to apply the GP model for the analysis of the sensitivity of the scattering dynamics of complex molecules to details of the PES. In the present work we assume that the scattering cross sections of a molecule A-H are known, either from an experiment or from a quantum dynamics calculation. The goal is to formulate the methodology to explore the effects of the substitution H X on the scattering cross sections without explicitly computing the PES for A-X. We begin with a general formulation and specialize the problem to the collision systems of benzene and benzonitrile with He in the next section.
We assume that X is a molecular group consisting of atoms. For simplicity, we consider the interaction of the molecules A-H and A-X with a structureless atom Rg and assume that the interaction of A-H and A-X with Rg is described by a single adiabatic PES. This implies that both of the molecules and the Rg atom have closed electronic shells. The substitution of a hydrogen atom with an -atom molecular group adds degrees of freedom. The change of the PES due to this substitution depends on the coordinates specifying the position of each atom in the group X relative to Rg. The PES for the interaction of the molecule A-X with Rg can be most generally written as
where is the PES for the interaction of the molecule A-H with Rg, the vector specifies the separation between the position of the center of mass of the molecule A-X and the position of the Rg atom, and is a vector with components that represent the positions of the individual atoms in the molecular group X relative to . The allowed range of the coordinates is restricted by the parameters of the individual bonds in the molecular group X.
A rigorous calculation of the scattering cross sections for collisions of A-X with Rg requires the evaluation of the PES (1) as a function of and , followed by the scattering calculations. This is generally a prohibitively difficult task that we want to avoid. Instead, we propose to determine the range of collision cross sections given physically reasonable variations of the PES as a function of the coordinates . The specific implementation of this procedure is illustrated in the next section.
The coordinates can be written as , where are the coordinates of the atoms in the equilibrium geometry of the molecular group X, and describe the deviations from the equilibrium geometry. The equilibrium geometry is generally unknown so the equilibrium coordinates must be treated as unknown parameters. The value of the PES change at each point of and can be restricted on physical grounds to be within . For the methodology proposed here, it is necessary to represent by a series of analytical functions characterized by a finite set of parameters , such that the variation of the individual parameters within a certain range changes within . The change of the PES in Eq. (1) thus becomes a function of these parameters: . Since are assumed to be unknown, of the parameters in the vector must correspond to .
With this formulation, the problem is reduced to determining the variation of the scattering cross sections as functions of . If the computation of the PES is to be avoided, one could – in principle – compute the scattering cross sections as functions of the individual parameters , thus producing a prediction interval of the scattering observables. Assuming that can be parametrized by ten parameters in addition to values of the equilibrium positions , the total number of parameters is . This ranges from 16 for a two-atom molecular group X to 22 for a four-atom molecular group. It is clearly unfeasible to perform dynamical scattering calculations, whether quantum or classical, on a grid of these 22 parameters. We propose to use the GP model to determine the dependence of the cross sections on , in order to find the prediction interval of the cross sections corresponding to the range of the PES variations .
We consider the scattering cross section as a function of parameters described by vector . The components of the vector are the collision energy, the internal energy of the molecule A-X and the parameters . The cross sections can be calculated by means of a classical trajectory method (16) at fixed values of , fixed values of the collision energy and well-defined internal energies of the molecule. Given the calculated values of at a small number of randomly chosen values of , we determine the correlations between different values of . These correlations are then used to train a GP model in order to make predictions of the cross sections for arbitrary values of within a given interval. The number of computed cross sections necessary to make accurate predictions can be estimated to be the number of parameters (14). With 22 parameters and the collision and internal energies as the additional free parameters, the number of values required to determine an accurate dependence of on should be about 240. Whether this number of calculations yields sufficient accuracy can be easily tested by computing the cross sections at new, arbitrary values of and comparing the computed values with the GP model predictions.
A GP can be viewed as a family of random functions normally distributed around a mean function. We denote the GP by . The realization of a GP at a particular site of the multi-parameter space is a value of a function randomly drawn from this family of functions and evaluated at . In the following, we will refer to as a GP or as a Gaussian random function. A GP is completely defined by its mean function and a covariance function . We assume a GP with a constant variance so that , where is a correlation function. The random outputs at fixed thus form a normal distribution with mean and variance and the multiple outputs at different jointly follow a multivariate normal distribution (17); (18).
The scattering cross section at each value of is assumed to be a realization of a Gaussian process . Given a set of cross sections computed at a finite number of values , the goal is to find and . Since the correlation function is fitted to calculated data, the choice of is somewhat flexible. We use the following analytical form for the correlation function, known as the Matérn correlation function (19); (20); (21); (22):
where , is the modified Bessel function of order and are the unknown parameters representing the characteristic length scales of the correlation variations. We fix the value of to , which reduces the correlation function to
We find that the Matérn correlation function yields more accurate results for the present application than the popular Gaussian correlation function (19); (20); (21); (22) we used in the previous work (15). With the Gaussian correlation function, the GP is differentiable to any order. With the Matérn correlation function with parameter , the process is differentiable to order (). Thus, with , the GP is twice differentiable.
The mean of the Gaussian random function can be modelled as
where is a vector of regression functions (9); (10) and is a vector of unknown coefficients. If the dependence of the cross sections on any of the parameters in is known, the regression functions can be chosen to mimic this dependence. This can make the GP model more efficient, i.e. fewer cross section values may be required to achieve the desired level of accuracy at arbitrary values of parameters. However, we note that the GP model does not rely on any specific form of the regression functions in Eq. (4). In the present work, we assume that and , which reduces Eq. (3) to a single unknown parameter . The problem is thus reduced to finding the parameters , and that provide the most accurate correlation function.
The GP model analysis begins with the computation of the cross sections at input vectors randomly chosen to cover the allowed interval of the parameters. These vectors are referred to as the training sites. The multiple outputs of a GP at the training sites follow a multivariate normal distribution
with the mean vector and the covariance matrix . Here, is an design matrix with regressors for each training site as the matrix elements
and is a matrix defined as
where the hat over the symbol denotes the MLE. To find the MLE of , we maximize the log-likelihood function
numerically by an iterative computation of the determinant and the matrix inverse .
The goal is to make a prediction of the cross section at an arbitrary input vector , given the values of the cross sections at the training sites. The values obtained by multiple realizations of the GP at and the multiple outputs of the GP at training sites are jointly distributed as
where is a column vector specified by the correlation function with the MLE of . This means that the conditional distribution of possible values given the values is a normal distribution
with the conditional mean and variance given by
Note that the conditional mean given by Eq. (13) is shifted from the unconditional mean equal, in our case, to due to the point-to-point correlations . Replacing in Eq. (13) with the vector of the known cross section values , we obtain the GP model prediction for the value of the cross section at
Iii From Benzene to Benzonitrile
In this section we apply the GP model to predict the cross sections for elastic scattering of benzonitrile (CHCN) with He, required for the interpretation of the experiments on cooling polyatomic molecules in a buffer gas of He (1). The PES for CHCN - He is currently unknown. The quantum scattering calculations of cross sections for collisions of CHCN with He are prohibitively difficult. At the same time, the interaction of benzene (CH) with He can be accurately described using a semi-empirical bond-additive method (23); (24) and the cross sections for collisions of CH with He can be computed using an accurate coupled states approach (25); (26) based on the time-independent quantum scattering theory (2). In the present section, we show how the GP model can be used to obtain the range of the cross sections for He - CHCN collisions, given the cross sections for He - CH collisions.
We use the semi-empirical approach of Pirani and coworkers (23); (24) for constructing the PES for the interactions of the polyatomic molecules with He. This method treats a polyatomic hydrocarbon molecule as an ensemble of C-C and C-H bonds and represents the PES for the molecule - He interaction as a sum of pairwise He - CH bond and He - CC bond interaction energies, optimized based on accurate ab initio calculations and measurements of bond polarizabilities. This method of representing the PES is particularly well suited for the analysis of the substitution on the molecular scattering properties using the GP model described in the previous section.
Within the approach of Pirani and coworkers (24), each He-CH and He-CC bond interaction is represented by the following analytical function:
is the distance between the atom and the center of the bond, and is the reduced distance , where is the position of the potential well and is the energy at the bottom of the potential well. The parameters and represent the well depth and location for the parallel and perpendicular approaches of the He atom to the corresponding bond and is the angle between the bond axis and the vector connecting the He atom to the center of the bond axis. The eight parameters specifying the PES for the He - CH interaction are given in Table 1.
In order to construct the PES for the He - CHCN interaction, we use the same approach. However, the PES for the He - CHCN system must include, in addition to the aromatic CC - He and CH - He interactions, the parameters describing the interaction of the He atom with the non-aromatic single CC bond and the CN bond fragment. The parameters for the He - non-aromatic, single CC bond interaction are known from the work in Refs. (24); (27) and are listed in Table I. The parameters for the He - CN bond interaction are, at present, unknown and we propose to treat them as variable parameters. The PES for He - CHCN collisions is thus parametrized by twelve parameters listed in Table 1 and four unknown parameters , , and characterizing the locations and energies of the well depths arising from the perpendicular and parallel approaches of the He atom to the CN bond fragment in the molecule. We treat these parameters as variables in the input vector of the GP model. The vector is thus a five-dimensional (5D) vector, containing four parameters , , and and the collision energy. In this work, we fix the internal energy of the molecule to correspond to the ground rotational state.
Given that the size of the N atom is smaller than the size of the C atom and that the C-N bond is more polarized than the C-C bond, one should expect that the He - CN interaction parameters should lead to a larger potential depth but smaller equilibrium distance than the parameters of the He - CC bond interactions. This puts the upper limit on the values of and and the lower limit on the values of and . In order to explore the dependence of the collision cross sections on the variation of these parameters, we construct 100 different PESs for the He - CHCN system with the variable parameters in the following ranges: Å, Å, meV, meV.
The substitution of the H atom in benzene with the molecular group CN is expected to introduce a local change of the global PES. To illustrate this, we plot in Figure 1 the cross sections of the PES for the He - CHCN and He - benzene interactions for various geometries of the atom - molecule approach. Figure 1 illustrates that the CN molecular group leads to a significant change of the potential energy only when the He atom approaches the CN side of the molecule. This justifies the representation (1) of the global PES and the approach adopted here.
Since the quantum scattering calculations of cross sections for He - CHCN collisions are prohibitively difficult, we use the classical trajectory method developed in Refs. (28); (2); (29) for the dynamical calculations for benzonitrile. Using the classical trajectory computations, as described in detail in Ref. (28), we compute the scattering cross sections for 100 combinations of different PESs and collision energies. Each computation represents an average of 5000 trajectories. The cross sections are computed using the corrected version of Eq. (23) in Ref. (2); (3). These cross sections are then used to train the GP model in order to provide the global dependence of the cross sections on the collision energy and the four unknown interaction potential parameters. The resulting ranges of the cross sections are then scaled by the ratio , where is the thermally averaged cross section for He - benzene collisions computed with the coupled states method (25); (26) and is the same cross section computed with the classical trajectory method used for benzonitrile.
Our first goal is to obtain the range of the cross sections corresponding to a physical variation of the He - CN bond interaction parameters , , , and the atom - molecule collision energy. Figure 2 shows the variation of the cross sections as a function of one of the PES parameters, with the remaining PES parameters and the collision energy chosen for each point at random. As evident from Figure 2, the scattered points cannot be assumed to follow any well-defined dependence.
We obtain the global dependence of the cross sections on the four PES parameters and the collision energy by computing the cross sections at 100 points, placed randomly and quasi-uniformly, in this 5D parameter space. These points are used to train the GP model as discussed in the previous section. Given the GP model, we compute the global 5D surface illustrated in Figure 3. In order to demonstrate the accuracy of the global 5D surface, we computed the cross sections at a different set of 100 points randomly chosen in the 5D parameter space. Figure 4 compares the predicted values with the computed values for these 100 points.
In order to quantify the error of the GP model, we compute the empirical root mean squared error (ERMSE) defined as
and the scaled root mean squared error (SRMSE) defined by . In Eq. (19), is the number of the cross section values, represents the values of the computed cross sections and – the value predicted by the GP model. For the model with only 50 scattering calculations used as training points, Å and . If the number of the scattering calculations is increased to 100, the errors decrease to Å and .
Given the GP model of the collision cross sections trained by the computations at random values of , , , and the atom - molecule collision energy, we can use Eq. (15) to examine the variation of the collision energy dependence of the cross sections as the four parameters are varied within the physical ranges specified in the previous section. This produces a band of the cross section values for each value of the collision energy presented in Figure 5. The results illustrate that the physical variation of the PES parameters describing the interaction of the He atom with the CN group changes the cross sections by less than %, with the percentage defined as the difference between the maximum and minimum values divided by the mean value. The computation of the cross sections for benzonitrile - He collisions with an accurate PES must fall within the grey area of Figure 5 with the probability 95 %.
Since the GP model provides the global dependence of the cross sections on the underlying parameters, it is possible to perform the analysis of the sensitivity of the cross section variations to the individual PES parameters by using the functional analysis of variance decomposition (30); (31); (32). The results shown in Figure 6 illustrate that, of the four unknown parameters, has the biggest effect on the variation of the cross sections. Figure 6 also shows that the effect of the individual parameters is largely correlated with the values of the other parameters. This means that the prediction interval in Figure 5 must be obtained by the simultaneous variation of each of the underlying parameters and cannot be obtained by a simple variation of one of the parameters with the other parameters kept fixed at a few random values.
The final results of this work are presented in Figure 7, showing the thermally averaged cross sections for benzene - He collisions computed with the quantum CS approach (blue squares) and the interval of the thermally averaged cross sections for benzonitrile - He collisions obtained by integrating the results in Figure 5 with the Maxwell-Boltzmann distribution of collision energies and scaling the thermally averaged cross sections by the ratio . Figure 7 shows that the presence of the molecular group CN enhances the elastic scattering cross sections by a factor of 1.1 – 1.5. The uncertainty of the parameters in the He - CN group interaction leads to the uncertainty in the final thermally averaged cross section. This percentage is defined as the difference between the maximum and minimum values of the grey band in Figure 7 divided by the corresponding mean value.
In the present work we consider a general problem of extrapolating the known collision properties of a complex polyatomic molecule A-H to make predictions for another molecule related to A-H by the substitution of the H atom with a small molecular group. Using the example of CH – He and CHCN – He collision systems, we show that the H CN substitution leads to a local modification of the global potential energy surface. While the quantitative effect of the H CN substitution on the PES is unknown, it can be parametrized by a finite number () of parameters , collectively denoted by . With this parametrization, the effect of the H CN substitution on the collision observable is embodied in a multidimensional function . Even if the parameters are unknown, once the function of on is determined, it can be used to obtain the prediction interval for given physically reasonable variation of .
We propose to apply the Gaussian Process model to determine . The model requires about 10 () calculations of the scattering observables to provide accurate results. Thus, with , an accurate five-dimensional dependence of on and on the collision energy can be obtained with about 50 dynamical calculations. We showed that this procedure can be used to obtain an accurate dependence of the cross sections for collisions of CHCN with He on the parameters describing the interaction of the He atom with the CN fragment of the molecule. The results are then compared with the cross section for He - CH collisions known from the quantum scattering calculations in order to provide a prediction interval of the thermally averaged cross sections for collisions of CHCN with He. This allowed us to obtain the prediction interval (Figure 7) of the collision cross sections for CHCN – He collisions without the knowledge of the potential energy surface.
This work illustrates that the Gaussian Process model can be used for a variety of applications in molecular dynamics research. For example, once the function is obtained, it can be used to perform the sensitivity analysis (such as in Figure 6) to determine which of the PES parameters is more or less important in determining the collision observable . The dependence of on the collision energy, which is treated as one of the unknown model parameters, can be used for an efficient integration of the collision properties to obtain thermally averaged observables. The function can be used to integrate out the dependence of the collision observables on the PES parameters, thereby minimizing the uncertainties associated with inaccuracies of the PES calculations and providing the error bars associated with the uncertainties of the PES calculations.
The computational effort associated with training the Gaussian Process model is determined by matrix inversion and scales as the third power of the number of training points. Given that the number of training points required for accurate predictions is typically 10 times the number of unknown parameters and that it has now become routine to invert matrices with the dimension of 1000 1000, the methodology proposed here can be applied to determine the dependence of the scattering observables on up to 100 unknown parameters. It is thus easy to envision an application where the Gaussian Process model is used to characterize the dependence of a dynamical observable on all PES parameters describing the interaction of two polyatomic molecules. This dependence can be used to determine which of the atoms in the two molecules are important for the detectable outcome of the molecule - molecule interaction and which of the atoms can be parametrized by simple functions irrelevant for the outcome of the interaction.
We thank Dr. Massimiliano Bartolomei for sending us the parameters of the interaction potentials ans for allowing us to include his unpublished data in Table I. This work is supported by NSERC of Canada.
- D. Patterson and J. M. Doyle, Phys. Chem. Chem. Phys. 17, 5372 (2015).
- Z. Li, R. V. Krems, and E. J. Heller, J. Chem. Phys. 141, 104317 (2014).
- Note that Eq. (23) of Ref. (2) incorrectly includes the factor 4. The correct equation is . The computations reported in Ref. (2) used the correct equation.
- J. R. Taylor, Scattering Theory: the Quantum Theory of Nonrelativistic Collisions (Dover Publications, 2006).
- J. Küpper, F. Filsinger, and G. Meijer, Faraday Discuss. 142, 155 (2009).
- C. Sommer, L. van Buuren, M. Motsch, S. Pohle, J. Bayerl, P. Pinkse, and G. Rempe, Faraday Discuss. 142, 203 (2009).
- D. Patterson, E. Tsikita, and J. M. Doyle, Phys. Chem. Chem. Phys. 12, 9736 (2010).
- J. Sacks, S. B. Schiller, and W. J. Welch, Technometrics 31, 41 (1989).
- T. J. Santner, B. J. Williams, and W. I. Notz, The Design and Analysis of Computer Experiments (Springer Science Bussiness Media, New York, 2003).
- C. E. Rasmussen and C. K. I. Williams, Gaussian Process for Machine Learning (The MIT Press, Cambridge, 2006).
- D. Higdon, M. Kennedy, J. C. Cavendish, J. A. Cafeo, and R. D. Ryne, SIAM J. Sci. Comput. 26, 448 (2004).
- D. Higdon, J. Gattiker, B. Williams, and M. Rightley, J. Am. Statist. Assoc. 103, 570 (2008).
- R. B. Gramacy and H. K. H. Lee, J. Am. Statist. Assoc. 103, 1119 (2008).
- J. L. Loeppky, J. Sacks, and W. J. Welch, Technometrics 51, 366 (2009).
- J. Cui and R. V. Krems, accepted for publication in Phys. Rev. Lett. (arXiv:1503.01432v2).
- R. B. Bernstein, Atom-Molecule Collision Theory: A Guide for the Experimentalist, (Plenum, New York, 1979).
- R. J. Adler, The Geometry of Random Fields, (SIAM,1981).
- H. Cramér and M. R. Leadbetter, Stationary and Related Stochastic Processes: Sample Function Properties and Their Applications, (Courier Corporation, 2013).
- T. Mitchell, M. Morris, and D. Ylvisaker, Stoch. Proc. Appl. 35, 109 (1990).
- N. A. C. Cressie, Statistics for Spatial Data, (Wiley-Interscience, New York, 1993).
- M. L. Stein, Interpolation of Spatial Data: Some Theory for Kriging, (Springer Science Business Media, 1999).
- M. Abt, Scand. J. Stat. 26, 563 (1999).
- F. Pirani, D. Cappelletti, G. Liuti, Chem. Phys. Lett. 350, 286 (2001).
- F. Pirani, M. Albertt́, A. Catro, M. Moix Teixdidor, and D. Cappelletti, Chem. Phys. Lett. 394, 37 (2004).
- B. J. Garrison, and W. A. Lester, Jr., J. Chem. Phys. 66, 531 (1977).
- S. Green, J. Chem. Phys. 64, 3463 (1976).
- M. Bartolomei, private communication (October 2012).
- Z. Li and E. J. Heller, J. Chem. Phys. 136, 054306 (2012).
- J. Cui, Z. Li, and R. V. Krems, J. Chem. Phys. 141, 164315 (2014).
- A. Saltelli, K. Chan, and E. M. Scott, Sensitivity Analysis (Wiley, New York, 2009).
- A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, and S. Tarantola, Global Sensitivity Analysis: the Primer (John Wiley Sons, 2008).
- O. Roustant, D. Ginsbourger, and Y. Deville, J. Stat. Softw. 51, 1 (2012).