TREX/qmc-lttc
1
0
Fork 0
mirror of https://github.com/TREX-CoE/qmc-lttc.git synced 2024-12-22 04:15:01 +01:00
Code Issues Projects Releases Wiki Activity

Fixed VMC

Browse Source
This commit is contained in:
Panadestein 2021-04-27 20:06:31 +02:00
parent 0bbfb469c9
commit c55a1bb1f1
3 changed files with 3000029 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3000002
pp2 Normal file
View File

File diff suppressed because it is too large Load Diff

10
qmc_stats.h
View File

@ -33,8 +33,6 @@ void random_gauss(double *z, const int n) {
const double two_pi = 2.0 * acos(-1.0); const double two_pi = 2.0 * acos(-1.0);
double u[n + 1]; double u[n + 1];
srand(time(NULL));
for (int i = 0; i < n + 1; ++i) { for (int i = 0; i < n + 1; ++i) {
u[i] = (double) rand() / RAND_MAX; u[i] = (double) rand() / RAND_MAX;
} }
@ -44,18 +42,18 @@ void random_gauss(double *z, const int n) {
// n is even // n is even
for (int i = 0; i < n; i+=2) { for (int i = 0; i < n; i+=2) {
z[i] = sqrt(-2.0 * log(u[i])); z[i] = sqrt(-2.0 * log(u[i]));
z[i] = z[i] * cos(two_pi * u[i + 1]);
z[i + 1] = z[i] * sin(two_pi * u[i + 1]); z[i + 1] = z[i] * sin(two_pi * u[i + 1]);
z[i] = z[i] * cos(two_pi * u[i + 1]);
} }
} }
else { else {
// n is odd // n is odd
for (int i = 0; i < n - 1; i+=2) { for (int i = 0; i < n - 1; i+=2) {
z[i] = sqrt(-2.0 * log(u[i])); z[i] = sqrt(-2.0 * log(u[i]));
z[i] = z[i] * cos(two_pi * u[i + 1]);
z[i + 1] = z[i] * sin(two_pi * u[i + 1]); z[i + 1] = z[i] * sin(two_pi * u[i + 1]);
z[i] = z[i] * cos(two_pi * u[i + 1]);
} }
z[n] = sqrt(-2.0 * log(u[n])); z[n - 1] = sqrt(-2.0 * log(u[n - 1]));
z[n] = z[n] * cos(two_pi * u[n + 1]); z[n - 1] = z[n - 1] * cos(two_pi * u[n]);
} }
} }

40
vmc_metropolis.c
View File

@ -6,17 +6,21 @@ void variational_montecarlo(double a, const int nmax, double dt,
int n_accept; int n_accept;
double r_old[3], r_new[3], psi_old, psi_new; double r_old[3], r_new[3], psi_old, psi_new;
double d_old[3], d_new[3], d2_old, d2_new; double d_old[3], d_new[3], d2_old, d2_new;
double rnd[3], aval, fact_a, fact_b, fact_exp, u; double rnd[3], fact_a, fact_b, sq_dt, q, u;
sq_dt = sqrt(dt);
// Initial position // Initial position
random_gauss(r_old, 3); random_gauss(r_old, 3);
psi_old = psi(a, r_old, 3) * psi(a, r_old, 3);
drift(a, r_old, d_old, 3); drift(a, r_old, d_old, 3);
d2_old = 0.0; d2_old = 0.0;
for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) {
d2_old += d_old[i] * d_old[i]; d2_old += d_old[j] * d_old[j];
} }
psi_old = psi(a, r_old, 3);
*energy = 0.0; *energy = 0.0;
*accept = 0.0; *accept = 0.0;
n_accept = 0; n_accept = 0;
@ -24,39 +28,41 @@ void variational_montecarlo(double a, const int nmax, double dt,
// Compute and accumulate the local energy // Compute and accumulate the local energy
*energy += e_loc(a, r_old, 3); *energy += e_loc(a, r_old, 3);
// Compute new position // Compute new position (correct variance of sampled Gaussian)
random_gauss(rnd, 3); random_gauss(rnd, 3);
for (int j = 0; j < 3; ++j) { for (int j = 0; j < 3; ++j) {
r_new[j] = r_old[j] + dt * d_old[j] + rnd[j]; r_new[j] = r_old[j] + dt * d_old[j] + rnd[j] * sq_dt;
} }
// New WF and acceptance probability // New WF and acceptance probability
psi_new = psi(a, r_new, 3) * psi(a, r_new, 3);
drift(a, r_new, d_new, 3); drift(a, r_new, d_new, 3);
d2_new = 0.0; d2_new = 0.0;
for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) {
d2_new += d_new[i] * d_new[i]; d2_new += d_new[j] * d_new[j];
} }
psi_new = psi(a, r_new, 3);
// Compute the ratio of probabilities q // Compute the ratio of probabilities q
fact_a = 0.5 * dt * (d2_new - d2_old);
fact_b = 0.0; fact_b = 0.0;
for (int i = 0; i < 3; ++i) { for (int j = 0; j < 3; ++j) {
fact_b += (d_new[i] + d_old[i]) * (r_new[i] - r_old[i]); fact_b += (d_new[j] + d_old[j]) * (r_new[j] - r_old[j]);
} }
fact_exp = exp(fact_b - fact_a); fact_a = 0.5 * dt * (d2_new - d2_old) + fact_b;
aval = fact_exp * psi_new / psi_old;
q = psi_new / psi_old;
q = exp(-fact_a) * q * q;
u = (double) rand() / RAND_MAX; u = (double) rand() / RAND_MAX;
if (u <= aval) { if (u <= q) {
n_accept += 1;
for (int j = 0; j < 3; ++j) { for (int j = 0; j < 3; ++j) {
r_old[j] = r_new[j]; r_old[j] = r_new[j];
d_old[j] = d_new[j]; d_old[j] = d_new[j];
} }
d2_old = d2_new;
psi_old = psi_new; psi_old = psi_new;
n_accept += 1;
} }
} }
*energy /= nmax; *energy /= nmax;