We show how a Monte Carlo procedure (based on random numbers) can generate a large sample of electron locations in any simple molecule. Based on this sampling, we can accurately estimate the molecule’s ground-state energy and other properties of interest. We demonstrate this using the LiH molecule.

Trial Solution to Schrödinger Equation

A mathematical description of a molecule (say LiH) is provided by a solution to the Schrödinger equation [1] (a differential eigenvalue problem with smallest eigenvalue ):

(1)

Here is the Laplacian operator , the summation is over all electrons, is a function of their positions, and is the molecule’s ground-state energy. The electrostatic potential in the molecule is given by

(2)

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, and the last one, analogously, is over all pairs of nuclei. The nuclear locations are kept fixed in accordance with the Born-Oppenheimer approximation [2]. We use atomic units, in which the electron’s charge, mass, and Planck’s constant are set equal to 1.

Thus, for the LiH molecule is computed by the following functions, assuming that the Li nucleus is located at the origin.

The argument Q consists of the four electron positions, each a list of three Cartesian coordinates.

The solution to (1) is subject to two boundary conditions: that whenever one or more electrons approach infinity, and that must change sign whenever two electrons of the same spin are interchanged, in accord with the Pauli exclusion principle [3].

There are several techniques for finding an approximate solution to (1), which we denote by , to distinguish it from the exact solution . Monte Carlo is a method that “borrows” one of these solutions, called, in this context, a trial function [1], and seeks to improve its accuracy. This is done by generating, with the help of , a random statistical sample, called an ensemble in this context, representing the exact solution. Based on this ensemble of electron locations, or configurations, one can then easily find, within the “standard” statistical error, the value of the molecule’s ground-state energy, and related properties such as dipole moment and polarizability, etc.

To build a trial solution for LiH, we start with two molecular orbitals, linear combinations of four simple atomic orbitals [2].

The argument q is a list of three coordinates that describe a location of a single electron. The parameters of the MO functions and those of J in the next expression have been obtained by minimizing the variational energy (7).

Secondly, we define the so-called Jastrow function [2], its arguments being locations of each pair of electrons.

The resulting trial function has the following form.

The two factors are determinants of the two molecular orbitals, each evaluated at the location of the four electrons. Electrons 1 and 2 have spin “up,” while 3 and 4 have spin “down.” The remaining factors are Jastrow functions for each pair of opposite-spin electrons.

Next we estimate the ground-state energy of LiH based on this trial function.

Variational Estimate of the Smallest Eigenvalue

Let us rewrite (1) more compactly as

(3)

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

(4)

The integration is over the three coordinates of each of the four electrons, altogether a 12-dimensional problem—no mean task—and is any trial solution to (1). 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 (4) 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. Well, Monte Carlo has the answer.

Let us first define the so-called drift function by

(5)

and local energy by

(6)

Here are the corresponding commands.

The coordinates of the four electrons are now called , understanding that these are the , , coordinates of the first electron, followed by the , , coordinates of the second electron, etc. The left-hand side of (4) can now be rewritten as

(7)

Next, we randomly generate a large sample of 1000 configurations of values denoted collectively as Ro and compute the corresponding , , and .

By averaging the 1000 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 (7).

To fix this, we move each configuration to a new location, specified by

(8)

where is the drift function evaluated at the old location , 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: the so-called grand mean).

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 [3]: for each proposed move (8) we compute a scalar quantity

(9)

where the subscripts and mean that and have been evaluated at the new or old location, respectively. The move is then accepted with a probability equal to . When , 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 , when this number exceeds a certain value, such as 10. 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 (confirming that 10 consecutive rejections was a good indication of a “stuck” configuration).

The following program carries this out (its execution will take a minute or two).

Note that moni keeps track, for each configuration, of how many of its last consecutive proposed moves have been rejected, its value being reset to 0 as soon as a move is accepted. The program returns (in res) the value of all sample averages of the local energy , together with the average acceptance rate.

This displays the former.

We see that about 50 iterations are necessary to reach equilibrium.

To get an accurate estimate of , we repeat the simulation with substantially more iterations, changing the Do loop’s “count” from 60 to 1000, and then computing the grand mean of the values by . In our case, this yields atomic units, with an average acceptance rate of about 85%.

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 . The “exact” ground-state energy of LiH is atomic units. 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.

Monte Carlo Estimate of the Smallest Eigenvalue

When, in (7), we replace by , the expression then yields “nearly” the exact value of , subject only to a small nodal error [1]. 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

(10)

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.

1. Occasionally (e.g., 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.
2. 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.

This can all be achieved by the following simple modifications of the program from the previous section. Monte Carlo techniques in general require a long time to execute (this one may take several hours).

For this reason, we have made it (and the subsequent command, which processes its output) non-evaluatable.

A ListPlot of the iteration averages of will show that equilibration now takes many more steps (about 500, 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.

Similarly, we can produce six such values with and , calling them r050 and r075, respectively.

It is now easy to find the resulting intercept.

This yields the value of for the corresponding intercept. This is in reasonable agreement, in view of the nodal error, with the exact value of atomic units.

This visualizes the regression fit.

In a follow-up article, we will show how this procedure can be extended to estimate other significant molecular properties, including geometry and polarizability, etc. and how to optimize parameters of a trial function, to make the Monte Carlo method more “self-sufficient.”

References

[1]	J. B. Anderson, “Quantum Chemistry by Random Walk: Higher Accuracy,” The Journal of Chemical Physics, 73(8), 1980 pp. 3897-3899. doi:10.1063/1.440575.
[2]	P. J. Reynolds, D. M. Ceperley, B. J. Alder, and W. A. Lester, Jr., “Fixed-Node Quantum Monte Carlo for Molecules,” The Journal of Chemical Physics, 77(12), 1982 pp. 5593-5603. doi:10.1063/1.443766.
[3]	D. M. Ceperley and B. J. Alder, “Quantum Monte Carlo,” Science, 231(4738), 1986 pp. 555-560. doi:10.1126/science.231.4738.555.
J. Vrbik, “Monte Carlo Simulation of Simple Molecules,” The Mathematica Journal, 2011. dx.doi.org/doi:10.3888/tmj.13-5.

About the Author

Jan Vrbik
Department of Mathematics, Brock University
500 Glenridge Ave., St. Catharines
Ontario, Canada, L2S 3A1
jvrbik@brocku.ca