Variational Monte Carlo for Ground State Energy

Recommended resources

The freely available resources this lesson will direct you to.

Hello! Welcome back to our study of numerical methods in quantum mechanics.

In our last lesson, we mastered the Metropolis-Hastings algorithm, a powerful technique for sampling from complex, high-dimensional classical probability distributions like the Boltzmann distribution for the Ising model. You learned how to construct a "smart" random walk that explores the state space according to a target probability, and you understood the core principles of detailed balance, acceptance rates, and burn-in.

Today, we will apply this powerful sampling machinery directly to a quantum problem. We'll combine the variational principle—a concept you're already familiar with from your coursework—with the Monte Carlo methods we just learned. This fusion creates one of the workhorse algorithms of computational many-body physics: Variational Monte Carlo (VMC).

Your learning outcome is to apply Variational Monte Carlo to calculate the ground state energy of a simple quantum system. By the end of this lesson, you will understand the VMC algorithm and be able to implement it to find a variational estimate for the ground state energy.

1. The Variational Principle Meets Monte Carlo

Let's begin with the familiar variational principle. For any trial wavefunction $Ψ_{T} (R)$ , where $R$ represents the coordinates of all particles in the system, the expectation value of the Hamiltonian $H^$ provides an upper bound to the true ground state energy $E_{0}$ :

$E [Ψ_{T}] = \frac{⟨ Ψ _{T} ∣ H ^ ∣ Ψ _{T} ⟩}{⟨ Ψ _{T} ∣ Ψ _{T} ⟩} \geq E_{0}$

When we write this out as an integral over the $3 N$ -dimensional configuration space, it becomes:

$E [Ψ_{T}] = \frac{\int Ψ _{T}^{*} ( R ) H ^ Ψ _{T} ( R ) d R}{\int ∣ Ψ _{T} ( R ) ∣ ^{2} d R}$

As we discussed in the last lesson, evaluating such multi-dimensional integrals with standard grid-based methods is impossible for more than a few particles due to the curse of dimensionality. This is where Monte Carlo integration becomes essential.

The key insight of VMC is to reinterpret this integral as the expectation value of a quantity called the local energy, $E_{L} (R)$ , over a specific probability distribution.

First, let's define the local energy:
$E_{L} (R) = \frac{H ^ Ψ _{T} ( R )}{Ψ _{T} ( R )}$
The local energy is the result of applying the Hamiltonian to the trial wavefunction at a single point in configuration space, $R$ , divided by the value of the wavefunction at that point.

Now, we can rewrite the energy expectation value:
$E [Ψ_{T}] = \frac{\int ∣ Ψ _{T} ( R ) ∣ ^{2} ( \frac{H ^ Ψ _{T} ( R )}{Ψ _{T} ( R )} ) d R}{\int ∣ Ψ _{T} ( R ) ∣ ^{2} d R} = \int E_{L} (R) P (R) d R$
where we've defined the probability distribution function (PDF) as $P (R) = \frac{∣ Ψ _{T} ( R ) ∣ ^{2}}{\int ∣ Ψ _{T} ( R ) ∣ ^{2} d R}$ .

This final expression, $⟨ E ⟩ = ⟨ E_{L} ⟩_{P (R)}$ , is exactly what we need. It tells us that the variational energy is simply the average of the local energy, sampled from the probability distribution defined by the square of our trial wavefunction, $∣ Ψ_{T} ∣^{2}$ .

And how do we sample from a complicated distribution like $P (R)$ ? With the Metropolis algorithm we learned in the last lesson!

5. Variational Monte Carlo methods

The following text, 'Variational Monte Carlo methods', provides an excellent and concise summary of this theoretical foundation. It formally connects the variational principle to the Monte Carlo framework.

Please read the section 'Quantum Monte Carlo Motivation'. Focus on understanding the definitions of the local energy and how the problem of calculating the ground state energy is transformed into computing the expectation value of this local energy.

2. The Variational Monte Carlo Algorithm

We can now lay out the complete VMC algorithm. Suppose we have a trial wavefunction $Ψ_{T} (R, α)$ that depends on a set of variational parameters $α$ . Our goal is to find the parameters $α$ that minimize the energy $E [Ψ_{T}]$ .

For a fixed set of parameters $α$ , the VMC procedure is:

