diff --git a/contents/split-operator_method/code/c++/split_op.cpp b/contents/split-operator_method/code/c++/split_op.cpp new file mode 100644 index 000000000..efdad46c1 --- /dev/null +++ b/contents/split-operator_method/code/c++/split_op.cpp @@ -0,0 +1,201 @@ +#include +#include +#include +#include +#include + +// Using fftw3 library. +#include + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +using complex = std::complex; +using vector_real = std::vector; +using vector_complex = std::vector; + +struct Params { + Params(double _xmax, unsigned int _res, double _dt, unsigned int _timesteps, bool im) { + xmax = _xmax; + res = _res; + dt = _dt; + timesteps = _timesteps; + dx = 2.0 * xmax / res; + x.reserve(res); + dk = M_PI / xmax; + k.reserve(res); + im_time = im; + + for (size_t i = 0; i < res; ++i) { + x.emplace_back(xmax / res - xmax + i * (2.0 * xmax / res)); + if (i < res / 2) { + k.push_back(i * M_PI / xmax); + } else { + k.push_back((static_cast(i) - res) * M_PI / xmax); + } + } + } + + double xmax; + unsigned int res; + double dt; + unsigned int timesteps; + double dx; + vector_real x; + double dk; + vector_real k; + bool im_time; +}; + +struct Operators { +public: + Operators(Params &par, double voffset, + double wfcoffset) { + size = par.res; + v.reserve(size); + pe.reserve(size); + ke.reserve(size); + wfc.reserve(size); + + for (size_t i = 0; i < size; ++i) { + v.push_back(0.5 * pow(par.x[i] - voffset, 2)); + wfc.push_back(exp(-pow(par.x[i] - wfcoffset, 2) / 2.0)); + + if (par.im_time) { + ke.push_back(exp(-0.5 * par.dt * pow(par.k[i], 2))); + pe.push_back(exp(-0.5 * par.dt * v[i])); + } else { + ke.push_back(exp(-0.5 * par.dt * pow(par.k[i], 2) * complex(0.0, 1.0))); + pe.push_back(exp(-0.5 * par.dt * v[i] * complex(0.0, 1.0))); + } + } + } + + size_t size; + vector_complex v; + vector_complex pe; + vector_complex ke; + vector_complex wfc; +}; + +void fft(vector_complex &x, int n, bool inverse) { + complex y[n]; + memset(y, 0, sizeof(y)); + fftw_plan p; + + fftw_complex *in = reinterpret_cast(x.data()); + fftw_complex *out = reinterpret_cast(y); + p = fftw_plan_dft_1d(n, in, out, + (inverse ? FFTW_BACKWARD : FFTW_FORWARD), FFTW_ESTIMATE); + + fftw_execute(p); + fftw_destroy_plan(p); + + for (size_t i = 0; i < n; ++i) { + x[i] = y[i] / sqrt(static_cast(n)); + } +} + +void split_op(Params &par, 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(abs(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); + } + } + + // Writing data into a file in the format of: + // index, density, real potential. + std::stringstream filename_stream; + filename_stream << "output" << i << ".dat"; + + std::ofstream fstream = std::ofstream(filename_stream.str()); + + if (fstream) { + for (int i = 0; i < opr.size; ++i) { + std::stringstream data_stream; + + data_stream << i << "\t" << density[i] << "\t" << real(opr.v[i]) << "\n"; + + fstream.write(data_stream.str().c_str(), data_stream.str().length()); + } + } + + fstream.close(); + } +} + +double calculate_energy(Params &par, Operators &opr) { + vector_complex wfc_r(opr.wfc); + vector_complex wfc_k(opr.wfc); + vector_complex wfc_c(opr.size); + fft(wfc_k, opr.size, false); + + for (size_t i = 0; i < opr.size; ++i) { + wfc_c[i] = conj(wfc_r[i]); + } + + vector_complex energy_k(opr.size); + vector_complex energy_r(opr.size); + + for (size_t i = 0; i < opr.size; ++i) { + energy_k[i] = wfc_k[i] * pow(complex(par.k[i], 0.0), 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 += real(energy_k[i] + energy_r[i]); + } + + return energy_final * par.dx; +} + +int main() { + Params par = Params(5.0, 256, 0.05, 100, true); + Operators opr = Operators(par, 0.0, -1.0); + + split_op(par, opr); + + std::cout << "The energy is " << calculate_energy(par, opr) << "\n"; + + return 0; +} diff --git a/contents/split-operator_method/split-operator_method.md b/contents/split-operator_method/split-operator_method.md index 60f2f33ee..ff74f99d9 100644 --- a/contents/split-operator_method/split-operator_method.md +++ b/contents/split-operator_method/split-operator_method.md @@ -102,6 +102,8 @@ Regardless, we first need to set all the initial parameters, including the initi {% sample lang="c" %} [import:11-21, lang:"c_cpp"](code/c/split_op.c) [import:52-73, lang:"c_cpp"](code/c/split_op.c) +{% sample lang="cpp" %} +[import:14-49, lang:"c_cpp"](code/c++/split_op.cpp) {% sample lang="py" %} [import:11-30, lang:"python"](code/python/split_op.py) {% sample lang="hs" %} @@ -121,6 +123,8 @@ Afterwards, we turn them into operators: {% sample lang="c" %} [import:23-29, lang:"c_cpp"](code/c/split_op.c) [import:75-96, lang:"c_cpp"](code/c/split_op.c) +{% sample lang="cpp" %} +[import:51-80, lang:"c_cpp"](code/c++/split_op.cpp) {% sample lang="py" %} [import:33-54, lang:"python"](code/python/split_op.py) {% sample lang="hs" %} @@ -140,6 +144,8 @@ The final step is to do the iteration, itself. [import:65-112, lang:"julia"](code/julia/split_op.jl) {% sample lang="c" %} [import:98-148, lang:"c_cpp"](code/c/split_op.c) +{% sample lang="cpp" %} +[import:100-157, lang:"c_cpp"](code/c++/split_op.cpp) {% sample lang="py" %} [import:57-95, lang:"python"](code/python/split_op.py) {% sample lang="hs" %} @@ -165,6 +171,8 @@ Checking to make sure your code can output the correct energy for a harmonic tra [import, lang:"julia"](code/julia/split_op.jl) {% sample lang="c" %} [import, lang:"c_cpp"](code/c/split_op.c) +{% sample lang="cpp" %} +[import, lang:"c_cpp"](code/c++/split_op.cpp) {% sample lang="py" %} [import:5-127, lang:"python"](code/python/split_op.py) {% sample lang="hs" %}