Quantum Simulation of non-Born-Oppenheimer dynamics in molecular systems by path integrals

TL;DR

This paper introduces a path integral-based numerical algorithm to accurately compute non-Born-Oppenheimer energies in molecular systems, providing a nonperturbative alternative to traditional methods.

Contribution

The authors develop a novel path integral approach combined with variational functions to treat non-Born-Oppenheimer molecular dynamics without relying on perturbation theory.

Findings

Calculated non-Born-Oppenheimer energies for hydrogen molecules and ions.

Results agree well with existing theoretical and experimental data.

Demonstrated the method's potential as a nonperturbative testing tool.

Abstract

A numerical algorithm based on the probabilistic path integral approach for solving Schroedinger equation has been devised to treat molecular systems without Born-Oppenheimer approximation in the non relativistic limit at zero temperature as an alternative to conventional Variational and perturbation methods. Using high quality variational trial functions and path integral method based on Generalized Feynman-Kac method, we have been able to calculate the non-Born-Oppenheimer energy for hydrogen molecule for the sigma state and hydrogen molecular ion. Combining these values and the value for ionization potential for atomic hydrogen, dissociation energy and ionization potential for hydrogen molecule have been determined to be 36 113.672(3) cm inverse and 124.446.066(10) cm inverse.respectively. Our results favorably compare with other theoretical and experimental results and thus show the promise of being a nonperturbative alternative for testing fundamental physical theories.

Tables

Table 1. Table 1: Notation Table

Notation/Phrase	Meaning
$r_{\beta }$	position vector of electrons inside the atom
$R_{\alpha }$	position vector of nuclei
$R$	distance between two nuclei
$r$	distance between two electrons
$X\left(t\right)$	Brownian motion with a non-ergodic probabilistic measure or
	Wiener Measure
$Y\left(t\right)$	A stochastic process with an ergodic or stationary measure
${\psi }_T$	The trial function corresponding to mathematical ground state
${\varphi }_T$	The trial function for the molecular ion
$\hat{T}$	Kinetic energy operator
$\hat{V}$	Potential Energy operator
	Virial theorem
$E_p\left(H\right)$	ionization potential of atomic hydrogen $=0.5a.u.$
$E\left(H_2\right)$	Total energy for the hydrogen molecule
$E\left(H_2^{}{}_{}{}^{+}\right)$	Total energy for the hydrogen molecular ion
$E_d\left(H_2\right)$	dissociation energy for the hydrogen molecule
$E_p\left(H_2\right)$	ionization potential for the hydrogen molecule

Table 2. Table 2: Born-Oppenheimer Energy of Hydrgen molecule for the ground state for different internuclear distances. The number in the parentheses is the statiscal error.

R	E	Refs.
0.4	-0.122 348(5)	This work
	-0.120 230	Sims et al[54]
	-0.120 228 2(9)	Alexander et al[55]
0.6	-0.771 535(1)	This Work
	-0.769 635	Sims et al[54]
	-0.769 635 1 (4)	Alexander et al[55]
0.8	-1.021 424(6)	This Work
	-1.020 056	Sims et al[54]
	-1.020 056 1(4)	Alexander et al[55]
1.0	-1.125 406	This work
	-1.124 539	Sims et al[54]
	-1.124 539 2(4)	Aleaxnder et al[55]
1.2	-1.165 377(1)	This work
	-1.164 935	Sims et al[54]
	-1.164 934 8(5)	Aleaxnder at al[55]
1.4	-1.174 564(2)	This work
	-1.174 475	Sims at al[54]
	-1.174 475 4(6)	Aleaxnder et al[55]
	-1.174 475 1(5)	Datta et al[56]
	-1.174 447 477	Kolos et al[17]
	-1.174 475 686	Kolos[57]
1.6	-1.168 371(1)	This work
	-1.168 583	Sims et al[54]
	-1.168 583 4(5)	Alexander et al[55]
2.0	-1.137 488(2)	This work
	-1.138 132	Sims et al[54]
	-1.138 132 0(4)	Aleaxnder et al[55]
3.0	-1.057 351(1)	This work
	-1.057 326	Sims et al[54]
	-1.057 324 9(3)	Alexander et al[55]
4.0	-1.019 750(4)	This work
	-1.016 390	Sims et al[54]
	-1.016 389 2(2)	Alexander et al[55]

Table 3. Table 3: Total nonrelativistic Born-Oppenheimer Energy of Hydrgen molecule ion and hydrogen molecule for the ground state at equilibrium distances. The number in the parentheses is the statiscal error.

molecule	R	E
$H_2^{+}$	2.0	-0.609 148(1)	This work
	2.01	-0.602 1	Swarwono et al[59]
