Adding split_op.c and energy.c

contents/quantum_systems/code/c/energy.c

@@ -0,0 +1,60 @@
+#include <complex.h>
+#include <string.h>
+#include <math.h>
+
+#include <fftw3.h>
+
+void fft(double complex *x, int n, bool inverse) {
+ double complex y[n];
+ memset(y, 0, sizeof(y));
+ fftw_plan p;
+
+ if (inverse) {
+ p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
+ FFTW_BACKWARD, FFTW_ESTIMATE);
+ } else {
+ p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
+ FFTW_FORWARD, FFTW_ESTIMATE);
+ }
+
+ fftw_execute(p);
+ fftw_destroy_plan(p);
+
+ for (size_t i = 0; i < n; ++i) {
+ x[i] = y[i] / sqrt((double)n);
+ }
+}
+
+double calculate_energy(double complex *wfc, double complex *h_r,
+ double complex *h_k, double dx, size_t size) {
+ double complex wfc_k[size];
+ double complex wfc_c[size];
+ memcpy(wfc_k, wfc, sizeof(wfc_k));
+ fft(wfc_k, size, false);
+
+ for (size_t i = 0; i < size; ++i) {
+ wfc_c[i] = conj(wfc[i]);
+ }
+
+ double complex energy_k[size];
+ double complex energy_r[size];
+
+ for (size_t i = 0; i < size; ++i) {
+ energy_k[i] = wfc_k[i] * h_k[i];
+ }
+
+ fft(energy_k, size, true);
+
+ for (size_t i = 0; i < size; ++i) {
+ energy_k[i] *= wfc_c[i];
+ energy_r[i] = wfc_c[i] * H_r[i] * wfc[i];
+ }
+
+ double energy_final = 0;
+
+ for (size_t i = 0; i < size; ++i) {
+ energy_final += creal(energy_k[i] + energy_r[i]);
+ }
+
+ return energy_final * dx;
+}

contents/quantum_systems/quantum_systems.md

@@ -227,6 +227,8 @@ This ultimately looks like this:
 {% method %}
 {% sample lang="jl" %}
 [import, lang:"julia"](code/julia/energy.jl)
+{% sample lang="c" %}
+[import:28-, lang:"c_cpp"](code/julia/energy.jl)
 {% endmethod %}
 
 This calculation will be used in many different simulations of quantum systems to check our results.

contents/split-operator_method/code/c/split_op.c

