Checking sampling distribution of velocities: the KS test

When generating multiple initial geometries or velocities, it is good practice to check that we are sampling in the correct way, i.e that we are uniformily sampling from the distribution we expect or desire to sample from.

In this post I will show an example where the Kolmogorov-Smirnov test (KS-test) is used to check that velocities are sampled uniformily from the Maxwell-Boltzmann distribution function.

From a practical point of view, the test consists in computing the maximum difference \sf D^{max}_n between the cumulative of the probability distribution function of your n-point dataset \sf C^{data}_n and the cumulative of the continuous distribution you believe you are sampling from \sf C^{distr}:

\sf D^{max}_{n}=sup_{x}(|C^{data}_n(x)-C^{distr}(x)|)

When the distributions match, as n tends to infinity, the distance should tend to zero.

Depending on how large this difference is respect to a critical value, one can establish whether their data was sampled from the continuous distribution.

Lets assume we have a set of 1000 values of velocities (norm of vectors) \sf [ v_0,v_1,...,v_i,...,v_n ] we can compute the probability and cumulative distribution of the velocities \sf C^{data}(v_i) . Here below I show a histogram plot of the probability distribution function at T=100K for a set of 1000 sampled velocities:


We then also compute the cumulative distribution of the distribution we expect to match, i.e the Maxwell-Boltzmann distribution,

\sf p(v) = \sqrt{(\frac{m}{ 2 \pi k_B T})^3} \,\,4 \pi v^2 e^{-\frac{mv^2}{2 k_B T}}

and here is a plot showing both cumulatives:


We see that the two cumulatives are very similar. The maximum distance in this case is D = 0.0238454225145 and the critical probability value for this number of points is \sf p_{critical} = 0.0430069761783 (computed as \sf 1.36/\sqrt{n} see e.g. for a 0.05 level of significance)

We have found D < \sf p_{critical} therefore the test confirms the fact that our data’s distribution is the Maxwell-Boltzmann distribution.



Randomly rotated molecules and bimolecular collision

In this post I will show an example where initial molecules are rotated to generate initial conditions for bimolecular collisions. The molecules which are considered are sulfanide and acetone cyanohydrin. These molecules most likely played an important role in the production of origins of life compounds (see B.H. Patel, C. Percivalle, D. J. Ritson, C. D. Duffy & J. D. Sutherland, Nature Chemistry 7, 301–307 (2015))

In Fig. 1, the minimized structure of each molecule is shown.

Fig 1: Geometry of two molecules, on the left acetone cyanohydrin, on the right sulfanide

In Fig. 2, a set of randomly rotated structures of the molecules of Fig. 1 is shown.

Fig 2: Two sets of rotated molecules. The molecules from Fig 1 were all randomly rotated and are shown here in the same plot.

