1
0
mirror of https://github.com/TREX-CoE/qmc-lttc.git synced 2024-07-16 08:00:45 +02:00
qmc-lttc/vmc_metropolis.py
Ravindra Shinde 65efa3a248
variable was incorrectly used
Runtime error   ::   NameError: name 'accep_rate' is not defined
2021-02-02 15:30:11 +01:00

61 lines
1.3 KiB
Python

#!/usr/bin/env python3
from hydrogen import *
from qmc_stats import *
def MonteCarlo(a,nmax,dt):
sq_dt = np.sqrt(dt)
energy = 0.
N_accep = 0
r_old = np.random.normal(loc=0., scale=1.0, size=(3))
d_old = drift(a,r_old)
d2_old = np.dot(d_old,d_old)
psi_old = psi(a,r_old)
for istep in range(nmax):
chi = np.random.normal(loc=0., scale=1.0, size=(3))
energy += e_loc(a,r_old)
r_new = r_old + dt * d_old + sq_dt * chi
d_new = drift(a,r_new)
d2_new = np.dot(d_new,d_new)
psi_new = psi(a,r_new)
# Metropolis
prod = np.dot((d_new + d_old), (r_new - r_old))
argexpo = 0.5 * (d2_new - d2_old)*dt + prod
q = psi_new / psi_old
q = np.exp(-argexpo) * q*q
if np.random.uniform() <= q:
N_accep += 1
r_old = r_new
d_old = d_new
d2_old = d2_new
psi_old = psi_new
return energy/nmax, N_accep/nmax
# Run simulation
a = 1.2
nmax = 100000
dt = 1.0
X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
# Energy
X = [ x for (x, _) in X0 ]
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
print(f"A = {A} +/- {deltaA}")