feat(librpa): add Delta-Sternheimer response solver

CMakeLists.txt

@@ -79,6 +79,14 @@ cmake_dependent_option(ENABLE_PEXSI "Enable PEXSI for LCAO" OFF
                        "ENABLE_LCAO;ENABLE_MPI" OFF)
 cmake_dependent_option(ENABLE_MLALGO "Enable the machine learning algorithms"
                        OFF "ENABLE_LCAO;ENABLE_MPI" OFF)
+# Optional external GreenX minimax API for the Delta-Sternheimer solver.  Off by
+# default so a stock build needs neither Fortran nor GreenX.
+cmake_dependent_option(
+  ENABLE_GREENX_MINIMAX
+  "Enable the external GreenX minimax API for Delta-Sternheimer" OFF
+  "ENABLE_LCAO;ENABLE_MPI;ENABLE_LIBRI" OFF)
+set(GREENX_DIR "" CACHE PATH
+    "GreenX installation prefix or CMake package directory")
 
 # Two-center FFT is only used in LCAO
 cmake_dependent_option(
@@ -635,6 +643,11 @@ if(ENABLE_LIBRI)
   find_package(cereal REQUIRED CONFIG)
 endif()
 
+if(ENABLE_GREENX_MINIMAX)
+  enable_language(Fortran)
+  find_package(GreenX REQUIRED)
+endif()
+
 if(ENABLE_LIBXC)
   find_package(Libxc CONFIG REQUIRED)
   if(Libxc_VERSION VERSION_LESS "5.1.7")

source/CMakeLists.txt

@@ -836,7 +836,7 @@ endif()
 if(ENABLE_LIBRI)
   target_link_libraries(${ABACUS_BIN_NAME} PRIVATE ri module_exx_symmetry)
   if(ENABLE_LCAO)
-    target_link_libraries(${ABACUS_BIN_NAME} PRIVATE rpa librpa_io)
+    target_link_libraries(${ABACUS_BIN_NAME} PRIVATE rpa delta_st librpa_io)
   endif()
 endif()
 if(ENABLE_JSON)

source/Makefile.Objects

@@ -82,6 +82,7 @@ VPATH=./src_global:\
 ./source_lcao/module_ri:\
 ./source_lcao/module_ri/module_exx_symmetry:\
 ./source_lcao/module_rpa:\
+./source_lcao/module_rpa/module_delta_st:\
 ./source_lcao/module_rpa/module_librpa_io:\
 ./source_lcao/module_bse:\
 ./source_io/module_parameter:\
@@ -937,6 +938,46 @@ OBJS_MODULE_RPA=rpa_abfs_preorth.o\
       rpa_lri_basis.o\
       rpa_lri_eigenvector_v1.o\
       rpa_lri_eigenvector.o\
+      sternheimer_abacus_dft_reference.o\
+      sternheimer_abacus_dft_zero_order.o\
+      sternheimer_abacus_fd_adapter.o\
+      sternheimer_abacus_fd_bindings.o\
+      sternheimer_abacus_fd_nonlocal.o\
+      sternheimer_abacus_fd_smoke.o\
+      sternheimer_abacus_st_smoke.o\
+      sternheimer_abfs_perturbation.o\
+      sternheimer_abfs_strict2d.o\
+      sternheimer_channel_resources.o\
+      sternheimer_chi0_mpi.o\
+      sternheimer_coulomb_whitening.o\
+      sternheimer_delta.o\
+      sternheimer_dft_zero_order.o\
+      sternheimer_fd_hamiltonian.o\
+      sternheimer_fd_nonlocal_projector.o\
+      sternheimer_fd_preconditioner.o\
+      sternheimer_fd_projector_sampler.o\
+      sternheimer_fd_solver.o\
+      sternheimer_galerkin_audit.o\
+      sternheimer_galerkin_operator.o\
+      sternheimer_grid_diagnostics.o\
+      sternheimer_grid_transfer.o\
+      sternheimer_kq.o\
+      sternheimer_periodic_solver.o\
+      sternheimer_response_grid.o\
+      sternheimer_rpa.o\
+      sternheimer_runtime_options.o\
+      sternheimer_siab_memory.o\
+      sternheimer_siab_mpi.o\
+      sternheimer_siab_overlap.o\
+      sternheimer_siab_provenance.o\
+      sternheimer_siab_writer.o\
+      sternheimer_supercell_perturbation.o\
+      sternheimer_supercell_sector.o\
+      sternheimer_wavefunction_diagnostic.o\
+      sternheimer_weak_augmented.o\
+      sternheimer_weak_grid.o\
+      sternheimer_weak_preconditioner.o\
+      sternheimer_weak_q_reference.o\
 
 OBJS_BSE=esolver_lr_lcao_bse.o\
     hamilt_bse.o\

source/source_io/module_bessel/numerical_basis.cpp

@@ -15,6 +15,10 @@
 #include <cstring>
 #include <functional>
 #include <vector>
+#include <cmath>
+#include <memory>
+#include <stdexcept>
+#include <string>
 Numerical_Basis::Numerical_Basis()
 {
 }
@@ -55,6 +59,269 @@ Numerical_Basis::~Numerical_Basis()
 //     this->numerical_atomic_wfc(ik, wfcpw, psi, sf, ucell);
 // }
 
+int Numerical_Basis::siab_abacus_m(const int conventional_m)
+{
+    if (conventional_m == 0)
+    {
+        return 0;
+    }
+    return conventional_m > 0 ? 2 * conventional_m - 1 : -2 * conventional_m;
+}
+
+Numerical_Basis::SIABPrimitiveParameters
+Numerical_Basis::siab_parameters_from_input(const SIABPrimitiveInputs& settings, const int rcut_index, const int lmax)
+{
+    if (rcut_index < 0 || rcut_index >= static_cast<int>(settings.rcuts.size()))
+    {
+        throw std::out_of_range("SIAB primitive rcut index is out of range");
+    }
+    const double primitive_ecut_ry
+        = settings.ecut_spec == "default" ? settings.ecutwfc : std::stod(settings.ecut_spec);
+    if (lmax < -1)
+    {
+        throw std::invalid_argument("SIAB primitive lmax must be -1 or non-negative");
+    }
+    return SIABPrimitiveParameters{
+        primitive_ecut_ry,
+        settings.rcuts[rcut_index],
+        settings.smooth,
+        settings.sigma,
+        settings.tolerance,
+        lmax,
+    };
+}
+
+void Numerical_Basis::initialize_siab_basis(const UnitCell& ucell,
+                                            const SIABPrimitiveParameters& parameters)
+{
+    if (!std::isfinite(parameters.ecut_ry) || !std::isfinite(parameters.rcut_bohr)
+        || !std::isfinite(parameters.sigma) || !std::isfinite(parameters.tolerance)
+        || !(parameters.ecut_ry > 0.0) || !(parameters.rcut_bohr > 0.0)
+        || (parameters.smooth && !(parameters.sigma > 0.0)) || !(parameters.tolerance > 0.0))
+    {
+        throw std::invalid_argument("SIAB primitive parameters are invalid");
+    }
+
+    const int target_lmax = parameters.lmax < 0 ? ucell.lmax : parameters.lmax;
+    const std::vector<int> basis_shape = Numerical_Basis::siab_basis_shape_signature(ucell);
+    const bool same_parameters = this->initialized_siab_parameters.ecut_ry == parameters.ecut_ry
+                                 && this->initialized_siab_parameters.rcut_bohr == parameters.rcut_bohr
+                                 && this->initialized_siab_parameters.smooth == parameters.smooth
+                                 && this->initialized_siab_parameters.sigma == parameters.sigma
+                                 && this->initialized_siab_parameters.tolerance == parameters.tolerance
+                                 && this->initialized_siab_parameters.lmax == parameters.lmax;
+    if (this->init_label && same_parameters && this->initialized_siab_ntype == ucell.ntype
+        && this->initialized_siab_lmax == target_lmax && this->initialized_siab_nmax == ucell.nmax
+        && this->initialized_siab_basis_shape == basis_shape)
+    {
+        return;
+    }
+
+    this->bessel_basis.init(false,
+                            parameters.ecut_ry,
+                            ucell.ntype,
+                            target_lmax,
+                            parameters.smooth,
+                            parameters.sigma,
+                            parameters.rcut_bohr,
+                            parameters.tolerance,
+                            ucell);
+    this->mu_index = this->init_mu_index(ucell);
+    this->initialized_siab_parameters = parameters;
+    this->initialized_siab_ntype = ucell.ntype;
+    this->initialized_siab_lmax = target_lmax;
+    this->initialized_siab_nmax = ucell.nmax;
+    this->initialized_siab_basis_shape = basis_shape;
+    this->init_label = true;
+}
+
+std::vector<int> Numerical_Basis::siab_basis_shape_signature(const UnitCell& ucell)
+{
+    std::vector<int> signature;
+    for (int type = 0; type < ucell.ntype; ++type)
+    {
+        signature.push_back(ucell.atoms[type].na);
+        signature.push_back(ucell.atoms[type].nwl);
+        signature.push_back(static_cast<int>(ucell.atoms[type].l_nchi.size()));
+        signature.insert(signature.end(), ucell.atoms[type].l_nchi.begin(), ucell.atoms[type].l_nchi.end());
+    }
+    return signature;
+}
+
+void Numerical_Basis::fill_reciprocal_primitive(
+    const int l,
+    const int abacus_m,
+    const int primitive_index,
+    const std::complex<double>* structure_factor,
+    const ModuleBase::realArray& flq,
+    const ModuleBase::matrix& ylm,
+    const std::vector<double>& gpow,
+    const double normalization,
+    std::vector<std::complex<double>>& values)
+{
+    const int npw = static_cast<int>(gpow.size());
+    if (values.size() < static_cast<std::size_t>(npw))
+    {
+        throw std::invalid_argument("SIAB reciprocal primitive storage is too small");
+    }
+
+    const int lm = l * l + abacus_m;
+    const std::complex<double> lphase = normalization * std::pow(ModuleBase::IMAG_UNIT, -l);
+    for (int ig = 0; ig < npw; ++ig)
+    {
+        values[ig] = lphase * structure_factor[ig] * ylm(lm, ig) * flq(l, primitive_index, ig) * gpow[ig];
+    }
+}
+
+std::vector<Numerical_Basis::SIABPrimitiveReciprocalBlock> Numerical_Basis::siab_primitive_reciprocal_values(
+    const int ik,
+    const ModulePW::PW_Basis_K* wfcpw,
+    const Structure_Factor& sf,
+    const UnitCell& ucell,
+    const SIABPrimitiveParameters& parameters)
+{
+    if (wfcpw == nullptr)
+    {
+        throw std::invalid_argument("SIAB reciprocal primitives require a PW basis");
+    }
+    if (ik < 0 || ik >= wfcpw->nks)
+    {
+        throw std::out_of_range("SIAB primitive k-point index is out of range");
+    }
+    if (!(ucell.omega > 0.0))
+    {
+        throw std::invalid_argument("SIAB reciprocal primitives require a positive cell volume");
+    }
+
+    this->initialize_siab_basis(ucell, parameters);
+    const int target_lmax = parameters.lmax < 0 ? ucell.lmax : parameters.lmax;
+
+    const int npw = wfcpw->npwk[ik];
+    std::unique_ptr<ModuleBase::realArray> flq;
+    std::unique_ptr<ModuleBase::matrix> ylm;
+    std::vector<double> gpow;
+    if (npw > 0)
+    {
+        std::vector<ModuleBase::Vector3<double>> gk(npw);
+        for (int ig = 0; ig < npw; ++ig)
+        {
+            gk[ig] = wfcpw->getgpluskcar(ik, ig) * ucell.tpiba;
+        }
+        flq.reset(new ModuleBase::realArray(this->cal_flq(gk, target_lmax)));
+        ylm.reset(new ModuleBase::matrix(Numerical_Basis::cal_ylm(gk, target_lmax)));
+        gpow.assign(npw, 1.0);
+    }
+    const double normalization = 4.0 * ModuleBase::PI / std::sqrt(ucell.omega);
+    const int nprimitive = this->bessel_basis.get_ecut_number();
+
+    std::vector<SIABPrimitiveReciprocalBlock> blocks;
+    int offset = 0;
+    int global_atom = 0;
+    for (int type = 0; type < ucell.ntype; ++type)
+    {
+        for (int atom = 0; atom < ucell.atoms[type].na; ++atom, ++global_atom)
+        {
+            std::unique_ptr<std::complex<double>[]> sk(sf.get_sk(ik, type, atom, wfcpw));
+            for (int l = 0; l <= target_lmax; ++l)
+            {
+                for (int conventional_m = -l; conventional_m <= l; ++conventional_m)
+                {
+                    SIABPrimitiveReciprocalBlock block;
+                    block.type_index = type;
+                    block.element = ucell.atoms[type].label;
+                    block.atom_index = global_atom;
+                    block.l = l;
+                    block.m = conventional_m;
+                    block.n_primitive = nprimitive;
+                    block.offset = offset;
+                    block.values.reserve(nprimitive);
+
+                    const int abacus_m = Numerical_Basis::siab_abacus_m(conventional_m);
+                    for (int ie = 0; ie < nprimitive; ++ie)
+                    {
+                        std::vector<std::complex<double>> reciprocal(static_cast<std::size_t>(npw), ModuleBase::ZERO);
+                        if (npw > 0)
+                        {
+                            Numerical_Basis::fill_reciprocal_primitive(l,
+                                                                      abacus_m,
+                                                                      ie,
+                                                                      sk.get(),
+                                                                      *flq,
+                                                                      *ylm,
+                                                                      gpow,
+                                                                      normalization,
+                                                                      reciprocal);
+                        }
+
+                        block.values.push_back(std::move(reciprocal));
+                    }
+                    blocks.push_back(std::move(block));
+                    offset += nprimitive;
+                }
+            }
+        }
+    }
+    return blocks;
+}
+
+std::vector<Numerical_Basis::SIABPrimitiveGridBlock> Numerical_Basis::siab_primitive_grid_values(
+    const int ik,
+    const ModulePW::PW_Basis_K* wfcpw,
+    const Structure_Factor& sf,
+    const UnitCell& ucell,
+    const SIABPrimitiveParameters& parameters)
+{
+    if (wfcpw == nullptr)
+    {
+        throw std::invalid_argument("SIAB primitive grid requires a PW basis");
+    }
+    const auto reciprocal_blocks
+        = this->siab_primitive_reciprocal_values(ik, wfcpw, sf, ucell, parameters);
+    const double real_space_scale = 1.0 / std::sqrt(ucell.omega);
+    std::vector<SIABPrimitiveGridBlock> blocks;
+    blocks.reserve(reciprocal_blocks.size());
+    for (const SIABPrimitiveReciprocalBlock& reciprocal_block: reciprocal_blocks)
+    {
+        SIABPrimitiveGridBlock block;
+        block.type_index = reciprocal_block.type_index;
+        block.element = reciprocal_block.element;
+        block.atom_index = reciprocal_block.atom_index;
+        block.l = reciprocal_block.l;
+        block.m = reciprocal_block.m;
+        block.n_primitive = reciprocal_block.n_primitive;
+        block.offset = reciprocal_block.offset;
+        block.values.reserve(reciprocal_block.values.size());
+        for (const std::vector<std::complex<double>>& reciprocal: reciprocal_block.values)
+        {
+            std::vector<std::complex<double>> local_grid(static_cast<std::size_t>(wfcpw->nrxx), ModuleBase::ZERO);
+            if (wfcpw->gamma_only)
+            {
+                std::vector<double> local_real(static_cast<std::size_t>(std::max(1, wfcpw->nrxx)), 0.0);
+                wfcpw->recip2real(reciprocal.data(), local_real.data(), ik);
+                for (int ir = 0; ir < wfcpw->nrxx; ++ir)
+                {
+                    local_grid[static_cast<std::size_t>(ir)]
+                        = std::complex<double>(local_real[static_cast<std::size_t>(ir)] * real_space_scale, 0.0);
+                }
+            }
+            else
+            {
+                std::vector<std::complex<double>> local_complex(
+                    static_cast<std::size_t>(std::max(1, wfcpw->nrxx)), ModuleBase::ZERO);
+                wfcpw->recip2real(reciprocal.data(), local_complex.data(), ik);
+                for (int ir = 0; ir < wfcpw->nrxx; ++ir)
+                {
+                    local_grid[static_cast<std::size_t>(ir)]
+                        = local_complex[static_cast<std::size_t>(ir)] * real_space_scale;
+                }
+            }
+            block.values.push_back(std::move(local_grid));
+        }
+        blocks.push_back(std::move(block));
+    }
+    return blocks;
+}
+
 std::vector<ModuleBase::IntArray> Numerical_Basis::init_mu_index(const UnitCell& ucell)
 {
     GlobalV::ofs_running << " Initialize the mu index" << std::endl;