Solving the Schrodinger Equation on the Basis of Finite-Difference and Monte-Carlo Approaches

Konstantin Eduardovich Plokhotnikov^1,2

¹Physics Department Lomonosov Moscow State University, Moscow, Russia.
²Financial University under the Government of the Russian Federation, Moscow, Russia.

DOI: 10.4236/jamp.2021.92024                

Abstract

The paper presents a method of numerical solution of the Schrodinger equation, which combines the finite-difference and Monte-Carlo approaches. The resulting method was effective and economical and, to a certain extent, not improved, i.e. optimal. The method itself is formalized as an algorithm for the numerical solution of the Schrodinger equation for a molecule with an arbitrary number of quantum particles. The method is presented and simultaneously illustrated by examples of solving the one-dimensional and multidimensional Schrodinger equation in such problems: linear one-dimensional oscillator, hydrogen atom, ion and hydrogen molecule, water, benzene and metallic hydrogen.

Keywords

Schrodinger Equation, Numerical Methods, Finite Difference and Monte-Carlo Methods

1. Introduction

The paper [1] draws attention to the fact that the Schrodinger equations describing the dynamics of most interesting quantum systems cannot be obtained constructively due to the large dimension of the wave function. Indeed, if the quantum system includes n particles (hereinafter the spin of the particles is not taken into account), then the wave function ψ will depend on 1 + 3n arguments, i.e. $\psi =\psi \left(t,r_1,\cdots ,r_n\right)$, where t is time, $r_1,\cdots ,r_n$ are the radii-vectors of the positions of particles in space, of dimension 3. In quantum chemistry [2], in the Monte-Carlo quantum method [3] [4], this problem is solved in different ways. We will especially note a group of variation methods and approaches, starting with the classical self-consistent Hartree-Fock field method [5] and subsequent variations and upgrades of the type of the functional density method [6].

In this paper, we consider a method for solving the Schrodinger equation, which will allow us to overcome this non-constructivity. The method is based on the idea of representing a set of quantum particles of an arbitrary as a finite set of subparticles (subunits). In terms of positioning a set of subparticles in a space of dimension 3n, a specially prepared random procedure is used, repeated use of this procedure allows to reconstruct the distribution of the average positions of the quantum particles of a molecule in the normal three-dimensional space. The proposed numerical method is efficient and computationally cost-effective and can be attributed to the intersection of finite difference and Monte-Carlo methods.

Some fragments of this work are published in a short article [7]. This article presents a numerical procedure for solving the Schrodinger equation in full.

2. One Quantum Particle in One-Dimensional Space

We write the Schrodinger equation for a single quantum particle in the one dimensional case, i.e.

$i\hslash {\psi }_t=-\frac{{\hslash }^2}{2m}{\psi }_{xx}+U\psi$, (1)

where $\hslash$ is the Planck’s constant, ${\psi }_t=\frac{\partial \psi }{\partial t},{\psi }_{xx}=\frac{{\partial }^2\psi }{\partial x^2}$, x is a space variable, $U=U\left(t,x\right)$ is a function of some potential, m is the mass of the quantum particle, $i^2=-1$.

Write the wave function in (1) as a decomposition into real and imaginary components, i.e. $\psi \left(t,x\right)=u\left(t,x\right)+iv\left(t,x\right)$, then

