Quantum Monte Carlo

From Wikipedia, the free encyclopedia

Electronic structure methods
v  d  e
Tight binding
Hartree-Fock
Configuration interaction
Coupled cluster
Multi-configurational self-consistent field
Density functional theory
Quantum Monte Carlo
Density matrix renormalization group

Quantum Monte Carlo is a large class of computer algorithms that simulate quantum systems with the idea of solving the many-body problem. They use, in one way or another, the Monte Carlo method to handle the many dimensional integrals that arise. Quantum Monte Carlo allows a direct representation of many-body effects in the wavefunction, at the cost of statistical uncertainty that can be reduced with more simultation time. For bosons, there exist numerically exact and polynomial-scaling algorithms. For fermions, there exist very good approximations and numerically exact exponentially scaling quantum Monte Carlo algorithms, but none that are both.

Contents

[edit] Background

In principle, any physical system can be described by the many-body Schrödinger equation, as long as the constituent particles are not moving 'too' fast; that is, they are not moving near the speed of light. This includes the electrons in almost every material in the world, so if we could solve the Schrödinger equation, we could predict the behavior of any electronic system, which has important applications in fields from computers to biology. This also includes the nuclei in Bose-Einstein condensates and superfluids like liquid helium. The difficulty is that the Schrödinger equation involves a function of three times the number of particles(in 3 dimensions), and is difficult(and impossible in the case of fermions) to solve in a reasonable amount of time. Traditionally, theorists have approximated the many-body wave function as a separable function of one-body wave functions: \Psi(x_1,x_2,\dots,x_n)=f(\Phi_1(x_1),\Phi_2(x_1),\dots, \Phi_n(x_1);\Phi_1(x_2) \Phi_2(x_2),\dots), for an example, see Hartree-Fock theory. This kind of formulation either limits the possible wave functions, as in the case of Hartree-Fock, or converges very slowly, as in configuration interaction.

Quantum Monte Carlo is a way around these problems, which allows us to work with the full many-body wave function directly. Most methods aim at computing the ground state wave function of the system, with the exception of Path Integral Monte Carlo and finite-temperature auxiliary field Monte Carlo, which calculate the density matrix. What follows is a sampling of some of the quantum Monte Carlo flavors, which approach the problem in different ways.

[edit] Variational Monte Carlo

Variational Monte Carlo takes advantage of the Variational method in quantum mechanics to approximate the ground state. The expectation value necessary can be written in the x representation as \frac{\langle \Psi(a) | H | \Psi(a) \rangle} {\langle \Psi(a) |  \Psi(a) \rangle } = \frac{\int \Psi(X,a)^2 \frac{H\Psi(X,a)}{\Psi(X,a)} dX} { \int \Psi(X,a)^2 dX}. Following the Monte Carlo method for evaluating integrals, we can interpret \frac{ \Psi(X,a)^2 } { \int \Psi(X,a) dX } as a probability distribution function, sample it, and evaluate the energy expectation value E(a) as the average of the local function \frac{H\Psi(X,a)}{\Psi(X,a)}, and minimize E(a).

This is no different from any other variational method, except that since the many-dimensional integrals are evaluated numerically, we only need to calculate the value of the possibly very complicated wave function, which gives a large amount of flexibility to the method. One of the largest gains in accuracy over writing the wave function separably comes from the introduction of the so-called Jastrow factor, where the wave function is written as exp(\sum{u(r_{ij})}), where rij is the distance between a pair of quantum particles. With this factor, we can explicitly account for particle-particle correlation, but the many-body integral becomes unseparable, so Monte Carlo is the only way to evaluate it efficiently. In chemical systems, slightly more sophisticated versions of this factor can obtain 80-90% of the correlation energy (see electronic correlation) with less than 30 parameters. In comparison, a configuration interaction calculation may require around 50,000 parameters to reach that accuracy, although it depends greatly on the particular case being considered. In addition, VMC usually scales as a small power of the number of particles in the simulation, usually something like N2-4 for calculation of the energy expectation value, depending on the form of the wave function.

