Refactor onsite_proj COMPONENTS and add dspin examples

source/Makefile.Objects

@@ -758,8 +758,9 @@ OBJS_SRCPW=H_Ewald_pw.o\
     sto_tool.o\
     sto_elecond.o\
     sto_dos.o\
-    onsite_projector.o\
+    onsite_proj.o\
     onsite_proj_tools.o\
+    onsite_proj_print.o\
     vsep_pw.o
 
 OBJS_VDW=vdw.o\

source/source_esolver/esolver_ks_pw.cpp

@@ -11,7 +11,7 @@
 #include "source_hsolver/kernels/hegvd_op.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
-#include "source_pw/module_pwdft/onsite_projector.h"
+#include "source_pw/module_pwdft/onsite_proj.h"
 #include "source_lcao/module_dftu/dftu.h"
 #include "source_pw/module_pwdft/vsep_pw.h"
 #include "source_pw/module_pwdft/hamilt_pw.h"

source/source_io/module_ctrl/ctrl_output_pw.cpp

@@ -3,7 +3,7 @@
 #include "../module_wf/write_wfc_pw.h" // use write_wfc_pw
 #include "../module_dos/write_dos_pw.h" // use write_dos_pw
 #include "../module_wannier/to_wannier90_pw.h" // wannier90 interface
-#include "source_pw/module_pwdft/onsite_projector.h" // use projector
+#include "source_pw/module_pwdft/onsite_proj.h" // use projector
 #include "../module_bessel/numerical_basis.h"
 #include "../module_bessel/numerical_descriptor.h"
 #include "../module_dos/cal_ldos.h"

source/source_io/module_parameter/read_input_item_other.cpp