$\{\begin{array}{l}\hslash u_t=-\frac{{\hslash }^2}{2m}{v}_{xx}+Uv,\\ -\hslash v_t=-\frac{{\hslash }^2}{2m}{u}_{xx}+Uu.\end{array}$ (2)

As you know, the Schrodinger Equation (1) admits the existence of a probability conservation law, which can be written as an equation:

$\frac{\partial }{\partial t}\left({\psi }^{\ast }\psi \right)+\frac{\partial }{\partial x}\left[i\frac{\hslash }{2m}\left(\psi {\psi }_x^{\ast }-{\psi }^{\ast }{\psi }_x\right)\right]=0$, (3)

where the star denotes the complex conjugation operation. In the notation of the real and imaginary parts of the wave function, Equation (3) can be rewritten as:

$\frac{\partial }{\partial t}\left(u^2+v^2\right)+\frac{\hslash }{m}\frac{\partial }{\partial x}\left(u{v}_x-u_{x}v\right)=0$. (3')

Let’s return to the pair of Equation (2), which are equivalent to the original Schrodinger equation in the form (1). According to (2) the quantum particle is described by a pair of functions $u\left(t,x\right)$ and $v\left(t,x\right)$. It is the system of Equation (2) that we will use to construct a numerical method for solving the Schrodinger equation.

We represent the original quantum particle in the form of a set of N subparticles. Let each of the subparticles have a corresponding positioning in space, $x_i,i=1,\cdots ,N$. In this case, each subparticle is described by a pair of functions $u_i\left(t\right),v_i\left(t\right),i=1,\cdots ,N$. One of the important aspects of the numerical method under consideration is the method of replacing the second space derivatives in the system (2) with the following expressions:

$u_{xx}\to -{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{u}_{i,j}}$, $v_{xx}\to -{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{v}_{i,j}}$, (4)

where $u_{i,j}=u_i-u_j$, $v_{i,j}=v_i-v_j$ ; $g_{i,j},i,j=1,\cdots ,N$, some matrix, the form of which we will specify later.

We rewrite the system of Equation (2) taking into account the substitutions (4), then

$\{\begin{array}{l}{\stackrel{\dot{}}{u}}_i=\frac{\hslash }{2m}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{v}_{i,j}}+\frac{1}{\hslash }{U}_i{v}_i,\\ {\stackrel{\dot{}}{v}}_i=-\frac{\hslash }{2m}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{u}_{i,j}}-\frac{1}{\hslash }{U}_i{u}_i,\end{array}$ (5)

where the point is the time derivative, and $U_i=U\left(t,x_i\right)$.

To satisfy the law of conservation of probability (3), (3') by the system of Equation (5), it is necessary:

$\frac{\text{d}}{\text{d}t}{\displaystyle {\sum }_{i=1}^N\left(u_i^2+v_i^2\right)}=0$. (6)

Multiply the first Equation in (5) by $u_i$, and the second, by $v_i$, add both equations and sum by index $i=1,\cdots ,N$, then we get

$\frac{\text{d}}{\text{d}t}{\displaystyle {\sum }_{i=1}^N\left(u_i^2+v_i^2\right)}=\frac{\hslash }{2m}{\displaystyle {\sum }_{i,j=1,j\ne i}^N\left(g_{i,j}-g_{j,i}\right)\left(u_j{v}_i-u_i{v}_j\right)}$. (7)

From (7), it clearly follows that to conserve the probability (in the format (6)), it is sufficient to assume the symmetry of the matrix g, i.e. require that there was the condition $g_{i,j}=g_{j,i},i,j=1,\cdots ,N$.

To select the species of matrix g, we compare the solution of the Schrodinger equation in the format (5) with the finite-difference representation of Equation (1). Without loss of generality, we further assume that there is no potential, i.e. $U\equiv 0$.

Let a certain scale of length $L_0$ be defined, which allows to introduce the characteristic time equal to the value $\frac{m{L}_0^2}{\hslash }$ and reduce the Equation (1) to a dimensionless form. We consider further that the spatial variable x changes in the range $\left[-L,L\right]$, and the value L is expressed in units of the length scale $L_0$.

We introduce a uniform grid $x_i=-L+h\left(i-1\right),i=1,\cdots ,N$, where the grid step $h=\frac{2L}{N-1}$. Using the finite-difference three-point approximation pattern of the second derivative in space, in dimensionless form, Equation (1) without potential can be rewritten as:

$i{\stackrel{\dot{}}{\psi }}_i=-\frac{1}{2}\frac{{\psi }_{i+1}-2{\psi }_i+{\psi }_{i-1}}{h^2}$, (8)

where $i=2,\cdots ,N-1$.

From Equation (8) it follows that as the matrix g should choose a symmetric matrix with non-negative elements of the form:

$g_{i,j}=\{\begin{array}{l}{h}^{-2},\left|i-j\right|=1;\\ 0,\left|i-j\right|\ne 1.\end{array}$ (9)

Note that the finite-difference version of the solution of the Schrodinger Equation (8) does not suit us for two reasons: 1) it is necessary to determine the boundary conditions, which, in general, is not very important when solving the Schrodinger equation in the context of quantum mechanical problems; 2) the finite difference method of the type (8) has no obvious prospect of generalization to the multidimensional case. Note that in (9) the elements of the matrix g that are far from the main diagonal and its two nearest diagonals are zero. Taking into account the above remarks, consider the following representation for the matrix g:

$g_{i,j}=\{\begin{array}{l}{\rho }_1{h}^{-2}{\text{e}}^{-{\rho }_2\left|i-j\right|},i\ne j;\\ 0,i=j;\end{array}$ (10)

where ${\rho }_1,{\rho }_2$, some non-negative parameters. Representation (10) defines a symmetric matrix with non-negative elements.

We rewrite the system of Equation (5) in dimensionless form and in the absence of potential, then

$\{\begin{array}{c}{\stackrel{\dot{}}{u}}_i=\frac{1}{2}{\displaystyle {\sum }_j{g}_{i,j}{v}_{i,j}},\\ {\stackrel{\dot{}}{v}}_i=-\frac{1}{2}{\displaystyle {\sum }_j{g}_{i,j}{u}_{i,j}.}\end{array}$ (11)

We introduce a vector column $y={\left(u_1,\cdots ,u_N,v_1,\cdots ,v_N\right)}^{\text{T}}$ and a block matrix $G=\left|\begin{array}{cc}0\left(N\times N\right)& \frac{1}{2}\left(g_2-g\right)\\ -\frac{1}{2}\left(g_2-g\right)& 0\left(N\times N\right)\end{array}\right|$, where g, matrix (10) and the $g_2$ matrix is diagonal, with the sums on the diagonal $g_{2,i,i}={\displaystyle {\sum }_{j,j\ne i}{g}_{i,j}},i=1,\cdots ,N$. Note that the unknown vector y has dimension 2N, and the matrix G is $2N\times 2N$ in size. Taking into account the vector y and the matrix G, the system of Equation (11) can be rewritten in the following compact form:

$\stackrel{\dot{}}{y}=Gy$ . (12)

Compare the solutions of Equations (8) and (12). To do this, substitute in (8) the solution in the form ${\psi }_i=\frac{1}{\sqrt{2L}}{\text{e}}^{i\left(-\omega t+k{x}_i\right)}$. We impose periodic boundary conditions on the solutions, ${\psi }_1=\frac{1}{\sqrt{2L}}{\text{e}}^{i\left(-\omega t-kL\right)}={\psi }_N=\frac{1}{\sqrt{2L}}{\text{e}}^{i\left(-\omega t+kL\right)}$, then we find $k=\frac{\pi l}{L},l=0,\pm 1,\pm 2,\cdots$ and after substitution in (8), we obtain an expression for the energy:

${\omega }_l=\frac{1}{h^2}\left[1-\cos \left(\frac{\pi h}{L}l\right)\right],l=0,\pm 1,\pm 2,\cdots$ (13)

Similarly, we consider the solution of the system of Equation (12) with matrix (10). To do this, we present the solution of Equation (12) as $y=z{\text{e}}^{i\Omega t}$, where z is some constant vector of dimension column $2N\times 1$. In this case, the value of $i\Omega$ is the eigenvalue of the matrix G. Such eigenvalues are in the general case 2N. It is easy to understand that the search for the eigenvalues of the matrix G is reduced to the solution of the equation:

$\det \Vert \frac{1}{2}\left(g_2-g\right)\pm \Omega e_N\Vert =0$, (14)

where $e_N$, unit matrix of size $N\times N$. We write the solutions of Equation (14) as a set $\pm {\Omega }_1,\cdots ,\pm {\Omega }_N$. Choose the top plus sign. In this case, we obtain a set of frequencies ${\Omega }_1,\cdots ,{\Omega }_N$, which can be compared with the set (13), when $l=0,\cdots ,N-1$.

To ensure coherence of the pair of energy spectra ${\omega }_0,\cdots ,{\omega }_{N-1}$ and ${\Omega }_1,\cdots ,{\Omega }_N$, it is necessary to select the appropriate values of the parameters h and L, which are included in Equation (13). Let’s call these values approximating and denote: $h_{app}$ and $L_{app}$. From the above values, one can also find the number of nodes $N_{app}$ of the finite-difference scheme of format (8), which is equivalent to our numerical scheme from the point of view of the maximum proximity of the spectra. The parameters $h_{app}$ and $L_{app}$ will be found by comparing the spectra, i.e. by minimizing the total comparison error of the type:

$S=S\left(a,b\right)={\displaystyle {\sum }_{i=1}^N{\left\{a-a\cos \left[b\left(i-1\right)\right]-{\Omega }_i\right\}}^2}$. (15)

Let function (15) reach a minimum when $a=a_{opt}$ and $b=b_{opt}$, then taking into account (13) we find $h_{app}=\frac{1}{\sqrt{a_{opt}}}$, $L_{app}=\frac{\pi }{b_{opt}\sqrt{a_{opt}}}$, $N_{app}=1+\frac{2L_{app}}{h_{app}}$. Figure 1 shows two graphs of the spectral decompositions of the pair of schemes, one of them is studied and the other is a finite difference, obtained at N = 10 (Figure 1(a)) and N = 100 (Figure 1(b)), respectively.

Note that the proximity of the spectra strongly depends on the coefficients ${\rho }_1,{\rho }_2$. In the calculations, the results of which are shown in Figure 1, the optimal values of the coefficients ${\rho }_{1,opt},{\rho }_{2,opt}$ were found. Table 1 shows the values of the parameters $h_{app},L_{app},N_{app}$ in calculations with different numbers of grid nodes N in our scheme. In all three cases, both spectra were close. From Table 1, it can be seen that, roughly speaking, the equivalent finite-difference scheme has nodes twice, and the localization segment is five times larger than the scheme under study.

Let us proceed to the numerical solution of the system of Equations (10)-(12). We solve the Cauchy problem by choosing the initial values of the functions $u_i\left(0\right),v_i\left(0\right),i=1,\cdots ,N$ random, provided that ${\displaystyle {\sum }_{i=1}^N\left[u_i^2\left(0\right)+v_i^2\left(0\right)\right]}=1$. Since the systems of Equations (11), (12) preserve the probability, so far the equation ${\displaystyle {\sum }_{i=1}^N\left[u_i^2\left(t\right)+v_i^2\left(t\right)\right]}=1$ will take place at all subsequent moments of the time ( $t>0$ ). The last equation means that we have a discrete random variable X, the position of the quantum particle, which takes the values $x_1,\cdots ,x_N$ with probabilities $w_1=u_1^2+v_1^2$, $\cdots$, $w_N=u_N^2+v_N^2$. You can also determine the average position of a quantum particle $\bar{X}$ and its dispersion $D_X$ according to classical formulas:

$\bar{X}={\displaystyle {\sum }_{i=1}^N{w}_i{x}_i}$, $D_X={\displaystyle {\sum }_{i=1}^N{w}_i{\left(x_i-\bar{X}\right)}^2}$. (16)

Since the probabilities in the set $w_1\left(t\right),\cdots ,w_N\left(t\right)$ depend on time, the average position $\bar{X}=\bar{X}\left(t\right)$ and the dispersion $D_X=D_X\left(t\right)$ of a quantum particle generally depend on time.

(a)(b)

Figure 1. (a) Comparison of the spectra of the studied and finite difference schemes for N = 10. (b) Comparison of the spectra of the studied and finite difference schemes for N = 100.

Table 1. List of parameter values when comparing spectra of the studied and finite-difference scheme.

On the left graph Figure 2 an example of solving the system of Equations (11), (12) is given in terms of dynamics of the components of the wave function $u_{i^{\prime }}\left(t\right)$ and $v_{i^{\prime }}\left(t\right)$ of a random i’-th subunit of a quantum particle placed within the interval $\left[-L,L\right]$. The right graph in Figure 2 shows the probability distribution of finding a quantum particle $w_i\left(t\right),i=1,\cdots ,N$ in space at some time point t. Other design parameters were chosen as follows:

$L=1,N={10}^2,h=\frac{2}{99},{\rho }_1=25,{\rho }_2=5$. From the right graph in Figure 2, it follows that at each instant of time the quantum particle is approximately uniformly “smeared” over the entire localization segment $\left[-L,L\right]$.

Figure 3 shows two dynamics: the average position of a quantum particle (left graph) and the dispersion of the position (right graph). Graphs in Figure 3 constructed according to the formulas (16) after finding solutions to the system of Equations (11), (12) with parameters having the following values:

$L=1,N={10}^2,h=\frac{2}{99},{\rho }_1=25,{\rho }_2=5$. According to the left graph in Figure 3 the average position of the quantum particle $\bar{X}$ vibrates, and the dispersion of vibrations $D_X$ also experiences vibrations.

3. One Quantum Particle in Three-Dimensional Space

We write the Schrodinger equation for a free quantum particle in a dimensionless form in three-dimensional space with coordinates x, y, z then

$i{\psi }_t=-\frac{1}{2}\left({\psi }_{xx}+{\psi }_{yy}+{\psi }_{zz}\right)$. (17)

In space we will define a cube ${\left[-L,L\right]}^3$ with the center at the origin of coordinates and place a quantum particle in it. We introduce finite-difference grids for each coordinate, considering the grid steps for each coordinate to be the same and equal to h, then $x_{i_1}=-L+h\left(i_1-1\right),i_1=1,\cdots ,N,N=1+\frac{2L}{h}$. Similarly, for the other two variables: $y_{i_2}=-L+h\left(i_2-1\right)$, $z_{i_3}=-L+h\left(i_3-1\right)$ ; $i_2,i_3=1,\cdots ,N$. In total, thus, in our grid that fills the selected cube, there are $N^3$ nodes.

Replace in (17) the second derivatives with the usual finite-difference expressions according to the three-point pattern in the same way as it was done in (8). In this case, the finite-difference approximation of Equation (17) is written as:

$i{\stackrel{\dot{}}{\psi }}_{i_1,i_2,i_3}=-\frac{{\psi }_{i_1+1}-2{\psi }_{i_1}+{\psi }_{i_1-1}}{h^2}-\frac{{\psi }_{i_2+1}-2{\psi }_{i_2}+{\psi }_{i_2-1}}{h^2}-\frac{{\psi }_{i_3+1}-2{\psi }_{i_3}+{\psi }_{i_3-1}}{h^2}$. (18)

Figure 2. The dynamics of the components of the wave function of a random subparticle (left schedule); probability distribution of particle localization in space at some point in time (right graph).

Figure 3. Dynamics of the average position of a quantum particle (left graph); the dynamics of the dispersion of the position of a quantum particle (right graph).

Looking for a solution to Equation (18) in the form:

${\psi }_{i_1,i_2,i_3}=\frac{1}{{\left(2L\right)}^{3/2}}{\text{e}}^{i\left(-\omega t+k_1{x}_{i_1}+k_2{y}_{i_2}+k_3{z}_{i_3}\right)}$. (19)

After substitution (19) in (18) we find the condition that (19) is the solution of Equation (18), namely $\omega =\frac{1}{h^2}\left(3-\cos k_{1}h-\cos k_{2}h-\cos k_{3}h\right)$. Let periodic boundary conditions be set on opposite faces of a three-dimensional cube ${\left[-L,L\right]}^3$, then the wave vector $k=\left(k_1,k_2,k_3\right)$ will take a discrete set of values, i.e. $k=\frac{\pi }{L}\left(j_1,j_2,j_3\right)$, $j_1,j_2,j_3=0,\pm 1,\cdots$ In this case, the dependence $\omega =\omega \left(k\right)$ can be rewritten as:

${\omega }_{j_1,j_2,j_3}=\frac{1}{h^2}\left(3-\cos \frac{\pi h}{L}{j}_1-\cos \frac{\pi h}{L}{j}_2-\cos \frac{\pi h}{L}{j}_3\right)$, (20)

where $j_1,j_2,j_3=0,\pm 1,\cdots$

We generalize the expression for the matrix g in (10) to the three-dimensional case. To do this, we introduce multicomponent numbers $i=\left(i_1,i_2,i_3\right)$, $j=\left(j_1,j_2,j_3\right)$, where $i_1,\cdots ,j_1,\cdots =1,\cdots ,N$. Suppose that

$g_{i,j}=\{\begin{array}{l}{\rho }_1{h}^{-2}{\text{e}}^{-{\rho }_2\left(\left|i_1-j_1\right|+\left|i_2-j_2\right|+\left|i_3-j_3\right|\right)},i\ne j;\\ 0,i=j.\end{array}$ (21)

Note that the finite-difference representation of (18) can be interpreted in the sense that the Laplace operator $\Delta =\frac{\partial }{\partial x^2}+\frac{\partial }{\partial y^2}+\frac{\partial }{\partial z^2}$ is approximated on a seven-point template. As in the one-dimensional case (4), replace the Laplace operator in the Schrodinger Equation (17) with the expression $-{\displaystyle {\sum }_{j,j\ne i}{g}_{i,j}{\psi }_j}$, where the summation extends over all index values $j=\left(j_1,j_2,j_3\right),j_1,j_2,j_3=1,\cdots ,N$, except $i$. After substituting the representation $\psi =u+iv$ into Equation (17), we obtain a system (11) where the indices $i$ and $j$ are multicomponent indices $i=\left(i_1,i_2,i_3\right)$, $j=\left(j_1,j_2,j_3\right)$.

It is possible to renumber the subparticles not by three subindexes, but by one. Let be $i=\left(i_1,i_2,i_3\right)$, $i_1,i_2,i_3=1,\cdots ,N$. We introduce a single (one-component) numbering i of all grid nodes in a three-dimensional cube using the formula:

$i=i_1+N\left(i_2-1\right)+N^2\left(i_3-1\right)=1,\cdots ,N^3$. (22)

After determining the single index of subparticles by the formula (22), the Schrodinger equation for the dynamics of a quantum particle in a three dimensional space (without interaction potential) can be reduced to the equations in the format (11), (12).

Let us compare the frequency spectrum of the finite-difference scheme (18)-(20) with the spectrum of the studied scheme of Equations (11), (12), when the transition from multi-index to monoindex is taken into account. Taking into account (22), solving the eigenvalue problem (14), we find the set of frequencies ${\Omega }_1,\cdots ,{\Omega }_{N^3}$. Compare the found set with the ordered set ${\omega }_1,\cdots ,{\omega }_{N^3}$, obtained from (20) according to the following procedure. We assume that $j_1,j_2,j_3=0,1,\cdots ,N-1$ then we introduce a single frequency numbering according to the algorithm of formula (22), i.e. $j=j_1+1+N{j}_2+N^2{j}_3=1,\cdots ,N^3$, then sort the obtained values in ascending order.

To ensure the proximity of a pair of spectra, the problem of minimizing the functional S was solved:

$S=S\left(a,b\right)={\displaystyle {\sum }_{i=1}^{N^3}{\left[{\omega }_i\left(a,b\right)-{\Omega }_i\right]}^2}$, (23)

where before of the numbering and sorting of frequencies (20) it was believed that ${\omega }_{j_1,j_2,j_3}=a\left(3-\cos b{j}_1-\cos b{j}_2-\cos b{j}_3\right)$. After obtaining the optimal values of parameters $a=a_{opt}$ и $b=b_{opt}$, the parameters $h_{app}=\frac{1}{\sqrt{a_{opt}}}$, $L_{app}=\frac{\pi }{b_{opt}\sqrt{a_{opt}}}$, $N_{app}=1+\frac{2L_{app}}{h_{app}}$ of the equivalent finite-difference scheme (18)-(20) were found.

Figure 4 shows several samples comparing a pair of frequency spectra for the three cases when $N=7,10,15$, the other parameters took the following values: $L=1,h=\frac{2}{N-1},{\rho }_1=25,{\rho }_2=5$. Analysis of the graphs in Figure 4 indicates a good correspondence of the frequency spectrum of the studied and equivalent finite difference schemes.

Table 2 summarizes the parameters of the studied and equivalent finite-difference scheme, the spectra of which are compared in Figure 4. Note that the parameters $h_{app}$ and $L_{app}$ were calculated by minimizing the functional (23).

4. Linear Oscillator in One-Dimensional Case

Let us write the potential function of a one-dimensional oscillator in the form $U\left(x\right)=\frac{1}{2}m{\omega }^2{x}^2$, where ω is the oscillation frequency of the oscillator. Introducing the characteristic values of time and length ${\omega }^{-1}$ and ${\left(\frac{\hslash }{m\omega }\right)}^{1/2}$, we rewrite the corresponding Schrodinger equation in dimensionless form. Based on the notation of the system of Equation (2), we write

$\{\begin{array}{l}{u}_t=-\frac{1}{2}{v}_{xx}+\frac{1}{2}{x}^{2}v,\\ v_t=\frac{1}{2}{u}_{xx}-\frac{1}{2}{x}^{2}u.\end{array}$ (24)

We choose a certain interval $\left[-L,L\right]$ of integration over the space of the system of Equation (24). We define at this interval generally speaking uneven grid with nodes: $-L=x_1<x_2<\cdots <x_{N-1}<x_N=L$. In accordance with the studied solution scheme, we rewrite Equation (24) as a set of the following differential equations:

$\{\begin{array}{l}{\stackrel{\dot{}}{u}}_i=\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{v}_{i,j}}+\frac{1}{2}{x}_i^2{v}_i,\\ {\stackrel{\dot{}}{v}}_i=-\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{u}_{i,j}}-\frac{1}{2}{x}_i^2{u}_i,\end{array}$ (25)

(a)(b)(c)

Figure 4. (a) Comparison of the frequency spectra of the studied and the finite difference schemes for N = 4; (b) Comparison of the frequency spectra of the studied and the finite difference schemes for N = 10; (c) Comparison of the frequency spectra of the studied and the finite difference schemes for N = 15.

Table 2. List of parameter values when comparing spectra of the studied and finite-difference scheme.

where $i=1,\cdots ,N$. The expression for the matrix $g_{i,j},i,j=1,\cdots ,N$ is taken from (10).

If you enter a column vector $y={\left(u_1,\cdots ,u_N,v_1,\cdots ,v_N\right)}^{\text{T}}$, then the system of Equation (25) can be rewritten in a form similar to (12), where $G=\left|\begin{array}{cc}0& Q\\ -Q& 0\end{array}\right|$, $Q=Q\left(N\times N\right)=\frac{1}{2}\left(g_2-g\right)$, $g_2=diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}+x_i^2\right)$ is the diagonal matrix of size $N\times N$, where on the diagonal are line-by-line sums of the elements of the matrix g plus double the value of the potential.

Find solutions to the linear system of Equation (12). To do this, we will look for a solution in the form: $y={\text{e}}^{\delta t}z$, where z is some constant column vector of size $2N\times 1$. Let vector z consist of two subvectors $a={\left(a_1,\cdots ,a_N\right)}^{\text{T}}$, $b={\left(b_1,\cdots ,b_N\right)}^{\text{T}}$, i.e. $z=\left|\begin{array}{c}a\\ b\end{array}\right|$. Substitute the selected solution into Equation (12), then we obtain the following problem for finding the eigenvalues (δ) and vectors (z) of the matrix G:

$\left|\begin{array}{cc}0& Q\\ -Q& 0\end{array}\right|\cdot \left|\begin{array}{c}a\\ b\end{array}\right|=\delta \left|\begin{array}{c}a\\ b\end{array}\right|$. (26)

The problem (26) is reduced to a simpler problem of finding eigenvalues (Ω) and vectors ( $c\left(N\times 1\right)$ ) of the matrix Q, i.e. $Qc=\Omega c$ with

$\delta =\pm i\Omega$, $\left|\begin{array}{c}a\\ b\end{array}\right|=\left|\begin{array}{c}c\\ \pm ic\end{array}\right|$. (27)

According to (27), there are 2N eigenvalues and vectors of the problem (26), since the solution of the problem of finding eigenvalues and vectors of the matrix Q implies N solutions.

Choose some eigenvalue Ω and the corresponding eigenvector c and write the real particular solution $y_{\Omega }\left(t\right)$ of Equation (12), i.e.

$y_{\Omega }=\left|\begin{array}{c}\left(C_1\cos \Omega t+C_2\sin \Omega t\right)c\\ \left(-C_1\sin \Omega t+C_2\cos \Omega t\right)c\end{array}\right|$, (28)

where $C_1,C_2$ arbitrary real constants.

Let us compare some approximate solutions (28) of the considered scheme with exact solutions for the quantum oscillator.

As is known [8], the energy of the quantum oscillator (in dimensionless form) is found by the formula $E_n=\frac{1}{2}+n,n=0,1,\cdots$ we write out the known first two exact wave functions:

${\psi }_0=\frac{1}{{\pi }^{1/4}}{\text{e}}^{-\frac{1}{2}it-\frac{1}{2}{x}^2}$, ${\psi }_1=\frac{\sqrt{2}}{{\pi }^{1/4}}x{\text{e}}^{-\frac{3}{2}it-\frac{1}{2}{x}^2}$. (29)

Taking into account (29), we find the probability density $W_0\left(x\right)={\psi }_0^{\ast }{\psi }_0$, $W_1\left(x\right)={\psi }_1^{\ast }{\psi }_1$ and compare it with the probability density distributions of solutions (28) for the first two eigenvalues Ω₁ and Ω₂, and that eigenvalues are ordered in ascending order. We denote the corresponding probability distribution, constructed using the solutions to (28), in the form: $w_{1,i},i=1,\cdots ,N$ and $w_{2,i},i=1,\cdots ,N$ for the eigenvalues Ω₁ and Ω₂ respectively.

Figure 5 shows a comparison of the approximate (solid line) and the exact solution (dashed line). In Figure 5(a), the distribution densities $w_{1,i},i=1,\cdots ,N$ and $W_0\left(x\right)$ were compared, and in Figure 5(b), $w_{2,i},i=1,\cdots ,N$ and $W_1\left(x\right)$. In the calculations of approximate solutions, the following parameter values were used: $L=3,N=2\times {10}^3,h=\frac{6}{1999},{\rho }_1=25,{\rho }_2=3.35$. Approximate solutions are constructed according to (28), with the values of $t,C_1,C_2$ chosen random.

The results of the comparison of the approximate and exact solutions presented in Figure 5 look quite satisfactory. It should be noted that a noticeable non-smooth nature of the approximate curves is due to the fact that the used grid $\left\{x_1,\cdots ,x_N\right\}$ was uneven. The uneven grid was deliberately chosen to ensure the effectiveness of this method of calculation. The advantages of the method are expressed in the absence of too strict requirements for the choice of grid uniformity, both in one-dimensional and multi-dimensional cases.

5. Hydrogen Atom

We construct an algorithm for calculating the dynamics of a pair of quantum particles, a proton and an electron, interacting according to the Coulomb law.

Let $m_1$ and $m_2$ be the masses of proton and electron. We write the Schrodinger equation describing the dynamics of the hydrogen atom as:

$i\hslash {\psi }_t=-\frac{{\hslash }^2}{2m_1}{\Delta }_1\psi -\frac{{\hslash }^2}{2m_2}{\Delta }_2\psi -\frac{e^2}{\left|r_1-r_2\right|}\psi$, (30)

where $\psi =\psi \left(t,r_1,r_2\right)$, wave function of the hydrogen atom, $r_1,r_2$, the spatial positions of the proton and electron, respectively, e, the value of the charge. Laplace operators for each of the quantum particles have the form:

${\Delta }_1=\frac{{\partial }^2}{\partial x_1^2}+\frac{{\partial }^2}{\partial y_1^2}+\frac{{\partial }^2}{\partial z_1^2}$, ${\Delta }_2=\frac{{\partial }^2}{\partial x_2^2}+\frac{{\partial }^2}{\partial y_2^2}+\frac{{\partial }^2}{\partial z_2^2}$.

We rewrite Equation (30) in dimensionless form. To do this, we introduce the characteristic length, mass, time and energy: $r_B=\frac{{\hslash }^2}{m_2{e}^2}\cong 5.2918\times {10}^{-9}\text{cm}$ ( $r_B$ Bohr radius), $m_2=9.1093\times {10}^{-28}\text{g}$, $\frac{{\hslash }^3}{m_2{e}^4}=2.4189\times {10}^{-17}\sec$, $\frac{m_2{e}^4}{{\hslash }^2}=4.3597\times {10}^{-11}\text{erg}$ respectively, then

$i{\psi }_t=-\frac{\mu }{2}{\Delta }_1\psi -\frac{1}{2}{\Delta }_2\psi -\frac{1}{\left|r_1-r_2\right|}\psi$, (31)

where $\mu =\frac{m_2}{m_1}\cong 5.4462\times {10}^{-4}$, the ratio of the mass of the electron to the mass of the proton.

Let us replace the variables $r_1\to \sqrt{\mu }{r}_1$, move on to the real and imaginary parts of the wave function by representing $\psi =u+iv$, then Equation (31) can be rewritten as:

(a)(b)

Figure 5. (a) Comparison of the approximate distribution $w_{1,i},i=1,\cdots ,N$ and exact $W_0\left(x\right)$ ; (b) Comparison of the approximate distribution $w_{2,i},i=1,\cdots ,N$ and exact $W_1\left(x\right)$.

$\{\begin{array}{l}{u}_t=-\frac{1}{2}\Delta v-\frac{1}{\left|\sqrt{\mu }{r}_1-r_2\right|}v,\\ v_t=\frac{1}{2}\Delta u+\frac{1}{\left|\sqrt{\mu }{r}_1-r_2\right|}u,\end{array}$ (32)

where $\Delta ={\Delta }_1+{\Delta }_2$, single Laplace operator for a pair of quantum particles.

Let’s define N points with coordinates $r_i=\left(r_{1,i},r_{2,i}\right),i=1,\cdots ,N$. The selected points have six coordinates, the first three of which characterize the first quantum particle, and the second three, the second quantum particle.

Following the previous sections in the development of a numerical algorithm, we replace the Laplace operators in the system of Equation (32) with expressions: $\Delta u\to -{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{u}_{i,j}}$, $u_{i,j}=u_i-u_j$ ; $\Delta v\to -{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{v}_{i,j}}$, $v_{i,j}=v_i-v_j$. After these substitutions, the system of Equation (32) will be rewritten as a system of linear ordinary differential equations of the form:

$\{\begin{array}{l}{\stackrel{\dot{}}{u}}_t=\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{v}_{i,j}}-\frac{v_i}{\left|\sqrt{\mu }{r}_{1,i}-r_{2,i}\right|},\\ {\stackrel{\dot{}}{v}}_t=-\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{u}_{i,j}}+\frac{u_i}{\left|\sqrt{\mu }{r}_{1,i}-r_{2,i}\right|},\end{array}$ (33)

where $i=1,\cdots ,N$.

To complete the algorithm for calculating the hydrogen atom according to the scheme (33), it is necessary to construct a matrix $g_{i,j},i,j=1,\cdots ,N$. Considering the examples of the construction of this matrix in (9), (10), (21), we choose the following type of matrix:

$g_{i,j}=\{\begin{array}{l}{\rho }_1{h}^{-2}{\text{e}}^{-{\rho }_2\left|r_i-r_j\right|/h},i\ne j;\\ 0,i=j;\end{array}$ (34)

where $i,j=1,\cdots ,N$ ; ${\rho }_1,{\rho }_2$, some non-negative constants, h, the step of the equivalent uniform grid, which we define below.

Taking into account the definition of $r_i=\left(r_{1,i},r_{2,i}\right),i=1,\cdots ,N$, we write the module in (34) and return to the original $r_1$ by replacing $r_1\to \frac{1}{\sqrt{\mu }}{r}_1$, then we get:

$g_{i,j}=\{\begin{array}{l}{\rho }_1{h}^{-2}{\text{e}}^{-\frac{{\rho }_2}{h}\sqrt{\frac{1}{\mu }{\left(r_{1,i}-r_{1,j}\right)}^2+{\left(r_{2,i}-r_{2,j}\right)}^2}},i\ne j;\\ 0,i=j.\end{array}$ (35)

In order for the two terms of the radicand in (35) to be the same in order of magnitude, we require the following representations for the positions of the first and second quantum particles:

$r_{1,i}=a_1+\sqrt{\mu }\left(-L+2L{\xi }_{1,i}\right),r_{2,i}=a_2-L+2L{\xi }_{2,i}$, (36)

where $i=1,\cdots ,N$ ; $a_1,a_2$, some constant vectors, which we will interpret later as the scattering centers of possible positions of proton and electron. We assume that the vectors ${\xi }_{1,i},{\xi }_{2,i}$ in (36) have independent uniformly random coordinate values from the interval [0, 1]. As the selection of step of the equivalent uniform grid, we assume that the cube ${\left(2L\right)}^3$ is divided into N subcubes, of volume h³, then $h=\frac{2L}{N^{1/3}}$.

The representation (36), in which independent uniformly random variables ${\xi }_{1,i},{\xi }_{2,i},i=1,\cdots ,N$ appear, means that the considered method of numerical solution of the Schrodinger equation becomes stochastic and can be referred to the class of Monte-Carlo methods. In what sense will become clear from further constructions.

Taking into account the considerations set out in (34)-(36), after the inverse replacement of $r_1\to \frac{1}{\sqrt{\mu }}{r}_1$, the system of Equation (33) can be rewritten as:

$\{\begin{array}{l}{\stackrel{\dot{}}{u}}_t=\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{v}_{i,j}}-\frac{v_i}{\left|r_{1,i}-r_{2,i}\right|},\\ {\stackrel{\dot{}}{v}}_t=-\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{u}_{i,j}}+\frac{u_i}{\left|r_{1,i}-r_{2,i}\right|},\end{array}$ (37)

where $i=1,\cdots ,N$.

The transition from Equation (33) to Equation (37) made it possible to naturally take into account the different mass of quantum particles, as well as to determine their scattering centers, while, as it will be clear further, the procedure described is universal and applicable to systems with many particles.

Since the system of Equation (37) is linear, we find the eigenvalues of the matrix of the right side, which act as the energy levels of the hydrogen atom. As noted in the previous section, the matrix of the right-hand side of a system of equations similar to (37) has 2N purely imaginary eigenvalues, which are denoted as $\pm i{\Omega }_1,\cdots ,\pm i{\Omega }_N$.

Taking into account the discrete nature of the energy spectrum of a hydrogen atom, which, as is known, is characterized by the dependence $E_n=-\frac{1}{2n^2},n=1,2,\cdots$, we compare it with the spectrum of negative imaginary parts of eigenvalues $-{\Omega }_1,\cdots ,-{\Omega }_N$ matrices of the right side of the system of Equations (35)-(37).

Figure 6(a) shows a sample comparison of a pair of spectra, when it was thought that $L=2,N={10}^4,{\rho }_1=0.075,{\rho }_2=1,a_1=\left(0,0,0\right),a_2=\left(0,0,0\right)$. When selecting the parameters, we sought to ensure that most of the spectrum of the matrix of the right part of the system (35)-(37) was negative.

According to Figure 6(a) the spectrum of the design scheme is significantly different from the exact on the edges. Thus, the value $-{\Omega }_1$ is significantly different from the exact value 0.5. In addition, part of the spectrum is often in the positive region, which corresponds not to the hydrogen atom, but to a pair of individual quantum particles, a proton and an electron. This, in general, is natural, since space in our numerical model appears in discrete form. In consequence of this circumstance, the spectra of the continuous and discrete models at the edges do not necessarily coincide.

Consider the solutions of the system of Equation (37), the expressions for which are given in (28). In the formula (28) Ω and c, eigenvalue and eigenvector of the matrix $Q=\frac{1}{2}\left(g_2-g\right)+diag\left(U_i\right)$, $g_2=diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}\right)$, $U_i=-\frac{1}{\left|r_{1,i}-r_{2,i}\right|}$. Since the matrix Q has dimensions $N\times N$, there are generally N eigenvalues ${\Omega }_1,\cdots ,{\Omega }_N$ and vectors $c_1,\cdots ,c_N$.

(a)(b)

Figure 6. (a) Comparison of the spectrum of the studied numerical scheme (35)-(37) and the exact energy distribution; (b) The dependence of the distribution of the number of localization points in the space of a hydrogen atom on the energy.

Assuming that the eigenvectors are normalized by one, we can find the probability of localization of $w_i,i=1,\cdots ,N$ of the hydrogen atom in one of the points $r_i=\left(r_{1,i},r_{2,i}\right),i=1,\cdots ,N$. To do this, you need to normalize the length of the vector (28) by one, what is provided by elementary transformation: $y_{\Omega }\to y_{\Omega }/\sqrt{C_1^2+C_2^2}$, resulting in ${\Vert y_{\Omega }\Vert }^2={\Vert c\Vert }^2=1$, where $\Vert \mathrm{...}\Vert$, the norm of the vector. Thus, for any of the eigenvalues of (28), $y_{{\Omega }_i}$, the probability $w_j$ of the localization of the hydrogen atom at the point $r_j$ coincides with the square of the coordinate of the eigenvector, $c_i$ i.e. $w_j=c_{i,j}^2,i,j=1,\cdots ,N$.

For a specific example, for each of the N eigenvalues we calculate the number of points $n_k,k=1,\cdots ,N$, in which the hydrogen atom is noticeably localized. As a criterion of noticeable localization, we will choose a certain probability threshold $p_{th}$. In this case, you can write the following formula to calculate the distribution of the number of localization points:

$n_k={\displaystyle {\sum }_{i=1}^N\theta \left(w_i^{\left(k\right)}>p_{th}\right)}$, (38)

where $k=1,\cdots ,N$ ; $w_i^{\left(k\right)}$, the probability of finding a hydrogen atom at the i-th point of space for the k-th eigenvalue, $\theta \left(\mathrm{...}\right)$, a Boolean function that returns one when the inequality is true and, conversely, zero when the inequality is not satisfied.

Taking into account the criterion (38) Figure 6(b) shows an example of calculating the distribution of the number of localization points of the solution in the space of each of N eigenvalues (28) is given. The following calculation parameter values were selected:

$L=2,N={10}^4,{\rho }_1=0.075,{\rho }_2=1,p_{th}=\frac{15}{N}=0.0015$.

In Figure 6(b), attention is drawn to the first and last numbers of its own solutions, which are localized within several or even one point of space. Localization at one point in space means that in this scheme of calculation at the hydrogen atom at low (high) energy levels there has one, far away, nearest neighbor. Note that in the equivalent uniform finite-difference scheme, each point in the six-dimensional space has 12 nearest neighbors. Since the grid in our calculations consists of a set of random points, the number of nearest neighbors of each point with probability close to one is reduced to one. In this circumstance, the role of stochastics is shown, which is included in the calculation scheme according to the algorithm (36) in connection with the random determination of the spatial positions of the proton and electron subparticles.

Taking into account the large variability of the number of localization points of the hydrogen atom depending on the energy (Figure 6(b)), we will conduct a series of calculations (M times), in each of them, new positions of proton and electron particles are randomly selected by formulas (36) and a proper solution with a given energy E is found. In total, this way, M pieces eigensolutions $y_{{\Omega }_{\alpha }}^{\left(\alpha \right)},\alpha =1,\cdots ,M$, will be found, whose eigenvalues selected from the set ${\Omega }_1,\cdots ,{\Omega }_N$, are close to the energy value E, i.e. ${\Omega }_{\alpha }\cong E,\alpha =1,\cdots ,M$.

Figure 7 shows typical examples of localization of the subparticles of a hydrogen atom, the probability of being in which was divided into three categories. The first category includes places where the probability of staying exceeds the threshold $p_{th}$, they are marked with markers in the form of black dots and blue stars. The second and third categories include places, the probabilities of being in which belong to the ranges $\left[\frac{1}{3}{p}_{th},p_{th}\right]$ and $\left[\frac{1}{9}{p}_{th},\frac{1}{3}{p}_{th}\right]$ respectively, these positions are indicated by markers in the form of a circle and rhombs. Markers in the form of stars indicate the position of the proton, and in the form of points, circles and rhombuses, the position of the electron. Other parameters of the computational experiment were chosen as follows:

$L=4,M=50,N={10}^3,{\rho }_1=0.075,{\rho }_2=1,p_{th}=\frac{5}{N}=0.005$.

In Figure 7(a), it is clear that the hydrogen atom itself is rather compact when E = –0.5. The locations of the electron (markers in the form of dots, circles and rhombuses) surround the proton, which is positioned in the center (marker in the form of a blue star) of the electron cluster. After increasing the energy of the hydrogen atom to the value of E = –0.25, the region of electron accumulation sites increases markedly (Figure 7(b)). Finally, when the energy is even closer to zero, E = –0.1, the cloud of electron clusters in Figure 7(c) noticeably moves away from the positions of the proton and fills the entire cube reserved for calculations.

Consider the dynamics of mixed states, i.e. the dynamics of some linear combination of pure solutions (28). To do this, we numerically solve a system of differential Equations (35)-(37), whereas the initial state vector we choose a random vector, $y_0$, normalized by one and prepared according to the algorithm: $y_0={\left(-L+2L{\xi }_1,\cdots ,-L+2L{\xi }_{2N}\right)}^{\text{T}}$, $y_0\to y_0/\Vert y_0\Vert$, ${\xi }_1,\cdots ,{\xi }_{2N}$, uniformly random numbers from the segment [0, 1]. For mixed solutions, the energy of the hydrogen atom is not determined, but it is possible to find the spectrum of energy values. The mixed solution acts as a mixture of pure solutions with energies from the available spectrum.

(a)(b)(c)

Figure 7. (a) Multiple places of localization of a hydrogen atom when E = –0.5; (b) Multiple places of localization of a hydrogen atom when E = –0.25; (c) Multiple places of localization of a hydrogen atom when E = –0.1.

Define kinetic spectra, E_ke, potential, E_pe, and total, E_te, energy by the formulas:

$\begin{array}{l}{E}_{ke}=eig\left[\frac{1}{2}diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}\right)-\frac{1}{2}g\right],\\ E_{pe}=\underset{1\le i\le N}{\text{sort}}{U}_i,\\ E_{te}=eig\left[\frac{1}{2}diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}\right)-\frac{1}{2}g+diag\left(U_i\right)\right],\end{array}$ (39)

