Refactor module_hamilt global dependencies

source/source_cell/cal_atoms_info.h

@@ -111,24 +111,10 @@ class CalAtomsInfo
             atoms[it].set_index();
         }
 
-        // calculate the total number of local basis
-        // 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 (nspin != 4)
-            {
-                result.nlocal += nlocal_it;
-            }
-            else
-            {
-                result.nlocal += nlocal_it * 2; // zhengdy-soc
-            }
-        }
+        // calculate the total number of local basis. This value is used by cal_nwfc()
+        // to initialize index arrays (iwt2iat, iwt2iw, itia2iat). The formula lives in
+        // unitcell::cal_nlocal() so that GintInfo::init_trace_lo_() shares it.
+        result.nlocal = unitcell::cal_nlocal(atoms, ntype, nspin);
 
         result.nelec = nelec;
         unitcell::cal_nelec(atoms, ntype, result.nelec, nelec_delta);

source/source_cell/cal_nelec_nband.cpp

@@ -8,6 +8,24 @@
 
 namespace unitcell {
 
+int cal_nlocal(const Atom* atoms, const int ntype, const int nspin)
+{
+    int nlocal = 0;
+    for (int it = 0; it < ntype; ++it)
+    {
+        const int nlocal_it = atoms[it].nw * atoms[it].na;
+        if (nspin != 4)
+        {
+            nlocal += nlocal_it;
+        }
+        else
+        {
+            nlocal += nlocal_it * 2; // zhengdy-soc
+        }
+    }
+    return nlocal;
+}
+
 void cal_nelec(const Atom* atoms, const int& ntype, double& nelec, const double nelec_delta)
 {
     ModuleBase::TITLE("UnitCell", "cal_nelec");

source/source_cell/cal_nelec_nband.h

@@ -14,6 +14,29 @@ namespace unitcell {
      */
     void cal_nelec(const Atom* atoms, const int& ntype, double& nelec, const double nelec_delta);
 
+    /**
+     * @brief Calculate the total number of local numerical atomic orbitals.
+     *
+     * nlocal = sum over all atom types of (atoms[it].nw * atoms[it].na).
+     * For nspin == 4 (non-collinear) each basis function carries 2 polarizations,
+     * so nlocal is doubled.
+     *
+     * Shared by cal_atoms_info() (which stores the result in PARAM.globalv.nlocal)
+     * and GintInfo::init_trace_lo_(), so those two can no longer drift apart.
+     * cal_wfc() still repeats the loop inline because it also needs the per-type
+     * prefix sums for Atom::stapos_wf, and asserts its own total against the value
+     * cal_atoms_info() produced.
+     *
+     * @note atoms[it].nw must already be populated, i.e. Atom::set_index() must have
+     *       run for every type before calling this.
+     *
+     * @param atoms [in] atom pointer
+     * @param ntype [in] number of atom types
+     * @param nspin [in] number of spin components
+     * @return total number of local basis functions
+     */
+    int cal_nlocal(const Atom* atoms, const int ntype, const int nspin);
+
     /**
      * @brief Calculate the number of bands.
      *

source/source_esolver/esolver_fp.cpp

@@ -46,6 +46,21 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp)
 
     this->inp_ = &inp;
 
+    SurchemParameters surchem_parameters;
+    surchem_parameters.eb_k = inp.eb_k;
+    surchem_parameters.tau = inp.tau;
+    surchem_parameters.sigma_k = inp.sigma_k;
+    surchem_parameters.nc_k = inp.nc_k;
+    this->solvent.set_parameters(surchem_parameters);
+
+    XCFunctionalParameters xc_parameters;
+    xc_parameters.xc_temperature = inp.xc_temperature;
+    xc_parameters.exx_fock_alpha = inp.exx_fock_alpha;
+    xc_parameters.exx_erfc_alpha = inp.exx_erfc_alpha;
+    xc_parameters.xc_exch_ext = inp.xc_exch_ext;
+    xc_parameters.xc_corr_ext = inp.xc_corr_ext;
+    XC_Functional::set_runtime_parameters(xc_parameters);
+
     ModuleBase::TITLE("ESolver_FP", "before_all_runners");
 
     //! 1) read pseudopotentials
@@ -224,7 +239,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     //! calculate ewald energy
     if (!this->inp_->test_skip_ewald)
     {
-        this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac);
+        this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac, this->inp_->test_energy, GlobalV::ofs_running);
     }
 
     //! set direction of magnetism, used in non-collinear case 

source/source_esolver/esolver_ks_lcao.cpp

@@ -137,7 +137,12 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         this->pw_rho->nx, this->pw_rho->ny, this->pw_rho->nz,
         0, 0, this->pw_big->nbzp_start,
         this->pw_big->nbx, this->pw_big->nby, this->pw_big->nbzp,
-        orb_.Phi, ucell, this->gd));
+        orb_.Phi, ucell, this->gd,
+        this->inp_->nspin,
+        PARAM.globalv.gamma_only_local,
+        PARAM.globalv.domag,
+        this->inp_->device == "gpu",
+        this->inp_->nstream));
     ModuleGint::Gint::set_gint_info(gint_info_.get());
 
     // 7) For each atom, calculate the adjacent atoms in different cells

source/source_esolver/esolver_lr_lcao_bse.cpp

@@ -123,7 +123,12 @@ void ESolver_BSE<T, TR>::before_all_runners(BaseCell& basecell, const Input_para
             this->pw_big->nbzp,
             this->orb_.Phi,
             ucell,
-            this->gd));
+            this->gd,
+            inp.nspin,
+            PARAM.globalv.gamma_only_local,
+            PARAM.globalv.domag,
+            inp.device == "gpu",
+            inp.nstream));
             ModuleGint::Gint::set_gint_info(this->gint_info_.get());
 
     this->pot.resize(this->nspin, nullptr);

source/source_esolver/esolver_lr_lcao_tddft.cpp

@@ -470,7 +470,12 @@ void ModuleESolver::ESolver_LR<T, TR>::initialize_from_unitcell_(UnitCell& ucell
         this->pw_big->nbzp,
         orb.Phi,
         ucell,
-        this->gd));
+        this->gd,
+        this->inp_->nspin,
+        PARAM.globalv.gamma_only_local,
+        PARAM.globalv.domag,
+        this->inp_->device == "gpu",
+        this->inp_->nstream));
     ModuleGint::Gint::set_gint_info(gint_info_.get());
     // if EXX from scratch, init 2-center integral and calculate Cs, Vs 
 #ifdef __EXX

source/source_esolver/lcao_others.cpp

@@ -99,7 +99,12 @@ void ESolver_KS_LCAO<TK, TR>::others(BaseCell& basecell, const int istep)
                                               this->pw_big->nbzp,
                                               orb_.Phi,
                                               ucell,
-                                              this->gd));
+                                              this->gd,
+                                              this->inp_->nspin,
+                                              gamma_only_local,
+                                              PARAM.globalv.domag,
+                                              this->inp_->device == "gpu",
+                                              this->inp_->nstream));
     ModuleGint::Gint::set_gint_info(gint_info_.get());
 
     // (2)For each atom, calculate the adjacent atoms in different cells

source/source_hamilt/module_ewald/h_ewald_pw.cpp

@@ -1,8 +1,7 @@
 #include "h_ewald_pw.h"
+#include "source_base/global_function.h"
 #include "source_base/parallel_comm.h"
-#include "source_io/module_parameter/parameter.h"
 #include "source_base/mymath.h" // use heapsort
-#include "source_io/module_parameter/parameter.h"
 #include "dnrm2.h"
 #include "source_base/parallel_reduce.h"
 #include "source_base/constants.h"
@@ -27,7 +26,9 @@ int H_Ewald_pw::estimate_mxr(const double &rmax, const ModuleBase::Matrix3 &bg)
 
 double H_Ewald_pw::compute_ewald(const UnitCell& cell,
                                  const ModulePW::PW_Basis* rho_basis,
-                                 const ModuleBase::ComplexMatrix& strucFac)
+                                 const ModuleBase::ComplexMatrix& strucFac,
+                                 const int test_energy,
+                                 std::ofstream& output_stream)
 {
     ModuleBase::TITLE("H_Ewald_pw","compute_ewald");
     ModuleBase::timer::start("H_Ewald_pw","compute_ewald");
@@ -73,9 +74,9 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
             charge += cell.atoms[it].na * cell.atoms[it].ncpp.zv;//mohan modify 2007-11-7
         }
     }
-    if(PARAM.inp.test_energy) 
+    if(test_energy)
     {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Total ionic charge",charge);
+        ModuleBase::GlobalFunc::OUT(output_stream,"Total ionic charge",charge);
     }
 
 	// (2) calculate the converged value: alpha
@@ -94,10 +95,10 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
                      erfc(sqrt(cell.tpiba2 * rho_basis->ggecut / 4.0 / alpha));
     }
     while (upperbound > 1.0e-7);
-    if(PARAM.inp.test_energy) 
+    if(test_energy)
     {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"alpha",alpha);
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Upper bound",upperbound);
+        ModuleBase::GlobalFunc::OUT(output_stream,"alpha",alpha);
+        ModuleBase::GlobalFunc::OUT(output_stream,"Upper bound",upperbound);
     }
 
     // G-space sum here.
@@ -123,7 +124,7 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
 	// but that's not the term "gamma_only" I want to use in LCAO,  
 	fact = 1.0;
 
-    //GlobalV::ofs_running << "\n pwb.gstart = " << pwb.gstart << std::endl;
+    //output_stream << "\n pwb.gstart = " << pwb.gstart << std::endl;
     const int ig0 = rho_basis->ig_gge0;
     for (int ig = 0; ig < rho_basis->npw; ig++)
     {
@@ -165,9 +166,9 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
     rmax = 4.0 / sqrt(alpha) / cell.lat0;
     mxr = H_Ewald_pw::estimate_mxr(rmax, cell.G);
 
-    if(PARAM.inp.test_energy) 
+    if(test_energy)
     {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"mxr",mxr);
+        ModuleBase::GlobalFunc::OUT(output_stream,"mxr",mxr);
     }
     std::vector<ModuleBase::Vector3<double>> vec_r(mxr);
     std::vector<double> vec_r2(mxr);
@@ -177,9 +178,9 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
     double* r2 = vec_r2.data();
 
 #ifdef __MPI
-    if(PARAM.inp.test_energy) 
+    if(test_energy)
     {
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"rmax(unit lat0)",rmax);
+        ModuleBase::GlobalFunc::OUT(output_stream,"rmax(unit lat0)",rmax);
     }
 
     int size = 0;
@@ -209,11 +210,11 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
             // calculate tau[na1]-tau[na2]
             dtau = cell.atoms[it1].tau[ia1] - cell.atoms[it2].tau[ia2];
             // generates nearest-neighbors shells
-            H_Ewald_pw::rgen(dtau, rmax, irr, cell.latvec, cell.G, r, r2, mxr, nrm);
+            H_Ewald_pw::rgen(dtau, rmax, irr, cell.latvec, cell.G, r, r2, mxr, nrm, test_energy);
             // at-->cell.latvec, bg-->G
             // and sum to the real space part
 
-            if(PARAM.inp.test_energy>1)
+            if(test_energy>1)
             {
                 ModuleBase::GlobalFunc::OUT("dtau.x",dtau.x);
                 ModuleBase::GlobalFunc::OUT("dtau.y",dtau.y);
@@ -228,7 +229,7 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
                             erfc(sqrt(alpha) * rr) / rr;
                 }
             }
-            if (PARAM.inp.test_energy>1) 
+            if (test_energy>1)
             {
                 ModuleBase::GlobalFunc::OUT("ewaldr",ewaldr);
             }
@@ -237,7 +238,7 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
 #else
     if (rho_basis->ig_gge0 >= 0)
     {	
-        if(PARAM.inp.test_energy) ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"rmax(unit lat0)",rmax);
+        if(test_energy) ModuleBase::GlobalFunc::OUT(output_stream,"rmax(unit lat0)",rmax);
         // with this choice terms up to ZiZj*erfc(4) are counted (erfc(4)=2x10^-8
         int nt1=0;
         int nt2=0;
@@ -253,11 +254,11 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
                         //calculate tau[na]-tau[nb]
                         dtau = cell.atoms[nt1].tau[na] - cell.atoms[nt2].tau[nb];
                         //generates nearest-neighbors shells
-                        H_Ewald_pw::rgen(dtau, rmax, irr, cell.latvec, cell.G, r, r2, mxr, nrm);
+                        H_Ewald_pw::rgen(dtau, rmax, irr, cell.latvec, cell.G, r, r2, mxr, nrm, test_energy);
                         // at-->cell.latvec, bg-->G
                         // and sum to the real space part
 
-                        if (PARAM.inp.test_energy>1)
+                        if (test_energy>1)
                         {
                             ModuleBase::GlobalFunc::OUT("dtau.x",dtau.x);
                             ModuleBase::GlobalFunc::OUT("dtau.y",dtau.y);
@@ -272,7 +273,7 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
                                         erfc(sqrt(alpha) * rr) / rr;
                             }
                         } // enddo
-                        if (PARAM.inp.test_energy>1) ModuleBase::GlobalFunc::OUT("ewaldr",ewaldr);
+                        if (test_energy>1) ModuleBase::GlobalFunc::OUT("ewaldr",ewaldr);
                     } // enddo
                 } // enddo
             } // nt2
@@ -285,7 +286,7 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
 	// mohan fix bug 2010-07-26
     Parallel_Reduce::reduce_pool(ewalds);
 
-    if (PARAM.inp.test_energy>1)
+    if (test_energy>1)
     {
         ModuleBase::GlobalFunc::OUT("ewaldg",ewaldg);
         ModuleBase::GlobalFunc::OUT("ewaldr",ewaldr);
@@ -306,7 +307,8 @@ void H_Ewald_pw::rgen(
     ModuleBase::Vector3<double> *r,
     double *r2,
     const int mxr,
-    int &nrm)
+    int &nrm,
+    const int test_energy)
 {
     //-------------------------------------------------------------------
     // generates neighbours shells (in units of alat) with length
@@ -377,7 +379,7 @@ void H_Ewald_pw::rgen(
 
     nm3 = (int)(dnrm2(3, bg1, 1) * rmax + 2);
 
-    if (PARAM.inp.test_energy>1)
+    if (test_energy>1)
     {
         ModuleBase::GlobalFunc::OUT("nm1",nm1);
         ModuleBase::GlobalFunc::OUT("nm2",nm2);

source/source_hamilt/module_ewald/h_ewald_pw.h

@@ -1,10 +1,11 @@
 #ifndef H_EWALD_PW_H
 #define H_EWALD_PW_H
 
-#include "source_base/global_function.h"
 #include "source_cell/unitcell.h"
 #include "source_basis/module_pw/pw_basis.h"
 
+#include <fstream>
+
 class H_Ewald_pw 
 {
   public:
@@ -14,7 +15,9 @@ class H_Ewald_pw
     // compute the Ewald energy
     static double compute_ewald(const UnitCell& cell,
                                 const ModulePW::PW_Basis* rho_basis,
-                                const ModuleBase::ComplexMatrix& strucFac);
+                                const ModuleBase::ComplexMatrix& strucFac,
+                                int test_energy,
+                                std::ofstream& output_stream);
 
   public:
     static int estimate_mxr(const double &rmax, const ModuleBase::Matrix3 &bg);
@@ -28,7 +31,8 @@ class H_Ewald_pw
         ModuleBase::Vector3<double> *r,
         double *r2,
       const int mxr,
-        int  &nrm
+        int  &nrm,
+        int test_energy
     );
 
 	// the coefficient of ewald method

source/source_hamilt/module_gint/gint.h

@@ -1,4 +1,5 @@
 #pragma once
+#include <cassert>
 #include <memory>
 #include "gint_info.h"
 #include "gint_type.h"
@@ -19,6 +20,14 @@ class Gint
         gint_info_ = gint_info;
     }
 
+    static const GintInfo& get_gint_info()
+    {
+        // set_gint_info() must have been called by the owning ESolver before any
+        // grid integration runs; dereferencing a null gint_info_ here would be UB.
+        assert(gint_info_ != nullptr && "Gint::set_gint_info() has not been called");
+        return *gint_info_;
+    }
+
     protected:
     static GintInfo* gint_info_;
 };