@@ -0,0 +1,210 @@
+#include <complex.h>
+#include <stdbool.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include <fftw3.h>
+
+struct params {
+ double xmax;
+ int res;
+ double dt;
+ int timesteps;
+ double dx;
+ double *x;
+ double dk;
+ double *k;
+ bool im_time;
+};
+
+struct operators {
+ size_t size;
+ double complex *v;
+ double complex *pe;
+ double complex *ke;
+ double complex *wfc;
+};
+
+void fft(double complex *x, int n, bool inverse) {
+ double complex y[n];
+ memset(y, 0, sizeof(y));
+ fftw_plan p;
+
+ if (inverse) {
+ p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
+ FFTW_BACKWARD, FFTW_ESTIMATE);
+ } else {
+ p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
+ FFTW_FORWARD, FFTW_ESTIMATE);
+ }
+
+ fftw_execute(p);
+ fftw_destroy_plan(p);
+
+ for (size_t i = 0; i < n; ++i) {
+ x[i] = y[i] / sqrt((double)n);
+ }
+}
+
+void init_params(struct params *par, double xmax, int res, double dt,
+ int timesteps, bool im) {
+
+ par->xmax = xmax;
+ par->res = res;
+ par->dt = dt;
+ par->timesteps = timesteps;
+ par->dx = 2.0 * xmax / res;
+ par->x = malloc(sizeof(double) * res);
+ par->dk = M_PI / xmax;
+ par->k = malloc(sizeof(double) * res);
+ par->im_time = im;
+
+ for (size_t i = 0; i < res; ++i) {
+ par->x[i] = xmax / res - xmax + i * (2.0 * xmax / res);
+ if (i < res / 2) {
+ par->k[i] = i * M_PI / xmax;
+ } else {
+ par->k[i] = ((double)i - res) * M_PI / xmax;
+ }
+ }
+}
+
+void init_operators(struct operators *opr, struct params par, double voffset,
+ double wfcoffset) {
+
+ opr->size = par.res;
+ opr->v = malloc(sizeof(double complex) * par.res);
+ opr->pe = malloc(sizeof(double complex) * par.res);
+ opr->ke = malloc(sizeof(double complex) * par.res);
+ opr->wfc = malloc(sizeof(double complex) * par.res);
+
+ for (size_t i = 0; i < par.res; ++i) {
+ opr->v[i] = 0.5 * cpow(par.x[i] - voffset, 2);
+ opr->wfc[i] = cexp(-cpow(par.x[i] - wfcoffset, 2) / 2.0);
+
+ if (par.im_time) {
+ opr->ke[i] = cexp(-0.5 * par.dt * cpow(par.k[i], 2));
+ opr->pe[i] = cexp(-0.5 * par.dt * opr->v[i]);
+ } else {
+ opr->ke[i] = cexp(-0.5 * par.dt * cpow(par.k[i], 2) * I);
+ opr->pe[i] = cexp(-0.5 * par.dt * opr->v[i] * I);
+ }
+ }
+}
+
+void split_op(struct params par, struct operators opr) {
+ double density[opr.size];
+
+ for (size_t i = 0; i < par.timesteps; ++i) {
+ for (size_t j = 0; j < opr.size; ++j) {
+ opr.wfc[j] *= opr.pe[j];
+ }
+
+ fft(opr.wfc, opr.size, false);
+
+ for (size_t j = 0; j < opr.size; ++j) {
+ opr.wfc[j] *= opr.ke[j];
+ }
+
+ fft(opr.wfc, opr.size, true);
+
+ for (size_t j = 0; j < opr.size; ++j) {
+ opr.wfc[j] *= opr.pe[j];
+ }
+
+ for (size_t j = 0; j < opr.size; ++j) {
+ density[j] = pow(cabs(opr.wfc[j]), 2);
+ }
+
+ if (par.im_time) {
+ double sum = 0;
+
+ for (size_t j = 0; j < opr.size; ++j) {
+ sum += density[j];
+ }
+
+ sum *= par.dx;
+
+ for (size_t j = 0; j < opr.size; ++j) {
+ opr.wfc[j] /= sqrt(sum);
+ }
+ }
+
+ char filename[256];
+ sprintf(filename, "output%lu.dat", i);
+ FILE *fp = fopen(filename, "w");
+
+ for (int i = 0; i < opr.size; ++i) {
+ fprintf(fp, "%d\t%f\t%f\n", i, density[i], creal(opr.v[i]));
+ }
+
+ fclose(fp);
+ }
+}
+
+double calculate_energy(struct params par, struct operators opr) {
+ double complex wfc_r[opr.size];
+ double complex wfc_k[opr.size];
+ double complex wfc_c[opr.size];
+ memcpy(wfc_r, opr.wfc, sizeof(wfc_r));
+
+ memcpy(wfc_k, opr.wfc, sizeof(wfc_k));
+ fft(wfc_k, opr.size, false);
+
+ for (size_t i = 0; i < opr.size; ++i) {
+ wfc_c[i] = conj(wfc_r[i]);
+ }
+
+ double complex energy_k[opr.size];
+ double complex energy_r[opr.size];
+
+ for (size_t i = 0; i < opr.size; ++i) {
+ energy_k[i] = wfc_k[i] * cpow(par.k[i] + 0.0*I, 2);
+ }
+
+ fft(energy_k, opr.size, true);
+
+ for (size_t i = 0; i < opr.size; ++i) {
+ energy_k[i] *= 0.5 * wfc_c[i];
+ energy_r[i] = wfc_c[i] * opr.v[i] * wfc_r[i];
+ }
+
+ double energy_final = 0;
+
+ for (size_t i = 0; i < opr.size; ++i) {
+ energy_final += creal(energy_k[i] + energy_r[i]);
+ }
+
+ return energy_final * par.dx;
+}
+
+void free_params(struct params par) {
+ free(par.x);
+ free(par.k);
+}
+
+void free_operators(struct operators opr) {
+ free(opr.v);
+ free(opr.pe);
+ free(opr.ke);
+ free(opr.wfc);
+}
+
+int main() {
+ struct params par;
+ struct operators opr;
+
+ init_params(&par, 5.0, 256, 0.05, 100, true);
+ init_operators(&opr, par, 0.0, -1.0);
+
+ split_op(par, opr);
+
+ printf("the energy is %f\n", calculate_energy(par, opr));
+
+ free_params(par);
+ free_operators(opr);
+
+ return 0;
+}

contents/split-operator_method/split-operator_method.md

@@ -99,6 +99,9 @@ Regardless, we first need to set all the initial parameters, including the initi
 {% method %}
 {% sample lang="jl" %}
 [import:9-32, lang:"julia"](code/julia/split_op.jl)
+{% sample lang="c" %}
+[import:10-20, lang:"c_cpp"](code/c/split_op.c)
+[import:51-72, lang:"c_cpp"](code/c/split_op.c)
 {% endmethod %}
 
 As a note, when we generate our grid in momentum space `k`, we need to split the grid into two lines, one that is going from `0` to `-kmax` and is then discontinuous and goes from `kmax` to `0`.
@@ -111,6 +114,9 @@ Afterwards, we turn them into operators:
 {% method %}
 {% sample lang="jl" %}
 [import:34-60, lang:"julia"](code/julia/split_op.jl)
+{% sample lang="c" %}
+[import:22-28, lang:"c_cpp"](code/c/split_op.c)
+[import:74-95, lang:"c_cpp"](code/c/split_op.c)
 {% endmethod %}
 
 Here, we use a standard harmonic potential for the atoms to sit in and a gaussian distribution for an initial guess for the probability distribution.
@@ -124,6 +130,8 @@ The final step is to do the iteration, itself.
 {% method %}
 {% sample lang="jl" %}
 [import:63-109, lang:"julia"](code/julia/split_op.jl)
+{% sample lang="c" %}
+[import:97-145, lang:"c_cpp"](code/c/split_op.c)
 {% endmethod %}
 
 And that's it.
@@ -143,6 +151,8 @@ Checking to make sure your code can output the correct energy for a harmonic tra
 {% method %}
 {% sample lang="jl" %}
 [import, lang:"julia"](code/julia/split_op.jl)
+{% sample lang="c" %}
+[import, lang:"c_cpp"](code/c/split_op.c)
 {% endmethod %}
 
 <script>