where $U_i=-\frac{1}{\left|r_{1,i}-r_{2,i}\right|}$, potential energy, $eig(\dots )$, a procedure that returns a set of eigenvalues of a matrix ranked in ascending order, sort operation of sorting a set of values in ascending order.

Figure 8 on the left shows an example of the graphs of the spectra of all three energies (kinetic, potential, and total). The following set of parameters was selected: $L=2,N=2\times {10}^2,{\rho }_1=0.075,{\rho }_2=1$. The graphs show that the kinetic energy spectrum E_ke is positive everywhere, the potential energy spectrum E_pe on the contrary, is negative everywhere.

The right graph in Figure 8 shows the trajectories of the middle positions of the proton (indicated by a marker in the form of a star) and an electron (indicated by a marker in the form of a dot), while smaller markers describe the beginning of the trajectory, and larger ones, the end of the trajectory. The average positions $R_1$, $R_2$ proton and electron were calculated accordingly by analogy with the formula (16), i.e.

$R_1\left(t\right)={\displaystyle {\sum }_{i=1}^N{r}_{1,i}{w}_i\left(t\right)}$, $R_2\left(t\right)={\displaystyle {\sum }_{i=1}^N{r}_{2,i}{w}_i\left(t\right)}$, (40)

where $w_i\left(t\right)=u_i^2\left(t\right)+v_i^2\left(t\right),i=1,\cdots ,N$, and the functions $u_i,v_i,i=1,\cdots ,N$ are obtained by solving the system of Equation (37) on a time interval [0.2 × 10²]. Note that the average trajectories of the proton and electron were found and displayed on the right graph Figure 8 in M = 10² computational experiments. When switching from one experiment to another according to the algorithm (36), the positions of the proton and electron were randomly changed.

