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

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.


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.
( ) 1 , , , n t = r r  ψ ψ , where t is time, 1 , , n r r  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.

One Quantum Particle in One-Dimensional Space
We write the Schrodinger equation for a single quantum particle in the one dimensional case, i.e.
where  is the Planck's constant, As you know, the Schrodinger Equation (1) admits the existence of a probability conservation law, which can be written as an equation: 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: 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 where , i j g i j 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 , , where the point is the time derivative, and ( ) 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 , , , , 1, , 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.
From Equation (8) it follows that as the matrix g should choose a symmetric matrix with non-negative elements of the form: 2 , , 1; 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: where 1 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 , , Taking into account the vector y and the matrix G, the system of Equation (11) can be rewritten in the following compact form: Compare the solutions of Equations (8) and (12). To do this, substitute in (8) the solution in the form  = ± ± π  and after substitution in (8), we obtain an expression for the energy: 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 e i t y z Ω = , where z is some constant vector of dimension column 2 1 N × . In this case, the value of iΩ 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: where N e , unit matrix of size N N × . We write the solutions of Equation (14) as a set 1 , , N ±Ω ±Ω  . Choose the top plus sign. In this case, we obtain a set of frequencies 1 , , N Ω Ω  , which can be compared with the set (13), when 0, , To ensure coherence of the pair of energy spectra 0 , 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: app h and app L . From the above values, one can also find the number of nodes app N 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 app h and app L will be found by comparing the spectra, i.e. by minimizing the total comparison error of the type: Let function (15) 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 1 2 , ρ ρ . In the calculations, the results of which are shown in Figure 1 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 Journal of Applied Mathematics and Physics ( ) ( ) Since the systems of Equations (11), (12) preserve the probability, so far the equation will take place at all subsequent moments of the time ( 0 t > ). The last equation means that we have a discrete random variable X, the position of the quantum particle, which takes the values 1 , , N x x  with probabilities You can also determine the average position of a quantum particle X and its dispersion X D according to classical formulas: Since the probabilities in the set depend on time, the average position ( ) of a quantum particle generally depend on time.  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 of a random i'-th subunit of a quantum particle placed with- . The right graph in Figure 2 shows the probability distribution of finding a quantum particle in space at some time point t. Other design parameters were chosen as follows: . 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 [ ] , L L − . 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: According to the left graph in Figure 3 the average position of the quantum particle X vibrates, and the dispersion of vibrations X D also experiences vibrations.

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 ( ) In space we will define a cube [ ] 3 , L L − 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 ( ) for the other two variables: In total, thus, in our grid that fills the selected cube, there are 3 N 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: We generalize the expression for the matrix g in (10) to the three-dimensional case. To do this, we introduce multicomponent numbers ( ) , , Note that the finite-difference representation of (18) can be interpreted in the sense that the Laplace operator 2 2 2 x y z , , , , j j j = j .
It is possible to renumber the subparticles not by three subindexes, but by one. Let be ( ) We introduce a single (one-component) numbering i of all grid nodes in a three-dimensional cube using the formula: 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  To ensure the proximity of a pair of spectra, the problem of minimizing the functional S was solved: where before of the numbering and sorting of frequencies (20) it was believed that ( ) found. Figure 4 shows several samples comparing a pair of frequency spectra for the three cases when 7,10,15 N = , the other parameters took the following values: . 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 app h and app L were calculated by minimizing the functional (23).

