PyQCAMS: Python Quasi-Classical Atom-Molecule Scattering

We present Python Quasi-classical atom-molecule scattering (PyQCAMS), a new Python package for atom-molecule scattering within the quasi-classical trajectory approach. The input consists of mass, collision energy, impact parameter, and pair-wise interactions to choose between Buckingham, generalized Lennard-Jones, and Morse potentials. As the output, the code provides the vibrational quenching, dissociation, and reactive cross sections along with the rovibrational energy distribution of the reaction products. Furthermore, we treat H$_2$ + Ca $\rightarrow$ CaH + H reactions as a prototypical example to illustrate the properties and performance of the software. Finally, we study the parallelization performance of the code by looking into the time per trajectory as a function of the number of CPUs used.


Introduction
For decades, one of the main approaches to studying molecular dynamics has been via the quasi-classical trajectory (QCT) method [1,2].This technique treats collisions semi-classically.The nuclear dynamics in the underlying potential energy surface is treated classically.However, the initial and final states are picked following the Bohr-Sommerfeld quantization rule, yielding accurate predictions at significantly less cost than quantum dynamics as long as they fall within given conditions regarding the collision energy -number of contributing partial waves [3,4].QCT has been used in multitude of scenarios relevant to chemical physics [5,6,7,8,9,10,11], and cold and ultracold chemistry [12,13,14], ranging from the ultracold to the hyperthermal regimes.In particular, It has been used to study the relaxation and reaction dynamics of cold atom-ionic molecule [13,12] and atom-molecule collisions [15].
We present an open-source, object-oriented program written in Python to perform QCT calculations on atom-molecule systems called Python Quasi-Classical Atom-Molecule Scattering (PyQCAMS).While chemical dynamics programs such as Gaussian [16] and VENUS [17], we introduce a more accessible, user-friendly and dedicated path to performing these simulations.Our program consists of completely open-source software, and relies mainly on the NumPy [18] and SciPy [19] packages.We use matplotlib [20] for data visualization, pandas [21] for data storage and analysis, and multiprocess [22,23] for parallel implementation.
The outline of this paper is as follows: In Section 2 we discuss the theory behind the quasi-classical trajectories method, including initial conditions, trajectory reactions, and analysis.In Section 3, we discuss the PyQCAMS program as separated into the inputs, main code, and outputs.We also discuss the implementation and performance of the program.In Section 4, we provide an example of a typical workflow to study the reaction (H 2 + Ca), where we outline how a user can obtain reaction rates and product distributions using the program.

Theoretical approach
In this section we describe the basics of QCT.The quasi-classical trajectory (QCT) method treats scattering processes semi-classically.First, by solving Newton's equations of motion of the colliding nuclei.Next, in the case of atom-molecule scattering, at the start of each trajectory, the internal degrees of freedom of the molecule are treated within the Wentzel-Kramers-Brillouin (WKB) approximation, such that the initial rovibrational state (v, j) satisfies a quantum-mechanically viable state.The classical Hamiltonian for a 3-particle, atom-molecule system with masses m i , i = 1, 2, 3 takes the form: where p i and r i represent the momentum and position vectors of each atom with respect to the origin.This Hamiltonian is expressed in Jacobi coordinates as: where ρ 1 is the Jacobi vector associated with the molecule and P 1 is its conjugate momentum.ρ 2 is the Jacobi vector connecting the atom with the center of mass of the molecule, and P 2 is its conjugate momentum, as shown in Fig. 1.Note that this Hamiltonian uses the reduced masses of the corresponding atoms: Defining the coordinates in this way separates the center-of-mass degree of freedom from the relative one.Then, the momentum associated with the center-ofmass degrees of freedom is a constant of motion since the interaction potential do not depend on the center-of-mass position, and neglected in the analysis of the dynamics..We trace the atoms' subsequent motion by solving Hamilton's equations of motion: ) Jacobi coordinates of a atom-molecule system.Here, the molecule is rotating in the x-y plane, with the angular momentum J along the z-axis and κ = ρ 1 × ẑ. η is the angle between κ and J, and φ and θ are defined as usual in spherical coordinates.

