Feature: enable get_s for nspin=4 (#5614)

Author: YuLiu98

source/module_esolver/esolver.cpp

@@ -186,39 +186,29 @@ ESolver* init_esolver(const Input_para& inp, UnitCell& ucell)
         }
     }
     else if (esolver_type == "ksdft_lcao")
-	{
-		if (PARAM.globalv.gamma_only_local)
-		{
-            if (PARAM.inp.calculation == "get_S")
+    {
+        if (PARAM.inp.calculation == "get_S")
+        {
+            if (PARAM.globalv.gamma_only_local)
             {
                 ModuleBase::WARNING_QUIT("ESolver", "get_S is not implemented for gamma_only");
             }
             else
             {
-                return new ESolver_KS_LCAO<double, double>();
+                return new ESolver_GetS();
             }
         }
+        if (PARAM.globalv.gamma_only_local)
+        {
+            return new ESolver_KS_LCAO<double, double>();
+        }
         else if (PARAM.inp.nspin < 4)
         {
-            if (PARAM.inp.calculation == "get_S")
-            {
-                return new ESolver_GetS();
-            }
-            else
-            {
-                return new ESolver_KS_LCAO<std::complex<double>, double>();
-            }
+            return new ESolver_KS_LCAO<std::complex<double>, double>();
         }
         else
         {
-            if (PARAM.inp.calculation == "get_S")
-            {
-                ModuleBase::WARNING_QUIT("ESolver", "get_S is not implemented for npsin=4");
-            }
-            else
-            {
-                return new ESolver_KS_LCAO<std::complex<double>, std::complex<double>>();
-            }
+            return new ESolver_KS_LCAO<std::complex<double>, std::complex<double>>();
         }
     }
     else if (esolver_type == "ksdft_lcao_tddft")

source/module_esolver/esolver_gets.cpp

@@ -103,11 +103,24 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
 
     if (this->p_hamilt == nullptr)
     {
-        this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(&this->pv,
-                                                                              this->kv,
-                                                                              *(two_center_bundle_.overlap_orb),
-                                                                              orb_.cutoffs());
-        dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, double>*>(this->p_hamilt->ops)->contributeHR();
+        if (PARAM.inp.nspin == 4)
+        {
+            this->p_hamilt
+                = new hamilt::HamiltLCAO<std::complex<double>, std::complex<double>>(&this->pv,
+                                                                                     this->kv,
+                                                                                     *(two_center_bundle_.overlap_orb),
+                                                                                     orb_.cutoffs());
+            dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>*>(this->p_hamilt->ops)
+                ->contributeHR();
+        }
+        else
+        {
+            this->p_hamilt = new hamilt::HamiltLCAO<std::complex<double>, double>(&this->pv,
+                                                                                  this->kv,
+                                                                                  *(two_center_bundle_.overlap_orb),
+                                                                                  orb_.cutoffs());
+            dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, double>*>(this->p_hamilt->ops)->contributeHR();
+        }
     }
 
     const std::string fn = PARAM.globalv.global_out_dir + "SR.csr";

source/module_io/write_HS_R.cpp

@@ -171,14 +171,28 @@ void ModuleIO::output_SR(Parallel_Orbitals& pv,
 
     const int istep = 0;
 
-    ModuleIO::save_sparse(HS_Arrays.SR_sparse,
-                          HS_Arrays.all_R_coor,
-                          sparse_thr,
-                          binary,
-                          SR_filename,
-                          pv,
-                          "S",
-                          istep);
+    if (PARAM.inp.nspin == 4)
+    {
+        ModuleIO::save_sparse(HS_Arrays.SR_soc_sparse,
+                              HS_Arrays.all_R_coor,
+                              sparse_thr,
+                              binary,
+                              SR_filename,
+                              pv,
+                              "S",
+                              istep);
+    }
+    else
+    {
+        ModuleIO::save_sparse(HS_Arrays.SR_sparse,
+                              HS_Arrays.all_R_coor,
+                              sparse_thr,
+                              binary,
+                              SR_filename,
+                              pv,
+                              "S",
+                              istep);
+    }
 
     sparse_format::destroy_HS_R_sparse(HS_Arrays);