ORIGINAL ARTICLE

Field of the paper

Numerical Solution of Schrodinger Equation for

Morse Potential using Matrix Methods with
Fourier Sine Basis

Aditi Sharma and O.S.K.S. Sastri1

1 Department of Physical and Astronomical
Sciences, Central University of Himachal In this paper, an elegant and easy to implement numerical
Pradesh(CUHP), Dharamshala, H.P.,176215, method using matrix mechanics approach is proposed, to
India(Bharat)
solve the time independent Schrodinger equation (TISE) for
Correspondence Morse potential. It is specifically applied to non-homogeneous
O.S.K.S. Sastri, Department of Physical and
Astronomical Sciences, Central University diatomic molecule HCl to obtain its rotating-vibrator spec-
of Himachal Pradesh(CUHP), Dharamshala, trum. While matrix diagonalization technique is utilised for
H.P.,176215, India(Bharat)
Email: [email protected] solving TISE, model parameters for Morse potential are op-
timized using variational Monte-Carlo (VMC) approach by
Funding information
minimizing χ 2 − value. Thus, validation with experimental
vibrational frequencies is completely numerical based with
no recourse to analytical solutions. The ro-vibrational spec-
tra of HCl molecule obtained using the optimized parame-
ters through VMC have resulted in least χ 2 − value as com-
pared to those determined using best parameters from mul-
tiple regression analysis of analytical expressions. Numer-
ical algorithm for solving the Hamiltonian matrix has been
implemented utilizing Free Open Source Software (FOSS)
Scilab and simulation results are matching well with those
obtained using analytical solutions from Nikiforov-Uvarov
(NU) method and asymptotic iteration method (AIM).
Keywords — Morse potential, Matrix method, Variational
Monte Carlo (VMC), Scilab.

Abbreviations: ABC, a black cat; DEF, doesn’t ever fret; GHI, goes home immediately.

1
2 Aditi and O.S.K.S. Sastri

1 | INTRODUCTION

Obtaining the bound state solutions of time independent Schrodinger equation (TISE) has remained one of the cen-
tral aspects of quantum mechanics and its application to atomic, molecular, nuclear and solid state physics. All these
physical systems are modeled using various potentials that typically have an attractive nature that dies down zero
reasonably fast and a repulsive core at short distances. A large number of exponential type potentials such as Rosen-
Morse((Zhang2012), (Ikhdair 2010)), Manning-Rosen ((Manning1933), (Manning1933), (Qiang2007), (Ikhdair2008)),
Kratzer(Berkdemir2006), (Berkdemir2007), (Ikhdair2008)) and other multi-parameter potentials ((Egrifes2000), (Jia2004))
are used in both molecular and nuclear physics. Among all these, Morse potential has proved to be one of the most
successful in explaining rotational and vibrational structure of diatomic molecules (Rosen1932), even till date.
Many analytical methods have been proposed to solve TISE for obtaining its exact energies and corresponding
wave-functions. The recent analytical techniques like super-symmetry quantum mechanics (SUSYQM) ((Levai1992),
(Morales2004)), Nikiforov-Uvarov ((Ikhdair1992), (Nikiforov2004), (Berkdemir2005)), point canonical transformations
(De1992), Pekeris-type approximation (Zhang2011), Taylor expansion(Chou2011) Lie algebraic method(Palma2011)
and Asymptotic Iteration Method(AIM) have been very successful in solving TISE for complex potential functions
which were beyond the reach of regular solutions using Frobenius method or special functions approach. The eigen
energies of Morse potential with the use of Perkeris approximation are obtained for Morse potential using AIM(Al-
Dossary 2007) and SUSYQM (Hassanabadi 2012). Also, bound state solutions of D-dimensional Schrodinger equation
are determined for class of exponential type Hulthen, Manning-Rosen, Eckart and Deng Fan potentials(Pena 2014).
Likewise, various numerical methods have been implemented to solve second order differential equation which con-
sists of 1/N expansion ((Qiang2007a),(Ikhdair2003)), finite element analysis (Xu2009), Runge–Kutta techniques ((Van
de Vyver 2005),(Anastassi 2005)), homotopy analysis((Alomari 2009),(Jia 2010)) and so on. Recently, Hassanabadi et.
al. followed both SUSYQM and finite difference methods to obtain reliable solution of radial Schrodinger equation for
Hua potential which is similar to Eckart, Manning–Rosen, or Morse potentials.
In this paper, TISE is solved for obtaining ro-vibrational spectrum of Morse potential using matrix method tech-
nique proposed by Marsglio et.al,. (Marsiglio 2009). They first implemented it for harmonic oscillator potential(Marsiglio
2009) and then extended it further for central potentials such as spherical well, Coulomb and Yukawa potentials (Jug-
dutt 2013). The central idea of the method is to embed the potential of interest within an infinite square potential,
whose solutions which are sine eigen-functions, are utilised as basis for obtaining the matrix form of TISE to be solved
using an eigen solver. Following this, we have recently applied it to obtain single particle neutron and proton energies
using nuclear shell model with Woods-Saxon and spin-orbit interactions (Sharma 2020). Also, a study of square well
potential (Sastri 2019), pure quartic and anharmonic oscillators (Sharma 2020a), and vibrational spectrum of diatomic
molecules (Sastri 2020) have been undertaken in simple worksheet environment in Gnumeric software. Here, we ex-
tend the study to obtaining the rotational fine structure of the vibrational spectrum. One of the important features of
this study involves optimising the model parameters using VMC technique. The χ 2 -value between the simulated and
experimental frequencies (which is akin to least squares minimisation) is reduced in a variational sense, by randomly
allowing each of the parameters to be varied successively in every iteration, till the required accuracy is achieved. The
paper is structured as follows.
First, in section 2, the need for modeling the interaction by Morse type function has been reinforced by obtaining
the electronic ground state energy for different distances from first principle calculations. This is followed by section
3 that deals with preparation of system for simulation, by briefly discussing the numerical technique being used and
rephrasing the problem in appropriate units. Then, the steps for implementation are clearly specified. Following which,
testing of code and optimising of algorithm parameters, for a known set of model parameters obtained from fitting
Aditi and O.S.K.S. Sastri 3

