diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp
index f9f493eaf1d..199be857f82 100644
--- a/source/source_lcao/FORCE_STRESS.cpp
+++ b/source/source_lcao/FORCE_STRESS.cpp
@@ -276,17 +276,41 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
         const int nks = (PARAM.inp.nspin == 1 || PARAM.inp.nspin == 2) ? 1 : kv.get_nks();
         if (PARAM.globalv.gamma_only_local)
         {
-            DeePKS_domain::cal_f_delta<double>(deepks.ld.dm_r, ucell, orb, gd,
-                                               *flk.ParaV, nks, deepks.ld.deepks_param,
-                                               kv.kvec_d, deepks.ld.phialpha, deepks.ld.gedm,
-                                               fvnl_dalpha, isstress, svnl_dalpha);
+            DeePKS_domain::cal_f_delta<double>(
+                ucell,
+                orb,
+                gd,
+                *flk.ParaV,
+                nks,
+                deepks.ld.deepks_param,
+                kv.kvec_d,
+                deepks.ld.phialpha,
+                fvnl_dalpha,
+                isstress,
+                svnl_dalpha,
+                deepks.ld.dm_r,
+                deepks.ld.gedm,
+                (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv) ? deepks.ld.dm_r_mag : nullptr,
+                (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv) ? deepks.ld.gedm_mag : nullptr);
         }
         else
         {
-            DeePKS_domain::cal_f_delta<std::complex<double>>(deepks.ld.dm_r, ucell, orb, gd,
-                                                              *flk.ParaV, nks, deepks.ld.deepks_param,
-                                                              kv.kvec_d, deepks.ld.phialpha, deepks.ld.gedm,
-                                                              fvnl_dalpha, isstress, svnl_dalpha);
+            DeePKS_domain::cal_f_delta<std::complex<double>>(
+                ucell,
+                orb,
+                gd,
+                *flk.ParaV,
+                nks,
+                deepks.ld.deepks_param,
+                kv.kvec_d,
+                deepks.ld.phialpha,
+                fvnl_dalpha,
+                isstress,
+                svnl_dalpha,
+                deepks.ld.dm_r,
+                deepks.ld.gedm,
+                (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv) ? deepks.ld.dm_r_mag : nullptr,
+                (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv) ? deepks.ld.gedm_mag : nullptr);
         }
 
         if (isforce)
diff --git a/source/source_lcao/FORCE_gamma.cpp b/source/source_lcao/FORCE_gamma.cpp
index 85059648986..379ceb29459 100644
--- a/source/source_lcao/FORCE_gamma.cpp
+++ b/source/source_lcao/FORCE_gamma.cpp
@@ -217,8 +217,7 @@ void Force_LCAO<double>::ftable(const bool isforce,
     {
         // No need to update E_delta here since it have been done in LCAO_Deepks_Interface in after_scf
         const int nks = 1;
-        DeePKS_domain::cal_f_delta<double>(deepks.ld.dm_r,
-                                           ucell,
+        DeePKS_domain::cal_f_delta<double>(ucell,
                                            orb,
                                            gd,
                                            *this->ParaV,
@@ -226,10 +225,11 @@ void Force_LCAO<double>::ftable(const bool isforce,
                                            deepks.ld.deepks_param,
                                            kv->kvec_d,
                                            deepks.ld.phialpha,
-                                           deepks.ld.gedm,
                                            fvnl_dalpha,
                                            isstress,
-                                           svnl_dalpha);
+                                           svnl_dalpha,
+                                           deepks.ld.dm_r,
+                                           deepks.ld.gedm);
     }
 #endif
 
diff --git a/source/source_lcao/FORCE_k.cpp b/source/source_lcao/FORCE_k.cpp
index 8d5a090a8dd..475d8368884 100644
--- a/source/source_lcao/FORCE_k.cpp
+++ b/source/source_lcao/FORCE_k.cpp
@@ -224,8 +224,7 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
     if (PARAM.inp.deepks_scf)
     {
         // No need to update E_delta since it have been done in LCAO_Deepks_Interface in after_scf
-        DeePKS_domain::cal_f_delta<std::complex<double>>(deepks.ld.dm_r,
-                                                         ucell,
+        DeePKS_domain::cal_f_delta<std::complex<double>>(ucell,
                                                          orb,
                                                          gd,
                                                          pv,
@@ -233,10 +232,11 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
                                                          deepks.ld.deepks_param,
                                                          kv->kvec_d,
                                                          deepks.ld.phialpha,
-                                                         deepks.ld.gedm,
                                                          fvnl_dalpha,
                                                          isstress,
-                                                         svnl_dalpha);
+                                                         svnl_dalpha,
+                                                         deepks.ld.dm_r,
+                                                         deepks.ld.gedm);
     }
 #endif
 
diff --git a/source/source_lcao/module_deepks/LCAO_deepks.cpp b/source/source_lcao/module_deepks/LCAO_deepks.cpp
index 41c7e13fbea..6ae7eae91be 100644
--- a/source/source_lcao/module_deepks/LCAO_deepks.cpp
+++ b/source/source_lcao/module_deepks/LCAO_deepks.cpp
@@ -11,6 +11,7 @@ LCAO_Deepks<T>::LCAO_Deepks()
 {
     deepks_param.inl_index = new ModuleBase::IntArray[1];
     gedm = nullptr;
+    gedm_mag = nullptr;
     this->phialpha.resize(1);
 }
 
@@ -33,6 +34,16 @@ LCAO_Deepks<T>::~LCAO_Deepks()
         }
         delete[] gedm;
     }
+    if (gedm_mag)
+    {
+        for (int inl = 0; inl < this->deepks_param.inlmax; inl++)
+        {
+            delete[] gedm_mag[inl];
+        }
+        delete[] gedm_mag;
+    }
+    pdm_mag.clear();
+    delete dm_r_mag;
 }
 
 template <typename T>
@@ -124,6 +135,15 @@ void LCAO_Deepks<T>::init(const LCAO_Orbitals& orb,
         }
     }
 
+    if (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv) // magnetization-channel PDM
+    {
+        this->pdm_mag.resize(this->deepks_param.inlmax);
+        for (int inl = 0; inl < this->deepks_param.inlmax; ++inl)
+        {
+            this->pdm_mag[inl] = torch::zeros_like(this->pdm[inl]);
+        }
+    }
+
     this->pv = &pv_in;
 
     ModuleBase::timer::end("LCAO_Deepks", "init");
@@ -202,6 +222,16 @@ void LCAO_Deepks<T>::allocate_V_delta(const int nat, const int nks)
         ModuleBase::GlobalFunc::ZEROS(this->gedm[inl], pdm_size);
     }
 
+    if (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv) // magnetization-channel gedm
+    {
+        this->gedm_mag = new double*[this->deepks_param.inlmax];
+        for (int inl = 0; inl < this->deepks_param.inlmax; inl++)
+        {
+            this->gedm_mag[inl] = new double[pdm_size];
+            ModuleBase::GlobalFunc::ZEROS(this->gedm_mag[inl], pdm_size);
+        }
+    }
+
     ModuleBase::timer::end("LCAO_Deepks", "allocate_V_delta");
     return;
 }
@@ -249,6 +279,11 @@ void LCAO_Deepks<T>::init_DMR(const UnitCell& ucell,
                                    this->dm_r->insert_pair(dm_pair);
                                });
     this->dm_r->allocate(nullptr, true);
+
+    if (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv) // magnetization-channel real-space DM
+    {
+        this->dm_r_mag = new hamilt::HContainer<double>(*this->dm_r);
+    }
 }
 
 template <typename T>
diff --git a/source/source_lcao/module_deepks/LCAO_deepks.h b/source/source_lcao/module_deepks/LCAO_deepks.h
index 75daa9ffad1..07751e7e021 100644
--- a/source/source_lcao/module_deepks/LCAO_deepks.h
+++ b/source/source_lcao/module_deepks/LCAO_deepks.h
@@ -79,6 +79,11 @@ class LCAO_Deepks
     /// dE/dD, autograd from loaded model(E: Ry)
     double** gedm = nullptr; //[tot_Inl][(2l+1)*(2l+1)]
 