From the right graph in Figure 8 it can be seen that the trajectories of the average position of the electron as a lighter particle have noticeably higher amplitude lines compared to the average trajectories of the proton.

6. Hydrogen Molecule Ion

As part of our calculation scheme, we will study the system of three quantum particles, the ion of a hydrogen molecule. Given the notation of Equation (31), we introduce the second proton as the third particle, then the Schrodinger equation can be rewritten for the hydrogen molecule ion after dimensionalization as:

$i{\psi }_t=-\frac{\mu }{2}{\Delta }_1\psi -\frac{1}{2}{\Delta }_2\psi -\frac{\mu }{2}{\Delta }_3\psi +\left(\frac{1}{\left|r_1-r_3\right|}-\frac{1}{\left|r_1-r_2\right|}-\frac{1}{\left|r_3-r_2\right|}\right)\psi$, (41)

Figure 8. Kinetic, potential and total energy spectra (left graph); dynamics of the average positions of the proton and electron (right graph).

where $\psi =\psi \left(t,r_1,r_2,r_3\right)$, and µ, the ratio of the mass of the electron to the mass of the proton. The expression in brackets in (41) describes all possible Coulomb contributions to the three-particle system.

In Equation (41), we will make the change of variables: $r_1\to \sqrt{\mu }{r}_1$, $r_3\to \sqrt{\mu }{r}_3$, then after switching to the real and imaginary parts of the wave function $\psi =u+iv$, we get