Initial Conditions
The initial orientation angles are randomly generated to initialize each trajectory.The momentum of the molecule, P 1 , is subsequently defined by the orientation angles.If we initialize the molecule at its classical outer turning point | ρ 1 | = r + , the momentum has no radial component.Since the angular momentum J = r × p, and the fact that ρ 1 ⊥ P 1 , we find the magnitude of P 1 = j(j + 1)/r + , with components [4]: The initial vibrational phase is randomly generated by choosing the initial distance between the atom and molecule as where R 0 is a fixed far-away distance where the interaction potential strength is negligible.The second term probes the molecule's vibrational phase since χ ∈ [0, 1] is randomly generated following a uniform distribution and τ v,j is the vibrational period of the molecule corresponding to the quantummechanical state (v, j).This is calculated as where V (r 12 ) is the molecular potential energy and r 12 is the molecular separation.E int is the internal energy of the molecule, which is calculated following a discrete variable representation (DVR) method using particle-in-a-box eigenfunctions [24].

Reaction Products
For the atom-molecule reaction AB + C, we expect three possible final products: The final product is determined by the relative energy between each atom.The condition for whether two atoms are bound is determined by the effective potential: where r 0 is the local maximum of the effective potential.Here, a, b ∈ [1,3] represent each of the atoms.Two atoms are considered bound only if this condition is satisfied for just one pair of atoms, so that the other two pairs can be deemed unbound.
In the same manner as the initial rovibrational states (v, j), the final states (v , j ) are calculated within the WKB approximation.The rotational quantum number is given by: where J = ρ i × P i is the effective conjugate angular momentum of the Jacobi vector ρ i .The vibrational quantum number is given by: Although these equations lead to continuous numbers, they need to be interpreted as integers before assigning them as quantum numbers.This is done through the Gaussian binning (GB) process, where each v is weighted against its nearest integer v t with a Gaussian [3]:

Analysis
The cross section and reaction rates for an atom-molecule collision are calculated from QCT simulations.The first step in doing so is to calculate the opacity function, which is formally defined as where dΩ = sin θdθdφdηdχ.This integral is evaluated via the Monte Carlo technique over the randomly generated variables θ, φ, η, and χ.
The final state distribution of, for example, a reaction is calculated via the Gaussian Binning method [25]: with the total weight W evaluated as where N r , N q , andN d are the number of reaction, quenching, and dissociation products, respectively.Each Gaussian weight is summed over both the total number of products and the total number of vibrational states produced by the scattering process.Each dissociation reaction has a weight W d = 1 since there is no vibrational state associated with this product.
For N calculations, the opacity function P q,r,d (E c , b) gives the probability of quenching (q), reacting (r), or dissociation (d) and δ is the error associated with the Monte Carlo technique:

The program
PyQCAMS is written in an object-oriented manner, containing three main classes; Potentials, Energy, and QCT.This allows for the addition of new methods within each step of the calculation outlined in Figure 3.The program takes an input file where the user specifies the details of the reaction of interest.The variables are passed into their respective classes, where the trajectory calculations are performed.The results of each trajectory are then output into a csv file.The details of each step are outlined in this section.

Input
The input file is a JSON file whose keywords are described in this section (Listing 1).Each trajectory is tuned by the initial rovibrational energy level of the molecule vi, ji, the collision energy Ec (K), impact parameter b, and initial distance between the atom and molecule r0.For reproducibility, the user can input a seed for the random number generator, or leave it as null to generate a new random number.The interaction potential between each atom is also tune-able, with choices from the Morse, generalized Lennard-Jones, and Buckingham potentials.The user should choose each of these internuclear potentials labeled as "potential AB (BC,CA)".Note that AB should always represent the initial molecule, so that the "C" represents the colliding atom.This convention should also be followed when specifying the masses.
The associated parameters of each of these potentials must be entered by the user, which are outlined in Section 3.2.Since the energy spectrum is calculated via the DVR method, the user should also enter the number of DVR points and specify the range for which the potential is well-defined, such that it captures the potential minimum, short-range, and long-range behavior.
Finally, the parameters for integration should be entered.The user can choose when to stop a trajectory using t stop, which is a multiplicative factor of the trajectory's time scale.Another stopping condition is the maximum distance between any two atoms, r stop.Finally, the absolute and relative tolerances of the integrator can be controlled with a tol, r tol.
Note that all values entered in the input file should be in atomic units except for collision energy, which is in K for the user's convenience.However, when importing the energy into the program, ensure that it is converted to atomic units, as is done in the start() method.

