Feature: DFT+U for noncollinear spin calculation with PW base

CMakeLists.txt

@@ -711,17 +711,17 @@ target_link_libraries(
   esolver
   vdw
   device
-  container)
+  container
+  dftu
+  deltaspin)
 if(ENABLE_LCAO)
   target_link_libraries(
     ${ABACUS_BIN_NAME}
     hamilt_lcao
     tddft
     orb
     gint
-    dftu
     hcontainer
-    deltaspin
     numerical_atomic_orbitals
     lr
     rdmft)

docs/advanced/input_files/input-main.md

@@ -1235,6 +1235,12 @@ Note: In new angle mixing, you should set `mixing_beta_mag >> mixing_beta`. The
 - **Description**: To determine the number of old iterations' `drho` used in slope calculations.
 - **Default**: `mixing_ndim`
 
+### sc_os_ndim
+
+- **Type**: int
+- **Description**: To determine the number of old iterations to judge oscillation, it occured,  more accurate lambda with DeltaSpin method would be calculated, only for PW base.
+- **Default**: 5
+
 ### chg_extrap
 
 - **Type**: String

python/pyabacus/src/ModuleNAO/CMakeLists.txt

@@ -12,7 +12,6 @@ list(APPEND _naos
     ${NAO_PATH}/two_center_bundle.cpp
     ${NAO_PATH}/two_center_integrator.cpp
     ${NAO_PATH}/two_center_table.cpp
-    ${NAO_PATH}/projgen.cpp
     # dependency
     ${ABACUS_SOURCE_DIR}/module_base/kernels/math_op.cpp
     # ${ABACUS_SOURCE_DIR}/module_psi/kernels/psi_memory_op.cpp

source/CMakeLists.txt

@@ -47,6 +47,7 @@ list(APPEND device_srcs
 
   module_hamilt_pw/hamilt_pwdft/kernels/force_op.cpp
   module_hamilt_pw/hamilt_pwdft/kernels/stress_op.cpp
+  module_hamilt_pw/hamilt_pwdft/kernels/onsite_op.cpp
   module_hamilt_pw/hamilt_pwdft/kernels/wf_op.cpp
   module_hamilt_pw/hamilt_pwdft/kernels/vnl_op.cpp
   module_base/kernels/math_op.cpp
@@ -60,6 +61,7 @@ if(USE_CUDA)
     module_hamilt_pw/hamilt_pwdft/kernels/cuda/ekinetic_op.cu
     module_hamilt_pw/hamilt_pwdft/kernels/cuda/meta_op.cu
     module_hamilt_pw/hamilt_stodft/kernels/cuda/hpsi_norm_op.cu
+    module_hamilt_pw/hamilt_pwdft/kernels/cuda/onsite_op.cu
     module_basis/module_pw/kernels/cuda/pw_op.cu
     module_hsolver/kernels/cuda/dngvd_op.cu
     module_hsolver/kernels/cuda/math_kernel_op.cu
@@ -83,6 +85,7 @@ if(USE_ROCM)
     module_hamilt_pw/hamilt_pwdft/kernels/rocm/veff_op.hip.cu
     module_hamilt_pw/hamilt_pwdft/kernels/rocm/ekinetic_op.hip.cu
     module_hamilt_pw/hamilt_pwdft/kernels/rocm/meta_op.hip.cu
+    module_hamilt_pw/hamilt_pwdft/kernels/rocm/onsite_op.hip.cu
     module_hamilt_pw/hamilt_stodft/kernels/rocm/hpsi_norm_op.hip.cu
     module_basis/module_pw/kernels/rocm/pw_op.hip.cu
     module_hsolver/kernels/rocm/dngvd_op.hip.cu

source/Makefile.Objects

@@ -635,10 +635,13 @@ OBJS_SRCPW=H_Ewald_pw.o\
     forces_nl.o\
     forces_cc.o\
     forces_scc.o\
+    forces_onsite.o\
+    onsite_proj_pw.o\
     fs_nonlocal_tools.o\
     fs_kin_tools.o\
     force_op.o\
     stress_op.o\
+    onsite_op.o\
     wf_op.o\
     vnl_op.o\
     global.o\
@@ -663,6 +666,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     stress_func_loc.o\
     stress_func_nl.o\
     stress_func_us.o\
+    stress_func_onsite.o\
     stress_pw.o\
     of_stress_pw.o\
     symmetry_rho.o\
@@ -673,7 +677,9 @@ OBJS_SRCPW=H_Ewald_pw.o\
     elecond.o\
     sto_tool.o\
     sto_elecond.o\
-    sto_dos.o
+    sto_dos.o\
+    onsite_projector.o\
+    onsite_proj_tools.o
 
 OBJS_VDW=vdw.o\
     vdwd2_parameters.o\
@@ -691,7 +697,8 @@ OBJS_DFTU=dftu.o\
       dftu_io.o\
       dftu_tools.o\
       dftu_occup.o\
-      dftu_hamilt.o
+      dftu_hamilt.o\
+      dftu_pw.o
 
 OBJS_DELTASPIN=basic_funcs.o\
       cal_mw_from_lambda.o\

source/module_base/CMakeLists.txt

@@ -54,6 +54,7 @@ add_library(
     spherical_bessel_transformer.cpp
     cubic_spline.cpp
     parallel_2d.cpp
+    projgen.cpp
     module_mixing/mixing_data.cpp
     module_mixing/mixing.cpp
     module_mixing/plain_mixing.cpp

source/module_basis/module_nao/CMakeLists.txt

@@ -14,7 +14,6 @@ if(ENABLE_LCAO)
     two_center_table.cpp
     two_center_integrator.cpp
     two_center_bundle.cpp
-    projgen.cpp
   )
 
   if(ENABLE_COVERAGE)

source/module_basis/module_nao/atomic_radials.cpp

@@ -7,7 +7,7 @@
 // FIXME: should update with pyabacus
 // #include "module_io/orb_io.h"
 
-#include "projgen.h"
+#include "module_base/projgen.h"
 
 #include <fstream>
 #include <iostream>

source/module_basis/module_nao/test/CMakeLists.txt

@@ -14,7 +14,6 @@ AddTest(
     ../atomic_radials.cpp
     ../radial_set.cpp
     ../numerical_radial.cpp
-    ../projgen.cpp
     ../../module_ao/ORB_atomic_lm.cpp
     ../../module_ao/ORB_atomic.cpp
     ../../../module_io/orb_io.cpp
@@ -84,7 +83,6 @@ AddTest(
     ../pswfc_radials.cpp
     ../radial_set.cpp
     ../numerical_radial.cpp
-    ../projgen.cpp
     ../sphbes_radials.cpp
     ../../module_ao/ORB_atomic_lm.cpp
     ../../module_ao/ORB_atomic.cpp
@@ -104,7 +102,6 @@ AddTest(
     ../pswfc_radials.cpp
     ../sphbes_radials.cpp
     ../radial_set.cpp
-    ../projgen.cpp
     ../numerical_radial.cpp
     ../two_center_bundle.cpp
     ../two_center_integrator.cpp
@@ -131,7 +128,6 @@ AddTest(
     ../real_gaunt_table.cpp
     ../radial_collection.cpp
     ../atomic_radials.cpp
-    ../projgen.cpp
     ../beta_radials.cpp
     ../hydrogen_radials.cpp
     ../pswfc_radials.cpp
@@ -158,7 +154,6 @@ AddTest(
     ../pswfc_radials.cpp
     ../sphbes_radials.cpp
     ../radial_set.cpp
-    ../projgen.cpp
     ../numerical_radial.cpp
     ../../../module_io/orb_io.cpp
   LIBS parameter ${math_libs} device base container orb 
@@ -179,7 +174,6 @@ AddTest(
     ../pswfc_radials.cpp
     ../sphbes_radials.cpp
     ../radial_set.cpp
-    ../projgen.cpp
     ../numerical_radial.cpp
     ../../../module_io/orb_io.cpp
   LIBS parameter ${math_libs} device base container orb 

source/module_basis/module_nao/test/projgen_test.cpp

@@ -1,4 +1,4 @@
-#include "module_basis/module_nao/projgen.h"
+#include "module_base/projgen.h"
 #include "gtest/gtest.h"
 
 #include "module_base/math_integral.h"

source/module_cell/read_atoms.cpp

@@ -455,7 +455,7 @@ bool UnitCell::read_atom_positions(std::ifstream &ifpos, std::ofstream &ofs_runn
             }
             else if(PARAM.inp.basis_type == "pw")
             {
-                if ((PARAM.inp.psi_initializer)&&(PARAM.inp.init_wfc.substr(0, 3) == "nao"))
+                if ((PARAM.inp.psi_initializer)&&(PARAM.inp.init_wfc.substr(0, 3) == "nao") || PARAM.inp.onsite_radius > 0.0)
                 {
                     std::string orbital_file = PARAM.inp.orbital_dir + orbital_fn[it];
                     this->read_orb_file(it, orbital_file, ofs_running, &(atoms[it]));

source/module_cell/unitcell.cpp

@@ -64,7 +64,6 @@ UnitCell::UnitCell() {
     atom_mass = nullptr;
     pseudo_fn = new std::string[1];
     pseudo_type = new std::string[1];
-    orbital_fn = new std::string[1];
 
     set_atom_flag = false;
 }
@@ -114,6 +113,15 @@ void UnitCell::bcast_unitcell() {
     Parallel_Common::bcast_int(lc[1]);
     Parallel_Common::bcast_int(lc[2]);
 
+    if(this->orbital_fn == nullptr)
+    {
+        this->orbital_fn = new std::string[ntype];
+    }
+    for (int i = 0; i < ntype; i++)
+    {
+        Parallel_Common::bcast_string(orbital_fn[i]);
+    }
+
     // distribute lattice vectors.
     Parallel_Common::bcast_double(a1.x);
     Parallel_Common::bcast_double(a1.y);

source/module_cell/unitcell.h

@@ -216,7 +216,7 @@ class UnitCell {
     std::string* pseudo_fn;
     std::string* pseudo_type; // pseudopotential types for each elements,
                               // sunliang added 2022-09-15.
-    std::string* orbital_fn;  // filenames of orbitals, liuyu add 2022-10-19
+    std::string* orbital_fn = nullptr;  // filenames of orbitals, liuyu add 2022-10-19
     std::string
         descriptor_file; // filenames of descriptor_file, liuyu add 2023-04-06
 

source/module_elecstate/elecstate.h

@@ -151,9 +151,7 @@ class ElecState
         return 0.0;
     }
 
-#ifdef __LCAO
     double get_dftu_energy();
-#endif
 
 #ifdef __DEEPKS
     double get_deepks_E_delta();

source/module_elecstate/elecstate_energy.cpp

@@ -287,7 +287,6 @@ void ElecState::cal_energies(const int type)
     }
 
     //! spin constrained energy
-#ifdef __LCAO
     if (PARAM.inp.sc_mag_switch)
     {
         this->f_en.escon = get_spin_constrain_energy();
@@ -298,7 +297,6 @@ void ElecState::cal_energies(const int type)
     {
         this->f_en.edftu = get_dftu_energy();
     }
-#endif
 
 #ifdef __DEEPKS
     // energy from deepks

source/module_elecstate/elecstate_energy_terms.cpp

@@ -34,12 +34,10 @@ double ElecState::get_solvent_model_Acav()
     return GlobalC::solvent_model.Acav;
 }
 
-#ifdef __LCAO
 double ElecState::get_dftu_energy()
 {
     return GlobalC::dftu.get_energy();
 }
-#endif
 
 #ifdef __DEEPKS
 double ElecState::get_deepks_E_delta()

source/module_esolver/esolver_ks_lcao.cpp

@@ -193,7 +193,8 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     // 7) initialize DFT+U
     if (PARAM.inp.dft_plus_u)
     {
-        GlobalC::dftu.init(ucell, &this->pv, this->kv.get_nks(), orb_);
+        auto* dftu = ModuleDFTU::DFTU::get_instance();
+        dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_);
     }
 
     // 8) initialize ppcell
@@ -1140,7 +1141,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
         //! Perform Mulliken charge analysis
         if (PARAM.inp.out_mul)
         {
-            ModuleIO::cal_mag(&(this->pv), this->p_hamilt, this->kv, this->pelec, ucell, istep, true);
+            ModuleIO::cal_mag(&(this->pv), this->p_hamilt, this->kv, this->pelec, this->two_center_bundle_, this->orb_, ucell, istep, true);
         }
     }