$\{\begin{array}{l}{u}_t=-\frac{1}{2}\Delta v+\left(\frac{1}{\sqrt{\mu }\left|r_1-r_3\right|}-\frac{1}{\left|\sqrt{\mu }{r}_1-r_2\right|}-\frac{1}{\left|\sqrt{\mu }{r}_3-r_2\right|}\right)v,\\ v_t=\frac{1}{2}\Delta u-\left(\frac{1}{\sqrt{\mu }\left|r_1-r_3\right|}-\frac{1}{\left|\sqrt{\mu }{r}_1-r_2\right|}-\frac{1}{\left|\sqrt{\mu }{r}_3-r_2\right|}\right)u,\end{array}$ (42)

where $\Delta ={\Delta }_1+{\Delta }_2+{\Delta }_3$ is the joint Laplace operator for molecular hydrogen ion.

Acting according to the procedure of derivation of Equations (35)-(37), replace in (42) the joint Laplace operator by the corresponding sum, perform the inverse replacement of variables $r_1\to \frac{1}{\sqrt{\mu }}{r}_1$, $r_3\to \frac{1}{\sqrt{\mu }}{r}_3$, then we obtain

$\{\begin{array}{l}{\stackrel{\dot{}}{u}}_i=\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{v}_{i,j}}+\left(\frac{1}{\left|r_{1,i}-r_{3,i}\right|}-\frac{1}{\left|r_{1,i}-r_{2,i}\right|}-\frac{1}{\left|r_{3,i}-r_{2,i}\right|}\right)v_i,\\ {\stackrel{\dot{}}{v}}_i=-\frac{1}{2}{\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}{u}_{i,j}}-\left(\frac{1}{\left|r_{1,i}-r_{3,i}\right|}-\frac{1}{\left|r_{1,i}-r_{2,i}\right|}-\frac{1}{\left|r_{3,i}-r_{2,i}\right|}\right)u_i,\end{array}$ (43)

$g_{i,j}=\{\begin{array}{l}{\rho }_1{h}^{-2}{\text{e}}^{-\frac{{\rho }_2}{h}\sqrt{\frac{1}{\mu }{\left(r_{1,i}-r_{1,j}\right)}^2+{\left(r_{2,i}-r_{2,j}\right)}^2+\frac{1}{\mu }{\left(r_{3,i}-r_{3,j}\right)}^2}},i\ne j;\\ 0,i=j;\end{array}$ (44)

where $i,j=1,\cdots ,N$.

In (43), (44) it is assumed that N points are introduced in a nine-dimensional space with radius vectors $r_i=\left(r_{1,i},r_{2,i},r_{3,i}\right),i=1,\cdots ,N$. In order for the terms under the root in (44) to be of the same order, just as in the case of (36), we put

$\begin{array}{l}{r}_{1,i}=a_1+\sqrt{\mu }\left(-L+2L{\xi }_{1,i}\right),\\ r_{2,i}=a_2-L+2L{\xi }_{2,i},\\ r_{3,i}=a_3+\sqrt{\mu }\left(-L+2L{\xi }_{3,i}\right),\end{array}$ (45)

where $i=1,\cdots ,N$ ; $a_1,a_2,a_3$, some constant vectors that are naturally interpreted as the scattering centers of the positions of two protons and an electron. Vectors ${\xi }_{1,i},{\xi }_{2,i},{\xi }_{3,i}$ in (45) are assumed to have independent uniformly random coordinate values from the segment [0, 1].

It is known [9] that in the ion of a hydrogen molecule a single-electron coupling has an energy of 62 kcal/mol, in the CGS system, it is equal to 4.243 × 10^–12 erg and, after passing to the dimensionless value in our problem, will be 0.0989. Let searched and sorted in ascending order all the eigenvalues

${\Omega }_1,\cdots ,{\Omega }_N$ of the matrix $Q=\frac{1}{2}\left(g_2-g\right)+diag\left(U_i\right)$, $g_2=diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}\right)$, $U_i=\frac{1}{\left|r_{1,i}-r_{3,i}\right|}-\frac{1}{\left|r_{1,i}-r_{2,i}\right|}-\frac{1}{\left|r_{3,i}-r_{2,i}\right|}$. Let’s choose among the set of eigenvalues the value, ${\Omega }_{\alpha }$, which is closest to the energy of the single-electron coupling in the ion of the hydrogen molecule, i.e. the value of 0.0989. In addition to the eigenvalue ${\Omega }_{\alpha }$, we find the eigenvector $c_{\alpha }$ of the matrix Q, which is normalized by one. Carrying out this procedure M times, we find the set $\left\{{\Omega }_{\alpha },c_{\alpha }\right\},\alpha =1,\cdots ,M$, wherein ${\Omega }_{\alpha }\cong -0.0989,\alpha =1,\cdots ,M$. In the previous section, it was already noted that the element square of the eigenvector acts as the probability of the particle’s presence in the selected set of spatial points.

The approximate structure of the ion of the hydrogen molecule is known [10], the protons are spaced at a distance of 1.06 Å, which corresponds to two Bohr radii, i.e. two units in a dimensionless form, with the center of scattering of electronic subparticles located exactly in the middle between the protons. This structural information was sufficient to build an algorithm for generating the entire set of points $r_1,\cdots ,r_N$ according to the procedure (45), when it is considered that, for example, $a_1=\left(-1,0,0\right),a_2=\left(0,0,0\right),a_3=\left(1,0,0\right)$.

Figure 9 shows: an example of calculating the kinetic energy, potential and total energies (left graph) of an ion of a hydrogen molecule, as well as the most likely locations of protons and electron localization of the hydrogen molecule ion (right diagram), when the binding energy is equal to $E\cong -0.0989$. All three energies are shown in the left Figure 9, calculated according to the formulas (39). Other parameter values were selected as follows:

Figure 9. Kinetic, potential and total energy spectra (left graph); the most likely localization of protons and electrons hydrogen molecule ion (right graph).

$L=3,M={10}^2,N={10}^3,{\rho }_1=0.5,{\rho }_2=1$.

On the right graph Figure 9 places of proton and electron localization were marked by markers in the form of rhombuses, pentagrams and points, circle, rhombic, respectively. The probabilities of the electron being in the appropriate places were divided into three categories: 1) $>p_{th}$ (markers in the form of black dots); 2) $\left[\frac{1}{3}{p}_{th},p_{th}\right]$ (markers circles); 3) $\left[\frac{1}{9}{p}_{th},\frac{1}{3}{p}_{th}\right]$ (markers rhombus), where $p_{th}=\frac{25}{N}=0.025$.

The study of the localization of the components of the ion of the hydrogen molecule suggests the following. First, both the proton experiences a “shiver”. Second, the electron “surrounds” both protons.

For control, the kinetic, potential and total energy spectra were found in the absence of variability in the positions of subparticles in space. Instead of formulas (45), it was assumed that $r_{1,i}=a_1,r_{2,i}=a_2,r_{3,i}=a_3,i=1,\cdots ,N$. In this case, it was found that the spectra of the kinetic and total energies have two different values. So when $L=3,N={10}^3,{\rho }_1=0.5,{\rho }_2=1$ it turned out that $E_{ke}=\left\{0,694.4444\right\}$, $E_{pe}=\frac{1}{\left|a_1-a_3\right|}-\frac{1}{\left|a_1-a_2\right|}-\frac{1}{\left|a_3-a_2\right|}=-\frac{3}{2}$, $E_{te}=\left\{-\frac{3}{2},692.9444\right\}$. Thus, when quantum particles are not divided into subparticles and are classically concentrated at separate points, they are in two states: 1) the kinetic energy of the particles is zero, potential and total energy is negative; 2) the kinetic and total energy is positive, i.e. the ion of the hydrogen molecule as a whole does not exist.

7. Hydrogen Molecule

Consider in our computational approach a hydrogen molecule, which consists of four quantum particles.

We note Equations (36), (45), which for the hydrogen atom and the ion of a hydrogen molecule provided the possibility of statistical generation of positions in the space of quantum subparticles. In these formulas, there are scattering centers of subparticles of each of the quantum particles $a_1,a_2,\cdots$, which were chosen on the basis of existing considerations of symmetry and knowledge of the average distances between the quantum particles. In this regard, we consider the distance between protons to be known, which is the distance 2b = 0.7412 Å [9] or in the dimensionless form 2b = 1.4014. Considering the obvious symmetry of the hydrogen molecule with respect to the line passing through the pair of protons, we assume that the electron scattering centers are concentrated in the middle of a straight line connecting the centers of protons.

We introduce designations for the scattering centers of two protons $a_1,a_4$ and two electrons $a_2,a_3$, then we can write: $a_1=\left(-b,0,0\right)$, $a_2=\left(0,0,0\right)$, $a_3=\left(0,0,0\right)$, $a_4=\left(b,0,0\right)$. Note that the diameter of the hydrogen molecule is $d_{{\text{H}}_{\text{2}}}=0.259\text{nm}$ [11] or in the dimensionless form $d_{{\text{H}}_{\text{2}}}=4.9$, which is slightly less than the length of the edge of the cube 2L = 6, which contained the ion of the hydrogen molecule in previous calculations.

For the purposes of numerical modeling, it remains to estimate the energy of the hydrogen molecule. Let us turn to the reference book [12], which presents the dissociation energy of a hydrogen molecule, equal to 103.267 kcal/mol. In our dimensionless quantities it will be 0.1647 units of energy. We believe that before the formation of the hydrogen molecule, hydrogen atoms were in the ground state with energy equals –0.5. In this case, we can assume that the energy of the hydrogen molecule will be $E=-2\times 0.5-0.1647=-1.1647$.

We introduce N points in twelve-dimensional space with radius vectors $r_i=\left(r_{1,i},r_{2,i},r_{3,i},r_{4,i}\right),i=1,\cdots ,N$, which determine the positions of hydrogen molecule subparticles. Taking into account (45), we write the algorithm for generating a set of vectors $r_i,i=1,\cdots ,N$, namely

$\begin{array}{l}{r}_{1,i}=a_1+\sqrt{\mu }\left(-L+2L{\xi }_{1,i}\right),\\ r_{2,i}=a_2+-L+2L{\xi }_{2,i},\\ r_{3,i}=a_3+-L+2L{\xi }_{3,i},\\ r_{4,i}=a_4+\sqrt{\mu }\left(-L+2L{\xi }_{4,i}\right),\end{array}$ (46)

where $i=1,\cdots ,N$ ; vectors ${\xi }_{1,i},{\xi }_{2,i},{\xi }_{3,i},{\xi }_{4,i}$ have independent uniformly random coordinate values from the segment [0, 1].

Define the matrix $Q=\frac{1}{2}\left(g_2-g\right)+diag\left(U_i\right)$, $g_2=diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}\right)$, $U_i=\frac{1}{\left|r_{1,i}-r_{4,i}\right|}+\frac{1}{\left|r_{2,i}-r_{3,i}\right|}-\frac{1}{\left|r_{1,i}-r_{2,i}\right|}-\frac{1}{\left|r_{1,i}-r_{3,i}\right|}-\frac{1}{\left|r_{4,i}-r_{1,i}\right|}-\frac{1}{\left|r_{4,i}-r_{2,i}\right|}$, describingthe interaction potential of protons and electrons in a hydrogen molecule. To find the pure states it is necessary to find the eigenvalues ${\Omega }_1,\cdots ,{\Omega }_N$ of the matrix Q and the corresponding eigenvectors. We will carry out this procedure in two stages. In the beginning, we assume that the kinetic energy is absent. This can be ensured either with ${\rho }_1=0$, or with ${\rho }_2\to \infty$. In the second stage, we choose the values of ${\rho }_1,{\rho }_2$ such that the kinetic and potential energies are connected within the framework of the well-known virial theorem.

