| 2404 | } |
| 2405 | |
| 2406 | XcEvaluationSpin evaluate_xc_spin_polarized( |
| 2407 | const PlaneWaveBasis& pw, |
| 2408 | const std::vector<double>& rho_up_r, |
| 2409 | const std::vector<double>& rho_dn_r, |
| 2410 | const std::string& xc_label |
| 2411 | ) { |
| 2412 | const int n_grid = pw.fft_size[0] * pw.fft_size[1] * pw.fft_size[2]; |
| 2413 | if (static_cast<int>(rho_up_r.size()) != n_grid || static_cast<int>(rho_dn_r.size()) != n_grid) { |
| 2414 | throw std::invalid_argument("evaluate_xc_spin_polarized: rho size mismatch with FFT grid"); |
| 2415 | } |
| 2416 | |
| 2417 | const std::string xc_kind = normalize_xc_label(xc_label); |
| 2418 | XcEvaluationSpin out; |
| 2419 | out.eps_xc_ha.assign(static_cast<std::size_t>(n_grid), 0.0); |
| 2420 | out.v_xc_up_ha.assign(static_cast<std::size_t>(n_grid), 0.0); |
| 2421 | out.v_xc_dn_ha.assign(static_cast<std::size_t>(n_grid), 0.0); |
| 2422 | |
| 2423 | if (xc_kind.empty() || xc_kind == "NONE" || xc_kind == "NOXC" || xc_kind == "NOX" || xc_kind == "HARTREE") { |
| 2424 | out.functional_name = "none"; |
| 2425 | return out; |
| 2426 | } |
| 2427 | |
| 2428 | #ifndef MaterialDFT_HAS_LIBXC |
| 2429 | (void)pw; |
| 2430 | (void)rho_up_r; |
| 2431 | (void)rho_dn_r; |
| 2432 | (void)xc_label; |
| 2433 | throw std::runtime_error("evaluate_xc_spin_polarized: Libxc support is disabled (MaterialDFT_USE_LIBXC=OFF)"); |
| 2434 | #else |
| 2435 | std::vector<double> rho_up_clamped(static_cast<std::size_t>(n_grid), 0.0); |
| 2436 | std::vector<double> rho_dn_clamped(static_cast<std::size_t>(n_grid), 0.0); |
| 2437 | std::vector<double> rho_pair(static_cast<std::size_t>(2 * n_grid), 0.0); |
| 2438 | for (int i = 0; i < n_grid; ++i) { |
| 2439 | const double up = std::max(0.0, rho_up_r[static_cast<std::size_t>(i)]); |
| 2440 | const double dn = std::max(0.0, rho_dn_r[static_cast<std::size_t>(i)]); |
| 2441 | rho_up_clamped[static_cast<std::size_t>(i)] = up; |
| 2442 | rho_dn_clamped[static_cast<std::size_t>(i)] = dn; |
| 2443 | rho_pair[static_cast<std::size_t>(2 * i)] = up; |
| 2444 | rho_pair[static_cast<std::size_t>(2 * i + 1)] = dn; |
| 2445 | } |
| 2446 | |
| 2447 | if (xc_kind.rfind("LDA", 0) == 0) { |
| 2448 | xc_func_type lda_x; |
| 2449 | if (xc_func_init(&lda_x, XC_LDA_X, XC_POLARIZED) != 0) { |
| 2450 | throw std::runtime_error("evaluate_xc_spin_polarized: failed to initialize LDA exchange"); |
| 2451 | } |
| 2452 | xc_func_type lda_c; |
| 2453 | #ifdef XC_LDA_C_PW |
| 2454 | const int corr_id = XC_LDA_C_PW; |
| 2455 | const char* corr_name = "LDA_C_PW"; |
| 2456 | #elif defined(XC_LDA_C_PZ) |
| 2457 | const int corr_id = XC_LDA_C_PZ; |
| 2458 | const char* corr_name = "LDA_C_PZ"; |
| 2459 | #else |
| 2460 | xc_func_end(&lda_x); |
| 2461 | throw std::runtime_error("evaluate_xc_spin_polarized: no supported LDA correlation functional"); |
| 2462 | #endif |
| 2463 | if (xc_func_init(&lda_c, corr_id, XC_POLARIZED) != 0) { |
no test coverage detected