1. refactor class H_XC_pw

ABACUS.develop/source/src_lcao/ELEC_scf.cpp

@@ -189,7 +189,7 @@ void ELEC_scf::scf(const int &istep)
 				// so be careful here, make sure
 				// rho1 and rho2 are the same rho.
 				// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-				pot.v_of_rho(CHR.rho, pot.vr);
+				pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core);
 				en.delta_escf();
 				if (vext == 0)	
 				{
@@ -417,7 +417,7 @@ void ELEC_scf::scf(const int &istep)
 		if(!conv_elec)
 		{
 			// option 1
-			pot.v_of_rho(CHR.rho, pot.vr);
+			pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core);
 			en.delta_escf();
 
 			// option 2
@@ -426,16 +426,16 @@ void ELEC_scf::scf(const int &istep)
 			// use real E_tot functional.
 			//------------------------------
 			/*
-			pot.v_of_rho(CHR.rho_save, pot.vr);
+			pot.vr = pot.v_of_rho(CHR.rho_save, CHR.rho);
 			en.calculate_etot();
 			en.print_etot(conv_elec, istep, iter, dr2, 0.0, ETHR, avg_iter,0);
-			pot.v_of_rho(CHR.rho, pot.vr);
+			pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core);
 			en.delta_escf();
 			*/
 		}
 		else
 		{
-			pot.v_of_rho(CHR.rho, pot.vnew);
+			pot.vnew = pot.v_of_rho(CHR.rho, CHR.rho_core);
 			//(used later for scf correction to the forces )
 			pot.vnew -= pot.vr;
 			en.descf = 0.0;

ABACUS.develop/source/src_pw/H_XC_pw.cpp

@@ -5,21 +5,21 @@
 double H_XC_pw::etxc;
 double H_XC_pw::vtxc;
 
-void H_XC_pw::v_xc
+std::tuple<double,double,matrix> H_XC_pw::v_xc
 (
 	const int &nrxx, // number of real-space grid
 	const int &ncxyz, // total number of charge grid
 	const double &omega, // volume of cell
-    double **rho_in,
-	double *rho_core, // core charge density
-    matrix &v)
+    const double*const*const rho_in,
+	const double*const rho_core) // core charge density
 {
     TITLE("H_XC_pw","v_xc");
     timer::tick("H_XC_pw","v_xc");
 
     //Exchange-Correlation potential Vxc(r) from n(r)
-    etxc = 0.0;
-    vtxc = 0.0;
+    double et_xc = 0.0;
+    double vt_xc = 0.0;
+	matrix v(NSPIN, nrxx);
 
     // the square of the e charge
     // in Rydeberg unit, so * 2.0.
@@ -51,9 +51,9 @@ void H_XC_pw::v_xc
                 XC_Functional::xc(arhox, ex, ec, vx[0], vc[0]);
                 v(0,ir) = e2 * (vx[0] + vc[0]);
 				// consider the total charge density
-                etxc += e2 * (ex + ec) * rhox;
+                et_xc += e2 * (ex + ec) * rhox;
 				// only consider rho_in
-                vtxc += v(0, ir) * rho_in[0][ir];
+                vt_xc += v(0, ir) * rho_in[0][ir];
             } // endif
         } //enddo
     }
@@ -97,9 +97,9 @@ void H_XC_pw::v_xc
                     v(is, ir) = e2 * (vx[is] + vc[is]);
                 }
 
-                etxc += e2 * (ex + ec) * rhox;
+                et_xc += e2 * (ex + ec) * rhox;
 
-                vtxc += v(0, ir) * rho_in[0][ir] + v(1, ir) * rho_in[1][ir];
+                vt_xc += v(0, ir) * rho_in[0][ir] + v(1, ir) * rho_in[1][ir];
             }
         }
 
@@ -129,18 +129,18 @@ void H_XC_pw::v_xc
 
                 XC_Functional::xc_spin( arhox, zeta, ex, ec, vx[0], vx[1], vc[0], vc[1] );
 
-                etxc += e2 * ( ex + ec ) * rhox;
+                et_xc += e2 * ( ex + ec ) * rhox;
 
                 v(0, ir) = e2*( 0.5 * ( vx[0] + vc[0] + vx[1] + vc[1] ) );
-                vtxc += v(0,ir) * rho_in[0][ir];
+                vt_xc += v(0,ir) * rho_in[0][ir];
 
                 double vs = 0.5 * ( vx[0] + vc[0] - vx[1] - vc[1] );
                 if ( amag > vanishing_charge )
                 {
                     for(int ipol = 1;ipol< 4;ipol++)
                     {
                         v(ipol, ir) = e2 * vs * rho_in[ipol][ir] / amag;
-                        vtxc += v(ipol,ir) * rho_in[ipol][ir];
+                        vt_xc += v(ipol,ir) * rho_in[ipol][ir];
                     }//end do
                 }//end if
             }//end if