Among the set of eigenvalues, choose ${\Omega }_{\alpha }$, the value of which is closest to the binding energy of the hydrogen molecule, i.e. the magnitude –1.1647. In addition to the eigenvalue ${\Omega }_{\alpha }$, we find the eigenvector $c_{\alpha }$ of the matrix Q, which is normalized by one. Conducting this procedure M times we find the set $\left\{{\Omega }_{\alpha },c_{\alpha }\right\},\alpha =1,\cdots ,M$, wherein ${\Omega }_{\alpha }\cong -1.1647,\alpha =1,\cdots ,M$. The element-wise square of the eigenvector acts as the probability of the presence of a particle in a selected set of spatial points.

For Figure 10 shows the result of the calculation of the spatial structure of the hydrogen molecule with account (46), as well as the algorithm, which is tested on the example of calculating the ion of the hydrogen molecule, provided that the kinetic energy, defined in (39), is absent. On the left graph Figure 10 shows the potential energy of the hydrogen molecule, $E_{pe}$.

Figure 10. Potential energy of the hydrogen molecule (left graph); localization of the middle positions of protons and electrons in space (graph in the middle) and in the projection on the plane (x, y) (right graph).

On the middle and right graphs the results of the search for a set of eigenvalues and eigenvectors $\left\{{\Omega }_{\alpha },c_{\alpha }\right\},\alpha =1,\cdots ,M$ such that ${\Omega }_{\alpha }\cong -1.1647,\alpha =1,\cdots ,M=5\times {10}^3$ are shown. The element squares of eigenvalues allowed us to find the places of localization of the average positions of the quantum particles of the hydrogen molecule. Other parameter values were selected as follows: $L=3,N={10}^3$.

The middle and right graphs of Figure 10 show the average positions of protons (markers in the form of pentagrams and hexagrams) and electrons (markers in the form of points and asterisks) of the hydrogen molecule in space (x, y, z) and in the form of a projection on the coordinate plane (x, y). The average positions of quantum particles were calculated using formulas similar to (40), i.e. $R_{k,\alpha }\left(t\right)={\displaystyle {\sum }_{i=1}^N{r}_{k,i}{c}_{\alpha ,i}^2}$, where $k=1,2,3,4$, numbers of quantum particles of a hydrogen molecule, and $\alpha =1,\cdots ,M$.

According to the middle and right graphs Figure 10 it turns out that a cloud of electron subparticles “envelops” both protons, and it is clear that the electron subparticles form something in the form of a three-dimensional figure of rotation.

Note that in terms of calculations, the search for eigenvalues and vectors of the matrix $Q=\frac{1}{2}diag\left(U_i\right)$ is a trivial problem, since the matrix Q in this case is diagonal. The task of finding the eigenvalues and vectors of the matrix Q becomes noticeably more expensive from the computational point of view in the second case, when the kinetic energy of the protons and electrons of the hydrogen molecule is different from zero. In this case, when determining the matrix g, the coefficients ${\rho }_1,{\rho }_2$ must be chosen appropriately.

To select the appropriate values of the coefficients ${\rho }_1,{\rho }_2$ we use the virial theorem known in both classical and quantum mechanics. Applied to our case, the theorem states that the average value of the kinetic energy is equal to minus half of the potential energy for a given value of $E=-1.1647$ of the energy of the hydrogen molecule.

Figure 11 shows the results of a computational experiment for the second case, when the kinetic, potential, and total energies of the hydrogen molecule are taken into account. Other parameter values were selected as follows:

$L=3,N={10}^3,M=1500,{\rho }_1=\frac{3}{N}=0.003,{\rho }_2=0.07$.

On the left graph Figure 11 it is clearly seen that the kinetic energy of the hydrogen molecule is available and it is comparable with the potential energy, and the virial theorem is fulfilled.

Analysis and comparison of Figure 10, Figure 11 shows that electrons (markers in the form of points and asterisks) of the hydrogen molecule condense in the vicinity of a somewhat elongated figure of rotation. If in Figure 10 this figure is only a little like a dumbbell, then in Figure 11 it is quite clearly visible. Note that the figure of rotation in the form of a dumbbell acts as a standard image of a covalent chemical bond.

Note that the topology of the average positions of electrons for the first case (Figure 10) and the second (Figure 11) are noticeably different, although in Figure 10 it is possible to identify the presence of a dumbbell prototype. This means that, to study the general case, when there are both kinetic and potential energies in the molecule, the first case can be used as heuristic calculations when the kinetic energy is considered to be zero. In the latter case, in terms of the amount of computation, the algorithm is greatly simplified.

8. Water Molecule

Consider a water molecule as the next example demonstrating the use of the proposed numerical algorithm for solving the Schrodinger equation. There are 13 quantum particles in a water molecule: an oxygen atom, two protons and ten electrons.

It is known that the angle between the lines directed from the oxygen atom to the protons in the water molecule is approximately 104.45^˚. In addition, the length of these lines is equal to the value 0.9584 Å.

We estimate the dissociation energy $E_{{\text{H}}_{\text{2}}\text{O}}$ of the water molecule on 13 independent quantum particles. It is obvious that $E_{{\text{H}}_{\text{2}}\text{O}}=2E_{\text{H}}+E_{\text{O}}+E_{\text{H}-\text{O}-\text{H}}$, where $E_{\text{H}},E_{\text{O}}$, energy of hydrogen and oxygen atoms formation; $E_{\text{H}-\text{O}-\text{H}}$, energy of water molecule formation. Taking into account section 5, we believe that hydrogen atoms have dimensionless energy equal to –0.5, i.e. $E_{\text{H}}=-0.5$. To estimate the energy of the oxygen atom $E_{\text{O}}$ it is necessary to refer to the experimental data on the energy of the first, second, etc. ionization. It is found [13] that the whole set of ionization energies of the oxygen atom is the following values (in eV): 13.6; 35.1; 54.9; 77.4; 113.9; 138.9; 739.1; 871.1. Summing the last eight values and moving to the dimensionless form, we find $E_{\text{O}}=-75.0826$. Finally, 216.87 kcal/mol, is required to break both covalent bonds in a water molecule [12], which in dimensionless form will be $E_{\text{H}-\text{O}-\text{H}}=-0.3459$. It is obvious that the main contribution to the water molecule gives the energy of formation of the oxygen atom. As a result, by elementary calculation we find the energy of complete dissociation of the water molecule, i.e. $E_{{\text{H}}_{\text{2}}\text{O}}=-77.4285$.

Figure 11. Kinetic, potential and total energy of the hydrogen molecule (left graph); places of localization of the average positions of protons and electrons in space (middle graph) and in the projection on the plane (x, y) (right graph).

We enumerate the quantum particles as follows: oxygen atom (#1); protons (#2, 3); electron (#4), associated with proton #2; electron (#5), associated with proton #3; the remaining 8 electrons (#6-#13) are associated with the oxygen atom.

According to the chosen numbering, we introduce a set of ${\mu }_1,\cdots ,{\mu }_{13}$ ratios of the electron mass to each of the 13 masses of quantum particles of the water molecule. In this case, we get:

${\mu }_1=3.4039\times {10}^{-5};{\mu }_2={\mu }_3=5.4462\times {10}^{-4};{\mu }_4=\cdots ={\mu }_{13}=1$. (47)

We also introduce a set of charges $q_1,\cdots ,q_{13}$ of each of the 13 particles of a water molecule in dimensionless form, then

$q_1=8,q_2=q_3=1,q_4=\cdots =q_{13}=-1$. (48)

Taking into account the spatial position of the particles of the hydrogen molecule, we determine their scattering centers $a_1,\cdots ,a_{13}$. Figure 12(a) shows the positioning in the (x, y) plane of the scattering centers of the oxygen atom and the pair of protons. The electron scattering centers #4, #5 are connected with the centers of scattering of protons, and the centers of scattering of electrons #6-#13 with the scattering center of the oxygen atom. As a result, we have:

$a_1=\left(0,0,0\right);a_2=\left(-b\cos \left(\pi -\gamma \right),b\sin \left(\pi -\gamma \right),0\right);a_3=\left(b,0,0\right)$ ;

$a_4=\left(-b\cos \left(\pi -\gamma \right),b\sin \left(\pi -\gamma \right),0\right);a_5=\left(b,0,0\right)$ ;

$a_6=\cdots =a_{13}=\left(0,0,0\right)$,

where $\gamma =104.45˚$, the angle between the lines going from the oxygen atom to the protons; b = 0.9584 Å = 1.8117, the distance between the oxygen atom and the proton.

We introduce N points in 39-dimensional space with radius vectors $r_i=\left(r_{1,i},\cdots ,r_{13,i}\right),i=1,\cdots ,N$, which determine the positions of water molecule subparticles. Considering (46), we write the algorithm for generating a set of vectors $r_i,i=1,\cdots ,N$, namely

$r_{k,i}=a_k+\sigma 2L\sqrt{{\mu }_k}{\eta }_{k.i}$, (49)

where $k=1,\cdots ,13$ ; $i=1,\cdots ,N$ ; s, some auxiliary non-negative numeric parameter; coordinates of vectors ${\eta }_{k,i},k=1,\cdots ,13;i=1,\cdots ,N$ have independent random values obeying the normal law with mean 0 and standard deviation 1.

Note that according to (49) the random selection of the k-th radius vector obeys the normal law with the mean $a_k$ and standard deviation $\sigma 2L\sqrt{{\mu }_k}$. The choice of the normal distribution instead of the uniform one, as it was in the previous section, is connected with the economy of computing resources, the deficit of which was connected with the provision in the calculations of the low energy value of the complete dissociation of the water molecule, equal to $E_{{\text{H}}_{\text{2}}\text{O}}=-77.4285$.

At the first stage, we find the average positions of all quantum particles of the water molecule in M = 10³ statistical experiments when the kinetic energy is absent. During each statistical experiment, the positions of all thirteen quantum particles are generated according to the algorithm (49). Figure 12(b) shows a typical sample of the potential energy spectrum, $E_{pe}$, of a water molecule.

(a)(b)

Figure 12. (a) Positioning the scattering sites of the oxygen atom and two protons; (b) Potential energy of water molecule.

Potential energy profile in Figure 12(b) is obtained after sorting the set potentials $U_i={\displaystyle {\sum }_{j,k=1,j<k}^{13}\frac{q_j{q}_k}{\left|r_{j,i}-r_{k,i}\right|}},i=1,\cdots ,N$ by ascending order. In addition, the eigenvectors $c_1,\cdots ,c_N$ of the matrix $Q=diag\left(U_i\right)$, were found, which allowed to calculate the average positions of all quantum particles of the water molecule. Figure 13 shows the results of the computational experiment.

Figure 13(a) shows the localization of the middle positions of the oxygen atom (a marker in the form of a circle), protons (a marker in the form of a pentagram) and electrons (markers in the form of points and asterisks) of the water molecule in the space (x, y, z) and in the projection on the plane (x, y). The positioning of individual electrons in the vicinity of oxygen and protons is clearly visible. The same can be seen in Figure 13(b), where the average positions of the electrons are connected by lines.

(a)(b)

Figure 13. (a) Places of localization of the average positions of quantum particles of a water molecule in space (x, y, z) (left graph) and in the projection on the plane (x, y) (right graph); (b) Localization of the middle positions of the electrons of the water molecule connected by lines.

Note that shown in Figure 13 the three-dimensional design of the water molecule corresponds to the concepts of the covalent bond of the oxygen atom with hydrogen atoms. In connection with the manipulation of the electron scattering centers $a_4,\cdots ,a_{13}$, it seems possible to construct other structures for positioning oxygen and hydrogen atoms into a single molecule.

Let us proceed to the general case, when the presence of kinetic energy is taken into account. We take into account its consideration with the virial theorem, i.e. it is believed that the kinetic energy is equal to minus half of the potential. Let us find the spectrum of the kinetic, potential and total energies, as well as the average positions of the oxygen atom, protons and electrons in M statistical tests by playing the positions of quantum particles of a water molecule according to the algorithm (49). To do this, following (44), we define the matrix g by the formula:

$g_{i,j}=\{\begin{array}{l}{\rho }_1{h}^{-2}{\text{e}}^{-\frac{{\rho }_2}{h}\sqrt{{\displaystyle {\sum }_{k=1}^{13}\frac{1}{{\mu }_k}{\left(r_{k,i}-r_{k,j}\right)}^2}}},i\ne j;\\ 0,i=j.\end{array}$ (50)

Let all eigenvalues ${\Omega }_1,\cdots ,{\Omega }_N$ of the matrix $Q=\frac{1}{2}\left(g_2-g\right)+diag\left(U_i\right)$, $g_2=diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}\right)$ be found and ordered by ascending order. Among the set of eigenvalues, choose ${\Omega }_{\alpha }$, the value of which is closest to the binding energy of the water molecule, i.e. the magnitude –77.4285. In addition to the eigenvalue ${\Omega }_{\alpha }$ we find the eigenvector $c_{\alpha }$ of the matrix Q, which is normalized by one. Conducting this procedure M times we find the set $\left\{{\Omega }_{\alpha },c_{\alpha }\right\},\alpha =1,\cdots ,M$, wherein ${\Omega }_{\alpha }\cong -77.4285,\alpha =1,\cdots ,M$.

