deepmodeling · mohanchen · Dec 18, 2024 · Dec 17, 2024 · Dec 17, 2024 · Dec 17, 2024
diff --git a/source/driver.cpp b/source/driver.cpp
@@ -43,9 +43,6 @@ void Driver::init()
     // (4) close all of the running logs
     ModuleBase::Global_File::close_all_log(GlobalV::MY_RANK, PARAM.inp.out_alllog,PARAM.inp.calculation);
 
-    // (5) output the json file
-    // Json::create_Json(&GlobalC::ucell.symm,GlobalC::ucell.atoms,&INPUT);
-    Json::create_Json(&GlobalC::ucell, PARAM);
 }
 
 void Driver::print_start_info()
@@ -180,7 +177,7 @@ void Driver::atomic_world()
     //--------------------------------------------------
 
     // where the actual stuff is done
-    this->driver_run(GlobalC::ucell);
+    this->driver_run();
 
     ModuleBase::timer::finish(GlobalV::ofs_running);
     ModuleBase::Memory::print_all(GlobalV::ofs_running);

diff --git a/source/driver.h b/source/driver.h
@@ -36,7 +36,7 @@ class Driver
     void atomic_world();
 
     // the actual calculations
-    void driver_run(UnitCell& ucell);
+    void driver_run();
 };
 
 #endif
diff --git a/source/driver_run.cpp b/source/driver_run.cpp
@@ -24,13 +24,12 @@
  * the configuration-changing subroutine takes force and stress and updates the
  * configuration
  */