+    // magnetization-channel (rho_up - rho_dn) counterparts, used only for nspin=2
+    hamilt::HContainer<double>* dm_r_mag = nullptr;
+    std::vector<torch::Tensor> pdm_mag;
+    double** gedm_mag = nullptr;
+
     // functions for hr status: 1. get value; 2. set value;
     int get_hr_cal()
     {
diff --git a/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp b/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp
index 928617f4cdb..21cc3b4e927 100644
--- a/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp
+++ b/source/source_lcao/module_deepks/LCAO_deepks_interface.cpp
@@ -100,6 +100,7 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
 
     const int nspin = PARAM.inp.nspin;
     const int nk = nks / nspin;
+    const bool deepks_spin2 = (nspin == 2 && !PARAM.inp.deepks_equiv);
 
     const bool is_after_scf = (iter == -1); // called in after_scf, not in electronic steps
     const bool output_base
@@ -130,13 +131,26 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
         DeePKS_domain::cal_descriptor(nat, deepks_param, pdm, descriptor); // final descriptor
         DeePKS_domain::check_descriptor(deepks_param, ucell, PARAM.globalv.global_out_dir, descriptor, rank);
 
+        // nspin=2 (traditional): also build the magnetization-channel descriptor,
+        // reused for the dm_eig label and the 2-channel model below.
+        std::vector<torch::Tensor> descriptor_mag;
+        if (deepks_spin2)
+        {
+            DeePKS_domain::update_dmr(kvec_d, dm->get_DMK_vector(), ucell, orb, *ParaV, GridD, ld->dm_r_mag, 2, true);
+            bool init_pdm_mag = false;
+            DeePKS_domain::cal_pdm<
+                TK>(init_pdm_mag, deepks_param, kvec_d, ld->dm_r_mag, phialpha, ucell, orb, GridD, *ParaV, ld->pdm_mag);
+            DeePKS_domain::cal_descriptor(nat, deepks_param, ld->pdm_mag, descriptor_mag);
+        }
+
         const std::string file_d = get_filename("dm_eig", PARAM.inp.deepks_out_labels, iter);
         LCAO_deepks_io::save_npy_d(nat,
                                    PARAM.inp.deepks_equiv,
                                    deepks_param,
-                                   descriptor,
                                    file_d,
-                                   rank); // libnpy needed
+                                   rank,
+                                   descriptor,
+                                   descriptor_mag); // libnpy needed
 
         if (PARAM.inp.deepks_scf)
         {
@@ -146,9 +160,31 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
             {
                 DeePKS_domain::cal_edelta_gedm_equiv(nat, deepks_param, descriptor, ld->model_deepks, ld->gedm, E_delta, rank);
             }
-            else
+            else // traditional version
             {
-                DeePKS_domain::cal_edelta_gedm(nat, deepks_param, descriptor, pdm, ld->model_deepks, ld->gedm, E_delta);
+                if (deepks_spin2)
+                {
+                    DeePKS_domain::cal_edelta_gedm(nat,
+                                                   deepks_param,
+                                                   ld->model_deepks,
+                                                   E_delta,
+                                                   descriptor,
+                                                   pdm,
+                                                   ld->gedm,
+                                                   descriptor_mag,
+                                                   ld->pdm_mag,
+                                                   ld->gedm_mag);
+                }
+                else
+                {
+                    DeePKS_domain::cal_edelta_gedm(nat,
+                                                   deepks_param,
+                                                   ld->model_deepks,
+                                                   E_delta,
+                                                   descriptor,
+                                                   pdm,
+                                                   ld->gedm);
+                }
             }
         }
     }
@@ -158,9 +194,14 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
     {
         // Used for deepks_scf == 1 or deepks_out_freq_elec!=0, for *precalc items, not for deepks_out_labels=2
         std::vector<torch::Tensor> gevdm;
+        std::vector<torch::Tensor> gevdm_mag;
         if (output_precalc)
         {
             DeePKS_domain::cal_gevdm(nat, deepks_param, pdm, gevdm);
+            if (deepks_spin2)
+            {
+                DeePKS_domain::cal_gevdm(nat, deepks_param, ld->pdm_mag, gevdm_mag);
+            }
         }
 
         //================================================================================
@@ -206,6 +247,18 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
 
                     torch::Tensor gvx;
                     DeePKS_domain::cal_gvx(ucell.nat, deepks_param, gevdm, gdmx, gvx, rank);
+                    if (deepks_spin2)
+                    {
+                        torch::Tensor gdmx_mag;
+                        DeePKS_domain::cal_gdmx<
+                            TK>(nks, deepks_param, kvec_d, phialpha, ld->dm_r_mag, ucell, orb, *ParaV, GridD, gdmx_mag);
+                        torch::Tensor gvx_mag;
+                        DeePKS_domain::cal_gvx(ucell.nat, deepks_param, gevdm_mag, gdmx_mag, gvx_mag, rank);
+                        if (rank == 0)
+                        {
+                            gvx = torch::stack({gvx, gvx_mag}, gvx.dim() - 1);
+                        }
+                    }
                     const std::string file_gradvx = get_filename("gradvx", PARAM.inp.deepks_out_labels, iter);
                     LCAO_deepks_io::save_tensor2npy<double>(file_gradvx, gvx, rank);
 
@@ -230,6 +283,26 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
 
                     torch::Tensor gvepsl;
                     DeePKS_domain::cal_gvepsl(ucell.nat, deepks_param, gevdm, gdmepsl, gvepsl, rank);
+                    if (deepks_spin2)
+                    {
+                        torch::Tensor gdmepsl_mag;
+                        DeePKS_domain::cal_gdmepsl<TK>(nks,
+                                                       deepks_param,
+                                                       kvec_d,
+                                                       phialpha,
+                                                       ld->dm_r_mag,
+                                                       ucell,
+                                                       orb,
+                                                       *ParaV,
+                                                       GridD,
+                                                       gdmepsl_mag);
+                        torch::Tensor gvepsl_mag;
+                        DeePKS_domain::cal_gvepsl(ucell.nat, deepks_param, gevdm_mag, gdmepsl_mag, gvepsl_mag, rank);
+                        if (rank == 0)
+                        {
+                            gvepsl = torch::stack({gvepsl, gvepsl_mag}, gvepsl.dim() - 1);
+                        }
+                    }
                     const std::string file_gvepsl = get_filename("gvepsl", PARAM.inp.deepks_out_labels, iter);
                     LCAO_deepks_io::save_tensor2npy<double>(file_gvepsl, gvepsl, rank);
 
@@ -346,6 +419,7 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
 
                 ModuleBase::matrix o_delta(nks, range);
                 torch::Tensor orbital_precalc;
+                torch::Tensor orbital_precalc_mag;
                 for (int ir = 0; ir < range; ++ir)
                 {
                     std::vector<TH> dm_bandgap(nks);
@@ -374,6 +448,31 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
                         orbital_precalc = torch::cat({orbital_precalc, orbital_precalc_temp}, 0);
                     }
 
+                    if (deepks_spin2)
+                    {
+                        torch::Tensor orbital_precalc_mag_temp;
+                        DeePKS_domain::cal_orbital_precalc<TK, TH>(dm_bandgap,
+                                                                   nat,
+                                                                   nks,
+                                                                   deepks_param,
+                                                                   kvec_d,
+                                                                   phialpha,
+                                                                   gevdm_mag,
+                                                                   ucell,
+                                                                   orb,
+                                                                   *ParaV,
+                                                                   GridD,
+                                                                   orbital_precalc_mag_temp);
+                        if (ir == 0)
+                        {
+                            orbital_precalc_mag = orbital_precalc_mag_temp;
+                        }
+                        else
+                        {
+                            orbital_precalc_mag = torch::cat({orbital_precalc_mag, orbital_precalc_mag_temp}, 0);
+                        }
+                    }
+
                     if (PARAM.inp.deepks_scf)
                     {
                         DeePKS_domain::cal_o_delta<TK, TH>(dm_bandgap, *h_delta, o_delta_temp, *ParaV, nks, nspin);
@@ -384,6 +483,14 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
                     }
                 }
                 // save obase and orbital_precalc