experimental data using analytical expressions, is presented. In section 4, the procedure for optimisation of model
parameters using VMC approach is deliberated and the results are generated in the process. The simulation results for
vibrational spectrum consisting of its rotational fine structure have been validated with experimental values obtained
from FTIR (Polik 1999). Finally, the optimized model parameters using VMC have been substituted into analytical
expressions for energies from NU and AIM methods, to obtain the rotational energies within each vibrational level. A
comparision with simulation results, shows perfect match between all the three approaches.

2 | MODELING DIATOMIC MOLECULE AS A VIBRATING ROTATOR:

In this section, we first discuss aspects of modeling such as description, formulation and ramification stages involved
in solving for vibrational and rotational energies of a diatomic molecule.

2.1 | Description

The rotational and vibrational motions of a linear non-homogeneous diatomic molecule of AB-type consisting of nuclei
N A & N B , and n electrons can be described by considering three types of Coulomb interactions:

1. e− -e− repulsion:

e2
n−1 Õ
Õ n
Ve − −e − = (1)
4πo r i j
i =1 j =i +1

where r i j = | r¯i − r¯j | is the distance between it h e− and jt h e− .

2. NA -NB repulsion:

ZA ZB e 2
VN −N = (2)
4πo R AB

where R AB = | R̄ A − R̄ B | is distance between the nuclei, ZA and ZB are their respective charges.
3. e− -N attraction:

−Z A e 2 −Z B e 2
N 
Õ 
Ve − −N = + (3)
i =1
4πo | r¯i − R̄ A | 4πo | r¯i − R̄ B |

and is dependent on the distance between it h e− and corresponding nuclei.

2.2 | Formulation:

The Hamiltonian for the diatomic molecule is given by

~2 ∂ 2 ~2 ∂ 2 ~2 Õ ∂ 2
! N
!
H = − − + − + VN −N + Ve − −N + Ve − −e − (4)
2m A ∂ R̄ 2 2m B ∂ R̄ 2 2m e ∂ r¯i2
A B i =1
4 Aditi and O.S.K.S. Sastri

Using Born Oppenheimer (BO) approximation, wave-function Ψ(r¯i , R̄ A , R̄ B ) as product of e− and nuclear wave-
functions φe (r¯i ; R̄ A , R̄ b ) and χN (R̄ A , R̄ b ) respectively, and using separation of variable, one obtains the TISE for e − s
and nuclear vibrational and rotational motion as

~2 Õ ∂ 2 φe
N
− + [Ve − N + Ve − e − + VN N (R AB )]φe = E e φe (5)
2m e
i =1
∂ r¯i2

and

1 ∂2 1 ∂2
!
− − + U e (R AB ) χN (R̄ A , R̄ B ) = E χN (R̄ A , R̄ B ) (6)
2m A ∂ R̄ 2 2m B ∂ R̄ 2
A B

where U e (R AB ) is the potential function obtained by solving electronic TISE for its ground state energy E e (r¯i ; R̄ A , R̄ B )
for different inter-nuclear distances R AB . When TISE is transformed, using (R̄ A , R̄ B ) → (R̄ C oM , R̄ ), into a one-body
problem describing the reduced mass µ of the molecule, R AB becomes the relative distance co-ordinate R = | R̄ | w.r.t
center of mass co-ordinate R̄ C oM .
Using FOSS such as Gamess, a code for obtaining ground state energy of the molecule using Hartree-Fock (HF) or
density functional theory (DFT), typical functional form of U e (R ) is obtained (see Fig.1) and it has a form as proposed
by Morse (Morse 1932), given by

 
U e (R ) = D e e −2α x − 2e −α x (7)

(R −R e )
Here x = Re and α = bR e , where R = R e is the equilibrium bond length at which the potential attains a minima.
Parameter α describes the shape of the curve and is a characteristic that is different for various molecules. With
increasing R , the potential reaches a constant value D e , as a result of dissociation of the molecule. Since the potential
is dependent only on distance R , it is central in character and spherical polar co-ordinates are appropriate choice for
reference system.
Once again invoking the BO approximation, in this case, to separate out the vibrational and rotational wave-functions,
one can obtain the radial equation describing the vibrations as (Struve 1989)

d2 S (R ) 2µ h J (J + 1)~2 i
+ 2 E − U e (R ) − S (R ) = 0 (8)
dR 2 ~ 2µR 2

