Merge branch 'develop' into ucell5

Author: A-006

source/module_elecstate/elecstate_pw.h

@@ -17,6 +17,7 @@ class ElecStatePW : public ElecState
 {
   private:
     using Real = typename GetTypeReal<T>::type;
+
   public:
     ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
                 Charge* chg_in,
@@ -26,46 +27,59 @@ class ElecStatePW : public ElecState
                 ModulePW::PW_Basis* rhodpw_in,
                 ModulePW::PW_Basis* rhopw_in,
                 ModulePW::PW_Basis_Big* bigpw_in);
-    // void init(Charge* chg_in):charge(chg_in){} override;
 
     ~ElecStatePW();
-    // interface for HSolver to calculate rho from Psi
+
+    //! interface for HSolver to calculate rho from Psi
     virtual void psiToRho(const psi::Psi<T, Device>& psi);
-    // return current electronic density rho, as a input for constructing Hamiltonian
-    // const double* getRho(int spin) const override;
-    virtual void cal_tau(const psi::Psi<T, Device>& psi);
 
-    // update charge density for next scf step
-    // void getNewRho() override;
+    virtual void cal_tau(const psi::Psi<T, Device>& psi);
 
     Real* becsum = nullptr;
-    // init rho_data and kin_r_data
+
+    //! init rho_data and kin_r_data
     void init_rho_data();
-    Real ** rho = nullptr, ** kin_r = nullptr; //[Device] [spin][nrxx] rho and kin_r
+    Real** rho = nullptr; 
+    Real** kin_r = nullptr; //[Device] [spin][nrxx] rho and kin_r
 
   protected:
+
     ModulePW::PW_Basis* rhopw_smooth = nullptr;
+
     ModulePW::PW_Basis_K* basis = nullptr;
+
     UnitCell* ucell = nullptr;
+
     const pseudopot_cell_vnl* ppcell = nullptr;
 
-    // calculate electronic charge density on grid points or density matrix in real space
-    // the consequence charge density rho saved into rho_out, preparing for charge mixing.
+    //! calculate electronic charge density on grid points or density matrix in real space
+    //! the consequence charge density rho saved into rho_out, preparing for charge mixing.
     void updateRhoK(const psi::Psi<T, Device>& psi); // override;
-    // sum over all pools for rho and ebands
+    
+    //! sum over all pools for rho and ebands
     void parallelK();
-    // calcualte rho for each k
+    
+    //! calcualte rho for each k
     void rhoBandK(const psi::Psi<T, Device>& psi);
-    // add to the charge density in reciprocal space the part which is due to the US augmentation.
+    
+    //! add to the charge density in reciprocal space the part which is due to the US augmentation.
     void add_usrho(const psi::Psi<T, Device>& psi);
-    // \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
+
+    //! Non-local pseudopotentials
+    //! \sum_lm Q_lm(r) \sum_i <psi_i|beta_l><beta_m|psi_i> w_i
     void addusdens_g(const Real* becsum, T* rhog);
 
     Device * ctx = {};
+
     bool init_rho = false;
+
     mutable T* vkb = nullptr;
-    Real * rho_data = nullptr, * kin_r_data = nullptr;
-    T * wfcr = nullptr, * wfcr_another_spin = nullptr;
+
+    Real* rho_data = nullptr; 
+    Real* kin_r_data = nullptr;
+    T* wfcr = nullptr; 
+    T* wfcr_another_spin = nullptr;
+
   private:
     using meta_op = hamilt::meta_pw_op<Real, Device>;
     using elecstate_pw_op = elecstate::elecstate_pw_op<Real, Device>;
@@ -85,4 +99,4 @@ class ElecStatePW : public ElecState
 
 } // namespace elecstate
 
-#endif
+#endif

source/module_elecstate/test/charge_extra_test.cpp

@@ -87,7 +87,7 @@ Structure_Factor::Structure_Factor()
 Structure_Factor::~Structure_Factor()
 {
 }
