<jupyter_start><jupyter_text>From Variational Monte Carlo to Boltzmann Machines and Machine Learning  **<NAME>  Email <EMAIL>  Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University, East Lansing, 48824 MI, USA**Date: **Notebook 1:  Variational Monte Carlo** Introduction Structure and AimsThese notebooks serve the aim of linking traditional variational MonteCarlo VMC calculations methods with recent progress on solvingmany-particle problems using Machine Learning algorithms.Furthermore, when linking with Machine Learning algorithms, in particularso-called Boltzmann Machines, there are interesting connections betweenthese algorithms and so-called [Shadow Wave functions (SWFs)](https://journals.aps.org/pre/abstract/10.1103/PhysRevE.90.053304) (and references therein). The implications of the latter have been explored in various Monte Carlo calculations. In total there are three notebooks:1. the one you are reading now on Variational Monte Carlo methods, 2. notebook 2 on Machine Learning and quantum mechanical problems and in particular on Boltzmann Machines, 3. and finally notebook 3 on the link between Boltzmann machines and SWFs.  This notebookIn this notebook the aim is to give you an introduction as well as anunderstanding of the basic elements that are needed in order todevelop a professional variational Monte Carlo code. We will focus ona simple system of two particles in an oscillator trap (oralternatively two fermions moving in a Coulombic potential) which caninteract via repulsive or attrative force.The advantage of these systems is that for two particles (boson orfermions) we have analytical solutions for the eigenpairs for thenon-interacting case. Furthermore, for a two- or three-dimensionalsystem of two electrons moving in a harmonic oscillator trap, we have[analytical solutions for the interacting case as well](https://iopscience.iop.org/article/10.1088/0305-4470/27/3/040/meta).  Having analytical eigenpairs is an invaluable feature that allow us to assess the physical relevance of the trial wave functions, bethese either from a standard VMC procedure, from Boltzmann Machines orfrom Shadow Wave functions.In this notebook we start with the basics of a VMC calculation andintroduce concepts like Markow Chain Monte Carlo methods and theMetropolis algorithm, importance sampling and Metropolis-Hastingsalgorithm, resampling methods to obtain better estimates of thestatistical errors and minimization of the expectation values of theenergy and the variance. The latter is done in order to obtain thebest possible variational parameters. Furthermore it will define theso-called **cost** function, a commonly encountered quantity in MachineLearning algorithms. Minimizing the latter is the one which leads tothe determination of the optimal parameters in basically all Machine Learning algorithms.For our purposes, it will serve as the first link between VMC methods and Machine Learning methods.Topics like Markov Chain Monte Carlo and various resampling techniquesare also central to Machine Learning methods. Presenting them in thecontext of VMC approaches leads hopefully to an easier starting pointfor the understanding of these methods.Finally, the reader may ask what do we actually want to achieve withcomplicating life with Machine Learning methods when we can easilystudy interacting systems with standard Monte Carlo approaches.  Ourhope is that by adding additional degrees of freedom via MachineLearning algorithms, we can let the algorithms we employ learn theparameters of the model via a given optimization algorithm. Instandard Monte Carlo calculations the practitioners end fine tuningthe trial wave function using all possible insights about the systemunderstudy. This may not always lead to the best possible ansatz andcan in the long run be rather time-consuming. In fields like nuclearmany-body physics with complicated interaction terms, guessing ananalytical form for the trial wave fuction can be difficult. Lettingthe machine learn the form of the trial function or find the optimalparameters may lead to insights about the problem which cannot beobtained by selecting various trial wave functions.The emerging and rapidly expanding fields of Machine Learning and Quantum Computing hold also great promise in tackling the dimensionality problems (the so-called dimensionality curse in many-body problems) we encounter when studying complicated many-body problems. The approach to Machine Learning we will focus on  is inspired by the idea of representing the wave function witha restricted Boltzmann machine (RBM), presented recently by [<NAME> and <NAME>, Science **355**, Issue 6325, pp. 602-606 (2017)](http://science.sciencemag.org/content/355/6325/602). Theynamed such a wave function/network a *neural network quantum state* (NQS). In their article they apply it to the quantum mechanicalspin lattice systems of the Ising model and Heisenberg model, withencouraging results.Machine learning (ML) is an extremely rich field, in spite of its young age. Theincreases we have seen during the last three decades in computationalcapabilities have been followed by developments of methods andtechniques for analyzing and handling large date sets, relying heavilyon statistics, computer science and mathematics.  The field is rathernew and developing rapidly. Machine learning is the science of giving computers the ability tolearn without being explicitly programmed.  The idea is that thereexist generic algorithms which can be used to find patterns in a broadclass of data sets without having to write code specifically for eachproblem. The algorithm will build its own logic based on the data.Machine learning is a subfield of computer science, and is closelyrelated to computational statistics.  It evolved from the study ofpattern recognition in artificial intelligence (AI) research, and hasmade contributions to AI tasks like computer vision, natural languageprocessing and speech recognition. It has also, especially in lateryears, found applications in a wide variety of other areas, includingbioinformatics, economy, physics, finance and marketing.An excellent reference we will come to back to is [Mehta *et al.*, arXiv:1803.08823](https://arxiv.org/abs/1803.08823).Our focus will first be on the basics of VMC calculations.  Basic Quantum Monte CarloWe start with the variational principle.Given a hamiltonian $H$ and a trial wave function $\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})$, the variational principle states that the expectation value of $\langle H \rangle$, defined through $$E[H]= \langle H \rangle =   \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})}        {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})},$$ is an upper bound to the ground state energy $E_0$ of the hamiltonian $H$, that is $$E_0 \le \langle H \rangle .$$ In general, the integrals involved in the calculation of various  expectation values  are multi-dimensional ones. Traditional integration methods such as the Gauss-Legendre will not be adequate for say the  computation of the energy of a many-body system.Here we have defined the vector $\boldsymbol{R} =[\boldsymbol{r}_1,\boldsymbol{r}_2,\dots,\boldsymbol{r}_n}]$  as an array that contains the positions of all particles $n$ while the vector $\boldsymbol{\alpha} = [\alpha_1,\alpha_2,\dots,\alpha_m]$ contains the variational parameters of the model, $m$ in total. The trial wave function can be expanded in the eigenstates of the hamiltonian since they form a complete set, viz., $$\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})=\sum_i a_i\Psi_i(\boldsymbol{R}),$$ and assuming the set of eigenfunctions to be normalized one obtains $$\frac{\sum_{nm}a^*_ma_n \int d\boldsymbol{R}\Psi^{\ast}_m(\boldsymbol{R})H(\boldsymbol{R})\Psi_n(\boldsymbol{R})}        {\sum_{nm}a^*_ma_n \int d\boldsymbol{R}\Psi^{\ast}_m(\boldsymbol{R})\Psi_n(\boldsymbol{R})} =\frac{\sum_{n}a^2_n E_n}        {\sum_{n}a^2_n} \ge E_0,$$ where we used that $H(\boldsymbol{R})\Psi_n(\boldsymbol{R})=E_n\Psi_n(\boldsymbol{R})$.In general, the integrals involved in the calculation of various  expectationvalues  are multi-dimensional ones. The variational principle yields the lowest state of a given symmetry.In most cases, a wave function has only small values in large parts of configuration space, and a straightforward procedure which useshomogenously distributed random points in configuration space will most likely lead to poor results. This may suggest that some kindof importance sampling combined with e.g., the Metropolis algorithm may be  a more efficient way of obtaining the ground state energy.The hope is then that those regions of configurations space wherethe wave function assumes appreciable values are sampled more efficiently. The tedious part in a VMC calculation is the search for the variationalminimum. A good knowledge of the system is required in order to carry outreasonable VMC calculations. This is not always the case, and often VMC calculations serve rather as the startingpoint for so-called diffusion Monte Carlo calculations (DMC). DMC is a way ofsolving exactly the many-body Schroedinger equation by means of a stochastic procedure. A good guess on the binding energyand its wave function is however necessary. A carefully performed VMC calculation can aid in this context. The basic procedure of a Variational Monte Carlo calculations consists thus of 1. Construct first a trial wave function $\psi_T(\boldsymbol{R};\boldsymbol{\alpha})$,  for a many-body system consisting of $n$ particles located at positions  $\boldsymbol{R}=(\boldsymbol{R}_1,\dots ,\boldsymbol{R}_n)$. The trial wave function depends on $\alpha$ variational parameters $\boldsymbol{\alpha}=(\alpha_1,\dots ,\alpha_M)$.2. Then we evaluate the expectation value of the hamiltonian $H$ $$\overline{E}[\boldsymbol{alpha}]=\frac{\int d\boldsymbol{R}\Psi^{\ast}_{T}(\boldsymbol{R},\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_{T}(\boldsymbol{R},\boldsymbol{\alpha})}        {\int d\boldsymbol{R}\Psi^{\ast}_{T}(\boldsymbol{R},\boldsymbol{\alpha})\Psi_{T}(\boldsymbol{R},\boldsymbol{\alpha})}.$$ 1. Thereafter we vary $\boldsymbol{\alpha}$ according to some minimization algorithm and return eventually to the first step if we are not satisfied with the results.Here we have used the notation $\overline{E}$ to label the expectation value of the energy.  Linking with standard statistical expressions for expectation valuesIn order to bring in the Monte Carlo machinery, we define first a likelihood distribution, or probability density distribution (PDF). Using our ansatz for the trial wave function $\psi_T(\boldsymbol{R};\boldsymbol{\alpha})$. $$P(\boldsymbol{R})= \frac{\left|\psi_T(\boldsymbol{R};\boldsymbol{\alpha})\right|^2}{\int \left|\psi_T(\boldsymbol{R};\boldsymbol{\alpha})\right|^2d\boldsymbol{R}}.$$ This is our new probability distribution function  (PDF).The approximation to the expectation value of the Hamiltonian is now $$\overline{E}[\boldsymbol{\alpha}] =    \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})}        {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})}.$$ Define a new quantity $$E_L(\boldsymbol{R};\boldsymbol{\alpha})=\frac{1}{\psi_T(\boldsymbol{R};\boldsymbol{\alpha})}H\psi_T(\boldsymbol{R};\boldsymbol{\alpha}),\label{eq:locale1} \tag{1}$$ called the local energy, which, together with our trial PDF yields $$\overline{E}[\boldsymbol{\alpha}]=\int P(\boldsymbol{R})E_L(\boldsymbol{R};\boldsymbol{\alpha}) d\boldsymbol{R}\approx \frac{1}{N}\sum_{i=1}^NE_L(\boldsymbol{R_i};\boldsymbol{\alpha})\label{eq:vmc1} \tag{2}$$ with $N$ being the number of Monte Carlo samples.The Algorithm for performing a variational Monte Carlo calculations runs thus as this   * Initialisation: Fix the number of Monte Carlo steps. Choose an initial $\boldsymbol{R}$ and variational parameters $\alpha$ and calculate $\left|\psi_T^{\alpha}(\boldsymbol{R})\right|^2$.    * Initialise the energy and the variance and start the Monte Carlo calculation.      * Calculate  a trial position  $\boldsymbol{R}_p=\boldsymbol{R}+r*step$ where $r$ is a random variable $r \in [0,1]$.      * Metropolis algorithm to accept or reject this move  $w = P(\boldsymbol{R}_p)/P(\boldsymbol{R})$.      * If the step is accepted, then we set $\boldsymbol{R}=\boldsymbol{R}_p$.       * Update averages   * Finish and compute final averages.Observe that the jumping in space is governed by the variable *step*. This is called brute-force sampling.Need importance sampling to get more relevant sampling, see lectures below. Simple example, the hydrogen atomThe radial Schroedinger equation for the hydrogen atom can bewritten as (when we have gotten rid of the first derivative term in the kinetic energy and used $rR(r)=u(r)$) $$-\frac{\hbar^2}{2m}\frac{\partial^2 u(r)}{\partial r^2}-\left(\frac{ke^2}{r}-\frac{\hbar^2l(l+1)}{2mr^2}\right)u(r)=Eu(r).$$ or with dimensionless variables $$-\frac{1}{2}\frac{\partial^2 u(\rho)}{\partial \rho^2}-\frac{u(\rho)}{\rho}+\frac{l(l+1)}{2\rho^2}u(\rho)-\lambda u(\rho)=0,\label{eq:hydrodimless1} \tag{3}$$ with the hamiltonian $$H=-\frac{1}{2}\frac{\partial^2 }{\partial \rho^2}-\frac{1}{\rho}+\frac{l(l+1)}{2\rho^2}.$$ Use variational parameter $\alpha$ in the trialwave function $$u_T^{\alpha}(\rho)=\alpha\rho e^{-\alpha\rho}. \label{eq:trialhydrogen} \tag{4}$$ Inserting this wave function into the expression for thelocal energy $E_L$ gives $$E_L(\rho)=-\frac{1}{\rho}-              \frac{\alpha}{2}\left(\alpha-\frac{2}{\rho}\right).$$ We note that at $\alpha=1$ we obtain the exactresult, and the variance is zero, as it should. The reason is that we then have the exact wave function, and the action of the hamiltionanon the wave function $$H\psi = \mathrm{constant}\times \psi,$$ yields just a constant. The integral which defines various expectation values involving moments of the hamiltonian becomes then $$\langle H^n \rangle =   \frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})H^n(\boldsymbol{R})\Psi_T(\boldsymbol{R})}        {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}=\mathrm{constant}\times\frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}        {\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}=\mathrm{constant}.$$ **This gives an important information: the exact wave function leads to zero variance!**Variation is then performed by minimizing both the energy and the variance. The Metropolis algorithmThe Metropolis algorithm , see [the original article](http://scitation.aip.org/content/aip/journal/jcp/21/6/10.1063/1.1699114) was invented by Metropolis et. aland is often simply called the Metropolis algorithm.It is a method to sample a normalized probabilitydistribution by a stochastic process. We define $\mathbf{P}_i^{(n)}$ tobe the probability for finding the system in the state $i$ at step $n$.The algorithm is then* Sample a possible new state $j$ with some probability $T_{i\rightarrow j}$.* Accept the new state $j$ with probability $A_{i \rightarrow j}$ and use it as the next sample. With probability $1-A_{i\rightarrow j}$ the move is rejected and the original state $i$ is used again as a sample.We wish to derive the required properties of $T$ and $A$ such that$\mathbf{P}_i^{(n\rightarrow \infty)} \rightarrow p_i$ so that startingfrom any distribution, the method converges to the correct distribution.Note that the description here is for a discrete probability distribution.Replacing probabilities $p_i$ with expressions like $p(x_i)dx_i$ willtake all of these over to the corresponding continuum expressions.The dynamical equation for $\mathbf{P}_i^{(n)}$ can be written directly fromthe description above. The probability of being in the state $i$ at step $n$is given by the probability of being in any state $j$ at the previous step,and making an accepted transition to $i$ added to the probability ofbeing in the state $i$, making a transition to any state $j$ andrejecting the move: $$\mathbf{P}^{(n)}_i = \sum_j \left [\mathbf{P}^{(n-1)}_jT_{j\rightarrow i} A_{j\rightarrow i} +\mathbf{P}^{(n-1)}_iT_{i\rightarrow j}\left ( 1- A_{i\rightarrow j} \right)\right ] \,.$$ Since the probability of making some transition must be 1,$\sum_j T_{i\rightarrow j} = 1$, and the above equation becomes $$\mathbf{P}^{(n)}_i = \mathbf{P}^{(n-1)}_i + \sum_j \left [\mathbf{P}^{(n-1)}_jT_{j\rightarrow i} A_{j\rightarrow i} -\mathbf{P}^{(n-1)}_iT_{i\rightarrow j}A_{i\rightarrow j}\right ] \,.$$ For large $n$ we require that $\mathbf{P}^{(n\rightarrow \infty)}_i = p_i$,the desired probability distribution. Taking this limit, gives thebalance requirement $$\sum_j \left [p_jT_{j\rightarrow i} A_{j\rightarrow i}-p_iT_{i\rightarrow j}A_{i\rightarrow j}\right ] = 0 \,.$$ The balance requirement is very weak. Typically the much stronger detailedbalance requirement is enforced, that is rather than the sum beingset to zero, we set each term separately to zero and use thisto determine the acceptance probabilities. Rearranging, the result is $$\frac{ A_{j\rightarrow i}}{A_{i\rightarrow j}}= \frac{p_iT_{i\rightarrow j}}{ p_jT_{j\rightarrow i}} \,.$$ The Metropolis choice is to maximize the $A$ values, that is $$A_{j \rightarrow i} = \min \left ( 1,\frac{p_iT_{i\rightarrow j}}{ p_jT_{j\rightarrow i}}\right ).$$ Other choices are possible, but they all correspond to multilplying$A_{i\rightarrow j}$ and $A_{j\rightarrow i}$ by the same constantsmaller than unity.\footnote{The penalty function method uses just sucha factor to compensate for $p_i$ that are evaluated stochasticallyand are therefore noisy.}Having chosen the acceptance probabilities, we have guaranteed thatif the  $\mathbf{P}_i^{(n)}$ has equilibrated, that is if it is equal to $p_i$,it will remain equilibrated. Next we need to find the circumstances forconvergence to equilibrium.The dynamical equation can be written as $$\mathbf{P}^{(n)}_i = \sum_j M_{ij}\mathbf{P}^{(n-1)}_j$$ with the matrix $M$ given by $$M_{ij} = \delta_{ij}\left [ 1 -\sum_k T_{i\rightarrow k} A_{i \rightarrow k}\right ] + T_{j\rightarrow i} A_{j\rightarrow i} \,.$$ Summing over $i$ shows that $\sum_i M_{ij} = 1$, and since$\sum_k T_{i\rightarrow k} = 1$, and $A_{i \rightarrow k} \leq 1$, theelements of the matrix satisfy $M_{ij} \geq 0$. The matrix $M$ is thereforea stochastic matrix.The Metropolis method is simply the power method for computing theright eigenvector of $M$ with the largest magnitude eigenvalue.By construction, the correct probability distribution is a right eigenvectorwith eigenvalue 1. Therefore, for the Metropolis method to convergeto this result, we must show that $M$ has only one eigenvalue with thismagnitude, and all other eigenvalues are smaller. The system: two electrons in a harmonic oscillator trap in two dimensionsThe Hamiltonian of the quantum dot is given by $$\hat{H} = \hat{H}_0 + \hat{V},$$ where $\hat{H}_0$ is the many-body HO Hamiltonian, and $\hat{V}$ is theinter-electron Coulomb interactions. In dimensionless units, $$\hat{V}= \sum_{i < j}^N \frac{1}{r_{ij}},$$ with $r_{ij}=\sqrt{\mathbf{r}_i^2 - \mathbf{r}_j^2}$.This leads to the  separable Hamiltonian, with the relative motion part given by ($r_{ij}=r$) $$\hat{H}_r=-\nabla^2_r + \frac{1}{4}\omega^2r^2+ \frac{1}{r},$$ plus a standard Harmonic Oscillator problem  for the center-of-mass motion.This system has analytical solutions in two and three dimensions ([<NAME> 1993 and 1994](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.48.3561)). We want to perform  a Variational Monte Carlo calculation of the ground state of two electrons in a quantum dot well with different oscillator energies, assuming total spin $S=0$.Our trial wave function has the following form $$\begin{equation}   \psi_{T}(\boldsymbol{r}_1,\boldsymbol{r}_2) =    C\exp{\left(-\alpha_1\omega(r_1^2+r_2^2)/2\right)}   \exp{\left(\frac{r_{12}}{(1+\alpha_2 r_{12})}\right)}, \label{eq:trial} \tag{5}\end{equation}$$ where the $\alpha$s represent our variational parameters, two in this case.Why does the trial function look like this? How did we get there? **This will be our main motivation** for switching toMachine Learning.To find an ansatz for the correlated part of the wave function, it is useful to rewrite the two-particlelocal energy in terms of the relative and center-of-mass motion.Let us denote the distance between the two electrons as$r_{12}$. We omit the center-of-mass motion since we are only interested in the case when$r_{12} \rightarrow 0$. The contribution from the center-of-mass (CoM) variable $\boldsymbol{R}_{\mathrm{CoM}}$gives only a finite contribution.We focus only on the terms that are relevant for $r_{12}$ and for three dimensions. The relevant local energy becomes then $$\lim_{r_{12} \rightarrow 0}E_L(R)=    \frac{1}{{\cal R}_T(r_{12})}\left(2\frac{d^2}{dr_{ij}^2}+\frac{4}{r_{ij}}\frac{d}{dr_{ij}}+\frac{2}{r_{ij}}-\frac{l(l+1)}{r_{ij}^2}+2E\right){\cal R}_T(r_{12}) = 0.$$ Set $l=0$ and we have the so-called **cusp** condition $$\frac{d {\cal R}_T(r_{12})}{dr_{12}} = -\frac{1}{2(l+1)}{\cal R}_T(r_{12})\qquad r_{12}\to 0$$ The above  results in $${\cal R}_T  \propto \exp{(r_{ij}/2)},$$ for anti-parallel spins and $${\cal R}_T  \propto \exp{(r_{ij}/4)},$$ for anti-parallel spins. This is the so-called cusp condition for the relative motion, resulting in a minimal requirementfor the correlation part of the wave fuction.For general systems containing more than say two electrons, we have thiscondition for each electron pair $ij$.<jupyter_code>python code here<jupyter_output><empty_output><jupyter_text>Importance samplingWe need to replace the brute forceMetropolis algorithm with a walk in coordinate space biased by the trial wave function.This approach is based on the Fokker-Planck equation and the Langevin equation for generating a trajectory in coordinate space.  The link between the Fokker-Planck equation and the Langevin equations are explained, only partly, in the slides below.An excellent reference on topics like Brownian motion, Markov chains, the Fokker-Planck equation and the Langevin equation is the text by  [<NAME>](http://www.elsevier.com/books/stochastic-processes-in-physics-and-chemistry/van-kampen/978-0-444-52965-7)Here we will focus first on the implementation part first.For a diffusion process characterized by a time-dependent probability density $P(x,t)$ in one dimension the Fokker-Planckequation reads (for one particle /walker) $$\frac{\partial P}{\partial t} = D\frac{\partial }{\partial x}\left(\frac{\partial }{\partial x} -F\right)P(x,t),$$ where $F$ is a drift term and $D$ is the diffusion coefficient. The new positions in coordinate space are given as the solutions of the Langevin equation using Euler's method, namely,we go from the Langevin equation $$\frac{\partial x(t)}{\partial t} = DF(x(t)) +\eta,$$ with $\eta$ a random variable,yielding a new position $$y = x+DF(x)\Delta t +\xi\sqrt{\Delta t},$$ where $\xi$ is gaussian random variable and $\Delta t$ is a chosen time step. The quantity $D$ is, in atomic units, equal to $1/2$ and comes from the factor $1/2$ in the kinetic energy operator. Note that $\Delta t$ is to be viewed as a parameter. Values of $\Delta t \in [0.001,0.01]$ yield in general rather stable values of the ground state energy.  The process of isotropic diffusion characterized by a time-dependent probability density $P(\mathbf{x},t)$ obeys (as an approximation) the so-called Fokker-Planck equation $$\frac{\partial P}{\partial t} = \sum_i D\frac{\partial }{\partial \mathbf{x_i}}\left(\frac{\partial }{\partial \mathbf{x_i}} -\mathbf{F_i}\right)P(\mathbf{x},t),$$ where $\mathbf{F_i}$ is the $i^{th}$ component of the drift term (drift velocity) caused by an external potential, and $D$ is the diffusion coefficient. The convergence to a stationary probability density can be obtained by setting the left hand side to zero. The resulting equation will be satisfied if and only if all the terms of the sum are equal zero, $$\frac{\partial^2 P}{\partial {\mathbf{x_i}^2}} = P\frac{\partial}{\partial {\mathbf{x_i}}}\mathbf{F_i} + \mathbf{F_i}\frac{\partial}{\partial {\mathbf{x_i}}}P.$$ The drift vector should be of the form $\mathbf{F} = g(\mathbf{x}) \frac{\partial P}{\partial \mathbf{x}}$. Then, $$\frac{\partial^2 P}{\partial {\mathbf{x_i}^2}} = P\frac{\partial g}{\partial P}\left( \frac{\partial P}{\partial {\mathbf{x}_i}}  \right)^2 + P g \frac{\partial ^2 P}{\partial {\mathbf{x}_i^2}}  + g \left( \frac{\partial P}{\partial {\mathbf{x}_i}}  \right)^2.$$ The condition of stationary density means that the left hand side equals zero. In other words, the terms containing first and second derivatives have to cancel each other. It is possible only if $g = \frac{1}{P}$, which yields $$\mathbf{F} = 2\frac{1}{\Psi_T}\nabla\Psi_T,$$ which is known as the so-called *quantum force*. This term is responsible for pushing the walker towards regions of configuration space where the trial wave function is large, increasing the efficiency of the simulation in contrast to the Metropolis algorithm where the walker has the same probability of moving in every direction.The Fokker-Planck equation yields a (the solution to the equation) transition probability given by the Green's function $$G(y,x,\Delta t) = \frac{1}{(4\pi D\Delta t)^{3N/2}} \exp{\left(-(y-x-D\Delta t F(x))^2/4D\Delta t\right)}$$ which in turn means that our brute force Metropolis algorithm $$A(y,x) = \mathrm{min}(1,q(y,x))),$$ with $q(y,x) = |\Psi_T(y)|^2/|\Psi_T(x)|^2$ is now replaced by the [Metropolis-Hastings algorithm](http://scitation.aip.org/content/aip/journal/jcp/21/6/10.1063/1.1699114) as well as [Hasting's article](http://biomet.oxfordjournals.org/content/57/1/97.abstract), $$q(y,x) = \frac{G(x,y,\Delta t)|\Psi_T(y)|^2}{G(y,x,\Delta t)|\Psi_T(x)|^2}$$<jupyter_code>add python code here<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11234195347415868, 'fraction_numerical': 0.012872913403873334, 'mean_word_length': 5.212624584717608, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 15, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '146395', 'license': 'CC0-1.0', 'n_tokens_mistral': 7248, 'n_tokens_neox': 6755, 'n_words': 3086, 'path': 'doc/pub/notebook1/ipynb/.ipynb_checkpoints/notebook1-checkpoint.ipynb', 'repo_name': 'GabrielSCabrera/ComputationalPhysics2', 'size': '26212'}