J (J +1)~2
So, an effective vibrational potential can be defined as U (R ) = U e (R ) + U cf (R ), where U cf (R ) = , represents
2µR 2
the centrifugal potential that gives rise to the rotational levels of the molecule.
Aditi and O.S.K.S. Sastri 5

F I G U R E 1 Graph of inter-nuclear distance R with ground state energy U e (R ) of HCl molecule, obtained using
FOSS Gamess

2.3 | Ramification:

The solution of TISE, as given by Eq. (8), can be obtained either analytically or numerically. The analytical solution is
available for only studying vibrational spectra and energy values of vibrator-rotator is given by (Banwell 1972)

1 1
ε(v , J ) = ω̄e (v + ) − ω̄e x e (v + )2 + B J (J + 1)cm −1 (9)
2 2

where ω̄e is oscillation frequency for anharmonic system at equilibrium and ω̄e x e is anharmonicity constant which is
always smaller than ω¯e and is positive. B is rotational constant given by

h
B= cm−1 . (10)
8π 2 µR e2 c

Emergent Properties:
The selection rules for transitions between the energy levels are obtained by utilising the eigen-functions of vibrator-
rotator, which are the product of eigen-functions of vibrator and rotator respectively. Therefore, selection rules will
be same as those for these two systems individually,(Struve 1989) i.e.

∆v = ±1, ±2, ±3, ...

∆J = ±1 (11)

Vibrational Spectrum:
The analytical expressions for the first three vibrational frequencies corresponding to absorption from v 00 = 0 and
their experimental values for HCl, are given by (Banwell 1972)

ν̄0→1 = εv =1 − εv =0 = ω¯e (1 − 2x e )cm −1 = 2886cm −1 (12)

6 Aditi and O.S.K.S. Sastri

ν̄0→2 = εv =2 − εv =0 = 2ω¯e (1 − 3x e )cm −1 = 5668cm −1 (13)

ν̄0→3 = εv =3 − εv =0 = 3ω¯e (1 − 4x e )cm −1 = 8347cm −1 (14)

Vibrating Rotor Spectrum:

Now, including the rotational fine structure, for the main absorbtion, that is, v 00 = 0 −→ v 0 = 1, the wavenumbers
corresponding to this transition are given by

ν̄ = ω¯0 + B 1 J 0 (J 0 + 1) − B 0 J 00 (J 00 + 1)cm −1 (15)

where ω¯0 = ω¯e (1 − 2x e )cm −1 is fundamental frequency and B 0 & B 1 are the rotational constants corresponding to
vibrational states v = 0 and v = 1 respectively, with B 0 > B 1 . The ν̄ can be determined from different values of J 0
and J 00 . So, the rotational transitions ∆J = J 0 − J 00 = +1 gives one set of lines called as R branch, whereas those
corresponding to ∆J = J 0 − J 00 = −1 are known as P branch.

3 | NUMERICAL SOLUTION OF VIBRATING ROTATOR

The main goal of this paper is to solve TISE given by Eq. (8) by numerical approach using matrix methods to obtain
the energy eigen values and then apply variational Monte-Carlo (VMC) approach to optimize the model parameters
such that frequencies obtained from simulation match the experimental ones by minimizing the χ 2 -value.

3.1 | Preparation of system for numerical solution

This stage involves three steps:

1. Choice of numerical technique: Schrodinger equation is often solved using Runge-kutta method, central divided
difference or Numerov methods. In this paper, we focus on matrix diagonalization technique suggested by Marsiglio
((Dauphinee T. 2015),(Felipe 2016)) which is stable, highly efficient and gives accurate results as compared with the
mentioned methods.
Matrix method: In this approach, the potential of interest Uv i b (R ) is embedded within an infinite square well potential
[0, a 0 ]. a 0 is the width, chosen to be certain cut-off radius, such that the characteristics of the potential are well
represented. In this case, a 0 can be chosen as the value of R beyond which the potential value saturates to D e and
does not change any further significantly. This is equivalent to writing Eq. (8) as

h −~2 d2 S (R ) i
+ Vi nf (R ) + U T (R ) S (R ) = E S (R ) (16)
2µ dR 2

which can be rewritten as

h i
H 0 + U T (R ) S (R ) = E S (R ) (17)
Aditi and O.S.K.S. Sastri 7

where H 0 is the Hamiltonian of the infinite square well potential and U T (R ) is the truncated potential defined in [0,a 0 ].
Expanding radial wave-function S (R ) in terms of infinite square eigen functions, chosen as basis functions

s
∞
2
 
Õ nπR
S (R ) = cn si n (18)
a0 a0
n=1

Eq. (17), can be transformed into the following matrix equation

∞
Õ
H mn c n = E c m for m = 1,2,3,. . . (19)
n=1

where H mn = E n0 δ mn + Vmn is an infinite dimensional square symmetric matrix.

E n0 are the energy eigen values of infinite square well Hamiltonian H 0 given by

n 2 π 2 ~2
E n0 = = n 2 E 10 (20)
2µa 02

and Vmn is given by

s s
a0 2 2
∫    
mπR nπR
Vmn = si n T
U (R ) si n dR (21)
0 a0 a0 a0 a0

We need to choose a finite number of basis functions, say N 0 for solving the problem numerically. The parameters a 0
and N 0 introduced into the technique are related to the algorithm, which are adjusted to ensure convergence and to
match with expected results from experiment.
2. Rephrasing of units: In atomic and molecular physics, energy is measured in eV, distances in Å and frequencies in
cm−1 . Therefore, these are chosen as the units for numerical calculations and the centrifugal and ground state energy
terms are rephrased in eV as follows:

