-
Notifications
You must be signed in to change notification settings - Fork 4
/
qmc_metropolis.c
73 lines (57 loc) · 1.67 KB
/
qmc_metropolis.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
#include "hydrogen.h"
#include "qmc_stats.h"
#include <stdlib.h>
#include <time.h>
void metropolis_montecarlo(double a, const int nmax, double dt,
double *energy, double *accept) {
int n_accept;
double r_old[3], r_new[3], psi_old, psi_new, rnd;
double aval, u;
// Initial position
for (int j = 0; j < 3; ++j) {
rnd = (double) rand() / RAND_MAX;
r_old[j] = dt * (2.0 * rnd - 1.0);
}
psi_old = psi(a, r_old, 3) * psi(a, r_old, 3);
*energy = 0.0;
*accept = 0.0;
n_accept = 0;
for (int i = 0; i < nmax; ++i) {
// Compute and accumulate the local energy
*energy += e_loc(a, r_old, 3);
// Compute new position
for (int j = 0; j < 3; ++j) {
rnd = (double) rand() / RAND_MAX;
r_new[j] = r_old[j] + dt * (2.0 * rnd - 1.0);
}
// New WF and acceptance probability
psi_new = psi(a, r_new, 3) * psi(a, r_new, 3);
aval = psi_new / psi_old;
u = (double) rand() / RAND_MAX;
if (u <= aval) {
for (int j = 0; j < 3; ++j) {
r_old[j] = r_new[j];
}
psi_old = psi_new;
n_accept += 1;
}
}
*energy /= nmax;
*accept = (double) n_accept / nmax;
}
int main() {
const double a = 1.2;
const double dt = 1.0;
const long nmax = 1e5;
int nruns = 30;
srand(time(NULL));
double x[nruns], y[nruns], obs[2];
for (int i = 0; i < nruns; ++i) {
metropolis_montecarlo(a, nmax, dt, &x[i], &y[i]);
}
ave_error(x, nruns, obs);
printf("E = %.5lf +/- %.5lf\n", obs[0], obs[1]);
ave_error(y, nruns, obs);
printf("A = %.5lf +/- %.5lf\n", obs[0], obs[1]);
return 0;
}