|
| 1 | +#include <vector> |
| 2 | +#include <complex> |
| 3 | + |
| 4 | +#include <fftw3.h> |
| 5 | + |
| 6 | +void fft(std::vector<std::complex<double>> x, bool inverse) { |
| 7 | + std::vector<std::complex<double>> y(x.size()); |
| 8 | + for (size_t i = 0; i < x.size(); i++) { |
| 9 | + y.push_back(std::complex(0.0, 0.0)); |
| 10 | + } |
| 11 | + |
| 12 | + fftw_plan p; |
| 13 | + |
| 14 | + p = fftw_plan_dft_1d(x.size(), reinterpret_cast<fftw_complex*>(x.data()), |
| 15 | + reinterpret_cast<fftw_complex*>(y.data()), |
| 16 | + (inverse ? FFTW_BACKWARD : FFTW_FORWARD), FFTW_ESTIMATE); |
| 17 | + |
| 18 | + fftw_execute(p); |
| 19 | + fftw_destroy_plan(p); |
| 20 | + |
| 21 | + for (size_t i = 0; i < x.size(); ++i) { |
| 22 | + x[i] = y[i] / sqrt((double)x.size()); |
| 23 | + } |
| 24 | +} |
| 25 | + |
| 26 | +double calculate_energy(std::vector<std::complex<double>> wfc, std::vector<std::complex<double>> h_r, |
| 27 | + std::vector<std::complex<double>> h_k, double dx, size_t size) { |
| 28 | + std::vector<std::complex<double>> wfc_k(wfc); |
| 29 | + std::vector<std::complex<double>> wfc_c(size); |
| 30 | + fft(wfc_k, false); |
| 31 | + |
| 32 | + for (size_t i = 0; i < size; ++i) { |
| 33 | + wfc_c.push_back(conj(wfc[i])); |
| 34 | + } |
| 35 | + |
| 36 | + std::vector<std::complex<double>> energy_k(size); |
| 37 | + std::vector<std::complex<double>> energy_r(size); |
| 38 | + |
| 39 | + for (size_t i = 0; i < size; ++i) { |
| 40 | + energy_k.push_back(wfc_k[i] * h_k[i]); |
| 41 | + } |
| 42 | + |
| 43 | + fft(energy_k, true); |
| 44 | + |
| 45 | + for (size_t i = 0; i < size; ++i) { |
| 46 | + energy_k[i] *= wfc_c[i]; |
| 47 | + energy_r.push_back(wfc_c[i] * h_r[i] * wfc[i]); |
| 48 | + } |
| 49 | + |
| 50 | + double energy_final = 0; |
| 51 | + |
| 52 | + for (size_t i = 0; i < size; ++i) { |
| 53 | + energy_final += real(energy_k[i] + energy_r[i]); |
| 54 | + } |
| 55 | + |
| 56 | + return energy_final * dx; |
| 57 | +} |
0 commit comments