The left graph of Figure 14 shows the typical profiles of the kinetic, potential and total energies of the water molecule. The graph shows that the desired value of the energy of the complete dissociation of the water molecule is included in the energy spectrum represented by the curve $E_{te}$.

Taking into account (48)-(50), the middle and right graphs in Figure 14 show the results of a typical calculation of the average positions of all 13 quantum particles of a water molecule. The large circle of blue indicates an oxygen atom, the pentagrams indicate the positions of a pair of protons, the black dots indicate the positions of electrons. Other values of the calculation parameters are as follows:

$L=3,N={10}^3,M={10}^3,{\rho }_1=\frac{60}{N}=0.06,{\rho }_2=0.01,\sigma =0.1$.

The middle graph in Figure 14 shows the three-dimensional image of the water molecule in the form of the averaged positions of oxygen, a pair of protons and a halo of electron positions. The right graph in Figure 14 shows the projection of the positions of all quantum particles on the plane (x, y). The increased concentration of electron positioning sites in the vicinity of the oxygen atom is clearly visible. Note that the locations of the oxygen atom and protons experience fluctuations with different amplitudes, remaining in approximately the same position relative to each other.

Figure 14. Energy spectra (left); three-dimensional (center) and two-dimensional image (right, projection on the (x, y) plane) positions of the oxygen atom (marker: blue circle), pairs of protons (marker: pentagram) and electrons (marker: point).

9. Benzene Molecule

Consider the application of the algorithm of numerical solution of the Schrodinger equation to the benzene molecule. In the benzene molecule, whose chemical formula is C₆H₆, 54 quantum particles: 6 carbon atoms, 6 protons and 42 electrons.

Using the example of a benzene molecule we will consider a new feature of the algorithm under study, it is associated with the matrix g, the typical form of which is given in (50). The mentioned feature has already manifested itself in the study of hydrogen and water molecules, where, due to the conditions of the virial theorem, it was necessary to significantly reduce the value of the coefficient ${\rho }_2$. The specified coefficient was chosen to be 0.07 and 0.01 for hydrogen and water molecules, respectively.

Let the algorithm for generating a set of radius vectors $r_i=\left(r_{1,i},\cdots ,r_{54,i}\right),i=1,\cdots ,N$, which determine the positions of the quantum particles of the benzene molecule, be similar to the algorithm of the water molecule (49), and the matrix g is constructed by analogy with the matrix (50). Substitute (49) in (50), then we get applied to the benzene molecule:

$g_{i,j}=\{\begin{array}{l}{\rho }_1{h}^{-2}{\text{e}}^{-\frac{\sigma 2L{\rho }_2}{h}\sqrt{{\displaystyle {\sum }_{k=1}^{54}\left[{\left({\eta }_{x,k,i}-{\eta }_{x,k,j}\right)}^2+{\left({\eta }_{y,k,i}-{\eta }_{y,k,j}\right)}^2+{\left({\eta }_{z,k,i}-{\eta }_{z,k,j}\right)}^2\right]}}},i\ne j;\\ 0,i=j.\end{array}$ (51)

where subindexes x, y, z denote projections of vectors ${\eta }_{k,i},k=1,\cdots ,54;i=1,\cdots ,N$ on the corresponding coordinate axes. Note that in (51) under the sign of the sum there are 54 × 3 = 162 terms, each of which represents the realization of the same random variable $\Phi ={\left({\eta }_{\ast ,\ast ,\ast }-{\eta }_{\ast ,\ast ,\ast }\right)}^2$. Since the number of terms in the sum of random numbers is quite large, we can use the central limit theorem of probability theory.

Find the mathematical expectation MΦ and the variance DΦ of a random variable Φ. Given that all values of ${\eta }_{\ast ,\ast ,\ast }$ are independent and have the standard normal distribution N(0, 1), it is easy to find $M\Phi =2$, $D\Phi =12$. We define the

sample mean ${\bar{\Phi }}_n=\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{\Phi }_k}$, then, according to the central limit theorem, we can write $\sqrt{n}\frac{{\bar{\Phi }}_n-M\Phi }{\sqrt{D\Phi }}=\sqrt{n}\frac{{\bar{\Phi }}_n-2}{12}\to N\left(0,1\right)$ by distribution, when $n\to \infty$.

Let η be some random variable obeying the normal law N(0, 1), then, according to the central limit theorem, we can write ${\bar{\Phi }}_n=2+\frac{12}{\sqrt{n}}\eta \to 2,n\to \infty$, where convergence is understood as convergence in probability. The last passage to the limit in expression (51) allows the sum to be approximately replaced by the constant 162 × 2, in this case it is assumed that n = 162. This is where the new feature of the algorithm under consideration lies, and the above-mentioned substitution will be the more accurate the larger the number of quantum particles in a molecule. We rewrite the expression for the matrix g taking into account this feature, then

$g_{i,j}=\{\begin{array}{l}\epsilon ,i\ne j;\\ 0,i=j;\end{array}$ (52)

where $\epsilon ={\rho }_1{h}^{-2}{\text{e}}^{-\frac{36L\sigma {\rho }_2}{h}}$.

Taking into account (52) the eigenvalue problem for the matrix $Q=\frac{1}{2}\left(g_2-g\right)$, $g_2=diag\left({\displaystyle {\sum }_{j=1,j\ne i}^N{g}_{i,j}}\right)$ can be solved analytically. In this case, the kinetic energy spectrum will have two non-negative values, ${\Omega }_1,{\Omega }_2$, one of which is zero. The corresponding characteristic equation has the form: $\Omega {\left(\Omega -\frac{\epsilon N}{2}\right)}^{N-1}=0$, whence it follows that ${\Omega }_1=0,{\Omega }_2=\frac{\epsilon N}{2}$, and the second eigenvalue has a degeneracy of order N – 1.

Note that the above-discussed feature of the matrix g remains in the case when the random variables of type $\Phi ={\left({\eta }_{\ast ,\ast ,\ast }-{\eta }_{\ast ,\ast ,\ast }\right)}^2$ they are subject to different distributions. Due to the well-known generalized Chebyshev’s theorem, only their independence, the existence of averages and the limitation of variances in the aggregate are required. In this case, convergence in probability is provided by the convergence to the mean of the mathematical expectations of random variables.

Figure 15(a) shows a block diagram of a benzene molecule in the form of a ring, which is well known in chemistry. Figure 15(a) shows the characteristic dimensions of the ring in nanometers (nm) and in dimensionless units (d.u.).

Let us calculate the energy of complete dissociation of the benzene ring, $E_{{\text{C}}_{\text{6}}{\text{H}}_{\text{6}}}$. So, the ionization energies of a carbon atom are (MJ/mol) [12]: 1.09; 2.35; 4.62; 6.22; 37.83; 47.28. Summing up, we find $E_{\text{C}}=-99.39\text{Mj}/\text{mol}=-37.8669\text{d}.\text{u}$. The ionization energy of the hydrogen atom will be $E_{\text{H}}=-0.\text{5}\text{d}.\text{u}$. Finally, the energy of the benzene ring [14] is. As a result, we find the total dissociation energy of the benzene ring:

.

(a)(b)

Figure 15. (a) Appearance of the benzene ring; (b) Kinetic, potential and total energy of the benzene molecule.

Taking into account (49), we write the algorithm for generating the set of vectors $r_i,i=1,\cdots ,N$, namely

$r_{k,i}=a_k+\sigma 2L\sqrt{{\mu }_k}{\eta }_{k,i}$, (53)

where $k=1,\cdots ,54;i=1,\cdots ,N$ ; $a_k,k=1,\cdots ,54$, the centers of the scattering of the quantum particles of the benzene molecule. The coordinates of the vectors ${\eta }_{k,i},k=1,\cdots ,54;i=1,\cdots ,N$ in (53) are assumed to have independent random values with zero expectation and finite variance.

Taking into account (53), we consider several schemes for choosing the electron scattering centers of the benzene molecule. This circumstance is caused by the fact that several schemes of delocalization of electrons of double bond in benzene molecule are considered in chemistry.

Note that to provide the virial theorem, it is sufficient to put that $E_{ke}=\frac{\epsilon N}{2}=-E_{te}=-E_{{\text{C}}_{\text{6}}{\text{H}}_{\text{6}}}=232.3103\text{d}.\text{u}$. Where it follows that $\epsilon =-\frac{2E_{{\text{C}}_{\text{6}}{\text{H}}_{\text{6}}}}{N}$.

Scheme No.1. In the beginning, we consider the case when the electron scattering centers coincide with the corresponding scattering centers of carbon and proton atoms. Figure 15(b) shows an example of a typical calculation of the graphs of the kinetic, potential and total energies of the benzene molecule. Figure 16 shows the localization of the average positions of all quantum particles of the benzene molecule in M = 500 statistical experiments. The left graph shows positioning in space, and the right graph shows positioning in a plane. Carbon atoms, protons and electrons are indicated by markers in the form of pentagrams, large dots and small dots, respectively.

Scheme No.2. The valence electrons of carbon and hydrogen are paired in two. Since only 6 × 4 + 6 × 1 = 30 valence electrons, 15 pairs of bonds can be determined. These bonds are related to C–C, C=C and C–H bonds in the amount of 3, 3 and 6, respectively. The centers of scattering of electron pairs will be placed in the centers of the lines defining the bonds C–C, C=C and C–H. Among the carbon bonds, there are three single C–C and three double C=C. The scattering centers of non-valent electrons (12 pieces) are associated with the scattering centers of the corresponding carbon atoms.

Figure 17 shows the result of calculating the localization sites of the average positions of all the quantum particles of the benzene molecule according to the scheme No. 2 for localizing the scattering centers of 30 electrons. Arrows in Figure 17 indicate double bonds C=C, which differ from single C–C with a slightly higher density of markers in the form of points that indicate the average positions of electrons. Markers in the form of pentagrams and large dots indicate the locations of carbon atoms and protons, respectively. Other parameters of calculation took the following values:

$L=7,N={10}^3,M=500,\epsilon =-\frac{2E_{{\text{C}}_{\text{6}}{\text{H}}_{\text{6}}}}{N}=0.4646,\sigma =0.0185$.

In scheme No.3 we consider the case when the electrons of double carbon bonds are delocalized.

Scheme No. 3. The electrons of hydrogen, four electrons from each of the carbon atoms are uniformly randomly distributed over the benzene ring. When positioning 30 valence electrons, it is considered that in each statistical experiment the electron scattering centers are uniformly randomly positioned over a length of 6 × 2.6465 + 6 × 2.0605 d.u. total benzene ring (Figure 15(a)), where 2.6465 d.u., distance between a pair of carbon atoms, 2.0605 d.u., the distance between carbon and hydrogen.

Figure 18 shows a typical result of the calculation of the average localization sites of quantum particles of benzene in the space (x, y, z) and in the projection on the plane (x, y). The calculation parameters were chosen as follows:

$L=7,N={10}^3,M=500,\epsilon =-\frac{2E_{{\text{C}}_{\text{6}}{\text{H}}_{\text{6}}}}{N}=0.4646,\sigma =0.025$.

Using the example of the scheme No.3, we construct a mixed solution of the Schrodinger equation describing the dynamics of the average trajectories of 30 valence electrons of the benzene molecule.

Figure 16. Scheme No.1. Locations of average positions of quantum particles of a molecule benzene in space (graph on the left) and in projection on a plane (graph on the right).

Figure 17. Scheme No.2. Locations of average positions of quantum particles of a molecule benzene in space (graph on the left) and in projection on a plane (graph on the right).

Figure 18. Scheme No.3. Locations of average positions of quantum particles of benzene molecule in space (graph on the left) and in projection on the plane (graph on the right).

Let’s make a matrix $G=\left|\begin{array}{cc}0& Q\\ -Q& 0\end{array}\right|$ and solve numerically on a time interval [0, 5] a system of linear differential Equation (12) starting from a random initial vector $y_0$ having a unit norm. An analogue of the procedure for constructing the average trajectory of each of the valence electrons is described in section 5 on the example of constructing the average trajectories of the proton and electron of the hydrogen atom. We repeat the indicated procedure for constructing a trajectory M = 50 times, each time changing the positions of the quantum particles of the benzene molecule according to algorithm (53). Other calculation parameters were reduced to the values: $L=7,N=25,\epsilon =0.4646,\sigma =0.025$. The result is shown in Figure 19.

Figure 19(a) shows the appearance of the average trajectories of motion of 30 valence electrons within the entire benzene molecule, while the trajectories of electrons are separated from each other by random colors in various statistical experiments. Due to the high frequency of oscillations, the trajectories merge into a single unit, while the characteristic geometric structure of the benzene molecule is clearly visible. On Figure 19(a) a fragment of the benzene molecule is indicated, which is enlarged and shown in Figure 19(b). On the enlarged fragment of the molecule, separate trajectories are visible. Markers in the form of black dots (small and large) in both Figure 19 indicate the beginning and end of individual trajectories.

The comparison between the patterns of positioning of electrons in Figures 16-19 shows that they are significantly different. According to scheme No.3, the scattering centers of valence electrons of the benzene molecule are uniformly randomly distributed along the benzene ring, including hydride bonds, while the variability corridor in the perpendicular direction is regulated by the parameter σ. In this scheme electrons of double and single bonds of carbon atoms are delocalized.

(a)(b)

Figure 19. (a) Average trajectories of 30 valence electrons of the benzene molecule in 50 calculations; (b) Enlarged fragment of the set of average trajectories of the motion of valence electrons in the benzene molecule.

Let us present in a generalized form the considered numerical algorithm for solving the Schrodinger equation in terms of constructing the pure states of a molecule with the number of quantum particles greater than several tens.

Let the molecule consist of n quantum particles and its energy is $E_{\Sigma }$. We consider that quantum particles of a molecule are characterized by ratios of electron mass to mass of atoms, including electrons, ${\mu }_1,\cdots ,{\mu }_n$ and charges $q_1,\cdots ,q_n$, respectively.

1) We introduce a set of N radius vectors $r_i=\left(r_{1,i},\cdots ,r_{n,i}\right),i=1,\cdots ,N$ in a space of dimension 3n, describing the initial data for determining the positions of quantum particles, which will then change according to point 4) of this algorithm.

