Skip to content

Commit 2fa7a20

Browse files
Gathrosleios
Gathros
authored and
leios
committed
Adding split_op.c and energy.c (algorithm-archivists#308)
* Adding split_op.c * updating split-operator_method.md * Adding energy.c * updating quantum_systems.md * Changing energy.c to match the julia code * fixing energy.c * fixing quantum_systems.md * Change ints to unsigned ints in split_op.c * adding comments to split_op.c
1 parent f69c197 commit 2fa7a20

File tree

4 files changed

+285
-0
lines changed

4 files changed

+285
-0
lines changed

‎contents/quantum_systems/code/c/energy.c

+60
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
#include <complex.h>
2+
#include <string.h>
3+
#include <math.h>
4+
5+
#include <fftw3.h>
6+
7+
void fft(double complex *x, int n, bool inverse) {
8+
double complex y[n];
9+
memset(y, 0, sizeof(y));
10+
fftw_plan p;
11+
12+
if (inverse) {
13+
p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
14+
FFTW_BACKWARD, FFTW_ESTIMATE);
15+
} else {
16+
p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
17+
FFTW_FORWARD, FFTW_ESTIMATE);
18+
}
19+
20+
fftw_execute(p);
21+
fftw_destroy_plan(p);
22+
23+
for (size_t i = 0; i < n; ++i) {
24+
x[i] = y[i] / sqrt((double)n);
25+
}
26+
}
27+
28+
double calculate_energy(double complex *wfc, double complex *h_r,
29+
double complex *h_k, double dx, size_t size) {
30+
double complex wfc_k[size];
31+
double complex wfc_c[size];
32+
memcpy(wfc_k, wfc, sizeof(wfc_k));
33+
fft(wfc_k, size, false);
34+
35+
for (size_t i = 0; i < size; ++i) {
36+
wfc_c[i] = conj(wfc[i]);
37+
}
38+
39+
double complex energy_k[size];
40+
double complex energy_r[size];
41+
42+
for (size_t i = 0; i < size; ++i) {
43+
energy_k[i] = wfc_k[i] * h_k[i];
44+
}
45+
46+
fft(energy_k, size, true);
47+
48+
for (size_t i = 0; i < size; ++i) {
49+
energy_k[i] *= wfc_c[i];
50+
energy_r[i] = wfc_c[i] * H_r[i] * wfc[i];
51+
}
52+
53+
double energy_final = 0;
54+
55+
for (size_t i = 0; i < size; ++i) {
56+
energy_final += creal(energy_k[i] + energy_r[i]);
57+
}
58+
59+
return energy_final * dx;
60+
}

‎contents/quantum_systems/quantum_systems.md

+2
Original file line numberDiff line numberDiff line change
@@ -227,6 +227,8 @@ This ultimately looks like this:
227227
{% method %}
228228
{% sample lang="jl" %}
229229
[import, lang:"julia"](code/julia/energy.jl)
230+
{% sample lang="c" %}
231+
[import:28-, lang:"c_cpp"](code/c/energy.c)
230232
{% sample lang="py" %}
231233
[import:4-17, lang:"python"](code/python/energy.py)
232234
{% endmethod %}

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

