ORIGINAL ARTICLE: George D. Xavier, M-CHIPR: A Mathematica program for constructing multi-state coupled adiabatic potential energy functions in triatomic molecule using many body partitioning approach, Comput. Phys. Comm. 278 (2022) 108419. DOI: https://doi.org/10.1016/j.cpc.2022.108419
Article 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
O
2
+
H
triatomic ion and its diatomic sub units upon dissociation with an unmatchable accuracy of fitting in comparison with the previous works (below 1
-1
cm
in the case of diatomic function and below 300
-1
cm
in the case of triatomic function). In general, this model is capable of modeling
A
3
(
H
3
,
O
3
etc.)
and ABC (OCS,
+
HCO
etc.) triatomic system in addition to the
A
2
B
(
O
2
+
H
)
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.

Introduction

This note book contains Mathematica code to fit the diatomic and triatomic ab initio data and seam trajectory using the eigen value expression of 2 by 2 diabatic matrix
​
M-CHIPR is a Mathematica note book that can be downloaded and used for the above mentioned purpose by any one. The result from using this code appeared in few articles in the past. The user can consult them if they so wish.
​
The following sections describe the details about the model function and Mathematica command and how to use the code with proper syntax. The result is presented in the form of 2d and 3d plot. Comprehensive documentation and to run in Parallel can be found from Mathematica’s Documentation Center.

Diatomic script and its Application

User can write their directory path where the file containing the fitting data lies. Example:/Users/a/b/c
SetDirectory["write users directory"];
Importing the data file named “OH+.dat” here with Import command and saved in another variable “dataohp” and the points not required can be deleted using Delete command, consult online manual
dataohp=Delete[Import["OH+.dat"],{{1},{2},{3},{4}}]
x is the linear combination mentioned in the manuscript and the number of terms can be altered to improve the fitting quality
x=
a
1
Tanh
β
1
r1-
ζ
1
0
(
R
oh
)

0
(
R
oh
)
Sech
μ
1
r1-
ζ
1
0
(
R
oh
)
+
a
2
Tanh
β
1
r1-
ζ
1
1
(
R
oh
)

1
(
R
oh
)
Sech
μ
2
r1-
ζ
1
1
(
R
oh
)
+
a
3
Tanh
β
1
r1-
ζ
1
2
(
R
oh
)

2
(
R
oh
)
Sech
μ
3
r1-
ζ
1
2
(
R
oh
)
+
a
4
Tanh
β
1
r1-
ζ
1
3
(
R
oh
)

