1
0
mirror of https://github.com/TREX-CoE/qmc-lttc.git synced 2024-10-05 15:56:17 +02:00
qmc-lttc/qmc_metropolis.py

45 lines
858 B
Python
Raw Normal View History

#!/usr/bin/env python3
2021-01-20 18:31:49 +01:00
from hydrogen import *
from qmc_stats import *
2021-01-29 13:23:00 +01:00
def MonteCarlo(a,nmax,dt):
energy = 0.
2021-01-26 10:02:38 +01:00
N_accep = 0
2021-01-29 13:23:00 +01:00
r_old = np.random.uniform(-dt, dt, (3))
2021-01-20 18:31:49 +01:00
psi_old = psi(a,r_old)
2021-01-29 13:23:00 +01:00
2021-01-20 18:31:49 +01:00
for istep in range(nmax):
2021-01-29 13:23:00 +01:00
energy += e_loc(a,r_old)
r_new = r_old + np.random.uniform(-dt,dt,(3))
2021-01-20 18:31:49 +01:00
psi_new = psi(a,r_new)
2021-01-29 13:23:00 +01:00
2021-01-20 18:31:49 +01:00
ratio = (psi_new / psi_old)**2
2021-01-29 13:23:00 +01:00
if np.random.uniform() <= ratio:
2021-01-26 10:02:38 +01:00
N_accep += 1
2021-01-29 13:23:00 +01:00
r_old = r_new
2021-01-20 18:31:49 +01:00
psi_old = psi_new
2021-01-29 13:23:00 +01:00
2021-01-26 10:02:38 +01:00
return energy/nmax, N_accep/nmax
2021-01-20 18:31:49 +01:00
2021-01-26 10:02:38 +01:00
# Run simulation
a = 1.2
2021-01-20 18:31:49 +01:00
nmax = 100000
dt = 1.0
2021-01-29 13:23:00 +01:00
X0 = [ MonteCarlo(a,nmax,dt) for i in range(30)]
2021-01-26 10:02:38 +01:00
# Energy
X = [ x for (x, _) in X0 ]
2021-01-20 18:31:49 +01:00
E, deltaE = ave_error(X)
print(f"E = {E} +/- {deltaE}")
2021-01-26 10:02:38 +01:00
# Acceptance rate
X = [ x for (_, x) in X0 ]
A, deltaA = ave_error(X)
2021-01-20 18:31:49 +01:00
print(f"A = {A} +/- {deltaA}")