deepmodeling · dyzheng · Jan 8, 2024 · Jan 2, 2024 · Jan 4, 2024 · Jan 5, 2024
diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md
@@ -141,6 +141,7 @@
     - [out\_mat\_hs2](#out_mat_hs2)
     - [out\_mat\_t](#out_mat_t)
     - [out\_mat\_dh](#out_mat_dh)
+    - [out\_mat\_xc](#out_mat_xc)
     - [out\_app\_flag](#out_app_flag)
     - [out\_ndigits](#out_ndigits)
     - [out\_interval](#out_interval)
@@ -1570,6 +1571,13 @@ These variables are used to control the output of properties.
 - **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix (in Ry/Bohr). The format will be the same as the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `data-dHRx-sparse_SPIN0.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
 - **Default**: False
 
+### out_mat_xc
+
+- **Type**: Boolean
+- **Availability**: Numerical atomic orbital basis
+- **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in **Kohn-Sham orbital representation** (unit: Ry): $\braket{\psi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\psi_j}$ for each k point into files in the directory `OUT.${suffix}`, which is useful for the subsequent GW calculation. (Note that currently DeePKS term is not included. ) The files are named `k-$k-Vxc`, the meaning of `$k`corresponding to k point and spin  is same as [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs).
+- **Default**: False
+
 ### out_app_flag
 
 - **Type**: Boolean

diff --git a/source/Makefile.Objects b/source/Makefile.Objects
@@ -457,7 +457,6 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       unk_overlap_lcao.o\
       read_wfc_nao.o\
       write_wfc_nao.o\
-      write_HS.o\
       write_HS_sparse.o\
       single_R_io.o\
       write_HS_R.o\

diff --git a/source/module_base/global_variable.cpp b/source/module_base/global_variable.cpp
@@ -277,6 +277,8 @@ double nelec = 0;
 bool out_bandgap = false; // QO added for bandgap printing
 int out_interval = 1;    // convert from out_hsR_interval liuyu 2023-04-18
 
+bool out_mat_xc = false; // output Vxc in KS-wfc representation for GW calculation
+
 //==========================================================
 // Deltaspin related
 //==========================================================

diff --git a/source/module_base/global_variable.h b/source/module_base/global_variable.h
@@ -309,6 +309,8 @@ extern double nelec;
 extern bool out_bandgap;
 extern int out_interval;
 
+extern bool out_mat_xc; // output Vxc in KS-wfc representation for GW calculation
+
 // Deltaspin related
 extern bool sc_mag_switch; // 0: no deltaspin; 1: constrain atomic magnetic moments;
 extern bool decay_grad_switch; // 0: decay grad will be set to zero; 1: with decay grad set for some elements

diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -9,6 +9,7 @@
 #include "module_io/write_proj_band_lcao.h"
 #include "module_io/output_log.h"
 #include "module_io/to_qo.h"
+#include "module_io/write_Vxc.hpp"
 
 //--------------temporary----------------------------
 #include "module_base/global_function.h"
@@ -691,14 +692,11 @@ namespace ModuleESolver
             {
                 hamilt::MatrixBlock<TK> h_mat, s_mat;
                 this->p_hamilt->matrix(h_mat, s_mat);
-                ModuleIO::saving_HS(istep,
-                                    h_mat.p,
-                                    s_mat.p,
-                                    bit,
-                    hsolver::HSolverLCAO<TK>::out_mat_hs,
-                                    "data-" + std::to_string(ik),
-                                    this->LOWF.ParaV[0],
-                                    1); // LiuXh, 2017-03-21
+                if (hsolver::HSolverLCAO<TK>::out_mat_hs)
+                {
+                    ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK);
+                    ModuleIO::save_mat(istep, s_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK);
+                }
             }
         }
     }
@@ -829,6 +827,12 @@ namespace ModuleESolver
         }
     }
 
+    bool out_exc = true;    // tmp, add parameter!
+    if (GlobalV::out_mat_xc)
+        ModuleIO::write_Vxc<TK, TR>(GlobalV::NSPIN, GlobalV::NLOCAL, GlobalV::DRANK,
+            *this->psi, GlobalC::ucell, this->sf, *this->pw_rho, *this->pw_rhod, GlobalC::ppcell.vloc,
+            *this->pelec->charge, this->UHM, this->LM, this->LOC, this->kv, this->pelec->wg, GlobalC::GridD);
+
 #ifdef __EXX
     if (GlobalC::exx_info.info_global.cal_exx) // Peize Lin add if 2022.11.14
     {

diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -258,14 +258,11 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
             {
                 hamilt::MatrixBlock<complex<double>> h_mat, s_mat;
                 this->p_hamilt->matrix(h_mat, s_mat);
-                ModuleIO::saving_HS(istep,
-                                    h_mat.p,
-                                    s_mat.p,
-                                    bit,
-                    hsolver::HSolverLCAO<std::complex<double>>::out_mat_hs,
-                                    "data-" + std::to_string(ik),
-                                    this->LOWF.ParaV[0],
-                                    1); // LiuXh, 2017-03-21
+                if (hsolver::HSolverLCAO<std::complex<double>>::out_mat_hs)
+                {
+                    ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK);
+                    ModuleIO::save_mat(istep, h_mat.p, GlobalV::NLOCAL, bit, 1, GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), *this->LOWF.ParaV, GlobalV::DRANK);
+                }
             }
         }
     }

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp
@@ -210,14 +210,8 @@ void Force_LCAO_gamma::allocate_gamma(const Parallel_Orbitals& pv)
         this->UHM->genH
             .build_ST_new('S', cal_deri, GlobalC::ucell, this->UHM->genH.LM->Sloc.data(), INPUT.cal_syns, INPUT.dmax);
         bool bit = false; // LiuXh, 2017-03-21
-        ModuleIO::saving_HS(0,
-                            this->UHM->genH.LM->Hloc.data(),
-                            this->UHM->genH.LM->Sloc.data(),
-                            bit,
-                            1,
-                            "data-" + std::to_string(0),
-                            this->ParaV[0],
-                            0); // LiuXh, 2017-03-21
+        ModuleIO::save_mat(0, this->UHM->genH.LM->Hloc.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "H", "data-" + std::to_string(0), *this->ParaV, GlobalV::DRANK);
+        ModuleIO::save_mat(0, this->UHM->genH.LM->Sloc.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "S", "data-" + std::to_string(0), *this->ParaV, GlobalV::DRANK);
     }
 
     ModuleBase::timer::tick("Force_LCAO_gamma", "allocate_gamma");

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp
@@ -243,14 +243,8 @@ void Force_LCAO_k::allocate_k(const Parallel_Orbitals& pv,
             this->UHM->genH.LM->zeros_HSk('S');
             this->UHM->genH.LM->folding_fixedH(ik, kvec_d, 1);
             bool bit = false; // LiuXh, 2017-03-21
-            ModuleIO::saving_HS(0,
-                                this->UHM->genH.LM->Hloc2.data(),
-                                this->UHM->genH.LM->Sloc2.data(),
-                                bit,
-                                1,
-                                "data-" + std::to_string(ik),
-                                this->ParaV[0],
-                                0); // LiuXh, 2017-03-21
+            ModuleIO::save_mat(0, this->UHM->genH.LM->Hloc2.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "H", "data-" + std::to_string(ik), *this->ParaV, GlobalV::DRANK);
+            ModuleIO::save_mat(0, this->UHM->genH.LM->Sloc2.data(), GlobalV::NLOCAL, bit, 0, GlobalV::out_app_flag, "S", "data-" + std::to_string(ik), *this->ParaV, GlobalV::DRANK);
         }
     }
 

diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt
@@ -53,7 +53,6 @@ if(ENABLE_LCAO)
       read_dm.cpp
       read_wfc_nao.cpp
       write_wfc_nao.cpp
-      write_HS.cpp
       write_dm.cpp
       dos_nao.cpp
       output_dm.cpp

diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp
@@ -339,6 +339,7 @@ void Input::Default(void)
     out_band = 0;
     out_proj_band = 0;
     out_mat_hs = 0;
+    out_mat_xc = 0;
     cal_syns = 0;
     dmax = 0.01;
     out_mat_hs2 = 0; // LiuXh add 2019-07-15
@@ -1400,6 +1401,10 @@ bool Input::Read(const std::string& fn)
         {
             read_bool(ifs, out_mat_dh);
         }
+        else if (strcmp("out_mat_xc", word) == 0)
+        {
+            read_bool(ifs, out_mat_xc);
+        }
         else if (strcmp("out_interval", word) == 0)
         {
             read_value(ifs, out_interval);
@@ -3302,6 +3307,7 @@ void Input::Bcast()
     Parallel_Common::bcast_bool(out_mat_hs2); // LiuXh add 2019-07-15
     Parallel_Common::bcast_bool(out_mat_t);
     Parallel_Common::bcast_bool(out_mat_dh);
+    Parallel_Common::bcast_bool(out_mat_xc);
     Parallel_Common::bcast_bool(out_mat_r); // jingan add 2019-8-14
     Parallel_Common::bcast_int(out_wfc_lcao);
     Parallel_Common::bcast_bool(out_alllog);

diff --git a/source/module_io/input.h b/source/module_io/input.h
@@ -265,6 +265,7 @@ class Input
     bool out_band; // band calculation pengfei 2014-10-13
     bool out_proj_band; // projected band structure calculation jiyy add 2022-05-11
     bool out_mat_hs; // output H matrix and S matrix in local basis.
+    bool out_mat_xc; // output exchange-correlation matrix in KS-orbital representation.
     bool cal_syns; // calculate asynchronous S matrix to output
     double dmax; // maximum displacement of all atoms in one step (bohr)
     bool out_mat_hs2; // LiuXh add 2019-07-16, output H(R) matrix and S(R) matrix in local basis.

diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp
@@ -639,6 +639,7 @@ void Input_Conv::Convert(void)
     hsolver::HSolverLCAO<std::complex<double>>::out_mat_hsR = INPUT.out_mat_hs2; // LiuXh add 2019-07-16
     hsolver::HSolverLCAO<std::complex<double>>::out_mat_t = INPUT.out_mat_t;
     hsolver::HSolverLCAO<std::complex<double>>::out_mat_dh = INPUT.out_mat_dh;
+    GlobalV::out_mat_xc = INPUT.out_mat_xc;
     if (GlobalV::GAMMA_ONLY_LOCAL)
     {
         elecstate::ElecStateLCAO<double>::out_wfc_lcao = INPUT.out_wfc_lcao;

diff --git a/source/module_io/parameter_pool.cpp b/source/module_io/parameter_pool.cpp
@@ -927,6 +927,10 @@ bool input_parameters_set(std::map<std::string, InputParameter> input_parameters
     {
         INPUT.out_mat_hs = *static_cast<bool*>(input_parameters["out_mat_hs"].get());
     }
+    else if (input_parameters.count("out_mat_xc") != 0)
+    {
+        INPUT.out_mat_xc = *static_cast<bool*>(input_parameters["out_mat_xc"].get());
+    }
     else if (input_parameters.count("cal_syns") != 0)
     {
         INPUT.cal_syns = *static_cast<bool*>(input_parameters["cal_syns"].get());

diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp
@@ -153,6 +153,7 @@ TEST_F(InputConvTest, Conv)
     EXPECT_EQ(hsolver::HSolverLCAO<std::complex<double>>::out_mat_t, false);
     EXPECT_EQ(hsolver::HSolverLCAO<double>::out_mat_dh, INPUT.out_mat_dh);
     EXPECT_EQ(hsolver::HSolverLCAO<std::complex<double>>::out_mat_dh, INPUT.out_mat_dh);
+    EXPECT_EQ(GlobalV::out_mat_xc, false);
     EXPECT_EQ(GlobalV::out_interval, 1);
     EXPECT_EQ(elecstate::ElecStateLCAO<double>::out_wfc_lcao, false);
     EXPECT_EQ(berryphase::berry_phase_flag, false);

diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp
@@ -180,6 +180,7 @@ TEST_F(InputTest, Default)
         EXPECT_EQ(INPUT.out_proj_band,0);
         EXPECT_EQ(INPUT.out_mat_hs,0);
         EXPECT_EQ(INPUT.out_mat_hs2,0);
+        EXPECT_EQ(INPUT.out_mat_xc, 0);
         EXPECT_EQ(INPUT.out_interval,1);
         EXPECT_EQ(INPUT.out_app_flag,1);
         EXPECT_EQ(INPUT.out_mat_r,0);
@@ -542,6 +543,7 @@ TEST_F(InputTest, Read)
         EXPECT_EQ(INPUT.out_proj_band,0);
         EXPECT_EQ(INPUT.out_mat_hs,0);
         EXPECT_EQ(INPUT.out_mat_hs2,0);
+        EXPECT_EQ(INPUT.out_mat_xc, 0);
         EXPECT_EQ(INPUT.out_interval,1);
         EXPECT_EQ(INPUT.out_app_flag,0);
         EXPECT_EQ(INPUT.out_mat_r,0);

diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp
@@ -184,6 +184,7 @@ TEST_F(InputParaTest, Bcast)
     EXPECT_EQ(INPUT.out_proj_band, 0);
     EXPECT_EQ(INPUT.out_mat_hs, 0);
     EXPECT_EQ(INPUT.out_mat_hs2, 0);
+    EXPECT_EQ(INPUT.out_mat_xc, 0);
     EXPECT_EQ(INPUT.out_interval, 1);
     EXPECT_EQ(INPUT.out_app_flag, 1);
     EXPECT_EQ(INPUT.out_mat_r, 0);

diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp
@@ -323,6 +323,7 @@ TEST_F(write_input, LCAO5)
     EXPECT_THAT(output,
                 testing::HasSubstr("lcao_rmax                      30 #max R for 1D two-center integration table"));
     EXPECT_THAT(output, testing::HasSubstr("out_mat_hs                     0 #output H and S matrix"));
+    EXPECT_THAT(output, testing::HasSubstr("out_mat_xc                     0 #output exchange-correlation matrix in KS-orbital representation"));
     EXPECT_THAT(output, testing::HasSubstr("out_mat_hs2                    0 #output H(R) and S(R) matrix"));
     EXPECT_THAT(output, testing::HasSubstr("out_mat_dh                     0 #output of derivative of H(R) matrix"));
     EXPECT_THAT(