$H_2$	1.4	-1.174 564(2)	This work

Table 4. Table 4: Total non-relativistic non Born-Oppenheimer energy in (a.u.) of hydrogen molecular ion for the ground state with Generalized Feynman-Kac(GFK) path integral method at different time with a stepsize of 1 / 30 1 30 1/30 and 5000 paths.

Time	$E_{mol}\left(GFK\right)$			Virial ratio
8	-0.554700(1)	0.59921	-0.92368	1.541
16	-0.553210(4)	0.78576	-1.4218	1.809
24	-0.458169(2)	0.3784	-0.6724	1.776
32	-0.431555(6)	0.7140	-1.2891	1.805
40	-0.420000(4)	0.616	-1.036	1.681
48	-0.376000(6)	0.417	-0.793	1.9

Table 5. Table 5: Total Nonrelativistic non-Born-Oppenheimer enegies E i ​ o ​ n subscript 𝐸 𝑖 𝑜 𝑛 E_{ion} of Hydrgen molecule ion in the ground state. The number in the parentheses is the statistical error.

Work	Method	$E_{ion}$ (a.u.)
Yuh et al	free iterative complement method(variational)	-0.597 139
Jeziorski et al	rel variational	-0.597 144
present work	GFK	-0.597 528(2)

Table 6. Table 6: Total non-relativistic non Born-Oppenheimer energy in (a.u.) of hydrogen molecule for the ground state with Generalized Feynman-Kac(GFK) path integral method at different time with a stepsize of 1 / 30 1 30 1/30 and 10000 paths.

Time	$E_{mol}\left(GFK\right)$			Virial ratio
8	-1.164376(5)	1.17594	-2.159756	1.836
16	-1.159098(4)	1.136411	-2.23381	1.97
24	-1.163043(1)	1.266482	-2.275504	1.8
32	-1.165153(2)	0.9991487	-1.90263	1.904
40	-1.164310(5)	1.235915	-2.584469	2.09
48	-1.15559(2)	1.120355	-2.27630	2.03

Table 7. Table 7: Total Nonrelativistic non-Born-Oppenheimer enegies E m ​ o ​ l subscript 𝐸 𝑚 𝑜 𝑙 E_{mol} of Hydrgen molecule in the ground state . The number in the parentheses is the statiscal error.

Work	Method	$E_{mol}$ (a.u.)
Tubman et al[58]	FN DMC-full	-1.164 01(5)
Bubin et al[36]	non-BO-Var	-1.164 0250
Alexander et al[40]	non-BO-Var	-1.164 02491(8)
present work	GFK	-1.164 546(3)

Table 8. Table 8: ionization energy E p subscript 𝐸 𝑝 E_{p} ,dissociation energy E d subscript 𝐸 𝑑 E_{d} of hydrogen molecule in c m − 1 ) {cm}^{-1}) for the ground state. The number in the parentheses is the statiscal error.

Work	method	$E_p$	$E_d$
Herzberg et al[18]	Expt		36 113.6$\pm 0.3$
Herzberg[19]	Expt		36 116.3
Wolneiwicz[57]	Theo(BO)	124 417.491
Zhang et al[26]	Expt		36 118.062(10)
Liu et al[30]	Hybrid
	Expt
	-Theo	124 417.491 13(37)	36 118.069 62(37)
Altmann et al[31]	Expt		36 118 06945(31)
Piszczatowski
et al[27]	Var
	(Theo)		36 118.0695(10)
Stwalley[20]	Expt		36 118.6
Pachuki et al[28]	Var
	(Theo)		36 118.797 746 3(2)
Puchalskii et al[29]	(Theo)		36 118.067 8(6)
Cheng et al[32]	Expt	124 357.238062(25)	35 999.582 894 (25)
Wang et al[38]	Var(Theo)		36 118.069 71(33)
Present work	GFK(Theo)(non rel)	124 446.066 (10)	36 113.672(3)
Present work	GFK with rel correction
	(energy nBO+rel corr BO)		36 116.072(10)

Equations

$$H=\sum _{\alpha =A,B}-\frac{{\nabla }_{\alpha }^2}{2\frac{M}{m}}-\sum _{\beta =1,2}-\frac{{\nabla }_{\beta }^2}{2}-\frac{1}{r_{1A}}-\frac{1}{r_{1B}}-\frac{1}{r_{2A}}-\frac{1}{r_{2B}}+\frac{1}{r_{12}}+\frac{1}{R}$$

$$H\Psi (r_A,r_B,R)=E\psi (r_A,r_B,R)$$