-void Structure_Factor::setup_structure_factor(UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
+void Structure_Factor::setup_structure_factor(const UnitCell* Ucell, const ModulePW::PW_Basis* rho_basis)
 {
 }
 

source/module_hamilt_general/module_surchem/cal_pseudo.cpp

@@ -8,7 +8,7 @@ void surchem::gauss_charge(const UnitCell& cell,
                            complex<double>* N,
                            Structure_Factor* sf)
 {
-    sf->setup_structure_factor(&GlobalC::ucell, rho_basis); // call strucFac(ntype,ngmc)
+    sf->setup_structure_factor(&cell, rho_basis); // call strucFac(ntype,ngmc)
     for (int it = 0; it < cell.ntype; it++)
     {
         double RCS = GetAtom.atom_RCS[cell.atoms[it].ncpp.psd];

source/module_hamilt_general/module_surchem/cal_vcav.cpp

@@ -2,7 +2,10 @@
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #include "surchem.h"
 
-void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW_Basis* rho_basis)
+void lapl_rho(const double& tpiba2,
+              const std::complex<double>* rhog, 
+              double* lapn, 
+              const ModulePW::PW_Basis* rho_basis)
 {
     std::complex<double> *gdrtmpg = new std::complex<double>[rho_basis->npw];
     ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);
@@ -22,9 +25,9 @@ void lapl_rho(const std::complex<double>* rhog, double* lapn, const ModulePW::PW
         // bring the gdr from G --> R
         rho_basis->recip2real(aux, aux);
         // remember to multily 2pi/a0, which belongs to G vectors.
-        for (int ir = 0; ir < rho_basis->nrxx; ir++) {
-            lapn[ir] -= aux[ir].real() * GlobalC::ucell.tpiba2;
-}
+        for (int ir = 0; ir < rho_basis->nrxx; ir++)
+            lapn[ir] -= aux[ir].real() * tpiba2;
+
     }
 
     delete[] gdrtmpg;
@@ -70,7 +73,7 @@ void surchem::createcavity(const UnitCell& ucell,
     ModuleBase::GlobalFunc::ZEROS(lapn, rho_basis->nrxx);
 
     // nabla n
-    XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, GlobalC::ucell.tpiba);
+    XC_Functional::grad_rho(PS_TOTN, nablan, rho_basis, ucell.tpiba);
 
     //  |\nabla n |^2 = nablan_2
     for (int ir = 0; ir < rho_basis->nrxx; ir++)
@@ -79,7 +82,7 @@ void surchem::createcavity(const UnitCell& ucell,
     }
 
     // Laplacian of n
-    lapl_rho(PS_TOTN, lapn, rho_basis);
+    lapl_rho(ucell.tpiba2,PS_TOTN, lapn, rho_basis);
 
     //-------------------------------------------------------------
     // add -Lap(n)/|\nabla n| to vwork and copy \sqrt(|\nabla n|^2)
@@ -131,7 +134,7 @@ void surchem::createcavity(const UnitCell& ucell,
 
     // \nabla(1 / |\nabla n|), ggn in real space
     ModuleBase::Vector3<double> *ggn = new ModuleBase::Vector3<double>[rho_basis->nrxx];
-    XC_Functional::grad_rho(inv_gn, ggn, rho_basis, GlobalC::ucell.tpiba);
+    XC_Functional::grad_rho(inv_gn, ggn, rho_basis, ucell.tpiba);
 
     //-------------------------------------------------------------
     // add -(\nabla n . \nabla(1/ |\nabla n|)) to Vcav in real space

source/module_hamilt_general/module_surchem/cal_vel.cpp

@@ -17,6 +17,7 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis
 }
 
 void eps_pot(const double* PS_TOTN_real,
+             const double& tpiba,
              const complex<double>* phi,
              const ModulePW::PW_Basis* rho_basis,
              double* d_eps,
@@ -36,7 +37,7 @@ void eps_pot(const double* PS_TOTN_real,
     double *phisq = new double[rho_basis->nrxx];
 
     // nabla phi
-    XC_Functional::grad_rho(phi, nabla_phi, rho_basis, GlobalC::ucell.tpiba);
+    XC_Functional::grad_rho(phi, nabla_phi, rho_basis, tpiba);
 
     for (int ir = 0; ir < rho_basis->nrxx; ir++)
     {
@@ -118,7 +119,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell& cell,
     this->Ael *= cell.omega / rho_basis->nxyz;
 
     // the 2nd item of tmp_Vel
-    eps_pot(PS_TOTN_real, Sol_phi, rho_basis, epsilon, epspot);
+    eps_pot(PS_TOTN_real, cell.tpiba, Sol_phi, rho_basis, epsilon, epspot);
 
     for (int i = 0; i < rho_basis->nrxx; i++)
     {

source/module_hamilt_general/module_surchem/sol_force.cpp

@@ -32,7 +32,7 @@ void force_cor_one(const UnitCell& cell,
                 vloc_at[ig] = vloc(it, rho_basis->ig2igg[ig]) * phase;
                 if(rho_basis->ig_gge0 == ig)
                 {
-                    N[ig] = GlobalC::ucell.atoms[it].ncpp.zv / GlobalC::ucell.omega;
+                    N[ig] = cell.atoms[it].ncpp.zv / cell.omega;
                 }
                 else
                 {
@@ -48,9 +48,9 @@ void force_cor_one(const UnitCell& cell,
                 forcesol(iat, 2) += rho_basis->gcar[ig][2] * imag(conj(delta_phi_g[ig]) * N[ig]);
             }
 
-                forcesol(iat, 0) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
-                forcesol(iat, 1) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
-                forcesol(iat, 2) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
+                forcesol(iat, 0) *= (cell.tpiba * cell.omega);
+                forcesol(iat, 1) *= (cell.tpiba * cell.omega);
+                forcesol(iat, 2) *= (cell.tpiba * cell.omega);
             //unit Ry/Bohr
                 forcesol(iat, 0) *= 2 ;
                 forcesol(iat, 1) *= 2 ;
@@ -128,9 +128,9 @@ void force_cor_two(const UnitCell& cell, const ModulePW::PW_Basis* rho_basis, Mo
                 forcesol(iat, 2) -= rho_basis->gcar[ig][2] * imag(conj(Vcav_g[ig]+Vel_g[ig]) * n_pseudo[ig]);
             }
 
-                forcesol(iat, 0) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
-                forcesol(iat, 1) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
-                forcesol(iat, 2) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
+                forcesol(iat, 0) *= (cell.tpiba * cell.omega);
+                forcesol(iat, 1) *= (cell.tpiba * cell.omega);
+                forcesol(iat, 2) *= (cell.tpiba * cell.omega);
             //eV/Ang
                 forcesol(iat, 0) *= 2 ;
                 forcesol(iat, 1) *= 2 ;
@@ -157,17 +157,17 @@ void surchem::cal_force_sol(const UnitCell& cell,
     ModuleBase::TITLE("surchem", "cal_force_sol");
     ModuleBase::timer::tick("surchem", "cal_force_sol");
 
-    int nat = GlobalC::ucell.nat;
+    int nat = cell.nat;
 	ModuleBase::matrix force1(nat, 3);
     ModuleBase::matrix force2(nat, 3);
 
     force_cor_one(cell, rho_basis, vloc, force1);
     force_cor_two(cell, rho_basis,force2);
 
     int iat = 0;
-    for (int it = 0;it < GlobalC::ucell.ntype;it++)
+    for (int it = 0;it < cell.ntype;it++)
 	{
-		for (int ia = 0;ia < GlobalC::ucell.atoms[it].na;ia++)
+		for (int ia = 0;ia < cell.atoms[it].na;ia++)
 		{
             for(int ipol = 0; ipol < 3; ipol++)
             {

source/module_hamilt_general/module_surchem/test/cal_pseudo_test.cpp

@@ -26,6 +26,7 @@ class cal_pseudo_test : public testing::Test
 {
   protected:
     surchem solvent_model;
+    UnitCell ucell;
 };
 
 TEST_F(cal_pseudo_test, gauss_charge)
@@ -41,7 +42,7 @@ TEST_F(cal_pseudo_test, gauss_charge)
     double wfcecut;
     bool gamma_only;
 
-    Setcell::setupcell(GlobalC::ucell);
+    Setcell::setupcell(ucell);
 
     wfcecut = 80;
     gamma_only = false;
@@ -60,7 +61,7 @@ TEST_F(cal_pseudo_test, gauss_charge)
 #endif
     // pwtest.initgrids(lat0,latvec,wfcecut);
 
-    GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
+    GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);
 
     GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
     GlobalC::rhopw->setuptransform();
@@ -75,7 +76,7 @@ TEST_F(cal_pseudo_test, gauss_charge)
     Structure_Factor sf;
     sf.nbspline = -1;
 
-    solvent_model.gauss_charge(GlobalC::ucell, GlobalC::rhopw, N, &sf);
+    solvent_model.gauss_charge(ucell, GlobalC::rhopw, N, &sf);
 
     EXPECT_NEAR(N[14].real(), 0.002, 1e-9);
     EXPECT_NEAR(N[16].real(), -0.001573534, 1e-9);
@@ -88,7 +89,7 @@ TEST_F(cal_pseudo_test, cal_pseudo)
     std::string precision_flag, device_flag;
     precision_flag = "double";
     device_flag = "cpu";
-
+    Setcell::setupcell(ucell);
     ModulePW::PW_Basis pwtest(device_flag, precision_flag);
     GlobalC::rhopw = &pwtest;
     ModuleBase::Matrix3 latvec;
@@ -115,7 +116,7 @@ TEST_F(cal_pseudo_test, cal_pseudo)
 #endif
     // pwtest.initgrids(lat0,latvec,wfcecut);
 
-    GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
+    GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);
     GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
     GlobalC::rhopw->setuptransform();
     GlobalC::rhopw->collect_local_pw();
@@ -135,7 +136,7 @@ TEST_F(cal_pseudo_test, cal_pseudo)
     }
 
     complex<double>* PS_TOTN = new complex<double>[npw];
-    solvent_model.cal_pseudo(GlobalC::ucell, GlobalC::rhopw, Porter_g, PS_TOTN, &sf);
+    solvent_model.cal_pseudo(ucell, GlobalC::rhopw, Porter_g, PS_TOTN, &sf);
 
     EXPECT_NEAR(PS_TOTN[16].real(), 0.098426466, 1e-9);
     EXPECT_NEAR(PS_TOTN[14].real(), 0.102, 1e-9);

source/module_hamilt_general/module_surchem/test/cal_totn_test.cpp

@@ -25,6 +25,7 @@ class cal_totn_test : public testing::Test
 {
   protected:
     surchem solvent_model;
+    UnitCell ucell;
 };
 
 TEST_F(cal_totn_test, cal_totn)
@@ -40,7 +41,7 @@ TEST_F(cal_totn_test, cal_totn)
     double wfcecut;
     bool gamma_only;
 
-    Setcell::setupcell(GlobalC::ucell);
+    Setcell::setupcell(ucell);
 
     wfcecut = 80;
     gamma_only = false;
@@ -59,7 +60,7 @@ TEST_F(cal_totn_test, cal_totn)
 #endif
     // pwtest.initgrids(lat0,latvec,wfcecut);
 
-    GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
+    GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);
 
     GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
     GlobalC::rhopw->setuptransform();
@@ -87,7 +88,7 @@ TEST_F(cal_totn_test, cal_totn)
         vloc[i] = 0.1;
     }
 
-    solvent_model.cal_totn(GlobalC::ucell, GlobalC::rhopw, Porter_g, N, TOTN, vloc);
+    solvent_model.cal_totn(ucell, GlobalC::rhopw, Porter_g, N, TOTN, vloc);
 
     EXPECT_NEAR(TOTN[0].real(), -0.0999496256, 1e-10);
     EXPECT_NEAR(TOTN[0].imag(), -1.299621928166352e-7, 1e-10);
@@ -103,7 +104,8 @@ TEST_F(cal_totn_test, induced_charge)
     std::string precision_flag, device_flag;
     precision_flag = "double";
     device_flag = "cpu";
-
+    Setcell::setupcell(ucell);
+    
     ModulePW::PW_Basis pwtest(device_flag, precision_flag);
     GlobalC::rhopw = &pwtest;
     ModuleBase::Matrix3 latvec;
@@ -128,15 +130,15 @@ TEST_F(cal_totn_test, induced_charge)
 #endif
     // pwtest.initgrids(lat0,latvec,wfcecut);
 
-    GlobalC::rhopw->initgrids(GlobalC::ucell.lat0, GlobalC::ucell.latvec, wfcecut);
+    GlobalC::rhopw->initgrids(ucell.lat0, ucell.latvec, wfcecut);
 
     GlobalC::rhopw->initparameters(gamma_only, wfcecut, distribution_type, xprime);
     GlobalC::rhopw->setuptransform();
     GlobalC::rhopw->collect_local_pw();
     GlobalC::rhopw->collect_uniqgg();
 
     double fac;
-    fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (GlobalC::ucell.tpiba2 * GlobalC::rhopw->gg[0]);
+    fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (ucell.tpiba2 * GlobalC::rhopw->gg[0]);
     complex<double> delta_phi{-2.0347933860e-05, 4.5900395826e-07};
     complex<double> induced_charge;
     induced_charge = -delta_phi / fac;