Estimation of Boron Ground State Energy by Monte Carlo Simulation
K. M. Ariful Kabir^{1}, Amal Halder^{2}
^{1}Department Mathematics, Bangladesh University of Engineering and Technology, Dhaka, Bangladesh
^{2}Department of Mathematics, University of Dhaka, Dhaka, Bangladesh
Email address:
K. M. Ariful Kabir, Amal Halder. Estimation of Boron Ground State Energy by Monte Carlo Simulation. American Journal of Applied Mathematics. Vol. 3, No. 3, 2015, pp. 106111. doi: 10.11648/j.ajam.20150303.15
Abstract: Quantum Monte Carlo (QMC) method is a powerful computational tool for finding accurate approximation solutions of the quantum many body stationary Schrödinger equations for atoms, molecules, solids and a variety of model systems. Using Variational Monte Carlo method we have calculated the ground state energy of the Boron atom. Our calculations are based on using a modified five parameters trial wave function which leads to good result comparing with fewer parameters trial wave functions presented before. Based on random Numbers we can generate a large sample of electron locations to estimate the ground state energy of Boron. Based on comparisons, the energy obtained in our simulation are in excellent agreement with experimental and other well established values.
Keywords: Monte Carlo Simulation, Boron, Ground State Energy, Schrödinger Equation
1. Introduction
Variational Monte Carlo (VMC) method has become a powerful tool in Quantum Chemistry calculations [13]. In most of its current applications the VMC method has become a valuable method because of a wide variety of wave function forms for which analytical integration would be impossible. The major advantage of this method is the possibility to freely choose the analytical form of the trial wave function which may contain highly sophisticated term in such a way that electron correlation is explicitly taken into account [45]. This is an important feature valid for quantum Monte Carlo methods, which are therefore extremely useful to study physical cases where the electron correlation plays a crucial role. For twoelectron system, which considered as the simplest fewbody systems, VMC method provides accurate estimations of the ground and excited state energies and properties of atomic and molecular systems [6]. Boron is a chemical element with symbol B and atomic number 5. The word boron was coined from borax, the mineral from which it was isolated, by analogy with carbon, which it resembles chemically. Boron is concentrated on Earth by the watersolubility of its more common naturally occurring compounds, the borate minerals. These are mined industrially as evaporates, such as borax and kernite. The largest proven boron deposits are in Turkey, which is also the largest producer of boron minerals.
The Hylleraas method is a variational method which introduces the correlation effects by including explicity the interelectronic distances in the wave function. In recent years efforts have been done to develop an approach for constructing trial wave functions in order to calculate the ground state energy of the Boron atom and achieve high level of accuracy. In 2006, the VMC method has been used to study the ground state energy of the atoms Li through Kr using explicitly correlated wave functions which consists of product of a Jastrow correlation factor times a model function with 17 variational parameters [7]. The obtained results were an improvement over all the previous results. Recently, calculations on the ground state of boron atom have been made using the single and 150 term wave functions constructed with Slater orbitals by M.B Ruiz. [8].
Despite the fact that there is no shortage of extremely accurate wave functions for the Boron atom and some fiveelectron atomic ions, most of these studies search for accurate results but nevertheless compact wave functions. From this point, [913] proposed a simple compact trial function for the ground state of the Boron atom which provided a very accurate energy in such a way that it could considered as the most accurate among existent fewparameter trial functions.
In the present paper we shall use VMC method to study five electron system (or Boron atom), which may be studied by using different method and different form of trial wave functions. Accordingly, we shall use a compact five parameters wave function to obtain the ground state energy for Boron ground state. The calculation will be done in the frame of VMC method.
2. ManyBodies Stationary Schrödinger Equation
Let us consider a system of quantum particles, such as electrons and ions interacting via Coulomb potentials. Since the masses of nuclei and electrons differ by three orders of magnitude or more and the Hamiltonian is given by
(1)
where i and j label the electrons and I runs over the ions with charges . Throughout the review, we employ the atomic units,, where is the electron mass, is the electron charge and is the permittivity of a vacuum. The Schrödinger equation is
(2)
Here, is the Laplacian operator, is a function of their positions, and is the molecule’s ground state energy. The electrostatic potential in the molecule is given by
(3)
Here, the first summation is over all electrons and nuclei, where and are the nuclear charges and locations, respectively. The second summation is over all pairs of electrons.
3. Variational Monte Carlo for the Boron Atom
Variational Monte Carlo (VMC) is based on a direct application of Monte Carlo integration to explicitly correlated manybody wave functions. The variational principle of quantum mechanics, states that the energy of a trial wave function will be greater than or equal to the energy of the exact wave function. Optimized forms for manybody wave functions enable the accurate determination of expectation values. Variational Monte Carlo is a method of computing the total energy
(4)
and its variance ( statistical error)
(5)
Here, H is the Hamiltonian, is the value of the trial wave function at the Monte Carlo integration point. The constant is fixed at a value close to the desired state in order to start the optimization in the proper region. The exact wave function is known to give both the lowest value of and a zero variance. If the adjustable parameters in the trial wave function are optimized so as to minimize the energy, instability often occurs. This happens when a set of parameters causes to be estimated a few sigmas too low. Although such parameters will produce a large variance, they are favored by the minimization. This problem can be avoided only by using a very large number of configurations during the optimization of the wave function so as to distinguish between those wave functions for which is truly low and those which are merely estimated to be low. In contrast, variance minimization favors those wave functions which have a constant local energy. Parameter values which do not produce this property will be eliminated by the optimization process. As a result, only a small fixed set of configurations is needed to accurately determine the variance. Previous studies have shown that the rate of convergence of a variational calculation can be tremendously accelerated by using basis functions which satisfy the twoelectron cusp condition and which have the correct asymptotic behavior [12]. Unfortunately, the integrals of such functions can rarely be evaluated analytically. Because our method uses Monte Carlo integration, we can easily build into the trial wave function many features which will accelerate convergence. Although, in principle, this flexibility leads to an enormous number of possible forms, in practice, the ideal trial wavefunction form must have a low variance, must add adjustable parameters in a straightforward manner, and must be easy to optimize.
4. The Trial Wave Function
The calculations for the Boron atom obtained previously using trial wave functions which takes the form [14]:
 (6) 

