Merge branch 'develop' into feat/issue-7726-dftu-pw-nspin-fixes

Author: mohanchen

docs/advanced/input_files/input-main.md

@@ -653,7 +653,7 @@
 ### bndpar
 
 - **Type**: Integer
-- **Description**: Divide all processors into bndpar groups, and bands (only stochastic orbitals now) will be distributed among each group. It should be larger than 0.
+- **Description**: Divide all processors into bndpar groups for SDFT or the BPCG solver. bndpar must be positive, no greater than the number of MPI processes, and kpar * bndpar must divide the number of MPI processes exactly.
 - **Default**: 1
 
 ### latname
@@ -1265,7 +1265,7 @@
 - **Description**: The number of spin components of wave functions.
   - 1: Spin degeneracy
   - 2: Collinear spin polarized.
-  - 4: For the case of noncollinear polarized, nspin will be automatically set to 4 without being specified by the user.
+  - 4: Noncollinear or spin-orbit calculations. Set nspin to 4 explicitly when noncolin or lspinorb is enabled.
 - **Default**: 1
 
 ### smearing_method
@@ -1469,7 +1469,7 @@
 - **Type**: Boolean
 - **Description**: Whether to consider spin-orbit coupling (SOC) effect in the calculation.
   - True: Consider spin-orbit coupling effect. When enabled:
-  - nspin is automatically set to 4 (noncollinear spin representation)
+  - nspin must be explicitly set to 4 (noncollinear spin representation)
   - Symmetry is automatically disabled (SOC breaks inversion symmetry)
   - Requires full-relativistic pseudopotentials with has_so=true in the UPF header
   - False: Do not consider spin-orbit coupling effect.
@@ -1481,7 +1481,7 @@
 - **Type**: Boolean
 - **Description**: Whether to allow non-collinear magnetic moments, where magnetization can point in arbitrary directions (x, y, z components) rather than being constrained to the z-axis.
   - True: Allow non-collinear polarization. When enabled:
-  - nspin is automatically set to 4
+  - nspin must be explicitly set to 4
   - Wave function dimension is doubled (npol=2), and the number of occupied states is doubled
   - Charge density has 4 components (Pauli spin matrices)
   - Cannot be used with gamma_only=true
@@ -1535,8 +1535,8 @@
 - **Type**: Integer or string
 - **Availability**: *esolver_type = sdft*
 - **Description**: The number of stochastic orbitals
-  - > 0: Perform stochastic DFT. Increasing the number of bands improves accuracy and reduces stochastic errors; To perform mixed stochastic-deterministic DFT, you should set nbands, which represents the number of KS orbitals.
-  - 0: Perform Kohn-Sham DFT.
+  - 1-1000000: Perform stochastic DFT. Increasing the number of bands improves accuracy and reduces stochastic errors; To perform mixed stochastic-deterministic DFT, you should set nbands, which represents the number of KS orbitals.
+  - 0: Invalid. Use all for the complete-basis SDFT mode.
   - all: All complete basis sets are used to replace stochastic orbitals with the Chebyshev method (CT), resulting in the same results as KSDFT without stochastic errors.
 - **Default**: 256
 

docs/advanced/scf/spin.md

@@ -30,11 +30,11 @@ If **"nupdown"** is set to non-zero, number of spin-up and spin-down electrons w
 
 ## Noncollinear Spin Polarized Calculations
 The spin non-collinear polarization calculation corresponds to setting **"noncolin 1"**, in which case the coupling between spin up and spin down will be taken into account.
-In this case, nspin is automatically set to 4, which is usually not required to be specified manually.
+In this case, **"nspin 4"** must also be specified. ABACUS reports an input error instead of silently changing an incompatible or omitted nspin value.
 The weight of each band will not change, but the number of occupied states will be double.
 If the nbands parameter is set manually, it is generally set to twice what it would be when nspin<4.
 
