Refactor:Remove GlobalC::ucell in module_pwdft

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_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

source/module_esolver/esolver_ks_pw.cpp

@@ -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.ntype, 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
@@ -619,7 +619,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,

source/module_esolver/esolver_of.cpp

@@ -112,7 +112,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS");
 
     // initialize local pseudopotential
-    this->locpp.init_vloc(pw_rho);
+    this->locpp.init_vloc(ucell,pw_rho);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
 
@@ -538,7 +538,7 @@ double ESolver_OF::cal_energy()
 void ESolver_OF::cal_force(UnitCell& ucell, ModuleBase::matrix& force)
 {
     Forces<double> ff(ucell.nat);
-    ff.cal_force(force, *pelec, this->pw_rho, &ucell.symm, &sf, &this->locpp);
+    ff.cal_force(ucell,force, *pelec, this->pw_rho, &ucell.symm, &sf, &this->locpp);
 }
 
 /**

source/module_esolver/esolver_sdft_pw.cpp

@@ -171,6 +171,7 @@ void ESolver_SDFT_PW<T, Device>::hamilt2density_single(UnitCell& ucell, int iste
                               this->psi[0],
                               this->pelec,
                               this->pw_wfc,
+                              ucell,
                               this->stowf,
                               istep,
                               iter,

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -869,11 +869,11 @@ void Force_Stress_LCAO<T>::calForcePwPart(const UnitCell& ucell,
     // local pseudopotential force:
     // use charge density; plane wave; local pseudopotential;
     //--------------------------------------------------------
-    f_pw.cal_force_loc(fvl_dvl, rhopw, nlpp.vloc, chr);
+    f_pw.cal_force_loc(ucell,fvl_dvl, rhopw, nlpp.vloc, chr);
     //--------------------------------------------------------
     // ewald force: use plane wave only.
     //--------------------------------------------------------
-    f_pw.cal_force_ew(fewalds, rhopw, &sf); // remain problem
+    f_pw.cal_force_ew(ucell,fewalds, rhopw, &sf); // remain problem
 
     //--------------------------------------------------------
     // force due to core correlation.
@@ -1009,17 +1009,17 @@ void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
     // local pseudopotential stress:
     // use charge density; plane wave; local pseudopotential;
     //--------------------------------------------------------
-    sc_pw.stress_loc(sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr);
+    sc_pw.stress_loc(ucell,sigmadvl, rhopw, nlpp.vloc, &sf, 0, chr);
 
     //--------------------------------------------------------
     // hartree term
     //--------------------------------------------------------
-    sc_pw.stress_har(sigmahar, rhopw, 0, chr);
+    sc_pw.stress_har(ucell,sigmahar, rhopw, 0, chr);
 
     //--------------------------------------------------------
     // ewald stress: use plane wave only.
     //--------------------------------------------------------
-    sc_pw.stress_ewa(sigmaewa, rhopw, 0); // remain problem
+    sc_pw.stress_ewa(ucell,sigmaewa, rhopw, 0); // remain problem
 
     //--------------------------------------------------------
     // stress due to core correlation.
@@ -1034,7 +1034,7 @@ void Force_Stress_LCAO<T>::calStressPwPart(const UnitCell& ucell,
         sigmaxc(i, i) = -etxc / ucell.omega;
     }
     // Exchange-correlation for PBE
-    sc_pw.stress_gga(sigmaxc, rhopw, chr);
+    sc_pw.stress_gga(ucell,sigmaxc, rhopw, chr);
 
     return;
 }

source/module_hamilt_pw/hamilt_ofdft/of_stress_pw.cpp

@@ -63,20 +63,20 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot,
     }
 
     // hartree contribution
-    stress_har(sigmahar, this->rhopw, true, pelec->charge);
+    stress_har(ucell,sigmahar, this->rhopw, true, pelec->charge);
 
     // ewald contribution
-    stress_ewa(sigmaewa, this->rhopw, true);
+    stress_ewa(ucell,sigmaewa, this->rhopw, true);
 
     // xc contribution: add gradient corrections(non diagonal)
     for (int i = 0; i < 3; i++)
     {
         sigmaxc(i, i) = -(pelec->f_en.etxc - pelec->f_en.vtxc) / ucell.omega;
     }
-    stress_gga(sigmaxc, this->rhopw, pelec->charge);
+    stress_gga(ucell,sigmaxc, this->rhopw, pelec->charge);
 
     // local contribution
-    stress_loc(sigmaloc, this->rhopw, locpp.vloc, p_sf, true, pelec->charge);
+    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);

source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp

@@ -17,7 +17,8 @@ pseudopot_cell_vl::~pseudopot_cell_vl()
 	delete[] zp;
 }
 
-void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis)
+void pseudopot_cell_vl::init_vloc(const UnitCell& ucell,
+								  const ModulePW::PW_Basis* rho_basis)
 {
 	if(PARAM.inp.use_paw) { return;
 }
@@ -30,19 +31,19 @@ void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis)
 	double *vloc1d = new double[rho_basis->ngg];
 	ModuleBase::GlobalFunc::ZEROS(vloc1d, rho_basis->ngg);
 
-	this->allocate(rho_basis->ngg);
+	this->allocate(ucell,rho_basis->ngg);
 
-	for (int it = 0; it < GlobalC::ucell.ntype; it++) 
+	for (int it = 0; it < ucell.ntype; it++) 
 	{
-		const Atom* atom = &GlobalC::ucell.atoms[it];
+		const Atom* atom = &ucell.atoms[it];
 
 		ModuleBase::GlobalFunc::ZEROS(vloc1d, rho_basis->ngg);
 
 		this->zp[it] = atom->ncpp.zv;
 		// compute V_loc(G) for a given type of atom
 		if(atom->coulomb_potential)
 		{
-			this->vloc_coulomb(this->zp[it], vloc1d, rho_basis);
+			this->vloc_coulomb(ucell,this->zp[it], vloc1d, rho_basis);
 		}
 		else if(numeric[it]==true)
 		{
@@ -53,6 +54,7 @@ void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis)
 					atom->ncpp.vloc_at.data(), // local potential in real space radial form.  
 		          	this->zp[it],
 					vloc1d,
+					ucell,
 					rho_basis);
 		}
 		else
@@ -69,26 +71,27 @@ void pseudopot_cell_vl::init_vloc(const ModulePW::PW_Basis* rho_basis)
 
 	delete[] vloc1d;
 
-	this->print_vloc(rho_basis);
+	this->print_vloc(ucell,rho_basis);
 
 	ModuleBase::timer::tick("ppcell_vl","init_vloc");
 	return;
 }
 
 
-void pseudopot_cell_vl::allocate(const int ngg)
+void pseudopot_cell_vl::allocate(const UnitCell& ucell,
+							     const int ngg)
 {
 	if(PARAM.inp.test_pp>0) { ModuleBase::TITLE("pseudopot_cell_vl","allocate");
 }
 	if(PARAM.inp.use_paw) { return;
 }
-	this->vloc.create(GlobalC::ucell.ntype, ngg);
+	this->vloc.create(ucell.ntype, ngg);
 
 	delete[] numeric;
-	this->numeric = new bool[GlobalC::ucell.ntype];
-	ModuleBase::GlobalFunc::ZEROS(numeric, GlobalC::ucell.ntype);
+	this->numeric = new bool[ucell.ntype];
+	ModuleBase::GlobalFunc::ZEROS(numeric, ucell.ntype);
 
-	for (int it = 0; it < GlobalC::ucell.ntype; it++)
+	for (int it = 0; it < ucell.ntype; it++)
 	{ 
 		this->numeric[it] = true; 
 	}
@@ -104,7 +107,10 @@ void pseudopot_cell_vl::allocate(const int ngg)
 	return;
 }
 
-void pseudopot_cell_vl::vloc_coulomb(const double& zp_in, double* vloc_1d, const ModulePW::PW_Basis* rho_basis) const
+void pseudopot_cell_vl::vloc_coulomb(const UnitCell& ucell,
+									 const double& zp_in, 
+									 double* vloc_1d, 
+									 const ModulePW::PW_Basis* rho_basis) const
 {
     int igl0 = 0;
     // start from |G|=0 or not.
@@ -117,14 +123,14 @@ void pseudopot_cell_vl::vloc_coulomb(const double& zp_in, double* vloc_1d, const
     {
         igl0 = 0;
     }
-    const double d_fpi_omega = ModuleBase::FOUR_PI / GlobalC::ucell.omega; // mohan add 2008-06-04
+    const double d_fpi_omega = ModuleBase::FOUR_PI / ucell.omega; // mohan add 2008-06-04
     double fac = -zp_in * ModuleBase::e2 * d_fpi_omega;
 #ifdef _OPENMP
 #pragma omp for
 #endif
     for (int ig = igl0; ig < rho_basis->ngg; ig++)
     {
-        double gx2 = rho_basis->gg_uniq[ig] * GlobalC::ucell.tpiba2;
+        double gx2 = rho_basis->gg_uniq[ig] * ucell.tpiba2;
         vloc_1d[ig] = fac / gx2;
     }
     return;
@@ -145,6 +151,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh,
                                   const double* vloc_at,
                                   const double& zp_in,
                                   double* vloc_1d,
+								  const UnitCell& ucell,
                                   const ModulePW::PW_Basis* rho_basis) const
 {
 	//----------------------------------------------------------------
@@ -173,7 +180,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh,
 	/*
 	for(ir=0; ir<msh; ir++)
 	{
-		aux[ir] = r[ir] * zp_in * e2 / GlobalC::ucell.omega;
+		aux[ir] = r[ir] * zp_in * e2 / ucell.omega;
 	}
 	ModuleBase::Integral::Simpson_Integral(msh, aux, rab, vloc_1d[0] );
 	vloc_1d[0] *= 4*3.1415926;
@@ -224,7 +231,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh,
 #endif
 	for (int ig = igl0;ig < rho_basis->ngg;ig++) 
 	{
-		double gx2= rho_basis->gg_uniq[ig] * GlobalC::ucell.tpiba2;
+		double gx2= rho_basis->gg_uniq[ig] * ucell.tpiba2;
 		double gx = std::sqrt(gx2);
 		for (int ir = 0;ir < msh;ir++) 
 		{
@@ -235,7 +242,7 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh,
 		vloc_1d[ig] -= fac * ModuleBase::libm::exp(- gx2 * 0.25)/ gx2;
 	} // enddo
 
-	const double d_fpi_omega = ModuleBase::FOUR_PI/GlobalC::ucell.omega;//mohan add 2008-06-04
+	const double d_fpi_omega = ModuleBase::FOUR_PI/ucell.omega;//mohan add 2008-06-04
 #ifdef _OPENMP
 #pragma omp for
 #endif
@@ -253,21 +260,22 @@ void pseudopot_cell_vl::vloc_of_g(const int& msh,
 	return;
 } // end subroutine vloc_of_g
 
-void pseudopot_cell_vl::print_vloc(const ModulePW::PW_Basis* rho_basis) const
+void pseudopot_cell_vl::print_vloc(const UnitCell& ucell,
+								   const ModulePW::PW_Basis* rho_basis) const
 {
 	if(GlobalV::MY_RANK!=0) { return; //mohan fix bug 2011-10-13
 }
 	bool check_vl = PARAM.inp.out_element_info;
 	if(check_vl)
 	{
-		for(int it=0; it<GlobalC::ucell.ntype; it++)
+		for(int it=0; it<ucell.ntype; it++)
 		{
 			std::stringstream ss ;
-			ss << PARAM.globalv.global_out_dir << GlobalC::ucell.atoms[it].label << "/v_loc_g.dat" ;
+			ss << PARAM.globalv.global_out_dir << ucell.atoms[it].label << "/v_loc_g.dat" ;
 			std::ofstream ofs_vg( ss.str().c_str() );
 			for(int ig=0;ig<rho_basis->ngg;ig++)
 			{
-				ofs_vg << std::setw(15) << rho_basis->gg_uniq [ig] * GlobalC::ucell.tpiba2 
+				ofs_vg << std::setw(15) << rho_basis->gg_uniq [ig] * ucell.tpiba2 
 				   	<< std::setw(15) << this->vloc(it, ig) << std::endl;
 			}
 			ofs_vg.close();

source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.h

@@ -5,6 +5,7 @@
 #include "module_base/global_variable.h"
 #include "module_base/matrix.h"
 #include "module_basis/module_pw/pw_basis.h"
+#include "module_cell/unitcell.h"
 
 class pseudopot_cell_vl
 {
@@ -19,7 +20,8 @@ class pseudopot_cell_vl
      * @param rho_basis pw basis
      * @return this->vloc 
      */