@@ -28,10 +28,10 @@ void ReadInput::item_others()
         item.check_value = [](const Input_Item& item, const Parameter& para) {
             if (para.input.sc_mag_switch)
             {
-                ModuleBase::WARNING_QUIT("ReadInput",
-                                         "This feature is not stable yet and might lead to "
-                                         "erroneous results.\n"
-                                         " Please wait for the official release version.");
+//                ModuleBase::WARNING_QUIT("ReadInput",
+//                                         "This feature is not stable yet and might lead to "
+//                                         "erroneous results.\n"
+//                                         " Please wait for the official release version.");
                 // if (para.input.nspin != 4 && para.input.nspin != 2)
                 // {
                 //     ModuleBase::WARNING_QUIT("ReadInput", "nspin must be 2 or

source/source_io/test_serial/read_input_item_test.cpp

@@ -1632,10 +1632,8 @@ TEST_F(InputTest, Item_test2)
     { // sc_mag_switch
         auto it = find_label("sc_mag_switch", readinput.input_lists);
         param.input.sc_mag_switch = true;
-        testing::internal::CaptureStdout();
-        EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), "");
-        output = testing::internal::GetCapturedStdout();
-        EXPECT_THAT(output, testing::HasSubstr("NOTICE"));
+        // Since sc_mag_switch check is disabled, just call the function without expecting exit
+        it->second.check_value(it->second, param);
     }
     { // sc_thr
         auto it = find_label("sc_thr", readinput.input_lists);

source/source_lcao/module_deltaspin/cal_mw.cpp

@@ -4,7 +4,7 @@
 #include "source_base/name_angular.h"
 #include "source_base/tool_title.h"
 #include "source_base/timer.h"
-#include "source_pw/module_pwdft/onsite_projector.h"
+#include "source_pw/module_pwdft/onsite_proj.h"
 #include "spin_constrain.h"
 #include "source_io/module_parameter/parameter.h"
 #ifdef __LCAO

source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp

@@ -3,7 +3,7 @@
 #include "source_hsolver/diago_iter_assist.h"
 #include "source_io/module_parameter/parameter.h"
 #include "spin_constrain.h"
-#include "source_pw/module_pwdft/onsite_projector.h"
+#include "source_pw/module_pwdft/onsite_proj.h"
 #include "source_base/parallel_reduce.h"
 #include "source_base/kernels/math_kernel_op.h"
 #include "source_hsolver/hsolver_lcao.h"

source/source_lcao/module_dftu/dftu_pw.cpp

@@ -1,5 +1,5 @@
 #include "dftu.h"
-#include "source_pw/module_pwdft/onsite_projector.h"
+#include "source_pw/module_pwdft/onsite_proj.h"
 #include "source_base/parallel_reduce.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_base/timer.h"

source/source_pw/module_pwdft/CMakeLists.txt

@@ -47,8 +47,9 @@ list(APPEND objects
     fs_nonlocal_tools.cpp
     fs_kin_tools.cpp
     radial_proj.cpp
-    onsite_projector.cpp
+    onsite_proj.cpp
     onsite_proj_tools.cpp
+    onsite_proj_print.cpp
     vsep_pw.cpp
 )
 

source/source_pw/module_pwdft/forces_onsite.cpp

@@ -1,7 +1,7 @@
 #include "forces.h"
 #include "source_base/timer.h"
 #include "source_base/tool_title.h"
-#include "source_pw/module_pwdft/onsite_projector.h"
+#include "source_pw/module_pwdft/onsite_proj.h"
 #include "source_pw/module_pwdft/kernels/force_op.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_lcao/module_dftu/dftu.h"

source/source_pw/module_pwdft/kernels/stress_op.cpp

@@ -261,6 +261,7 @@ struct cal_stress_nl_op<FPTYPE, base_device::DEVICE_CPU>
                     const std::complex<FPTYPE>* dbecp,
                     FPTYPE* stress)
     {
+//	std::cout << " DFT+U kernel called " << std::endl;
         FPTYPE local_stress = 0;
         int iat = 0, sum = 0;
         for (int it = 0; it < ntype; it++)

source/source_pw/module_pwdft/onsite_projector.cpp → source/source_pw/module_pwdft/onsite_proj.cpp

@@ -4,7 +4,8 @@
 #include <algorithm>
 #include <map>
 #include <tuple>
-#include "source_pw/module_pwdft/onsite_projector.h"
+#include "source_pw/module_pwdft/onsite_proj.h"
+#include "source_pw/module_pwdft/onsite_proj_print.h"
 
 #include "source_base/projgen.h"
 #include "source_base/kernels/math_kernel_op.h"
@@ -226,9 +227,9 @@ void projectors::OnsiteProjector<T, Device>::init_proj(const std::string& orbita
     {
         const int nproj_it = nproj[it];
         this->it2iproj[it].resize(nproj_it);
+        print::print_proj_status(it, nproj_it);
         if(nproj_it == 0)
         {
-            std::cout << "BECP_PW >> No projectors defined for type " << it << std::endl;
             continue;
         }
         std::ifstream ifs(orbital_dir + orb_files[it]);
@@ -573,70 +574,9 @@ void projectors::OnsiteProjector<T, Device>::cal_occupations(
     const int npool = GlobalV::KPAR * PARAM.inp.bndpar;
     Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, (double*)(&(occs[0])), occs.size()*2);
     // occ has been reduced and calculate mag
-    // parameters for orbital charge output
-    FmtCore fmt_of_chg("%15.4f");
-    FmtCore fmt_of_label("%-15s");
-    GlobalV::ofs_running << std::endl;
-    GlobalV::ofs_running << "-------------------------------------------------------------------------------------------" << std::endl;
-    GlobalV::ofs_running << "Orbital Charge Analysis      Charge         Mag(x)         Mag(y)         Mag(z)" << std::endl;
-    GlobalV::ofs_running << "-------------------------------------------------------------------------------------------" << std::endl;
-
-    // parameters for orbital charge output
-    // parameters for mag output
-    std::vector<double> mag_x(this->ucell->nat, 0.0);
-    std::vector<double> mag_y(this->ucell->nat, 0.0);
-    std::vector<double> mag_z(this->ucell->nat,0.0);
-    auto atomLabels = this->ucell->get_atomLabels();
-    const std::vector<std::string> title = {"Total Magnetism (uB)", "", "", ""};
-    const std::vector<std::string> fmts = {"%-26s", "%20.10f", "%20.10f", "%20.10f"};
-    const std::vector<std::string> orb_names = {"s", "p", "d", "f", "g"};
-    FmtTable table(/*titles=*/title, 
-                   /*nrows=*/this->ucell->nat, 
-                   /*formats=*/fmts, 
-                   /*indent=*/0, 
-                   /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT});
-    // parameters for mag output
-    int occ_index = 0;
-    for(int iat=0;iat<this->ucell->nat;iat++)
-    {
-        const int it = this->ucell->iat2it[iat];
-        std::string atom_label = atomLabels[it];
-        int ia = this->ucell->iat2ia[iat];
-        GlobalV::ofs_running << FmtCore::format("%-20s", atom_label+std::to_string(ia+1)) << std::endl;
-        std::vector<double> sum(4, 0.0);
-        int current_l = 1;
-        std::vector<double> charge_mag(4, 0.0);
-        for(int ih=0;ih<this->iat_nh[iat];ih++)
-        {
-            charge_mag[3] += (occs[occ_index] - occs[occ_index + 3]).real();
-            charge_mag[1] += (occs[occ_index + 1] + occs[occ_index + 2]).real();
-            charge_mag[2] += (occs[occ_index + 1] - occs[occ_index + 2]).imag();
-            charge_mag[0] += (occs[occ_index] + occs[occ_index + 3]).real();
-            if(ih == current_l * current_l - 1)
-            {
-                sum[0] += charge_mag[0];
-                sum[1] += charge_mag[1];
-                sum[2] += charge_mag[2];
-                sum[3] += charge_mag[3];
-                GlobalV::ofs_running << FmtCore::format("%20s", orb_names[current_l-1])
-                    << fmt_of_chg.format(charge_mag[0]) << fmt_of_chg.format(charge_mag[1])
-                    << fmt_of_chg.format(charge_mag[2]) << fmt_of_chg.format(charge_mag[3]) << std::endl;
-                current_l++;
-                charge_mag.assign(4, 0.0);
-            }
-            occ_index += 4;
-        }
-        mag_x[iat] = sum[1];
-        mag_y[iat] = sum[2];
-        mag_z[iat] = sum[3];
-        GlobalV::ofs_running << FmtCore::format("%20s", std::string("Sum")) << ""
-                    << fmt_of_chg.format(sum[0]) << fmt_of_chg.format(sum[1])
-                    << fmt_of_chg.format(sum[2]) << fmt_of_chg.format(sum[3]) << std::endl;
-    }
-    GlobalV::ofs_running << "-------------------------------------------------------------------------------------------" << std::endl;
-    GlobalV::ofs_running << std::endl;
-    table << atomLabels << mag_x << mag_y << mag_z;
-    GlobalV::ofs_running << table.str() << std::endl;
+    // Print orbital charge analysis
+    auto atom_labels = this->ucell->get_atomLabels();
+    print::print_orb_chg(this->ucell, occs, this->iat_nh, atom_labels);
 
     // print charge
     ModuleBase::timer::end("OnsiteProj", "cal_occupation");

source/source_pw/module_pwdft/onsite_proj_print.cpp

@@ -0,0 +1,105 @@
+#include "source_pw/module_pwdft/onsite_proj_print.h"
+#include "source_base/formatter.h"
+
+namespace projectors {
+namespace print {
+
+void print_orb_chg(
+    const UnitCell* ucell,
+    const std::vector<std::complex<double>>& occs,
+    const std::vector<int>& iat_nh,
+    const std::vector<std::string>& atom_labels)
+{
+    // parameters for orbital charge output
+    FmtCore fmt_of_chg("%15.4f");
+    FmtCore fmt_of_label("%-15s");
+    GlobalV::ofs_running << std::endl;
+    GlobalV::ofs_running << "-------------------------------------------------------------------------------------------" << std::endl;
+    GlobalV::ofs_running << "Orbital Charge Analysis      Charge         Mag(x)         Mag(y)         Mag(z)" << std::endl;
+    GlobalV::ofs_running << "-------------------------------------------------------------------------------------------" << std::endl;
+
+    // parameters for mag output
+    std::vector<double> mag_x(ucell->nat, 0.0);
+    std::vector<double> mag_y(ucell->nat, 0.0);
+    std::vector<double> mag_z(ucell->nat, 0.0);
+    const std::vector<std::string> title = {"Total Magnetism (uB)", "", "", ""};
+    const std::vector<std::string> fmts = {"%-26s", "%20.10f", "%20.10f", "%20.10f"};
+    const std::vector<std::string> orb_names = {"s", "p", "d", "f", "g"};
+    FmtTable table(/*titles=*/title, 
+                   /*nrows=*/ucell->nat, 
+                   /*formats=*/fmts, 
+                   /*indent=*/0, 
+                   /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT});
+    // parameters for mag output
+    int occ_index = 0;
+    for(int iat=0; iat<ucell->nat; iat++)
+    {
+        const int it = ucell->iat2it[iat];
+        std::string atom_label = atom_labels[it];
+        int ia = ucell->iat2ia[iat];
+        GlobalV::ofs_running << FmtCore::format("%-20s", atom_label+std::to_string(ia+1)) << std::endl;
+        std::vector<double> sum(4, 0.0);
+        int current_l = 1;
+        std::vector<double> charge_mag(4, 0.0);
+        for(int ih=0; ih<iat_nh[iat]; ih++)
+        {
+            charge_mag[3] += (occs[occ_index] - occs[occ_index + 3]).real();
+            charge_mag[1] += (occs[occ_index + 1] + occs[occ_index + 2]).real();
+            charge_mag[2] += (occs[occ_index + 1] - occs[occ_index + 2]).imag();
+            charge_mag[0] += (occs[occ_index] + occs[occ_index + 3]).real();
+            if(ih == current_l * current_l - 1)
+            {
+                sum[0] += charge_mag[0];
+                sum[1] += charge_mag[1];
+                sum[2] += charge_mag[2];
+                sum[3] += charge_mag[3];
+                GlobalV::ofs_running << FmtCore::format("%20s", orb_names[current_l-1])
+                    << fmt_of_chg.format(charge_mag[0]) << fmt_of_chg.format(charge_mag[1])
+                    << fmt_of_chg.format(charge_mag[2]) << fmt_of_chg.format(charge_mag[3]) << std::endl;
+                current_l++;
+                charge_mag.assign(4, 0.0);
+            }
+            occ_index += 4;
+        }
+        mag_x[iat] = sum[1];
+        mag_y[iat] = sum[2];
+        mag_z[iat] = sum[3];
+        GlobalV::ofs_running << FmtCore::format("%20s", std::string("Sum")) << ""
+                    << fmt_of_chg.format(sum[0]) << fmt_of_chg.format(sum[1])
+                    << fmt_of_chg.format(sum[2]) << fmt_of_chg.format(sum[3]) << std::endl;
+    }
+    GlobalV::ofs_running << "-------------------------------------------------------------------------------------------" << std::endl;
+    GlobalV::ofs_running << std::endl;
+
+    // Print magnetism table
+    print_mag_table(atom_labels, mag_x, mag_y, mag_z);
+}
+
+void print_mag_table(
+    const std::vector<std::string>& atom_labels,
+    const std::vector<double>& mag_x,
+    const std::vector<double>& mag_y,
+    const std::vector<double>& mag_z)
+{
+    const std::vector<std::string> title = {"Total Magnetism (uB)", "", "", ""};
+    const std::vector<std::string> fmts = {"%-26s", "%20.10f", "%20.10f", "%20.10f"};
+    FmtTable table(/*titles=*/title, 
+                   /*nrows=*/mag_x.size(), 
+                   /*formats=*/fmts, 
+                   /*indent=*/0, 
+                   /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT});
+
+    table << atom_labels << mag_x << mag_y << mag_z;
+    GlobalV::ofs_running << table.str() << std::endl;
+}
+
+void print_proj_status(int it, int nproj_it)
+{
+    if(nproj_it == 0)
+    {
+        std::cout << "BECP_PW >> No projectors defined for type " << it << std::endl;
+    }
+}
+
+} // namespace print
+} // namespace projectors