@@ -150,16 +150,16 @@ void H_XC_pw::v_xc
 
     // add gradient corrections (if any)
     // mohan modify 2009-12-15
-    GGA_PW::gradcorr(etxc, vtxc, v);
+    GGA_PW::gradcorr(et_xc, vt_xc, v);
 
-    // parallel code : collect vtxc,etxc
+    // parallel code : collect vt_xc,et_xc
     // mohan add 2008-06-01
-    Parallel_Reduce::reduce_double_pool( etxc );
-    Parallel_Reduce::reduce_double_pool( vtxc );
+    Parallel_Reduce::reduce_double_pool( et_xc );
+    Parallel_Reduce::reduce_double_pool( vt_xc );
 
-    etxc *= omega / ncxyz;
-    vtxc *= omega / ncxyz;
+    et_xc *= omega / ncxyz;
+    vt_xc *= omega / ncxyz;
 
     timer::tick("H_XC_pw","v_xc");
-    return;
+    return std::make_tuple(et_xc,vt_xc,v);
 }

ABACUS.develop/source/src_pw/H_XC_pw.h

@@ -25,13 +25,12 @@ class H_XC_pw
 	static double vtxc;
 
 	// compute the exchange-correlation energy 
-    static void v_xc(
+    static std::tuple<double,double,matrix> v_xc(
 		const int &nrxx, // number of real-space grid
 		const int &ncxyz, // total number of charge grid
 		const double &omega, // volume of cell
-		double** rho_in, 
-		double *rho_core, // core charge density
-		matrix &v);
+		const double*const*const rho_in, 
+		const double*const rho_core); // core charge density
 
 };
 

ABACUS.develop/source/src_pw/efield.cpp

@@ -13,7 +13,7 @@ double Efield::etotefield; // energy correction due to the field
 double Efield::bvec[3];
 double Efield::bmod;
 
-void Efield::add_efield(double* rho, double* v_in)
+void Efield::add_efield(const double*const rho, double* v_in)
 {
 	TITLE("Efield","add_efield");
 
@@ -207,7 +207,7 @@ void Efield::add_efield(double* rho, double* v_in)
 }
 
 void Efield::compute_el_dip(const double &emaxpos, const double &eopreg,
-	const int &edir, const double *rho, double &e_dipole)const
+	const int &edir, const double*const rho, double &e_dipole)const
 {
 	TITLE("Efield","compute_el_dip");
 

ABACUS.develop/source/src_pw/efield.h

@@ -10,7 +10,7 @@ class Efield
 	Efield();
 	~Efield();
 
-	void add_efield(double* rho, double* v_in);
+	void add_efield(const double*const rho, double* v_in);
 
 	// efield.
     static int edir;         // direction of the field
@@ -26,7 +26,7 @@ class Efield
 	private:
 
 	void compute_el_dip(const double& emaxpos, const double &eopreg, 
-		const int &edir, const double* rho, double& e_dipole)const;
+		const int &edir, const double*const rho, double& e_dipole)const;
 
 	void compute_ion_dip(const double& emaxpos, const double& eoperg, const int &edir, 
 		double& ion_dipole, const double& bmod, const double* bg)const;

ABACUS.develop/source/src_pw/electrons.cpp

@@ -286,7 +286,7 @@ void Electrons::self_consistent(const int &istep)
         if (!conv_elec)
         {
             // not converged yet, calculate new potential from mixed charge density
-            pot.v_of_rho(CHR.rho, pot.vr);
+            pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core);
 
             // because <T+V(ionic)> = <eband+deband> are calculated after sum
             // band, using output charge density.
@@ -308,7 +308,7 @@ void Electrons::self_consistent(const int &istep)
 
             // mohan fix bug 2012-06-05,
             // the new potential V(PL)+V(H)+V(xc)
-            pot.v_of_rho(CHR.rho, pot.vr);
+            pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core);
             //cout<<"Exc = "<<en.etxc<<endl;
             //( vnew used later for scf correction to the forces )
             pot.vnew = pot.vr - pot.vnew;

ABACUS.develop/source/src_pw/forces.cpp

@@ -551,8 +551,10 @@ void Forces::cal_force_ew(matrix& forceion)
 void Forces::cal_force_cc(matrix& forcecc)
 {
 	// recalculate the exchange-correlation potential.
-    matrix vxc(NSPIN, pw.nrxx);
-    H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core, vxc);
+    const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core);
+	H_XC_pw::etxc    = std::get<0>(etxc_vtxc_v);			// may delete?
+	H_XC_pw::vtxc    = std::get<1>(etxc_vtxc_v);			// may delete?
+	const matrix vxc = std::get<2>(etxc_vtxc_v);
 
     complex<double> * psiv = new complex<double> [pw.nrxx];
     ZEROS(psiv, pw.nrxx);