Initialization: Place the $N$ particles (our "walkers") at some initial random configuration $R_{0}$ .
Markov Chain Generation: Use the Metropolis algorithm to generate a sequence of $M$ configurations $R_{1}, R_{2}, \dots, R_{M}$ that are sampled from the distribution $P (R) = ∣ Ψ_{T} (R) ∣^{2}$ .
- For each step $t$ , propose a move from $R_{t}$ to a new configuration $R^{'}$ (e.g., by randomly displacing one particle).
- Calculate the Metropolis acceptance ratio: $A = \frac{P ( R ^{'} )}{P ( R _{t} )} = \frac{∣ Ψ _{T} ( R ^{'} ) ∣ ^{2}}{∣ Ψ _{T} ( R _{t} ) ∣ ^{2}}$ .
- Accept the move with probability $min (1, A)$ . If rejected, the next state in the chain is just $R_{t + 1} = R_{t}$ .
Local Energy Calculation: For each configuration $R_{i}$ in the chain (after an initial burn-in period), calculate the local energy $E_{L} (R_{i}) = \frac{H ^ Ψ _{T} ( R _{i} )}{Ψ _{T} ( R _{i} )}$ . The Hamiltonian is typically $H^= T^+ V^$ , so the local energy is:
$E_{L} (R) = \frac{- \frac{ℏ ^{2}}{2 m} \sum _{j} \nabla _{j}^{2} Ψ _{T} ( R )}{Ψ _{T} ( R )} + V (R)$
Note that the potential energy $V (R)$ is usually simple to calculate. The kinetic energy term requires calculating the Laplacian of the trial wavefunction, which must be derived analytically beforehand.
Energy Estimation: The variational energy is estimated as the average of the local energies over the generated samples:
$E (α) \approx E ˉ_{L} = \frac{1}{M} \sum_{i = 1}^{M} E_{L} (R_{i})$
The statistical error of this estimate scales as $σ / M$ , where $σ$ is the standard deviation of the local energy values.
Optimization: Repeat steps 2-4 for different values of the parameters $α$ to find the minimum energy.

Python Code for Variational Monte Carlo Setup — An example of the core Python functions for a VMC simulation of a 1D harmonic oscillator. You can see the definitions for the trial wavefunction (`psi_trial`) and the local energy (`local_energy`), which includes the analytically derived kinetic energy term and the potential energy.

A powerful feature of VMC is related to the variance of the local energy. If our trial wavefunction $Ψ_{T}$ happens to be the exact ground state eigenfunction $Ψ_{0}$ , then $H^Ψ_{0} = E_{0} Ψ_{0}$ . In this case, the local energy is constant everywhere:
$E_{L} (R) = \frac{H ^ Ψ _{0} ( R )}{Ψ _{0} ( R )} = \frac{E _{0} Ψ _{0} ( R )}{Ψ _{0} ( R )} = E_{0}$
This means the variance of the local energy, $σ^{2} = ⟨(E_{L} - ⟨ E_{L} ⟩)^{2} ⟩$ , will be zero. A smaller variance generally indicates a better trial wavefunction. Therefore, a common strategy is to optimize $α$ by minimizing not just the energy, but a combination of the energy and its variance.

Test your understanding!

You are performing a VMC calculation for the 1D simple harmonic oscillator ( $V (x) = \frac{1}{2} m ω^{2} x^{2}$ ). You choose the trial wavefunction $Ψ_{T} (x) = e^{- α x^{2}}$ . The exact ground state wavefunction has this form with $α_{e x a c t} = \frac{mω}{2ℏ}$ .

What do you expect will happen to the list of local energies you calculate during your Monte Carlo run as your variational parameter $α$ gets very close to $α_{e x a c t}$ ?

Show answer

As $α$ approaches $α_{e x a c t}$ , your trial wavefunction gets closer to the true ground state eigenfunction. Consequently, the local energy $E_{L} (x)$ will become nearly constant for all positions $x$ . The list of local energies calculated at each step of your Monte Carlo walk will show very little fluctuation, and their standard deviation (and variance) will approach zero. This "zero variance principle" is a unique and powerful feature of VMC.

3. Implementation: VMC for the Hydrogen Molecule

Now for the practical part. Your background in Python and Numpy makes you well-equipped to build a VMC simulation. We will follow a guided tutorial to calculate the ground state energy of a simple diatomic molecule: H₂.

The system consists of two protons (which we'll fix in space) and two electrons. Our configuration vector is $R = (r_{1}, r_{2})$ . The Hamiltonian in atomic units ( $ℏ = m_{e} = e = 1$ ) for protons fixed at positions $R_{A}$ and $R_{B}$ is:
$H^= Kinetic (e⁻) - \frac{1}{2} \nabla_{1}^{2} - \frac{1}{2} \nabla_{2}^{2} Potential (e⁻-p⁺) - \frac{1}{∣ r _{1} - R _{A} ∣} - \frac{1}{∣ r _{1} - R _{B} ∣} - \frac{1}{∣ r _{2} - R _{A} ∣} - \frac{1}{∣ r _{2} - R _{B} ∣} + Potential (e⁻-e⁻) \frac{1}{∣ r _{1} - r _{2} ∣} + Constant (p⁺-p⁺) \frac{1}{∣ R _{A} - R _{B} ∣}$

We will use a very simple (and frankly, not very accurate, but good for learning) trial wavefunction that ignores electron correlation:
$Ψ_{T} (r_{1}, r_{2}; α) = e^{- α (∣ r_{1} ∣^{2} + ∣ r_{2} ∣^{2})}$
This function has one variational parameter, $α$ .

In this tutorial, we will write a variational Monte Carlo code ...

The following resource provides a step-by-step guide to writing a VMC code in Python to solve this exact problem. We will work through its main sections to build our simulation.

This is a hands-on coding task. Follow the tutorial to build your VMC program. Read 'Trial Wave functions for a hydrogen molecule' and 'Building the WaveFunction'. Implement the psi function in Python and verify it with the provided test case. Read 'Building a Monte Carlo Code'. Implement the main VMC loop that proposes moves and uses the Metropolis criterion to accept/reject them. Initially, just compute the acceptance ratio and check it against the expected value. Read 'Computing the local energy'. This is the most physics-intensive step. You need to analytically derive the expression for the local energy for our given Hamiltonian and trial wavefunction. Then, implement the LocalEnergy function and verify it with its test case. Integrate the pieces. Add the local energy calculation to your main VMC loop, as described in the paragraph starting 'Now add the calculation of the local energy...'. Collect a list of energies and compute their mean and standard error.

This implementation will give you a working VMC code. To make it more robust and extensible—essential for tackling the many-body problems you're interested in—it's good practice to use an object-oriented approach. The tutorial also discusses an important optimization called importance sampling, which modifies the proposal step of the Metropolis algorithm to achieve much higher acceptance rates.

In this tutorial, we will write a variational Monte Carlo code ...

To see how to structure your code more effectively and how to improve the sampling efficiency, continue with the same tutorial.

Read the sections 'Classes all around' and 'A better acceptance ratio'. You don't need to re-implement everything right now, but focus on understanding: How a class can encapsulate the wavefunction and its parameters (alpha, ion positions). The concept of the 'quantum force' (F(R) = abla \Psi / \Psi) and how it's used to bias the random walk towards regions of higher probability, improving efficiency. This connects to the Metropolis-Hastings algorithm (asymmetric proposal) from our last lesson.

The concept of using the "quantum force" for importance sampling is a standard technique in advanced QMC methods, connecting the random walk to a drift-diffusion process described by the Langevin equation. The text 5. Variational Monte Carlo methods (resource LINK) provides a deeper theoretical treatment of this in its section on 'Importance sampling' if you wish to explore it further.

Conclusion

Congratulations! You have now bridged the gap between the variational principle, a cornerstone of quantum theory, and the powerful Monte Carlo methods from computational science.

Key Takeaways:

Variational Monte Carlo (VMC) transforms the problem of finding the ground state energy into a statistical sampling problem.
The method relies on calculating the local energy $E_{L} (R) = H^Ψ_{T} / Ψ_{T}$ at configurations $R$ sampled from the probability distribution $∣ Ψ_{T} (R) ∣^{2}$ .
The Metropolis algorithm is the engine used to generate the configurations that correctly sample the $∣ Ψ_{T} ∣^{2}$ distribution.
The quality of a trial wavefunction can be judged not only by the resulting energy but also by the variance of the local energy, which is zero for an exact eigenstate.
Implementing VMC involves a blend of analytical work (deriving the local energy) and programming (implementing the Metropolis sampler and statistical averaging).

Preview of the Next Lesson:

While VMC is powerful, its accuracy is fundamentally limited by the quality of the trial wavefunction you choose. How do we deal with systems where we don't have a good guess for $Ψ_{T}$ ? In the next lesson, we will explore Diffusion Monte Carlo (DMC), a more advanced QMC method that can, in principle, find the exact ground state energy, stochastically evolving an ensemble of walkers in imaginary time to project out the ground state component. VMC often serves as the crucial first step for initializing a DMC simulation.