-In general, non-collinear magnetic moment settings are often used in calculations considering [SOC effects](#soc-effects). When **"lspinorb 1"** in INPUT file, "nspin" is also automatically set to 4.
+In general, non-collinear magnetic moment settings are often used in calculations considering [SOC effects](#soc-effects). When **"lspinorb 1"** is set in INPUT, **"nspin 4"** is also required.
 
 Note: different settings for "noncolin" and "lspinorb" correspond to different calculations:
 
@@ -119,22 +119,22 @@ Example from a full-relativistic UPF file:
 - **PseudoDOJO**: Provides both scalar and full-relativistic versions
 - **ABACUS official**: [abacus.ustc.edu.cn](http://abacus.ustc.edu.cn/pseudo/list.htm)
 
-## Automatic Parameter Settings
+## Parameter Requirements and Automatic Settings
 
-When using SOC or non-collinear calculations, ABACUS automatically adjusts several parameters:
+When using SOC or non-collinear calculations, set the required spin representation explicitly. ABACUS still derives internal spin state and some related settings after validating the input:
 
 ### When `lspinorb=true`:
-1. **nspin**: Automatically set to 4 (noncollinear spin representation)
+1. **nspin**: Must be explicitly set to 4 (noncollinear spin representation)
 2. **Symmetry**: Automatically disabled (`symm_flag=-1`) because SOC breaks inversion symmetry
 3. **Magnetization**: NOT automatically set when `noncolin=0` (implies non-magnetic material with SOC)
 
 ### When `noncolin=true`:
-1. **nspin**: Automatically set to 4
+1. **nspin**: Must be explicitly set to 4
 2. **npol**: Set to 2 (wave function has two spinor components)
 3. **Magnetization**: Automatically set if user provides zero values (unless `lspinorb=1` and `noncolin=0`)
 
 ### Important Notes:
-- You do NOT need to manually set `nspin=4` when using `lspinorb=1` or `noncolin=1`
+- You must set `nspin=4` when using `lspinorb=1` or `noncolin=1`; missing or incompatible values are rejected during input validation
 - Symmetry operations are incompatible with SOC, so they are automatically turned off
 - For `lspinorb=1, noncolin=0`: This is a special case for non-magnetic materials with SOC, where magnetization is not initialized
 
@@ -172,7 +172,7 @@ basis_type          pw
 ecutwfc             50
 lspinorb            1        # Enable SOC
 noncolin            0        # No non-collinear magnetism
-# nspin will be automatically set to 4
+nspin               4        # Required spinor representation
 # symmetry will be automatically disabled
 ```
 
@@ -185,7 +185,7 @@ calculation         scf
 basis_type          lcao
 lspinorb            0        # No SOC
 noncolin            1        # Enable non-collinear magnetism
-# nspin will be automatically set to 4
+nspin               4        # Required spinor representation
 # Magnetization directions should be specified in STRU file
 ```
 
@@ -199,7 +199,7 @@ basis_type          pw
 ecutwfc             60
 lspinorb            1        # Enable SOC
 noncolin            1        # Enable non-collinear magnetism
-# nspin will be automatically set to 4
+nspin               4        # Required spinor representation
 # symmetry will be automatically disabled
 # Magnetization directions should be specified in STRU file
 ```
@@ -213,6 +213,7 @@ calculation         scf
 basis_type          pw
 ecutwfc             50
 lspinorb            1        # Enable SOC
+nspin               4        # Required spinor representation
 soc_lambda          0.5      # 50% SOC strength
 # Useful when full SOC overestimates or underestimates experimental results
 ```

docs/community/faq.md

@@ -50,7 +50,7 @@ To perform SOC calculations in ABACUS, follow these steps:
 2. **Use full-relativistic pseudopotentials**: SOC calculations require pseudopotentials with `has_so=true` in the UPF header
    - Download full-relativistic versions of SG15_ONCV pseudopotentials from [quantum-simulation.org](http://quantum-simulation.org/potentials/sg15_oncv/upf/)
    - Check the UPF file header for `relativistic="full"` and `has_so="T"`
-3. **Verify automatic settings**: When `lspinorb=1` is set, `nspin` is automatically set to 4 and symmetry is automatically disabled
+3. **Set the spin representation**: When `lspinorb=1` is set, explicitly set `nspin=4`; symmetry is automatically disabled
 
 **Basis set support**: Both `basis_type=pw` (plane wave) and `basis_type=lcao` (numerical atomic orbitals) support SOC calculations for both SCF and NSCF.
 

docs/parameters.yaml

@@ -108,7 +108,7 @@ parameters:
     category: System variables
     type: Integer
     description: |
-      Divide all processors into bndpar groups, and bands (only stochastic orbitals now) will be distributed among each group. It should be larger than 0.
+      Divide all processors into bndpar groups for SDFT or the BPCG solver. bndpar must be positive, no greater than the number of MPI processes, and kpar * bndpar must divide the number of MPI processes exactly.
     default_value: "1"
     unit: ""
     availability: ""
@@ -661,7 +661,7 @@ parameters:
       The number of spin components of wave functions.
       * 1: Spin degeneracy
       * 2: Collinear spin polarized.
-      * 4: For the case of noncollinear polarized, nspin will be automatically set to 4 without being specified by the user.
+      * 4: Noncollinear or spin-orbit calculations. Set nspin to 4 explicitly when noncolin or lspinorb is enabled.
     default_value: "1"
     unit: ""
     availability: ""
@@ -906,7 +906,7 @@ parameters:
     description: |
       Whether to consider spin-orbit coupling (SOC) effect in the calculation.
       * True: Consider spin-orbit coupling effect. When enabled:
-        * nspin is automatically set to 4 (noncollinear spin representation)
+        * nspin must be explicitly set to 4 (noncollinear spin representation)
         * Symmetry is automatically disabled (SOC breaks inversion symmetry)
         * Requires full-relativistic pseudopotentials with has_so=true in the UPF header
       * False: Do not consider spin-orbit coupling effect.
@@ -920,7 +920,7 @@ parameters:
     description: |
       Whether to allow non-collinear magnetic moments, where magnetization can point in arbitrary directions (x, y, z components) rather than being constrained to the z-axis.
       * True: Allow non-collinear polarization. When enabled:
-        * nspin is automatically set to 4
+        * nspin must be explicitly set to 4
         * Wave function dimension is doubled (npol=2), and the number of occupied states is doubled
         * Charge density has 4 components (Pauli spin matrices)
         * Cannot be used with gamma_only=true
@@ -2198,8 +2198,8 @@ parameters:
     type: Integer or string
     description: |
       The number of stochastic orbitals
-      * > 0: Perform stochastic DFT. Increasing the number of bands improves accuracy and reduces stochastic errors; To perform mixed stochastic-deterministic DFT, you should set nbands, which represents the number of KS orbitals.
-      * 0: Perform Kohn-Sham DFT.
+      * 1-1000000: Perform stochastic DFT. Increasing the number of bands improves accuracy and reduces stochastic errors; To perform mixed stochastic-deterministic DFT, you should set nbands, which represents the number of KS orbitals.
+      * 0: Invalid. Use all for the complete-basis SDFT mode.
       * all: All complete basis sets are used to replace stochastic orbitals with the Chebyshev method (CT), resulting in the same results as KSDFT without stochastic errors.
     default_value: "256"
     unit: ""

examples/21_deepks/03_lcao_CsPbI3/INPUT

@@ -29,6 +29,7 @@ deepks_model         model.ptg
 
 #Parameters (7.SOC)
 lspinorb             1
+nspin                4
 
 
 

source/Makefile.Objects

@@ -343,18 +343,28 @@ OBJS_HAMILT_OF=kedf_tf.o\
 OBJS_HAMILT_LCAO=hamilt_lcao.o\
     operator_lcao.o\
     ekinetic.o\
+    ekinetic_fs.o\
+    ekinetic_dh.o\
     nonlocal.o\
+    nonlocal_dh.o\
+    nonlocal_fs.o\
     overlap.o\
+    overlap_fs.o\
     td_ekinetic_lcao.o\
     td_nonlocal_lcao.o\
     td_pot_hybrid.o\
+    td_pot_hybrid_fs.o\
     veff_lcao.o\
+    veff_dh.o\
     meta_lcao.o\
     op_dftu_lcao.o\
     deepks_lcao.o\
     op_exx_lcao.o\
     dspin_lcao.o\
+    dspin_fs.o\
     dftu_lcao.o\
+    dftu_fs.o\
+    operator_fs_utils.o\
 
 OBJS_HCONTAINER=base_matrix.o\
     atom_pair.o\
@@ -721,7 +731,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     setup_pwrho.o\
     setup_pwwfc.o\
     update_cell_pw.o\
-    dftu_pw.o\
+    setup_dftu_pw.o\
     deltaspin_pw.o\
     forces.o\
     forces_us.o\
@@ -793,7 +803,7 @@ OBJS_DFTU=dftu.o\
       dftu_tools.o\
       dftu_occup.o\
       dftu_hamilt.o\
-      dftu_pw.o
+      setup_dftu_pw.o
 
 OBJS_DELTASPIN=basic_funcs.o\
       cal_mw_from_lambda.o\

source/source_esolver/esolver_double_xc.cpp

@@ -1,7 +1,6 @@
 #include "esolver_double_xc.h"
 
 #include "source_hamilt/module_ewald/H_Ewald_pw.h"
-#include "source_hamilt/module_vdw/vdw.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 #ifdef __MLALGO
 #include "source_lcao/module_deepks/LCAO_deepks.h"
@@ -121,13 +120,9 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);
 
     //----------------------------------------------------------
-    //! calculate D2 or D3 vdW
+    //! Reuse the vdW correction prepared by ESolver_FP::before_scf.
     //----------------------------------------------------------
-    auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
-    if (vdw_solver != nullptr)
-    {
-        this->pelec_base->f_en.evdw = vdw_solver->get_energy();
-    }
+    this->pelec_base->f_en.evdw = this->pelec->f_en.evdw;
 
     //----------------------------------------------------------
     //! calculate ewald energy
@@ -398,6 +393,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&
     this->deepks.dpks_out_type = "base"; // for deepks method
 
     fsl.getForceStress(ucell,
+                       this->get_vdw_result(),
                        PARAM.inp.cal_force,
                        PARAM.inp.cal_stress,
                        PARAM.inp.test_force,

source/source_esolver/esolver_fp.cpp

@@ -190,11 +190,18 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
                                     GlobalV::ofs_running, GlobalV::ofs_warning);
     }
 
-    //! calculate D2 or D3 vdW
+    //! Evaluate the vdW correction once for this ionic configuration.
+    this->vdw_result_.reset();
     auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running));
     if (vdw_solver != nullptr)
     {
-        this->pelec->f_en.evdw = vdw_solver->get_energy();
+        const vdw::VdwRequest request(PARAM.inp.cal_force, PARAM.inp.cal_stress);
+        this->vdw_result_.reset(new vdw::VdwResult(vdw_solver->evaluate(request)));
+        this->pelec->f_en.evdw = this->vdw_result_->energy;
+    }
+    else
+    {
+        this->pelec->f_en.evdw = 0.0;
     }
 
     //! calculate ewald energy