[edit] Wave Function Optimization in VMC

QMC calculations crucially depend on the quality of the trial-function, and so it is essential to have an optimized wave-function as close as possible to the ground state. The problem of function optimization is a very important research topic in numerical simulation. In QMC, in addition to the usual difficulties to find the minimum of multidimensional parametric function, the statistical noise is present in the estimate of the cost function (usually the energy), and its derivatives , required for an efficient optimization.

Different cost functions and different strategies were used to optimize a many-body trial-function. Usually three cost functions were used in QMC optimization energy, variance or a linear combination of them. In this thesis we always used energy optimization. The variance optimization have the advantage to be bounded by below, to be positive defined and its minimum is known, but different authors recently showed that the energy optimization is more effective than the variance one.

There are different motivations for this: first, usually one is interested in the lowest energy rather than in the lowest variance in both variational and diffusion Monte Carlo; second, variance optimization takes many iterations to optimize determinant parameters and often the optimization can get stuck in multiple local minimum and it suffers of the "false convergence" problem; third energy-minimized wave functions on average yield more accurate values of other expectation values than variance minimized wave functions do.

The optimization strategies can be divided into three categories. The first strategy is based on correlated sampling together with deterministic optimization methods. Even if this idea yielded very accurate results for the first-row atoms, this procedure can have problems if parameters affect the nodes, and moreover density ratio of the current and initial trial-function increases exponentially with the size of the system. In the second strategy one use a large bin to evaluate the cost function and its derivatives in such way that the noise can be neglected and deterministic methods can be used.

The third approach, is based on an iterative technique to handle directly with noise functions. The first example of these methods is the so called Stochastic Gradient Approximation (SGA), that was used also for structure optimization. Recently an improved and faster approach of this kind was proposed the so called Stochastic Reconfiguration (SR) method.

[edit] Diffusion Monte Carlo

Diffusion Monte Carlo (DMC) is potentially numerically exact, meaning that it can find the exact ground state energy within a given error for any quantum system. When actually attempting the calculation, one finds that for bosons, the algorithm is scales as a polynomial with the system size, but for fermions, DMC is exponentially scaling with the system size. This makes exact large-scale DMC simulations for fermions impossible; however, with a clever approximation known as fixed-node, very accurate results can be obtained. What follows is an explanation of the basic algorithm, how it works, why fermions cause a problem, and how the fixed-node approximation resolves this problem.

The Projector Method

To motivate the algorithm, let's look at the Schoedinger equation for a particle in some potential in one dimension: i\frac{d\Psi(x,t)}{dt}=\frac{d^2 \Psi(x,t)}{dx^2} + V(x)\Psi(x,t). We can condense the notation a bit by writing it in terms of an operator equation, with H=\frac{d^2 }{dx^2} + V(x). So then we have i\frac{d\Psi(x,t)}{dt}=H\Psi(x,t), where we have to keep in mind that H is an operator, not a simple number or function. There are special functions, called eigenfunctions, for which HΨ(x) = EΨ(x), where E is a number. These functions are special because no matter where we evaluate the action of the H operator on the wave function, we always get the same number E. These functions are called stationary states, because the time derivative at any point x is always the same, so the amplitude of the wave function never changes in time. Since the overall phase of a wave function is not measurable, the system does not change in time.

