The Full Monte Carlo Code


The purpose of this work is to evaluate Gaussian integrals both analytically and by Biased Selection Monte Carlo methods in order to illustrate and test the Monte-Carlo methods.


The biased selection method began as a method to evaluate partition functions[1].  By using it to estimate and minimize the variance in the local energy it can be used to solve the Schroedinger and Dirac equations.  It has been used extensively to estimate quantum mechanical expectation values.  At its heart it is an integration scheme which should seriously be considered for evaluating any integral of more than about 6 dimensions.  The best description of the method in its current form is in the book Recent Advances in Quantum Monte Carlo Methods[2]In the following sections our implementation of information sampling know as biased as random, the mixed extension that allows one to choose with respect to two centers and  the single_average method which allows for multistage selection are described and tested against analytical evaluations of Gaussian integrals. 

The integrals

The general form of the integrals tested is .

The analytical evaluation of this integral is discussed in analytical4.doc .htm.

The first integral is normalization Integral of the form


Biased As Random

                This is importance sampling[3] which is nothing more than a change of variables so that in (2)

 A new variable is given by


So that

  . (4)

This then makes


  It is convenient to define



So that all integrals are from 0 to 1 yielding


The Monte Carlo estimate of I is


The error estimate is given by


The weighting function g(r) was one left over from an atomic calculation on neon.

Figure 1 Guiding function used in this work with r2 exp(-r) for comparison purposes.

Figure 2 Integral of  g(r) versus r

                Two important items need to be mentioned here.  The actual g(r) is the straight line connecting the points representing it, which is exactly what is plotted in figure 1.  The subroutines setup and bli enable this to be accomplished.  A set of points to best represent the function desired a guiding function is found by minimizing the difference between the function and a linear interpolation to the function using 65 points.  These values of ri, g(ri) and the integrals from 0 to g(ri) are stored, setup.for.  Since the actual guiding function is the straight line connecting these points the integral of g(r) is exactly the trap rule up to each point and exactly quadratic between points.  A random y is picked and multiplied times the largest integral of g(r).  This value is then used in a quadratic equation to find the value of r(y) as indicated in figure 2. The details are shown in testing\BiasAsR.for. In the six electron example of (1), the above was applied to each electron, so that


And then (1) is evaluated as a single 18 dimensional integral so that


                In the last step in (12), the rki refers to a r for each electron





Err-errfrom large/Err

R=(Γ-Γexact)/errfrom large


21794 +- 3462





23942 +- 1848





27320 +- 1152





28479 +-   617





29662 +-   327





29802 +-   163





29703 +-     81





29771 +-     41





29800 +-     20




                With this many random numbers, the "randomness" of the set can be a problem for some calculations[4].  The method used here, first suggested by MacLaren and Marsaglia[5], is to insert numbers into a table with a multiplicative congruential method with one multiplier and seed and to extract them from the table with a multiplicative congruential method with a different multiplier and seed.  These are in Fortran as testing\random.for and in C as the pair C\random.cpp, and C\random.h.

Optimizing the Guiding Function

                A detailed account of optimizing the Guiding function for wave function minimization given by Alexander et al[6] can also be used for simple Monte Carlo evaluations.  The variance as given by equation (10) can is explicitly given by


The probability w serves to take this from a sum to an integral, so this can also be written as


At first glance, it seems that something must be wrong owing to the fact that the first term in the variance still contains the weighting function.  This is not wrong, the standard deviation per configuration depends on the method used to select the configurations.  As shown in reference 4, the practical implementation of this is to imagine a wi that was used to select a set of configurations and a  which can be varied to minimize the variance with respect to the set of points previously selected.   Thus equation (13) becomes


which can be minimized with respect to parameters in 

The electron-electron guiding function

                The expectation value of 1/rij is dominated by the low probability of finding ri next to rj, so it does not actually give the problems that might be expected from its singular nature.  However, just as there is an r2 in equation (6) that makes evaluation of singular nuclear electron potentials very straightforward, it is possible to put an rij2 in the weights, if configurations are selected with respect to the anti-parallel electrons.  In order to do this it is necessary to recognize that the weights in equation (6) are the probabilities relative to uniform and random of having chosen the points in the manner that we did.  As first proved in equation 28 of an early paper on the He-He potential[7] by imagining a large lattice of T points and taking the appropriate limits.  If there are two ways of choosing a point x that proceed by probabilities w1(x) and by w2(x), these two probabilities can be replaced by the average probability


This means that for any integral , so that w is first averaged over all the ways of selecting the configuration and then the inverse of this average multiplied by f is averaged.  If it is not possible to get to x with one guiding function, but is with another, obviously the non-contributor puts a zero in equation (16). In the case of the electron-electron guiding function for carbon, what is desired is a guiding function that explicitly removes the 1/rij between anti-parallel electrons with very little other effect.  This is described in more detail in eegfun.doc .htm and the code is in xnefin.for.  The spin up electrons have been selected before the first call to xnefin.  For testing purposes the probability of selecting the next spin down electron with respect to a spin up electron, pele, has been set in a data statement to 0.5.  On entry to the code a random number is used to decide if the selection is to be with the nuclear electron guiding function or with the electron-electron guiding function.  Then the electron's location is picked with the appropriate guiding function.  The average weight is given by the probability of selecting the electron with respect to the nuclear electron guiding function plus the probability of selecting with respect to any of the up electrons (17)