source/source_esolver/esolver_fp.h

@@ -12,6 +12,12 @@
 #include "source_pw/module_pwdft/vl_pw.h"             // local pseudopotential
 
 #include <fstream>
+#include <memory>
+
+namespace vdw
+{
+struct VdwResult;
+}
 
 //! The First-Principles (FP) Energy Solver Class
 /**
@@ -42,6 +48,11 @@ class ESolver_FP : public ESolver
 
     virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver);
 
+    const vdw::VdwResult* get_vdw_result() const { return this->vdw_result_.get(); }
+
+    //! vdW correction evaluated once for the current ionic configuration.
+    std::unique_ptr<const vdw::VdwResult> vdw_result_;
+
     //! These pointers will be deleted in the free_pointers() function every ion step.
     elecstate::ElecState* pelec = nullptr; ///< Electronic states
 

source/source_esolver/esolver_ks_lcao.cpp

@@ -5,7 +5,7 @@
 #include "source_lcao/module_deltaspin/spin_constrain.h"
 #include "source_lcao/module_deltaspin/deltaspin_lcao.h"
 #include "source_lcao/dftu_lcao.h"
-#include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp
+#include "source_hamilt/hs_matrix_k.h"
 #include "source_estate/module_charge/symmetry_rho.h"
 #include "source_lcao/LCAO_domain.h" // need DeePKS_init
 #include "source_lcao/FORCE_STRESS.h"
@@ -248,7 +248,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(BaseCell& basecell, ModuleBase::matrix&
 
     deepks.dpks_out_type = "tot";  // for deepks method
 
-    fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress, 
+    fsl.getForceStress(ucell, this->get_vdw_result(), PARAM.inp.cal_force, PARAM.inp.cal_stress,
                        PARAM.inp.test_force, PARAM.inp.test_stress,
                        this->gd, this->pv, this->pelec, this->dmat, this->psi,
                        two_center_bundle_, orb_, force, this->scs,