Linear Oscillator in One-Dimensional Case
Let us write the potential function of a one-dimensional oscillator in the form ( ) We choose a certain interval [ ] , L L − of integration over the space of the system of Equation (24). We define at this interval generally speaking uneven grid with nodes: . In accordance with the studied solution scheme, we rewrite Equation (24) as a set of the following differential equations: is taken from (10).
If you enter a column vector ( )  Find solutions to the linear system of Equation (12). To do this, we will look for a solution in the form: Substitute the selected solution into Equation (12), then we obtain the following problem for finding the eigenvalues (δ) and vectors (z) of the matrix G: The problem (26) is reduced to a simpler problem of finding eigenvalues (Ω) and vectors ( ( ) According to (27) where 1 2 , C C 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 1 , 0,1, 2 n E n n = + =  we write out the known first two Journal of Applied Mathematics and Physics exact wave functions: Taking into account (29), we find the probability density =ψ ψ and compare it with the probability density distributions of solutions (28) for the first two eigenvalues Ω 1 and Ω 2 , and that eigenvalues are ordered in ascending order. We denote the corresponding probability distribution, constructed using the solutions to (28), in the form: 1, , 1, , for the eigenvalues Ω 1 and Ω 2 respectively. Figure 5 shows a comparison of the approximate (solid line) and the exact solution (dashed line). In Figure 5(a), the distribution densities 1, , 1, , x were compared, and in Figure 5(b), 2, , 1, , In the calculations of approximate solutions, the following parameter values were used: 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 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.

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 1 m and 2 m be the masses of proton and electron. We write the Schrodinger equation describing the dynamics of the hydrogen atom as: where ( ) , wave function of the hydrogen atom, 1 2 , r r , 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: We rewrite Equation (30) in dimensionless form. To do this, we introduce the characteristic length, mass, time and energy: 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: ter these substitutions, the system of Equation (32) will be rewritten as a system of linear ordinary differential equations of the form: , , 1 , 2 To complete the algorithm for calculating the hydrogen atom according to the scheme (33), it is necessary to construct a matrix , , , 1, , Considering the examples of the construction of this matrix in (9), (10), (21), we choose the following type of matrix: where , 1, , i j N =  ; 1 2 , ρ ρ , some non-negative constants, h, the step of the equivalent uniform grid, which we define below.
Taking into account the definition of ( ) 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: , a a , some constant vectors, which we will interpret later as the scattering centers of possible positions of proton and electron. We assume Journal of Applied Mathematics and Physics The representation (36), in which independent uniformly random variables 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 1 , , ± Ω  . Taking into account the discrete nature of the energy spectrum of a hydrogen atom, which, as is known, is characterized by the dependence . 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 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 For a specific example, for each of the N eigenvalues we calculate the number of points , 1, , k n k N =  , in which the hydrogen atom is noticeably localized. As a criterion of noticeable localization, we will choose a certain probability threshold th p . In this case, you can write the following formula to calculate the distribution of the number of localization points: w , the probability of finding a hydrogen atom at the i-th point of space for the k-th eigenvalue, ( ) ... θ , 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:  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 th p , 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 Define kinetic spectra, E ke , potential, E pe , and total, E te , energy by the formulas: 1 , 2 2 sort ,  . 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 1 R , 2 R proton and electron were calculated accordingly by analogy with the formula (16), i.e.
( ) ( ) 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.

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: 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 ( ) 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 Let's choose among the set of eigenvalues the value, Ω α , 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 Ω α , we find the eigenvector c α of the matrix Q, which is normalized by one. Carrying out this procedure M times, we find the set 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  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). Journal of Applied Mathematics and Physics 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) 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 1,

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 1 2 , , a a  , 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 1 4 , a a and two electrons 2 3 , a a , then we can write:

d =
, 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 2 0.5 0.1647 We introduce N points in twelve-dimensional space with radius vectors ( )   Define the matrix  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 1 0 = ρ , or with 2 → ∞ ρ . In the second stage, we choose the values of 1 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 Ω α , 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 Ω α , we find the eigenvector c α of the matrix Q, which is normalized by one. Conducting this procedure M times we find the set 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, pe E . Journal of Applied Mathematics and Physics 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 ( ) 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 1 2 , ρ ρ must be chosen appropriately.
To select the appropriate values of the coefficients 1 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 1.1647 E = − 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: 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.

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 . Finally, 216.87 kcal/mol, is required to break both covalent bonds in a water molecule [12], which in dimensionless form will be H O 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.  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 1 13 , ,  µ µ ratios of the electron mass to each of the 13 masses of quantum particles of the water molecule. In this case, we get: We also introduce a set of charges 1 13 , , q q  of each of the 13 particles of a water molecule in dimensionless form, then 0,0,0 ; cos , sin ,0 ; , cos , sin ,0 ; , ( )  Figure 13 shows the results of the computational experiment. 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

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 6 H 6 , 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 .
where subindexes x, y, z denote projections of vectors Taking into account (52) the eigenvalue problem for the matrix ( ) can be solved analytically. In this case, the kinetic energy spectrum will have two non-negative values, 1 2 , Ω Ω , one of which is zero. The corresponding characteristic equation has the form:  Let us calculate the energy of complete dissociation of the benzene ring, 6 6 C H E . So, the ionization energies of a carbon atom are (MJ/mol) [12] As a result, we find the total dissociation energy of the benzene ring: Taking into account (49), we write the algorithm for generating the set of vectors , 1, , where 1, , 54; 1, , 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 Journal of Applied Mathematics and Physics . Where it follows that 6 6 C H 2E 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. 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: 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:   . 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 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 Σ . We consider that quantum particles of a molecule are characterized by ratios of electron mass to mass of atoms, including electrons, 1 , , n  µ µ and charges 1 , , n q q  , respectively. 1) We introduce a set of N radius vectors

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 hydro-gen 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 mh E 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, bond E , 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 8 bar, the density of metallic hydrogen [16] is 5.90 g/cm 3 , 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 bond E we assume that it is equal to the change in the potential energy 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:

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.