Variational Monte Carlo method

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

1D harmonic oscillator

Trial function which is exact for α = 0.5:

Clear[α];
ψ[x_]:=Exp[-α*x*x];
exactEnergy=Assuming[α>0,Integrate[ψ[x]*(-D[ψ[x],x,x]+x^2 ψ[x])/2,{x,-Infinity,Infinity}]/Integrate[ψ[x]^2,{x,-Infinity,Infinity}]]

n=100000;
a=-8.0;b=8.0;
Do[
α=0.2+0.05*i;
ρ[x_]:=ψ[x]^2;
x0=0.5;δ=0.5;
points=MyRandomMetropolis1D[n,ρ,a,b,x0,δ,1];
ϵ[x_]:=(α+x*x*(0.5-2.0*α*α));
{energy,error}=MyMCIntegration[points,ϵ];
Print["Alpha = ", α, ", energy = ", energy,", exact energy = ", exactEnergy, ", probable error = ", error],
{i,0,11}
];

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

Hydrogen atom

Trial function which is exact for α = 1.0:

Clear[α];
ψ[x_,y_,z_,α_]:=Exp[-α*Sqrt[x^2+y^2+z^2]];
ϵloc[x_,y_,z_,α_]:=Evaluate[(-1/2(D[ψ[x,y,z,α],x,x]+D[ψ[x,y,z,α],y,y]+D[ψ[x,y,z,α],z,z])-ψ[x,y,z,α]/Sqrt[x^2+y^2+z^2])/ψ[x,y,z,α]//Simplify];
ϵloc[x,y,z,α]

n=100000;
dim=3;
xmax=8.0;xmin=-xmax;
a={xmin,xmin,xmin};b={xmax,xmax,xmax};
Do[
αα=0.8+0.1*i;
ρ[x_]:=ψ[x[[1]],x[[2]],x[[3]],αα]^2;
(* Print[{N[ρ[a]],N[ρ[b]]}]; *)
x0=b/10.0;δ=b/10;
points=MyRandomMetropolisND[n,dim,ρ,a,b,x0,δ,1];
ϵ[x_]:=ϵloc[x[[1]],x[[2]],x[[3]],αα];
{energy,error}=MyMCIntegration[points,ϵ];
Print["Alpha = ", αα, ", energy = ", energy, ", probable error = ", error],
{i,0,4}
];

Trial function consisting of one Gaussian function:

Clear[α];
ψ[x_,y_,z_,α_]:=Exp[-α*(x^2+y^2+z^2)];
ϵloc[x_,y_,z_,α_]:=Evaluate[(-1/2(D[ψ[x,y,z,α],x,x]+D[ψ[x,y,z,α],y,y]+D[ψ[x,y,z,α],z,z])-ψ[x,y,z,α]/Sqrt[x^2+y^2+z^2])/ψ[x,y,z,α]//Simplify];
ϵloc[x,y,z,α]

n=10000;
dim=3;
xmax=8.0;xmin=-xmax;
a={xmin,xmin,xmin};b={xmax,xmax,xmax};
nalpha=21;
alpha=ConstantArray[0.0,nalpha];
energy=ConstantArray[0.0,nalpha];
Do[
alpha[[i]]=0.1+0.02*(i-1);
ρ[x_]:=ψ[x[[1]],x[[2]],x[[3]],alpha[[i]]]^2;
(* Print[{N[ρ[a]],N[ρ[b]]}]; *)
x0=b/10.0;δ=b/10;
points=MyRandomMetropolisND[n,dim,ρ,a,b,x0,δ,1];
ϵ[x_]:=ϵloc[x[[1]],x[[2]],x[[3]],alpha[[i]]];
{energy[[i]],error}=MyMCIntegration[points,ϵ];
Print["Alpha = ", alpha[[i]], ", energy = ", energy[[i]], ", probable error = ", error],
{i,1,nalpha}
];

plot1=ListPlot[Transpose[{alpha,energy}]]

Finding minimum using the least-square method:

fToFit[α_,par_]:=par[[1]]+par[[2]]*(α-par[[3]])^2;
f[par_]:=Module[
{i,fx,squares},
squares=0.0;
Do[
fx=fToFit[alpha[[i]],par]-energy[[i]];
squares+=fx*fx,
{i,1,nalpha}
];
Return[squares]
];
{minpar,value,niter}=MyMinimization[{-0.42,1.0,0.3}, {0.01,0.01,0.01}, f, 10^(-6)]
plot2=Plot[fToFit[x,minpar],{x,0.1,0.5}];
Show[plot1,plot2]

alphaSTO=1;
alphaGTO=minpar[[3]];
norm1=1/Sqrt[Integrate[Exp[-2*alphaSTO*r],{r,0,Infinity}]];
norm2=1/Sqrt[Integrate[Exp[-2*alphaGTO*r^2],{r,0,Infinity}]];
Plot[{Exp[-alphaSTO*r]*norm1,Exp[-alphaGTO*r^2]*norm2},{r,0,10},PlotRange->All]

Trial function consisting of two Gaussian function with minimization:

ψ[x_,y_,z_,α1_,α2_,c2_]:=Exp[-α1*(x^2+y^2+z^2)]+c2 *Exp[-α2*(x^2+y^2+z^2)];
ϵloc[x_,y_,z_,α1_,α2_,c2_]:=Evaluate[(-1/2(D[ψ[x,y,z,α1,α2,c2],x,x]+D[ψ[x,y,z,α1,α2,c2],y,y]+D[ψ[x,y,z,α1,α2,c2],z,z])-ψ[x,y,z,α1,α2,c2]/Sqrt[x^2+y^2+z^2])/ψ[x,y,z,α1,α2,c2]//Simplify];
ϵloc[x,y,z,α1,α2,c2]

(* par1 = α1, par2 = α2, par3 = c2 *)
minFunc[par_] :=
  Module[{n,dim,xmin,xmax,a,b,ρ,x0,δ,points,ϵ,energy},
n=100000;
dim=3;
xmax=8.0;xmin=-xmax;
a={xmin,xmin,xmin};b={xmax,xmax,xmax};
ρ[x_]:=ψ[x[[1]],x[[2]],x[[3]],par[[1]],par[[2]],par[[3]]]^2;
x0=b/10.0;δ=b/10;
points=MyRandomMetropolisND[n,dim,ρ,a,b,x0,δ,1];
ϵ[x_]:=ϵloc[x[[1]],x[[2]],x[[3]],par[[1]],par[[2]],par[[3]]];
energy=MyMCIntegration[points,ϵ];
Print["Parameters = ", par, ", energy = ", energy[[1]], ", probable error = ", energy[[2]]];
Return[energy[[1]]];
];

minFunc[{0.25,1.50,1.00}]

MyMinimization[{0.5,2.00,0.50}, {0.1,0.1,0.1}, minFunc, 0.01]

f1[r_]:=Exp[-25/100r^2]+Exp[-150/100r^2];f2[r_]:=Exp[-0.25r^2]+0.55Exp[-2.2r^2];
norm=1/Sqrt[Integrate[Exp[-2r],{r,0,Infinity}]];
norm1=1/Sqrt[Integrate[f1[r]^2,{r,0,Infinity}]];
norm2=1/Sqrt[Integrate[f2[r]^2,{r,0,Infinity}]];
Plot[{Exp[-r]*norm,f1[r]*norm1,f2[r]*norm2},{r,0,4},PlotRange->All]