where, the symbol denotes the angular variables of and the function is the hydrogenlike wave function in the nlstate with effective charge .
This wave function was used to calculate the energy of Beryllium ground state with quite accurate results. In the present paper we introduce some modifications to this trial wave function in order to get more accurate results. Firstly, we will consider the correlation between each two electrons. In order to include the electronelectro correlation we have to modify the form of this trial wave function by multiplying by the following factor:
(7)
which expresses the correlation between the two electrons i and j due to their Coulomb repulsion. That is, we expect f to be small when is small and to approach a large constant value as the electrons become well separated. Then, for the groundstate of the BoronEq. (6) reduces to the following form:
 (8) 

where the indices i, j run over the number of the electrons. The variational parameters and will determined for each value of Z by minimizing the energy. The function of equation (6) is constant for the groundstate of Beryllium.
5. Estimate of the Smallest Eigenvalue
To estimate the smallest eigenvalue, Let us rewrite equation (2) more compactly as
(9)
where represents the sum of both operators on the lefthand side of equation (9). The variational principle tells us [17] that
(10)
The integration is over the three coordinates of each of the four electrons, altogether a 12 dimensional problem and is any trial solution to equation (2). The limit of equality holds only for the exact solution, but for approximate solutions, called variational estimates of the groundstate energy, the lefthand side of equation (10) is usually quite close to. The main problem is how to evaluate the two 12dimensional integrals; this is impossible to do analytically and not feasible even numerically. But, Monte Carlo is the process that can easily simulate 12dimensional integrals.
Let us first define the local Energy by,
(11)
The left hand side of (10) can now be written as
(12)
Next, we randomly generate a large sample of 10000 configurations of 12 variables denoted collectively as and compute the corresponding, D and. By averaging the 10000 values of , we get an estimate of . Unfortunately, this estimate will be very inaccurate since our random sample of configurations bears, at this point, no relationship to of equation (12).
To fix this, we move each configuration to a new location, specified by
(13)
where is the drift function evaluated at the old location , N is a random vector of 12 independent components from the normal distribution (with mean 0 and standard deviation 1), and is an extra parameter called the step size, which controls the speed of this motion. This will bring us a step closer to the desired distribution of configurations whose probability density function is proportional to, but it will take dozens of such moves to reach it. Monitoring the consecutive sample averages of, we find no systematic change but only random fluctuations after reaching a socalled equilibration. Once in equilibrium, we continue advancing our configuration for as many steps (called iterations) as feasible, to reduce the statistical error of the final estimate. This is computed by combining all the individual sample averages into one. There is only one little snag: the result will still have an error proportional to the step size. To correct for this, we would have to make impractically small and equilibration would take forever. Fortunately, there is another way, called Metropolis sampling [19], for each proposed move (13) we compute a scalar quantity
(14)
where the subscripts n and o mean that and D have been evaluated at the new or old location, respectively. The move is then accepted with a probability equal to T. When T > 1, the move is accepted automatically. When a move is rejected, the configuration simply remains at its old location. The step size should be adjusted to yield a reasonable proportion of rejections, say between 10% and 30%. Rejecting configurations in this manner creates the last small problem, in our original random sample there is usually a handful of configurations which, because they have landed at "wrong" locations, just would not move. To fix this, we have to monitor, for each configuration, the number of consecutive times a move has been rejected, and let it move, regardless of T, when this number exceeds a certain value. After this is done and the sample equilibrates, the problem automatically disappears, and no configuration is ever refused its move more than six consecutive times.
To get an accurate estimate of, we repeat the simulation with substantially more iterations, and then computing the grand mean of the values. In our case, this yields atomic units, with an average acceptance rate of about 90%. The easiest way to find the corresponding statistical error is to execute the same program, independently, 5 to 10 times, and then to combine the individual results. This improves the estimate to atomic units, with the standard error of ±0.05326.
The obvious discrepancy, well beyond the statistical error, between our estimate and this value is due to our use of a rather primitive trial function. In accordance with the variational principle, our estimate remains higher than the exact value.
6. Monte Carlo Estimation
Replacing by in equation (12), the expression then yields "nearly" the exact value of, subject only to a small nodal error [16]. So, all we need to do is to modify our simulation program accordingly, to get a sample from a distribution whose probability density function is proportional to instead of . This can be achieved by assuming that each configuration carries a different weight, computed from
(15)
where is the local energy of the configuration as computed iterations ago, the summation is over all past iterations, and is a rough estimate of (the variational result will do). The sum in (10) "depreciates" the past values at a rate that should resemble the decrease in serial correlation of the sequence, which can be easily monitored during the variational simulation.
The new estimates of are then the correspondingly weighted averages, computed at each step and then combined in the usual grandmean fashion. There are two slight problems with this algorithm, but both can be easily alleviated.
Firstly, when an electron moves too close to a nucleus, may acquire an unusually low value, making the corresponding rather large, sometimes larger than all the remaining weights combined. We must eliminate "outliers" outside the range. It is better to do this in a symmetrical way by truncating the value to the nearest boundary of the interval.
Secondly, The final (grandmean) estimate may have a small, proportional bias. This can be removed only by repeating the simulation, preferably more than once, at several (say 3 to 5) distinct values of , and getting an unbiased estimate of by performing a simple polynomial regression. It is the intercept of the resulting regression line (corresponding to ) that yields the final answer.
7. Result and Discussion
In this paper we have calculated the ground state energy of the Boron atom using a modified wave function. Moreover, we have succeeded to use this trial wave function to obtain the energy of the Boron ground state by using VMC Method. The iteration averages of will show that equilibration now takes many more steps (about 1000, when ) than in the case of variational simulation. We have thus decided to discard the first 1000 results and partition the remaining 6000 into six blocks of 1000. We can produce six such values with and . It is now easy to find the resulting intercept.
Estimate  Standard Error  PValue  
1 



