1
0
mirror of https://github.com/TREX-CoE/qmc-lttc.git synced 2024-11-04 05:04:01 +01:00
This commit is contained in:
scemama 2021-02-04 10:25:57 +00:00
parent 99539f2445
commit 133703df94

View File

@ -3,7 +3,7 @@
"http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<!-- 2021-02-04 Thu 07:47 -->
<!-- 2021-02-04 Thu 10:25 -->
<meta http-equiv="Content-Type" content="text/html;charset=utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>Quantum Monte Carlo</title>
@ -329,75 +329,92 @@ for the JavaScript code in this tag.
<h2>Table of Contents</h2>
<div id="text-table-of-contents">
<ul>
<li><a href="#org5ef0488">1. Introduction</a>
<li><a href="#org01a29e9">1. Introduction</a>
<ul>
<li><a href="#org900f1d4">1.1. Energy and local energy</a></li>
<li><a href="#orgd30e15c">1.1. Energy and local energy</a></li>
</ul>
</li>
<li><a href="#orgeb42fdf">2. Numerical evaluation of the energy of the hydrogen atom</a>
<li><a href="#org08f395b">2. Numerical evaluation of the energy of the hydrogen atom</a>
<ul>
<li><a href="#org4a6f916">2.1. Local energy</a>
<li><a href="#orgb5021b6">2.1. Local energy</a>
<ul>
<li><a href="#orgd3ba702">2.1.1. Exercise 1</a></li>
<li><a href="#org0939e15">2.1.2. Exercise 2</a></li>
<li><a href="#org8603375">2.1.3. Exercise 3</a></li>
<li><a href="#org4d5cd01">2.1.4. Exercise 4</a></li>
<li><a href="#org3af2b6e">2.1.5. Exercise 5</a></li>
<li><a href="#org493ceaf">2.1.1. Exercise 1</a></li>
<li><a href="#org32df485">2.1.2. Exercise 2</a></li>
<li><a href="#org4c17cf7">2.1.3. Exercise 3</a></li>
<li><a href="#orgc43c8fc">2.1.4. Exercise 4</a></li>
<li><a href="#org0da2733">2.1.5. Exercise 5</a></li>
</ul>
</li>
<li><a href="#org5ab9959">2.2. Plot of the local energy along the \(x\) axis</a>
<li><a href="#org107b849">2.2. Plot of the local energy along the \(x\) axis</a>
<ul>
<li><a href="#orga7ce1d2">2.2.1. Exercise</a></li>
<li><a href="#orgce43b8f">2.2.1. Exercise</a></li>
</ul>
</li>
<li><a href="#orgd0d619f">2.3. Numerical estimation of the energy</a>
<li><a href="#orgcc5a700">2.3. Numerical estimation of the energy</a>
<ul>
<li><a href="#org1229a69">2.3.1. Exercise</a></li>
<li><a href="#org33238fa">2.3.1. Exercise</a></li>
</ul>
</li>
<li><a href="#orgad1ee8c">2.4. Variance of the local energy</a>
<li><a href="#org24aa290">2.4. Variance of the local energy</a>
<ul>
<li><a href="#org0afc4b1">2.4.1. Exercise (optional)</a></li>
<li><a href="#org10f6a19">2.4.2. Exercise</a></li>
<li><a href="#orgbd8786f">2.4.1. Exercise (optional)</a></li>
<li><a href="#orgdf837df">2.4.2. Exercise</a></li>
</ul>
</li>
</ul>
</li>
<li><a href="#org0de2b52">3. Variational Monte Carlo</a>
<li><a href="#orgf24e4c1">3. Variational Monte Carlo</a>
<ul>
<li><a href="#org4abd3b8">3.1. Computation of the statistical error</a>
<li><a href="#org5a7850e">3.1. Computation of the statistical error</a>
<ul>
<li><a href="#org5a1c95b">3.1.1. Exercise</a></li>
<li><a href="#org86ece31">3.1.1. Exercise</a></li>
</ul>
</li>
<li><a href="#org8a01ad5">3.2. Uniform sampling in the box</a>
<li><a href="#orgf3538bc">3.2. Uniform sampling in the box</a>
<ul>
<li><a href="#orgae0d75a">3.2.1. Exercise</a></li>
<li><a href="#org5ea7de2">3.2.1. Exercise</a></li>
</ul>
</li>
<li><a href="#org23b5713">3.3. Metropolis sampling with \(\Psi^2\)</a>
<li><a href="#orgc5a1ec2">3.3. Metropolis sampling with \(\Psi^2\)</a>
<ul>
<li><a href="#org06c5905">3.3.1. Optimal step size</a></li>
<li><a href="#org654010b">3.3.2. Exercise</a></li>
<li><a href="#org5e266bf">3.3.1. Optimal step size</a></li>
<li><a href="#org7580d40">3.3.2. Exercise</a></li>
</ul>
</li>
<li><a href="#org215f308">3.4. Generalized Metropolis algorithm</a>
<li><a href="#orgcaabae7">3.4. Generalized Metropolis algorithm</a>
<ul>
<li><a href="#org31cb975">3.4.1. Gaussian random number generator</a></li>
<li><a href="#org48820db">3.4.2. Exercise 1</a></li>
<li><a href="#org1ce793c">3.4.3. Exercise 2</a></li>
<li><a href="#org5445700">3.4.1. Gaussian random number generator</a></li>
<li><a href="#org1b2244b">3.4.2. Exercise 1</a></li>
<li><a href="#org0368ce7">3.4.3. Exercise 2</a></li>
</ul>
</li>
</ul>
</li>
<li><a href="#orgc4c9346">4. Project</a></li>
<li><a href="#org79362c8">5. Acknowledgments</a></li>
<li><a href="#org82635bc">4. Diffusion Monte Carlo</a>
<ul>
<li><a href="#org1ddc667">4.1. Schrödinger equation in imaginary time</a></li>
<li><a href="#org7695533">4.2. Relation to diffusion</a></li>
<li><a href="#orgff86247">4.3. Importance sampling</a>
<ul>
<li><a href="#orgc628347">4.3.1. Appendix : Details of the Derivation</a></li>
</ul>
</li>
<li><a href="#org429350d">4.4. Pure Diffusion Monte Carlo</a></li>
<li><a href="#org7b3d04e">4.5. Hydrogen atom</a>
<ul>
<li><a href="#org1a74409">4.5.1. Exercise</a></li>
</ul>
</li>
</ul>
</li>
<li><a href="#org60565bd">5. Project</a></li>
<li><a href="#orga85900b">6. Acknowledgments</a></li>
</ul>
</div>
</div>
<div id="outline-container-org5ef0488" class="outline-2">
<h2 id="org5ef0488"><span class="section-number-2">1</span> Introduction</h2>
<div id="outline-container-org01a29e9" class="outline-2">
<h2 id="org01a29e9"><span class="section-number-2">1</span> Introduction</h2>
<div class="outline-text-2" id="text-1">
<p>
This website contains the QMC tutorial of the 2021 LTTC winter school
@ -437,8 +454,8 @@ coordinates, etc).
</p>
</div>
<div id="outline-container-org900f1d4" class="outline-3">
<h3 id="org900f1d4"><span class="section-number-3">1.1</span> Energy and local energy</h3>
<div id="outline-container-orgd30e15c" class="outline-3">
<h3 id="orgd30e15c"><span class="section-number-3">1.1</span> Energy and local energy</h3>
<div class="outline-text-3" id="text-1-1">
<p>
For a given system with Hamiltonian \(\hat{H}\) and wave function \(\Psi\), we define the local energy as
@ -521,8 +538,8 @@ energy computed over these configurations:
</div>
</div>
<div id="outline-container-orgeb42fdf" class="outline-2">
<h2 id="orgeb42fdf"><span class="section-number-2">2</span> Numerical evaluation of the energy of the hydrogen atom</h2>
<div id="outline-container-org08f395b" class="outline-2">
<h2 id="org08f395b"><span class="section-number-2">2</span> Numerical evaluation of the energy of the hydrogen atom</h2>
<div class="outline-text-2" id="text-2">
<p>
In this section, we consider the hydrogen atom with the following
@ -551,8 +568,8 @@ To do that, we will compute the local energy and check whether it is constant.
</p>
</div>
<div id="outline-container-org4a6f916" class="outline-3">
<h3 id="org4a6f916"><span class="section-number-3">2.1</span> Local energy</h3>
<div id="outline-container-orgb5021b6" class="outline-3">
<h3 id="orgb5021b6"><span class="section-number-3">2.1</span> Local energy</h3>
<div class="outline-text-3" id="text-2-1">
<p>
You will now program all quantities needed to compute the local energy of the H atom for the given wave function.
@ -579,8 +596,8 @@ to catch the error.
</div>
</div>
<div id="outline-container-orgd3ba702" class="outline-4">
<h4 id="orgd3ba702"><span class="section-number-4">2.1.1</span> Exercise 1</h4>
<div id="outline-container-org493ceaf" class="outline-4">
<h4 id="org493ceaf"><span class="section-number-4">2.1.1</span> Exercise 1</h4>
<div class="outline-text-4" id="text-2-1-1">
<div class="exercise">
<p>
@ -626,8 +643,8 @@ and returns the potential.
</div>
</div>
<div id="outline-container-org0939e15" class="outline-4">
<h4 id="org0939e15"><span class="section-number-4">2.1.2</span> Exercise 2</h4>
<div id="outline-container-org32df485" class="outline-4">
<h4 id="org32df485"><span class="section-number-4">2.1.2</span> Exercise 2</h4>
<div class="outline-text-4" id="text-2-1-2">
<div class="exercise">
<p>
@ -663,8 +680,8 @@ input arguments, and returns a scalar.
</div>
</div>
<div id="outline-container-org8603375" class="outline-4">
<h4 id="org8603375"><span class="section-number-4">2.1.3</span> Exercise 3</h4>
<div id="outline-container-org4c17cf7" class="outline-4">
<h4 id="org4c17cf7"><span class="section-number-4">2.1.3</span> Exercise 3</h4>
<div class="outline-text-4" id="text-2-1-3">
<div class="exercise">
<p>
@ -746,8 +763,8 @@ Therefore, the local kinetic energy is
</div>
</div>
<div id="outline-container-org4d5cd01" class="outline-4">
<h4 id="org4d5cd01"><span class="section-number-4">2.1.4</span> Exercise 4</h4>
<div id="outline-container-orgc43c8fc" class="outline-4">
<h4 id="orgc43c8fc"><span class="section-number-4">2.1.4</span> Exercise 4</h4>
<div class="outline-text-4" id="text-2-1-4">
<div class="exercise">
<p>
@ -807,8 +824,8 @@ are calling is yours.
</div>
</div>
<div id="outline-container-org3af2b6e" class="outline-4">
<h4 id="org3af2b6e"><span class="section-number-4">2.1.5</span> Exercise 5</h4>
<div id="outline-container-org0da2733" class="outline-4">
<h4 id="org0da2733"><span class="section-number-4">2.1.5</span> Exercise 5</h4>
<div class="outline-text-4" id="text-2-1-5">
<div class="exercise">
<p>
@ -820,8 +837,8 @@ Find the theoretical value of \(a\) for which \(\Psi\) is an eigenfunction of \(
</div>
</div>
<div id="outline-container-org5ab9959" class="outline-3">
<h3 id="org5ab9959"><span class="section-number-3">2.2</span> Plot of the local energy along the \(x\) axis</h3>
<div id="outline-container-org107b849" class="outline-3">
<h3 id="org107b849"><span class="section-number-3">2.2</span> Plot of the local energy along the \(x\) axis</h3>
<div class="outline-text-3" id="text-2-2">
<p>
The program you will write in this section will be written in
@ -852,8 +869,8 @@ In Fortran, you will need to compile all the source files together:
</div>
</div>
<div id="outline-container-orga7ce1d2" class="outline-4">
<h4 id="orga7ce1d2"><span class="section-number-4">2.2.1</span> Exercise</h4>
<div id="outline-container-orgce43b8f" class="outline-4">
<h4 id="orgce43b8f"><span class="section-number-4">2.2.1</span> Exercise</h4>
<div class="outline-text-4" id="text-2-2-1">
<div class="exercise">
<p>
@ -949,8 +966,8 @@ plot './data' index 0 using 1:2 with lines title 'a=0.1', \
</div>
</div>
<div id="outline-container-orgd0d619f" class="outline-3">
<h3 id="orgd0d619f"><span class="section-number-3">2.3</span> Numerical estimation of the energy</h3>
<div id="outline-container-orgcc5a700" class="outline-3">
<h3 id="orgcc5a700"><span class="section-number-3">2.3</span> Numerical estimation of the energy</h3>
<div class="outline-text-3" id="text-2-3">
<p>
If the space is discretized in small volume elements \(\mathbf{r}_i\)
@ -980,8 +997,8 @@ The energy is biased because:
</div>
<div id="outline-container-org1229a69" class="outline-4">
<h4 id="org1229a69"><span class="section-number-4">2.3.1</span> Exercise</h4>
<div id="outline-container-org33238fa" class="outline-4">
<h4 id="org33238fa"><span class="section-number-4">2.3.1</span> Exercise</h4>
<div class="outline-text-4" id="text-2-3-1">
<div class="exercise">
<p>
@ -1054,8 +1071,8 @@ To compile the Fortran and run it:
</div>
</div>
<div id="outline-container-orgad1ee8c" class="outline-3">
<h3 id="orgad1ee8c"><span class="section-number-3">2.4</span> Variance of the local energy</h3>
<div id="outline-container-org24aa290" class="outline-3">
<h3 id="org24aa290"><span class="section-number-3">2.4</span> Variance of the local energy</h3>
<div class="outline-text-3" id="text-2-4">
<p>
The variance of the local energy is a functional of \(\Psi\)
@ -1082,8 +1099,8 @@ energy can be used as a measure of the quality of a wave function.
</p>
</div>
<div id="outline-container-org0afc4b1" class="outline-4">
<h4 id="org0afc4b1"><span class="section-number-4">2.4.1</span> Exercise (optional)</h4>
<div id="outline-container-orgbd8786f" class="outline-4">
<h4 id="orgbd8786f"><span class="section-number-4">2.4.1</span> Exercise (optional)</h4>
<div class="outline-text-4" id="text-2-4-1">
<div class="exercise">
<p>
@ -1094,8 +1111,8 @@ Prove that :
</div>
</div>
</div>
<div id="outline-container-org10f6a19" class="outline-4">
<h4 id="org10f6a19"><span class="section-number-4">2.4.2</span> Exercise</h4>
<div id="outline-container-orgdf837df" class="outline-4">
<h4 id="orgdf837df"><span class="section-number-4">2.4.2</span> Exercise</h4>
<div class="outline-text-4" id="text-2-4-2">
<div class="exercise">
<p>
@ -1174,8 +1191,8 @@ To compile and run:
</div>
</div>
<div id="outline-container-org0de2b52" class="outline-2">
<h2 id="org0de2b52"><span class="section-number-2">3</span> Variational Monte Carlo</h2>
<div id="outline-container-orgf24e4c1" class="outline-2">
<h2 id="orgf24e4c1"><span class="section-number-2">3</span> Variational Monte Carlo</h2>
<div class="outline-text-2" id="text-3">
<p>
Numerical integration with deterministic methods is very efficient
@ -1191,8 +1208,8 @@ interval.
</p>
</div>
<div id="outline-container-org4abd3b8" class="outline-3">
<h3 id="org4abd3b8"><span class="section-number-3">3.1</span> Computation of the statistical error</h3>
<div id="outline-container-org5a7850e" class="outline-3">
<h3 id="org5a7850e"><span class="section-number-3">3.1</span> Computation of the statistical error</h3>
<div class="outline-text-3" id="text-3-1">
<p>
To compute the statistical error, you need to perform \(M\)
@ -1232,8 +1249,8 @@ And the confidence interval is given by
</p>
</div>
<div id="outline-container-org5a1c95b" class="outline-4">
<h4 id="org5a1c95b"><span class="section-number-4">3.1.1</span> Exercise</h4>
<div id="outline-container-org86ece31" class="outline-4">
<h4 id="org86ece31"><span class="section-number-4">3.1.1</span> Exercise</h4>
<div class="outline-text-4" id="text-3-1-1">
<div class="exercise">
<p>
@ -1275,8 +1292,8 @@ input array.
</div>
</div>
<div id="outline-container-org8a01ad5" class="outline-3">
<h3 id="org8a01ad5"><span class="section-number-3">3.2</span> Uniform sampling in the box</h3>
<div id="outline-container-orgf3538bc" class="outline-3">
<h3 id="orgf3538bc"><span class="section-number-3">3.2</span> Uniform sampling in the box</h3>
<div class="outline-text-3" id="text-3-2">
<p>
We will now perform our first Monte Carlo calculation to compute the
@ -1337,8 +1354,8 @@ compute the statistical error.
</p>
</div>
<div id="outline-container-orgae0d75a" class="outline-4">
<h4 id="orgae0d75a"><span class="section-number-4">3.2.1</span> Exercise</h4>
<div id="outline-container-org5ea7de2" class="outline-4">
<h4 id="org5ea7de2"><span class="section-number-4">3.2.1</span> Exercise</h4>
<div class="outline-text-4" id="text-3-2-1">
<div class="exercise">
<p>
@ -1442,8 +1459,8 @@ well as the index of the current step.
</div>
</div>
<div id="outline-container-org23b5713" class="outline-3">
<h3 id="org23b5713"><span class="section-number-3">3.3</span> Metropolis sampling with \(\Psi^2\)</h3>
<div id="outline-container-orgc5a1ec2" class="outline-3">
<h3 id="orgc5a1ec2"><span class="section-number-3">3.3</span> Metropolis sampling with \(\Psi^2\)</h3>
<div class="outline-text-3" id="text-3-3">
<p>
We will now use the square of the wave function to sample random
@ -1562,8 +1579,8 @@ All samples should be kept, from both accepted <i>and</i> rejected moves.
</div>
<div id="outline-container-org06c5905" class="outline-4">
<h4 id="org06c5905"><span class="section-number-4">3.3.1</span> Optimal step size</h4>
<div id="outline-container-org5e266bf" class="outline-4">
<h4 id="org5e266bf"><span class="section-number-4">3.3.1</span> Optimal step size</h4>
<div class="outline-text-4" id="text-3-3-1">
<p>
If the box is infinitely small, the ratio will be very close
@ -1598,8 +1615,8 @@ the same variable later on to store a time step.
</div>
<div id="outline-container-org654010b" class="outline-4">
<h4 id="org654010b"><span class="section-number-4">3.3.2</span> Exercise</h4>
<div id="outline-container-org7580d40" class="outline-4">
<h4 id="org7580d40"><span class="section-number-4">3.3.2</span> Exercise</h4>
<div class="outline-text-4" id="text-3-3-2">
<div class="exercise">
<p>
@ -1710,8 +1727,8 @@ Can you observe a reduction in the statistical error?
</div>
</div>
<div id="outline-container-org215f308" class="outline-3">
<h3 id="org215f308"><span class="section-number-3">3.4</span> Generalized Metropolis algorithm</h3>
<div id="outline-container-orgcaabae7" class="outline-3">
<h3 id="orgcaabae7"><span class="section-number-3">3.4</span> Generalized Metropolis algorithm</h3>
<div class="outline-text-3" id="text-3-4">
<p>
One can use more efficient numerical schemes to move the electrons by choosing a smarter expression for the transition probability.
@ -1832,8 +1849,8 @@ The algorithm of the previous exercise is only slighlty modified as:
</ol>
</div>
<div id="outline-container-org31cb975" class="outline-4">
<h4 id="org31cb975"><span class="section-number-4">3.4.1</span> Gaussian random number generator</h4>
<div id="outline-container-org5445700" class="outline-4">
<h4 id="org5445700"><span class="section-number-4">3.4.1</span> Gaussian random number generator</h4>
<div class="outline-text-4" id="text-3-4-1">
<p>
To obtain Gaussian-distributed random numbers, you can apply the
@ -1897,8 +1914,8 @@ In Python, you can use the <a href="https://numpy.org/doc/stable/reference/rando
</div>
<div id="outline-container-org48820db" class="outline-4">
<h4 id="org48820db"><span class="section-number-4">3.4.2</span> Exercise 1</h4>
<div id="outline-container-org1b2244b" class="outline-4">
<h4 id="org1b2244b"><span class="section-number-4">3.4.2</span> Exercise 1</h4>
<div class="outline-text-4" id="text-3-4-2">
<div class="exercise">
<p>
@ -1941,8 +1958,8 @@ Write a function to compute the drift vector \(\frac{\nabla \Psi(\mathbf{r})}{\P
</div>
</div>
<div id="outline-container-org1ce793c" class="outline-4">
<h4 id="org1ce793c"><span class="section-number-4">3.4.3</span> Exercise 2</h4>
<div id="outline-container-org0368ce7" class="outline-4">
<h4 id="org0368ce7"><span class="section-number-4">3.4.3</span> Exercise 2</h4>
<div class="outline-text-4" id="text-3-4-3">
<div class="exercise">
<p>
@ -2041,10 +2058,542 @@ Modify the previous program to introduce the drift-diffusion scheme.
</div>
</div>
<div id="outline-container-orgc4c9346" class="outline-2">
<h2 id="orgc4c9346"><span class="section-number-2">4</span> Project</h2>
<div id="outline-container-org82635bc" class="outline-2">
<h2 id="org82635bc"><span class="section-number-2">4</span> Diffusion Monte Carlo</h2>
<div class="outline-text-2" id="text-4">
<p>
As we have seen, Variational Monte Carlo is a powerful method to
compute integrals in large dimensions. It is often used in cases
where the expression of the wave function is such that the integrals
can't be evaluated (multi-centered Slater-type orbitals, correlation
factors, etc).
</p>
<p>
Diffusion Monte Carlo is different. It goes beyond the computation
of the integrals associated with an input wave function, and aims at
finding a near-exact numerical solution to the Schrödinger equation.
</p>
</div>
<div id="outline-container-org1ddc667" class="outline-3">
<h3 id="org1ddc667"><span class="section-number-3">4.1</span> Schrödinger equation in imaginary time</h3>
<div class="outline-text-3" id="text-4-1">
<p>
Consider the time-dependent Schrödinger equation:
</p>
<p>
\[
i\frac{\partial \Psi(\mathbf{r},t)}{\partial t} = (\hat{H} -E_{\rm ref}) \Psi(\mathbf{r},t)\,.
\]
</p>
<p>
where we introduced a shift in the energy, \(E_{\rm ref}\), for reasons which will become apparent below.
</p>
<p>
We can expand a given starting wave function, \(\Psi(\mathbf{r},0)\), in the basis of the eigenstates
of the time-independent Hamiltonian, \(\Phi_k\), with energies \(E_k\):
</p>
<p>
\[
\Psi(\mathbf{r},0) = \sum_k a_k\, \Phi_k(\mathbf{r}).
\]
</p>
<p>
The solution of the Schrödinger equation at time \(t\) is
</p>
<p>
\[
\Psi(\mathbf{r},t) = \sum_k a_k \exp \left( -i\, (E_k-E_{\rm ref})\, t \right) \Phi_k(\mathbf{r}).
\]
</p>
<p>
Now, if we replace the time variable \(t\) by an imaginary time variable
\(\tau=i\,t\), we obtain
</p>
<p>
\[
-\frac{\partial \psi(\mathbf{r}, \tau)}{\partial \tau} = (\hat{H} -E_{\rm ref}) \psi(\mathbf{r}, \tau)
\]
</p>
<p>
where \(\psi(\mathbf{r},\tau) = \Psi(\mathbf{r},-i\,\tau)\)
and
</p>
\begin{eqnarray*}
\psi(\mathbf{r},\tau) &=& \sum_k a_k \exp( -(E_k-E_{\rm ref})\, \tau) \Phi_k(\mathbf{r})\\
&=& \exp(-(E_0-E_{\rm ref})\, \tau)\sum_k a_k \exp( -(E_k-E_0)\, \tau) \Phi_k(\mathbf{r})\,.
\end{eqnarray*}
<p>
For large positive values of \(\tau\), \(\psi\) is dominated by the
\(k=0\) term, namely, the lowest eigenstate. If we adjust \(E_{\rm ref}\) to the running estimate of \(E_0\),
we can expect that simulating the differetial equation in
imaginary time will converge to the exact ground state of the
system.
</p>
</div>
</div>
<div id="outline-container-org7695533" class="outline-3">
<h3 id="org7695533"><span class="section-number-3">4.2</span> Relation to diffusion</h3>
<div class="outline-text-3" id="text-4-2">
<p>
The <a href="https://en.wikipedia.org/wiki/Diffusion_equation">diffusion equation</a> of particles is given by
</p>
<p>
\[
\frac{\partial \psi(\mathbf{r},t)}{\partial t} = D\, \Delta \psi(\mathbf{r},t)
\]
</p>
<p>
where \(D\) is the diffusion coefficient. When the imaginary-time
Schrödinger equation is written in terms of the kinetic energy and
potential,
</p>
<p>
\[
\frac{\partial \psi(\mathbf{r}, \tau)}{\partial \tau} =
\left(\frac{1}{2}\Delta - [V(\mathbf{r}) -E_{\rm ref}]\right) \psi(\mathbf{r}, \tau)\,,
\]
</p>
<p>
it can be identified as the combination of:
</p>
<ul class="org-ul">
<li>a diffusion equation (Laplacian)</li>
<li>an equation whose solution is an exponential (potential)</li>
</ul>
<p>
The diffusion equation can be simulated by a Brownian motion:
</p>
<p>
\[ \mathbf{r}_{n+1} = \mathbf{r}_{n} + \sqrt{\delta t}\, \chi \]
</p>
<p>
where \(\chi\) is a Gaussian random variable, and the potential term
can be simulated by creating or destroying particles over time (a
so-called branching process) or by simply considering it as a
cumulative multiplicative weight along the diffusion trajectory
(pure Diffusion Monte Carlo):
</p>
<p>
\[
\prod_i \exp \left( - (V(\mathbf{r}_i) - E_{\text{ref}}) \delta t \right).
\]
</p>
<p>
We note that the ground-state wave function of a Fermionic system is
antisymmetric and changes sign. Therefore, its interpretation as a probability
distribution is somewhat problematic. In fact, mathematically, since
the Bosonic ground state is lower in energy than the Fermionic one, for
large \(\tau\), the system will evolve towards the Bosonic solution.
</p>
<p>
For the systems you will study, this is not an issue:
</p>
<ul class="org-ul">
<li>Hydrogen atom: You only have one electron!</li>
<li>Two-electron system (\(H_2\) or He): The ground-wave function is
antisymmetric in the spin variables but symmetric in the space ones.</li>
</ul>
<p>
Therefore, in both cases, you are dealing with a "Bosonic" ground state.
</p>
</div>
</div>
<div id="outline-container-orgff86247" class="outline-3">
<h3 id="orgff86247"><span class="section-number-3">4.3</span> Importance sampling</h3>
<div class="outline-text-3" id="text-4-3">
<p>
In a molecular system, the potential is far from being constant
and, in fact, diverges at the inter-particle coalescence points. Hence,
it results in very large fluctuations of the erm weight associated with
the potental, making the calculations impossible in practice.
Fortunately, if we multiply the Schrödinger equation by a chosen
<i>trial wave function</i> \(\Psi_T(\mathbf{r})\) (Hartree-Fock, Kohn-Sham
determinant, CI wave function, <i>etc</i>), one obtains
</p>
<p>
\[
-\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) =
\left[ -\frac{1}{2} \Delta \psi(\mathbf{r},\tau) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \right] \Psi_T(\mathbf{r})
\]
</p>
<p>
Defining \(\Pi(\mathbf{r},\tau) = \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})\), (see appendix for details)
</p>
<p>
\[
-\frac{\partial \Pi(\mathbf{r},\tau)}{\partial \tau}
= -\frac{1}{2} \Delta \Pi(\mathbf{r},\tau) +
\nabla \left[ \Pi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})}
\right] + (E_L(\mathbf{r})-E_{\rm ref})\Pi(\mathbf{r},\tau)
\]
</p>
<p>
The new "kinetic energy" can be simulated by the drift-diffusion
scheme presented in the previous section (VMC).
The new "potential" is the local energy, which has smaller fluctuations
when \(\Psi_T\) gets closer to the exact wave function.
This term can be simulated by
\[
\prod_i \exp \left( - (E_L(\mathbf{r}_i) - E_{\text{ref}}) \delta t \right).
\]
where \(E_{\rm ref}\) is the constant we had introduced above, which is adjusted to
an estimate of the average energy to keep the weights close to one.
</p>
<p>
This equation generates the <i>N</i>-electron density \(\Pi\), which is the
product of the ground state solution with the trial wave
function. You may then ask: how can we compute the total energy of
the system?
</p>
<p>
To this aim, we use the <i>mixed estimator</i> of the energy:
</p>
\begin{eqnarray*}
E(\tau) &=& \frac{\langle \psi(\tau) | \hat{H} | \Psi_T \rangle}{\langle \psi(\tau) | \Psi_T \rangle}\\
&=& \frac{\int \psi(\mathbf{r},\tau) \hat{H} \Psi_T(\mathbf{r}) d\mathbf{r}}
{\int \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) d\mathbf{r}} \\
&=& \frac{\int \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) E_L(\mathbf{r}) d\mathbf{r}}
{\int \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) d\mathbf{r}} \,.
\end{eqnarray*}
<p>
For large \(\tau\), we have that
</p>
<p>
\[
\Pi(\mathbf{r},\tau) =\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \rightarrow \Phi_0(\mathbf{r}) \Psi_T(\mathbf{r})\,,
\]
</p>
<p>
and, using that \(\hat{H}\) is Hermitian and that \(\Phi_0\) is an
eigenstate of the Hamiltonian, we obtain for large \(\tau\)
</p>
<p>
\[
E(\tau) = \frac{\langle \psi_\tau | \hat{H} | \Psi_T \rangle}
{\langle \psi_\tau | \Psi_T \rangle}
= \frac{\langle \Psi_T | \hat{H} | \psi_\tau \rangle}
{\langle \Psi_T | \psi_\tau \rangle}
\rightarrow E_0 \frac{\langle \Psi_T | \Phi_0 \rangle}
{\langle \Psi_T | \Phi_0 \rangle}
= E_0
\]
</p>
<p>
Therefore, we can compute the energy within DMC by generating the
density \(\Pi\) with random walks, and simply averaging the local
energies computed with the trial wave function.
</p>
</div>
<div id="outline-container-orgc628347" class="outline-4">
<h4 id="orgc628347"><span class="section-number-4">4.3.1</span> Appendix : Details of the Derivation</h4>
<div class="outline-text-4" id="text-4-3-1">
<p>
\[
-\frac{\partial \psi(\mathbf{r},\tau)}{\partial \tau} \Psi_T(\mathbf{r}) =
\left[ -\frac{1}{2} \Delta \psi(\mathbf{r},\tau) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \right] \Psi_T(\mathbf{r})
\]
</p>
<p>
\[
-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
= -\frac{1}{2} \Big( \Delta \big[
\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] -
\psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) - 2
\nabla \psi(\mathbf{r},\tau) \nabla \Psi_T(\mathbf{r}) \Big) + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
\]
</p>
<p>
\[
-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] +
\frac{1}{2} \psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) +
\Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + V(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
\]
</p>
<p>
\[
-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] +
\psi(\mathbf{r},\tau) \Delta \Psi_T(\mathbf{r}) +
\Psi_T(\mathbf{r})\nabla \psi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})} + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
\]
\[
-\frac{\partial \big[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big]}{\partial \tau}
= -\frac{1}{2} \Delta \big[\psi(\mathbf{r},\tau) \Psi_T(\mathbf{r}) \big] +
\nabla \left[ \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
\frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})}
\right] + E_L(\mathbf{r}) \psi(\mathbf{r},\tau) \Psi_T(\mathbf{r})
\]
</p>
<p>
Defining \(\Pi(\mathbf{r},t) = \psi(\mathbf{r},\tau)
\Psi_T(\mathbf{r})\),
</p>
<p>
\[
-\frac{\partial \Pi(\mathbf{r},\tau)}{\partial \tau}
= -\frac{1}{2} \Delta \Pi(\mathbf{r},\tau) +
\nabla \left[ \Pi(\mathbf{r},\tau) \frac{\nabla \Psi_T(\mathbf{r})}{\Psi_T(\mathbf{r})}
\right] + E_L(\mathbf{r}) \Pi(\mathbf{r},\tau)
\]
</p>
</div>
</div>
</div>
<div id="outline-container-org429350d" class="outline-3">
<h3 id="org429350d"><span class="section-number-3">4.4</span> Pure Diffusion Monte Carlo</h3>
<div class="outline-text-3" id="text-4-4">
<p>
Instead of having a variable number of particles to simulate the
branching process as in the <i>Diffusion Monte Carlo</i> (DMC) algorithm, we
use variant called <i>pure Diffusion Monte Carlo</i> (PDMC) where
the potential term is considered as a cumulative product of weights:
</p>
\begin{eqnarray*}
W(\mathbf{r}_n, \tau) = \prod_{i=1}^{n} \exp \left( -\delta t\,
(E_L(\mathbf{r}_i) - E_{\text{ref}}) \right) =
\prod_{i=1}^{n} w(\mathbf{r}_i)
\end{eqnarray*}
<p>
where \(\mathbf{r}_i\) are the coordinates along the trajectory and
we introduced a <i>time-step variable</i> \(\delta t\) to discretize the
integral.
</p>
<p>
The PDMC algorithm is less stable than the DMC algorithm: it
requires to have a value of \(E_\text{ref}\) which is close to the
fixed-node energy, and a good trial wave function. Moreover, we
can't let \(\tau\) become too large because the weight whether
explode or vanish: we need to have a fixed value of \(\tau\)
(projection time).
The big advantage of PDMC is that it is rather simple to implement
starting from a VMC code:
</p>
<ol class="org-ol">
<li>Start with \(W(\mathbf{r}_0)=1, \tau_0 = 0\)</li>
<li>Evaluate the local energy at \(\mathbf{r}_{n}\)</li>
<li>Compute the contribution to the weight \(w(\mathbf{r}_n) =
\exp(-\delta t(E_L(\mathbf{r}_n)-E_\text{ref}))\)</li>
<li>Update \(W(\mathbf{r}_{n}) = W(\mathbf{r}_{n-1}) \times w(\mathbf{r}_n)\)</li>
<li>Accumulate the weighted energy \(W(\mathbf{r}_n) \times
E_L(\mathbf{r}_n)\),
and the weight \(W(\mathbf{r}_n)\) for the normalization</li>
<li>Update \(\tau_n = \tau_{n-1} + \delta t\)</li>
<li>If \(\tau_{n} > \tau_\text{max}\), the long projection time has
been reached and we can start an new trajectory from the current
position: reset \(W(r_n) = 1\) and \(\tau_n
= 0\)</li>
<li>Compute a new position \(\mathbf{r'} = \mathbf{r}_n +
\delta t\, \frac{\nabla \Psi(\mathbf{r})}{\Psi(\mathbf{r})} + \chi\)</li>
<li>Evaluate \(\Psi(\mathbf{r}')\) and \(\frac{\nabla \Psi(\mathbf{r'})}{\Psi(\mathbf{r'})}\) at the new position</li>
<li>Compute the ratio \(A = \frac{T(\mathbf{r}' \rightarrow \mathbf{r}_{n}) P(\mathbf{r}')}{T(\mathbf{r}_{n} \rightarrow \mathbf{r}') P(\mathbf{r}_{n})}\)</li>
</ol>
<ol class="org-ol">
<li>Draw a uniform random number \(v \in [0,1]\)</li>
<li>if \(v \le A\), accept the move : set \(\mathbf{r}_{n+1} = \mathbf{r'}\)</li>
<li>else, reject the move : set \(\mathbf{r}_{n+1} = \mathbf{r}_n\)</li>
</ol>
<p>
Some comments are needed:
</p>
<ul class="org-ul">
<li><p>
You estimate the energy as
</p>
\begin{eqnarray*}
E = \frac{\sum_{k=1}^{N_{\rm MC}} E_L(\mathbf{r}_k) W(\mathbf{r}_k, k\delta t)}{\sum_{k=1}^{N_{\rm MC}} W(\mathbf{r}_k, k\delta t)}
\end{eqnarray*}</li>
<li><p>
The result will be affected by a time-step error
(the finite size of \(\delta t\)) due to the discretization of the
integral, and one has in principle to extrapolate to the limit
\(\delta t \rightarrow 0\). This amounts to fitting the energy
computed for multiple values of \(\delta t\).
</p>
<p>
Here, you will be using a small enough time-step and you should not worry about the extrapolation.
</p></li>
<li>The accept/reject step (steps 9-12 in the algorithm) is in principle not needed for the correctness of
the DMC algorithm. However, its use reduces significantly the time-step error.</li>
</ul>
</div>
</div>
<div id="outline-container-org7b3d04e" class="outline-3">
<h3 id="org7b3d04e"><span class="section-number-3">4.5</span> Hydrogen atom</h3>
<div class="outline-text-3" id="text-4-5">
</div>
<div id="outline-container-org1a74409" class="outline-4">
<h4 id="org1a74409"><span class="section-number-4">4.5.1</span> Exercise</h4>
<div class="outline-text-4" id="text-4-5-1">
<div class="exercise">
<p>
Modify the Metropolis VMC program into a PDMC program.
In the limit \(\delta t \rightarrow 0\), you should recover the exact
energy of H for any value of \(a\), as long as the simulation is stable.
We choose here a time step of 0.05 a.u. and a fixed projection
time \(\tau\) =100 a.u.
</p>
</div>
<p>
<b>Python</b>
</p>
<div class="org-src-container">
<pre class="src src-python"><span style="color: #a020f0;">from</span> hydrogen <span style="color: #a020f0;">import</span> *
<span style="color: #a020f0;">from</span> qmc_stats <span style="color: #a020f0;">import</span> *
<span style="color: #a020f0;">def</span> <span style="color: #0000ff;">MonteCarlo</span>(a, nmax, dt, Eref):
# <span style="color: #b22222;">TODO</span>
# <span style="color: #b22222;">Run simulation</span>
<span style="color: #a0522d;">a</span> = 1.2
<span style="color: #a0522d;">nmax</span> = 100000
<span style="color: #a0522d;">dt</span> = 0.05
<span style="color: #a0522d;">tau</span> = 100.
<span style="color: #a0522d;">E_ref</span> = -0.5
<span style="color: #a0522d;">X0</span> = [ MonteCarlo(a, nmax, dt, E_ref) <span style="color: #a020f0;">for</span> i <span style="color: #a020f0;">in</span> <span style="color: #483d8b;">range</span>(30)]
# <span style="color: #b22222;">Energy</span>
<span style="color: #a0522d;">X</span> = [ x <span style="color: #a020f0;">for</span> (x, _) <span style="color: #a020f0;">in</span> X0 ]
<span style="color: #a0522d;">E</span>, <span style="color: #a0522d;">deltaE</span> = ave_error(X)
<span style="color: #a020f0;">print</span>(f<span style="color: #8b2252;">"E = {E} +/- {deltaE}"</span>)
# <span style="color: #b22222;">Acceptance rate</span>
<span style="color: #a0522d;">X</span> = [ x <span style="color: #a020f0;">for</span> (_, x) <span style="color: #a020f0;">in</span> X0 ]
<span style="color: #a0522d;">A</span>, <span style="color: #a0522d;">deltaA</span> = ave_error(X)
<span style="color: #a020f0;">print</span>(f<span style="color: #8b2252;">"A = {A} +/- {deltaA}"</span>)
</pre>
</div>
<p>
<b>Fortran</b>
</p>
<div class="org-src-container">
<pre class="src src-f90"><span style="color: #a020f0;">subroutine</span> <span style="color: #0000ff;">pdmc</span>(a, dt, nmax, energy, accep, tau, E_ref)
<span style="color: #a020f0;">implicit</span> <span style="color: #228b22;">none</span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">intent</span>(in) ::<span style="color: #a0522d;"> a, dt, tau</span>
<span style="color: #228b22;">integer</span>*8 , <span style="color: #a020f0;">intent</span>(in) ::<span style="color: #a0522d;"> nmax </span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">intent</span>(out) ::<span style="color: #a0522d;"> energy, accep</span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">intent</span>(in) ::<span style="color: #a0522d;"> E_ref</span>
<span style="color: #228b22;">integer</span>*8 ::<span style="color: #a0522d;"> istep</span>
<span style="color: #228b22;">integer</span>*8 ::<span style="color: #a0522d;"> n_accep</span>
<span style="color: #228b22;">double precision</span> ::<span style="color: #a0522d;"> sq_dt, chi(3)</span>
<span style="color: #228b22;">double precision</span> ::<span style="color: #a0522d;"> psi_old, psi_new</span>
<span style="color: #228b22;">double precision</span> ::<span style="color: #a0522d;"> r_old(3), r_new(3)</span>
<span style="color: #228b22;">double precision</span> ::<span style="color: #a0522d;"> d_old(3), d_new(3)</span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">external</span> ::<span style="color: #a0522d;"> e_loc, psi</span>
! <span style="color: #b22222;">TODO</span>
<span style="color: #a020f0;">end subroutine</span> <span style="color: #0000ff;">pdmc</span>
<span style="color: #a020f0;">program</span> <span style="color: #0000ff;">qmc</span>
<span style="color: #a020f0;">implicit</span> <span style="color: #228b22;">none</span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">parameter</span> ::<span style="color: #a0522d;"> a = 1.2d0</span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">parameter</span> ::<span style="color: #a0522d;"> dt = 0.05d0</span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">parameter</span> ::<span style="color: #a0522d;"> E_ref = -0.5d0</span>
<span style="color: #228b22;">double precision</span>, <span style="color: #a020f0;">parameter</span> ::<span style="color: #a0522d;"> tau = 100.d0</span>
<span style="color: #228b22;">integer</span>*8 , <span style="color: #a020f0;">parameter</span> ::<span style="color: #a0522d;"> nmax = 100000</span>
<span style="color: #228b22;">integer</span> , <span style="color: #a020f0;">parameter</span> ::<span style="color: #a0522d;"> nruns = 30</span>
<span style="color: #228b22;">integer</span> ::<span style="color: #a0522d;"> irun</span>
<span style="color: #228b22;">double precision</span> ::<span style="color: #a0522d;"> X(nruns), accep(nruns)</span>
<span style="color: #228b22;">double precision</span> ::<span style="color: #a0522d;"> ave, err</span>
<span style="color: #a020f0;">do</span> irun=1,nruns
<span style="color: #a020f0;">call</span> <span style="color: #0000ff;">pdmc</span>(a, dt, nmax, X(irun), accep(irun), tau, E_ref)
<span style="color: #a020f0;">enddo</span>
<span style="color: #a020f0;">call</span> <span style="color: #0000ff;">ave_error</span>(X,nruns,ave,err)
<span style="color: #a020f0;">print</span> *, <span style="color: #8b2252;">'E = '</span>, ave, <span style="color: #8b2252;">'+/-'</span>, err
<span style="color: #a020f0;">call</span> <span style="color: #0000ff;">ave_error</span>(accep,nruns,ave,err)
<span style="color: #a020f0;">print</span> *, <span style="color: #8b2252;">'A = '</span>, ave, <span style="color: #8b2252;">'+/-'</span>, err
<span style="color: #a020f0;">end program</span> <span style="color: #0000ff;">qmc</span>
</pre>
</div>
<div class="org-src-container">
<pre class="src src-sh">gfortran hydrogen.f90 qmc_stats.f90 pdmc.f90 -o pdmc
./pdmc
</pre>
</div>
</div>
</div>
</div>
</div>
<div id="outline-container-org60565bd" class="outline-2">
<h2 id="org60565bd"><span class="section-number-2">5</span> Project</h2>
<div class="outline-text-2" id="text-5">
<p>
Change your PDMC code for one of the following:
</p>
<ul class="org-ul">
@ -2059,9 +2608,9 @@ And compute the ground state energy.
</div>
<div id="outline-container-org79362c8" class="outline-2">
<h2 id="org79362c8"><span class="section-number-2">5</span> Acknowledgments</h2>
<div class="outline-text-2" id="text-5">
<div id="outline-container-orga85900b" class="outline-2">
<h2 id="orga85900b"><span class="section-number-2">6</span> Acknowledgments</h2>
<div class="outline-text-2" id="text-6">
<div class="figure">
<p><img src="https://trex-coe.eu/sites/default/files/inline-images/euflag.jpg" alt="euflag.jpg" />
@ -2080,7 +2629,7 @@ Union is not responsible for any use that might be made of such content.
</div>
<div id="postamble" class="status">
<p class="author">Author: Anthony Scemama, Claudia Filippi</p>
<p class="date">Created: 2021-02-04 Thu 07:47</p>
<p class="date">Created: 2021-02-04 Thu 10:25</p>
<p class="validation"><a href="http://validator.w3.org/check?uri=referer">Validate</a></p>
</div>
</body>