M-CHIPR: A Mathematica program for constructing multi-state coupled adiabatic potential energy functions in triatomic molecule using many body partitioning approach☆,☆☆
Introduction
A spectroscopic accurate potential energy function describing all the usual features present in a (stationary point, saddle point, seam trajectory, dissociation limits/products, unified atom limit, optimal conical intersection geometry etc.) typical adiabatic potential energy surface of any bound molecular system is absolutely necessary to proceed further with the calculations of reaction dynamics to estimate the rate constants of any state specific reactive processes evolving over the single or multiple surfaces. Such functions (analytical model) often involve multiple variables with the total number of such variables varying as with N being the number of atoms. Therefore, depending upon the size of the molecular system as well as the overall size (length) of the model function, a single point evaluation time at a specific geometry varies considerably from being the least for the smaller (triatomic) systems to larger values for the poly atomic one involving more than 10 atoms. Over the last three decades, various groups had developed Fortran subroutines for the purpose of fitting the raw ab initio data to obtain a suitable model function and they have all been well documented in the literature [1], [2], [3], [4], [5], [6], [7]. The model and methodology used in the present work is almost similar to the one proposed earlier [2] which used morse function as the basis where else hyperbolic and function are used as basis in the present work.
The system we have tested using this basis (called collectively as CHIPR function standing for combined hyperbolic inverse power representation) is ion and its diatomic sub units. Earlier, the same function was used to fit the raw ab initio points of carbon clusters such as and and the Fortran code to achieve this task was published in this journal [8], [9]. Regarding this ion, it was extensively studied from computational point of view since eighties by many people [10], [11], [12], [13], [14], [15] and few experimental works [16], [17], [18], [19], [20], [21], [22] also appeared in the last few years. This ion upon the addition of one more neutral has the ability to form highly stable protonated dimer [23], [24], [25]. Before this, the weakly bound unprotonated dimer was also researched well [26], [27], [28], [29], [30], [31], [32], [33], [31].
In the next section, theory behind this model will be explained, followed by the explicit working details of the present Mathematica code with proper explanation in the next section. In the last section, the results obtained by running the code will be presented and its superiority over the Fortran code will be highlighted [8], [9]. Finally, concluding remarks will be given.
Section snippets
Analytical model
According to many body expansion theory, the total interaction potential energy for a N-atom molecule can be cast in the form of a generalized summation using the theory of combinatronics as given below.
In Eq. (1), , a combination (permutationally invariant) counting formula, gives the number of elements in k-subset of length k. It is a positive integer greater than zero. X denotes the set of complete list of atoms in the N-atom molecule,
Input details
The code can accept the input data either from an external file mentioned at the top of the small script given in the subsection 3.2 or from the manually typed list of ab initio data in the notebook cell after saving it in a variable to be used later and will look like in the format, {{1,-0.01},{2,-0.2},{3,-0.3},{4,-0.25}⋯{10,-0.001}}. Typing is practical only when the number of data is small enough to do. In the case of diatom, this input must be a set of pair of numbers {x,V} where x is
Conclusion
In this work, a working Mathematica code is presented with the objective of fitting the diatomic and triatomic ab initio data and obtaining a suitable potential energy function capable of predicting accurate molecular properties in a user friendly way. The quality of the fit is always superior when compared with the previous works [8], [9] and the accuracy is well below 1 cm−1 in the case of potential energy function and below 300 cm−1 in the case of triatomic surface.(see the text) The model
Declaration of Competing Interest
George D Xavier reports financial support was provided by University of Coimbra Faculty of Sciences and Technology. George D Xavier reports a relationship with University of Coimbra Faculty of Sciences and Technology that includes: employment.
References (36)
- et al.
J. Mol. Spectrosc.
(2007) - et al.
Comput. Phys. Commun.
(2020) - et al.
Comput. Phys. Commun.
(2021) - et al.
Chem. Phys.
(2000) - et al.
J. Mol. Spectrosc.
(1991) - et al.
Chem. Phys. Lett.
(1997) - et al.
Chem. Rev.
(1987) - et al.
International Reviews in Physical Chemistry
(2009) - et al.
Chem. Rev.
(2021) - et al.
Journal of Chemical Theory and Computation
(2010)
Annu. Rev. Phys. Chem.
J. Chem. Theory Comput.
Astrophys. J.
Chem. Phys. Lett.
J. Chem. Phys.
J. Phys. Chem. A
J. Phys. Chem.
J. Chem. Phys.
Cited by (1)
A global CHIPR potential energy surface of PH<inf>2</inf>(X<sup>2</sup>B<inf>1</inf>) via extrapolation to the complete basis set limit and the dynamics of P(<sup>2</sup>D) + H<inf>2</inf>(X<sup>1</sup>Σ<sup>+</sup><inf>g</inf>) → PH(X<sup>3</sup>Σ<sup>−</sup>) + H(<sup>2</sup>S)
2022, Physical Chemistry Chemical Physics
- ☆
The review of this paper was arranged by Prof. Stephan Fritzsche.
- ☆☆
This paper and its associated computer program are available via the Computer Physics Communications homepage on ScienceDirect (http://www.sciencedirect.com/science/journal/00104655).