Merge branch 'hamilt' of https://github.com/ErjieWu/abacus-develop into hamilt

Author: ErjieWu

source/module_elecstate/elecstate_pw.cpp

@@ -263,7 +263,7 @@ void ElecStatePW<T, Device>::add_usrho(const psi::Psi<T, Device>& psi)
         // get |beta>
         if (this->ppcell->nkb > 0)
         {
-            this->ppcell->getvnl(this->ctx, ik, this->vkb);
+            this->ppcell->getvnl(this->ctx, *ucell,ik, this->vkb);
         }
 
         // becp = <beta|psi>

source/module_elecstate/test/elecstate_pw_test.cpp

@@ -115,12 +115,14 @@ std::complex<double>* pseudopot_cell_vnl::get_vkb_data<double>() const
 }
 template <>
 void pseudopot_cell_vnl::getvnl<float, base_device::DEVICE_CPU>(base_device::DEVICE_CPU*,
+                                                                const UnitCell&,
                                                                 int const&,
                                                                 std::complex<float>*) const
 {
 }
 template <>
 void pseudopot_cell_vnl::getvnl<double, base_device::DEVICE_CPU>(base_device::DEVICE_CPU*,
+                                                                 const UnitCell&,
                                                                  int const&,
                                                                  std::complex<double>*) const
 {

source/module_esolver/esolver_fp.cpp

@@ -281,7 +281,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
             this->pw_rhod->collect_uniqgg();
         }
 
-        this->p_locpp->init_vloc(this->pw_rhod);
+        this->p_locpp->init_vloc(ucell,this->pw_rhod);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
         this->pelec->omega = ucell.omega;

source/module_esolver/esolver_gets.cpp

@@ -136,8 +136,8 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
     if (PARAM.inp.out_mat_r)
     {
         cal_r_overlap_R r_matrix;
-        r_matrix.init(pv, orb_);
-        r_matrix.out_rR(istep);
+        r_matrix.init(ucell,pv, orb_);
+        r_matrix.out_rR(ucell,istep);
     }
 
     ModuleBase::timer::tick("ESolver_GetS", "runner");

source/module_esolver/esolver_ks_lcao.cpp

@@ -197,7 +197,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     }
 
     // 8) initialize ppcell
-    this->ppcell.init_vloc(this->pw_rho);
+    this->ppcell.init_vloc(ucell,this->pw_rho);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
     // 9) inititlize the charge density
@@ -920,6 +920,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
                         PARAM.inp.out_dm1,
                         false,
                         PARAM.inp.out_app_flag,
+                        ucell.get_iat2iwt(),
+                        &ucell.nat,
                         istep);
 
     //! 4) write density matrix
@@ -1228,7 +1230,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
                                              PARAM.inp.nnkpfile,
                                              PARAM.inp.wannier_spin);
 
-            myWannier.calculate(this->pelec->ekb, 
+            myWannier.calculate(ucell,
+                                this->pelec->ekb, 
                                 this->pw_wfc, 
                                 this->pw_big, 
                                 this->sf, 
@@ -1247,7 +1250,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
                                        PARAM.inp.wannier_spin,
                                        orb_);
 
-            myWannier.calculate(this->pelec->ekb, this->kv, *(this->psi), &(this->pv));
+            myWannier.calculate(ucell,this->pelec->ekb, this->kv, *(this->psi), &(this->pv));
         }
         std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wave function to Wannier90");
     }