+213
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,213 @@
1+
#include <complex.h>
2+
#include <stdbool.h>
3+
#include <stdio.h>
4+
#include <stdlib.h>
5+
#include <string.h>
6+
#include <math.h>
7+
8+
// Using fftw3 library.
9+
#include <fftw3.h>
10+
11+
struct params {
12+
double xmax;
13+
unsigned int res;
14+
double dt;
15+
unsigned int timesteps;
16+
double dx;
17+
double *x;
18+
double dk;
19+
double *k;
20+
bool im_time;
21+
};
22+
23+
struct operators {
24+
size_t size;
25+
double complex *v;
26+
double complex *pe;
27+
double complex *ke;
28+
double complex *wfc;
29+
};
30+
31+
void fft(double complex *x, int n, bool inverse) {
32+
double complex y[n];
33+
memset(y, 0, sizeof(y));
34+
fftw_plan p;
35+
36+
if (inverse) {
37+
p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
38+
FFTW_BACKWARD, FFTW_ESTIMATE);
39+
} else {
40+
p = fftw_plan_dft_1d(n, (fftw_complex*)x, (fftw_complex*)y,
41+
FFTW_FORWARD, FFTW_ESTIMATE);
42+
}
43+
44+
fftw_execute(p);
45+
fftw_destroy_plan(p);
46+
47+
for (size_t i = 0; i < n; ++i) {
48+
x[i] = y[i] / sqrt((double)n);
49+
}
50+
}
51+
52+
void init_params(struct params *par, double xmax, unsigned int res, double dt,
53+
unsigned int timesteps, bool im) {
54+
55+
par->xmax = xmax;
56+
par->res = res;
57+
par->dt = dt;
58+
par->timesteps = timesteps;
59+
par->dx = 2.0 * xmax / res;
60+
par->x = malloc(sizeof(double) * res);
61+
par->dk = M_PI / xmax;
62+
par->k = malloc(sizeof(double) * res);
63+
par->im_time = im;
64+
65+
for (size_t i = 0; i < res; ++i) {
66+
par->x[i] = xmax / res - xmax + i * (2.0 * xmax / res);
67+
if (i < res / 2) {
68+
par->k[i] = i * M_PI / xmax;
69+
} else {
70+
par->k[i] = ((double)i - res) * M_PI / xmax;
71+
}
72+
}
73+
}
74+
75+
void init_operators(struct operators *opr, struct params par, double voffset,
76+
double wfcoffset) {
77+
78+
opr->size = par.res;
79+
opr->v = malloc(sizeof(double complex) * par.res);
80+
opr->pe = malloc(sizeof(double complex) * par.res);
81+
opr->ke = malloc(sizeof(double complex) * par.res);
82+
opr->wfc = malloc(sizeof(double complex) * par.res);
83+
84+
for (size_t i = 0; i < par.res; ++i) {
85+
opr->v[i] = 0.5 * cpow(par.x[i] - voffset, 2);
86+
opr->wfc[i] = cexp(-cpow(par.x[i] - wfcoffset, 2) / 2.0);
87+
88+
if (par.im_time) {
89+
opr->ke[i] = cexp(-0.5 * par.dt * cpow(par.k[i], 2));
90+
opr->pe[i] = cexp(-0.5 * par.dt * opr->v[i]);
91+
} else {
92+
opr->ke[i] = cexp(-0.5 * par.dt * cpow(par.k[i], 2) * I);
93+
opr->pe[i] = cexp(-0.5 * par.dt * opr->v[i] * I);
94+
}
95+
}
96+
}
97+
98+
void split_op(struct params par, struct operators opr) {
99+
double density[opr.size];
100+
101+
for (size_t i = 0; i < par.timesteps; ++i) {
102+
for (size_t j = 0; j < opr.size; ++j) {
103+
opr.wfc[j] *= opr.pe[j];
104+
}
105+
106+
fft(opr.wfc, opr.size, false);
107+
108+
for (size_t j = 0; j < opr.size; ++j) {
109+
opr.wfc[j] *= opr.ke[j];
110+
}
111+
112+
fft(opr.wfc, opr.size, true);
113+
114+
for (size_t j = 0; j < opr.size; ++j) {
115+
opr.wfc[j] *= opr.pe[j];
116+
}
117+
118+
for (size_t j = 0; j < opr.size; ++j) {
119+
density[j] = pow(cabs(opr.wfc[j]), 2);
120+
}
121+
122+
if (par.im_time) {
123+
double sum = 0;
124+
125+
for (size_t j = 0; j < opr.size; ++j) {
126+
sum += density[j];
127+
}
128+
129+
sum *= par.dx;
130+
131+
for (size_t j = 0; j < opr.size; ++j) {
132+
opr.wfc[j] /= sqrt(sum);
133+
}
134+
}
135+
136+
// Writing data into a file in the format of:
137+
// index, density, real potential.
138+
char filename[256];
139+
sprintf(filename, "output%lu.dat", i);
140+
FILE *fp = fopen(filename, "w");
141+
142+
for (int i = 0; i < opr.size; ++i) {
143+
fprintf(fp, "%d\t%f\t%f\n", i, density[i], creal(opr.v[i]));
144+
}
145+
146+
fclose(fp);
147+
}
148+
}
149+
150+
double calculate_energy(struct params par, struct operators opr) {
151+
double complex wfc_r[opr.size];
152+
double complex wfc_k[opr.size];
153+
double complex wfc_c[opr.size];
154+
memcpy(wfc_r, opr.wfc, sizeof(wfc_r));
155+
156+
memcpy(wfc_k, opr.wfc, sizeof(wfc_k));
157+
fft(wfc_k, opr.size, false);
158+
159+
for (size_t i = 0; i < opr.size; ++i) {
160+
wfc_c[i] = conj(wfc_r[i]);
161+
}
162+
163+
double complex energy_k[opr.size];
164+
double complex energy_r[opr.size];
165+
166+
for (size_t i = 0; i < opr.size; ++i) {
167+
energy_k[i] = wfc_k[i] * cpow(par.k[i] + 0.0*I, 2);
168+
}
169+
170+
fft(energy_k, opr.size, true);
171+
172+
for (size_t i = 0; i < opr.size; ++i) {
173+
energy_k[i] *= 0.5 * wfc_c[i];
174+
energy_r[i] = wfc_c[i] * opr.v[i] * wfc_r[i];
175+
}
176+
177+
double energy_final = 0;
178+
179+
for (size_t i = 0; i < opr.size; ++i) {
180+
energy_final += creal(energy_k[i] + energy_r[i]);
181+
}
182+
183+
return energy_final * par.dx;
184+
}
185+
186+
void free_params(struct params par) {
187+
free(par.x);
188+
free(par.k);
189+
}
190+
191+
void free_operators(struct operators opr) {
192+
free(opr.v);
193+
free(opr.pe);
194+
free(opr.ke);
195+
free(opr.wfc);
196+
}
197+
198+
int main() {
199+
struct params par;
200+
struct operators opr;
201+
202+
init_params(&par, 5.0, 256, 0.05, 100, true);
203+
init_operators(&opr, par, 0.0, -1.0);
204+
205+
split_op(par, opr);
206+
207+
printf("the energy is %f\n", calculate_energy(par, opr));
208+
209+
free_params(par);
210+
free_operators(opr);
211+
212+
return 0;
213+
}