+                if (deepks_spin2)
+                {
+                    if (rank == 0)
+                    {
+                        orbital_precalc
+                            = torch::stack({orbital_precalc, orbital_precalc_mag}, orbital_precalc.dim() - 1);
+                    }
+                }
                 const std::string file_orbpre = get_filename("orbpre", PARAM.inp.deepks_out_labels, iter);
                 LCAO_deepks_io::save_tensor2npy<double>(file_orbpre, orbital_precalc, rank);
 
@@ -478,7 +585,27 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
                                                            *ParaV,
                                                            GridD,
                                                            vdr_precalc);
-
+                            if (deepks_spin2)
+                            {
+                                torch::Tensor vdr_precalc_mag;
+                                DeePKS_domain::cal_vdr_precalc(nlocal,
+                                                               nat,
+                                                               nks,
+                                                               R_size,
+                                                               deepks_param,
+                                                               kvec_d,
+                                                               phialpha,
+                                                               gevdm_mag,
+                                                               ucell,
+                                                               orb,
+                                                               *ParaV,
+                                                               GridD,
+                                                               vdr_precalc_mag);
+                                if (rank == 0)
+                                {
+                                    vdr_precalc = torch::stack({vdr_precalc, vdr_precalc_mag}, vdr_precalc.dim() - 1);
+                                }
+                            }
                             const std::string file_vdrpre = PARAM.globalv.global_out_dir + "deepks_vdrpre.npy";
                             LCAO_deepks_io::save_tensor2npy<double>(file_vdrpre, vdr_precalc, rank);
                         }
@@ -486,6 +613,15 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
                         {
                             torch::Tensor gevdm_out;
                             DeePKS_domain::prepare_gevdm(nat, deepks_param, orb, gevdm, gevdm_out);
+                            if (deepks_spin2)
+                            {
+                                torch::Tensor gevdm_out_mag;
+                                DeePKS_domain::prepare_gevdm(nat, deepks_param, orb, gevdm_mag, gevdm_out_mag);
+                                if (rank == 0)
+                                {
+                                    gevdm_out = torch::stack({gevdm_out, gevdm_out_mag}, gevdm_out.dim() - 1);
+                                }
+                            }
                             const std::string file_gevdm = PARAM.globalv.global_out_dir + "deepks_gevdm.npy";
                             LCAO_deepks_io::save_tensor2npy<double>(file_gevdm, gevdm_out, rank);
 
@@ -563,7 +699,27 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
                                                            *ParaV,
                                                            GridD,
                                                            v_delta_precalc);
-
+                    if (deepks_spin2)
+                    {
+                        torch::Tensor v_delta_precalc_mag;
+                        DeePKS_domain::cal_v_delta_precalc<TK>(nlocal,
+                                                               nat,
+                                                               nks,
+                                                               deepks_param,
+                                                               kvec_d,
+                                                               phialpha,
+                                                               gevdm_mag,
+                                                               ucell,
+                                                               orb,
+                                                               *ParaV,
+                                                               GridD,
+                                                               v_delta_precalc_mag);
+                        if (rank == 0)
+                        {
+                            v_delta_precalc
+                                = torch::stack({v_delta_precalc, v_delta_precalc_mag}, v_delta_precalc.dim() - 1);
+                        }
+                    }
                     const std::string file_vdpre = get_filename("vdpre", PARAM.inp.deepks_out_labels, iter);
                     LCAO_deepks_io::save_tensor2npy<TK>(file_vdpre, v_delta_precalc, rank);
                 }
@@ -577,6 +733,15 @@ void LCAO_Deepks_Interface<TK, TR>::out_deepks_labels(const double& etot,
 
                     torch::Tensor gevdm_out;
                     DeePKS_domain::prepare_gevdm(nat, deepks_param, orb, gevdm, gevdm_out);
+                    if (deepks_spin2)
+                    {
+                        torch::Tensor gevdm_out_mag;
+                        DeePKS_domain::prepare_gevdm(nat, deepks_param, orb, gevdm_mag, gevdm_out_mag);
+                        if (rank == 0)
+                        {
+                            gevdm_out = torch::stack({gevdm_out, gevdm_out_mag}, gevdm_out.dim() - 1);
+                        }
+                    }
                     const std::string file_gevdm = get_filename("gevdm", PARAM.inp.deepks_out_labels, iter);
                     LCAO_deepks_io::save_tensor2npy<double>(file_gevdm, gevdm_out, rank);
                 }
diff --git a/source/source_lcao/module_deepks/LCAO_deepks_io.cpp b/source/source_lcao/module_deepks/LCAO_deepks_io.cpp
index 12c38dfcf4f..d13cd2c43f9 100644
--- a/source/source_lcao/module_deepks/LCAO_deepks_io.cpp
+++ b/source/source_lcao/module_deepks/LCAO_deepks_io.cpp
@@ -77,9 +77,10 @@ void LCAO_deepks_io::load_npy_gedm(const int nat,
 void LCAO_deepks_io::save_npy_d(const int nat,
                                 const bool deepks_equiv,
                                 const DeePKS_Param& deepks_param,
-                                const std::vector<torch::Tensor>& descriptor,
                                 const std::string& dm_eig_file,
-                                const int rank)
+                                const int rank,
+                                const std::vector<torch::Tensor>& descriptor,
+                                const std::vector<torch::Tensor>& descriptor_mag)
 {
     ModuleBase::TITLE("LCAO_deepks_io", "save_npy_d");
 
@@ -102,12 +103,42 @@ void LCAO_deepks_io::save_npy_d(const int nat,
                 npy_des.push_back(accessor[im]);
             }
         }
-        const long unsigned dshape[]
-            = {static_cast<unsigned long>(nat), static_cast<unsigned long>(deepks_param.des_per_atom)};
-        if (rank == 0)
+        if (descriptor_mag.empty())
         {
+            const long unsigned dshape[]
+                = {static_cast<unsigned long>(nat), static_cast<unsigned long>(deepks_param.des_per_atom)};
             npy::SaveArrayAsNumpy(dm_eig_file, false, 2, dshape, npy_des);
         }