@@ -1259,12 +1262,13 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep)
     {
         std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase calculation");
         berryphase bp(&(this->pv));
-        bp.lcao_init(this->kv,
+        bp.lcao_init(ucell,
+                     this->kv,
                      this->GridT,
                      orb_); // additional step before calling
                             // macroscopic_polarization (why capitalize
                             // the function name?)
-        bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv);
+        bp.Macroscopic_polarization(ucell,this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv);
         std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase calculation");
     }
 }

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -287,15 +287,16 @@ void ESolver_KS_LCAO_TDDFT::after_scf(UnitCell& ucell, const int istep)
         {
             std::stringstream ss_dipole;
             ss_dipole << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_DIPOLE";
-            ModuleIO::write_dipole(pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str());
+            ModuleIO::write_dipole(ucell,pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str());
         }
     }
     if (TD_Velocity::out_current == true)
     {
         elecstate::DensityMatrix<std::complex<double>, double>* tmp_DM
             = dynamic_cast<elecstate::ElecStateLCAO<std::complex<double>>*>(this->pelec)->get_DM();
 
-        ModuleIO::write_current(istep,
+        ModuleIO::write_current(ucell,
+                                istep,
                                 this->psi,
                                 pelec,
                                 kv,

source/module_esolver/esolver_ks_lcaopw.cpp

@@ -62,9 +62,9 @@ namespace ModuleESolver
     }
 
     template <typename T>
-    void ESolver_KS_LIP<T>::allocate_hamilt()
+    void ESolver_KS_LIP<T>::allocate_hamilt(const UnitCell& ucell)
     {
-        this->p_hamilt = new hamilt::HamiltLIP<T>(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell
+        this->p_hamilt = new hamilt::HamiltLIP<T>(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell, &ucell
 #ifdef __EXX
             , *this->exx_lip
 #endif

source/module_esolver/esolver_ks_lcaopw.h

@@ -33,7 +33,7 @@ namespace ModuleESolver
                                            const int iter,
                                            const double ethr) override;
 
-        virtual void allocate_hamilt() override;
+        virtual void allocate_hamilt(const UnitCell& ucell) override;
         virtual void deallocate_hamilt() override;
 
 #ifdef __EXX

source/module_esolver/esolver_ks_pw.cpp

@@ -115,9 +115,9 @@ ESolver_KS_PW<T, Device>::~ESolver_KS_PW()
 }
 
 template <typename T, typename Device>
-void ESolver_KS_PW<T, Device>::allocate_hamilt()
+void ESolver_KS_PW<T, Device>::allocate_hamilt(const UnitCell& ucell)
 {
-    this->p_hamilt = new hamilt::HamiltPW<T, Device>(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell);
+    this->p_hamilt = new hamilt::HamiltPW<T, Device>(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell, &ucell);
 }
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::deallocate_hamilt()
@@ -202,10 +202,10 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
     }
 
     //! init pseudopotential
-    this->ppcell.init(ucell.ntype, &this->sf, this->pw_wfc);
+    this->ppcell.init(ucell,&this->sf, this->pw_wfc);
 
     //! initalize local pseudopotential
-    this->ppcell.init_vloc(this->pw_rhod);
+    this->ppcell.init_vloc(ucell,this->pw_rhod);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
     //! Initalize non-local pseudopotential
@@ -279,7 +279,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     this->deallocate_hamilt();
 
     // allocate HamiltPW
-    this->allocate_hamilt();
+    this->allocate_hamilt(ucell);
 
     //----------------------------------------------------------
     // about vdw, jiyy add vdwd3 and linpz add vdwd2
@@ -585,8 +585,8 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep)
                            PARAM.inp.out_wannier_wvfn_formatted,
                            PARAM.inp.nnkpfile,
                            PARAM.inp.wannier_spin);
-
-        wan.calculate(this->pelec->ekb, this->pw_wfc, this->pw_big, this->kv, this->psi);
+        wan.set_tpiba_omega(ucell.tpiba, ucell.omega);
+        wan.calculate(ucell,this->pelec->ekb, this->pw_wfc, this->pw_big, this->kv, this->psi);
         std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Wannier functions calculation");
     }
 
@@ -595,7 +595,7 @@ void ESolver_KS_PW<T, Device>::after_scf(UnitCell& ucell, const int istep)
     {
         std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "Berry phase polarization");
         berryphase bp;
-        bp.Macroscopic_polarization(this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv);
+        bp.Macroscopic_polarization(ucell,this->pw_wfc->npwk_max, this->psi, this->pw_rho, this->pw_wfc, this->kv);
         std::cout << FmtCore::format(" >> Finish %s.\n * * * * * *\n", "Berry phase polarization");
     }
 }
@@ -621,7 +621,8 @@ void ESolver_KS_PW<T, Device>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                            : reinterpret_cast<psi::Psi<std::complex<double>, Device>*>(this->kspw_psi);
 
     // Calculate forces
-    ff.cal_force(force,
+    ff.cal_force(ucell,
+                 force,
                  *this->pelec,
                  this->pw_rhod,
                  &ucell.symm,
@@ -782,7 +783,7 @@ void ESolver_KS_PW<T, Device>::after_all_runners(UnitCell& ucell)
     //! 6) Print out electronic wave functions in real space
     if (PARAM.inp.out_wfc_r == 1) // Peize Lin add 2021.11.21
     {
-        ModuleIO::write_psi_r_1(this->psi[0], this->pw_wfc, "wfc_realspace", true, this->kv);
+        ModuleIO::write_psi_r_1(ucell,this->psi[0], this->pw_wfc, "wfc_realspace", true, this->kv);
     }
 
     //! 7) Use Kubo-Greenwood method to compute conductivities

source/module_esolver/esolver_ks_pw.h

@@ -46,7 +46,7 @@ class ESolver_KS_PW : public ESolver_KS<T, Device>
 
     virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr) override;
 
-    virtual void allocate_hamilt();
+    virtual void allocate_hamilt(const UnitCell& ucell);
     virtual void deallocate_hamilt();
 
     //! hide the psi in ESolver_KS for tmp use