-void Driver::driver_run(UnitCell& ucell)
+void Driver::driver_run()
 {
     ModuleBase::TITLE("Driver", "driver_line");
     ModuleBase::timer::tick("Driver", "driver_line");
 
     //! 1: setup cell and atom information
-
     // this warning should not be here, mohan 2024-05-22
 #ifndef __LCAO
     if (PARAM.inp.basis_type == "lcao_in_pw" || PARAM.inp.basis_type == "lcao") {
@@ -40,6 +39,13 @@ void Driver::driver_run(UnitCell& ucell)
 #endif
 
     // the life of ucell should begin here, mohan 2024-05-12
+    UnitCell ucell;
+    ucell.setup(PARAM.inp.latname,
+                    PARAM.inp.ntype,
+                    PARAM.inp.lmaxmax,
+                    PARAM.inp.init_vel,
+                    PARAM.inp.fixed_axes);
+
     ucell.setup_cell(PARAM.globalv.global_in_stru, GlobalV::ofs_running);
     Check_Atomic_Stru::check_atomic_stru(ucell, PARAM.inp.min_dist_coef);
 
@@ -86,6 +92,8 @@ void Driver::driver_run(UnitCell& ucell)
 
     ModuleESolver::clean_esolver(p_esolver);
 
+    //! 6: output the json file
+    Json::create_Json(&ucell, PARAM);
     ModuleBase::timer::tick("Driver", "driver_line");
     return;
 }
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
@@ -302,11 +302,11 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
 
     Force_Stress_LCAO<TK> fsl(this->RA, ucell.nat);
 
-    fsl.getForceStress(PARAM.inp.cal_force,
+    fsl.getForceStress(ucell,
+                       PARAM.inp.cal_force,
                        PARAM.inp.cal_stress,
                        PARAM.inp.test_force,
                        PARAM.inp.test_stress,
-                       ucell,
                        this->gd,
                        this->pv,
                        this->pelec,

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp
@@ -31,11 +31,11 @@ Force_Stress_LCAO<T>::~Force_Stress_LCAO()
 {
 }
 template <typename T>
-void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
+void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
+                                          const bool isforce,
                                           const bool isstress,
                                           const bool istestf,
                                           const bool istests,
-                                          const UnitCell& ucell,
                                           const Grid_Driver& gd,
                                           Parallel_Orbitals& pv,
                                           const elecstate::ElecState* pelec,
@@ -830,7 +830,7 @@ void Force_Stress_LCAO<T>::getForceStress(const bool isforce,
 
 // local pseudopotential, ewald, core correction, scc terms in force
 template <typename T>
-void Force_Stress_LCAO<T>::calForcePwPart(const UnitCell& ucell,
+void Force_Stress_LCAO<T>::calForcePwPart(UnitCell& ucell,
                                           ModuleBase::matrix& fvl_dvl,
                                           ModuleBase::matrix& fewalds,
                                           ModuleBase::matrix& fcc,
@@ -857,7 +857,7 @@ void Force_Stress_LCAO<T>::calForcePwPart(const UnitCell& ucell,
     //--------------------------------------------------------
     // force due to core correlation.
     //--------------------------------------------------------
-    f_pw.cal_force_cc(fcc, rhopw, chr, nlpp.numeric, GlobalC::ucell);
+    f_pw.cal_force_cc(fcc, rhopw, chr, nlpp.numeric, ucell);
     //--------------------------------------------------------
     // force due to self-consistent charge.
     //--------------------------------------------------------
@@ -975,7 +975,7 @@ void Force_Stress_LCAO<std::complex<double>>::integral_part(const bool isGammaOn
 
 // vlocal, hartree, ewald, core correction, exchange-correlation terms in stress
 template <typename T>
-void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
+void Force_Stress_LCAO<T>::calStressPwPart(UnitCell& ucell,
                                            ModuleBase::matrix& sigmadvl,
                                            ModuleBase::matrix& sigmahar,
                                            ModuleBase::matrix& sigmaewa,
@@ -1007,7 +1007,7 @@ void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
     //--------------------------------------------------------
     // stress due to core correlation.
     //--------------------------------------------------------
-    sc_pw.stress_cc(sigmacc, rhopw, &sf, 0, nlpp.numeric, chr);
+    sc_pw.stress_cc(sigmacc, rhopw, ucell, &sf, 0, nlpp.numeric, chr);
 
     //--------------------------------------------------------
     // stress due to self-consistent charge.

diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h
@@ -29,11 +29,11 @@ class Force_Stress_LCAO
     Force_Stress_LCAO(Record_adj& ra, const int nat_in);
     ~Force_Stress_LCAO();
 
-    void getForceStress(const bool isforce,
+    void getForceStress(UnitCell& ucell,
+                        const bool isforce,
                         const bool isstress,
                         const bool istestf,
                         const bool istests,
-                        const UnitCell& ucell,
                         const Grid_Driver& gd,
                         Parallel_Orbitals& pv,
                         const elecstate::ElecState* pelec,
@@ -66,7 +66,7 @@ class Force_Stress_LCAO
                        ModuleBase::matrix& fcs, 
                        ModuleSymmetry::Symmetry* symm);
 
-    void calForcePwPart(const UnitCell& ucell,
+    void calForcePwPart(UnitCell& ucell,
                         ModuleBase::matrix& fvl_dvl,
                         ModuleBase::matrix& fewalds,
                         ModuleBase::matrix& fcc,
@@ -105,7 +105,7 @@ class Force_Stress_LCAO
                        const Parallel_Orbitals& pv,
                        const K_Vectors& kv);
 
-    void calStressPwPart(const UnitCell& ucell,
+    void calStressPwPart(UnitCell& ucell,
                          ModuleBase::matrix& sigmadvl,
                          ModuleBase::matrix& sigmahar,
                          ModuleBase::matrix& sigmaewa,

diff --git a/source/module_hamilt_lcao/module_dftu/dftu.cpp b/source/module_hamilt_lcao/module_dftu/dftu.cpp
@@ -59,6 +59,7 @@ void DFTU::init(UnitCell& cell, // unitcell class
     {
         orb_cutoff_ = orb->cutoffs();
     }
+    ucell = &cell;
 #endif
 
     // needs reconstructions in future

diff --git a/source/module_hamilt_lcao/module_dftu/dftu.h b/source/module_hamilt_lcao/module_dftu/dftu.h
@@ -300,6 +300,7 @@ class DFTU
     void set_dmr(const elecstate::DensityMatrix<std::complex<double>, double>* dm_in_dftu_cd);
 
   private:
+    const UnitCell* ucell = nullptr;
     const elecstate::DensityMatrix<double, double>* dm_in_dftu_d = nullptr;
     const elecstate::DensityMatrix<std::complex<double>, double>* dm_in_dftu_cd = nullptr;
 #endif

diff --git a/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp b/source/module_hamilt_lcao/module_dftu/dftu_tools.cpp
@@ -12,48 +12,51 @@ void DFTU::cal_VU_pot_mat_complex(const int spin, const bool newlocale, std::com
     ModuleBase::TITLE("DFTU", "cal_VU_pot_mat_complex");
     ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc);
 
-    for (int it = 0; it < GlobalC::ucell.ntype; ++it)
+    for (int it = 0; it < this->ucell->ntype; ++it)
     {
-        if (PARAM.inp.orbital_corr[it] == -1) {
+        if (PARAM.inp.orbital_corr[it] == -1) 
+        {
             continue;
-}
-        for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++)
+        }
+        for (int ia = 0; ia < this->ucell->atoms[it].na; ia++)
         {
-            const int iat = GlobalC::ucell.itia2iat(it, ia);
-            for (int L = 0; L <= GlobalC::ucell.atoms[it].nwl; L++)
+            const int iat = this->ucell->itia2iat(it, ia);
+            for (int L = 0; L <= this->ucell->atoms[it].nwl; L++)
             {
-                if (L != PARAM.inp.orbital_corr[it]) {
+                if (L != PARAM.inp.orbital_corr[it]) 
+                {
                     continue;
-}
+                }
 
-                for (int n = 0; n < GlobalC::ucell.atoms[it].l_nchi[L]; n++)
+                for (int n = 0; n < this->ucell->atoms[it].l_nchi[L]; n++)
                 {
-                    if (n != 0) {
+                    if (n != 0) 
+                    {
                         continue;
-}
+                    }
 
                     for (int m1 = 0; m1 < 2 * L + 1; m1++)
                     {
                         for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++)
                         {
                             const int mu = this->paraV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]);
-                            if (mu < 0) {
+                            if (mu < 0) 
+                            {
                                 continue;
-}
+                            }
 
                             for (int m2 = 0; m2 < 2 * L + 1; m2++)
                             {
                                 for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++)
                                 {
                                     const int nu
                                         = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]);
-                                    if (nu < 0) {
+                                    if (nu < 0) 
+                                    {
                                         continue;
-}
-
+                                    }
                                     int m1_all = m1 + (2 * L + 1) * ipol1;
                                     int m2_all = m2 + (2 * L + 1) * ipol2;
-
                                     double val = get_onebody_eff_pot(it, iat, L, n, spin, m1_all, m2_all, newlocale);
                                     VU[nu * this->paraV->nrow + mu] = std::complex<double>(val, 0.0);
                                 } // ipol2
@@ -73,43 +76,47 @@ void DFTU::cal_VU_pot_mat_real(const int spin, const bool newlocale, double* VU)
     ModuleBase::TITLE("DFTU", "cal_VU_pot_mat_real");
     ModuleBase::GlobalFunc::ZEROS(VU, this->paraV->nloc);
 
-    for (int it = 0; it < GlobalC::ucell.ntype; ++it)
+    for (int it = 0; it < this->ucell->ntype; ++it)
     {
-        if (PARAM.inp.orbital_corr[it] == -1) {
+        if (PARAM.inp.orbital_corr[it] == -1) 
+        {
             continue;
-}
-        for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ia++)
+        }
+        for (int ia = 0; ia < this->ucell->atoms[it].na; ia++)
         {
-            const int iat = GlobalC::ucell.itia2iat(it, ia);
-            for (int L = 0; L <= GlobalC::ucell.atoms[it].nwl; L++)
+            const int iat = this->ucell->itia2iat(it, ia);
+            for (int L = 0; L <= this->ucell->atoms[it].nwl; L++)
             {
-                if (L != PARAM.inp.orbital_corr[it]) {
+                if (L != PARAM.inp.orbital_corr[it]) 
+                {
                     continue;
-}
+                }
 
-                for (int n = 0; n < GlobalC::ucell.atoms[it].l_nchi[L]; n++)
+                for (int n = 0; n < this->ucell->atoms[it].l_nchi[L]; n++)
                 {
-                    if (n != 0) {
+                    if (n != 0) 
+                    {
                         continue;
-}
-
+                    }
                     for (int m1 = 0; m1 < 2 * L + 1; m1++)
                     {
                         for (int ipol1 = 0; ipol1 < PARAM.globalv.npol; ipol1++)
                         {
                             const int mu = this->paraV->global2local_row(this->iatlnmipol2iwt[iat][L][n][m1][ipol1]);
-                            if (mu < 0) {
+                            if (mu < 0) 
+                            {
                                 continue;
-}
+                            }
                             for (int m2 = 0; m2 < 2 * L + 1; m2++)
                             {
                                 for (int ipol2 = 0; ipol2 < PARAM.globalv.npol; ipol2++)
                                 {
                                     const int nu
                                         = this->paraV->global2local_col(this->iatlnmipol2iwt[iat][L][n][m2][ipol2]);
-                                    if (nu < 0) {
+                                    if (nu < 0) 
+                                    {
                                         continue;
-}
+                                    }
 
                                     int m1_all = m1 + (2 * L + 1) * ipol1;
                                     int m2_all = m2 + (2 * L + 1) * ipol2;
@@ -157,20 +164,21 @@ double DFTU::get_onebody_eff_pot(const int T,
         {
             if (Yukawa)
             {
-                if (m0 == m1) {
+                if (m0 == m1) 
+                {
                     VU = (this->U_Yukawa[T][L][N] - this->J_Yukawa[T][L][N])
                          * (0.5 - this->locale[iat][L][N][spin](m0, m1));
                 } else {
                     VU = -(this->U_Yukawa[T][L][N] - this->J_Yukawa[T][L][N]) * this->locale[iat][L][N][spin](m0, m1);
-}
+                }
             }
             else
             {
                 if (m0 == m1) {
                     VU = (this->U[T]) * (0.5 - this->locale[iat][L][N][spin](m0, m1));
                 } else {
                     VU = -(this->U[T]) * this->locale[iat][L][N][spin](m0, m1);
-}
+                }
             }
         }
         else
@@ -183,15 +191,15 @@ double DFTU::get_onebody_eff_pot(const int T,
                 } else {
                     VU = -(this->U_Yukawa[T][L][N] - this->J_Yukawa[T][L][N])
                          * this->locale_save[iat][L][N][spin](m0, m1);
-}
+                }
             }
             else
             {
                 if (m0 == m1) {
                     VU = (this->U[T]) * (0.5 - this->locale_save[iat][L][N][spin](m0, m1));
                 } else {
                     VU = -(this->U[T]) * this->locale_save[iat][L][N][spin](m0, m1);
-}
+                }
             }
         }
 

diff --git a/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp b/source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp
@@ -79,7 +79,7 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
     stress_loc(ucell,sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge);
 
     // nlcc
-    stress_cc(sigmaxcc, this->rhopw, p_sf, true, locpp.numeric, pelec->charge);
+    stress_cc(sigmaxcc, this->rhopw, ucell, p_sf, true, locpp.numeric, pelec->charge);
 
     // vdw term
     stress_vdw(sigmavdw, ucell);

diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces.cpp
@@ -24,7 +24,7 @@
 #endif
 
 template <typename FPTYPE, typename Device>
-void Forces<FPTYPE, Device>::cal_force(const UnitCell& ucell,
+void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
                                        ModuleBase::matrix& force,
                                        const elecstate::ElecState& elec,
                                        ModulePW::PW_Basis* rho_basis,
@@ -161,15 +161,15 @@ void Forces<FPTYPE, Device>::cal_force(const UnitCell& ucell,
         // DFT+U and DeltaSpin
         if(PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch)
         {
-            this->cal_force_onsite(forceonsite, wg, wfc_basis, GlobalC::ucell, psi_in);
+            this->cal_force_onsite(forceonsite, wg, wfc_basis, ucell, psi_in);
         }
     }
 
     // non-linear core correction
     // not relevant for PAW
     if (!PARAM.inp.use_paw)
     {
-        Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, GlobalC::ucell);
+        Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, ucell);
     }
     else
     {

diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces.h b/source/module_hamilt_pw/hamilt_pwdft/forces.h
@@ -33,7 +33,7 @@ class Forces
     Forces(const int nat_in) : nat(nat_in){};
     ~Forces(){};
 
-    void cal_force(const UnitCell& ucell,
+    void cal_force(UnitCell& ucell,
                    ModuleBase::matrix& force,
                    const elecstate::ElecState& elec,
                    ModulePW::PW_Basis* rho_basis,
-Original file line number
+Diff line change
@@ Expand Up / @@ -59,6 +59,7 @@ void DFTU::init(UnitCell& cell, // unitcell class @@
         {
             orb_cutoff_ = orb->cutoffs();
         }
+        ucell = &cell;
     #endif
         // needs reconstructions in future
@@ Expand Down @@