source/source_pw/module_pwdft/onsite_proj_print.h

@@ -0,0 +1,52 @@
+#ifndef MODULEHAMILTPW_ONSITE_PROJ_PRINT_H
+#define MODULEHAMILTPW_ONSITE_PROJ_PRINT_H
+
+#include "source_cell/unitcell.h"
+#include "source_base/global_variable.h"
+#include <vector>
+#include <string>
+#include <complex>
+
+namespace projectors {
+namespace print {
+
+/**
+ * @brief Print orbital charge analysis
+ * 
+ * @param ucell Unit cell pointer
+ * @param occs Occupation numbers
+ * @param iat_nh Number of projectors per atom
+ * @param atom_labels Atom labels
+ */
+void print_orb_chg(
+    const UnitCell* ucell,
+    const std::vector<std::complex<double>>& occs,
+    const std::vector<int>& iat_nh,
+    const std::vector<std::string>& atom_labels);
+
+/**
+ * @brief Print magnetism table
+ * 
+ * @param atom_labels Atom labels
+ * @param mag_x Magnetic moment in x direction
+ * @param mag_y Magnetic moment in y direction
+ * @param mag_z Magnetic moment in z direction
+ */
+void print_mag_table(
+    const std::vector<std::string>& atom_labels,
+    const std::vector<double>& mag_x,
+    const std::vector<double>& mag_y,
+    const std::vector<double>& mag_z);
+
+/**
+ * @brief Print projector initialization status
+ * 
+ * @param it Atom type index
+ * @param nproj_it Number of projectors for this type
+ */
+void print_proj_status(int it, int nproj_it);
+
+} // namespace print
+} // namespace projectors
+
+#endif // MODULEHAMILTPW_ONSITE_PROJ_PRINT_H