1
0
mirror of https://github.com/TREX-CoE/qmc-lttc.git synced 2024-12-11 23:14:09 +01:00

Solutions in C up to 3.4

This commit is contained in:
Panadestein 2021-04-25 20:02:47 +02:00
parent b15ea4b02f
commit 6eb6eed6db
8 changed files with 230 additions and 3 deletions

7
.gitignore vendored
View File

@ -1,2 +1,7 @@
__pycache__
plt_dat
data
data
get_energ
get_var
qmc_uni
qmc_met

35
energy_hydrogen.c Normal file
View File

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

View File

@ -12,7 +12,7 @@ double potential(double *r, const int l) {
if (pot > 0.0) {
return -1.0 / sqrt(pot);
}
return -1;
return 1e6;
}

73
qmc_metropolis.c Normal file
View File

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

27
qmc_stats.h Normal file
View File

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

44
qmc_uniform.c Normal file
View File

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

41
variance_hydrogen.c Normal file
View File

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

View File

@ -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}")