diff --git a/.gitignore b/.gitignore index d32491f..c912f89 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,7 @@ +__pycache__ plt_dat -data \ No newline at end of file +data +get_energ +get_var +qmc_uni +qmc_met 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 index c73632f..905b482 100644 --- a/hydrogen.h +++ b/hydrogen.h @@ -12,7 +12,7 @@ double potential(double *r, const int l) { if (pot > 0.0) { return -1.0 / sqrt(pot); } - return -1; + return 1e6; } 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..3ac432d --- /dev/null +++ b/qmc_stats.h @@ -0,0 +1,27 @@ +#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); + } + +} 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/variance_hydrogen.py b/variance_hydrogen.py index 07b7d17..0892f79 100644 --- a/variance_hydrogen.py +++ b/variance_hydrogen.py @@ -4,6 +4,8 @@ import numpy as np from hydrogen import e_loc, psi interval = np.linspace(-5,5,num=50) +print(interval + ) delta = (interval[1]-interval[0])**3 @@ -34,5 +36,5 @@ for a in [0.1, 0.2, 0.5, 0.9, 1., 1.5, 2.]: E = E / norm E2 = E2 / norm - s2 = E2 - E**2 + s2 = E2 #- E**2 print(f"a = {a} \t E = {E:10.8f} \t \sigma^2 = {s2:10.8f}")