We are usually interested in the wave function with the lowest energy eigenvalue, the ground state. We're going to write a slightly different version of the Schroedinger equation that will have the same energy eigenvalue, but, instead of being oscillatory, it will be convergent. Here it is: \frac{-d\Psi(x,t)}{dt}=(H-E_0)\Psi(x,t). We've removed the imaginary number from the time derivative and added in a constant offset of E0, which is the ground state energy. We don't actually know the ground state energy, but there will be a way to determine it self-consistently which we'll introduce later. Our modified equation(some people call it the imaginary-time Schrodinger equation) has some nice properties. The first thing to notice is that if we happen to guess the ground state wave function, then HΦ0(x) = E0Φ0(x) and the time derivative is zero. Now suppose that we start with another wave function(Ψ), which is not the ground state but is not orthogonal to it. Then we can write it as a linear sum of eigenfunctions: \Psi=c_0\Phi_0+\sum_{i=1}^\infty c_i\Phi_i Since this is a linear differential equation, we can look at the action of each part separately. We already determined that Φ0 is stationary. Suppose we take Φ1. Since Φ0 is the lowest-energy eigenfunction, the associate eigenvalue of Φ1 satisfies the property E1 > E0. Thus the time derivative of c1 is negative, and will eventually go to zero, leaving us with only the ground state. This observation also gives us a way to determine E0. We watch the amplitude of the wave function as we propagate through time. If it increases, then decrease the estimation of the offset energy. If the amplitude decreases, then increase the estimate of the offset energy.

Stochastic Implementation

Now we have an equation that, as we propagate it forward in time and adjust E0 appropriately, we find the ground state of any given Hamiltonian. This is still a harder problem than classical mechanics, though, because instead of propagating single positions of particles, we must propagate entire functions. In classical mechanics, we could simulate the motion of the particles by setting x(t + τ) = x(t) + τv(t) + 0.5F(t2, if we assume that the force is constant over the time span of τ. For the imaginary time Schroedinger equation, instead, we propagate forward in time using a convolution integral with a special function called a Green's function. So we get \Psi(x,t+\tau)=\int G(x,x',\tau) \Psi(x',t) dx'. Similarly to classical mechanics, we can only propagate for small slices of time; otherwise the Green's function is inaccurate. As the number of particles increases, the dimensionality of the integral increases as well, since we have to integrate over all coordinates of all particles. We can do these integrals by Monte Carlo.

[edit] Flavors of quantum Monte Carlo

[edit] See also

[edit] References

  • Diffusion Monte Carlo
    • R.C. Grimm and R.G. Storer, J. Comput. Phys. 7,134(1971)
    • [1] J. Anderson, J. Chem. Phys. 63, 1499(1975)
    • B.L. Hammond, W.A Lester, Jr. & P.J. Reynolds "Monte Carlo Methods in Ab Initio Quantum Chemistry" (World Scientific, 1994)
  • Path Integral Monte Carlo
    • [2]D.M. Ceperley, Rev. Mod. Phys. 67,279(1995)
  • Reptation Monte Carlo
    • [3] S. Baroni and S. Moroni, Phys. Rev. Lett. 82,4745(1999)
  • Variational Monte Carlo
    • [4] W.L. McMillan, Phys. Rev. 138,A442(1964)
    • [5] D. Ceperley, G.V. Chester and M.H. Kalos, Phys. Rev. B 16,3081(1977)
  • Wave Function Optimization in VMC
    • [6] M. Snajdr. S.M. Rothstein., J. Chem. Phys. 112,4935(2000)
    • [7] D. Bressanini et al., J. Chem. Phys. 116,5345(2002)
    • [8] J. W. Wilkins C.J. Umrigar, K. G. Wilson, Phys. Rev. Lett. 60,1719(1988)
    • [9] P. R. C. Kent R. J. Needs and G. Rajagopal, Phys. Rev. B, 59,12344(1999)
    • [10] Xi Lin Hongkai Zhang and Andrew M. Rappe, J. Chem. Phys., 112,2650(2000)
    • [11] A. Harju, B. Barbiellini, S. Siljamaki, R. M. Nieminen, and G. Ortiz, Phys. Rev. Lett. 79, 1173(1997)
    • [12] S. Tanaka, J. Chem. Phys., 100,7416(1994)
    • [13] M. Casula, C. Attaccalite, and S. Sorella, J. Chem. Phys., 121,7110(2004)

[edit] External links

[edit] Lecture notes

[edit] Computer programs