+        else
+        {
+            // nspin=2: store charge and magnetization channels -> (nat, 2, des)
+            std::vector<double> npy_mag;
+            for (int inl = 0; inl < deepks_param.inlmax; ++inl)
+            {
+                auto accessor = descriptor_mag[inl].accessor<double, 1>();
+                int nm = 2 * deepks_param.inl2l[inl] + 1;
+                for (int im = 0; im < nm; im++)
+                {
+                    npy_mag.push_back(accessor[im]);
+                }
+            }
+            const int des = deepks_param.des_per_atom;
+            std::vector<double> npy_both;
+            npy_both.reserve(2 * nat * des);
+            for (int iat = 0; iat < nat; ++iat)
+            {
+                for (int i = 0; i < des; ++i)
+                {
+                    npy_both.push_back(npy_des[iat * des + i]);
+                }
+                for (int i = 0; i < des; ++i)
+                {
+                    npy_both.push_back(npy_mag[iat * des + i]);
+                }
+            }
+            const long unsigned dshape[] = {static_cast<unsigned long>(nat), 2ul, static_cast<unsigned long>(des)};
+            npy::SaveArrayAsNumpy(dm_eig_file, false, 3, dshape, npy_both);
+        }
     }
     else
     {
diff --git a/source/source_lcao/module_deepks/LCAO_deepks_io.h b/source/source_lcao/module_deepks/LCAO_deepks_io.h
index 46e9272508f..b7b5987c3fc 100644
--- a/source/source_lcao/module_deepks/LCAO_deepks_io.h
+++ b/source/source_lcao/module_deepks/LCAO_deepks_io.h
@@ -47,9 +47,10 @@ void load_npy_gedm(const int nat, const int des_per_atom, double** gedm, double&
 void save_npy_d(const int nat,
                 const bool deepks_equiv,
                 const DeePKS_Param& deepks_param,
-                const std::vector<torch::Tensor>& descriptor,
                 const std::string& dm_eig_file,
-                const int rank);
+                const int rank,
+                const std::vector<torch::Tensor>& descriptor,
+                const std::vector<torch::Tensor>& descriptor_mag = {});
 
 // save energy
 void save_npy_e(const double& e, /**<[in] \f$E_{base}\f$ or \f$E_{tot}\f$, in Ry*/
diff --git a/source/source_lcao/module_deepks/deepks_basic.cpp b/source/source_lcao/module_deepks/deepks_basic.cpp
index f15996be682..b9e6d63f4a0 100644
--- a/source/source_lcao/module_deepks/deepks_basic.cpp
+++ b/source/source_lcao/module_deepks/deepks_basic.cpp
@@ -217,20 +217,35 @@ void DeePKS_domain::cal_edelta_gedm_equiv(const int nat,
 // E_delta is also calculated here
 void DeePKS_domain::cal_edelta_gedm(const int nat,
                                     const DeePKS_Param& deepks_param,
+                                    torch::jit::script::Module& model_deepks,
+                                    double& E_delta,
                                     const std::vector<torch::Tensor>& descriptor,
                                     const std::vector<torch::Tensor>& pdm,
-                                    torch::jit::script::Module& model_deepks,
                                     double** gedm,
-                                    double& E_delta)
+                                    const std::vector<torch::Tensor>& descriptor_mag,
+                                    const std::vector<torch::Tensor>& pdm_mag,
+                                    double** gedm_mag)
 {
     ModuleBase::TITLE("DeePKS_domain", "cal_edelta_gedm");
     ModuleBase::timer::start("DeePKS_domain", "cal_edelta_gedm");
 
+    const bool two_channel = !descriptor_mag.empty();
+
     // forward
     std::vector<torch::jit::IValue> inputs;
-
     // input_dim:(natom, des_per_atom)
-    inputs.push_back(torch::cat(descriptor, 0).reshape({1, nat, deepks_param.des_per_atom}));
+    if (!two_channel)
+    {
+        // nspin=1: only charge channel -> (1, natom, des_per_atom)
+        inputs.push_back(torch::cat(descriptor, 0).reshape({1, nat, deepks_param.des_per_atom}));
+    }
+    else
+    {
+        // nspin=2: charge and magnetization channels -> (1, natom, 2, des_per_atom)
+        torch::Tensor d_charge = torch::cat(descriptor, 0).reshape({1, nat, deepks_param.des_per_atom});
+        torch::Tensor d_mag = torch::cat(descriptor_mag, 0).reshape({1, nat, deepks_param.des_per_atom});
+        inputs.push_back(torch::stack({d_charge, d_mag}, 2));
+    }
     std::vector<torch::Tensor> ec;
     try
     {
@@ -250,8 +265,13 @@ void DeePKS_domain::cal_edelta_gedm(const int nat,
     // cal gedm
     std::vector<torch::Tensor> gedm_shell;
     gedm_shell.push_back(torch::ones_like(ec[0]));
+    std::vector<torch::Tensor> grad_inputs(pdm.begin(), pdm.end());
+    if (two_channel)
+    {
+        grad_inputs.insert(grad_inputs.end(), pdm_mag.begin(), pdm_mag.end());
+    }
     std::vector<torch::Tensor> gedm_tensor = torch::autograd::grad(ec,
-                                                                   pdm,
+                                                                   grad_inputs,
                                                                    gedm_shell,
                                                                    /*retain_grad=*/true,
                                                                    /*create_graph=*/false,
@@ -261,6 +281,16 @@ void DeePKS_domain::cal_edelta_gedm(const int nat,
     for (int inl = 0; inl < deepks_param.inlmax; ++inl)
     {
         int nm = 2 * deepks_param.inl2l[inl] + 1;
+        // allow_unused=true may return an undefined gradient when the model
+        // output does not depend on this PDM block; treat it as zero
+        if (!gedm_tensor[inl].defined())
+        {
+            for (int index = 0; index < nm * nm; ++index)
+            {
+                gedm[inl][index] = 0.0;
+            }
+            continue;
+        }
         auto accessor = gedm_tensor[inl].accessor<double, 2>();
         for (int m1 = 0; m1 < nm; ++m1)
         {
@@ -271,6 +301,30 @@ void DeePKS_domain::cal_edelta_gedm(const int nat,
             }
         }
     }
+    if (two_channel)
+    {
+        for (int inl = 0; inl < deepks_param.inlmax; ++inl)
+        {
+            int nm = 2 * deepks_param.inl2l[inl] + 1;
+            const torch::Tensor& grad_mag = gedm_tensor[deepks_param.inlmax + inl];
+            if (!grad_mag.defined())
+            {
+                for (int index = 0; index < nm * nm; ++index)
+                {
+                    gedm_mag[inl][index] = 0.0;
+                }
+                continue;
+            }
+            auto accessor = grad_mag.accessor<double, 2>();
+            for (int m1 = 0; m1 < nm; ++m1)
+            {
+                for (int m2 = 0; m2 < nm; ++m2)
+                {
+                    gedm_mag[inl][m1 * nm + m2] = accessor[m1][m2] * 2;
+                }
+            }
+        }
+    }
     ModuleBase::timer::end("DeePKS_domain", "cal_edelta_gedm");
     return;
 }
diff --git a/source/source_lcao/module_deepks/deepks_basic.h b/source/source_lcao/module_deepks/deepks_basic.h
index 2a1d3180591..fef75e9179b 100644
--- a/source/source_lcao/module_deepks/deepks_basic.h
+++ b/source/source_lcao/module_deepks/deepks_basic.h
@@ -37,11 +37,14 @@ void cal_gevdm(const int nat,
 /// calculate partial of energy correction to descriptors
 void cal_edelta_gedm(const int nat,
                      const DeePKS_Param& deepks_param,
+                     torch::jit::script::Module& model_deepks,
+                     double& E_delta,
                      const std::vector<torch::Tensor>& descriptor,
                      const std::vector<torch::Tensor>& pdm,
-                     torch::jit::script::Module& model_deepks,
                      double** gedm,
-                     double& E_delta);
+                     const std::vector<torch::Tensor>& descriptor_mag = {},
+                     const std::vector<torch::Tensor>& pdm_mag = {},
+                     double** gedm_mag = nullptr);
 void check_gedm(const DeePKS_Param& deepks_param, double** gedm);
 void cal_edelta_gedm_equiv(const int nat,
                            const DeePKS_Param& deepks_param,
diff --git a/source/source_lcao/module_deepks/deepks_force.cpp b/source/source_lcao/module_deepks/deepks_force.cpp
index eb230df4d8a..a132cb6726a 100644
--- a/source/source_lcao/module_deepks/deepks_force.cpp
+++ b/source/source_lcao/module_deepks/deepks_force.cpp
@@ -11,8 +11,7 @@
 #include "source_lcao/module_hcontainer/atom_pair.h"
 
 template <typename TK>
-void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
-                                const UnitCell& ucell,
+void DeePKS_domain::cal_f_delta(const UnitCell& ucell,
                                 const LCAO_Orbitals& orb,
                                 const Grid_Driver& GridD,
                                 const Parallel_Orbitals& pv,
@@ -20,13 +19,17 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
                                 const DeePKS_Param& deepks_param,
                                 const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                                 std::vector<hamilt::HContainer<double>*> phialpha,
-                                double** gedm,
                                 ModuleBase::matrix& f_delta,
                                 const bool isstress,
-                                ModuleBase::matrix& svnl_dalpha)
+                                ModuleBase::matrix& svnl_dalpha,
+                                const hamilt::HContainer<double>* dmr,
+                                double** gedm,
+                                const hamilt::HContainer<double>* dmr_mag,
+                                double** gedm_mag)
 {
     ModuleBase::TITLE("DeePKS_domain", "cal_f_delta");
     ModuleBase::timer::start("DeePKS_domain", "cal_f_delta");
+    const bool spin2 = (dmr_mag != nullptr && gedm_mag != nullptr); // nspin=2 magnetization channel
     f_delta.zero_out();
 
     const int lmaxd = orb.get_lmax_d();
@@ -108,6 +111,11 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
                 }
                 ModuleBase::Vector3<double> dR(dRx, dRy, dRz);
                 const double* dm_current = dmr->find_matrix(ibt1, ibt2, dR.x, dR.y, dR.z)->get_pointer();
+                const double* dm_current_mag = nullptr;
+                if (spin2)
+                {
+                    dm_current_mag = dmr_mag->find_matrix(ibt1, ibt2, dR.x, dR.y, dR.z)->get_pointer();
+                }
 
                 hamilt::BaseMatrix<double>* overlap_1 = phialpha[0]->find_matrix(iat, ibt1, dR1);
                 hamilt::BaseMatrix<double>* overlap_2 = phialpha[0]->find_matrix(iat, ibt2, dR2);
@@ -131,6 +139,8 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
                     {
                         double nlm[3] = {0, 0, 0};
                         double nlm_t[3] = {0, 0, 0}; // for stress
+                        double nlm_mag[3] = {0, 0, 0};
+                        double nlm_t_mag[3] = {0, 0, 0};
 
                         if (!PARAM.inp.deepks_equiv)
                         {
@@ -158,6 +168,21 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
                                                                     * grad_overlap_1[dim]->get_value(row_indexes[iw1],
                                                                                                     ib + m2);
                                                 }
+                                                if (spin2)
+                                                {
+                                                    nlm_mag[dim]
+                                                        += gedm_mag[inl][m1 * nm + m2]
+                                                           * overlap_1->get_value(row_indexes[iw1], ib + m1)
+                                                           * grad_overlap_2[dim]->get_value(col_indexes[iw2], ib + m2);
+                                                    if (isstress)
+                                                    {
+                                                        nlm_t_mag[dim]
+                                                            += gedm_mag[inl][m1 * nm + m2]
+                                                               * overlap_2->get_value(col_indexes[iw2], ib + m1)
+                                                               * grad_overlap_1[dim]->get_value(row_indexes[iw1],
+                                                                                                ib + m2);
+                                                    }
+                                                }
                                             }
                                         }
                                     }
@@ -193,15 +218,27 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
                             }
                         }
 
-                        // HF term is minus, only one projector for each atom force.
-                        f_delta_local(iat, 0) -= 2.0 * *dm_current * nlm[0];
-                        f_delta_local(iat, 1) -= 2.0 * *dm_current * nlm[1];
-                        f_delta_local(iat, 2) -= 2.0 * *dm_current * nlm[2];
+                        // combine charge and (for nspin=2) magnetization channels
+                        double wnlm[3];
+                        double wnlm_t[3];
+                        for (int dim = 0; dim < 3; ++dim)
+                        {
+                            wnlm[dim] = (*dm_current) * nlm[dim];
+                            wnlm_t[dim] = (*dm_current) * nlm_t[dim];
+                            if (spin2)
+                            {
+                                wnlm[dim] += (*dm_current_mag) * nlm_mag[dim];
+                                wnlm_t[dim] += (*dm_current_mag) * nlm_t_mag[dim];
+                            }
+                        }
 
-                        // Pulay term is plus, only one projector for each atom force.
-                        f_delta_local(ibt2, 0) += 2.0 * *dm_current * nlm[0];
-                        f_delta_local(ibt2, 1) += 2.0 * *dm_current * nlm[1];
-                        f_delta_local(ibt2, 2) += 2.0 * *dm_current * nlm[2];
+                        // HF term is minus, Pulay term is plus
+                        f_delta_local(iat, 0) -= 2.0 * wnlm[0];
+                        f_delta_local(iat, 1) -= 2.0 * wnlm[1];
+                        f_delta_local(iat, 2) -= 2.0 * wnlm[2];
+                        f_delta_local(ibt2, 0) += 2.0 * wnlm[0];
+                        f_delta_local(ibt2, 1) += 2.0 * wnlm[1];
+                        f_delta_local(ibt2, 2) += 2.0 * wnlm[2];
 
                         if (isstress)
                         {
@@ -209,11 +246,14 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
                             {
                                 for (int jpol = ipol; jpol < 3; jpol++)
                                 {
-                                    svnl_dalpha_local(ipol, jpol)
-                                        += *dm_current * (nlm[ipol] * r2[jpol] + nlm_t[ipol] * r1[jpol]);
+                                    svnl_dalpha_local(ipol, jpol) += wnlm[ipol] * r2[jpol] + wnlm_t[ipol] * r1[jpol];
                                 }
                             }
                         }
+                        if (spin2)
+                        {
+                            dm_current_mag++;
+                        }
                         dm_current++;
                     } // iw2
                 }     // iw1