Main Code
The Potentials class houses the interaction potentials.Each method in this class represents a different potential, and returns a tuple of two function objects: the potential and its derivative.This class is useful for studying and manipulating the effect of different potentials on different pairs of the atommolecule system.Results can be sensitive to the chosen potential, so it is important to test different potentials when studying a system.The available potentials and associated parameters are: Parameters required: D e , r e , w e 2. Generalized Lennard-Jones: Parameters required: m, n, C m , C n , and r e (guess left) 3. Buckingham: V (r) = ae −br − C 6 /r 6  Parameters required: a, b, C 6 , and r e (guess left)."max": Guess of where the maximum is.At short range, Buckingham potentials can reach a maximum and collapse.Enter your nearest r value to this maximum.
Note that while r e is not a parameter in the analytical expressions of the generalized Lennard-Jones and Buckingham potentials, it is required as part of a root solver method The Energy class houses the DVR method and turning point calculations.This is a separate class since these quantities should be computed once at the beginning of each trajectory.The spectrum and turning points for different E (v,j) in different potentials can be stored separately for future use.
The QCT class performs the main QCT calculation, and is outlined in Figure 4.The iCond() method uses the calculated parameters from the Energy and Potentials class to generate a random set of initial conditions, yielding ( ρ 1 , ρ 2 , P 1 , P 1 ).The hamEq() function writes the Hamilton's equations of motion, and serves as an input function to solve ivp(), a SciPy [19] integrator utilizing the adaptive Runge-Kutta 5(4) intregration method [26].The vPrime() method serves to calculate Equation 10, yielding the final state vector s f .The integrator is housed in the run T() method, which processes the results and assigns the relevant outputs as attributes of the QCT class.
The hamiltonian() method defines the Hamiltonian of the system, and outputs the total energy and momentum, useful for checking conserved quantities.
Figure 4: Outline of the QCT class.The components within the dashed line are contained in the run T() method of the QCT class.The outputs are the final position, momentum, state, and product count vectors.The product count vector is defined by the final product bins, filled according to which pair of atoms is bound.

Output
For a single trajectory, the outputs are: 1. Product count n f = (n q , n r 1 , n r 2 , n d ) where n q,r 1 ,r 2 ,d is either 0 or 1.
For atom-diatomic molecule scattering, the two possible reactions are represented by n r 1 and n r 2 .2. Final state s f = (v t , w(v , v t ), j) where w(v , v t ) is the Gaussian weight given in equation 11.For trajectories yielding n d = 1, the final state outputs s f = (0,0,0).

Final positions
Each trajectory is labeled by the input parameters (E c , b).

