Diffusion quantum Monte Carlo method

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

1D harmonic oscillator

Starting wave function, local energy, diffusion coefficient

ϕ0[x_,α_]:=Exp[-α*x*x];
exactEnergy[α_]=Assuming[α>0,Integrate[ϕ0[x,α]*(-D[ϕ0[x,α],x,x]+x^2 ϕ0[x,α])/2,{x,-Infinity,Infinity}]/Integrate[ϕ0[x,α]^2,{x,-Infinity,Infinity}]];
locEn[x_,α_]=Simplify[1/2(-D[ϕ0[x,α],x,x]+x^2 ϕ0[x,α])/ϕ0[x,α]];
Dif[x_,α_]=D[ϕ0[x,α],x]/ϕ0[x,α];
Print["ϕ0(x,α) = ", ϕ0[x,α],", ϵ(x) = ", locEn[x,α],", D(x) = ", Dif[x,α]]

Get initial n points with distribution given by |ϕ0(x)|^2

αini=1.0; (* exact for α = 1/2 *)
nini=1000000;
a=-5.0;b=5.0;
ρ[x_]:=ϕ0[x,αini]^2;
x0=0.0;δ=4.0;
(* use in Matropolis each 10th point of random walk to avoid correlation*)
points=MyRandomMetropolis1D[nini,ρ,a,b,x0,δ,10]; 
ϵ[x_]:=locEn[x,αini];
energy=MyMCIntegration[points,ϵ];
Print["Alpha = ", αini, ", mean energy = ", energy[[1]],", exact energy = ", exactEnergy[αini]]

Histogram[points[[1;;nini]],{a,b,(b-a)/100},"PDF",Epilog->First@Plot[ρ[x],{x,a,b},PlotRange->All,PlotStyle->Red]]

Time evolution of random points:

nt=200; (* number of iterations *)
dt=0.01; (* time step *)
En=energy[[1]];

xx=ConstantArray[0.0,2*nini];
xx[[1;;nini]]=points;
nn=nini;
meanEnergy=ConstantArray[0.0,nt];
Do[
j=1;
While[j≤nn,
xx[[j]]=xx[[j]]+Dif[xx[[j]],αini]*dt+MyRandomNormal[{0.0,Sqrt[dt]}];
expEn=Exp[-(locEn[xx[[j]],αini]-En)*dt];
If[
expEn<1.0,
(* then *)
If[expEn<RandomReal[], (* delete the point (replace with the last one) *)
xx[[j]]=xx[[nn]];
nn--;
j--
],
(* else *)
If[expEn-1.0>RandomReal[], (* add another point *)
xx[[nn+1]]=xx[[j]];
nn++
]
];
j++
];
{En,error}=MyMCIntegration[xx[[1;;nn]],ϵ];
meanEnergy[[i]]=En;
If[Mod[i,10]==0,Print["Iteration ", i,", En = ", En,", number of points: ",nn]],
{i,1,nt}
];

Print["Average energy ",Mean[meanEnergy[[nt/2;;nt]]]];
ListPlot[meanEnergy]

Histogram[{xx[[1;;nn]],points[[1;;nini]]},{a,b,(b-a)/100},"PDF",Epilog->First@Plot[ρ[x],{x,a,b},PlotRange->All,PlotStyle->Red]]