‎contents/split-operator_method/split-operator_method.md

+10
Original file line numberDiff line numberDiff line change
@@ -99,6 +99,9 @@ Regardless, we first need to set all the initial parameters, including the initi
9999
{% method %}
100100
{% sample lang="jl" %}
101101
[import:9-32, lang:"julia"](code/julia/split_op.jl)
102+
{% sample lang="c" %}
103+
[import:10-20, lang:"c_cpp"](code/c/split_op.c)
104+
[import:51-72, lang:"c_cpp"](code/c/split_op.c)
102105
{% sample lang="py" %}
103106
[import:11-30, lang:"python"](code/python/split_op.py)
104107
{% endmethod %}
@@ -113,6 +116,9 @@ Afterwards, we turn them into operators:
113116
{% method %}
114117
{% sample lang="jl" %}
115118
[import:34-60, lang:"julia"](code/julia/split_op.jl)
119+
{% sample lang="c" %}
120+
[import:22-28, lang:"c_cpp"](code/c/split_op.c)
121+
[import:74-95, lang:"c_cpp"](code/c/split_op.c)
116122
{% sample lang="py" %}
117123
[import:33-54, lang:"python"](code/python/split_op.py)
118124
{% endmethod %}
@@ -128,6 +134,8 @@ The final step is to do the iteration, itself.
128134
{% method %}
129135
{% sample lang="jl" %}
130136
[import:63-109, lang:"julia"](code/julia/split_op.jl)
137+
{% sample lang="c" %}
138+
[import:97-145, lang:"c_cpp"](code/c/split_op.c)
131139
{% sample lang="py" %}
132140
[import:57-95, lang:"python"](code/python/split_op.py)
133141
{% endmethod %}
@@ -149,6 +157,8 @@ Checking to make sure your code can output the correct energy for a harmonic tra
149157
{% method %}
150158
{% sample lang="jl" %}
151159
[import, lang:"julia"](code/julia/split_op.jl)
160+
{% sample lang="c" %}
161+
[import, lang:"c_cpp"](code/c/split_op.c)
152162
{% sample lang="py" %}
153163
[import:5-127, lang:"python"](code/python/split_op.py)
154164
{% endmethod %}

0 commit comments

Comments
 (0)