@@ -263,8 +303,7 @@ void DeePKS_domain::cal_f_delta(const hamilt::HContainer<double>* dmr,
     return;
 }
 
-template void DeePKS_domain::cal_f_delta<double>(const hamilt::HContainer<double>* dmr,
-                                                 const UnitCell& ucell,
+template void DeePKS_domain::cal_f_delta<double>(const UnitCell& ucell,
                                                  const LCAO_Orbitals& orb,
                                                  const Grid_Driver& GridD,
                                                  const Parallel_Orbitals& pv,
@@ -272,13 +311,15 @@ template void DeePKS_domain::cal_f_delta<double>(const hamilt::HContainer<double
                                                  const DeePKS_Param& deepks_param,
                                                  const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                                                  std::vector<hamilt::HContainer<double>*> phialpha,
-                                                 double** gedm,
                                                  ModuleBase::matrix& f_delta,
                                                  const bool isstress,
-                                                 ModuleBase::matrix& svnl_dalpha);
+                                                 ModuleBase::matrix& svnl_dalpha,
+                                                 const hamilt::HContainer<double>* dmr,
+                                                 double** gedm,
+                                                 const hamilt::HContainer<double>* dmr_mag,
+                                                 double** gedm_mag);
 
-template void DeePKS_domain::cal_f_delta<std::complex<double>>(const hamilt::HContainer<double>* dmr,
-                                                               const UnitCell& ucell,
+template void DeePKS_domain::cal_f_delta<std::complex<double>>(const UnitCell& ucell,
                                                                const LCAO_Orbitals& orb,
                                                                const Grid_Driver& GridD,
                                                                const Parallel_Orbitals& pv,
@@ -286,9 +327,12 @@ template void DeePKS_domain::cal_f_delta<std::complex<double>>(const hamilt::HCo
                                                                const DeePKS_Param& deepks_param,
                                                                const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                                                                std::vector<hamilt::HContainer<double>*> phialpha,
-                                                               double** gedm,
                                                                ModuleBase::matrix& f_delta,
                                                                const bool isstress,
-                                                               ModuleBase::matrix& svnl_dalpha);
+                                                               ModuleBase::matrix& svnl_dalpha,
+                                                               const hamilt::HContainer<double>* dmr,
+                                                               double** gedm,
+                                                               const hamilt::HContainer<double>* dmr_mag,
+                                                               double** gedm_mag);
 
 #endif
diff --git a/source/source_lcao/module_deepks/deepks_force.h b/source/source_lcao/module_deepks/deepks_force.h
index ef5c0d8a949..4f0335e17c9 100644
--- a/source/source_lcao/module_deepks/deepks_force.h
+++ b/source/source_lcao/module_deepks/deepks_force.h
@@ -25,8 +25,7 @@ namespace DeePKS_domain
 // 1. cal_f_delta, which is used for F_delta calculation
 
 template <typename TK>
-void cal_f_delta(const hamilt::HContainer<double>* dmr,
-                 const UnitCell& ucell,
+void cal_f_delta(const UnitCell& ucell,
                  const LCAO_Orbitals& orb,
                  const Grid_Driver& GridD,
                  const Parallel_Orbitals& pv,
@@ -34,10 +33,13 @@ void cal_f_delta(const hamilt::HContainer<double>* dmr,
                  const DeePKS_Param& deepks_param,
                  const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                  std::vector<hamilt::HContainer<double>*> phialpha,
-                 double** gedm,
                  ModuleBase::matrix& f_delta,
                  const bool isstress,
-                 ModuleBase::matrix& svnl_dalpha);
+                 ModuleBase::matrix& svnl_dalpha,
+                 const hamilt::HContainer<double>* dmr,
+                 double** gedm,
+                 const hamilt::HContainer<double>* dmr_mag = nullptr,
+                 double** gedm_mag = nullptr);
 } // namespace DeePKS_domain
 
 #endif
diff --git a/source/source_lcao/module_deepks/deepks_pdm.cpp b/source/source_lcao/module_deepks/deepks_pdm.cpp
index 04bf4091550..cad6683a72e 100644
--- a/source/source_lcao/module_deepks/deepks_pdm.cpp
+++ b/source/source_lcao/module_deepks/deepks_pdm.cpp
@@ -92,7 +92,9 @@ void DeePKS_domain::update_dmr(const std::vector<ModuleBase::Vector3<double>>& k
                                const LCAO_Orbitals& orb,
                                const Parallel_Orbitals& pv,
                                const Grid_Driver& GridD,
-                               hamilt::HContainer<double>* dmr_deepks)
+                               hamilt::HContainer<double>* dmr_deepks,
+                               const int nspin,
+                               const bool mag)
 {
     dmr_deepks->set_zero();
     // save whether the pair with R has been calculated
@@ -157,6 +159,10 @@ void DeePKS_domain::update_dmr(const std::vector<ModuleBase::Vector3<double>>& k
                     const double arg = -(kvec_d[ik] * ModuleBase::Vector3<double>(dR)) * ModuleBase::TWO_PI;
                     kphase = std::complex<double>(cos(arg), sin(arg));
                 }
+                if (mag && nspin == 2 && ik >= (int)(dmk.size() / nspin))
+                {
+                    kphase *= -1.0; // spin-down block enters with a minus sign
+                }
                 TK* kphase_ptr = reinterpret_cast<TK*>(&kphase);
                 if (ModuleBase::GlobalFunc::IS_COLUMN_MAJOR_KS_SOLVER(PARAM.inp.ks_solver))
                 {
@@ -480,7 +486,9 @@ template void DeePKS_domain::update_dmr<double>(const std::vector<ModuleBase::Ve
                                                 const LCAO_Orbitals& orb,
                                                 const Parallel_Orbitals& pv,
                                                 const Grid_Driver& GridD,
-                                                hamilt::HContainer<double>* dmr_deepks);
+                                                hamilt::HContainer<double>* dmr_deepks,
+                                                const int nspin,
+                                                const bool mag);
 
 template void DeePKS_domain::update_dmr<std::complex<double>>(const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                                                               const std::vector<std::vector<std::complex<double>>>& dmk,
@@ -488,7 +496,9 @@ template void DeePKS_domain::update_dmr<std::complex<double>>(const std::vector<
                                                               const LCAO_Orbitals& orb,
                                                               const Parallel_Orbitals& pv,
                                                               const Grid_Driver& GridD,
-                                                              hamilt::HContainer<double>* dmr_deepks);
+                                                              hamilt::HContainer<double>* dmr_deepks,
+                                                              const int nspin,
+                                                              const bool mag);
 
 template void DeePKS_domain::cal_pdm<double>(bool& init_pdm,
                                              const DeePKS_Param& deepks_param,
diff --git a/source/source_lcao/module_deepks/deepks_pdm.h b/source/source_lcao/module_deepks/deepks_pdm.h
index 3339dc75e90..0a22ef86f7f 100644
--- a/source/source_lcao/module_deepks/deepks_pdm.h
+++ b/source/source_lcao/module_deepks/deepks_pdm.h
@@ -44,7 +44,9 @@ void update_dmr(const std::vector<ModuleBase::Vector3<double>>& kvec_d,
                 const LCAO_Orbitals& orb,
                 const Parallel_Orbitals& pv,
                 const Grid_Driver& GridD,
-                hamilt::HContainer<double>* dmr_deepks);
+                hamilt::HContainer<double>* dmr_deepks,
+                const int nspin = 1,
+                const bool mag = false);
 
 // calculate projected density matrix: pdm = sum_i,occ <phi_i|alpha1><alpha2|phi_k>
 // 3 cases to skip calculation of pdm:
diff --git a/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp b/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp
index 071605bf0a3..38a0cac574d 100644
--- a/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp
+++ b/source/source_lcao/module_deepks/test/LCAO_deepks_test.cpp
@@ -340,11 +340,11 @@ void test_deepks<T>::check_edelta(std::vector<torch::Tensor>& descriptor)
     {
         DeePKS_domain::cal_edelta_gedm(ucell.nat,
                                        this->ld.deepks_param,
+                                       this->ld.model_deepks,
+                                       this->ld.E_delta,
                                        descriptor,
                                        this->ld.pdm,
-                                       this->ld.model_deepks,
-                                       this->ld.gedm,
-                                       this->ld.E_delta);
+                                       this->ld.gedm);
     }
 
     std::ofstream ofs("E_delta.dat");
@@ -399,8 +399,7 @@ void test_deepks<T>::check_f_delta_and_stress_delta()
     svnl_dalpha.create(3, 3);
     const int cal_stress = 1;
     const int nks = kv.nkstot;
-    DeePKS_domain::cal_f_delta<T>(this->ld.dm_r,
-                                  ucell,
+    DeePKS_domain::cal_f_delta<T>(ucell,
                                   ORB,
                                   Test_Deepks::GridD,
                                   ParaO,
@@ -408,10 +407,11 @@ void test_deepks<T>::check_f_delta_and_stress_delta()
                                   this->ld.deepks_param,
                                   kv.kvec_d,
                                   this->ld.phialpha,
-                                  this->ld.gedm,
                                   fvnl_dalpha,
                                   cal_stress,
-                                  svnl_dalpha);
+                                  svnl_dalpha,
+                                  this->ld.dm_r,
+                                  this->ld.gedm);
     std::ofstream ofs_f("F_delta.dat");
     std::ofstream ofs_s("stress_delta.dat");
     ofs_f << std::setprecision(10);
diff --git a/source/source_lcao/module_operator_lcao/deepks_lcao.cpp b/source/source_lcao/module_operator_lcao/deepks_lcao.cpp
index 8c93965cdb0..7255645cbb5 100644
--- a/source/source_lcao/module_operator_lcao/deepks_lcao.cpp
+++ b/source/source_lcao/module_operator_lcao/deepks_lcao.cpp
@@ -147,13 +147,12 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
 {
 #ifdef __MLALGO
     ModuleBase::TITLE("DeePKS", "contributeHR");
+    ModuleBase::timer::start("DeePKS", "contributeHR");
     // if DM changed, HR of DeePKS need to refresh.
     // the judgement is based on the status of HR in ld
     // this operator should be informed that DM has changed and HR need to recalculate.
     if (this->ld->get_hr_cal())
     {
-        ModuleBase::timer::start("DeePKS", "contributeHR");
-
         DeePKS_domain::cal_pdm<TK>(this->ld->init_pdm,
                                    this->ld->deepks_param,
                                    this->kvec_d,
@@ -170,8 +169,27 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
                                       this->ld->deepks_param,
                                       this->ld->pdm,
                                       descriptor);
+        // nspin=2 (traditional): prepare the magnetization-channel descriptor
+        std::vector<torch::Tensor> descriptor_mag;
+        if (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv)
+        {
+            bool init_pdm_mag = false;
+            DeePKS_domain::cal_pdm<TK>(init_pdm_mag,
+                                       this->ld->deepks_param,
+                                       this->kvec_d,
+                                       this->ld->dm_r_mag,
+                                       this->ld->phialpha,
+                                       *this->ucell,
+                                       *ptr_orb_,
+                                       *(this->gd),
+                                       *(this->hR->get_paraV()),
+                                       this->ld->pdm_mag);
+            DeePKS_domain::cal_descriptor(this->ucell->nat, this->ld->deepks_param, this->ld->pdm_mag, descriptor_mag);
+        }
+
         if (PARAM.inp.deepks_equiv)
         {
+            // equivariant version: an independent path; spin is not handled here
             DeePKS_domain::cal_edelta_gedm_equiv(this->ucell->nat,
                                                  this->ld->deepks_param,
                                                  descriptor,
@@ -180,38 +198,56 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::contributeHR()
                                                  this->ld->E_delta,
                                                  GlobalV::MY_RANK);
         }
-        else
+        else // traditional version (descriptor_mag is empty for nspin=1 -> single channel)
         {
             DeePKS_domain::cal_edelta_gedm(this->ucell->nat,
                                            this->ld->deepks_param,
+                                           this->ld->model_deepks,
+                                           this->ld->E_delta,
                                            descriptor,
                                            this->ld->pdm,
-                                           this->ld->model_deepks,
                                            this->ld->gedm,
-                                           this->ld->E_delta);
+                                           descriptor_mag,
+                                           this->ld->pdm_mag,
+                                           this->ld->gedm_mag);
         }
 
-        // // recalculate the V_delta_R
-        // if (this->V_delta_R == nullptr)
-        // {
-        //     this->V_delta_R = new hamilt::HContainer<std::complex<double>>(*this->hR);
-        // }
-        this->V_delta_R->set_zero();
-        this->calculate_HR();
+        // For nspin=1 (and the equivariant version) V_delta_R is spin-independent;
+        // build it once here and reuse it across calls. For nspin=2 it is
+        // spin-dependent and is rebuilt per spin below.
+        if (!(PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv))
+        {
+            this->V_delta_R->set_zero();
+            this->calculate_HR(this->ld->gedm);
+        }
 
         this->ld->set_hr_cal(false);
+    }
 
-        ModuleBase::timer::end("DeePKS", "contributeHR");
+    // nspin=2 (traditional): the correction for spin sigma is
+    // |alpha>(gedm + sigma * gedm_mag)<alpha|; the combination is folded into
+    // calculate_HR, so no per-call gedm buffer is allocated.
+    if (PARAM.inp.nspin == 2 && !PARAM.inp.deepks_equiv)
+    {
+        const double sign = (this->current_spin == 0) ? 1.0 : -1.0;
+        this->V_delta_R->set_zero();
+        this->calculate_HR(this->ld->gedm, this->ld->gedm_mag, sign);
+        this->current_spin = 1 - this->current_spin;
     }
+
     // save V_delta_R to hR
     this->hR->add(*this->V_delta_R);
+
+    ModuleBase::timer::end("DeePKS", "contributeHR");
 #endif
 }
 
 #ifdef __MLALGO
 
 template <typename TK, typename TR>