2) Create the matrix $Q=\frac{1}{2}\epsilon N{e}_N-\frac{1}{2}\epsilon o_N+diag\left(U_i\right)$, where $\epsilon =-\frac{2E_{\Sigma }}{N}$ is a non-negative parameter ensuring the fulfillment of the virial theorem for the molecule, $e_N$ is the N × N unit matrix, $o_N$ is a special N × N matrix, all elements of which are units, $diag\left(U_i\right)$ is a diagonal matrix, on the diagonal of which are the potential energy values of the molecule at the points $r_i,i=1,\cdots ,N$. The potential energy is calculated by the formula:

$U_i={\displaystyle {\sum }_{j,k=1;j<k}^N\frac{q_j{q}_k}{\left|r_{j,i}-r_{k,i}\right|}},i=1,\cdots ,N$.

3) Find eigenvalues ${\Omega }_1,\cdots ,{\Omega }_N$ and eigenvectors $c_1,\cdots ,c_N$ of the matrix Q. Choose among a set of eigenvalues that of them, ${\Omega }_{\alpha }\cong E_{\Sigma }$, which is closest to the total energy of the molecule. Assuming that the eigenvectors are normalized by one, we find the localization of $R_{k,\alpha },k=1,\cdots ,n$ average positions of quantum particles included in the molecule by the formula: $R_{k,\alpha }={\displaystyle {\sum }_{i=1}^N{r}_{k,i}{c}_{\alpha ,i}^2}$.

4) The procedure in items 1)-3) is repeated M times, assuming that the radius vectors $r_i=\left(r_{1,i},\cdots ,r_{n,i}\right),i=1,\cdots ,N$ are selected each time according to the scheme of the form: $r_{k,i}=a_k+\sigma 2L\sqrt{{\mu }_k}{\eta }_{k,i},k=1,\cdots ,n;i=1,\cdots ,N$, where $a_k,k=1,\cdots ,n$ the so-called scattering centers of quantum particles, which in general can vary from one statistical test to another; σ is some non-negative fitting coefficient; L is the characteristic size of the task or, otherwise, the size of the three-dimensional box, ${\left[-L,L\right]}^3$, in which the molecule is placed; ${\eta }_{k,i},k=1,\cdots ,n;i=1,\cdots ,N$ are vectors of random variables whose coordinates have zero expectation and dispersion of the order of unity.

5) We construct the final graphs of the average positions $R_{k,\alpha },k=1,\cdots ,n;\alpha =1,\cdots ,M$ of quantum particles of the molecule in three dimensional physical space.

10. Metallic Hydrogen

As the final stage of testing the proposed scheme of numerical solution of the Schrodinger equation, we consider the so-called “metallic hydrogen”. It is believed [15] that metallic hydrogen can be obtained under high pressure. In this case, electrons are separated from specific protons and become free, and hydrogen acquires the properties of metal, in particular, becomes a conductor of electric current. In astrophysics [16], metallic hydrogen is mentioned in connection with the giant planets (Jupiter, etc.), where, due to high gravity, conditions are formed for the formation of metallic hydrogen.

Consider a molecule consisting of n hydrogen atoms. Let us estimate the total energy $E_{mh}$ of the molecule, when it is considered that there is a phase called metallic hydrogen. The energy of the molecule can be collected by summing the internal energy of n hydrogen atoms equal to –0.5n (d.u.), plus the binding energy of the metal hydrogen molecule, $E_{bond}$, which is not known to us.

For certainty, we assume that the proton scattering centers are concentrated in the nodes of a simple cubic lattice fragment, with a characteristic distance between the nearest nodes b. In this case, it is considered that the number of quantum particles can be represented as a product of two on the cube of some natural number. At each step of the statistical experiment, each electron scattering center is uniformly randomly connected to one of the proton scattering centers, then the electrons become delocalized and indistinguishable.

Thus, at a pressure of 10⁸ bar, the density of metallic hydrogen [16] is 5.90 g/cm³, where you can find the average distance between protons equal to b = 0.6569 Å = 1.2413 d.u. Recall that the distance between a pair of protons in a hydrogen molecule is considered to be equal to b' = 0.7416 Å = 1.4014 d.u.

To estimate the binding energy $E_{bond}$ we assume that it is equal to the change in the potential energy $\Delta U_p=U_p\left(b\right)-U_p\left(b^{\prime }\right)$ when n protons are closer together, where $U_p\left(b\right)$ is the Coulomb interaction potential of protons placed in the nodes of the cubic lattice, whose step takes the value b. As a result, we will write down $E_{mh}=-0.5n+E_{bond}=-0.5n+\Delta U_p$. The energy $E_{bond}=\Delta U_p$ defines the assumed part of the total energy, which is characteristic for the phase of metallic hydrogen. Considering this energy as negative, we believe that hydrogen atoms are more closely related to each other in metallic hydrogen.

We apply the calculation algorithm formulated at the end of the previous section M = 500 times and construct the final graphs of the average positions of the quantum particles of metallic hydrogen. A typical example of the final calculations is shown in Figure 20. Other parameter values were selected as follows:

$\begin{array}{l}L=4,N={10}^3,n=5^3=125,E_{mh}=-322.8112,\\ E_{bond}=-260.3112,\epsilon =-\frac{2E_{mh}}{N}=0.6456,\sigma =\mathrm{0.02.}\end{array}$

As in the previous examples, the parameter ε was chosen on the basis of the conditions of the virial theorem.

Figure 20 shows the characteristic results of the calculation of a small-sized metallic hydrogen molecule. It should be noted that due to the enormous value of the modulus of the binding energy $E_{bond}$, the electrons are noticeably locked in the vicinity of protons, and the value of the parameter σ is chosen to be rather small. The reality status of metallic hydrogen constructed in this example is

Figure 20. Kinetic, potential and total energy (graph on the left); locations of the average positions of protons (markers in the form of large points) and electrons (markers in the form of small points) of the metal hydrogen molecule on the right graph.

questionable, because the specific binding energy per hydrogen atom is $\frac{E_{mh}}{n}=-\frac{322.8112}{125}=-2.5825$, which is considerably less than the minimum of the binding energy of the hydrogen atom is –0.5.

11. Conclusion

The paper presents a numerical method for solving the Schrodinger equation for a molecule with an arbitrary set of quantum particles. This method combines the finite-difference and Monte-Carlo approaches, which allowed to remove the problem of multidimensionality of the Schrodinger equation. The resulting method was effective and economical and, to a certain extent, not improved, i.e. optimal. The method is presented and simultaneously outlined on the example of solving a number of problems, namely, a linear one-dimensional oscillator, a hydrogen atom, an ion and a molecule of hydrogen, water, benzene, and metallic hydrogen. The method itself is formalized as an algorithm for the numerical solution of the Schrodinger equation for a molecule with an arbitrary number of quantum particles. The method can be a useful addition to other methods of calculation of molecular systems in physics and chemistry.

Acknowledgements

The author is grateful for help in translating article V.N. Nikolenko.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Ozhigov, Y.I. (2010) Constructive Physics. Research Center “Regular and chaotic dynamics”, Moscow, Izhevsk, 424 p.
[2]	Stepanov, N.F. (2001) Quantum mechanics and quantum chemistry. Mir, Moscow, 519 p.
[3]	Motta, M., Ceperley, D.M., Chan, G.K.-L., Gomez, J.A., Gull, E., Guo, S., et al. (2017) Towards the Solution of the Many-Electron Problem in Real Materials: Equation of State of the Hydrogen Chain with State-of-the-Art Many-Body Methods. Physical Review X, 7, Article ID: 031059.
[4]	Kim, J., Baczewski, A.T., Beaudet, T.D., Benali, A., Chandler Bennett, M., Berrill, M.A., et al. (2018) QMCPACK: An Open Source ab Initio Quantum Monte Carlo Package for the Electronic Structure of Atoms, Molecules and Solids. Journal of Physics Condensed Matter, 30, Article ID: 195901. https://doi.org/10.1088/1361-648X/aab9c3
[5]	Hartree, D. (1960) Calculations of Atomic Structures. Foreign Literature Publishing House, Moscow, 271 p.
[6]	Kohn, W. (1999) Nobel Lecture: Electronic Structure of Matter-Wave Functions and Density Functionals. Reviews of Modern Physics, 71, 1253-1266. https://doi.org/10.1103/RevModPhys.71.1253
[7]	Plokhotnikov, K.E. (2020) About One Method of Numerical Solution of Schrodinger’s Equation. Mathematical Models and Computer Simulations, 12, 221-231. https://doi.org/10.1134/S2070048220020106
[8]	Landau, L.D. and Livshits, E.M. (1974) Quantum Mechanics. Non-Relativistic Theory. Vol. 3, Pergamon Press, Moscow, 752 p.
[9]	Nekrasov, B.V. (1973) Fundamentals of General Chemistry. Vol. 1. Moscow, Khimia, 656 p.
[10]	Nikolsky, B.P. (1966) Chemist Handbook. Vol. 1, Khimia, Moscow, Leningrad, 1072 p.
[11]	Ravdel, A.A. and Ponomareva, A.M. (2003) Quick Reference Physico-Chemical Variables/SPb. Ivan Fedorov, Moscow, 240 p.
[12]	Gurvich, L.V., Karachevtsev, G.V., Kondratiev, V.N., et al. (1974) The Energies of Breaking the Chemical Bond. Ionization Potentials and Electron Affinity. Science, Moscow, 351 p.
[13]	Akhmetov, N.S. (2001) General and Inorganic Chemistry. Higher. School, Ed. center “Academy”, Moscow, 743 p.
[14]	Poling, L. (1960) The Nature of the Chemical Bond. Cornell University Press, Ithaca, 655 p.
[15]	Dias, R.P. and Silvera, I.F. (2017) Observation of the Wigner-Huntington Transition to Metallic Hydrogen. Science, 355, 715-718.
[16]	Zharkov, V.N. (1983) The Internal Structure of the Earth and Planets. Gordon and Breach, New York, 416 p.