diff --git a/ABACUS.develop/source/src_global/global_function.h b/ABACUS.develop/source/src_global/global_function.h
index 5117468cb5..f54048f236 100644
--- a/ABACUS.develop/source/src_global/global_function.h
+++ b/ABACUS.develop/source/src_global/global_function.h
@@ -257,4 +257,63 @@ inline const void* MAP_EXIST( const T_map &ms, const T_key1 &key1, const T_key_t
 //==========================================================
 size_t MemAvailable();
 
+
+//==========================================================
+// GLOBAL FUNCTION :
+// NAME : DELETE_MUL_PTR
+// delete Multi-dimensional array pointer
+// example:
+//		int*** v;
+//		DELETE_MUL_PTR(v,N1,N2);
+//	->	for(int i1=0; i1<N1; ++i1){
+//			for(int i2=0; i2<N2; ++i2){
+//				delete[] v[i1][i2];	v[i1][i2]=nullptr;	}
+//			delete[] v[i1];	v[i1]=nullptr;	}
+//		delete[] v;	v=nullptr;
+// Peize Lin add 2021-05-09
+//==========================================================
+template <typename T_element>
+static inline void DELETE_MUL_PTR(T_element* v)
+{
+	delete[] v;
+	v = nullptr;
+}
+template <typename T_element, typename T_N_first, typename... T_N_tail>
+static inline void DELETE_MUL_PTR(T_element* v, const T_N_first N_first, const T_N_tail... N_tail)
+{
+	for(T_N_first i=0; i<N_first; ++i)
+		DELETE_MUL_PTR(v[i],N_tail...);
+	delete[] v;
+	v = nullptr;
+}
+
+//==========================================================
+// GLOBAL FUNCTION :
+// NAME : FREE_MUL_PTR
+// delete Multi-dimensional array pointer
+// example:
+//		int*** v;
+//		DELETE_MUL_PTR(v,N1,N2);
+//	->	for(int i1=0; i1<N1; ++i1){
+//			for(int i2=0; i2<N2; ++i2){
+//				free(v[i1][i2]);	v[i1][i2]=nullptr;	}
+//			free(v[i1]);	v[i1]=nullptr;	}
+//		free(v);	v=nullptr;
+// Peize Lin add 2021-05-09
+//==========================================================
+template <typename T_element>
+static inline void FREE_MUL_PTR(T_element* v)
+{
+	free(v);
+	v = nullptr;
+}
+template <typename T_element, typename T_N_first, typename... T_N_tail>
+static inline void FREE_MUL_PTR(T_element* v, const T_N_first N_first, const T_N_tail... N_tail)
+{
+	for(T_N_first i=0; i<N_first; ++i)
+		FREE_MUL_PTR(v[i],N_tail...);
+	free(v);
+	v = nullptr;
+}
+
 #endif
diff --git a/ABACUS.develop/source/src_global/ylm.cpp b/ABACUS.develop/source/src_global/ylm.cpp
index 0626954319..504b77af2e 100644
--- a/ABACUS.develop/source/src_global/ylm.cpp
+++ b/ABACUS.develop/source/src_global/ylm.cpp
@@ -6,7 +6,7 @@
 
 using namespace std;
 int Ylm::nlm = 0;
-vector<double> Ylm::ylmcoef(100);
+std::vector<double> Ylm::ylmcoef(100);
 
 // here Lmax == max angular momentum + 1
 void Ylm::get_ylm_real( const int &Lmax, const Vector3<double> &vec, double ylmr[] )
@@ -486,9 +486,11 @@ void Ylm::sph_harm
  	const double& xdr,
 	const double& ydr,
 	const double& zdr,
-	double rly[]
+	std::vector<double> &rly
 )
 {
+	rly.resize( (Lmax+1)*(Lmax+1) );
+
 	//begin calculation
 	/***************************
 			 L = 0
@@ -623,7 +625,7 @@ void Ylm::rl_sph_harm
  	const double& x,
 	const double& y,
 	const double& z,
-	vector<double>& rly
+	std::vector<double>& rly
 )
 {
 	rly.resize( (Lmax+1)*(Lmax+1) );
@@ -762,12 +764,12 @@ void Ylm::grad_rl_sph_harm
  	const double& x,
 	const double& y,
 	const double& z,
-	vector<double>& rly,
-	vector<vector<double>>& grly
+	std::vector<double>& rly,
+	std::vector<std::vector<double>>& grly
 )
 {
 	rly.resize( (Lmax+1)*(Lmax+1) );
-	grly.resize( (Lmax+1)*(Lmax+1), vector<double>(3) );
+	grly.resize( (Lmax+1)*(Lmax+1), std::vector<double>(3) );
 
 	double radius2 = x*x+y*y+z*z;
 	double tx = 2.0*x;
@@ -1105,7 +1107,7 @@ void Ylm::test1 (void)
 	int nu = 100;
 	
 	// Peize Lin change rlya 2016-08-26
-	vector<double> rlya;
+	std::vector<double> rlya;
 	double rlyb[400];
 	ZEROS( rlyb, 400);
 	
@@ -1136,10 +1138,10 @@ void Ylm::test2 (void)
 	
 	//int nu = 100;
 
-	vector<double> rlya;
+	std::vector<double> rlya;
 	double rlyb[400];
 	
-	vector<vector<double>> grlya;
+	std::vector<std::vector<double>> grlya;
 	double grlyb[400][3];
 	
 	Ylm::grad_rl_sph_harm (9, R.x, R.y, R.z, rlya, grlya);
diff --git a/ABACUS.develop/source/src_global/ylm.h b/ABACUS.develop/source/src_global/ylm.h
index f0c1df9109..9014397b80 100644
--- a/ABACUS.develop/source/src_global/ylm.h
+++ b/ABACUS.develop/source/src_global/ylm.h
@@ -48,7 +48,7 @@ class Ylm
 			const double& xdr,
 			const double& ydr,
 			const double& zdr,
-			double rly[]);
+			std::vector<double> &rly);
 	
 	// (6) used in getting overlap.
 		// Peize Lin change rly 2016-08-26
@@ -57,7 +57,7 @@ class Ylm
 			const double& x,
 			const double& y,
 			const double& z,
-			vector<double>& rly);
+			std::vector<double>& rly);
 	
 	// (6) used in getting derivative of overlap.
 		// Peize Lin change rly, grly 2016-08-26
@@ -66,8 +66,8 @@ class Ylm
 			const double& x,
 			const double& y,
 			const double& z,
-			vector<double>& rly,
-			vector<vector<double>>& grly);
+			std::vector<double>& rly,
+			std::vector<std::vector<double>>& grly);
 			
 	static void set_coefficients ();
 	static std::vector<double> ylmcoef;
diff --git a/ABACUS.develop/source/src_lcao/gint_gamma.h b/ABACUS.develop/source/src_lcao/gint_gamma.h
index e9633b9971..12eaf6fadf 100644
--- a/ABACUS.develop/source/src_lcao/gint_gamma.h
+++ b/ABACUS.develop/source/src_lcao/gint_gamma.h
@@ -96,7 +96,7 @@ class Gint_Gamma : public Grid_Base_Beta
 	void setVindex(const int ncyz, const int ibx, const int jby, const int kbz, int* vindex) const;
 
 	void cal_psir_ylm_rho(int size, int grid_index, double delta_r,
-        double** distance, double* ylma,
+        double** distance,
         int* at, int* block_index, int* block_iw, int* block_size, 
         int** cal_flag, double** psir_ylm);
 
diff --git a/ABACUS.develop/source/src_lcao/gint_gamma_env.cpp b/ABACUS.develop/source/src_lcao/gint_gamma_env.cpp
index b150aaf5dc..4453556487 100644
--- a/ABACUS.develop/source/src_lcao/gint_gamma_env.cpp
+++ b/ABACUS.develop/source/src_lcao/gint_gamma_env.cpp
@@ -63,9 +63,6 @@ void Gint_Gamma::gamma_envelope(const double* wfc, double* rho)
 		}
 	}
 
-	double* ylma = new double[nnnmax]; // Ylm for each atom: [bxyz, nnnmax]
-	ZEROS(ylma, nnnmax);
-
 	double mt[3]={0,0,0};
 	double *vldr3 = new double[pw.bxyz];
 	double v1 = 0.0;
@@ -136,6 +133,7 @@ void Gint_Gamma::gamma_envelope(const double* wfc, double* rho)
 							continue;
 						}
 						
+						std::vector<double> ylma;
 						//if(distance[id] > GridT.orbital_rmax) continue;
 						//	Ylm::get_ylm_real(this->nnn[it], this->dr[id], ylma);
 						if (distance[ib][id] < 1.0E-9) distance[ib][id] += 1.0E-9;
@@ -244,7 +242,6 @@ void Gint_Gamma::gamma_envelope(const double* wfc, double* rho)
     }// i
 
 	delete[] vindex;
-	delete[] ylma;
 	delete[] vldr3;
 	
 	if(max_size!=0) 
diff --git a/ABACUS.develop/source/src_lcao/gint_gamma_mull.cpp b/ABACUS.develop/source/src_lcao/gint_gamma_mull.cpp
index f19990fa4e..bc595e860f 100644
--- a/ABACUS.develop/source/src_lcao/gint_gamma_mull.cpp
+++ b/ABACUS.develop/source/src_lcao/gint_gamma_mull.cpp
@@ -68,9 +68,6 @@ void Gint_Gamma::gamma_mulliken(double** mulliken)
 		}
 	}
 
-	double* ylma = new double[nnnmax]; // Ylm for each atom: [bxyz, nnnmax]
-	ZEROS(ylma, nnnmax);
-
 	double mt[3]={0,0,0};
 	double *vldr3 = new double[pw.bxyz];
 	double v1 = 0.0;
@@ -141,6 +138,7 @@ void Gint_Gamma::gamma_mulliken(double** mulliken)
 							continue;
 						}
 						
+						std::vector<double> ylma;
 						//if(distance[id] > GridT.orbital_rmax) continue;
 						//	Ylm::get_ylm_real(this->nnn[it], this->dr[id], ylma);
 						if (distance[ib][id] < 1.0E-9) distance[ib][id] += 1.0E-9;
@@ -291,7 +289,6 @@ void Gint_Gamma::gamma_mulliken(double** mulliken)
         }// j
     }// i
 
-	delete[] ylma;
 	delete[] vldr3;
 	
 	if(max_size!=0) 
diff --git a/ABACUS.develop/source/src_lcao/gint_gamma_rho.cpp b/ABACUS.develop/source/src_lcao/gint_gamma_rho.cpp
index 02337a4b7b..56782a60ff 100644
--- a/ABACUS.develop/source/src_lcao/gint_gamma_rho.cpp
+++ b/ABACUS.develop/source/src_lcao/gint_gamma_rho.cpp
@@ -31,8 +31,7 @@ void Gint_Gamma::cal_psir_ylm_rho(
 	int size, // number of atoms on this grid 
 	int grid_index,  
 	double delta_r, // delta_r of uniform grid
-	double** distance, 
-	double* ylma, 
+	double** distance,
 	int* at, 
 	int* block_index, 
 	int* block_iw, 
@@ -90,6 +89,7 @@ void Gint_Gamma::cal_psir_ylm_rho(
 
 			cal_flag[ib][id]=1;
 			
+			std::vector<double> ylma;
 			//if(distance[id] > GridT.orbital_rmax) continue;
 			//	Ylm::get_ylm_real(this->nnn[it], this->dr[id], ylma);
 			if (distance[ib][id] < 1.0E-9) distance[ib][id] += 1.0E-9;
@@ -363,10 +363,7 @@ double Gint_Gamma::gamma_charge(void)					// Peize Lin update OpenMP 2020.09.28
 			}
 
 			int *block_iw=new int[max_size]; // index of wave functions of each block;
-		
-			double* ylma = new double[nnnmax]; // Ylm for each atom: [bxyz, nnnmax]
-			ZEROS(ylma, nnnmax);
-		
+				
 			double v1 = 0.0;
 			int* vindex=new int[pw.bxyz];
 			ZEROS(vindex, pw.bxyz);
@@ -399,7 +396,7 @@ double Gint_Gamma::gamma_charge(void)					// Peize Lin update OpenMP 2020.09.28
 						if(size==0) continue;
 						this->setVindex(ncyz, ibx, jby, kbz, vindex);
 						
-						this->cal_psir_ylm_rho(size, grid_index_thread, delta_r, distance, ylma,
+						this->cal_psir_ylm_rho(size, grid_index_thread, delta_r, distance,
 								at, block_index, block_iw, block_size, 
 								cal_flag, psir_ylm);
 						
@@ -411,7 +408,6 @@ double Gint_Gamma::gamma_charge(void)					// Peize Lin update OpenMP 2020.09.28
 			}// i
 			
 			delete[] vindex;
-			delete[] ylma;
 			
 			delete[] block_iw;
 			//OUT(ofs_running, "block_iw deleted");
diff --git a/ABACUS.develop/source/src_lcao/gint_gamma_vl.cpp b/ABACUS.develop/source/src_lcao/gint_gamma_vl.cpp
index 70da46fb66..9c45526711 100644
--- a/ABACUS.develop/source/src_lcao/gint_gamma_vl.cpp
+++ b/ABACUS.develop/source/src_lcao/gint_gamma_vl.cpp
@@ -54,7 +54,6 @@ inline void cal_psir_ylm(
 	double*** dr, // dr[ bxyz ; atoms_on_this_big_cell; xyz ] 
 	double** distance, // [ bxyz ; atoms_on_this_big_cell]
 	const Numerical_Orbital_Lm* pointer, // pointer for ORB.Phi[it].PhiLN 
-	double* ylma, // spherical harmonic functions 
 	int* colidx,  // count total number of atomis orbitals 
 	int* block_iw, // seems not belong to this subroutine
 	int* bsize,  // ??
@@ -118,6 +117,7 @@ inline void cal_psir_ylm(
             
 			cal_flag[ib][id]=1;
 
+			std::vector<double> ylma;
             //if(distance[id] > GridT.orbital_rmax) continue;
             //    Ylm::get_ylm_real(this->nnn[it], this->dr[id], ylma);
             if (distance[ib][id] < 1.0E-9) distance[ib][id] += 1.0E-9;
@@ -511,8 +511,6 @@ void Gint_Gamma::gamma_vlocal(void)						// Peize Lin update OpenMP 2020.09.27
 			//------------------------------------------------------
 			// spherical harmonic functions Ylm
 			//------------------------------------------------------
-			double* ylma=new double[nnnmax];
-			ZEROS(ylma, nnnmax);
 			double *vldr3=new double[pw.bxyz];
 			ZEROS(vldr3, pw.bxyz);
 			int* vindex=new int[pw.bxyz];
@@ -628,7 +626,7 @@ void Gint_Gamma::gamma_vlocal(void)						// Peize Lin update OpenMP 2020.09.27
 						// compute atomic basis phi(r) with both radial and angular parts
 						//------------------------------------------------------------------
 						cal_psir_ylm(size, grid_index_thread, delta_r, phi, mt, dr, 
-						distance, pointer, ylma, colidx, block_iw, bsize,  psir_ylm, cal_flag);
+						distance, pointer, colidx, block_iw, bsize,  psir_ylm, cal_flag);
 
 						//------------------------------------------------------------------
 						// calculate <phi_i|V|phi_j>
@@ -675,7 +673,6 @@ void Gint_Gamma::gamma_vlocal(void)						// Peize Lin update OpenMP 2020.09.27
 			}
 			delete[] cal_flag;
 			delete[] vindex;
-			delete[] ylma;
 			delete[] vldr3;    
 			delete[] block_iw;
 			delete[] psir_vlbr3;
diff --git a/ABACUS.develop/source/src_lcao/gint_k.cpp b/ABACUS.develop/source/src_lcao/gint_k.cpp
index 149d52f109..b557cf62a9 100644
--- a/ABACUS.develop/source/src_lcao/gint_k.cpp
+++ b/ABACUS.develop/source/src_lcao/gint_k.cpp
@@ -1150,8 +1150,7 @@ void Gint_k::set_ijk_atom(
 	double*** psir_ylm, 
 	double*** dr, 
 	bool** cal_flag, 
-	double** distance, 
-	double* ylma, 
+	double** distance,
 	const double &delta_r)
 {
 	const Numerical_Orbital_Lm* pointer;
@@ -1201,7 +1200,7 @@ void Gint_k::set_ijk_atom(
 
 			if (distance[ib][id] < 1.0E-9) distance[ib][id] += 1.0E-9;
 
-			
+			std::vector<double> ylma;
 			Ylm::sph_harm ( ucell.atoms[it].nwl,
 					dr[ib][id][0] / distance[ib][id],
 					dr[ib][id][1] / distance[ib][id],
diff --git a/ABACUS.develop/source/src_lcao/gint_k.h b/ABACUS.develop/source/src_lcao/gint_k.h
index 444819f459..f9a369c947 100644
--- a/ABACUS.develop/source/src_lcao/gint_k.h
+++ b/ABACUS.develop/source/src_lcao/gint_k.h
@@ -85,7 +85,7 @@ class Gint_k : public Gint_k_init
 	// set the orbital/Ylm information on each real space grid.
 	void set_ijk_atom(const int &grid_index, const int &size,
 		double*** psir_ylm, double*** dr, bool** cal_flag, 
-		double** distance, double* ylma, const double &delta_r);
+		double** distance, const double &delta_r);
 
 	//------------------------------------------------------
 	// in gint_k_vl.cpp 
diff --git a/ABACUS.develop/source/src_lcao/gint_k_rho.cpp b/ABACUS.develop/source/src_lcao/gint_k_rho.cpp
index 8214a69b4c..8f31d7f1bf 100644
--- a/ABACUS.develop/source/src_lcao/gint_k_rho.cpp
+++ b/ABACUS.develop/source/src_lcao/gint_k_rho.cpp
@@ -28,7 +28,7 @@ inline void setVindex(const int ncyz, const int ibx, const int jby, const int kb
 }
 
 
-inline void cal_psir_ylm(int size, int grid_index, double delta_r, double* ylma,
+inline void cal_psir_ylm(int size, int grid_index, double delta_r,
 						int* at, int* uc, int* block_index, int* block_iw, int* block_size, 
 						bool** cal_flag, double** psir_ylm)
 {
@@ -84,6 +84,7 @@ inline void cal_psir_ylm(int size, int grid_index, double delta_r, double* ylma,
 
 			cal_flag[ib][id]=true;
 			
+			std::vector<double> ylma;
 			//if(distance[id] > GridT.orbital_rmax) continue;
 			//	Ylm::get_ylm_real(this->nnn[it], this->dr[id], ylma);
 			if (distance < 1.0E-9) distance += 1.0E-9;
@@ -381,7 +382,6 @@ void Gint_k::calculate_charge(void)
 	int *at;
 	int *uc;
 	int *block_index;
-	double* ylma;
 	int* vindex;	
 	bool** cal_flag;
 	//int** AllOffset;
@@ -406,8 +406,6 @@ void Gint_k::calculate_charge(void)
 		{
 			nn = max(nn,(ucell.atoms[it].nwl+1)*(ucell.atoms[it].nwl+1));
 		}
-		ylma = new double[nn];
-		ZEROS(ylma, nn);
 
 		vindex = new int[pw.bxyz];
 		ZEROS(vindex, pw.bxyz);
@@ -442,7 +440,7 @@ void Gint_k::calculate_charge(void)
 				if(size==0) continue;
 				setVindex(ncyz, ibx, jby, kbz, vindex);
 				//timer::tick("Gint_k","cal_psir_ylm",'G');
-				cal_psir_ylm(size, grid_index, delta_r, ylma,
+				cal_psir_ylm(size, grid_index, delta_r,
 						at, uc, block_index, block_iw, block_size, 
 						cal_flag, psir_ylm);
 				//timer::tick("Gint_k","cal_psir_ylm",'G');
@@ -468,7 +466,6 @@ void Gint_k::calculate_charge(void)
 		delete[] block_index;
         delete[] cal_flag;
 
-		delete[] ylma;
 		delete[] vindex;
     }	
 
diff --git a/ABACUS.develop/source/src_lcao/gint_k_vl.cpp b/ABACUS.develop/source/src_lcao/gint_k_vl.cpp
index 49636e3b93..016da57c3b 100644
--- a/ABACUS.develop/source/src_lcao/gint_k_vl.cpp
+++ b/ABACUS.develop/source/src_lcao/gint_k_vl.cpp
@@ -128,7 +128,7 @@ inline int find_offset(const int size, const int grid_index,
 }
 
 inline void cal_psir_ylm(int size, int grid_index, double delta_r,
-						double** distance, double* ylma,
+						double** distance,
 						int* at, int* block_index, int* block_iw, int* block_size, 
 						bool** cal_flag, double** psir_ylm)
 {
@@ -182,6 +182,7 @@ inline void cal_psir_ylm(int size, int grid_index, double delta_r,
 
 			cal_flag[ib][id]=true;
 			
+			std::vector<double> ylma;
 			//if(distance[id] > GridT.orbital_rmax) continue;
 			//	Ylm::get_ylm_real(this->nnn[it], this->dr[id], ylma);
 			if (distance[ib][id] < 1.0E-9) distance[ib][id] += 1.0E-9;
@@ -368,7 +369,6 @@ void Gint_k::cal_vlocal_k(const double *vrs1, const Grid_Technique &GridT, const
 	int *block_size = nullptr; //band size: number of columns of a band
 	int *at = nullptr;
 	int *block_index = nullptr;
-	double* ylma = nullptr;
 	if(max_size!=0)
 	{
 		// save the small box information for a big box.
@@ -393,10 +393,6 @@ void Gint_k::cal_vlocal_k(const double *vrs1, const Grid_Technique &GridT, const
 			nn = std::max(nn, (ucell.atoms[it].nwl+1)*(ucell.atoms[it].nwl+1));
 		}
 
-		// ylma 
-		ylma = new double[nn];
-		ZEROS(ylma, nn);
-
 		for(int i=0; i<bxyz; i++)
 		{
 			// possible max atom number in a big box.
@@ -437,7 +433,7 @@ void Gint_k::cal_vlocal_k(const double *vrs1, const Grid_Technique &GridT, const
 				// get the wave functions in this
 				// grid.
 				//--------------------------------- 
-				cal_psir_ylm(size, grid_index, delta_r, distance, ylma,
+				cal_psir_ylm(size, grid_index, delta_r, distance,
 						at, block_index, block_iw, block_size, 
 						cal_flag, psir_ylm);
 
@@ -493,7 +489,6 @@ void Gint_k::cal_vlocal_k(const double *vrs1, const Grid_Technique &GridT, const
 		delete[] block_iw;
 		delete[] block_size;
 		delete[] block_index;
-		delete[] ylma;
 	}	
 
 	timer::tick("Gint_k","vlocal", 'K');
diff --git a/ABACUS.develop/source/src_pw/H_XC_pw.cpp b/ABACUS.develop/source/src_pw/H_XC_pw.cpp
index 2475ae17b0..c6914e4aa6 100644
--- a/ABACUS.develop/source/src_pw/H_XC_pw.cpp
+++ b/ABACUS.develop/source/src_pw/H_XC_pw.cpp
@@ -5,6 +5,7 @@
 double H_XC_pw::etxc;
 double H_XC_pw::vtxc;
 
+// [etxc, vtxc, v] = H_XC_pw::v_xc(...)
 std::tuple<double,double,matrix> H_XC_pw::v_xc
 (
 	const int &nrxx, // number of real-space grid
@@ -161,5 +162,5 @@ std::tuple<double,double,matrix> H_XC_pw::v_xc
     vt_xc *= omega / ncxyz;
 
     timer::tick("H_XC_pw","v_xc");
-    return std::make_tuple(et_xc,vt_xc,v);
+    return std::make_tuple(et_xc, vt_xc, std::move(v));
 }
diff --git a/ABACUS.develop/source/src_pw/H_XC_pw.h b/ABACUS.develop/source/src_pw/H_XC_pw.h
index bef5b032d7..baa15e5b7d 100644
--- a/ABACUS.develop/source/src_pw/H_XC_pw.h
+++ b/ABACUS.develop/source/src_pw/H_XC_pw.h
@@ -25,6 +25,7 @@ class H_XC_pw
 	static double vtxc;
 
 	// compute the exchange-correlation energy 
+	// [etxc, vtxc, v] = 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
diff --git a/ABACUS.develop/source/src_pw/potential.cpp b/ABACUS.develop/source/src_pw/potential.cpp
index a08a109e52..b16837e3db 100644
--- a/ABACUS.develop/source/src_pw/potential.cpp
+++ b/ABACUS.develop/source/src_pw/potential.cpp
@@ -290,9 +290,9 @@ matrix Potential::v_of_rho(
 //----------------------------------------------------------
 	
 	#ifdef USE_LIBXC
-    const auto etxc_vtxc_v = Potential_Libxc::v_xc(rho_in, CHR.rho_core);
+    const std::tuple<double,double,matrix> etxc_vtxc_v = Potential_Libxc::v_xc(rho_in, CHR.rho_core);
 	#else
-    const auto etxc_vtxc_v = H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, rho_in, CHR.rho_core);
+    const std::tuple<double,double,matrix> 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);
diff --git a/ABACUS.develop/source/src_pw/potential_libxc.cpp b/ABACUS.develop/source/src_pw/potential_libxc.cpp
index dc0cd45b5e..76c824cecb 100644
--- a/ABACUS.develop/source/src_pw/potential_libxc.cpp
+++ b/ABACUS.develop/source/src_pw/potential_libxc.cpp
@@ -13,6 +13,7 @@
 #include "src_global/global_function.h"
 #include "src_pw/xc_gga_pw.h"
 
+// [etxc, vtxc, v] = Potential_Libxc::v_xc(...)
 std::tuple<double,double,matrix> Potential_Libxc::v_xc(
 	const double * const * const rho_in,
 	const double * const rho_core_in)
@@ -30,13 +31,13 @@ std::tuple<double,double,matrix> Potential_Libxc::v_xc(
 	// use can check on website, for example:
 	// https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/
 	//----------------------------------------------------------
-	vector<xc_func_type> funcs = init_func();
+	std::vector<xc_func_type> funcs = init_func();
 
-	// the type of input_tmp is automatically set to 'tuple'
+	// the type of rho_sigma_gdr is automatically set to 'tuple'
 	// [rho, sigma, gdr] = cal_input( funcs, rho_in );
-	const auto input_tmp = cal_input( funcs, rho_in, rho_core_in );
-	const vector<double> &rho = std::get<0>(input_tmp);
-	const vector<double> &sigma = std::get<1>(input_tmp);
+	const auto rho_sigma_gdr = cal_input( funcs, rho_in, rho_core_in );
+	const std::vector<double> &rho = std::get<0>(rho_sigma_gdr);
+	const std::vector<double> &sigma = std::get<1>(rho_sigma_gdr);
 
 	for( xc_func_type &func : funcs )
 	{
@@ -45,9 +46,9 @@ std::tuple<double,double,matrix> Potential_Libxc::v_xc(
 		constexpr double grho_threshold = 1E-10;
 		xc_func_set_dens_threshold(&func, rho_threshold);
 		// sgn for threshold mask
-		const vector<double> sgn = [&]() -> vector<double>
+		const std::vector<double> sgn = [&]() -> std::vector<double>
 		{
-			vector<double> sgn( pw.nrxx * nspin0(), 1.0);
+			std::vector<double> sgn( pw.nrxx * nspin0(), 1.0);
 			if(nspin0()==2 && func.info->family != XC_FAMILY_LDA && func.info->kind==XC_CORRELATION)
 			{
 				for( size_t ir=0; ir!=pw.nrxx; ++ir )
@@ -61,9 +62,9 @@ std::tuple<double,double,matrix> Potential_Libxc::v_xc(
 			return sgn;
 		}();
 
-		vector<double> exc   ( pw.nrxx                       );
-		vector<double> vrho  ( pw.nrxx * nspin0()            );
-		vector<double> vsigma( pw.nrxx * ((1==nspin0())?1:3) );
+		std::vector<double> exc   ( pw.nrxx                       );
+		std::vector<double> vrho  ( pw.nrxx * nspin0()            );
+		std::vector<double> vsigma( pw.nrxx * ((1==nspin0())?1:3) );
 		
 		// cal etxc from rho, exc
 		auto process_exc = [&]()
@@ -93,7 +94,7 @@ std::tuple<double,double,matrix> Potential_Libxc::v_xc(
 				constexpr double vanishing_charge = 1.0e-12;
 				for( size_t ir=0; ir!=pw.nrxx; ++ir )
 				{
-					vector<double> v_tmp(4);
+					std::vector<double> v_tmp(4);
 					v_tmp[0] = e2 * (0.5 * (vrho[ir*2] + vrho[ir*2+1]));
 					const double vs = 0.5 * (vrho[ir*2] - vrho[ir*2+1]);
 					const double amag = sqrt( pow(rho_in[1][ir],2) 
@@ -120,9 +121,9 @@ std::tuple<double,double,matrix> Potential_Libxc::v_xc(
 		// cal vtxc, v from rho_in, rho, gdr, vsigma
 		auto process_vsigma = [&]()
 		{
-			const std::vector<std::vector<Vector3<double>>> &gdr = std::get<2>(input_tmp);
+			const std::vector<std::vector<Vector3<double>>> &gdr = std::get<2>(rho_sigma_gdr);
 			
-			vector<vector<Vector3<double>>> h( nspin0(), vector<Vector3<double>>(pw.nrxx) );
+			std::vector<std::vector<Vector3<double>>> h( nspin0(), std::vector<Vector3<double>>(pw.nrxx) );
 			if( 1==nspin0() )
 			{
 				for( size_t ir=0; ir!=pw.nrxx; ++ir )
@@ -142,7 +143,7 @@ std::tuple<double,double,matrix> Potential_Libxc::v_xc(
 			}
 
 			// define two dimensional array dh [ nspin, pw.nrxx ]
-			vector<vector<double>> dh(nspin0(), vector<double>(pw.nrxx));
+			std::vector<std::vector<double>> dh(nspin0(), std::vector<double>(pw.nrxx));
 			for( size_t is=0; is!=nspin0(); ++is )
 				GGA_PW::grad_dot( VECTOR_TO_PTR(h[is]), VECTOR_TO_PTR(dh[is]) );
 
@@ -215,7 +216,7 @@ std::tuple<double,double,matrix> Potential_Libxc::v_xc(
     vtxc *= ucell.omega / pw.ncxyz;
 
     timer::tick("Potential_Libxc","v_xc");
-	return std::make_tuple( etxc, vtxc, v );
+	return std::make_tuple( etxc, vtxc, std::move(v) );
 }
 
 
@@ -307,7 +308,7 @@ std::vector<xc_func_type> Potential_Libxc::init_func()
 
 
 
-// [rho, sigma, gdr] = cal_input( funcs, rho_in )
+// [rho, sigma, gdr] = Potential_Libxc::cal_input(...)
 std::tuple< std::vector<double>, 
 			std::vector<double>, 
 			std::vector<std::vector<Vector3<double>>> > 
@@ -363,7 +364,7 @@ Potential_Libxc::cal_input(
 			gdr.resize( nspin0() );
 			for( size_t is=0; is!=nspin0(); ++is )
 			{
-				vector<double> rhor(pw.nrxx);
+				std::vector<double> rhor(pw.nrxx);
 				for(int ir=0; ir<pw.nrxx; ++ir)
 					rhor[ir] = rho[ir*nspin0()+is];
 
@@ -371,7 +372,7 @@ Potential_Libxc::cal_input(
 				// initialize the charge density arry in reciprocal space
 				// bring electron charge density from real space to reciprocal space
 				//------------------------------------------
-				vector<complex<double>> rhog(pw.ngmc);
+				std::vector<std::complex<double>> rhog(pw.ngmc);
 				CHR.set_rhog(rhor.data(), rhog.data());
 
 				//-------------------------------------------
@@ -434,7 +435,7 @@ Potential_Libxc::cal_input(
 		}
 	}
 	
-	return std::make_tuple( rho, sigma, gdr );
+	return std::make_tuple( std::move(rho), std::move(sigma), std::move(gdr) );
 }
 
 #endif	//ifdef USE_LIBXC
diff --git a/ABACUS.develop/source/src_pw/potential_libxc.h b/ABACUS.develop/source/src_pw/potential_libxc.h
index 0328546af9..5ad20c79ec 100644
--- a/ABACUS.develop/source/src_pw/potential_libxc.h
+++ b/ABACUS.develop/source/src_pw/potential_libxc.h
@@ -23,8 +23,9 @@ class Potential_Libxc
 
 	//------------------------------------------------
 	// evaluate the exchange-correlation (XC) energy
-	// by using the input charge density rho_in
+	// by using the input charge density rho_in and rho_core_in
 	//------------------------------------------------
+	// [etxc, vtxc, v] = v_xc(...)
 	static std::tuple<double,double,matrix> v_xc(
 		const double * const * const rho_in,
 		const double * const rho_core_in);
@@ -41,7 +42,7 @@ class Potential_Libxc
 	// evaluate three quantities: rho, sigma, and gdr
 	// according to the input types of XC functionals
 	//------------------------------------------------
-	// [rho, sigma, gdr] = cal_input( funcs, rho_in )
+	// [rho, sigma, gdr] = cal_input(...)
 	static std::tuple< 
 		std::vector<double>, 
 		std::vector<double>,