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 MonteCarlo 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 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
This is importance sampling[3] which is nothing more than a change of variables so that in (2)
A new variable is given by
(3)
So that
. (4)
This then makes
(5).
It is convenient to define
(7)
So that all integrals are from 0 to 1 yielding
(8)
The
(9)
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 r^{2} 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 r_{i}, g(r_{i}) and the integrals from 0 to g(r_{i}) 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
(11).
And then (1) is evaluated as a single 18 dimensional integral so that
In the last step in (12), the r_{ki} refers to a r for each electron
Configurations 
Γ 
R=(ΓΓ_{exact})/err 
Errerr_{from large/Err } 
R=(ΓΓ_{exact})/err_{from large} 
1024 
21794 + 3462 
2.3 
0.328 
1.6 
4096 
23942 + 1848 
3.2 
0.278 
2.3 
16384 
27320 + 1152 
2.2 
0.100 
1.9 
64536 
28479 + 617 
2.2 
0.036 

262144 
29662 + 327 
0.4 
0.022 

1048576 
29802 + 163 
0.0 
0.019 

4194304 
29703 + 81 
1.3 
0.012 

16777216 
29771 + 41 
0.9 
0.025 

67108864 
29800 + 20 
0.4 
0 

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.
A
detailed account of optimizing the Guiding function for wave function
minimization given by Alexander et al[6] can
also be used for simple
The probability w serves to take this from a sum to an integral, so this can also be written as
(14)
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 w_{i} 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
(15)
which can be minimized with respect to parameters in .
The expectation value of 1/r_{ij} is dominated by the low probability of finding r_{i} next to r_{j}, so it does not actually give the problems that might be expected from its singular nature. However, just as there is an r^{2} in equation (6) that makes evaluation of singular nuclear electron potentials very straightforward, it is possible to put an r_{ij}^{2} in the weights, if configurations are selected with respect to the antiparallel 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 HeHe 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 w_{1}(x) and by w_{2}(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 noncontributor puts a zero in equation (16). In the case of the electronelectron guiding function for carbon, what is desired is a guiding function that explicitly removes the 1/r_{ij} between antiparallel 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 electronelectron 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
Configurations 
Value of integral 
Residual 
Seed1 
Seed2 
1600 
25900 + 9827 
0.4 
248 
131 
256000 
31214 + 2845 
+0.5 
248 
131 
4096000 
28957 + 731 
+1.2 
248 
131 
It is possible to select an important configuration from a subset of M configurations[8].
1. Generate M sets of up electron coordinates x_{i}^{up} and weights w_{i}^{up} values using the biased as random method above.
2. For each set of x^{up} 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 x_{j}^{up } is given by
The Fortran code config.for organizes this. It begins by getting N_{up} 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 phi^{2}(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.
Configurations 
M 
Value 
Residual 
1048576 
1 
29539 + 162 
 1.7 
262144 
4 
29596 + 117 
 1.8 
65536 
16 
29893 + 103 
+0.8 
16384 
64 
29846 + 100 
+0.4 
4096 
256 
29954 + 100 
+1.4 
1024 
1024 
29870 + 102 
+0.6 
In the second test I use M=64 and 16384 configurations, but vary pele from 0.1 to 0.999
Pele 
Value 
Residual 
0.1 
29722 + 102 
 0.9 
0.3 
29898 + 112 
+0.8 
0.7 
29651 + 188 
 0.8 
0.95 
30249 + 1021 
+0.4 
0.99 
30986 + 5730 
+0.2 
0.999 
5096 + 551 
44.9 
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 a_{B}. 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 a_{B} 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.
[1] R.
L. Coldwell,"
[2] “Atomic
Calculations Using Variational MonteCarlo” , S.A. Alexander and R.L.
Coldwell, pp 3964 in "Recent Advances in Quantum Monte Carlo
Methods" edited by W.A. Lester,World Scientific Publishing Co,
[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 BiasedSelection Monte Carlo Calculations" J. Chem. Phys. 97, 8407 (1992).
[7] R.L. Coldwell and R.E. Lowther,"MonteCarlo Calculations of the BornOppenheimer 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