Function compute_density_gradient

cpp_core/src/scf.cpp:2132–2184

};

DensityGradient compute_density_gradient(
    const PlaneWaveBasis& pw,
    const std::vector<double>& rho_r
) {
    const int nx = pw.fft_size[0];
    const int ny = pw.fft_size[1];
    const int nz = pw.fft_size[2];
    const int n_grid = nx * ny * nz;
    if (static_cast<int>(rho_r.size()) != n_grid) {
        throw std::invalid_argument("compute_density_gradient: rho_r size mismatch with FFT grid");
    }

    FFT3D fft(nx, ny, nz);
    std::vector<Vec3d> gcart_by_idx(static_cast<std::size_t>(n_grid), Vec3d{0.0, 0.0, 0.0});
    for (int ix = 0; ix < nx; ++ix) {
        const int gx = signed_miller_from_fft_index(ix, nx);
        for (int iy = 0; iy < ny; ++iy) {
            const int gy = signed_miller_from_fft_index(iy, ny);
            for (int iz = 0; iz < nz; ++iz) {
                const int gz = signed_miller_from_fft_index(iz, nz);
                const int idx = flat_index_3d(ix, iy, iz, ny, nz);
                gcart_by_idx[static_cast<std::size_t>(idx)] = matvec(pw.lattice.G, Vec3i{gx, gy, gz});
            }
        }
    }

    std::vector<Complex> rho_c(static_cast<std::size_t>(n_grid), Complex{0.0, 0.0});
    std::vector<Complex> rho_g(static_cast<std::size_t>(n_grid), Complex{0.0, 0.0});
    for (int i = 0; i < n_grid; ++i) {
        rho_c[static_cast<std::size_t>(i)] = Complex{rho_r[static_cast<std::size_t>(i)], 0.0};
    }
    fft.forward(rho_c.data(), rho_g.data());

    const Complex imag_unit{0.0, 1.0};
    std::vector<Complex> grad_g(static_cast<std::size_t>(n_grid), Complex{0.0, 0.0});
    std::vector<Complex> grad_r_c(static_cast<std::size_t>(n_grid), Complex{0.0, 0.0});
    std::vector<Vec3d> grad(static_cast<std::size_t>(n_grid), Vec3d{0.0, 0.0, 0.0});
    for (int comp = 0; comp < 3; ++comp) {
        for (int idx = 0; idx < n_grid; ++idx) {
            const double g_comp = gcart_by_idx[static_cast<std::size_t>(idx)][comp];
            grad_g[static_cast<std::size_t>(idx)] = imag_unit * g_comp * rho_g[static_cast<std::size_t>(idx)];
        }
        fft.backward(grad_g.data(), grad_r_c.data());
        for (int idx = 0; idx < n_grid; ++idx) {
            grad[static_cast<std::size_t>(idx)][comp] = std::real(grad_r_c[static_cast<std::size_t>(idx)]);
        }
    }

    DensityGradient out;
    out.grad = std::move(grad);
    out.gcart_by_idx = std::move(gcart_by_idx);
    return out;
}

XcEvaluation evaluate_gga_xc_unpolarized(
    const PlaneWaveBasis& pw,