modifying C++ energy calculation to match split-op code. (#533)

Author: leios

contents/quantum_systems/code/c++/energy.cpp

@@ -3,48 +3,50 @@
 
 #include <fftw3.h>
 
-void fft(std::vector<std::complex<double>> x, bool inverse) {
-    std::vector<std::complex<double>> y(x.size());
-    for (size_t i = 0; i < x.size(); i++) {
-        y.push_back(std::complex(0.0, 0.0));
-    }
-    
+void fft(std::vector<std::complex<double>> &x, bool inverse) {
+    std::vector<std::complex<double>> y(x.size(), std::complex<double>(0.0, 0.0));
+
     fftw_plan p;
 
-    p = fftw_plan_dft_1d(x.size(), reinterpret_cast<fftw_complex*>(x.data()),
-                         reinterpret_cast<fftw_complex*>(y.data()),
-                         (inverse ? FFTW_BACKWARD : FFTW_FORWARD), FFTW_ESTIMATE);
+    fftw_complex *in = reinterpret_cast<fftw_complex*>(x.data());
+    fftw_complex *out = reinterpret_cast<fftw_complex*>(y.data());
+
+    p = fftw_plan_dft_1d(x.size(), in, out,
+                        (inverse ? FFTW_BACKWARD : FFTW_FORWARD), FFTW_ESTIMATE);
+
 
     fftw_execute(p);
     fftw_destroy_plan(p);
 
     for (size_t i = 0; i < x.size(); ++i) {
-        x[i] = y[i] / sqrt((double)x.size());
+        x[i] = y[i] / sqrt(static_cast<double>(x.size()));
     }
 }
 
-double calculate_energy(std::vector<std::complex<double>> wfc, std::vector<std::complex<double>> h_r,
-                        std::vector<std::complex<double>> h_k, double dx, size_t size) {
+double calculate_energy(std::vector<std::complex<double>> wfc,
+                        std::vector<std::complex<double>> h_r,
+                        std::vector<std::complex<double>> h_k,
+                        double dx, size_t size) {
     std::vector<std::complex<double>> wfc_k(wfc);
     std::vector<std::complex<double>> wfc_c(size);
     fft(wfc_k, false);
 
     for (size_t i = 0; i < size; ++i) {
-        wfc_c.push_back(conj(wfc[i]));
+        wfc_c[i] = conj(wfc[i]);
     }
 
     std::vector<std::complex<double>> energy_k(size);
     std::vector<std::complex<double>> energy_r(size);
 
     for (size_t i = 0; i < size; ++i) {
-        energy_k.push_back(wfc_k[i] * h_k[i]);
+        energy_k[i] = wfc_k[i] * pow(h_k[i], 2);
     }
 
     fft(energy_k, true);
 
     for (size_t i = 0; i < size; ++i) {
-        energy_k[i] *= wfc_c[i];
-        energy_r.push_back(wfc_c[i] * h_r[i] * wfc[i]);
+        energy_k[i] *= 0.5 * wfc_c[i];
+        energy_r[i] = wfc_c[i] * h_r[i] * wfc[i];
     }
 
     double energy_final = 0;

contents/quantum_systems/quantum_systems.md

@@ -232,7 +232,7 @@ This ultimately looks like this:
 {% sample lang="c" %}
 [import:29-, lang:"c_cpp"](code/c/energy.c)
 {% sample lang="cpp" %}
-[import:26-57, lang:"c_cpp"](code/c++/energy.cpp)
+[import:26-, lang:"c_cpp"](code/c++/energy.cpp)
 {% sample lang="py" %}
 [import:4-17, lang:"python"](code/python/energy.py)
 {% endmethod %}

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

@@ -79,21 +79,20 @@ struct Operators {
     vector_complex wfc;
 };
 
-void fft(vector_complex &x, int n, bool inverse) {
-    complex y[n];
-    memset(y, 0, sizeof(y));
+void fft(vector_complex &x, bool inverse) {
+    std::vector<std::complex<double>> y(x.size(), std::complex<double>(0.0, 0.0));
     fftw_plan p;
 
     fftw_complex *in = reinterpret_cast<fftw_complex*>(x.data());
-    fftw_complex *out = reinterpret_cast<fftw_complex*>(y);
-    p = fftw_plan_dft_1d(n, in, out,
+    fftw_complex *out = reinterpret_cast<fftw_complex*>(y.data());
+    p = fftw_plan_dft_1d(x.size(), 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<double>(n));
+    for (size_t i = 0; i < x.size(); ++i) {
+        x[i] = y[i] / sqrt(static_cast<double>(x.size()));
     }
 }
 
@@ -105,13 +104,13 @@ void split_op(Params &par, Operators &opr) {
             opr.wfc[j] *= opr.pe[j];
         }
 
-        fft(opr.wfc, opr.size, false);
+        fft(opr.wfc, false);
 
         for (size_t j = 0; j < opr.size; ++j) {
             opr.wfc[j] *= opr.ke[j];
         }
 
-        fft(opr.wfc, opr.size, true);
+        fft(opr.wfc, true);
 
         for (size_t j = 0; j < opr.size; ++j) {
             opr.wfc[j] *= opr.pe[j];
@@ -160,7 +159,7 @@ 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);
+    fft(wfc_k, false);
 
     for (size_t i = 0; i < opr.size; ++i) {
         wfc_c[i] = conj(wfc_r[i]);
@@ -173,7 +172,7 @@ double calculate_energy(Params &par, Operators &opr) {
         energy_k[i] = wfc_k[i] * pow(complex(par.k[i], 0.0), 2);
     }
 
-    fft(energy_k, opr.size, true);
+    fft(energy_k, true);
 
     for (size_t i = 0; i < opr.size; ++i) {
         energy_k[i] *= 0.5 * wfc_c[i];

contents/split-operator_method/split-operator_method.md

@@ -145,7 +145,7 @@ The final step is to do the iteration, itself.
 {% 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)
+[import:99-156, lang:"c_cpp"](code/c++/split_op.cpp)
 {% sample lang="py" %}
 [import:57-95, lang:"python"](code/python/split_op.py)
 {% sample lang="hs" %}