With these initial positions and by defining some initial atomic velocities we can run dynamics for bimolecular collisions. Here this was done using TeraChem ( In the first video here below, I show a set of multiple trajectories (each of these was run separately) where the molecules collide at 300K.

In the second video a single trajectory is shown. From this we can see that the angle at which reactants approach is an important variable in determining whether a reaction will occur. In fact for this trajectory no reaction occurs.


Initial conditions – uniform random rotations

In theoretical chemistry we often simulate dynamics by propagating trajectories on potential energy surfaces. From these trajectories one can derive, e.g. a reaction rate constant. For each trajectories, initial conditions, i.e velocities and coordinates, need to be defined for the molecules. One interesting part of the process is defining uniformly distributed random rotations in three dimensions (uniform distribution on SO(3)).  In fact usually we want to give the molecules a random rotation. For a detailed theoretical description of this, please see e.g. Ref [1].

At first, one may be tempted to extract three random Euler angles \alpha, \beta and \gamma , each from a uniform distribution and to compute the rotation matrix R as \sf{R=R_x(\alpha)\cdot R_y(\beta) \cdot R_z(\gamma)} . If this procedure is repeated, a set of random matrices will be obtained. However, this set will not be uniformly distributed. In fact, if we take a vector v=(1,0,0) and rotate it 1000 times using the above procedure, where angles are extracted from uniform distributions \alpha,\gamma\in[-\pi,\pi] and \beta\in[-\pi/2,\pi/2] we get

We can see a higher density of points around the north and south poles. This happens because in three dimensions, the Haar measure of a matrix R in SO(3) can be written as [dR]=cos\beta d\beta d\alpha d\gamma . The cosine term is fundamental. If we do not take it into account when sampling  \beta , as we get closer to the poles, circles have a smaller circumference than near the equator thus if the \beta angle is sampled uniformly one gets more rotations in those locations (see e.g. Ref. [2]).

One of the solutions was proposed by Arvos (Ref. [3]): “To generate uniformly distributed random rotations of a unit sphere, first perform a random rotation about the vertical axis, then rotate the north pole to a random position.”   This is straightforward to implement and we can see (image below) that indeed we get a uniform distribution of rotations


Also, we can extract the angles from the random matrices and see that the distributions are as expected

Now we can apply this to a molecule and generate uniformly distributed rotations. I will show an image of this in the next post!

[1] David Kirk, Graphic Gems III, San Diego: Academic Press Professional (1992).

[2] John T. Kent and Asaad M. Ganeiber: Simulation of random rotation matrices, Events and Meetings of Italian Statistical Society, 46th scientific meeting of the Italian statistical society (2012).

[3] James Arvo, Fast random rotation matrices, Graphics Gems III, p117-120 (2012).

NEB: Computing a minimum energy path

Nudged Elastic Band (NEB) is a method which can be used to find the minimum energy path (MEP) and the transition state which connects reactants to products. Reactants and products are considered in their local minima and their states need to be provided by the user.

The algorithm works as following. A set of images (i.e a specific geometric configuration of the atoms) connecting the initial (reactant) and final (product) states are linearly interpolated.

Each image is connected to the other by a spring. Then, the energy of this string of images is minimized. The perpendicular component of the spring force is removed to avoid interference with relaxation to the MEP when minimizing. If the minimization converges, one obtains a guess of the minimum energy path.

For example, for the \sf{H+H_2}\rightarrow\sf{H_2+H} reaction, the blue images (see below) would be interpolations between start and end geometries. During minimization the intermediate configurations will change until the MEP is found.


NEB path for a reaction

Here I will describe how to write input files to submit an NEB calculation using TeraChem (

To begin one minimizes the reactant and products to get a converged geometry in their local minima.

Here is an example of the input file for TeraChem minimization:

charge          0
method          rhf
basis           6-31g*
run             minimize
min_coordinates cartesian

When the simulation converges, the output file in scr/ will contain the geometry of reactants at each optimization step. The last frame is the optimized frame.The same type of minimization is repeated for the products. The two optimized geometries are inserted into a file, first the reactant geometry then the product geometry. This file is used as input for the NEB calculation, here is the NEB input file:

charge          0
method          rhf
basis           6-31g*
run             ts
ts_method       neb_frozen
maxit           200
nstep           200
min_image       10

The main difference is the keyword ts_method used to select the NEB method with frozen end points. When the algorithm converges, the TS optimized geometry can be found in the file of the scr/ directory.


Finding the minimum energy path

To run reaction dynamics or compute reaction rate constants ab-initio, one needs a potential energy surface. In the case of transition state theory the minimum energy path needs to be identified and computed and subsequently the transition state needs to be determined. If the transition state is known then the full path need not be determined to implement TST.

The computation of a full potential energy surface is generally time-consuming and complex. Therefore, methods have been developed to try to interpolate between known fixed geometries to compute an energy pathway while avoiding to compute the full potential energy surface. These include the Nudge Elastic Band method (NEB) [1], the String Method [2],  the Growing String Method (GSM), the Freezing String Method (FSM) [4].

In the next post I will describe how to implement the NEB method for a reaction.

[1] G. Henkelman and H. Jónsson,  J. Chem. Phys. 113, 9901-9904 (2000).

[2] E. Weinan, W. Ren, and E. Vanden-Eijnden, Phys. Rev. B 66, 052301 (2002).

[3] B. Peters, A. Heyden, A. T. Bell, and A. Chakraborty, J. Chem. Phys. 120, 7877 (2004).

[4] A. Behn, P. M. Zimmerman, A. T. Bell, and M. Head-Gordon, J. Chem. Phys. 135, 224108 (2011)

Part 2 – Reaction rate constants with TST: rates

We can now insert our partition functions values in the expression for the rate constant at T=300K:

\sf k(T)= \sf \frac{k_BT}{h}\frac{Q^{\ddagger}(T)}{Q_R(T)}e^{-\beta E_a}

To find,

\sf k(T=300K)=  \sf \text{\sf 8.229957e+03 } cm^3 / molecule / s

Note that the overall reaction rate reactant partition function is taken as the product of each reactant partition function. This procedure can be repeated at each temperature to get a profile of the reaction rate constant as a function of temperature. Here is a plot in logarithmic scale where we compare our results (blue curve) to published results (red triangles)


The main role of the temperature dependence is in the exponential, this is why we see a linear dependence of the logarithm of the rate. Further, as expected the rate increases as temperature rises (i.e in the limit of the x axis going to zero) and decreases at lower temperatures. When quantum effects are included, one might see a large change at low temperatures which is due to the presence of tunneling. Other differences arise from the fact that the rate as is comes from a static picture and no dynamics has been carried out.

Part 1 – Reaction rate constant with transition state theory (TST): partition functions

As part of the work carried out to understand origins of life, it is crucial to determine which processes/reactions were most likely to occur. The thermal rate constant of a reaction defines the efficiency of the process at a given temperature. To get an exact rate, one needs an accurate potential energy surface and quantum time-dependent dynamics of reactants on the potential energy surface. This is computationally very demanding.

One simple theoretical approach which avoids having to propagate in time is classical Transition State Theory (TST). For more details on TST see the original papers by Wigner, Eyring, Evans and Polanyi:

  • Pelzer and Wigner, Z. Phys. Chem. B 15, 445 (1932)
  • Eyring, J. Chem. Phys. 3, 105 (1935)
  • Evans and Polyani, Trans. Faraday Soc. 31, 857 (1935)

The main idea of TST is that the reaction dynamics can be described by a one dimensional minimum energy path on the potential energy surface which connects reactants to products through a transition state. The transition state is found at the maximum energy along the one dimensional path and it is a first order saddle point in the multidimensional potential energy surface. The imaginary frequency of the transition state is not included in the computation of the transition state partition function as it corresponds to the reaction coordinate degree of freedom.


The approximations required to derive TST are

  • the Born-Oppenheimer approximation is valid
  • reactants are at equilibrium and in the Boltzmann distribution
  • no recrossing: molecules which have crossed the transition state may not return no recrossing of transition state and classical dynamics of the nuclei
  • at the transition state, the motion along the one dimensional reaction coordinate can be separated from other motion and treated classically as translational motion

The expression for the TST rate constant is then

\sf k(T)= \sf \frac{k_BT}{h}\frac{Q^{\ddagger}(T)}{Q_R(T)}e^{-\beta E_a}

Here we see the ratio of the canonical partition function for the transition state, \sf Q^{\ddagger}, and the canonical partition function for the reactants, \sf Q_R, as well as the exponential of the reaction energy barrier \sf E_a. As long as the ratio of partition functions stays constant as a function of temperature, the rate constant will decay linearly with temperature in log scale.


Let’s see how it works for a simple reaction

CH4 + H → [CH5+]→H2+CH3

Here for simplicity the transition state is indicated in the square brackets. We’ll begin by some general expressions and then implement these for the above reaction.

  • Compute partition functions of reactants and transition state

The most accurate way to compute the partition function of the molecules in their ground state would be to diagonalize the reactant and transition state Hamiltonians and sum over all eigenstates, i.e

\sf Q(T)=\sum_i g_i  e^{-\beta E_i}

With \sf \beta=1/(k_BT) and \sf g_i is the degeneracy of the i-th level. This expression is difficult to evaluate as one needs a full potential energy surface. Diagonalizing the Hamiltonian is expensive and gets worse as the number of degrees of freedom increases. Therefore, further approximations are made.

Assuming we have one of each of these molecules, the overall canonical partition function can be expressed as the product of the partition function for electronic, vibrational, translational and rotational degrees of freedom:

\sf Q(T)=Q_{\textsf{el}}(T)\cdot Q_{\textsf{vib}}(T)\cdot Q_{\textsf{tras}}(T)\cdot Q_{\textsf{rot}}(T)

For N identical molecules this expression needs to be elevated to the power of N and rescaled by N!.

Each partition function can then expressed as

\sf Q_\textsf{tras}(T)= \sf \left(\frac{2\pi \textit{\textsf {m}} \sf k_B T}{h^2}\right)^{3/2}\cdot V%s=2

\sf Q_\textsf{vib}(T)= \sf \prod_i^{n_{\textsf{vib}}}\frac{e^{\sf -\beta\hslash\omega_i/2}}{1-e^{\sf -\beta\hslash\omega_i}}%s=2

\sf Q_\textsf{rot}(T)= \sf \frac{\pi^{1/2}}{\sigma_r}\left(\frac{T^3}{\Theta_A\cdot\Theta_B\cdot\Theta_C}\right)^{1/2} %s=2

with m the mass and V the volume, \sf n_{\textsf{vib}}=3N-6 for non linear molecules and 3N-5 for linear molecules. The rotational constants are obtained from the three moments of inertia \sf I_j  as \sf \Theta_j=\frac{h^2}{8\pi^2 I_j k_b}. The paramerter \sf \sigma_r is the rotational constant which can be obtained from the point group of the molecule. In general if your molecule is not at equilibrium it may be complicated to estimate this value, for more details see e.g. M. K. Gilson and K. K. Irikura, J. Phys. Chem. B, 114:16304 (2010). For a linear molecule the rotational partition function is simply \sf Q_{\textsf{rot}}=\frac{T}{\sigma_r\Theta_r}.

The electronic partition function is approximated to be \sf Q_{\textsf{el}}=g_{el}^{\textsf{ground}} i.e the ground state electronic degeneracy.

Note that in the above expressions for the partition functions we have assumed that

  • the harmonic approximation is valid
  • there is no roto-vibrational coupling
  • the molecules can be seen as linear or non linear rigid rotor
  • the main contribution to the electronic partition function comes from the ground electronic state

For this reaction at T=300K, we find

\sf{ Q_{CH_4}=7.93e+12\text{; }Q_{H}=9.77e+29 \textsf{ and } Q_{CH_5^+}=1.27e+15}

and we have used


\sf{\omega(CH_4)=} \sf{\left(\textsf{1341 1341 1341 1569 1569 3031 3151 3151 3151}\right) cm^{-1}}


\sf{\omega(CH_5^+)=} \sf{\left(\textsf{534 534 1074 1125 1125 1442 1442 1795 3076 3223 3223}\right) cm^{-1}}

excluding the imaginary frequency at the saddle point and

\sf{I(CH_4)}= \sf{ \left(\textsf{3.173   3.173   3.173}\right) au\cdot\AA^2}


\sf{I(CH_5^+)}= \sf{ \left(\textsf{4.401   4.436   4.460}\right) au\cdot\AA^2}

Some of these values (e.g. frequencies) were taken from J. Chem. Phys. 124:164307 (2006). For the electronic partition function we used the spin multiplicity of 2 for hydrogen, 1 for methane and 2 for the transition state.

Next post we’ll put this all in and get the rate constant as a function of temperature.