If we assume a discrete set of events, our initial probability distribution function can be given by
wi(0)=δi,0,
and its time-development after a given time step Δt=ϵ is
wi(t)=∑jW(j→i)wj(t=0).
The continuous analog to wi(0) is
w(x)→δ(x),
where we now have generalized the one-dimensional position x to a generic-dimensional vector x. The Kroenecker δ function is replaced by the δ distribution function δ(x) at t=0.
The transition from a state j to a state i is now replaced by a transition to a state with position y from a state with position x. The discrete sum of transition probabilities can then be replaced by an integral and we obtain the new distribution at a time t+Δt as
w(y,t+Δt)=∫W(y,t+Δt|x,t)w(x,t)dx,
and after m time steps we have
w(y,t+mΔt)=∫W(y,t+mΔt|x,t)w(x,t)dx.
When equilibrium is reached we have
w(y)=∫W(y|x,t)w(x)dx,
that is no time-dependence. Note our change of notation for W
We can solve the equation for w(y,t) by making a Fourier transform to momentum space. The PDF w(x,t) is related to its Fourier transform ˜w(k,t) through
w(x,t)=∫∞−∞dkexp(ikx)˜w(k,t),
and using the definition of the δ-function
δ(x)=12π∫∞−∞dkexp(ikx),
we see that
˜w(k,0)=1/2π.
We can then use the Fourier-transformed diffusion equation
∂˜w(k,t)∂t=−Dk2˜w(k,t),
with the obvious solution
˜w(k,t)=˜w(k,0)exp[−(Dk2t))=12πexp[−(Dk2t)].
With the Fourier transform we obtain
w(x,t)=∫∞−∞dkexp[ikx]12πexp[−(Dk2t)]=1√4πDtexp[−(x2/4Dt)],
with the normalization condition
∫∞−∞w(x,t)dx=1.
The solution represents the probability of finding our random walker at position x at time t if the initial distribution was placed at x=0 at t=0.
There is another interesting feature worth observing. The discrete transition probability W itself is given by a binomial distribution. The results from the central limit theorem state that transition probability in the limit n→∞ converges to the normal distribution. It is then possible to show that
W(il−jl,nϵ)→W(y,t+Δt|x,t)=1√4πDΔtexp[−((y−x)2/4DΔt)],
and that it satisfies the normalization condition and is itself a solution to the diffusion equation.
Let us now assume that we have three PDFs for times t0<t′<t, that is w(x0,t0), w(x′,t′) and w(x,t). We have then
w(x,t)=∫∞−∞W(x.t|x′.t′)w(x′,t′)dx′,
and
w(x,t)=∫∞−∞W(x.t|x0.t0)w(x0,t0)dx0,
and
w(x′,t′)=∫∞−∞W(x′.t′|x0,t0)w(x0,t0)dx0.
We can combine these equations and arrive at the famous Einstein-Smoluchenski-Kolmogorov-Chapman (ESKC) relation
W(xt|x0t0)=∫∞−∞W(x,t|x′,t′)W(x′,t′|x0,t0)dx′.
We can replace the spatial dependence with a dependence upon say the velocity (or momentum), that is we have
W(v,t|v0,t0)=∫∞−∞W(v,t|v′,t′)W(v′,t′|v0,t0)dx′.
We will now derive the Fokker-Planck equation. We start from the ESKC equation
W(x,t|x0,t0)=∫∞−∞W(x,t|x′,t′)W(x′,t′|x0,t0)dx′.
Define s=t′−t0, τ=t−t′ and t−t0=s+τ. We have then
W(x,s+τ|x0)=∫∞−∞W(x,τ|x′)W(x′,s|x0)dx′.
Assume now that τ is very small so that we can make an expansion in terms of a small step ξ, with x′=x−ξ, that is
W(x,s|x0)+∂W∂sτ+O(τ2)=∫∞−∞W(x,τ|x−ξ)W(x−ξ,s|x0)dx′.
We assume that W(x,τ|x−ξ) takes non-negligible values only when ξ is small. This is just another way of stating the Master equation!!
We say thus that x changes only by a small amount in the time interval τ. This means that we can make a Taylor expansion in terms of ξ, that is we expand
W(x,τ|x−ξ)W(x−ξ,s|x0)=∞∑n=0(−ξ)nn!∂n∂xn[W(x+ξ,τ|x)W(x,s|x0)].
We can then rewrite the ESKC equation as
∂W∂sτ=−W(x,s|x0)+∞∑n=0(−ξ)nn!∂n∂xn[W(x,s|x0)∫∞−∞ξnW(x+
ξ,τ|x)dξ].
We have neglected higher powers of τ and have used that for n=0 we get simply W(x,s|x0) due to normalization.
We say thus that x changes only by a small amount in the time interval τ. This means that we can make a Taylor expansion in terms of ξ, that is we expand
W(x,τ|x−ξ)W(x−ξ,s|x0)=∞∑n=0(−ξ)nn!∂n∂xn[W(x+ξ,τ|x)W(x,s|x0)].
We can then rewrite the ESKC equation as
∂W(x,s|x0)∂sτ=−W(x,s|x0)+∞∑n=0(−ξ)nn!∂n∂xn
[W(x,s|x0)∫∞−∞ξnW(x+ξ,τ|x)dξ].
We have neglected higher powers of τ and have used that for n=0 we get simply W(x,s|x0) due to normalization.
We simplify the above by introducing the moments
Mn=1τ∫∞−∞ξnW(x+ξ,τ|x)dξ=⟨[Δx(τ)]n⟩τ,
resulting in
∂W(x,s|x0)∂s=∞∑n=1(−ξ)nn!∂n∂xn[W(x,s|x0)Mn].
When τ→0 we assume that ⟨[Δx(τ)]n⟩→0 more rapidly than τ itself if n>2. When τ is much larger than the standard correlation time of system then Mn for n>2 can normally be neglected. This means that fluctuations become negligible at large time scales.
If we neglect such terms we can rewrite the ESKC equation as
∂W(x,s|x0)∂s=−∂M1W(x,s|x0)∂x+12∂2M2W(x,s|x0)∂x2.
In a more compact form we have
∂W∂s=−∂M1W∂x+12∂2M2W∂x2,
which is the Fokker-Planck equation! It is trivial to replace position with velocity (momentum).
Consider a particle suspended in a liquid. On its path through the liquid it will continuously collide with the liquid molecules. Because on average the particle will collide more often on the front side than on the back side, it will experience a systematic force proportional with its velocity, and directed opposite to its velocity. Besides this systematic force the particle will experience a stochastic force F(t). The equations of motion are
From hydrodynamics we know that the friction constant ξ is given by
ξ=6πηa/m
where η is the viscosity of the solvent and a is the radius of the particle .
Solving the second equation in the previous slide we get
v(t)=v0e−ξt+∫t0dτe−ξ(t−τ)F(τ).
If we want to get some useful information out of this, we have to average over all possible realizations of F(t), with the initial velocity as a condition. A useful quantity for example is
⟨v(t)⋅v(t)⟩v0=v−ξ2t0+2∫t0dτe−ξ(2t−τ)v0⋅⟨F(τ)⟩v0
+∫t0dτ′∫t0dτe−ξ(2t−τ−τ′)⟨F(τ)⋅F(τ′)⟩v0.
In order to continue we have to make some assumptions about the conditional averages of the stochastic forces. In view of the chaotic character of the stochastic forces the following assumptions seem to be appropriate
⟨F(t)⟩=0,
and
⟨F(t)⋅F(t′)⟩v0=Cv0δ(t−t′).
We omit the subscript v0, when the quantity of interest turns out to be independent of v0. Using the last three equations we get
⟨v(t)⋅v(t)⟩v0=v20e−2ξt+Cv02ξ(1−e−2ξt).
For large t this should be equal to 3kT/m, from which it follows that
⟨F(t)⋅F(t′)⟩=6kTmξδ(t−t′).
This result is called the fluctuation-dissipation theorem .
Integrating
v(t)=v0e−ξt+∫t0dτe−ξ(t−τ)F(τ),
we get
r(t)=r0+v01ξ(1−e−ξt)+∫t0dτ∫τ0τ′e−ξ(τ−τ′)F(τ′),
from which we calculate the mean square displacement
⟨(r(t)−r0)2⟩v0=v20ξ(1−e−ξt)2+3kTmξ2(2ξt−3+4e−ξt−e−2ξt).
For very large t this becomes
⟨(r(t)−r0)2⟩=6kTmξt
from which we get the Einstein relation
D=kTmξ
where we have used ⟨(r(t)−r0)2⟩=6Dt.
The general derivative formula of the Jastrow factor (or the ansatz for the correlated part of the wave function) is (the subscript C stands for Correlation)
1ΨC∂ΨC∂xk=k−1∑i=1∂gik∂xk+N∑i=k+1∂gki∂xk
However, with our written in way which can be reused later as
ΨC=∏i<jg(rij)=exp{∑i<jf(rij)},
the gradient needed for the quantum force and local energy is easy to compute. The function f(rij) will depends on the system under study. In the equations below we will keep this general form.
In the Metropolis/Hasting algorithm, the acceptance ratio determines the probability for a particle to be accepted at a new position. The ratio of the trial wave functions evaluated at the new and current positions is given by (OB for the onebody part)
R≡ΨnewTΨoldT=ΨnewOBΨoldOBΨnewCΨoldC
Here ΨOB is our onebody part (Slater determinant or product of boson single-particle states) while ΨC is our correlation function, or Jastrow factor. We need to optimize the ∇ΨT/ΨT ratio and the second derivative as well, that is the ∇2ΨT/ΨT ratio. The first is needed when we compute the so-called quantum force in importance sampling. The second is needed when we compute the kinetic energy term of the local energy.
∇ΨΨ=∇(ΨOBΨC)ΨOBΨC=ΨC∇ΨOB+ΨOB∇ΨCΨOBΨC=∇ΨOBΨOB+∇ΨCΨC
The expectation value of the kinetic energy expressed in scaled units for particle i is
⟨ˆKi⟩=−12⟨Ψ|∇2i|Ψ⟩⟨Ψ|Ψ⟩,
ˆKi=−12∇2iΨΨ.
The second derivative which enters the definition of the local energy is
∇2ΨΨ=∇2ΨOBΨOB+∇2ΨCΨC+2∇ΨOBΨOB⋅∇ΨCΨC
We discuss here how to calculate these quantities in an optimal way.
We have defined the correlated function as
ΨC=∏i<jg(rij)=N∏i<jg(rij)=N∏i=1N∏j=i+1g(rij),
with rij=|ri−rj|=√(xi−xj)2+(yi−yj)2+(zi−zj)2 in three dimensions or rij=|ri−rj|=√(xi−xj)2+(yi−yj)2 if we work with two-dimensional systems.
In our particular case we have
ΨC=∏i<jg(rij)=exp{∑i<jf(rij)}.
The total number of different relative distances rij is N(N−1)/2. In a matrix storage format, the relative distances form a strictly upper triangular matrix
r≡(0r1,2r1,3⋯r1,N⋮0r2,3⋯r2,N⋮⋮0⋱⋮⋮⋮⋮⋱rN−1,N000⋯0).
This applies to g=g(rij) as well.
In our algorithm we will move one particle at the time, say the kth-particle. This sampling will be seen to be particularly efficient when we are going to compute a Slater determinant.
We have that the ratio between Jastrow factors RC is given by
RC=ΨnewCΨcurC=k−1∏i=1gnewikgcurikN∏i=k+1gnewkigcurki.
For the Pade-Jastrow form
RC=ΨnewCΨcurC=expUnewexpUcur=expΔU,
where
ΔU=k−1∑i=1(fnewik−fcurik)+N∑i=k+1(fnewki−fcurki)
One needs to develop a special algorithm that runs only through the elements of the upper triangular matrix g and have k as an index.
The expression to be derived in the following is of interest when computing the quantum force and the kinetic energy. It has the form
∇iΨCΨC=1ΨC∂ΨC∂xi,
for all dimensions and with i running over all particles.
For the first derivative only N−1 terms survive the ratio because the g-terms that are not differentiated cancel with their corresponding ones in the denominator. Then,
1ΨC∂ΨC∂xk=k−1∑i=11gik∂gik∂xk+N∑i=k+11gki∂gki∂xk.
An equivalent equation is obtained for the exponential form after replacing gij by exp(fij), yielding:
1ΨC∂ΨC∂xk=k−1∑i=1∂gik∂xk+N∑i=k+1∂gki∂xk,
with both expressions scaling as O(N).
Using the identity
∂∂xigij=−∂∂xjgij,
we get expressions where all the derivatives acting on the particle are represented by the second index of g:
1ΨC∂ΨC∂xk=k−1∑i=11gik∂gik∂xk−N∑i=k+11gki∂gki∂xi,
and for the exponential case:
1ΨC∂ΨC∂xk=k−1∑i=1∂gik∂xk−N∑i=k+1∂gki∂xi.
For correlation forms depending only on the scalar distances rij we can use the chain rule. Noting that
∂gij∂xj=∂gij∂rij∂rij∂xj=xj−xirij∂gij∂rij,
we arrive at
1ΨC∂ΨC∂xk=k−1∑i=11gikrikrik∂gik∂rik−N∑i=k+11gkirkirki∂gki∂rki.
Note that for the Pade-Jastrow form we can set gij≡g(rij)=ef(rij)=efij and
∂gij∂rij=gij∂fij∂rij.
Therefore,
1ΨC∂ΨC∂xk=k−1∑i=1rikrik∂fik∂rik−N∑i=k+1rkirki∂fki∂rki,
where
rij=|rj−ri|=(xj−xi)e1+(yj−yi)e2+(zj−zi)e3
is the relative distance.
The second derivative of the Jastrow factor divided by the Jastrow factor (the way it enters the kinetic energy) is
[∇2ΨCΨC]x= 2N∑k=1k−1∑i=1∂2gik∂x2k + N∑k=1(k−1∑i=1∂gik∂xk−N∑i=k+1∂gki∂xi)2
But we have a simple form for the function, namely
ΨC=∏i<jexpf(rij),
and it is easy to see that for particle k we have
∇2kΨCΨC=∑ij≠k(rk−ri)(rk−rj)rkirkjf′(rki)f′(rkj)+∑j≠k(f″
Our aim with this part of the project is to be able to
To achieve this will look at methods like the simplest possible gradient descent, Steepest descent, the conjugate gradient method and stochastic gradient descent. These methods allow us to find the minima of a multivariable function like our energy (function of several variational parameters). Alternatively, you can always use Newton's method. In particular, since we will normally have one variational parameter, Newton's method can be easily used in finding the minimum of the local energy.
Let us illustrate what is needed in our calculations using a simple example, the harmonic oscillator in one dimension. For the harmonic oscillator in one-dimension we have a trial wave function and probability
\psi_T(x;\alpha) = \exp{-(\frac{1}{2}\alpha^2x^2)},
which results in a local energy
\frac{1}{2}\left(\alpha^2+x^2(1-\alpha^4)\right).
We can compare our numerically calculated energies with the exact energy as function of \alpha
\overline{E}[\alpha] = \frac{1}{4}\left(\alpha^2+\frac{1}{\alpha^2}\right).
The derivative of the energy with respect to \alpha gives
\begin{equation*}
\frac{d\langle E_L[\alpha]\rangle}{d\alpha} = \frac{1}{2}\alpha-\frac{1}{2\alpha^3}
\end{equation*}
and a second derivative which is always positive (meaning that we find a minimum)
\begin{equation*}
\frac{d^2\langle E_L[\alpha]\rangle}{d\alpha^2} = \frac{1}{2}+\frac{3}{2\alpha^4}
\end{equation*}
The condition
\begin{equation*}
\frac{d\langle E_L[\alpha]\rangle}{d\alpha} = 0,
\end{equation*}
gives the optimal \alpha=1 , as expected.
a) Derive the local energy for the harmonic oscillator in one dimension and find its expectation value.
b) Show also that the optimal value of optimal \alpha=1
c) Repeat the above steps in two dimensions for N bosons or electrons. What is the optimal value of \alpha ?
We can also minimize the variance. In our simple model the variance is
\sigma^2[\alpha]=\frac{1}{4}\left(1+(1-\alpha^4)^2\frac{3}{4\alpha^4}\right)-\overline{E}^2.
which yields a second derivative which is always positive.
In general we end up computing the expectation value of the energy in terms of some parameters \alpha_0,\alpha_1,\dots,\alpha_n and we search for a minimum in this multi-variable parameter space. This leads to an energy minimization problem where we need the derivative of the energy as a function of the variational parameters.
In the above example this was easy and we were able to find the expression for the derivative by simple derivations. However, in our actual calculations the energy is represented by a multi-dimensional integral with several variational parameters. How can we can then obtain the derivatives of the energy with respect to the variational parameters without having to resort to expensive numerical derivations?
To find the derivatives of the local energy expectation value as function of the variational parameters, we can use the chain rule and the hermiticity of the Hamiltonian.
Let us define
\bar{E}_{\alpha}=\frac{d\langle E_L[\alpha]\rangle}{d\alpha}.
as the derivative of the energy with respect to the variational parameter \alpha (we limit ourselves to one parameter only). In the above example this was easy and we obtain a simple expression for the derivative. We define also the derivative of the trial function (skipping the subindex T ) as
\bar{\psi}_{\alpha}=\frac{d\psi[\alpha]\rangle}{d\alpha}.
The elements of the gradient of the local energy are then (using the chain rule and the hermiticity of the Hamiltonian)
\bar{E}_{\alpha} = 2\left( \langle \frac{\bar{\psi}_{\alpha}}{\psi[\alpha]}E_L[\alpha]\rangle -\langle \frac{\bar{\psi}_{\alpha}}{\psi[\alpha]}\rangle\langle E_L[\alpha] \rangle\right).
From a computational point of view it means that you need to compute the expectation values of
\langle \frac{\bar{\psi}_{\alpha}}{\psi[\alpha]}E_L[\alpha]\rangle,
and
\langle \frac{\bar{\psi}_{\alpha}}{\psi[\alpha]}\rangle\langle E_L[\alpha]\rangle
a) Show that
\bar{E}_{\alpha} = 2\left( \langle \frac{\bar{\psi}_{\alpha}}{\psi[\alpha]}E_L[\alpha]\rangle -\langle \frac{\bar{\psi}_{\alpha}}{\psi[\alpha]}\rangle\langle E_L[\alpha] \rangle\right).
b) Find the corresponding expression for the variance.
# 2-electron VMC code for 2dim quantum dot with importance sampling
# Using gaussian rng for new positions and Metropolis- Hastings
# Added energy minimization with gradient descent using fixed step size
# To do: replace with optimization codes from scipy and/or use stochastic gradient descent
from math import exp, sqrt
from random import random, seed, normalvariate
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import sys
# Trial wave function for the 2-electron quantum dot in two dims
def WaveFunction(r,alpha,beta):
r1 = r[0,0]**2 + r[0,1]**2
r2 = r[1,0]**2 + r[1,1]**2
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = r12/(1+beta*r12)
return exp(-0.5*alpha*(r1+r2)+deno)
# Local energy for the 2-electron quantum dot in two dims, using analytical local energy
def LocalEnergy(r,alpha,beta):
r1 = (r[0,0]**2 + r[0,1]**2)
r2 = (r[1,0]**2 + r[1,1]**2)
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = 1.0/(1+beta*r12)
deno2 = deno*deno
return 0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha + 1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
# Derivate of wave function ansatz as function of variational parameters
def DerivativeWFansatz(r,alpha,beta):
WfDer = np.zeros((2), np.double)
r1 = (r[0,0]**2 + r[0,1]**2)
r2 = (r[1,0]**2 + r[1,1]**2)
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = 1.0/(1+beta*r12)
deno2 = deno*deno
WfDer[0] = -0.5*(r1+r2)
WfDer[1] = -r12*r12*deno2
return WfDer
# Setting up the quantum force for the two-electron quantum dot, recall that it is a vector
def QuantumForce(r,alpha,beta):
qforce = np.zeros((NumberParticles,Dimension), np.double)
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = 1.0/(1+beta*r12)
qforce[0,:] = -2*r[0,:]*alpha*(r[0,:]-r[1,:])*deno*deno/r12
qforce[1,:] = -2*r[1,:]*alpha*(r[1,:]-r[0,:])*deno*deno/r12
return qforce
# Computing the derivative of the energy and the energy
def EnergyMinimization(alpha, beta):
NumberMCcycles= 10000
# Parameters in the Fokker-Planck simulation of the quantum force
D = 0.5
TimeStep = 0.05
# positions
PositionOld = np.zeros((NumberParticles,Dimension), np.double)
PositionNew = np.zeros((NumberParticles,Dimension), np.double)
# Quantum force
QuantumForceOld = np.zeros((NumberParticles,Dimension), np.double)
QuantumForceNew = np.zeros((NumberParticles,Dimension), np.double)
# seed for rng generator
seed()
energy = 0.0
DeltaE = 0.0
EnergyDer = np.zeros((2), np.double)
DeltaPsi = np.zeros((2), np.double)
DerivativePsiE = np.zeros((2), np.double)
#Initial position
for i in range(NumberParticles):
for j in range(Dimension):
PositionOld[i,j] = normalvariate(0.0,1.0)*sqrt(TimeStep)
wfold = WaveFunction(PositionOld,alpha,beta)
QuantumForceOld = QuantumForce(PositionOld,alpha, beta)
#Loop over MC MCcycles
for MCcycle in range(NumberMCcycles):
#Trial position moving one particle at the time
for i in range(NumberParticles):
for j in range(Dimension):
PositionNew[i,j] = PositionOld[i,j]+normalvariate(0.0,1.0)*sqrt(TimeStep)+\
QuantumForceOld[i,j]*TimeStep*D
wfnew = WaveFunction(PositionNew,alpha,beta)
QuantumForceNew = QuantumForce(PositionNew,alpha, beta)
GreensFunction = 0.0
for j in range(Dimension):
GreensFunction += 0.5*(QuantumForceOld[i,j]+QuantumForceNew[i,j])*\
(D*TimeStep*0.5*(QuantumForceOld[i,j]-QuantumForceNew[i,j])-\
PositionNew[i,j]+PositionOld[i,j])
GreensFunction = exp(GreensFunction)
ProbabilityRatio = GreensFunction*wfnew**2/wfold**2
#Metropolis-Hastings test to see whether we accept the move
if random() <= ProbabilityRatio:
for j in range(Dimension):
PositionOld[i,j] = PositionNew[i,j]
QuantumForceOld[i,j] = QuantumForceNew[i,j]
wfold = wfnew
DeltaE = LocalEnergy(PositionOld,alpha,beta)
DerPsi = DerivativeWFansatz(PositionOld,alpha,beta)
DeltaPsi += DerPsi
energy += DeltaE
DerivativePsiE += DerPsi*DeltaE
# We calculate mean values
energy /= NumberMCcycles
DerivativePsiE /= NumberMCcycles
DeltaPsi /= NumberMCcycles
EnergyDer = 2*(DerivativePsiE-DeltaPsi*energy)
return energy, EnergyDer
#Here starts the main program with variable declarations
NumberParticles = 2
Dimension = 2
# guess for variational parameters
alpha = 0.9
beta = 0.2
# Set up iteration using gradient descent method
Energy = 0
EDerivative = np.zeros((2), np.double)
eta = 0.01
Niterations = 50
#
for iter in range(Niterations):
Energy, EDerivative = EnergyMinimization(alpha,beta)
alphagradient = EDerivative[0]
betagradient = EDerivative[1]
alpha -= eta*alphagradient
beta -= eta*betagradient
print(alpha, beta)
print(Energy, EDerivative[0], EDerivative[1])
The following function uses the above described BFGS algorithm. Here we have defined a function which calculates the energy and a function which computes the first derivative.
# 2-electron VMC code for 2dim quantum dot with importance sampling
# Using gaussian rng for new positions and Metropolis- Hastings
# Added energy minimization using the BFGS algorithm, see p. 136 of https://www.springer.com/it/book/9780387303031
from math import exp, sqrt
from random import random, seed, normalvariate
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from scipy.optimize import minimize
import sys
# Trial wave function for the 2-electron quantum dot in two dims
def WaveFunction(r,alpha,beta):
r1 = r[0,0]**2 + r[0,1]**2
r2 = r[1,0]**2 + r[1,1]**2
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = r12/(1+beta*r12)
return exp(-0.5*alpha*(r1+r2)+deno)
# Local energy for the 2-electron quantum dot in two dims, using analytical local energy
def LocalEnergy(r,alpha,beta):
r1 = (r[0,0]**2 + r[0,1]**2)
r2 = (r[1,0]**2 + r[1,1]**2)
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = 1.0/(1+beta*r12)
deno2 = deno*deno
return 0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha + 1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
# Derivate of wave function ansatz as function of variational parameters
def DerivativeWFansatz(r,alpha,beta):
WfDer = np.zeros((2), np.double)
r1 = (r[0,0]**2 + r[0,1]**2)
r2 = (r[1,0]**2 + r[1,1]**2)
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = 1.0/(1+beta*r12)
deno2 = deno*deno
WfDer[0] = -0.5*(r1+r2)
WfDer[1] = -r12*r12*deno2
return WfDer
# Setting up the quantum force for the two-electron quantum dot, recall that it is a vector
def QuantumForce(r,alpha,beta):
qforce = np.zeros((NumberParticles,Dimension), np.double)
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = 1.0/(1+beta*r12)
qforce[0,:] = -2*r[0,:]*alpha*(r[0,:]-r[1,:])*deno*deno/r12
qforce[1,:] = -2*r[1,:]*alpha*(r[1,:]-r[0,:])*deno*deno/r12
return qforce
# Computing the derivative of the energy and the energy
def EnergyDerivative(x0):
# Parameters in the Fokker-Planck simulation of the quantum force
D = 0.5
TimeStep = 0.05
NumberMCcycles= 10000
# positions
PositionOld = np.zeros((NumberParticles,Dimension), np.double)
PositionNew = np.zeros((NumberParticles,Dimension), np.double)
# Quantum force
QuantumForceOld = np.zeros((NumberParticles,Dimension), np.double)
QuantumForceNew = np.zeros((NumberParticles,Dimension), np.double)
energy = 0.0
DeltaE = 0.0
alpha = x0[0]
beta = x0[1]
EnergyDer = 0.0
DeltaPsi = 0.0
DerivativePsiE = 0.0
#Initial position
for i in range(NumberParticles):
for j in range(Dimension):
PositionOld[i,j] = normalvariate(0.0,1.0)*sqrt(TimeStep)
wfold = WaveFunction(PositionOld,alpha,beta)
QuantumForceOld = QuantumForce(PositionOld,alpha, beta)
#Loop over MC MCcycles
for MCcycle in range(NumberMCcycles):
#Trial position moving one particle at the time
for i in range(NumberParticles):
for j in range(Dimension):
PositionNew[i,j] = PositionOld[i,j]+normalvariate(0.0,1.0)*sqrt(TimeStep)+\
QuantumForceOld[i,j]*TimeStep*D
wfnew = WaveFunction(PositionNew,alpha,beta)
QuantumForceNew = QuantumForce(PositionNew,alpha, beta)
GreensFunction = 0.0
for j in range(Dimension):
GreensFunction += 0.5*(QuantumForceOld[i,j]+QuantumForceNew[i,j])*\
(D*TimeStep*0.5*(QuantumForceOld[i,j]-QuantumForceNew[i,j])-\
PositionNew[i,j]+PositionOld[i,j])
GreensFunction = exp(GreensFunction)
ProbabilityRatio = GreensFunction*wfnew**2/wfold**2
#Metropolis-Hastings test to see whether we accept the move
if random() <= ProbabilityRatio:
for j in range(Dimension):
PositionOld[i,j] = PositionNew[i,j]
QuantumForceOld[i,j] = QuantumForceNew[i,j]
wfold = wfnew
DeltaE = LocalEnergy(PositionOld,alpha,beta)
DerPsi = DerivativeWFansatz(PositionOld,alpha,beta)
DeltaPsi += DerPsi
energy += DeltaE
DerivativePsiE += DerPsi*DeltaE
# We calculate mean values
energy /= NumberMCcycles
DerivativePsiE /= NumberMCcycles
DeltaPsi /= NumberMCcycles
EnergyDer = 2*(DerivativePsiE-DeltaPsi*energy)
return EnergyDer
# Computing the expectation value of the local energy
def Energy(x0):
# Parameters in the Fokker-Planck simulation of the quantum force
D = 0.5
TimeStep = 0.05
# positions
PositionOld = np.zeros((NumberParticles,Dimension), np.double)
PositionNew = np.zeros((NumberParticles,Dimension), np.double)
# Quantum force
QuantumForceOld = np.zeros((NumberParticles,Dimension), np.double)
QuantumForceNew = np.zeros((NumberParticles,Dimension), np.double)
energy = 0.0
DeltaE = 0.0
alpha = x0[0]
beta = x0[1]
NumberMCcycles= 10000
#Initial position
for i in range(NumberParticles):
for j in range(Dimension):
PositionOld[i,j] = normalvariate(0.0,1.0)*sqrt(TimeStep)
wfold = WaveFunction(PositionOld,alpha,beta)
QuantumForceOld = QuantumForce(PositionOld,alpha, beta)
#Loop over MC MCcycles
for MCcycle in range(NumberMCcycles):
#Trial position moving one particle at the time
for i in range(NumberParticles):
for j in range(Dimension):
PositionNew[i,j] = PositionOld[i,j]+normalvariate(0.0,1.0)*sqrt(TimeStep)+\
QuantumForceOld[i,j]*TimeStep*D
wfnew = WaveFunction(PositionNew,alpha,beta)
QuantumForceNew = QuantumForce(PositionNew,alpha, beta)
GreensFunction = 0.0
for j in range(Dimension):
GreensFunction += 0.5*(QuantumForceOld[i,j]+QuantumForceNew[i,j])*\
(D*TimeStep*0.5*(QuantumForceOld[i,j]-QuantumForceNew[i,j])-\
PositionNew[i,j]+PositionOld[i,j])
GreensFunction = exp(GreensFunction)
ProbabilityRatio = GreensFunction*wfnew**2/wfold**2
#Metropolis-Hastings test to see whether we accept the move
if random() <= ProbabilityRatio:
for j in range(Dimension):
PositionOld[i,j] = PositionNew[i,j]
QuantumForceOld[i,j] = QuantumForceNew[i,j]
wfold = wfnew
DeltaE = LocalEnergy(PositionOld,alpha,beta)
energy += DeltaE
# We calculate mean values
energy /= NumberMCcycles
return energy
#Here starts the main program with variable declarations
NumberParticles = 2
Dimension = 2
# seed for rng generator
seed()
# guess for variational parameters
x0 = np.array([0.9,0.2])
# Using Broydens method
res = minimize(Energy, x0, method='BFGS', jac=EnergyDerivative, options={'gtol': 1e-4,'disp': True})
print(res.x)
Note that the minimize function returns the finale values for the variable \alpha=x0[0] and \beta=x0[1] in the array x .
Let us quickly remind ourselves how we derive the above method.
Perhaps the most celebrated of all one-dimensional root-finding routines is Newton's method, also called the Newton-Raphson method. This method requires the evaluation of both the function f and its derivative f' at arbitrary points. If you can only calculate the derivative numerically and/or your function is not of the smooth type, we normally discourage the use of this method.
The Newton-Raphson formula consists geometrically of extending the tangent line at a current point until it crosses zero, then setting the next guess to the abscissa of that zero-crossing. The mathematics behind this method is rather simple. Employing a Taylor expansion for x sufficiently close to the solution s , we have
f(s)=0=f(x)+(s-x)f'(x)+\frac{(s-x)^2}{2}f''(x) +\dots.
\tag{1}
For small enough values of the function and for well-behaved functions, the terms beyond linear are unimportant, hence we obtain
f(x)+(s-x)f'(x)\approx 0,
yielding
s\approx x-\frac{f(x)}{f'(x)}.
Having in mind an iterative procedure, it is natural to start iterating with
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}.
The above is Newton-Raphson's method. It has a simple geometric interpretation, namely x_{n+1} is the point where the tangent from (x_n,f(x_n)) crosses the x -axis. Close to the solution, Newton-Raphson converges fast to the desired result. However, if we are far from a root, where the higher-order terms in the series are important, the Newton-Raphson formula can give grossly inaccurate results. For instance, the initial guess for the root might be so far from the true root as to let the search interval include a local maximum or minimum of the function. If an iteration places a trial guess near such a local extremum, so that the first derivative nearly vanishes, then Newton-Raphson may fail totally
Newton's method can be generalized to systems of several non-linear equations and variables. Consider the case with two equations
\begin{array}{cc} f_1(x_1,x_2) &=0\\
f_2(x_1,x_2) &=0,\end{array}.
We Taylor expand to obtain
\begin{array}{cc} 0=f_1(x_1+h_1,x_2+h_2)=&f_1(x_1,x_2)+h_1
\partial f_1/\partial x_1+h_2
\partial f_1/\partial x_2+\dots\\
0=f_2(x_1+h_1,x_2+h_2)=&f_2(x_1,x_2)+h_1
\partial f_2/\partial x_1+h_2
\partial f_2/\partial x_2+\dots
\end{array}.
Defining the Jacobian matrix \hat{J} we have
\hat{J}=\left( \begin{array}{cc}
\partial f_1/\partial x_1 & \partial f_1/\partial x_2 \\
\partial f_2/\partial x_1 &\partial f_2/\partial x_2
\end{array} \right),
we can rephrase Newton's method as
\left(\begin{array}{c} x_1^{n+1} \\ x_2^{n+1} \end{array} \right)=
\left(\begin{array}{c} x_1^{n} \\ x_2^{n} \end{array} \right)+
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right),
where we have defined
\left(\begin{array}{c} h_1^{n} \\ h_2^{n} \end{array} \right)=
-{\bf \hat{J}}^{-1}
\left(\begin{array}{c} f_1(x_1^{n},x_2^{n}) \\ f_2(x_1^{n},x_2^{n}) \end{array} \right).
We need thus to compute the inverse of the Jacobian matrix and it is to understand that difficulties may arise in case \hat{J} is nearly singular.
It is rather straightforward to extend the above scheme to systems of more than two non-linear equations. In our case, the Jacobian matrix is given by the Hessian that represents the second derivative of cost function.