M-CHIPR: A Mathematica program for constructing multi-state coupled adiabatic potential energy functions in triatomic molecule using many body partitioning approach,☆☆

https://doi.org/10.1016/j.cpc.2022.108419Get rights and content

Abstract

The permutationally invariant combined hyperbolic inverse power representation (CHIPR) polynomial expansion in terms of hyperbolic secant basis proves to be an efficient model in describing the multi state coupled potential energy surface (MSCPES) and this fact was tested on the O2H+ triatomic ion and its diatomic sub units upon dissociation with an unmatchable accuracy of fitting in comparison with the previous works (below 1cm1 in the case of diatomic function and below 300cm1 in the case of triatomic function). In general, this model is capable of modeling A3 (H3, O3 etc.) and ABC (OCS, HCO+ etc.) triatomic system in addition to the A2B (O2H+) system. This paper presents a Mathematica notebook (M-CHIPR.nb file) to perform the fitting of di and triatomic potentials in a user friendly manner with proper explanations written along side wherever required in the .nb file supplied. The advantage of this code is that the computing clusters are not required to use it and can be run in normal pc with good processor and memory, apart from its numerical accuracy and the system independent consistency over the similar programs written in other languages.

Program summary

Program title: M-CHIPR

CPC Library link to program files: https://doi.org/10.17632/4mbydkzyts.1

Licensing provisions: GPLv3

Programming language: Wolfram Mathematica version12 or higher

Nature of problem: Finding an multi-variable Mathematica combined hyperbolic inverse power representation (M-CHIPR) function to describe a multi state coupled potential energy surface for a typical diatomic and triatomic molecule via extensive nonlinear regression method involving hundreds of optimizable variables with minimal fitting error with respect to the raw data points.

Solution method: Default Mathematica implementation of nonlinear least square fitting of ab initio points to obtain a well optimized model function having ability to accurately predicting the data behavior.

Additional comments including restrictions and unusual features: The entire work was tested on Apple Mac pc with Intel i7 processor having 4 cores. The version of Mathematica is 12 and therefore, the users are advised to test the supplied notebook with the same version and beyond to avoid incompatibility issues.

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 N(N1)/2 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 sech and tanh 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 O2H+ ion and its diatomic sub units. Earlier, the same function was used to fit the raw ab initio points of carbon clusters such as C3 and C4 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 O2 has the ability to form highly stable protonated O2 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.Vtotal=k=1Ni=1(Nk)Vi(k)[(Xk)i]

In Eq. (1), (Nk)=N!k!(Nk)!, 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, X={A1,A2,A

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)

  • S.N. Yurchenko et al.

    J. Mol. Spectrosc.

    (2007)
  • C. Rocha et al.

    Comput. Phys. Commun.

    (2020)
  • C. Rocha et al.

    Comput. Phys. Commun.

    (2021)
  • J.M. Robbe et al.

    Chem. Phys.

    (2000)
  • W.C. Ho et al.

    J. Mol. Spectrosc.

    (1991)
  • S.A. Nizkorodov et al.

    Chem. Phys. Lett.

    (1997)
  • D.G. Truhlar et al.

    Chem. Rev.

    (1987)
  • B.J. Braams et al.

    International Reviews in Physical Chemistry

    (2009)
  • S. Manzhos et al.

    Chem. Rev.

    (2021)
  • Z. Xie et al.

    Journal of Chemical Theory and Computation

    (2010)
  • C. Qu et al.

    Annu. Rev. Phys. Chem.

    (2018)
  • J. Li et al.

    J. Chem. Theory Comput.

    (2020)
  • S.L. Widicus Weaver et al.

    Astrophys. J.

    (2009)
  • F. George et al.

    Chem. Phys. Lett.

    (2018)
  • F.G.D. Xavier et al.

    J. Chem. Phys.

    (2010)
  • F.G.D. Xavier

    J. Phys. Chem. A

    (2010)
  • F. George et al.

    J. Phys. Chem.

    (2019)
  • H. Kohguchi et al.

    J. Chem. Phys.

    (2018)
  • 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).

    View full text