The following were run with pele=0.9


Value of integral





25900 +- 9827





31214 +- 2845





28957 +- 731




Choosing from a "single average"

                It is possible to select an important configuration from a subset of M configurations[8]

1. Generate M sets of up electron coordinates xiup and weights wiup values using the biased as random method above. 

2. For each set of xup values calculate the secondary weight function  

3. Generate a random number R such that 0 < R <  

4 Select a single configuration from the original set of M by finding the one at which  first becomes larger than R.  The probability of generating xjup  is given by



The Fortran code config.for organizes this.  It begins by getting Nup  electron coordinates from the biased as random generator xnfind.  These are saved in N temporary locations.  The determinant is calculated by a separate subroutine WUP.  The value of phi2(k) for each determinant and the sum up to k for each determinant is saved.  A random number is generated and locate called to determine which of the configurations is to actually be used.  Then the electron set and its total weight are set for the spin up electrons.  The locations of the spin down electrons are found using xnefin which a probability set in a data statement of finding each spin down electron with respect to the spin up electrons rather than with respect to the nucleus.  This is also repeated M times and the final w for this is multiplied times the w for the spin up electrons to compute the entire w.  In the following tests,  for the spin up electrons and  for the spin down electrons.  It is expected to gain a small amount due to the fact that the spin down electrons are placed preferentially with those spin up electrons that contribute most significantly to the integral. The functions phi used in the single averaging are calculated in deter.for.

                In the first test I use pele=0, so that only the single averaging is being tested.  The random seeds for all of the following tests are 45247, 89139.







29539 +- 162

- 1.7



29596 +- 117

- 1.8



29893 +- 103




29846 +- 100




29954 +- 100




29870 +- 102


                In the second test I use M=64 and 16384 configurations, but vary pele from 0.1 to 0.999





29722 +-   102

- 0.9


29898 +-   112



29651 +-   188

- 0.8


30249 +- 1021



30986 +- 5730



  5096 +-   551



                The standard deviations increase rapidly beyond 0.7.  This comes about because there are regions of the 6 dimensional space in which all electrons are separated by more than 0.5 aB.  The electron electron selection does not sample this region of space.  Thus this part of the integration region is sampled only by those rare events where the two spin down electrons are both chosen with respect to the nuclear electron guiding function and in addition are selected more than 0.5 aB from all of the spin up electrons.  When one of these events happens, the probability w of its happening is small relative to the value of the integrand and thus 3/w for this point becomes large enough to compensate for the usually missed region.  This large value causes the square to increase even more rapidly and gives rise to the increasing standard deviations from pele = 0.95 to 0.99.  At pele = 0.999, there are regions that are simply not sampled.  Since the bulk of the regions that are sampled appear correct, a standard deviation is given that compounds the problem by seeming to be accurate.  This has been called the epsilon connection where a region of space with probability epsilon of being selected actually contributes a value of R times epsilon which is significant.

The rij tests.doc .htm

[1] R. L. Coldwell,"Monte Carlo Evaluation of the Partition Function for a Hard Disk System", Phys. Rev. A 7, 270 (1973).

[2] “Atomic Calculations Using Variational Monte-Carlo” , S.A. Alexander and R.L. Coldwell,  pp 39-64 in  "Recent Advances in Quantum Monte Carlo Methods" edited          by  W.A. Lester,World Scientific Publishing Co, Singapore. (1997)

[3] J. M. Hammersley and D. C. Handscomb, Monte Carlo Methods, Methune & Co. Ltd, London, John Wiley & Sons Inc, New York.  pp. 57 -59

[4] R. L. Coldwell, "Correlational Defects in the Standard IBM 360 Random Number Generator and the Classical Ideal Gas Correlation Function", J.         Computational Phys. 14, 223 (1974).

[5] M. D. MacLaren and G. Marsaglia, J. Assoc. Comput. Mach. 12 (1965),83

[6] S.A. Alexander, R.L. Coldwell and J.D. Morgan III, "Guiding Function Optimization in Biased-Selection Monte Carlo Calculations" J. Chem. Phys. 97, 8407 (1992).

[7] R.L. Coldwell and R.E. Lowther,"Monte-Carlo Calculations of the Born-Oppenheimer Potential Between Two Helium Atoms Using Hylleraas Type Electronic Wave Functions", Int. J. of Quantum Chem. S12, 329 (1978).


[8] S. A. Alexander and R.L. Coldwell, "Calculating atomic properties using variational Monte Carlo", J. Chem. Phys. 103, 2572 (1995)