$$H=\sum _{\alpha =A,B}-\frac{{\nabla }_{\alpha }^2}{2\frac{M}{m}}-\frac{{\nabla }^2}{2}-\frac{1}{r_A}-\frac{1}{r_B}+\frac{1}{R}$$

$$i\frac{\partial \psi (x,t)}{\partial t}=(-\frac{{\nabla }^2}{2}+V)\psi (x,t)$$

$$\psi (x,0)=f(x)$$

$$\frac{\partial \psi (x,t)}{\partial t}=(\frac{{\nabla }^2}{2}-V)\psi (x,t)$$

$$\psi (x,0)=f(x)$$

$$\psi (x,t)=E_{x}exp\{-{\int }_0^{t}V(X(s))ds\}$$

$${\lambda }_1=-\underset{t\to \infty }{\lim }\frac{1}{t}ln{E}_x[e^{-{\int }_0^{t}V(X(s))ds})]$$

$${\lambda }_1={\lambda }_T-\underset{t\to \infty }{\lim }\frac{1}{t}ln{E}_x[e^{-{\int }_0^{t}V(Y(s))ds})]$$

$$dY(t)=\frac{\nabla {\Psi }_T(Y(t))}{{\Phi }_T(Y(t))}+dX(t)$$

$$⟨Y|A|Y⟩=\frac{{\lim }_{t\to \infty }\int dY(t)A(Y(t))e^{-\int [V_p(Y(s)]ds}}{\int dY(t)e^{-\int [V_p(Y(s)]ds}}.$$

$$=-\frac{{\hslash }^2}{2M}\sum _{\alpha }{\nabla }_{\alpha }^2-\frac{{\hslash }^2}{2m}\sum _{\beta }{\nabla }_{\beta }^2-\sum _{\alpha ,\beta }\frac{e^2}{4\pi {ϵ}_0{r}_{\alpha \beta }}+\sum _{A>B}\frac{e^2}{4\pi {ϵ}_0{R}_{AB}}+\sum _{\beta >\gamma }\frac{e^2}{4\pi {ϵ}_0{r}_{\beta \gamma }}$$

$$H=-\frac{1}{2M}\sum _{\alpha }{\nabla }_{\alpha }^2-\frac{1}{2}\sum _{\beta }{\nabla }_{\beta }^2-\sum _{\alpha ,\beta }\frac{1}{r_{\alpha \beta }}+\sum _{A>B}\frac{1}{R_{AB}}+\sum _{\beta >\gamma }\frac{1}{r_{\beta \gamma }}$$

$$H=-\frac{1}{2}\sum _{\beta }{\nabla }_{\beta }^2-\sum _{\alpha ,\beta }\frac{1}{r_{\alpha \beta }}+\sum _{A>B}\frac{1}{R_{AB}}+\sum _{\beta >\gamma }\frac{1}{r_{\beta \gamma }}$$

$$H=-\frac{1}{2m}{\nabla }^2(r_1,r_2)-\frac{1}{2M}{\nabla }^2(R_A,R_B)$$

$$+\frac{1}{|r_1-r_2|}+\frac{1}{|R_A-R_B|}$$

$$-\frac{1}{|r_1-R_A|}-\frac{1}{|r_1-R_B|}-\frac{1}{|r_2-R_A|}-\frac{1}{|r_2-R_B|}$$

$$=-\frac{1}{2}{\nabla }^2+V(r_1,r_2,R_A,R_B)$$

$$[-\sum _i\frac{{\hslash }^2}{2m_i}{\nabla }^2({\vec{r}}_i)-\sum _j\frac{{\hslash }^2}{2M_j}{\nabla }^2({\vec{R}}_j)$$

$$+\sum _{ik}\frac{1}{\sqrt{(r_i-r_k{)}^2}}+\sum _{jl}\frac{1}{\sqrt{(R_j-R_l{)}^2}}-\sum _{ij}\frac{1}{\sqrt{(R_j-r_i{)}^2}}]\psi ({\vec{r}}_i,{\vec{R}}_j)=\mu \psi ({\vec{r}}_i,{\vec{R}}_j)$$

$$[-\sum _i\frac{{\nabla }^2({\overrightarrow{r_i}}^{\prime })}{2}-\sum _j\frac{{\nabla }^2({\overrightarrow{R_j}}^{\prime })}{2}$$

$$+\sum _{ik}\frac{1}{\sqrt{({\overrightarrow{r_i}}^{\prime }-{\overrightarrow{r_k}}^{\prime }{)}^2}}+\sum _{jl}\frac{1}{\sqrt{({\overrightarrow{R_j}}^{\prime }-{\overrightarrow{R_l}}^{\prime }{)}^2}}-\sum _{ij}\frac{1}{\sqrt{(\frac{m_i}{M_j}{\overrightarrow{R_j}}^{\prime }-{\overrightarrow{r_i}}^{\prime }{)}^2}}]\psi ({{\vec{r}}^{\prime }}_i,{{\vec{R}}^{\prime }}_j)=\mu \psi ({{\vec{r}}^{\prime }}_i,{{\vec{R}}^{\prime }}_j)$$

$$H=-\frac{1}{2}{\nabla }^2(r_1,r_2)-\frac{1}{2}{\nabla }^2(R_A,R_B)$$

$$-\frac{1}{\sqrt{(\sqrt{\frac{m}{M}}{X}_A-x_1{)}^2+(\sqrt{\frac{m}{M}}{Y}_A-y_1{)}^2+(\sqrt{\frac{m}{M}}{Z}_A-z_1{)}^2}}$$

$$-\frac{1}{\sqrt{(\sqrt{\frac{m}{M}}{X}_B-x_1{)}^2+(\sqrt{\frac{m}{M}}{Y}_B-y_1{)}^2+(\sqrt{\frac{m}{M}}{Z}_B-z_1{)}^2}}$$

$$-\frac{1}{\sqrt{(\sqrt{\frac{m}{M}}{X}_A-x_2{)}^2+(\sqrt{\frac{m}{M}}{Y}_A-y_2{)}^2+(\sqrt{\frac{m}{M}}{Z}_A-z_2{)}^2}}$$

$$-\frac{1}{\sqrt{(\sqrt{\frac{m}{M}}{X}_B-x_2{)}^2+(\sqrt{\frac{m}{M}}{Y}_A-y_2{)}^2+(\sqrt{\frac{m}{M}}{Z}_A-z_2{)}^2}}$$

$$+\frac{1}{\sqrt{(X_A-X_B{)}^2+(Y_A-Y_B{)}^2+(Z_A-Z_B{)}^2}}$$

$$+\frac{1}{\sqrt{(x_1-x_2{)}^2+(y_1-y_2{)}^2+(z_1-z_2{)}^2}}$$

$$E_P=E(H_2^{+})-E(H_2)$$

$$\begin{gathered} E_d=-\underset{R\to \infty }{\lim }E(H_2)-E(H_2) \\ =-2E(H)-E(H_2) \end{gathered}$$

$${\psi }_T=(1+P_{12})(1+P_{AB})exp({\sum }_{k=0}{a}_k{q_{1A}}^u{q_{1B}}^v{q_{2A}}^w{q_{2B}}^n{q_{12}}^g{q_{AB}}^h-\chi r_{1A}-\delta r_{2B})$$

$${\varphi }_T=exp(-\sigma r^2)$$

$$W(l)\equiv W(t,n,l)$$

$${w_j}^i(t,n,l)=\sum _{k=1}^l\frac{{ϵ}_{jk}^i}{\sqrt{n}}$$

$$\underset{n\to \infty }{\lim }P(\frac{1}{n}\sum _{l=1}^{nt}V(W(l)))\le \nu$$

Full text

1 Introducton

Sovling eigenvalue problems for molecular systems are in general complicated and people quite often make several simplications to handle the difficulties associted with it, namely Born-Oppenheimer(BO) approximation or adiabatic Born-Oppenheimer(ABO)[1,2,3] and rovibrational approach[4]. In the Born-Oppenheimer(BO) approximation in a chemical system, the coupling between the nuclear and electronic movements is neglected. The large mass disparity in the nuclei and electrons justifies the decoupling of their different time scale motions and provides a very practical way to model a lot of chemical system adequetely. In quantum mechanics, a solution of Schrödinger equation in the adiabatic approximation is defined to have a time scale separation in the fast and slow degrees of motion. Sometimes the BO approximation is also referred to as ’Adiabatic Born-Oppenheimer’ approximation as the lighter electrons follow the motion of the heavier nulei adiabatically. The ABO approximation breaks down when two or more potential energy surfaces approach each other or cross and one must take resort to the coupled equations. It is justifiable when energy gap between the ground and excited states is larger than the energy scale of nucleus. In metals, the applicability of ABO is quaestionble as this energy gap turns out to be zero. For example, ABO fails in the case of Graphene[5].

With the typical BO approximation in a diatomic molecule the non-relativistic ground state energy can be evaluated by solving the relevant Schrödinger equation neglecting the nuclear kinetic energy. One can choose some fixed value for the nuclear confugurations and solve for the electronic wavefunction which depends parametrically on this fixed value of nuclear configuration which we will describe in the next section. Since the nuclear configuration is considered as a parameter and not a quantum mechanical variable the nuclear motion is modelled only classically. To get the full quantum dynamics of the diatomic molecule one needs to treat molecular systems as a whole, including electrons and atom nuclei on the same footing. Or in ther words one needs to consider the motion of all the constituents of the molecule simultaneously assuming nuclei have finite masses and they move in the configuration space as freely as the electrons do. In view of the above reasons in this paper we have taken a relook at the sigma state of hydrogen molecule and hydrogen ion molecule as a testbed so explore the fully non-Born-Oppenheimer(nBO) or non-adiabatic efects in diatomic molecular systems in general. As a matter of fact we have adopted a quantum Monte Carlo method based on Generalized Feynman-Kac(GFK) method[6-11] to calculate the energies for the sigma state of hydrogen molecule and molecular hydrogen. Since GFK is a non-perterbative approach, it is easier to study motion of all the particles in the molecular system even in a fully quantum mechanical scenario. In the framework of GFK the non-BO study of hydrogen molecule now boils down to solving a four particle Schrödinger equation. In the Stochastic scenario, the all four particles execute twelve dimensional random walk scaled differently due to mass disparity of electrons and nuclei. In the case of a hydrogen molecular ion it turns out to be a nine dimensional walk of three particles(one electron and two nuclei). Adopting this idea we prescribe the nBO model for any aribitrary many body system with more than one electron and nucleus. Now to test the power of our theory we calculate the total energy of hydrogen molecular system(hydrogen molecule and hydrogen molecular ion) using path integral Monte Carlo technique use those to determine the ionization potential $E_{p}$ and dissociation energy $E_{d}$ of hydrogen molecule. At this point we need to review the theoretical and experimental endeavors and their outcomes so far. Hydrogen molcule and hydrogen molecular ion are well explored topics in quantum mechanics. Their long history dates back to the first theoretical work of Heitler and London[12] followed by the work of James and Coolidge[13,14] and extension by Kolos and Wolniewicz[15-17]. The controversy of [16,17,18] was apparently resolved by the experimental results of Hertzberg[19] and Stwalley[20]. Subsequent theoretical[21,22,23] and experimental endeavors[24,25,26] seemed to reduce the discrepency between the theory and experiment. Up to 2017 the best accepted theoretical value of the dissociation energy of hydrogen molecule in the non-BO basis is 36,118.0695(10)$cm^{-1}$[Piszczatowski et al [27],Pachuki et al [28] and Puchalski[29] and the best experimental value was 36,118.0696(4)$cm^{-1}$[Liu et al[30] and Altmann et al[31]. One of the most recent experimental works[32] has reported the values for ionization potentials($E_{p}$) and dissociation energies($E_{d}$)for hydrogen molecule which are significantly lower than the accpted values so far. Most of the theoretical approaches[33-39] were based on variational principle and provided the sophistcated upper bound to the total energies, dissociation energy and ionization energy etc. In this work using Non-BO basis functions as the trial functions[40] our path integral approach yielded the nonrelativistic non-BO energies for both hydrogen molecule and molecular ion which are lower than the previously accepted values. These ground state energies for the hydrogen molecular system were used in calculating the ionization potentials($E_{p}$) and dissociation energies($E_{d}$) for hydrogen molecule and we get new benchmarks for those.

The orgnization of the paper is as follows: In Section 1, we introduce the problem and describe the contents of the different sections of the paper. In Section 2.1, we discuss the general aspects of BO and nBO approaches to Quantum Mechanical problems. In Section 2.2 we discuss our path integral approach for calculating the eigen energies of hydrogen molecular system. In Section 3, we show how it can be generalized to any N particle molecular systems. In Section 4 we discuss our results and we conclude in Section 5. Fig 1 and 2 are schematic diagrams for the overall motion of the molecule in the nBO scenario. Fig 3 show how our numerical code works in the BO limit.

2 Theory

2.1 General theoretical considerations in BO and nBO approach in connection with hydrogen molecular system

In the BO approximation the nuclei are clamped at a fixed position and only electrons are moving in the configuration space. In the following Figs 1 and 2 the arrows beside the electrons as well as the nuclei signify that in the non- BO approach electrons and nuclei all are moving in the configuration space and they are being treated in the same footing.

The exact time dependent Schrödinger equation for hydrogen molecule can be represented by $H\Psi(r_{1},r_{2},R_{1},R_{2})=E\Psi(r_{1},r_{2},R_{1},R_{2})$ with

[TABLE]

where $r_{\beta},{\beta}=1,2,R_{\alpha},{\alpha}=A,B$ are the positions of the electrons and the nuclei respectively. As a matter of fact $r_{1A}=r_{1}-R_{A},r_{1B}=r_{1}-R_{B}$ and so on $r=r_{1}-r_{2},R=R_{A}-R_{B}$. In defining the above Hamiltonian we have used atomic units and m and M denote the mass of electrons and nuclei respectively. Also $\frac{M}{m}=1836.152701$. The above Schrödinger equation describes the quantum motion of the four constituents of hydrogen molecule(two electrons and 2 nuclei). For the hydrogem molecular ion the Schrödinger equation can read as

[TABLE]

[TABLE]

In the case of Born-Oppenheimer approximation, corresponding to the general Hamiltonian for the electron-nuclei system of hydrogen molecule $H(R,r)=T_{R}+H_{e}(R,r)$ there will be only one potential surface defined by $V(R)=<{\psi}_{0}|H_{e}(R,r)|{\psi}_{0}>$(as evident from Fig 3), whereas in the case of a non-BO calculation there will be multiple potential surfaces[41] corresponding to $V_{ii}(R)=<{\psi}_{i}|H_{e}(R,r)|{\psi}_{i}>$ and $V_{ij}(R)=<{\psi}_{i}|H_{e}(R,r)|{\psi}_{j}>$.

2.2 General Path Integral Theory for the Energy and other Properties for many body systems

Let us consider the time-dependent Schroedinger equation for a system of N particles with Hamiltonian $H=-{{\nabla}^{2}}/2+V(x)$ as follows:

[TABLE]

where $x\in R^{d}$. The above initial value problem in imaginary time can be represented as

[TABLE]

The solution of Eq(5) for $f(x)=1$ can be written in the following Feynaman-Kac representation[42-43] and it provides a rigorous justification unlike ordinary Feynman path integration[44-47].

[TABLE]

for $V\in K_{\nu}$, the Kato class of potential[48] where $x(t)$ is a Brownian motion trajectory and E is the average value of the exponetial term with respect to these trajectories. The lowest eigenvalue for a given symmetry can be obtained by large deviation principle of Donsker and Varadhan[49],

[TABLE]

The above representation(Eq[6]) suffers from poor convergence rate as the underlying diffusion process-Brownian motion(Wiener Process[50] is non-recurrent. So it is necessary to use a representation which employs a diffusion which unlike Brownian motion, has a stationary distributions. To speed up the calculation we use generalized Feynman-Kac method in which the Hamiltonian is rewritten as $H=H_{0}+V_{p}$ where $H_{0}=-{{\nabla}^{2}}/2+{\lambda}_{T}+{\nabla}{\psi}_{T}/2{\psi}_{T}$. Here ${\psi}_{T}$ is a twice differentiable nonnegative reference function and $H{\psi}_{T}={\lambda}_{T}{\psi}_{T}$. The expression for the energy can now be rewritten as

[TABLE]

where $Y(t)$ is the diffusion process which solves the stochastic differential equation and is known as Ornstein-Uhlenbeck process.

[TABLE]

$V(Y(s))$ is summed over all the time steps and $e^{-V(Y(s))}$ is summed over all the trajectories. The presence of both drift and diffusion terms in the above expression enables the trajectory to be highly localized. As a result, the important regions of the potential are frequently sampled and Eq(4) converges rapidly. The expectation value for the other properties can be evaluated as follows[6,51]:

[TABLE]

3 Calculation BO and nBO energies and properties

The nonrelativistic Hamiltonian for hydrogen molecule can be written as a sum of five terms:

[TABLE]

In atomic units the above expression is just

[TABLE]

For Born-Oppenheimer calculations the kinetic energy of nuclei is neglected due to its small contribution towards the total energy as it gets divided by ’M’, the mass of nuclei. So within BO approximation,the above equation reads as follows:

[TABLE]

Here ${R_{A}}$ and ${R_{B}}$ are not variables hence $R_{AB}$ is treated as a parameter.

In general the Hamiltonian for all the moving electrons and the nuclei has the form

[TABLE]

where $r_{\beta},{\beta}=1,2,R_{\alpha},{\alpha}=A,B$ are the positions of the electrons and the ’quantum nuclei’(moving nuclei) respectively(Fig 2). $|.|$ is the 3 dimensional Euclidean distance and $V$ is the Coulombic interaction. For non-Born-Oppenheimer calculations the hydrogen molecules can be treated as general physical systems with inequivalent masses(electrons and nuclei) and can be represented with the above Hamiltonian. For any physical system with N ineqivalent masses the above can be generalized as follows: The Schrödinger equation for the above system can be written as

[TABLE]

i,j refer to number of elctrons and number of nuclei respectively. Now using $\vec{r_{i}}=s_{i}\vec{r_{i}}^{\prime}$ $\vec{R_{j}}=s_{j}\vec{R_{j}}^{\prime}$, ${\nabla}^{2}(\vec{r}_{i})={\nabla}^{2}(s_{i}\vec{r_{i}}^{\prime})=\frac{1}{{s_{i}}^{2}}{\nabla}^{2}(\vec{r_{i}}^{\prime})$

${\nabla}^{2}(\vec{R}_{j})={\nabla}^{2}(s_{j}\vec{R_{j}}^{\prime})=\frac{1}{{s_{j}}^{2}}{\nabla}^{2}(\vec{R_{j}}^{\prime})$ and putting $m_{i}{s_{i}}^{2}=M_{j}{s_{j}}^{2}$, multiplying throughout by $\frac{m_{i}{s_{i}}^{2}}{\hbar^{2}}$ the Schrödinger equation in the dimensionless form reads as (in $m_{i}=\hbar=1$ units)

[TABLE]

In general let $m_{i}{s_{i}}^{2}=M_{j}{s_{j}}^{2}$ and set $s_{i}=\sqrt{\frac{m_{i}}{M_{j}}}s_{j}$ No matter the number of distinct masses, the scale for the corresponding random walker will always be the square root of the ratio of their masses. Using this result, the Generalized Feynman-Kac path integral can be used for non adiabatic treatment of the molecules. Now for hydrogen molecule we need to simulate the random ealk associated with the following Hamiltonian

[TABLE]

4 Results and discusssions

Now by setting $x=\{x_{k}\},k=1,2,....,12,{x_{k}}{\bf{{\in}R^{12}}}$ we can write $V(x)=V(r_{1},r_{2},R_{A},R_{B})$ To calculate energies, we use Eq(8) of Sec 2.2 whereas the other properties are calculated using Eq(10) with $V$ as defined in Eq(16). In our program the stepsize is fixed and the direction of the path is chosen randomly. In Eq(9),the first and second terms represent the drift and diffusion respectively. For each system a number of paths are generated with a specific path length. These are then summed to produce an average value and a statistical error. In order to examine the behavior of our energies and other properties as a function of path length, we compute several different path lengths - from 8 to 48 units of time. To implement Eq(8) numerically we replace 12 dimensional Brownian motion with 12 dimensional Ornstein-Uhlenbeck process and simulate them by 12 independent, properly scaled one dimensional drifted random walk. For details plesae see Appendix A . The Hamiltonian in Eq(17) is represented as twelve dimensional random walk where six dimensions involve stepsizes $\frac{m}{M}$ of the size of other six as shown above. We calculate $E(H_{2})$ and $E(H_{2}^{+})$ using the formula given Eqn(8). Using these values and the value of ionization potential of hydrogen one can calculate the Ionization Potential and Dissociation Energy according to the following expressions[52].

**Ionization Potential:

**

[TABLE]

**Dissociation Energy:

**

[TABLE]

For the BO and NBO calculations for the hydrogen molecule we use the trial function of the following form:

[TABLE]

Here $P_{12}$ is the operator that interchanges the two electrons, $P_{AB}$ is the operator that interchages the two nuclei, and $q_{x}={r_{x}}/(1+c{r_{x}})$ is a coordinate transformation that allows terms in the exponent to go smoothly to the separated atom limit. The exponents u, v, w, n, g,and h are integers(0,1,……..) and all possible terms adding up to $N=u+v+w+n+g+h$ are selected. In Table 2, we show the variation of BO energies with internuclear distances and its agreement with the best nonrealtivistic estimates for this system. We show this agreement for BO case just to establish that the same code works when the nuclei are frozen as well. In Table 3, we show extrapolated values for the BO energies for hydrogen molecule and hydrogen molecular ion. In Table 4, we run simulations for ${H_{2}}^{+}$ for different time 8-48. In Table 5, we show the comparison of extrapolated energy from the data in Table 4 with the best non relativistic variational calculatios. Table 6 contains the simlations for nbo dynamics of $H_{2}$ molecule at different time (8-48). In Table 7 we show the extrapolated value of the nbo energy of hydrogen moecule and its agreement with other theoretical values. As can be seen from Eq(8) the most accurate estimate of the energy is obtained when we extrapolate our results to infinite time. We do this by performing a least square fit. We have verified that the 5000 path lengths selected with runtime selected from 8 to 48 are more than enough to perform an accurate fit. Oher tests have confirmed that a stepsize of 1/30 has little influence on the value of extrapolated energy. The final value obtained for the nBO energy for the hydrogen molecule,-1.164 546(3), is in excellent agreement with the best non relativistic value for this system[36,58]. It is, better than the value obtained from the Variational Monte Carlo calculation. This agreement can, however, be attributed to the quality of the original trial wavefunction. Other tests have confirmed that the stepsize used has little influence on these values. Unlike the energy there is no need to extrapolate any of the properties to infinite time In Table 8, we have compared our Ionization potential $E_{p}$ and dissociation energy $E_{d}$ with other theoretical and experimental results.

For the BO and NBO calculations for the hydrogen molecular ion the following trial function function is used:

[TABLE]

where $\sigma$ is a variational parameter. Using the calculated values for $E(H_{2}^{+})$(The extrapolated nBO energy for hydrogen molecular ion from Table 5), $E(H_{2})$ (The extrapolated nBO energy for hydrogen molecule Table 7)and Eq(18), we calculate the value for the ionization potential $E_{p}(H_{2})=0.567018a.u.=124446.0663cm^{-1}$. Using ionization energy for atomic hydrogen $E_{p}(H)$, and $E_{(}H_{2})$ and Eq(19), we calculate the dissociation energy $E_{d}(H_{2})=0.164546a.u.=36113.672cm^{-1}$. Now adding the value of previously calculated relativistic corrections[53] for sigma state to our nBO energy we determine the $E_{d}(H_{2})$ to be equal to 36,116.672(10).

Acknowledgements: The author would like to thank Alliance University for providing partial support for carrying the research work and The University of Texas at Arlington,USA where the idea behind the research work was conceived.

5 Conclusions and Outlook

In this paper using a probabilistic approach a solution to time dependent Schrödinger equation has been constructed in the form of a path integral. By simulating an approximation to Ornstein-Uhenlenbech process through trial function(drift term) and a toss of a unbiased coin(diffusion term), we have determined the ionization potential and dissociation energy of hydrogen molecule with a high accuracy. It looks very promising to observe that if we guide our random walk using a non-BO variational trial function and perform the numerical simulation for the path integarl solution for the molecular system, we already improve the variational energy for the hydrogen molecule. For implementing the simulation of our path integral solution we make only onefold approximation as the probalistic representation to Schrödinger equation can be written in a closed form. To be precise to calculate energy we approximate an exact solution (i. e. the GFK representation of it) to the Schrödinger equation, whereas most of the other numerical procedures approximate a solution to an approximate Schrödinger equation. Also since the path integral solution is based on fully quantum mechanical approach it can improve the variational energy to a graet extent provided our trial function has the right symmetry of the Hamiltonian. At this point our dissociation energy a little less than established theoretical and experimental values. We believe we need to add a non-BO relativistic correction to our present nBO non-relativistic energy to have a better agreement with the experimental values. We see a better agreement if we add our BO relativistic corrections[53] as a rough estimate of non-BO relativistic corrections. Also we need to improve the quality of the trial funtions particularly in the case of hydrogen molecular ion. The procedre can be applied to more complex systems for which energies are known up to a few significat figures from variational calculations and accuracy can be increased to include more significant figures. Our benchmarks for the dissociation energy and ionization potential can be a useful input for the other work for evaluating non-BO relativistic corrections for the molecular systems. We observe that adopting this Monte Carlo method and taking advantage of modern computer technology solving eigenvalue calculations in Quantum mechanics can be simplified to a great deal and hope it will inspire other people to carry out research along the same line.

Appendix A Details of Numerical Calculations

The formalism described in section 2 can include any generalized potential [60] and valid for any arbitrary dimension d (d=3N). To implement Eq(3) numerically, the 3N dimensional Brownian motion can be replaced by properly scaled one dimensional random walks as follows [9, 43, 61]:

[TABLE]

where

[TABLE]

with ${w_{j}}^{i}(0,n,l)=0$ for $i=1,2,....,N$;$j=1,2,3$ and $l=1,2,.....,nt$. Here $\epsilon$ denotes the binomially distributed random variables which are chosen independently and randomly with probability P for all i,j,k such that $P({\epsilon}^{i}_{jk}=1)$=$P({\epsilon}^{i}_{jk}=-1)$=$\frac{1}{2}$. It is known by an invariance principle [62] that for every $\nu$ and W(l) defined in Eq.(A1) and Eq(A2)

[TABLE]

Consequently for large n,

[TABLE]

Finally, by generating $N_{rep}$ independent realization $Z_{1}$,$Z_{2}$,….$Z_{N_{rep}}$ of

[TABLE]

and using the law of large numbers,with regard to Eq(A3), we conclude that

[TABLE]

is an approximation to Eq.(6) Here $W^{m}(l),m=1,2,N_{rep}$ denotes the $m^{th}$ realization of W(l) out of $N_{rep}$ independently run simulations. In the limit of large t and $N_{rep}$ this approximation approaches an equality, and forms the basis of a computational scheme for the lowest energy of a many particle system with a prescribed symmetry.