Revert the dependency between unitcell in source_cell and ORB files in source_basis

python/pyabacus/src/ModuleDriver/py_driver.cpp

@@ -425,7 +425,10 @@ CalculationResult PyDriver::run(
         );
 
         // Read structure
-        impl_->ucell_->setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
+        impl_->ucell_->setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running, PARAM.inp.symmetry_prec, PARAM.inp.dfthalf_type, PARAM.inp.pseudo_dir, PARAM.inp.nspin,
+            PARAM.inp.basis_type, PARAM.inp.orbital_dir, PARAM.inp.init_wfc,
+            PARAM.inp.onsite_radius, PARAM.globalv.deepks_setorb, PARAM.inp.rpa,
+            PARAM.inp.fixed_atoms, PARAM.inp.noncolin, PARAM.inp.calculation, PARAM.inp.esolver_type);
 
         // Check atomic structure
         unitcell::check_atomic_stru(*impl_->ucell_, PARAM.inp.min_dist_coef);

source/Makefile.Objects

@@ -205,6 +205,7 @@ OBJS_CELL=atom_pseudo.o\
     read_atom_species.o\
     sep.o\
     sep_cell.o\
+    cal_nelec_nband.o\
 
 OBJS_DEEPKS=LCAO_deepks.o\
         deepks_basic.o\
@@ -247,7 +248,6 @@ OBJS_ELECSTAT=elecstate.o\
     pot_xc.o\
     cal_ux.o\
     read_orb.o\
-    cal_nelec_nband.o\
     read_pseudo.o\
     cal_wfc.o\
     setup_estate_pw.o\

source/source_cell/CMakeLists.txt