3
(
R
oh
)
Sech
μ
4
r1-
ζ
1
3
(
R
oh
)
;
model is polynomial of degree 7(can be of any degree)in x basis defined in the previous line and saved in the variable anamod1. polynomial degree can be altered to improve the fit
anamod1=
8
ω
r1
Sum
A
i
i
x
,{i,7};
The line Mathematica does the fitting with Nonlinearmodelfit command on the data stored in datohp using the model anamod1 to find the unknown coefficients with the running variable named r1, the value 1.86 is the starting value given to
R
oh
fitoh=NonlinearModelFit[dataohp,{anamod1},{{
R
oh
,(*1.97*)1.86},{ω},
β
1
,
a
1
,
a
2
,
a
3
,
a
4
,
A
1
,
A
2
,
A
3
,
A
4
,
A
5
,
A
6
,
A
7
,
ζ
1
,
μ
1
,
μ
2
,
μ
3
,
μ
4
},r1,MaxIterations32000];
The optimized coefficients after the fitting run
fitoh["BestFitParameters"]
calculate root mean square error
Print"RMSD=",RootMeanSquare[Table[dataohp[[i]][[2]]-fitoh[dataohp[[i]][[1]]],{i,Length[dataohp]}]]219474.63,"
-1
cm
"
The user can adapt the following script to make their plot which will look nearly the same after making necessary changes to the script and also consulting mathematical manual regarding how to plot
Show[Plot[{fitoh[r],fito2p[r],fito2[r]},{r,0,50},PlotRange{{1.2,8},{-0.26,0.05}},FrameTicks{{2,4,6,8},{-0.25,-0.20,-0.15,-0.10,-0.05,0,0.05}},Epilog{Text[Style["
+
O
2
",FontSize14],{6,0.1}]},Axes{False,True},FrameStyleDirective[16],LabelStyle{FontSize16},FrameLabel{"bond distance/bohr","energy/hartree"},Frame->True,PlotLegendsPlaced[{Style["
+
OH
fit",White,14],Style["
+
O
2
fit",White,14],Style["
O
2
fit",White,14]},{Scaled[{0,0.5}],{-4.5,0.5}}],PlotTheme"Marketing"(*,GridLinesNone*)],ListPlot[{dataohp,datao2p,datao2},PlotLegendsPlaced[{Style["
+
OH
data",White,14],Style["
+
O
2
data",White,14],Style["
O
2
data",White,14]},{Scaled[{0,0.5}],{-3.95,1.5}}],PlotTheme"Marketing"],ImageSizeLarge]
This plot is obtained by running the above script and rmsd obtained is below
1
-1
cm
and the length of the contracted basis is not same for all diatomic species. For OH+, there are 4 terms and polynomial of degree 7 and for other two species, there are six terms with polynomial of degree 10. The user can adapt the above script for fitting the molecular data of their choice

Triatomic script and its Application

Analytical model to fit the data of three body potential given in this script.​
​
y
1
,
y
2
,
y
3
are the linear combination for A2B system, each one with five terms, but it can be changed to any number. Duration of fitting depends upon this number of terms, more the terms more the quality and heavier the computation. A minimum of 3 term is required to obtain meaningful fit. For other two types of molecular system, to write the expression for
y
1
,
y
2
,
y
3
please look at the manuscript. The expansion of polynomial is written up to the degree of 15, each degree with the correct number of monomial and their associated coefficients applicable to A2B system. For each degree, list of monomial are closed in parenthesis indicating the degree in pink. There is no need for user to use full expansion to fit initially. They can start with smaller expansion first and if not satisfied with the result, big expansion may be considered, sometimes overriding limit set in the code 15 to 16,17 etc by extending the expansion. Lastly, Mathematica key word to perform nonlinear fit is given and saved the arbitrary function called “fit” which is again called in the sequence fit[u,v,w] with replacement of variable, r1--->u,r2---->v,r3----->w. For example, fit[2,3,4] gives the result.
y
1
=p11
Tanh
β
1
u-
ζ
1
0
(
R
ref1
)

0
(
R
ref1
)
Sech
γ
1
u-
ζ
1
0
(
R
ref1
)
+p12
Tanh
β
1
u-
ζ
1
1
(
R
ref1
)

1
(
R
ref1
)
Sech
γ
2
u-
ζ
1
1
(
R
ref1
)
+p13
Tanh
β
1
u-
ζ
1
2
(
R
ref1
)

2
(
R
ref1
)
Sech
γ
3
u-
ζ
1
2
(
R
ref1
)
+p14
Tanh
β
1
u-
ζ
1
3
(
R
ref1
)

3
(
R
ref1
)
Sech
γ
4
u-
ζ
1
3
(
R
ref1
)
+p15
Tanh
β
1
u-
ζ
1
4
(
R
ref1
)

4
(
R
ref1
)
Sech
γ
17
u-
ζ
1
4
(
R
ref1
)
;​​​​
y
2
=p21
Tanh
β
2
v-
ζ
2
0
(
R
ref2
)

0
(
R
ref2
)
Sech
γ
5
v-
ζ
2
0
(
R
ref2
)
+p22
Tanh
β
2
v-
ζ
2
1
(
R
ref2
)

1
(
R
ref2
)
Sech
γ
6
v-
ζ
2
1
(
R
ref2
)
+p23
Tanh
β
2
v-
ζ
2
2
(
R
ref2
)

2
(
R
ref2
)
Sech
γ
7
v-
ζ
2
2
(
R
ref2
)
+p24
Tanh
β
2
v-
ζ
2
3
(
R
ref2
)

3
(
R
ref2
)
Sech
γ
8
v-
ζ
2
3
(
R
ref2
)
+p25
Tanh
β
2
v-
ζ
2
4
(
R
ref2
)

4
(
R
ref2
)
Sech
γ
9
v-
ζ
2
4
(
R
ref2
)
;​​​​
y
3
=p21
Tanh
β
2
w-
ζ
2
0
(
R
ref2
)

0
(
R
ref2
)
Sech
γ
5
w-
ζ
2
0
(
R
ref2
)
+p22
Tanh
β
2
w-
ζ
2
1
(
R
ref2
)

1
(
R
ref2
)
Sech
γ
6
w-
ζ
2
1
(
R
ref2
)
+p23
Tanh
β
2
w-
ζ
2
2
(
R
ref2
)

2
(
R
ref2
)
Sech
γ
7
w-
ζ
2
2
(
R
ref2
)
+p24
Tanh
β
2
w-
ζ
2
3
(
R
ref2
)

3
(
R
ref2
)
Sech
γ
8
w-
ζ
2
3
(
R
ref2
)
+p25
Tanh
β
2
w-
ζ
2
4
(
R
ref2
)

4
(
R
ref2
)
Sech
γ
9
w-
ζ
2
4
(
R
ref2
)
;​​​​VR=(*2*)
Plus[
]
(*15*);​​​​Block{iter=0},fit=NonlinearModelFitTotint,
1
8
(1+Tanh[xu(u-Abs[u0])])(1+Tanh[xv(v-Abs[v0])])(1+Tanh[xw(w-Abs[w0])])VR+fito2p[u]+fitoh[v]+fitoh[w],(*2*)
,{u,v,w},(*Weightswts,*)MaxIterations1000,StepMonitoriter++,iter​​fit["BestFitParameters"]
A 2d cut of model function for various Jacobi angle as a function of
O
2
bond distance
A relaxed Plot of mode function by relaxing
O
2
as a function of distance between cm and
+
H
expressed in Cartesian coordinates

Two state cusp fitting

Hyperspherical surface of two states showing cusp and respective geometires

Diatomic Curves (Tested Examples)

Calculation of Three Body Potentials in Valence Coordinates/Conversion from Jacobi to Valence Coordinates

CITE THIS NOTEBOOK

M-CHIPR: A Mathematica program for constructing multi-state coupled adiabatic potential energy functions in triatomic molecule using many body partitioning approach​
by George Xavier​
Wolfram Community, STAFF PICKS, November 6, 2023
​https://community.wolfram.com/groups/-/m/t/3061322