x 







This yields the value of for the corresponding intercept. This is in Reasonable agreement, in view of the nodal error, with the experimental value of 24.541246 atomic units. This visualizes the regression fit is
Experiment  Energy 
Hylleraas Experiment 

Monte Carlo Method 

Table1 shows the final results of the energy using a very large of MC points together with the standard deviation given by equation (4) and equation (5). It is clear from Table2 that, the obtained result is of good agreement with the recent experimental Hylleraas value. Also, the associated standard deviations have very small value; this is due to the large number of MC points. The electronelectron correlation factor, which has been included in the trial wave function, played a crucial role in improving the results. In fact, the ground state energy of the Boron atom which was obtained previously using the same trial wave function, but without introducing the factor, was . It is clear from Table1 that our result for the Boron atom (Z=5) is, which is very accurate.
Reference  Year  Method  Energy 
Clementi  1989  HF  24.495670 
Clementi, Roetti  1974  HF  24.498369 
Ruiz  2004  Ref Hy  24.498369 
Clementi, Raimondi  1963  HF  24.498370 
Huzinaga  1977  HF  24.528709 
FroeseFisheret La  1993  Numerical HF  24.529036 
Mayer  2004  FullCI  24.530874 
Ruiz  2004  Hylleraas  24.541246 
FroeseFisher La  1994  Multiref. CI  24.560354 
Estimated nonrelativistic  1993  FullCI  24.65391 
Ruiz[8]  2004  Hydrogen Like Orbital  24.501187 
Present work  Monte Carlo Method  24.546707 
Finally, in Table 3, we compare our energy results with the ones of other methods. The ground state energy obtained with Variational Monte Carlo Method was 24.546707 a.u. This energy lies 0.056% above the FroeseFisher et la result and is comparable in accuracy to Ruiz result (Table3)^{}.
Only the energy obtained by Ruiz [8] lies closer to our result and the error was only 0.02%. By comparing with the values of Clementi, Raimondi, FroeseFisheret La, Mayer and Huzinaga (Table3), we see that our result is an improved one when compared with experimental result. In a followup article, it showed that how Monte Carlo procedure can be extended to estimate atomic properties, including geometry and polarizability, etc. and how to optimize Parameters of a trial function, to make the Monte Carlo method more "selfsufficient. This means that VMC method introduced a very well description for the Beryllium ground state using a trial wave function expressed in hydrogenlike orbital’s with a rather simple factor, describing the correlations between the electrons.
Acknowledgments
I would like to express my very great appreciation to my teacher and colleague for their willingness to support me. I also like to acknowledge the financial, academic and technical support of Bangladesh University of Engineering and Technology, and its staff, that provided the necessary financial and logistic support.
References