diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e7ebbb4 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +__pycache__ +plt_dat +data +get_energ +get_var +qmc_uni +qmc_met +vmc diff --git a/energy_hydrogen.c b/energy_hydrogen.c new file mode 100644 index 0000000..2bb9ec3 --- /dev/null +++ b/energy_hydrogen.c @@ -0,0 +1,35 @@ +#include "hydrogen.h" + +int main() { + const double a[6] = {0.1, 0.2, 0.5, 1., 1.5, 2.}; + const int sizex = 50; + double energy, r[3]; + double x[sizex], dx, w, wsum, dV; + + dx = 10.0 / (sizex - 1); + for (int i = 0; i < sizex; ++i) { + x[i] = -5.0 + (i - 1) * dx; + } + dV = dx * dx * dx; + + for (int l = 0; l < 6; ++l) { + energy = 0.0; + wsum = 0.0; + for (int i = 0; i < sizex; ++i) { + r[0] = x[i]; + for (int j = 0; j < sizex; ++j) { + r[1] = x[j]; + for (int k = 0; k < sizex; ++k) { + r[2] = x[k]; + w = psi(a[l], r, 3) * psi(a[l], r, 3) * dV; + wsum += w; + energy += w * e_loc(a[l], r, 3); + } + } + } + printf("a = %.2f E = %.5f \n", a[l], energy / wsum); + } + + + return 0; +} diff --git a/hydrogen.h b/hydrogen.h new file mode 100644 index 0000000..40da3b1 --- /dev/null +++ b/hydrogen.h @@ -0,0 +1,64 @@ +#include +#include + +double potential(double *r, const int l) { + double pot; + + pot = 0.0; + for (int i = 0; i < l; ++i) { + pot += r[i] * r[i]; + } + + if (pot > 0.0) { + return -1.0 / sqrt(pot); + } + return 1e6; + +} + +double psi(double a, double *r, const int l) { + double psival, rnorm; + + psival = 0.0; + rnorm = 0.0; + for (int i = 0; i < l; ++i) { + rnorm += r[i]*r[i]; + } + rnorm = sqrt(rnorm); + + psival = exp(-a * rnorm); + + return psival; +} + +double kinetic(double a, double *r, const int l) { + double rnorm; + + rnorm = 0.0; + for (int i = 0; i < l; ++i) { + rnorm += r[i]*r[i]; + } + rnorm = sqrt(rnorm); + + return a * (1 / rnorm - 0.5 * a); + +} + +double e_loc(double a, double *r, const int l) { + return kinetic(a, r, l) + potential(r, l); +} + +void drift(double a, double *r, double *d, const int l) { + double rnorm, fact; + + rnorm = 0.0; + for (int i = 0; i < l; ++i) { + rnorm += r[i]*r[i]; + } + rnorm = sqrt(rnorm); + + fact = -a / rnorm; + for (int i = 0; i < l; ++i) { + d[i] = r[i] * fact; + } +} diff --git a/pdmc.c b/pdmc.c new file mode 100644 index 0000000..87a8f76 --- /dev/null +++ b/pdmc.c @@ -0,0 +1,111 @@ +#include "hydrogen.h" +#include "qmc_stats.h" + +void pdmc(double a, const int nmax, double dt, + double *energy, double *accept, double tau, double e_ref) { + int n_accept; + double r_old[3], r_new[3], psi_old, psi_new; + double d_old[3], d_new[3], d2_old, d2_new; + double rnd[3], fact_a, fact_b, sq_dt, q, u; + double e, w, tau_current, norma; + + sq_dt = sqrt(dt); + + // Initial position + random_gauss(r_old, 3); + + drift(a, r_old, d_old, 3); + d2_old = 0.0; + for (int j = 0; j < 3; ++j) { + d2_old += d_old[j] * d_old[j]; + } + + psi_old = psi(a, r_old, 3); + + *energy = 0.0; + *accept = 0.0; + n_accept = 0; + w = 1.0; + tau_current = 0.0; + norma = 0.0; + for (int i = 0; i < nmax; ++i) { + // Compute local energy and weight + e = e_loc(a, r_old, 3); + w = w * exp(-dt * (e - e_ref)); + + // Accumulate local energy and norm + norma = norma + w; + *energy += w * e; + + // Increase time and reset if threshold reached + tau_current += dt; + if (tau_current > tau) { + w = 1.0; + tau_current = 0.0; + } + + // Compute new position (correct variance of sampled Gaussian) + random_gauss(rnd, 3); + for (int j = 0; j < 3; ++j) { + r_new[j] = r_old[j] + dt * d_old[j] + rnd[j] * sq_dt; + } + + // New WF and acceptance probability + drift(a, r_new, d_new, 3); + d2_new = 0.0; + for (int j = 0; j < 3; ++j) { + d2_new += d_new[j] * d_new[j]; + } + + psi_new = psi(a, r_new, 3); + + // Compute the ratio of probabilities q + fact_b = 0.0; + for (int j = 0; j < 3; ++j) { + fact_b += (d_new[j] + d_old[j]) * (r_new[j] - r_old[j]); + } + fact_a = 0.5 * dt * (d2_new - d2_old) + fact_b; + + q = psi_new / psi_old; + q = exp(-fact_a) * q * q; + + u = (double) rand() / RAND_MAX; + + if (u <= q) { + n_accept += 1; + for (int j = 0; j < 3; ++j) { + r_old[j] = r_new[j]; + d_old[j] = d_new[j]; + } + d2_old = d2_new; + psi_old = psi_new; + } + } + *energy /= norma; + *accept = (double) n_accept / nmax; +} + +int main() { + const double a = 1.2; + const double dt = 0.05; + const double tau = 100.0; + const double e_ref = -0.5; + const long nmax = 1e5; + int nruns = 30; + + srand(time(NULL)); + + double ene[nruns], acc[nruns], obs[2]; + + for (int i = 0; i < nruns; ++i) { + pdmc(a, nmax, dt, &ene[i], &acc[i], tau, e_ref); + } + + ave_error(ene, nruns, obs); + printf("E = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + ave_error(acc, nruns, obs); + printf("A = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + return 0; +} diff --git a/plot_hydrogen.c b/plot_hydrogen.c new file mode 100644 index 0000000..074b8bc --- /dev/null +++ b/plot_hydrogen.c @@ -0,0 +1,27 @@ +#include "hydrogen.h" + +int main() { + const double a[6] = {0.1, 0.2, 0.5, 1., 1.5, 2.}; + const int sizex = 50; + double x[sizex], dx; + + dx = 10.0 / (sizex - 1); + for (int i = 0; i < sizex; ++i) { + x[i] = -5.0 + (i - 1) * dx; + } + + FILE * fil; + + fil = fopen("./data", "w+"); + + for (int i; i < 6; ++i) { + for (int j = 0; j < sizex; ++j) { + fprintf(fil, "%lf %lf\n", x[j], e_loc(a[i], &x[j], 1)); + } + fprintf(fil, "\n\n"); + } + + fclose(fil); + + return 0; +} diff --git a/qmc_gaussian.c b/qmc_gaussian.c new file mode 100644 index 0000000..69e4936 --- /dev/null +++ b/qmc_gaussian.c @@ -0,0 +1,44 @@ +#include "qmc_stats.h" +#include "hydrogen.h" + +double gaussian(double *r) { + const double norm_gauss = 1.0 / pow(2.0 * acos(-1.0), 1.5); + return norm_gauss * exp(-0.5 * (r[0]*r[0] + r[1]*r[1] + r[2]*r[2])); +} + +double gaussian_montecarlo(double a, const int nmax) { + double r[3], energy, norm, w; + + norm = 0.0; + energy = 0.0; + for (int i = 0; i < nmax; ++i) { + random_gauss(r, 3); + w = psi(a, r, 3) * psi(a, r, 3) / gaussian(r); + norm += w; + energy += w * e_loc(a, r, 3); + } + + return energy / norm; + +} + +int main() { + const double a = 1.2; + const long nmax = 1e5; + int nruns = 30; + + srand(time(NULL)); + + double x[nruns], obs[2]; + + for (int i = 0; i < nruns; ++i) { + x[i] = gaussian_montecarlo(a, nmax); + } + + ave_error(x, nruns, obs); + + printf("E = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + return 0; +} + diff --git a/qmc_metropolis.c b/qmc_metropolis.c new file mode 100644 index 0000000..ecfd997 --- /dev/null +++ b/qmc_metropolis.c @@ -0,0 +1,73 @@ +#include "hydrogen.h" +#include "qmc_stats.h" +#include +#include + +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; +} diff --git a/qmc_stats.h b/qmc_stats.h new file mode 100644 index 0000000..0feb898 --- /dev/null +++ b/qmc_stats.h @@ -0,0 +1,59 @@ +#include +#include +#include +#include + +void ave_error(double *x, const int n, double obs[2]) { + // obs[1] --> mean and obs[2] --> err + double var; + + if (n == 1) { + obs[0] = x[0]; + obs[1] = 0.0; + } + else { + for (int i = 0; i < n; ++i) { + obs[0] += x[i]; + } + obs[0] = obs[0] / n; + + var = 0.0; + for (int i = 0; i < n; ++i) { + var += (x[i] - obs[0]) * (x[i] - obs[0]); + } + var = var / (n - 1); + + obs[1] = sqrt(var / n); + } + +} + +void random_gauss(double *z, const int n) { + // Box Muller method + const double two_pi = 2.0 * acos(-1.0); + double u[n + 1]; + + for (int i = 0; i < n + 1; ++i) { + u[i] = (double) rand() / RAND_MAX; + } + + + if (!(n & 1)) { + // n is even + for (int i = 0; i < n; i+=2) { + z[i] = sqrt(-2.0 * log(u[i])); + z[i + 1] = z[i] * sin(two_pi * u[i + 1]); + z[i] = z[i] * cos(two_pi * u[i + 1]); + } + } + else { + // n is odd + for (int i = 0; i < n - 1; i+=2) { + z[i] = sqrt(-2.0 * log(u[i])); + z[i + 1] = z[i] * sin(two_pi * u[i + 1]); + z[i] = z[i] * cos(two_pi * u[i + 1]); + } + z[n - 1] = sqrt(-2.0 * log(u[n - 1])); + z[n - 1] = z[n - 1] * cos(two_pi * u[n]); + } +} diff --git a/qmc_uniform.c b/qmc_uniform.c new file mode 100644 index 0000000..ebdb891 --- /dev/null +++ b/qmc_uniform.c @@ -0,0 +1,44 @@ +#include "hydrogen.h" +#include "qmc_stats.h" +#include +#include + +double uniform_montecarlo(double a, const int nmax) { + double r[3], energy, norm, w, rnd; + + norm = 0.0; + energy = 0.0; + for (int i = 0; i < nmax; ++i) { + for (int j = 0; j < 3; ++j) { + rnd = (double) rand() / RAND_MAX; + r[j] = -5.0 + 10.0 * rnd; + } + w = psi(a, r, 3) * psi(a, r, 3); + norm += w; + energy += w * e_loc(a, r, 3); + } + + return energy / norm; + +} + +int main() { + const double a = 1.2; + const long nmax = 1e5; + int nruns = 30; + + srand(time(NULL)); + + double x[nruns], obs[2]; + + for (int i = 0; i < nruns; ++i) { + x[i] = uniform_montecarlo(a, nmax); + } + + ave_error(x, nruns, obs); + + printf("E = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + return 0; +} + diff --git a/variance_hydrogen.c b/variance_hydrogen.c new file mode 100644 index 0000000..a02c3e0 --- /dev/null +++ b/variance_hydrogen.c @@ -0,0 +1,41 @@ +#include "hydrogen.h" + +int main() { + const double a[6] = {0.1, 0.2, 0.5, 1., 1.5, 2.}; + const int sizex = 50; + double energy, r[3], sigma; + double x[sizex], dx, w, wsum, dV, e_tmp; + + dx = 10.0 / (sizex - 1); + for (int i = 0; i < sizex; ++i) { + x[i] = -5.0 + i * dx; + } + dV = dx * dx *dx; + + for (int l = 0; l < 6; ++l) { + energy = 0.0; + sigma = 0.0; + wsum = 0.0; + for (int i = 0; i < sizex; ++i) { + r[0] = x[i]; + for (int j = 0; j < sizex; ++j) { + r[1] = x[j]; + for (int k = 0; k < sizex; ++k) { + r[2] = x[k]; + w = psi(a[l], r, 3) * psi(a[l], r, 3) * dV; + wsum += w; + e_tmp = e_loc(a[l], r, 3); + energy += w * e_tmp; + sigma += w * e_tmp * e_tmp; + } + } + } + energy = energy / wsum; + sigma = sigma / wsum; + sigma -= energy * energy; + printf("a = %.2lf E = %.5lf Sigma^2 = %.5lf\n", a[l], energy, sigma); + } + + + return 0; +} diff --git a/vmc_metropolis.c b/vmc_metropolis.c new file mode 100644 index 0000000..651b71c --- /dev/null +++ b/vmc_metropolis.c @@ -0,0 +1,93 @@ +#include "hydrogen.h" +#include "qmc_stats.h" + +void variational_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; + double d_old[3], d_new[3], d2_old, d2_new; + double rnd[3], fact_a, fact_b, sq_dt, q, u; + + sq_dt = sqrt(dt); + + // Initial position + random_gauss(r_old, 3); + + drift(a, r_old, d_old, 3); + d2_old = 0.0; + for (int j = 0; j < 3; ++j) { + d2_old += d_old[j] * d_old[j]; + } + + psi_old = 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 (correct variance of sampled Gaussian) + random_gauss(rnd, 3); + for (int j = 0; j < 3; ++j) { + r_new[j] = r_old[j] + dt * d_old[j] + rnd[j] * sq_dt; + } + + // New WF and acceptance probability + drift(a, r_new, d_new, 3); + d2_new = 0.0; + for (int j = 0; j < 3; ++j) { + d2_new += d_new[j] * d_new[j]; + } + + psi_new = psi(a, r_new, 3); + + // Compute the ratio of probabilities q + fact_b = 0.0; + for (int j = 0; j < 3; ++j) { + fact_b += (d_new[j] + d_old[j]) * (r_new[j] - r_old[j]); + } + fact_a = 0.5 * dt * (d2_new - d2_old) + fact_b; + + q = psi_new / psi_old; + q = exp(-fact_a) * q * q; + + u = (double) rand() / RAND_MAX; + + if (u <= q) { + n_accept += 1; + for (int j = 0; j < 3; ++j) { + r_old[j] = r_new[j]; + d_old[j] = d_new[j]; + } + d2_old = d2_new; + psi_old = psi_new; + } + } + *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 ene[nruns], acc[nruns], obs[2]; + + for (int i = 0; i < nruns; ++i) { + variational_montecarlo(a, nmax, dt, &ene[i], &acc[i]); + } + + ave_error(ene, nruns, obs); + printf("E = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + ave_error(acc, nruns, obs); + printf("A = %.5lf +/- %.5lf\n", obs[0], obs[1]); + + return 0; +}