-void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
+void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR(double** gedm_use,
+                                                                double** gedm_mag_use,
+                                                                const double sign)
 {
     ModuleBase::TITLE("DeePKS", "calculate_HR");
     ModuleBase::timer::start("DeePKS", "calculate_HR");
@@ -245,7 +281,8 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
                 for (int N0 = 0; N0 < ptr_orb_->Alpha[0].getNchi(L0); ++N0)
                 {
                     const int inl = this->ld->deepks_param.inl_index[T0](I0, L0, N0);
-                    const double* pgedm = this->ld->gedm[inl];
+                    const double* pgedm = gedm_use[inl];
+                    const double* pgedm_mag = (gedm_mag_use != nullptr) ? gedm_mag_use[inl] : nullptr;
                     const int nm = 2 * L0 + 1;
 
                     for (int m1 = 0; m1 < nm; ++m1) // m1 = 1 for s, 3 for p, 5 for d
@@ -254,7 +291,8 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
                         {
                             trace_alpha_row.push_back(ib + m1);
                             trace_alpha_col.push_back(ib + m2);
-                            gedms.push_back(pgedm[m1 * nm + m2]);
+                            gedms.push_back(pgedm_mag != nullptr ? pgedm[m1 * nm + m2] + sign * pgedm_mag[m1 * nm + m2]
+                                                                 : pgedm[m1 * nm + m2]);
                         }
                     }
                     ib += nm;
@@ -263,7 +301,7 @@ void hamilt::DeePKS<hamilt::OperatorLCAO<TK, TR>>::calculate_HR()
         }
         else
         {
-            const double* pgedm = this->ld->gedm[iat0];
+            const double* pgedm = gedm_use[iat0];
             int nproj = 0;
             for (int il = 0; il < this->ld->deepks_param.lmaxd + 1; il++)
             {
diff --git a/source/source_lcao/module_operator_lcao/deepks_lcao.h b/source/source_lcao/module_operator_lcao/deepks_lcao.h
index d184e6cb09f..3f790c006f7 100644
--- a/source/source_lcao/module_operator_lcao/deepks_lcao.h
+++ b/source/source_lcao/module_operator_lcao/deepks_lcao.h
@@ -94,7 +94,7 @@ class DeePKS<OperatorLCAO<TK, TR>> : public OperatorLCAO<TK, TR>
      * @brief calculate the DeePKS correction matrix with specific <I,J,R> atom-pairs
      * use the adjs_all to calculate the HR matrix
      */
-    void calculate_HR();
+    void calculate_HR(double** gedm_use, double** gedm_mag_use = nullptr, const double sign = 0.0);
 
     /**
      * @brief calculate the HR local matrix of <I,J,R> atom pair
diff --git a/source/source_lcao/setup_deepks.cpp b/source/source_lcao/setup_deepks.cpp
index 88f84bbc8ad..d5cde5dd0a1 100644
--- a/source/source_lcao/setup_deepks.cpp
+++ b/source/source_lcao/setup_deepks.cpp
@@ -86,6 +86,10 @@ void Setup_DeePKS<TK>::delta_e(const UnitCell& ucell,
     {
         this->ld.dpks_cal_e_delta_band(dm_vec, kv.get_nks());
         DeePKS_domain::update_dmr(kv.kvec_d, dm_vec, ucell, orb, pv, gd, this->ld.dm_r);
+        if (inp.nspin == 2 && !inp.deepks_equiv) // magnetization-channel real-space DM
+        {
+            DeePKS_domain::update_dmr(kv.kvec_d, dm_vec, ucell, orb, pv, gd, this->ld.dm_r_mag, 2, true);
+        }
         f_en.edeepks_scf = this->ld.E_delta - this->ld.e_delta_band;
         f_en.edeepks_delta = this->ld.E_delta;
     }
diff --git a/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/INPUT b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/INPUT
new file mode 100644
index 00000000000..7e31886dc22
--- /dev/null
+++ b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/INPUT
@@ -0,0 +1,35 @@
+INPUT_PARAMETERS 
+#Parameters (1.General)
+suffix            autotest
+calculation       scf
+dft_functional    lda
+
+nbands 8
+symmetry          0
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc           20
+scf_thr 1e-8
+scf_nmax 200
+
+#Parameters (3.Basis)
+basis_type        lcao
+gamma_only        1
+
+#Parameters (4.Smearing)
+smearing_method   gaussian
+smearing_sigma    0.02
+
+#Parameters (5.Mixing)
+mixing_type       broyden
+mixing_beta 0.2
+
+cal_force         1
+cal_stress        1
+deepks_scf        1
+deepks_out_labels  0
+deepks_model      ../Model_ProjOrb/model_nspin2_demo.ptg
+nspin 2
+nupdown 2
diff --git a/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/KPT b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/KPT
new file mode 100644
index 00000000000..c289c0158aa
--- /dev/null
+++ b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0
diff --git a/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/README b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/README
new file mode 100644
index 00000000000..2afda875f2e
--- /dev/null
+++ b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/README
@@ -0,0 +1 @@
+deepks scf nspin=2 (net moment) for gamma_only CH4
diff --git a/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/STRU b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/STRU
new file mode 100644
index 00000000000..e85ed6df0f9
--- /dev/null
+++ b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/STRU
@@ -0,0 +1,33 @@
+ATOMIC_SPECIES
+C 12.000 C_ONCV_PBE-1.0.upf 	#Element, Mass, Pseudopotential
+H 1.008  H_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+C_gga_8au_100Ry_2s2p1d.orb
+H_gga_8au_100Ry_2s1p.orb
+
+NUMERICAL_DESCRIPTOR
+../Model_ProjOrb/2au_20Ry_jle.orb
+
+LATTICE_CONSTANT
+20   			#Lattice constant
+
+LATTICE_VECTORS
+1 0 0
+0 1 0
+0 0 1
+
+ATOMIC_POSITIONS
+Cartesian 		#Cartesian(Unit is LATTICE_CONSTANT)
+C 			#Name of element	
+0.0			#Magnetic for this element.
+1			#Number of atoms
+0.999619 0.000169648 0.00131031
+
+H
+0.0
+4
+0.998911 6.75673e-05 0.0374835
+0.0418908 0.000157703 0.987531
+0.980183 0.965348 0.9869
+0.98072 0.0340401 0.987153
diff --git a/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/result.ref b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/result.ref
new file mode 100644
index 00000000000..21bfb6d4a41
--- /dev/null
+++ b/tests/09_DeePKS/29_NO_GO_deepks_scf_nspin2/result.ref
@@ -0,0 +1,7 @@
+etotref -47.1202591348382853
+etotperatomref -9.4240518270
+totalforceref 1512.567633
+totalstressref 594.856569
+deepks_desc 7.310751
+deepks_dm_eig 31.477251302098573
+totaltimeref 5.69
diff --git a/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/INPUT b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/INPUT
new file mode 100644
index 00000000000..f51cc7d79c3
--- /dev/null
+++ b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/INPUT
@@ -0,0 +1,24 @@
+INPUT_PARAMETERS 
+calculation       scf
+suffix            autotest
+
+ecutwfc           50
+scf_thr 1e-8
+scf_nmax 200
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+basis_type        lcao
+dft_functional    lda
+gamma_only        0
+mixing_type       broyden
+mixing_beta 0.2
+smearing_method   gaussian
+smearing_sigma    0.02
+
+cal_force         1
+cal_stress        1
+deepks_out_labels  0
+deepks_scf        1
+deepks_model      ../Model_ProjOrb/model_nspin2_demo.ptg
+nspin 2
+nupdown 2
diff --git a/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/KPT b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/KPT
new file mode 100644
index 00000000000..7ae3c03fd5c
--- /dev/null
+++ b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+2 1 1 0 0 0.5
\ No newline at end of file
diff --git a/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/README b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/README
new file mode 100644
index 00000000000..6c03a353e20
--- /dev/null
+++ b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/README
@@ -0,0 +1 @@
+deepks scf nspin=2 (net moment) for multi-k system
diff --git a/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/STRU b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/STRU
new file mode 100644
index 00000000000..6b3e51ea3e3
--- /dev/null
+++ b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/STRU
@@ -0,0 +1,31 @@
+ATOMIC_SPECIES
+H 1.00 H_ONCV_PBE-1.0.upf
+O 1.00 O_ONCV_PBE-1.0.upf
+
+LATTICE_CONSTANT
+1
+
+LATTICE_VECTORS
+7 0.0 0.0 
+0.0 7 0.0 
+0.0 0.0 7
+
+ATOMIC_POSITIONS
+Cartesian    # Cartesian(Unit is LATTICE_CONSTANT)
+
+H
+0.0
+2
+-23.196290411877 12.932693996280 6.850457923460 
+-24.006711594306 11.850811692381 9.324876689886
+O
+0.0
+1
+-22.515669775626 12.546754901616 8.620909072323 
+
+NUMERICAL_ORBITAL
+1_H_gga_100Ry_7au_2s1p.orb
+8_O_gga_100Ry_7au_2s2p1d.orb
+
+NUMERICAL_DESCRIPTOR
+../Model_ProjOrb/2au_20Ry_jle.orb
diff --git a/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/result.ref b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/result.ref
new file mode 100644
index 00000000000..3326c3a07a9
--- /dev/null
+++ b/tests/09_DeePKS/30_NO_KP_deepks_scf_nspin2/result.ref
@@ -0,0 +1,7 @@
+etotref -455.8915396166352139
+etotperatomref -151.9638465389
+totalforceref 28.633572
+totalstressref 651.237311
+deepks_desc 2.173539
+deepks_dm_eig 12.023629818086746
+totaltimeref 2.80
diff --git a/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/INPUT b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/INPUT
new file mode 100644
index 00000000000..f86109903f9
--- /dev/null
+++ b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/INPUT
@@ -0,0 +1,36 @@
+INPUT_PARAMETERS 
+#Parameters (1.General)
+suffix            autotest
+calculation       scf
+dft_functional    lda
+
+nbands 8
+symmetry          0
+pseudo_dir        ../../PP_ORB
+orbital_dir       ../../PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc           20
+scf_thr           1e-9
+scf_nmax 200
+
+#Parameters (3.Basis)
+basis_type        lcao
+gamma_only        1
+
+#Parameters (4.Smearing)
+smearing_method   gaussian
+smearing_sigma    0.02
+
+#Parameters (5.Mixing)
+mixing_type       broyden
+mixing_beta 0.2
+
+cal_force         1
+cal_stress        1
+deepks_scf        1
+deepks_out_labels 1
+deepks_model      ../Model_ProjOrb/model_nspin2_demo.ptg
+nspin 2
+deepks_bandgap 1
+nupdown 2
diff --git a/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/KPT b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/KPT
new file mode 100644
index 00000000000..c289c0158aa
--- /dev/null
+++ b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0
diff --git a/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/README b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/README
new file mode 100644
index 00000000000..bf6b9ae9e56
--- /dev/null
+++ b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/README
@@ -0,0 +1 @@
+deepks nspin=2 with output labels (precalc) and bandgap for gamma_only CH4
diff --git a/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/STRU b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/STRU
new file mode 100644
index 00000000000..e85ed6df0f9
--- /dev/null
+++ b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/STRU
@@ -0,0 +1,33 @@
+ATOMIC_SPECIES
+C 12.000 C_ONCV_PBE-1.0.upf 	#Element, Mass, Pseudopotential
+H 1.008  H_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+C_gga_8au_100Ry_2s2p1d.orb
+H_gga_8au_100Ry_2s1p.orb
+
+NUMERICAL_DESCRIPTOR
+../Model_ProjOrb/2au_20Ry_jle.orb
+
+LATTICE_CONSTANT
+20   			#Lattice constant
+
+LATTICE_VECTORS
+1 0 0
+0 1 0
+0 0 1
+
+ATOMIC_POSITIONS
+Cartesian 		#Cartesian(Unit is LATTICE_CONSTANT)
+C 			#Name of element	
+0.0			#Magnetic for this element.
+1			#Number of atoms
+0.999619 0.000169648 0.00131031
+
+H
+0.0
+4
+0.998911 6.75673e-05 0.0374835
+0.0418908 0.000157703 0.987531
+0.980183 0.965348 0.9869
+0.98072 0.0340401 0.987153
diff --git a/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/result.ref b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/result.ref
new file mode 100644
index 00000000000..1706c5c9d9b
--- /dev/null
+++ b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/result.ref
@@ -0,0 +1,18 @@
+etotref -47.1202591336836605
+etotperatomref -9.4240518267
+totalforceref 1512.567631
+totalstressref 594.856569
+deepks_desc 7.310751
+deepks_dm_eig 31.477251301681157
+deepks_e_label 1.7316369631930555
+deepks_edelta 0.001382689937213355
+deepks_o_label 1.6494261075960732
+deepks_odelta 0.0006502451785370678
+deepks_oprec -12.94781242717467
+deepks_f_label 29.414734973926272
+deepks_fdelta 0.000884781727152761
+deepks_s_label 15.943259296616286
+deepks_sdelta 0.0004366517141523793
+deepks_fpre 113.57368130652188
+deepks_spre 46.835304309900394
+totaltimeref 8.83
diff --git a/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/threshold b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/threshold
new file mode 100644
index 00000000000..6b17b7c67f5
--- /dev/null
+++ b/tests/09_DeePKS/31_NO_GO_deepks_bandgap_nspin2/threshold
@@ -0,0 +1 @@
+threshold 1e-5
diff --git a/tests/09_DeePKS/CASES_CPU.txt b/tests/09_DeePKS/CASES_CPU.txt
index 6d0393047d5..73c5880b9b1 100644
--- a/tests/09_DeePKS/CASES_CPU.txt
+++ b/tests/09_DeePKS/CASES_CPU.txt
@@ -26,3 +26,6 @@
 26_NO_KP_deepks_out_freq_elec
 27_NO_GO_deepks_out_2
 28_NO_KP_deepks_out_2
+29_NO_GO_deepks_scf_nspin2
+30_NO_KP_deepks_scf_nspin2
+31_NO_GO_deepks_bandgap_nspin2
diff --git a/tests/09_DeePKS/Model_ProjOrb/model_nspin2_demo.ptg b/tests/09_DeePKS/Model_ProjOrb/model_nspin2_demo.ptg
new file mode 100644
index 00000000000..c417d9ac382
Binary files /dev/null and b/tests/09_DeePKS/Model_ProjOrb/model_nspin2_demo.ptg differ