| 25 | } // namespace |
| 26 | |
| 27 | std::vector<Complex> hartree_potential_g( |
| 28 | const std::vector<Complex>& rho_g, |
| 29 | const std::vector<double>& g2_bohr_inv2 |
| 30 | ) { |
| 31 | if (rho_g.size() != g2_bohr_inv2.size()) { |
| 32 | throw std::invalid_argument("hartree_potential_g: rho_g and g2 must have the same length"); |
| 33 | } |
| 34 | |
| 35 | // Compute in eV-Ang unit. |
| 36 | constexpr double FOUR_PI_HA_TO_EV = 4.0 * PI * HA_TO_EV; |
| 37 | |
| 38 | std::vector<Complex> out(rho_g.size(), Complex{0.0, 0.0}); |
| 39 | for (std::size_t i = 0; i < rho_g.size(); ++i) { |
| 40 | const double g2 = g2_bohr_inv2[i]; |
| 41 | if (g2 < 0.0) { |
| 42 | throw std::invalid_argument("hartree_potential_g: g2 entries must be non-negative"); |
| 43 | } |
| 44 | if (g2 <= 1e-14) { |
| 45 | out[i] = Complex{0.0, 0.0}; |
| 46 | continue; |
| 47 | } |
| 48 | out[i] = (FOUR_PI_HA_TO_EV / g2) * rho_g[i]; |
| 49 | } |
| 50 | return out; |
| 51 | } |
| 52 | |
| 53 | bool libxc_available() { |
| 54 | #ifdef MaterialDFT_HAS_LIBXC |