ABACUS.develop/source/src_pw/potential.cpp

@@ -197,7 +197,7 @@ void Potential::init_pot(
     //----------------------------------------------------------
     // (3) compute Hartree and XC potentials saves in vr 
     //----------------------------------------------------------
-    this->v_of_rho(CHR.rho, vr);
+    this->vr = this->v_of_rho(CHR.rho, CHR.rho_core);
 
     //----------------------------------------------------------
     // (4) total potentials 
@@ -276,48 +276,44 @@ void Potential::set_local_pot(
 // The XC potential is computed in real space, while the
 // Hartree potential is computed in reciprocal space.
 //==========================================================
-void Potential::v_of_rho
-(
-    double **rho_in,
-    matrix &v_in
-)
+matrix Potential::v_of_rho(
+	const double*const*const rho_in,
+	const double * const rho_core_in)
 {
     TITLE("Potential","v_of_rho");
-    v_in.zero_out();
-
     timer::tick("Potential","v_of_rho",'E');
 
+    matrix v(NSPIN,pw.nrxx);
+
 //----------------------------------------------------------
 //  calculate the exchange-correlation potential
 //----------------------------------------------------------
 
 	#ifdef USE_LIBXC
-	H_XC_pw::etxc = 0.0;
-	H_XC_pw::vtxc = 0.0;
     const auto etxc_vtxc_v = Potential_Libxc::v_xc(rho_in, CHR.rho_core);
-	H_XC_pw::etxc += std::get<0>(etxc_vtxc_v);
-	H_XC_pw::vtxc += std::get<1>(etxc_vtxc_v);
-	v_in          += std::get<2>(etxc_vtxc_v);
 	#else
-    H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, rho_in, CHR.rho_core, v_in);
+    const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, rho_in, CHR.rho_core);
 	#endif
+	H_XC_pw::etxc = std::get<0>(etxc_vtxc_v);
+	H_XC_pw::vtxc = std::get<1>(etxc_vtxc_v);
+	v            += std::get<2>(etxc_vtxc_v);
 
 //----------------------------------------------------------
 //  calculate the Hartree potential
 //----------------------------------------------------------
-	v_in += H_Hartree_pw::v_hartree(ucell, pw, NSPIN, rho_in);
+	v += H_Hartree_pw::v_hartree(ucell, pw, NSPIN, rho_in);
 
     // mohan add 2011-06-20
     if(EFIELD && DIPOLE)
     {
         Efield EFID;
         for (int is = 0;is < NSPIN;is++)
         {
-            EFID.add_efield(rho_in[is], &v_in.c[is*pw.nrxx]);
+            EFID.add_efield(rho_in[is], &v.c[is*pw.nrxx]);
         }
     }
     timer::tick("Potential","v_of_rho",'E');
-    return;
+    return v;
 } //end subroutine v_of_rho
 
 

ABACUS.develop/source/src_pw/potential.h

@@ -41,7 +41,9 @@ class Potential
 		ComplexMatrix &sf // structure factors
 	);
 
-    void v_of_rho(double** rho_in, matrix &v);
+    matrix v_of_rho(
+		const double*const*const rho_in,
+		const double * const rho_core_in);
 
     void set_vr_eff(void);
 

ABACUS.develop/source/src_pw/sto_elec.cpp

@@ -215,7 +215,7 @@ void Stochastic_Elec::scf_stochastic(const int &istep)
         if (!conv_elec)
         {
 			// not converged yet, calculate new potential from mixed charge density
-            pot.v_of_rho(CHR.rho, pot.vr);
+            pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core);
 
 			// because <T+V(ionic)> = <eband+deband> are calculated after sum
 			// band, using output charge density.
@@ -235,7 +235,7 @@ void Stochastic_Elec::scf_stochastic(const int &istep)
 			}
 
 			// the new potential V(PL)+V(H)+V(xc)
-            pot.v_of_rho(CHR.rho, pot.vr);
+            pot.vr = pot.v_of_rho(CHR.rho, CHR.rho_core);
 
             //( vnew used later for scf correction to the forces )
 	    	pot.vnew = pot.vr - pot.vnew;

ABACUS.develop/source/src_pw/stress_func_cc.cpp

@@ -34,8 +34,10 @@ void Stress_Func::stress_cc(matrix& sigma, const bool is_pw)
 	}
 
 	//recalculate the exchange-correlation potential
-    matrix vxc(NSPIN, pw.nrxx);
-    H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core, vxc);
+    const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, CHR.rho, CHR.rho_core);
+	H_XC_pw::etxc    = std::get<0>(etxc_vtxc_v);			// may delete?
+	H_XC_pw::vtxc    = std::get<1>(etxc_vtxc_v);			// may delete?
+	const matrix vxc = std::get<2>(etxc_vtxc_v);
 
 	complex<double> * psic = new complex<double> [pw.nrxx];