J (J + 1)~2 J (J + 1)(~c)2
U c.f (R ) = 2
= eV (22)
2µR 2(µc 2 )R 2 × 106

where µ is reduced mass of diatomic molecule expressed in eV /c 2 and ~c = 1973.29eV Å. and

π 2 ~2 π 2 (~c)2
E 10 = 2
= eV (23)
2µa 0 2(µc 2 )a 02 × 106

3. Discretizing the continuous variable: For doing simulation in computers, it is essential to discretize the continuous
variable R with appropriate step-size and restrict it to a finite region of interest [0,a 0 ].

3.2 | Implementation of Matrix Methods Algorithm:

1. Initialisation: The model parameters of Morse potential for HCl molecule to be initialised. In the NIST data
compiled by Huber & Herzberg as well as in other references (Huber 1972), instead of Morse potential parameter
values, the spectroscopic information of various diatomic molecules is given which includes constants like ω̄e ,
8 Aditi and O.S.K.S. Sastri

ω̄e x e , R e . The formula to calculate D e (Bernath 2020) and b (Zhang 2012) model values from data available is

ω̄e2 ω̄e2
De = cm −1 = ∗ 1.2397 × 10−4 eV (24)
4ω̄e x e 4ω̄e x e

and
s s
ke µ ω̄e2 −1
b= = Å (25)
2D e 2D e

where R e can be calculated by using expression (Bernath 2020)

s
h
Re = Å (26)
8π 2 µcB e

parameters from multiple regression analysis of analytical expressions. Numerical algorithm for solving the Hamil-
tonian matrix has been implemented utilizing Free Open Source Software(FOSS) Scilab and simulation results are
matching well with those obtained using analytical solutions from Nikiforov-Uvarov (NU) method and asymptotic
iteration method (AIM).
where B e is rotational constant corresponding to the separation R e at the minimum of potential curve expressed
in cm −1 .
2. Potential definition: The anharmonic Morse potential is chosen to be embedded within an infinite square well
potential of width a 0 = 6 (where the potential saturates to constant value), and is given by

J (J + 1)~2 c 2
UvTi b (R ) = D e (e −2α x − 2e −α x ) + (27)
2µc 2 R 2

3. Defining H mn matrix: Hamiltonian matrix defined in Eq. (19) will be a square symmetric matrix consisting of
diagonal and non-diagonal elements. Even though integrals can be solved analytically for J = 0 case, for the
more general J , 0 that results in the rotational fine-structure require the matrix elements to be obtained using
numerical integration. The diagonal and non-diagonal elements of H nm matrix are obtained as

1−cos(2nπR /a 0 )
 
 E n0 + 0 0 a2
∫a
UvTi b (R )dR , f or m = n

2

0
= ∫

