diff --git a/.gitignore b/.gitignore
index c912f89..e7ebbb4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,3 +5,4 @@ get_energ
 get_var
 qmc_uni
 qmc_met
+vmc
diff --git a/hydrogen.h b/hydrogen.h
index 905b482..c061238 100644
--- a/hydrogen.h
+++ b/hydrogen.h
@@ -46,3 +46,17 @@ double kinetic(double a, double *r, const int l) {
 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;
+
+      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/pp b/pp
new file mode 100755
index 0000000..30dfb73
Binary files /dev/null and b/pp differ
diff --git a/qmc_gaussian.c b/qmc_gaussian.c
new file mode 100644
index 0000000..1d537df
--- /dev/null
+++ b/qmc_gaussian.c
@@ -0,0 +1,7 @@
+#include <math.h>
+#include <stdio.h>
+
+double gaussian(double *r) {
+      // Pending
+      return 0;
+}
diff --git a/qmc_stats.h b/qmc_stats.h
index 3ac432d..0bd615f 100644
--- a/qmc_stats.h
+++ b/qmc_stats.h
@@ -1,5 +1,7 @@
 #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
@@ -25,3 +27,35 @@ void ave_error(double *x, const int n, double obs[2]) {
       }
 
 }
+
+void random_gauss(double *z, const int n) {
+      // Box Muller method
+      const double two_pi = 2.0 * acos(-1.0);
+      double u[n + 1];
+
+      srand(time(NULL));
+
+      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] = z[i] * cos(two_pi * u[i + 1]);
+		  z[i + 1] = z[i] * sin(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] = z[i] * cos(two_pi * u[i + 1]);
+		  z[i + 1] = z[i] * sin(two_pi * u[i + 1]);
+	    }
+	    z[n] = sqrt(-2.0 * log(u[n]));
+	    z[n] = z[n] * cos(two_pi * u[n + 1]);
+      }
+}
diff --git a/tt.c b/tt.c
new file mode 100644
index 0000000..20f4d12
--- /dev/null
+++ b/tt.c
@@ -0,0 +1,11 @@
+#include "qmc_stats.h"
+
+int main() {
+      double rnd[3];
+      random_gauss(rnd, 3);
+      for (int i = 0; i < 3; ++i) {
+	    printf("val %lf\n", rnd[i]);
+      }
+
+      return 0;
+}
diff --git a/vmc_metropolis.c b/vmc_metropolis.c
new file mode 100644
index 0000000..2175e34
--- /dev/null
+++ b/vmc_metropolis.c
@@ -0,0 +1,87 @@
+#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], aval, fact_a, fact_b, fact_exp, u;
+
+      // Initial position
+      random_gauss(r_old, 3);
+      psi_old = psi(a, r_old, 3) * psi(a, r_old, 3);
+      drift(a, r_old, d_old, 3);
+      d2_old = 0.0;
+      for (int i = 0; i < 3; ++i) {
+	    d2_old += d_old[i] * d_old[i];
+      }
+
+      *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
+	    random_gauss(rnd, 3);
+	    for (int j = 0; j < 3; ++j) {
+		  r_new[j] = r_old[j] + dt * d_old[j] + rnd[j];
+	    }
+
+	    // New WF and acceptance probability
+	    psi_new = psi(a, r_new, 3) * psi(a, r_new, 3);
+
+            drift(a, r_new, d_new, 3);
+	    d2_new = 0.0;
+            for (int i = 0; i < 3; ++i) {
+            	    d2_new += d_new[i] * d_new[i];
+            }
+
+	    // Compute the ratio of probabilities q
+	    fact_a = 0.5 * dt * (d2_new - d2_old);
+	    fact_b = 0.0;
+	    for (int i = 0; i < 3; ++i) {
+		  fact_b += (d_new[i] + d_old[i]) * (r_new[i] - r_old[i]);
+	    }
+	    fact_exp = exp(fact_b - fact_a);
+	    aval = fact_exp * 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];
+			d_old[j] = d_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 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;
+}