Merge pull request #4 from Panadestein/solutions_c

Solutions in C

.gitignore

@@ -0,0 +1,8 @@
+__pycache__
+plt_dat
+data
+get_energ
+get_var
+qmc_uni
+qmc_met
+vmc

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;
+}

hydrogen.h

@@ -0,0 +1,64 @@
+#include <math.h>
+#include <stdio.h>
+
+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;
+      }
+}

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;
+}

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;
+}

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;
+}
+

qmc_metropolis.c

@@ -0,0 +1,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;
+}

qmc_stats.h

@@ -0,0 +1,59 @@
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+
+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]);
+      }
+}

qmc_uniform.c

@@ -0,0 +1,44 @@
+#include "hydrogen.h"
+#include "qmc_stats.h"
+#include <stdlib.h>
+#include <time.h>
+
+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;
+}
+

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;
+}

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;
+}