-    void init_vloc(const ModulePW::PW_Basis* rho_basis);
+    void init_vloc(const UnitCell& ucell,
+                   const ModulePW::PW_Basis* rho_basis);
 
     ModuleBase::matrix vloc;   //(ntype,ngl),the local potential for each atom type(ntype,ngl)
 	bool *numeric; //[ntype], =true
@@ -28,22 +30,28 @@ class pseudopot_cell_vl
 
 	double *zp;   // (npsx),the charge of the pseudopotential
 
-	void allocate(const int ngg);
+	void allocate(const UnitCell& ucell,
+                  const int ngg);
     /**
      * @brief compute the coulomb potential in reciprocal space
      *        v(g) = -\frac{4pi}{V} * zp*e^2 / G^2
      */
-    void vloc_coulomb(const double& zp, double* vloc_1d, const ModulePW::PW_Basis* rho_basis) const;
+    void vloc_coulomb(const UnitCell& ucell,
+                      const double& zp, 
+                      double* vloc_1d, 
+                      const ModulePW::PW_Basis* rho_basis) const;
 	// generate vloc for a particular atom type.
     void vloc_of_g(const int& msh,
                    const double* rab,
                    const double* r,
                    const double* vloc_at,
                    const double& zp,
                    double* vloc,
+                   const UnitCell& ucell,
                    const ModulePW::PW_Basis* rho_basis) const;
 
-    void print_vloc(const ModulePW::PW_Basis* rho_basis) const;
+    void print_vloc(const UnitCell& ucell,
+                    const ModulePW::PW_Basis* rho_basis) const;
 };
 
 #endif // VL_IN_PW