H nm (28)
 a 0 2 (1−cos((n−m)πR )/a 0 )−(1−cos((n+m)πR /a 0 ) U T (R )dR ,
 
 0 a0

2 vi b f or m,n

For doing numerical integration, ‘intspline’ function is used.

4. Obtaining eigen values: In Scilab, ‘spec’ command is used for getting eigen values but the generated eigen values
are in random order. So, ‘gsort’ is used for sorting the data in descending order. The obtained energies (in eV)
are converted into wavenumbers by dividing all of them them with a factor hc = 1239.8419 ∗ 10−7 eV cm. The
vibrational frequencies can be obtained by taking appropriate energy differences as per selection rules.

The code for obtaining the vibrational frequencies from the energies obtained as eigen values of the h-matrix has
been written in Scilab and is given in Appendix 1.
Aditi and O.S.K.S. Sastri 9

3.3 | Testing the correctness of code

Using the data obtained through multiple regression analysis,(Struve 1989) and the formulae given in Eqs. (24) to (26),
the values of D e , b and R e are obtained as 5.346799eV , 1.734985Å−1 and 1.274808Å respectively. On running the code,
with a 0 = 6 and N 0 = 50, the first three pure vibrational frequencies (J = 0) corresponding to transitions from v 00 = 0 to
v 0 = 1, 2, 3 are obtained as 2911, 5824, 8730 in cm −1 . These are not very close to the experimental values as specified
in Eqs. (12)-(14). So, the number of basis functions, N 0 are increased in steps of 20 till 90, beyond which the resulting
frequencies are, ω̄10 = 2886, ω̄20 = 5668 & ω̄30 = 8346 in cm −1 , converging to experimental values. This confirms that
the code is implemented correctly to solve the TISE using matrix methods approach. Here, the model parameters
for Morse potential are chosen based on fitting experimental data with analytical expressions obtained from solving
the theoretical model. This is the usual procedure followed for performing numerical simulations. But, in this paper,
the main intent is to establish numerical simulation as an independent methodology for solving the theoretical model
without recourse to analytical solutions. Hence, we introduce in the next section, VMC optimization procedure for
obtaining the model parameters that best validate the experimental observations.

3.4 | Optimization of Morse potential parameters using variational Monte Carlo approach

While Monte-Carlo approach involves varying the model parameters randomly by a small amount in each iteration,
variational principle is invoked to minimize the χ 2 − value. So, over all VMC approach is a parallel methodology to least
square minimisation procedure that is often employed with analytical expressions. The implementation algorithm for
VMC is as follows:
1. Initialization: To begin with, the model parameters- D e , b and R e could be chosen based on some guess. Alterna-
tively, the values could be chosen randomly or based on inputs from other techniques.
Here, for HCl molecule, equilibrium bond length R e = 1.2745 known from rotational spectra is taken as its initial value.
The parameters D e = 5 and b = 1 are initial guesses, which are chosen some what close to data available from litera-
ture D e = 5.3 and b = 1.7439.
2. Determination of simulation frequencies: Using these initial values, the TISE is solved using the matrix methods
technique elucidated in the previous section and the energy eigen values are obtained. Since, each of the frequencies
are a result of transitions to different energy levels, they can be regarded as independent degrees of freedom.
χ 2 − value: Hence, χ 2 can be calculated using following formula

3 (ν si m − ν ex p )2
χ2 =
Õ
i i
(29)
3
i =1

This is akin to the idea of mean-square error. Store the value in a variable named chi sq r ol d . 3. Monte-Carlo step:
Generate a random number ‘r 0 within an interval, say [-I,I]. Choose one parameter among the three, say D e and change
it to a new value by adding r to it. That is, say, D e new = D e + r . Redetermine the simulated frequencies using this
D e new and obtain χ 2 − value. Store it in a new variable chi sq r new .
4. Variational Step: If chisqrnew is less than chisqrold, then the new value is accepted. That is D e is updated as D e new ,
else old value of D e is retained. That is nothing is done. If the condition is satisfied then one has to update the value
of χ 2 as well. That is, chi sq r ol d = chi sq r new .
5. Iterative Minimization: When the above steps are repeated for all parameters involved, then one iteration is
completed. Then, the process is repeated over many iterations till the change in χ 2 is not remarkable, at which instance,
size of interval is decreased by decreasing value of I. This allows for even smaller variations in the parameters. In
10 Aditi and O.S.K.S. Sastri

this way, one can reduce χ 2 till one achieves the desired level of accuracy to the appropriate decimal place in the
parameters involved.
The VMC code, for minimization of χ 2 value for an interval [-I,I], has been written in FOSS Scilab. It consists of inline
documentation that makes it easy to understand and is available with authors.

4 | SIMULATION AND DISCUSSION OF RESULTS

The simulation results, obtained using model parameters obtained from multiple regression analysis data during the
testing phase of the code, give a χ 2 -value of 0.245121(Polik 1999).
One of the major contributions in this paper is to solve the TISE for Morse potential numerically without recourse to
analytical solutions. Towards this goal, we now optimize the model parameters for the Morse potential to yield the
best convergence to the experimental vibrational frequencies using the VMC approach. This has been implemented
as three step process.
Step 1 Experimentally, from microwave spectra of molecules, the inter-nuclear distance, R e is available to an accuracy
of five decimal places as 1.27455Å(Huber 2013) and is kept constant initially. Now, one can begin with any random
values of D e and b as initial guesses but there is a possibility that it could lead the process either to get into a local
minima or may result in a set of parameters that may not be those which are close to expected values obtained from
analytical solutions. Also, starting with completely random guesses can increase the number of iterations to extremely
large values thus increasing the computational time. It is always important to make calculated guesses from graphical
procedures.
From observing Fig. 1, we have chosen the starting values for D e and b to be 5 and 1 respectively. Now, the first three
vibrational transition frequencies using matrix method are determined by converting the energy values obtained in
eV for non-rotating Morse potential (J = 0) to wavenumbers (by dividing them with a factor of hc = 1239.84193 ∗
10−7 eV cm) and χ 2 − value is calculated with respect to experimental data available. This is stored in variable, named
as chisqrold.
In each iteration, one by one, the parameters D e and R e are changed by a random value chosen in an interval [-0.1,
0.1]. Each time the program is run to determine the new frequencies and the corresponding χ 2 − value. This is stored
in chisrqnew. The values of the parameters are updated if chisqrnew is less than chisqrold. This is repeated for 20,000
iterations. Using the obtained values of D e and b as inputs, that is, as the new starting values, the program is run again
for another 10,000 iterations. But, this time, the random value for making the change to the parameters is chosen
from a smaller interval of [-0.01, 0.01]. This process of reducing the interval limits by a decimal place and running the
program for 10,000 iterations is continued till the χ 2 − values has reached a value to less than 0.02. The results are
shown in Table 1.

Step 2 Just to ensure that the variations in R e are also accounted for, optimization procedure is repeated as before
by allowing for changes in all three parameters in each iteration. That is, by starting from the last obtained parameters
in the previous run, that is, D e = 5.347819eV , b = 1.734941Å−1 , and R e = 1.27455Å. Each parameter is made to change
by a random value picked in an interval [-0.00001,0.00001], (same as in the last run of previous step), in every iteration.
At the end of 30,000 iterations, the χ 2 −value reduced to 0.019219, a change in the fourth decimal place, resulting in
model parameter values as D e = 5.348363eV , b = 1.734841Å−1 , and R e = 1.274187Å. This completes obtaining one
single set of parameters using VMC technique.
Step 3 By proceeding in a similar fashion as discussed in previous two steps, four more sets have been completed
and model parameters are determined for each set with their corresponding chi-square values which are displayed in
Aditi and O.S.K.S. Sastri 11

TA B L E 1 Determination of parameters D e (eV ) and b(Å−1 ) keeping R e (= 1.27455Å) as constant, for reducing
interval sizes

Interval No. of iterations Old parameters New parameters χ2

De b De b

(-0.1,0.1) 20000 5 1 5.346809 1.735122 0.022031

(-0.01,0.01) 10000 5.346809 1.735122 5.346446 1.735122 0.079440

(-0.001,0.001) 10000 5.346446 1.735122 5.346853 1.735120 0.021558

(-0.0001,0.0001) 10000 5.346853 1.735120 5.346953 1.735102 0.021375

(-0.00001,0.00001) 10000 5.346953 1.735102 5.347819 1.734941 0.019944

TA B L E 2 Final optimized Morse potential parameters D e (eV ), b(Å−1 ) and R e (Å) obtained using VMC for five
different runs:

Sr. No. De b Re χ 2 -value

1 5.348363 1.734841 1.274187 0.019219

2 5.353239 1.733939 1.264900 0.018536

3 5.353913 1.733814 1.274379 0.019265

4 5.353953 1.733807 1.274438 0.019315

5 5.352160 1.734139 1.274553 0.017783

Table 2. Even though, one could take mean of the obtained values and determine the uncertainty, it is best to consider
the parameter values as those for which minimum χ 2 −value has been obtained.
When number of decimal places in final parameter values obtained was rounded to five, four and three places, the
χ 2 -value is obtained as 0.0178044, 0.0281803 and 0.319032 respectively. Hence, six decimal places are retained for
performing further calculations.

4.1 | Ro-vibrational spectrum of HCl

To validate the goodness of fit of model parameters obtained using VMC, the study of rotation-vibration spectra of HCl
molecule has been undertaken for obtaining the rotational fine structure in the fundamental transition. To accomplish
this, effect of rotation due to presence of centrifugal term needs to be included in solving the TISE in Eq. (8). Even
though, for each J -value, one will obtain many vibrational states, only the first two energy values corresponding to
v = 0 and v = 1 are recorded as shown in Table 3. Using selection rule ∆J = J 0 − J 00 = ±1, differences in energies
are accordingly taken to obtain R- and P- branch lines respectively as shown in columns 4 and 5 of Table 3. The same
lines are also obtained using analytical expressions as discussed below.
Analytical treatment: For ∆J = J 0 − J 00 = +1 in Eq. (15), we obtain R- branch lines with wavenumbers given by

ν̄R = ω¯0 + (B 1 + B 0 )(J 00 + 1) + (B 1 − B 0 )(J 00 + 1)2 (30)

12 Aditi and O.S.K.S. Sastri

where lower rotational quantum number, J 00 = 0, 1, 2, 3, ...

and for getting P-branch lines, substitute ∆J = J 0 − J 00 = −1 in Eq. (15) with wavenumbers

ν̄R = ω¯0 − (B 1 + B 0 )(J 0 + 1) + (B 1 − B 0 )(J 0 + 1)2 (31)

where J 0 = 0, 1, 2, 3, .... Also, rotational constant, B corresponding to the vibrational state v is given by

Bv = B e − β e (v + 1/2) (32)

where β e reflects how rotational constant B varies with vibrational quantum number v . So, for v = 0 state,

1
B0 = Be − βe
2

and for v = 1 state

3
B1 = Be − βe (33)
2

Using the values of β e = 0.30167 cm −1 , B e = 10.58919 cm −1 from multiple regression data[? ] and substituting in
above equations, B 0 and B 1 are obtained as 10.438355 cm −1 and 10.136685 cm −1 respectively. Now, putting these
values in Eqs. (30) and (31), wavenumbers corresponding to R-branch and P-branch lines are calculated which are
shown in column 6 and 7 of Table 3 respectively. In the final two columns, the experimental data recorded using a
Nicolet 730 FTIR spectrometer is shown. Since, the experimental results are available to 3 decimal places, accordingly
the simulation and analytical outcomes are rounded to 3 decimal places as well.

From Table 3, one can observe that results obtained through simulation by optimizing the potential parameters
through VMC are better than those obtained from analytical technique using multiple regression data, as these values
are much closer to experimental ones.

4.2 | Comparative analysis of Matrix method with NU and AIM methods:

While experimentally only frequencies are available, analytical and numerical solutions result in the bound state energy
levels, from which one obtains the transition frequencies. So as to check the effectiveness of matrix methods based
numerical technique that has been implemented in this paper, a comparative study has been undertaken with two
analytical techniques. The rotational energy levels for each of the first three vibrational levels are obtained by choosing
the optimized Morse potential model parameters using VMC that resulted in minimum χ 2 − value, as given in Table
2. The energy eigen values so obtained are compared with analytical results determined explicitly from approximate
techniques of Nikiforov-Uvarov(NU) and Asymptotic Iteration method(AIM).
The analytical expression to determine energy values through NU method (Berkdemir 2005) is given as

J (J + 1)~2 3 3 ~2 a 2 h 2 1 i2
E nJ = 2
(1 − + 2
)− √ − (n + ) (34)
2µR e 2
a0 a R e 2µ 2 3 2
Aditi and O.S.K.S. Sastri 13

TA B L E 3 Comparative validation of simulated P- and R-branch lines corresponding to fundamental vibrational

band in fine-spectrum of HCl molecule with those obtained using analytical and experimental techniques

J 00  J 00,v 00 =0  J 00,v 00 =1 Simulation(MM) Analytical(Regression) Experiment(FTIR)

(Polik (Polik
1999) 1999)

(cm −1 ) (cm −1 ) ν̄R (cm −1 ) ν̄P (cm −1 ) ν̄R (cm −1 ) ν̄P (cm −1 ) ν̄R (cm −1 ) ν̄P (cm −1 )

0 -41686.3551 -38800.5308 2906.185 2864.913 2906.566 2865.416 2905.995 2864.834

1 -41665.4436 -38780.1701 2925.983 2843.462 2928.046 2845.746 2925.581 2843.315

2 -41623.6326 -38739.4601 2945.205 2821.489 2950.129 2826.679 2944.577 2821.249

3 -41560.9488 -38678.4273 2963.839 2799.002 2799.002 2808.2156 2962.955 2798.641

4 -41477.4293 -38597.1097 2981.870 2776.016 2996.106 2790.356 2980.689 2775.499

5 -41373.1257 -38495.5589 2999.288 2752.543 3019.999 2773.099 2997.788 2751.817

6 -41248.1017 -38373.8377 3016.078 2728.595 3044.496 2756.446 3014.202 2727.624

7 -41102.4331 -38232.0237 3032.228 2704.185 3069.596 2740.396 3029.941 2702.907

8 -40936.2087 -38070.2054 3047.726 2679.324 3095.300 2724.949 3044.965 2677.697

9 -40749.5293 -37888.4831 3062.558 2654.025 3121.606 2710.106 3059.234 2651.932

10 -40542.5085 -37686.9711 3076.713 2628.300 3148.517 2695.866 3072.771 2625.689

11 -40315.2715 -37465.7953
14 Aditi and O.S.K.S. Sastri

TA B L E 4 Comparison of bound state energy eigen values of HCl molecule for Morse potential obtained by
matrix method with those from NU and AIM methods

Energy values(eV)

J MM NU AIM

v =0 0 -5.1684 -5.1684 -5.1684

5 -5.1296 -5.1296 -5.1296

10 -5.0266 -5.0266 -5.0266

v =1 0 -4.8106 -4.8106 -4.8106

5 -4.7728 -4.7728 -4.7728

10 -4.6726 -4.6725 -4.6725

v =2 0 -4.4657 -4.4657 -4.4657

5 -4.4289 -4.4289 -4.4289

10 -4.3314 -4.3313 -4.3313

where

2 1 h 2µD e J (J + 1) 2 2 3 i
√ = 2√ 2
− 2
√ ( − ) (35)
2 3 b 3 ~ Re 2 3 bR e b R e2
2

2µR e2 (D e + γD 2 ) −1 3
3 = ; D2 = + 2 and α = bR e (36)
~2 α 2 α α

and using AIM, analytical expression (Bayrak 2006) is

2
~2 h β 1 1 i
E nJ = 2
− (n + )α + γc 0 (37)
2µR e 2β 2 2

where

2µR e2 2µR e2
β 12 = (2D e − γc 1 ); β 22 = (γc 2 + 2D e ) (38)
~2 ~2

3 3 4 6 3 J (J +1)~2
and c 0 = 1 − + ; c1 = − ; c2 = −1
+ ; and γ = .
α α2 α α2 α α2 2µR e2
In Table 4, comparision of energy levels has been made between matrix method(MM), NU and AIM techniques for
J = 0, 5and 10 values, for v = 0, 1 and 2 values. It can be seen that the results obtained by matrix method are in good
agreement with those obtained from AIM and NU methods up to four decimal accuracy in the vibrational ground
state. A small variation of 1 unit in the fourth decimal place could be observed for J = 10 values in first and second
excited vibrational states in comparison to the analytical results.
Aditi and O.S.K.S. Sastri 15

5 | CONCLUSION:

The matrix methods numerical technique has been implemented to solve the Schrodinger equation for Morse po-
tential. The vibrational frequencies of HCl, obtained from simulation are matched with the experimental ones in the
least square sense by minimizing χ 2 −value using variational Monte-Carlo approach for the first time. The best op-
timized model parameters obtained within few runs of VMC algorithm are D e = 5.352160eV , b = 1.734139Å−1 and
R e = 1.274553Å. These are utilised for determining the spectral lines in the rotational fine structure of HCl and are
found to be matching with experimental data better than those resulting from analytical expressions fitted using mul-
tiple regression analysis. The energies obtained from simulation match very well with those resulting from analytical
expressions of NU and AIM methods confirming the effectiveness of the numerical technique implemented.

References

Modified Rosen-Morse potential-energy model for diatomic molecules. (2012). Physical Review A, 86.
Approximate solutions of the Dirac equation for the Rosen-Morse potential including the Spin-orbit centrifugal term.
(2010).Journal of Mathematical Physics, 51.
Minutes of the Middletown meeting, October 14, 1933. (1933).Physical Review, 44
Exact solutions of the Schrödinger equation. (1935).Physical Review, 48
Analytical approximations to the solutions of the Manning-Rosen potential with centrifugal term. (2007). Physics Letters
A, 368.
Approximate l-state solutions of the D-dimensional Schrödinger equation for Manning-Rosen potential. (2008). Annalen
der Physik, 17.
Bound state solutions of the Schrödinger equation for modified Kratzer’s molecular potential. (2006). Chemical Physics
Letters, 417.
Relativistic treatment of a spin-zero particle subject to a Kratzer-type potential,” American Journal of Physics. (2007).
American Journal of Physics, 75.
Exact solutions of the radial Schrödinger equation for some physical potentials. (2008). Central European Journal of
Physics, 6.
Exact solutions of the Schrodinger equation for the deformed hyperbolic potential well and the deformed four-parameter
exponential type potential. (2000). Physics Letters A, 275
Mapping of the five-parameter exponential-type potential model into trigonometric-type potentials. (2004). Journal of
Physics A, 37.
On the vibrations of polyatomic molecules. (1932). Physical Review, 42.
On some exactly solvable potentials derived from supersymmetric quantum mechanics. (1992). Journal of Physics A,
25.
Supersymmetric improvement of the Pekeris approximation for the rotating Morse potential. (2004). Chemical physics
letters, 394.
Exact solutions of the modified Kratzer potential plus ring-shaped potential in the D-dimensional Schrödinger equation
by the Nikiforov–Uvarov method. (2008). International Journal of Modern Physics C, 19.
Special functions of mathematical physics. (1988). Basel: Birkhäuser.
Rotational correction on the Morse potential through the Pekeris approximation and Nikiforov-Uvarov method. (2005).
arXiv preprint quant-ph/0502182.
16 Aditi and O.S.K.S. Sastri

Mapping of shape invariant potentials under point canonical transformations. (1992). Journal of Physics A: Mathemat-
ical and General. 25.
Approximate solutions of the Schrödinger equation with the generalized Morse potential model including the centrifugal
term. (2011). International Journal of Quantum Chemistry. 111.
Applications of automatic differentiation to the time-dependent Schrödinger equation. (2011). International Journal of
Quantum Chemistry. 111.
On the time-dependent solutions of the Schrödinger equation. I. The linear time-dependent potential. (2011). Interna-
tional Journal of Quantum Chemistry. 111.
Morse potential eigen-energies through the asymptotic iteration method. (2007). International Journal of Quantum
Chemistry. 107.
Arbitrary l-state solutions of the rotating Morse potential through the exact quantization. (2007). Physics Letters A.
363.
Bc meson spectrum and hyperfine splittings in the shifted large-N- expansion technique. (2003). International Journal
of Modern Physics A. 18.
Solving the vibrational Schrödinger equation on an arbitrary multidimensional potential energy surface by the finite
element method. (2009). Computer Physics Communications. 180.
Comparison of some special optimized fourth-order Runge–Kutta methods for the numerical solution of the Schrödinger
equation. (2005). Computer physics communications. 166.
Trigonometrically fitted fifth-order runge-kutta methods for the numerical solution of the schrödinger equation. (2005).
Mathematical and computer modelling. 42.
Explicit series solutions of some linear and nonlinear Schrodinger equations via the homotopy analysis method. (2009).
Communications in Nonlinear Science and Numerical Simulation. 14.
Approximate solution for the Klein–Gordon–Schrödinger equation by the homotopy analysis method. (2010). Chinese
Physics B . 19.
Approximate analytical versus numerical solutions of Schrödinger equation under molecular Hua potential. (2012).
International Journal of Quantum Chemistry. 112.
Bound State Solutions of D-Dimensional Schrodinger Equation with Exponential-Type Potentials. (2014). International
Journal of Quantum Chemistry. 115.
The harmonic oscillator in quantum mechanics: A third way. (2009). American Journal of Physics. 77.
Solving for three-dimensional central potentials using numerical matrix methods. (2009). American Journal of Physics.
81.
Simulation study of nuclear shell model using sine basis. (2020). American Journal of Physics. 88.
Numerical simulation of quantum anharmonic oscillator, embedded within an infinite square well potential, by matrix
methods using gnumeric spreadsheet. (2020). American Journal of Physics. 41.
[Colin N Banwell]Fundamentals of molecular spectroscopy. (1972). McGraw-Hill, New York.
Analysis of the infrared spectra of diatomic molecules. (1999). Journal of Chemical Education, 76.
Fundamentals of molecular spectroscopy. (1989). John Wiley and Sons, New York.
Diatomic molecules according to the wave mechanics. II. Vibrational levels. (1929). Phys. Rev., 34. Asymmetric
wave functions from tiny perturbations. (2015). Am. J. Phys., 83.
Numerical matrix method for quantum periodic potentials. (2016). Am. J. Phys., 84. Molecular spectra and molecular
structure: IV. Constants of diatomic molecules. (2013). Spectra of atoms and molecules. (2020). Oxford university press.
Modified Rosen-Morse potential-energy model for diatomic molecules. (2012). Physical Review A, 86.
Aditi and O.S.K.S. Sastri 17

Rotational correction on the Morse potential through the Pekeris approximation and Nikiforov-Uvarov method.
(2005). ArXiv Preprint Quant-Ph/0502182(.
Arbitrary -state solutions of the rotating Morse potential by the asymptotic iteration method. (2006). Journal of
Physics A: Mathematical and General, 39.
Simulation of Vibrational Spectrum of Diatomic Molecules Using Morse Potential by Matrix Methods in Gnumeric
Worksheet. (2020). Physics Education.
Numerical Solution of Square Well Potential With Matrix Method Using Worksheets. (2019). Physics Education,
36.