Parallel Implementation & Performance
The code is best used in a parallel implementation, which dramatically speeds up the time per trajectory as the number of CPUs is increased (Figure 5.The save short() and save long() methods from util.py file uses the Python package multiprocess to do this.As shown in the sim.py file, the user can choose to save a long output, where the program outputs a new line per trajectory containing all of the output data, or a short output, which creates a new line per (E c , b) by summing the product count vector n f .The long output is required for final state distribution analysis.Figure 5 shows the time per trajectory as the number of CPUs varies, averaged over 1000 trajectories.This was performed on Stony Brook's Seawulf cluster.It is clear that parallel processing yields a significant decrease in the total time.The majority of the calculation time is spent during the solution of the Hamilton's equations of motion.Note that for parallel processing, the multiprocess package is required.This package uses dill, which gives more flexibility in what can be serialized for parallel computation.The solve ivp API allows the user to control absolute and relative tolerances to control local error estimates.These can be controlled via the input file, and will have a large influence on the energy conservation and time per trajectory.It is highly recommended to study the effect that these tolerances have on a system before running large calculations.Each QCT object has a delta e attribute which yields the total change in energy over the trajectory.

Visualization
The file plotters.pycontains several methods for visualizing the results of a trajectory.Fig. 2 was obtained using the traj plt function, which requires a trajectory object as input.There is also a 3d plot generator traj 3d which provides a trace of the event.Finally, there is a method to generate an animation of the trajectory, traj gif.The usage of these plotters is shown in the example Jupyter notebook.

Example
As an example, we demonstrate the calculation of CaH formation rate as a result of the reaction H 2 + Ca → CaH + H.The input file inputs.json.All of the simulation and potential parameters are listed.The potential ranges and parameters for H 2 and CaH were obtained from [27,28] and [29,30], respectively.We include the parameters of the three potentials we might be interested to study this system in; Morse, Lennard-Jones ('lj'), and Buckingham ('buck').Here we choose the Morse potential to describe the interaction of both H 2 and CaH: In this example, we run 10 Listing 2: Parallel implementation of PyQCAMS with comments.The "utils.saveshort()" runs the quasiclassical trajectories in parallel and outputs the summed data to the specified "out file".
The output data now has ∼ 10000 × 20 lines, each corresponding to one trajectory.Each of the 20 impact parameters leads to a different opacity function for a reaction, P r (E c , b).The opacity function yields the scattering cross section at different collision energies E c : and the rate For this example, we repeated the code in Listing 2 over 10 different collision energies E c .For each E c , we sum over the count vector n to obtain the opacity function, cross section, and rate as described.The rates for CaH formation are shown in Figure 6, where it is noticed that higher collision energies gives k (cm 3 /s) Figure 6: Rate of reaction H 2 + Ca → CaH + H. 9 different collision energies were considered, and 10 4 trajectories were run at each energy.We use 20 evenly spaced impact parameters between 0 and 5 a 0 for each collision energy.Here, the H 2 molecule was initiated at v = 1, j = 0 and each pairwise interaction was defined by a Morse potential.
From here, we see that CaH is most likely formed at a collision energy E c = 30000 K.
rise to larger reaction rates, as expected in endothermic reactions.However, at very high collision energies such a trend changes due to the dominance of molecular dissociation processes.A more detailed study of this reaction will be published elsewhere.We can also calculate the distribution of states of CaH and H 2 , using the final state vector s. Figure 7 shows these distributions at four different collision energies.We can see that the higher vibrational states fill up for both CaH and H 2 as the collision energy is increased, as is typical in endothermic reactions.All details of these calculations can be found in a Jupyter notebook example file.

Conclusions
We have presented a Python quasi-classical atom-molecule scattering program, PyQCAMS.The PyQCAMS program aims to provide an easyto-use platform for calculating quasi-classical trajectories for atom-diatomic  molecule systems, including the three most relevant potentials for diatomic molecules: Morse, Buckingham and the generalized Lennard-Jones.We discussed the underlying theory behind the program, and the methods of the program as they pertain to the theory.As output, the user obtains the reaction probability per energy and impact parameter.Then, with this information is possible to calculate the cross section and energy-dependent rate constant.Its object-oriented approach allows the user to study different properties of a given trajectory.The plotter tools make it easy to visualize the results of a trajectory, making it ideal for a new researcher studying trajectories or for presenting the topic in a classroom setting.

Acknowledgments
The authors acknowledge the generous support of the Simons Foundation.

Figure 1 :
Figure1: (Color online.)Jacobi coordinates of a atom-molecule system.Here, the molecule is rotating in the x-y plane, with the angular momentum J along the z-axis and κ = ρ 1 × ẑ. η is the angle between κ and J, and φ and θ are defined as usual in spherical coordinates.

Figure 2 :
Figure 2: (Color online.)Different product outcomes for a given QCT calculation for H 2 + CaH at E c = 50000 K, b = 0 a 0 .r 12 represents the internuclear distance between each H, whereas r 32 and r 31 represent the internuclear distance between the colliding atom Ca (3) and each molecular atom (H(1) and H(2)).These trajectories show dissociation (a), reaction (b), and quenching (c).In panel (b), r 31 oscillates around a fixed distance, an indicator that a new molecule is formed.

Figure 3 :
Figure 3: The general structure of the PyQCAMS program.

Figure 5 :
Figure 5: Time per trajectory as a function of the number of CPUs used in the parallel calculation.These trajectories were run at a collision energy of 40000 K for H 2 + Ca reactions, with a relative tolerance of 10 −12 and absolute tolerance of 10 −11 .

Figure 7 :
Figure 7: (Color online.)Probability distribution of final vibrational states of CaH (a) and H 2 (b) as a result of the reaction H 2 + Ca, calculated by PyQCAMS.H 2 was initiated at v = 1, j = 0, with each pairwise interaction defined by a Morse potential.Each color represents a different collision energy E c , and the impact parameter of each trajectory is fixed to b = 0.