American Journal of Applied Mathematics
Volume 3, Issue 3, June 2015, Pages: 106-111

Estimation of Boron Ground State Energy by Monte Carlo Simulation

K. M. Ariful Kabir1, Amal Halder2

1Department Mathematics, Bangladesh University of Engineering and Technology, Dhaka, Bangladesh

2Department of Mathematics, University of Dhaka, Dhaka, Bangladesh

Email address:

(K. M. A. Kabir)

To cite this article:

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. 106-111. 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 [1-3]. 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 [4-5]. 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 two-electron system, which considered as the simplest few-body 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 water-solubility 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 inter-electronic 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 five-electron atomic ions, most of these studies search for accurate results but nevertheless compact wave functions. From this point, [9-13] 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 few-parameter 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. Many-Bodies 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


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


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


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 many-body 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 many-body wave functions enable the accurate determination of expectation values. Variational Monte Carlo is a method of computing the total energy


and its variance ( statistical error)


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 two-electron 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 wave-function 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]:


where, the symbol denotes the angular variables of and the function  is the hydrogen-like wave function in the nl-state 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 electron-electro correlation we have to modify the form of this trial wave function by multiplying by the following factor:


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 ground-state of the BoronEq. (6) reduces to the following form:


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 ground-state of Beryllium.

5. Estimate of the Smallest Eigenvalue

To estimate the smallest eigenvalue, Let us rewrite equation (2) more compactly as


where  represents the sum of both operators on the left-hand side of equation (9). The variational principle tells us [17] that


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 ground-state energy, the left-hand side of equation (10) is usually quite close to. The main problem is how to evaluate the two 12-dimensional integrals; this is impossible to do analytically and not feasible even numerically. But, Monte Carlo is the process that can easily simulate 12-dimensional integrals.

Let us first define the local Energy by,


The left hand side of (10) can now be written as


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


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 so-called 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


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.

Figure 1. 1000 iterations of VMC to reach equilibrium.

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


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 grand-mean 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 (grand-mean) 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.

Table 1. Monte Carlo Result.

  Estimate Standard Error P-Value


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

Table 2. Energy Values.

Experiment Energy
Hylleraas Experiment

Monte Carlo Method

Table-1 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 Table-2 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 electron-electron 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 Table-1 that our result for the Boron atom (Z=5) is, which is very accurate.

Figure 2. Presenting the energy level as a function of Boron and using a set of Monte Carlo points.

Table 3. Comparison of the energies for boron atom calculated with different methods and the nonrelativistic energy (all values in a.u.)[8].

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
Froese-Fisheret La 1993 Numerical HF -24.529036
Mayer 2004 Full-CI -24.530874
Ruiz 2004 Hylleraas -24.541246
Froese-Fisher La 1994 Multiref. CI -24.560354
Estimated nonrelativistic 1993 Full-CI -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 Froese-Fisher et la result and is comparable in accuracy to Ruiz result (Table-3).

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, Froese-Fisheret La, Mayer and Huzinaga (Table-3), we see that our result is an improved one when compared with experimental result. In a follow-up 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 "self-sufficient. This means that VMC method introduced a very well description for the Beryllium ground state using a trial wave function expressed in hydrogen-like orbital’s with a rather simple factor, describing the correlations between the electrons.


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.


  1. B. L. Hammond, W. A. Lester Jr and P. J. Reynolds, Monte Carlo methods in ab initio quantum chemistry, pp.1-10,(1994).
  2. E. Buendía, F. J. Gálvez, A. Sarsa, Chem. Phys. Lett., Explicitly correlated energies for neutral atoms and cations ,465: pp.190-192 (2008).
  3. S. A. Alexander, R. L. Coldwell, Int. J. Quantum Chem, Atomic wave function forms, 63,pp.1001-1022 (1997).
  4. K. E. Riley, J. B. Anderson, Chem. Phys. Lett., A new variational Monte Carlo trial wave function with directional Jastrow functions 366: pp.153-156 (2002).
  5. S. A. Alexander, R. L. Coldwell, Int. J. Quantum Chem, Ro-vibrationally averaged properties of H2 using Monte Carlo methods, 107: pp.345-352 (2007).
  6. S. Datta, S. A. Alexander, R. L. Coldwell, Int. J. Quantum Chem. The lowest order relativistic corrections for the hydrogen molecule, pp.731-739 (2012).
  7. Buendı´a, F.J. Ga´lvez, A. Sarsa, Chem. Phys. Lett., Correlated wave functions for the ground state of the atoms Li through Kr, 428: pp.241-244 (2006).
  8. Ruiz, M.B and Rojas,M. Computational Method in Science and Technology, Variational Calculations on the 2P1/2 Ground State of Boron Atom using Hydrogenlike orbitals.9(1-2), 101-112 (2003)
  9. A. M. Frolov, Eur. Phys. J. D, Atomic Excitations During the Nuclear β− Decay in Light Atoms, 61: pp. 571-577 (2011).
  10. N. L. Guevara, F. E. Harris, A. Turbiner, Int. J. Quantum Chem., An Accurate Few-parameter Ground State Wave Function for the Lithium Atom, 109: pp.3036-3040 (2009).
  11. S. B. DomaandF. El-Gamal, Monte Carlo Variational Method and the Ground-State of Helium, pp-78-83.
  12. C. Schwarz, Many-body methods in quantum chemistry: proceedings of the symposium (1989).
  13. W. Kutzelnigg and J. D. Morgan III,Explicitly Correlated Wave Functions in Chemistry andPhysics. Phys. pp. 96-104,(1992).
  14. D. Ruenn Su, Theory of a Wave Function within a Band and a New Process for Calculation of the Band Energy, pp. 17-28 (1976)
  15. P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, Jr., The Journal of Chemical Physics, 77(12), 1982 pp. 5593–5603. doi:10.1063/1.443766.
  16. J. B. Anderson, The Journal of Chemical Physics, 73(8), 1980 pp. 3897–3899. doi:10.1063/1.440575
  17. D. M. Ceperley and B. J. Alder, "Quantum Monte Carlo," Science, 231(4738), 1986 pp. 555–560. doi:10.1126/science.231.4738.555.

Article Tools
Follow on us
Science Publishing Group
NEW YORK, NY 10018
Tel: (001)347-688-8931