Feature: write Vxc(R) (#6061)

* fix the R range of data-HR

* feature: write Vxc(R)

* add parameter out_mat_xc2 and update tests

* fix compile

Author: maki49

docs/advanced/input_files/input-main.md

@@ -153,6 +153,7 @@
     - [out\_mat\_t](#out_mat_t)
     - [out\_mat\_dh](#out_mat_dh)
     - [out\_mat\_xc](#out_mat_xc)
+    - [out\_mat\_xc2](#out_mat_xc2)
     - [out\_eband\_terms](#out_eband_terms)
     - [out\_hr\_npz/out\_dm\_npz](#out_hr_npzout_dm_npz)
     - [dm\_to\_rho](#dm_to_rho)
@@ -1799,6 +1800,13 @@ These variables are used to control the output of properties.
 The band (KS orbital) energy for each (k-point, spin, band) will be printed in the file `OUT.${suffix}/vxc_out.dat`. If EXX is calculated, the local and EXX part of band energy will also be printed in `OUT.${suffix}/vxc_local_out.dat`and `OUT.${suffix}/vxc_exx_out.dat`, respectively. All the `vxc*_out.dat` files contains 3 integers (nk, nspin, nband) followed by nk\*nspin\*nband lines of energy Hartree and eV.
 - **Default**: False
 
+### out_mat_xc2
+
+- **Type**: Boolean
+- **Availability**: Numerical atomic orbital (NAO) basis
+- **Description**: Whether to print the exchange-correlation matrices in **numerical orbital representation** (unit: Ry): $\braket{\phi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\phi_j}(\mathbf{R})$ in CSR format (the same format as [out_mat_hs2](../elec_properties/hs_matrix.md#out_mat_hs2)) in the directory `OUT.${suffix}`. (Note that currently DeePKS term is not included. ) The files are named `Vxc_R_spin$s`.
+- **Default**: False
+
 ### out_eband_terms
 
 - **Type**: Boolean

source/module_esolver/esolver_ks_lcao.cpp

@@ -33,6 +33,7 @@
 //be careful of hpp, there may be multiple definitions of functions, 20250302, mohan
 #include "module_io/write_eband_terms.hpp"
 #include "module_io/write_vxc.hpp"
+#include "module_io/write_vxc_r.hpp"
 #include "module_hamilt_lcao/hamilt_lcaodft/hs_matrix_k.hpp"
 
 //--------------temporary----------------------------
@@ -488,27 +489,49 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
     if (PARAM.inp.out_mat_xc)
     {
         ModuleIO::write_Vxc<TK, TR>(PARAM.inp.nspin,
-                                    PARAM.globalv.nlocal,
-                                    GlobalV::DRANK,
-                                    &this->pv,
-                                    *this->psi,
-                                    ucell,
-                                    this->sf,
-                                    this->solvent,
-                                    *this->pw_rho,
-                                    *this->pw_rhod,
-                                    this->locpp.vloc,
-                                    this->chr,
-                                    this->GG,
-                                    this->GK,
-                                    this->kv,
-                                    orb_.cutoffs(),
-                                    this->pelec->wg,
-                                    this->gd
+            PARAM.globalv.nlocal,
+            GlobalV::DRANK,
+            &this->pv,
+            *this->psi,
+            ucell,
+            this->sf,
+            this->solvent,
+            *this->pw_rho,
+            *this->pw_rhod,
+            this->locpp.vloc,
+            this->chr,
+            this->GG,
+            this->GK,
+            this->kv,
+            orb_.cutoffs(),
+            this->pelec->wg,
+            this->gd
 #ifdef __EXX
-                                    ,
-                                    this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr,
-                                    this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr
+            ,
+            this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr,
+            this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr
+#endif
+        );
+    }
+    if (PARAM.inp.out_mat_xc2)
+    {
+        ModuleIO::write_Vxc_R<TK, TR>(PARAM.inp.nspin,
+            &this->pv,
+            ucell,
+            this->sf,
+            this->solvent,
+            *this->pw_rho,
+            *this->pw_rhod,
+            this->locpp.vloc,
+            this->chr,
+            this->GG,
+            this->GK,
+            this->kv,
+            orb_.cutoffs(),
+            this->gd
+#ifdef __EXX
+            , this->exx_lri_double ? &this->exx_lri_double->Hexxs : nullptr,
+            this->exx_lri_complex ? &this->exx_lri_complex->Hexxs : nullptr
 #endif
         );
     }

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.cpp

@@ -4,6 +4,22 @@
 
 namespace hamilt
 {
+    RI::Cell_Nearest<int, int, 3, double, 3> init_cell_nearest(const UnitCell& ucell, const std::array<int, 3>& Rs_period)
+    {
+        RI::Cell_Nearest<int, int, 3, double, 3> cell_nearest;
+        std::map<int, std::array<double, 3>> atoms_pos;
+        for (int iat = 0; iat < ucell.nat; ++iat) {
+            atoms_pos[iat] = RI_Util::Vector3_to_array3(
+                ucell.atoms[ucell.iat2it[iat]]
+                .tau[ucell.iat2ia[iat]]);
+        }
+        const std::array<std::array<double, 3>, 3> latvec
+            = { RI_Util::Vector3_to_array3(ucell.a1),
+               RI_Util::Vector3_to_array3(ucell.a2),
+               RI_Util::Vector3_to_array3(ucell.a3) };
+        cell_nearest.init(atoms_pos, latvec, Rs_period);
+        return cell_nearest;
+    }
 
 template<>
 void OperatorEXX<OperatorLCAO<double, double>>::add_loaded_Hexx(const int ik)

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h

@@ -73,6 +73,22 @@ class OperatorEXX<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
       std::vector<std::vector<std::complex<double>>> Hexxc_k_load;
 };
 
+using TAC = std::pair<int, std::array<int, 3>>;
+
+RI::Cell_Nearest<int, int, 3, double, 3> init_cell_nearest(const UnitCell& ucell, const std::array<int, 3>& Rs_period);
+
+// allocate according to the read-in HexxR, used in nscf
+template <typename Tdata, typename TR>
+void reallocate_hcontainer(const std::vector<std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& Hexxs,
+    HContainer<TR>* hR,
+    const RI::Cell_Nearest<int, int, 3, double, 3>* const cell_nearest = nullptr);
+
+/// allocate according to BvK cells, used in scf
+template <typename TR>
+void reallocate_hcontainer(const int nat, HContainer<TR>* hR,
+    const std::array<int, 3>& Rs_period,
+    const RI::Cell_Nearest<int, int, 3, double, 3>* const cell_nearest = nullptr);
+
 } // namespace hamilt
 #endif // __EXX
 #include "op_exx_lcao.hpp"

source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.hpp

@@ -12,10 +12,12 @@
 namespace hamilt
 {
     using TAC = std::pair<int, std::array<int, 3>>;
+
     // allocate according to the read-in HexxR, used in nscf
     template <typename Tdata, typename TR>
-    inline void reallocate_hcontainer(const std::vector<std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& Hexxs,
-        HContainer<TR>* hR)
+    void reallocate_hcontainer(const std::vector<std::map<int, std::map<TAC, RI::Tensor<Tdata>>>>& Hexxs,
+        HContainer<TR>* hR,
+        const RI::Cell_Nearest<int, int, 3, double, 3>* const cell_nearest)
     {
         auto* pv = hR->get_paraV();
         bool need_allocate = false;
@@ -27,7 +29,10 @@ namespace hamilt
                 const int& iat1 = Htmp2.first.first;
                 if (pv->get_row_size(iat0) > 0 && pv->get_col_size(iat1) > 0)
                 {
-                    const Abfs::Vector3_Order<int>& R = RI_Util::array3_to_Vector3(Htmp2.first.second);
+                    const Abfs::Vector3_Order<int>& R = RI_Util::array3_to_Vector3(
+                        (cell_nearest ?
+                            cell_nearest->get_cell_nearest_discrete(iat0, iat1, Htmp2.first.second)
+                            : Htmp2.first.second));
                     BaseMatrix<TR>* HlocR = hR->find_matrix(iat0, iat1, R.x, R.y, R.z);
                     if (HlocR == nullptr)
                     { // add R to HContainer
@@ -40,11 +45,12 @@ namespace hamilt
         }
         if (need_allocate) { hR->allocate(nullptr, true); }
     }
+
     /// allocate according to BvK cells, used in scf
     template <typename TR>
-    inline void reallocate_hcontainer(const int nat, HContainer<TR>* hR,
+    void reallocate_hcontainer(const int nat, HContainer<TR>* hR,
         const std::array<int, 3>& Rs_period,
-        const RI::Cell_Nearest<int, int, 3, double, 3>* const cell_nearest = nullptr)
+        const RI::Cell_Nearest<int, int, 3, double, 3>* const cell_nearest)
     {
         auto* pv = hR->get_paraV();
         auto Rs = RI_Util::get_Born_von_Karmen_cells(Rs_period);
@@ -154,18 +160,7 @@ OperatorEXX<OperatorLCAO<TK, TR>>::OperatorEXX(HS_Matrix_K<TK>* hsk_in,
             const std::array<int, 3> Rs_period = { this->kv.nmp[0], this->kv.nmp[1], this->kv.nmp[2] };
             if (this->use_cell_nearest)
             {
-                // set cell_nearest
-                std::map<int, std::array<double, 3>> atoms_pos;
-                for (int iat = 0; iat < ucell.nat; ++iat) {
-                    atoms_pos[iat] = RI_Util::Vector3_to_array3(
-                        ucell.atoms[ucell.iat2it[iat]]
-                        .tau[ucell.iat2ia[iat]]);
-                }
-                const std::array<std::array<double, 3>, 3> latvec
-                    = { RI_Util::Vector3_to_array3(ucell.a1),
-                       RI_Util::Vector3_to_array3(ucell.a2),
-                       RI_Util::Vector3_to_array3(ucell.a3) };
-                this->cell_nearest.init(atoms_pos, latvec, Rs_period);
+                this->cell_nearest = init_cell_nearest(ucell, Rs_period);
                 reallocate_hcontainer(ucell.nat, this->hR, Rs_period, &this->cell_nearest);
             }
             else { reallocate_hcontainer(ucell.nat, this->hR, Rs_period); }

source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp

@@ -23,7 +23,7 @@ void sparse_format::cal_HSR(const UnitCell& ucell,
 ) {
     ModuleBase::TITLE("sparse_format", "cal_HSR");
 
-    sparse_format::set_R_range(HS_Arrays.all_R_coor, grid);
+    // sparse_format::set_R_range(HS_Arrays.all_R_coor, grid);
 
     const int nspin = PARAM.inp.nspin;
 
@@ -33,6 +33,8 @@ void sparse_format::cal_HSR(const UnitCell& ucell,
             = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>, double>*>(
                 p_ham);
 
+        HS_Arrays.all_R_coor = get_R_range(*(p_ham_lcao->getHR()));
+
         if (TD_Velocity::tddft_velocity) {
             sparse_format::cal_HContainer_td(
                 pv,
@@ -57,7 +59,9 @@ void sparse_format::cal_HSR(const UnitCell& ucell,
         hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>*
             p_ham_lcao
             = dynamic_cast<hamilt::HamiltLCAO<std::complex<double>,
-                                              std::complex<double>>*>(p_ham);
+            std::complex<double>>*>(p_ham);
+
+        HS_Arrays.all_R_coor = get_R_range(*(p_ham_lcao->getHR()));
 
         sparse_format::cal_HContainer_cd(pv,
                                          current_spin,

source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h

@@ -4,6 +4,23 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
 
 namespace sparse_format {
+    template<typename T>
+    std::set<Abfs::Vector3_Order<int>> get_R_range(const hamilt::HContainer<T>& hR)
+    {
+        std::set<Abfs::Vector3_Order<int>> all_R_coor;
+        for (int iap = 0;iap < hR.size_atom_pairs(); ++iap)
+        {
+            const hamilt::AtomPair<T>& atom_pair = hR.get_atom_pair(iap);
+            for (int iR = 0; iR < atom_pair.get_R_size(); ++iR)
+            {
+                const auto& r_index = atom_pair.get_R_index(iR);
+                Abfs::Vector3_Order<int> dR(r_index.x, r_index.y, r_index.z);
+                all_R_coor.insert(dR);
+            }
+        }
+        return all_R_coor;
+    };
+
     using TAC = std::pair<int, std::array<int, 3>>;
     void cal_HSR(const UnitCell& ucell,
                  const Parallel_Orbitals& pv,

source/module_io/read_input_item_output.cpp

@@ -298,6 +298,12 @@ void ReadInput::item_output()
         read_sync_bool(input.out_mat_xc);
         this->add_item(item);
     }
+    {
+        Input_Item item("out_mat_xc2");
+        item.annotation = "output exchange-correlation matrix in NAO representation";
+        read_sync_bool(input.out_mat_xc2);
+        this->add_item(item);
+    }
     {
         Input_Item item("out_eband_terms");
         item.annotation = "output the band energy terms separately";

source/module_io/test/read_input_ptest.cpp

@@ -208,6 +208,7 @@ TEST_F(InputParaTest, ParaRead)
     EXPECT_EQ(param.inp.out_mat_hs[1], 8);
     EXPECT_EQ(param.inp.out_mat_hs2, 0);
     EXPECT_FALSE(param.inp.out_mat_xc);
+    EXPECT_FALSE(param.inp.out_mat_xc2);
     EXPECT_FALSE(param.inp.out_eband_terms);
     EXPECT_EQ(param.inp.out_interval, 1);
     EXPECT_EQ(param.inp.out_app_flag, 0);

source/module_io/write_vxc.hpp

@@ -30,16 +30,6 @@ struct TGint<std::complex<double>>
 namespace ModuleIO
 {
 
-inline void gint_vl(Gint_Gamma& gg, Gint_inout& io)
-{
-    gg.cal_vlocal(&io, false);
-};
-
-inline void gint_vl(Gint_k& gk, Gint_inout& io, ModuleBase::matrix& wg)
-{
-    gk.cal_gint(&io);
-};
-
 inline void set_para2d_MO(const Parallel_Orbitals& pv, const int nbands, Parallel_2D& p2d)
 {
     std::ofstream ofs;
@@ -134,6 +124,8 @@ std::vector<double> orbital_energy(const int ik, const int nbands, const std::ve
     return e;
 }
 
+#ifndef SET_GINT_POINTER_H
+#define SET_GINT_POINTER_H
 // mohan update 2024-04-01
 template <typename T>
 void set_gint_pointer(Gint_Gamma& gint_gamma, Gint_k& gint_k, typename TGint<T>::type*& gint);
@@ -153,6 +145,7 @@ void set_gint_pointer<std::complex<double>>(Gint_Gamma& gint_gamma,
 {
     gint = &gint_k;
 }
+#endif
 
 inline void write_orb_energy(const K_Vectors& kv,
     const int nspin0, const int nbands,
@@ -223,7 +216,7 @@ void write_Vxc(const int nspin,
     // 2. allocate AO-matrix
     // R (the number of hR: 1 for nspin=1, 4; 2 for nspin=2)
     int nspin0 = (nspin == 2) ? 2 : 1;
-    std::vector<hamilt::HContainer<TR>> vxcs_R_ao(nspin0, hamilt::HContainer<TR>(pv));
+    std::vector<hamilt::HContainer<TR>> vxcs_R_ao(nspin0, hamilt::HContainer<TR>(ucell, pv));
     for (int is = 0; is < nspin0; ++is) {
         vxcs_R_ao[is].set_zero();
         if (std::is_same<TK, double>::value) { vxcs_R_ao[is].fix_gamma(); }