@@ -18,7 +18,6 @@ add_library(
     read_atoms.cpp
     read_atoms_helper.cpp
     read_orb.cpp
-    setup_nonlocal.cpp
     klist.cpp
     parallel_kpoints.cpp
     cell_index.cpp
@@ -34,6 +33,7 @@ add_library(
     sep_cell.cpp
     qlist.cpp
     qlist.h
+    cal_nelec_nband.cpp
 )
 
 if(ENABLE_COVERAGE)

source/source_cell/atom_pseudo.cpp

@@ -1,7 +1,5 @@
 #include "atom_pseudo.h"
-#include "source_io/module_parameter/parameter.h"
 
-#include "source_io/module_parameter/parameter.h"
 Atom_pseudo::Atom_pseudo()
 {
 }
@@ -14,8 +12,12 @@ Atom_pseudo::~Atom_pseudo()
 void Atom_pseudo::set_d_so(ModuleBase::ComplexMatrix& d_so_in,
                            const int& nproj_in,
                            const int& nproj_in_so,
-                           const bool has_so)
+                           const bool has_so,
+                           const bool lspinorb,
+                           const int nspin)
 {
+    const bool lspinorb_ = lspinorb;
+    const int nspin_ = nspin;
     if (this->lmax < -1 || this->lmax > 20)
     {
         ModuleBase::WARNING_QUIT("Numerical_Nonlocal", "bad input of lmax : should be between -1 and 20");
@@ -69,7 +71,7 @@ void Atom_pseudo::set_d_so(ModuleBase::ComplexMatrix& d_so_in,
 
         if (this->lmax > -1)
         {
-            if (PARAM.inp.lspinorb)
+            if (lspinorb_)
             {
                 int is = 0;
                 for (int is1 = 0; is1 < 2; is1++)
@@ -107,7 +109,7 @@ void Atom_pseudo::set_d_so(ModuleBase::ComplexMatrix& d_so_in,
                 {
                     for (int is2 = 0; is2 < 2; is2++)
                     {
-                        if (is >= PARAM.inp.nspin) {
+                        if (is >= nspin_) {
                             break;
 }
                         for (int L1 = 0; L1 < nproj_soc; L1++)

source/source_cell/atom_pseudo.h

@@ -27,7 +27,9 @@ class Atom_pseudo : public pseudo
 		ModuleBase::ComplexMatrix &d_so_in,
 		const int &nproj_in,
 		const int &nproj_in_so,
-		const bool has_so);
+		const bool has_so,
+		const bool lspinorb,
+		const int nspin);
 
 
 	inline void get_d(const int& is, const int& p1, const int& p2, const std::complex<double>*& tmp_d)

source/source_cell/bcast_cell.cpp

@@ -1,6 +1,6 @@
 #include "unitcell.h"
 #include "source_base/parallel_common.h"
-#include "source_io/module_parameter/parameter.h"
+
 #include "source_hamilt/module_xc/exx_info.h" // use GlobalC::exx_info
 
 #include <string>
@@ -99,32 +99,32 @@ namespace unitcell
     #endif
     }
 
-    void bcast_magnetism(Magnetism& magnet, const int ntype)
+    void bcast_magnetism(Magnetism& magnet, const int ntype, const int nspin)
     {
-    #ifdef __MPI
+#ifdef __MPI
         MPI_Barrier(MPI_COMM_WORLD);
         if (GlobalV::MY_RANK != 0)
         {
             magnet.start_mag.resize(ntype, 0.0);
         }
         Parallel_Common::bcast_double(magnet.start_mag.data(), ntype);
-        if (PARAM.inp.nspin == 4) 
+        if (nspin == 4) 
         {
             Parallel_Common::bcast_double(magnet.ux_[0]);
             Parallel_Common::bcast_double(magnet.ux_[1]);
             Parallel_Common::bcast_double(magnet.ux_[2]);
         }
-    #endif
+#endif
     }
 
-    void bcast_unitcell(UnitCell& ucell)
+    void bcast_unitcell(UnitCell& ucell, const int nspin)
     {
-    #ifdef __MPI
+#ifdef __MPI
         const int ntype = ucell.ntype;
         Parallel_Common::bcast_int(ucell.nat);
 
         bcast_Lattice(ucell.lat);
-        bcast_magnetism(ucell.magnet,ntype);
+        bcast_magnetism(ucell.magnet,ntype, nspin);
         bcast_atoms_tau(ucell.atoms,ntype);
 
         for (int i = 0; i < ntype; i++)

source/source_cell/bcast_cell.h

@@ -35,14 +35,17 @@ namespace unitcell
      * @param nytpe: the number of types of the atoms [in]
      */
     void bcast_magnetism(Magnetism& magnet,
-                         const int ntype);
-
+                         const int ntype,
+                         const int nspin);
+
     /**
      * @brief broadcast the unitcell
      * 
      * @param ucell: the unitcell to be broadcasted [in/out]
+     * @param nspin: the number of spin components
      */
-    void bcast_unitcell(UnitCell& ucell);
+    void bcast_unitcell(UnitCell& ucell,
+                        const int nspin);
 
 
 }

source/source_cell/cal_atoms_info.h

@@ -1,90 +1,124 @@
 #ifndef CAL_ATOMS_INFO_H
 #define CAL_ATOMS_INFO_H
-#include "source_io/module_parameter/parameter.h"
-#include "source_estate/cal_nelec_nband.h"
+#include "source_cell/cal_nelec_nband.h"
 #include "source_base/global_function.h"
+
+struct AtomsInfoResult
+{
+    int nlocal = 0;
+    double nelec = 0.0;
+    int nbands = 0;
+    double nupdown = 0.0;
+    bool use_uspp = false;
+    int nbands_l = 0;
+    bool ks_run = false;
+};
+
 class CalAtomsInfo
 {
   public:
     CalAtomsInfo(){};
     ~CalAtomsInfo(){};
 
     /**
-     * @brief Calculate the atom information from pseudopotential to set Parameter
+     * @brief Calculate the atom information from pseudopotential
      *
      * @param atoms [in] Atom pointer
      * @param ntype [in] number of atom types
-     * @param para [out] Parameter object
+     * @param nspin [in] number of spin components
+     * @param two_fermi [in] two fermi level flag
+     * @param nelec_delta [in] electron number delta
+     * @param esolver_type [in] solver type
+     * @param lspinorb [in] spin-orbit coupling flag
+     * @param basis_type [in] basis type
+     * @param smearing_method [in] smearing method
+     * @param ks_solver [in] KS solver type
+     * @param bndpar [in] band parallel parameter
+     * @return AtomsInfoResult containing calculated atom information
      */
-    void cal_atoms_info(const Atom* atoms, const int& ntype, Parameter& para)
+    AtomsInfoResult cal_atoms_info(const Atom* atoms, const int& ntype,
+                                   const int nspin, const bool two_fermi,
+                                   const double nelec_delta,
+                                   const std::string& esolver_type,
+                                   const bool lspinorb,
+                                   const std::string& basis_type,
+                                   const std::string& smearing_method,
+                                   const std::string& ks_solver,
+                                   const int bndpar)
     {
+        AtomsInfoResult result;
+
         // calculate initial total magnetization when NSPIN=2
-        if (para.inp.nspin == 2 && !para.globalv.two_fermi)
+        if (nspin == 2 && !two_fermi)
         {
             for (int it = 0; it < ntype; ++it)
             {
                 for (int ia = 0; ia < atoms[it].na; ++ia)
                 {
-                    para.input.nupdown  += atoms[it].mag[ia];
+                    result.nupdown += atoms[it].mag[ia];
                 }
             }
             GlobalV::ofs_running << std::endl;
-            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "The readin total magnetization", para.inp.nupdown);
+            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "The readin total magnetization", result.nupdown);
         }
 
-
         // decide whether to be USPP
         for (int it = 0; it < ntype; ++it)
         {
             if (atoms[it].ncpp.tvanp)
             {
-                para.sys.use_uspp = true;
+                result.use_uspp = true;
             }
         }
 
         // calculate the total number of local basis
-        para.sys.nlocal = 0;
+        // nlocal = sum over all atom types of (atoms[it].nw * atoms[it].na)
+        // For nspin == 4 (non-collinear), each basis function has 2 polarizations,
+        // so nlocal is doubled. This value is used by cal_nwfc() to initialize
+        // index arrays (iwt2iat, iwt2iw, itia2iat).
+        result.nlocal = 0;
         for (int it = 0; it < ntype; ++it)
         {
             const int nlocal_it = atoms[it].nw * atoms[it].na;
-            if (para.inp.nspin != 4)
+            if (nspin != 4)
             {
-                para.sys.nlocal += nlocal_it;
+                result.nlocal += nlocal_it;
             }
             else
             {
-                para.sys.nlocal += nlocal_it * 2; // zhengdy-soc
+                result.nlocal += nlocal_it * 2; // zhengdy-soc
             }
         }
 
         // calculate the total number of electrons
-        elecstate::cal_nelec(atoms, ntype, para.input.nelec);
+        unitcell::cal_nelec(atoms, ntype, result.nelec, nelec_delta);
 
         // autoset and check GlobalV::NBANDS
         std::vector<double> nelec_spin(2, 0.0);
-        if (para.inp.nspin == 2)
+        if (nspin == 2)
         {
-            nelec_spin[0] = (para.inp.nelec + para.inp.nupdown ) / 2.0;
-            nelec_spin[1] = (para.inp.nelec - para.inp.nupdown ) / 2.0;
+            nelec_spin[0] = (result.nelec + result.nupdown) / 2.0;
+            nelec_spin[1] = (result.nelec - result.nupdown) / 2.0;
         }
-        elecstate::cal_nbands(para.inp.nelec, para.sys.nlocal, nelec_spin, para.input.nbands);
+        unitcell::cal_nbands(static_cast<int>(result.nelec), result.nlocal, nelec_spin, result.nbands,
+                              esolver_type, lspinorb, nspin, basis_type, smearing_method);
 
         // calculate the number of nbands_local
-        para.sys.nbands_l = para.inp.nbands;
-        if (para.inp.ks_solver == "bpcg") // only bpcg support band parallel
+        result.nbands_l = result.nbands;
+        if (ks_solver == "bpcg")
         {
-            para.sys.nbands_l = para.inp.nbands / para.inp.bndpar;
-            if (GlobalV::MY_BNDGROUP < para.inp.nbands % para.inp.bndpar)
+            result.nbands_l = result.nbands / bndpar;
+            if (GlobalV::MY_BNDGROUP < result.nbands % bndpar)
             {
-                para.sys.nbands_l++;
+                result.nbands_l++;
             }
         }
         // temporary code
-        if (GlobalV::MY_BNDGROUP == 0 || para.inp.ks_solver == "bpcg")
+        if (GlobalV::MY_BNDGROUP == 0 || ks_solver == "bpcg")
         {
-            para.sys.ks_run = true;
+            result.ks_run = true;
         }
-        return;
+        return result;
     }
 };
 #endif