diff --git a/ABACUS.develop/source/Makefile.Objects b/ABACUS.develop/source/Makefile.Objects
index caef758e51..90577331a7 100644
--- a/ABACUS.develop/source/Makefile.Objects
+++ b/ABACUS.develop/source/Makefile.Objects
@@ -30,9 +30,11 @@ dftd3_subroutine.o\
 vdwd3.o\
 vdwd3_parameters.o\
 pseudopot_upf.o\
+pseudopot_upf201.o\
+pseudopot_vwr.o\
 pseudo_nc.o \
-pseudopot_cell_vl.o \
-pseudopot_cell_vnl.o \
+VL_in_pw.o\
+VNL_in_pw.o\
 atom_pseudo.o\
 unitcell_pseudo.o\
 threshold_elec.o\
@@ -230,7 +232,7 @@ H_XC_pw.o \
 H_TDDFT_pw.o\
 read_rho.o\
 read_atoms.o\
-read_pseudopot.o\
+read_cell_pseudopots.o\
 read_dm.o\
 write_pot.o\
 write_rho.o\
diff --git a/ABACUS.develop/source/src_global/matrix3.h b/ABACUS.develop/source/src_global/matrix3.h
index 966e209030..8989a3ad38 100644
--- a/ABACUS.develop/source/src_global/matrix3.h
+++ b/ABACUS.develop/source/src_global/matrix3.h
@@ -1,7 +1,3 @@
-//==========================================================
-// AUTHOR: Lixin He, Mohan Chen
-// LAST UPDATE : 2009-03-24 add operator == and !=
-//==========================================================
 #ifndef MatriX3_H
 #define MatriX3_H
 
diff --git a/ABACUS.develop/source/src_io/read_pseudopot.cpp b/ABACUS.develop/source/src_io/read_cell_pseudopots.cpp
similarity index 97%
rename from ABACUS.develop/source/src_io/read_pseudopot.cpp
rename to ABACUS.develop/source/src_io/read_cell_pseudopots.cpp
index 8ddb8c1b47..88db649f47 100644
--- a/ABACUS.develop/source/src_io/read_pseudopot.cpp
+++ b/ABACUS.develop/source/src_io/read_cell_pseudopots.cpp
@@ -7,9 +7,9 @@
 //==========================================================
 // Read pseudopotential according to the dir
 //==========================================================
-void UnitCell_pseudo::read_pseudopot(const string &pp_dir)
+void UnitCell_pseudo::read_cell_pseudopots(const string &pp_dir)
 {
-	TITLE("UnitCell_pseudo","read_pseudopot");
+	TITLE("UnitCell_pseudo","read_cell_pseudopots");
 	// setup reading log for pseudopot_upf
 	stringstream ss;
 	ss << global_out_dir << "atom_pseudo.log";
diff --git a/ABACUS.develop/source/src_lcao/ORB_control.cpp b/ABACUS.develop/source/src_lcao/ORB_control.cpp
index b1dc36c8aa..a094b42411 100644
--- a/ABACUS.develop/source/src_lcao/ORB_control.cpp
+++ b/ABACUS.develop/source/src_lcao/ORB_control.cpp
@@ -1,11 +1,11 @@
-#include "../src_pw/tools.h"
+//#include "../src_pw/tools.h"
 #include "ORB_control.h"
 //#include "src_global/sltk_atom_arrange.h"
 #include "ORB_gen_tables.h"
 #include "build_st_pw.h"
 #include "../src_pdiag/pdiag_double.h"
 
-#include "../src_pw/global.h"
+//#include "../src_pw/global.h"
 
 ORB_control::ORB_control()
 {}
diff --git a/ABACUS.develop/source/src_lcao/ORB_read.cpp b/ABACUS.develop/source/src_lcao/ORB_read.cpp
index 6a9b1eadaa..653eb01540 100644
--- a/ABACUS.develop/source/src_lcao/ORB_read.cpp
+++ b/ABACUS.develop/source/src_lcao/ORB_read.cpp
@@ -102,7 +102,7 @@ void LCAO_Orbitals::bcast_files(const int &ntype_in)
 #endif
 
 
-#include "../src_pw/global.h"
+//#include "../src_pw/global.h"
 void LCAO_Orbitals::Read_Orbitals(
 	const int &ntype_in, 
 	const int &lmax_in,
@@ -813,7 +813,12 @@ void LCAO_Orbitals::Read_Descriptor(void)	//read descriptor basis
 }
 
 
-void LCAO_Orbitals::read_orb_file(ifstream &ifs,const int &it, int &lmax, int &nchimax, Numerical_Orbital* ao)
+void LCAO_Orbitals::read_orb_file(
+	ifstream &ifs,
+	const int &it, 
+	int &lmax, 
+	int &nchimax, 
+	Numerical_Orbital* ao)
 {
 	TITLE("LCAO_Orbitals","read_orb_file");
 	char word[80];
@@ -1048,4 +1053,4 @@ void LCAO_Orbitals::read_orb_file(ifstream &ifs,const int &it, int &lmax, int &n
 	
 	delete[] nchi;
 	return;
-}
\ No newline at end of file
+}
diff --git a/ABACUS.develop/source/src_pw/pseudopot_cell_vl.cpp b/ABACUS.develop/source/src_pw/VL_in_pw.cpp
similarity index 99%
rename from ABACUS.develop/source/src_pw/pseudopot_cell_vl.cpp
rename to ABACUS.develop/source/src_pw/VL_in_pw.cpp
index 8c879bdde9..0bc3d9f18b 100644
--- a/ABACUS.develop/source/src_pw/pseudopot_cell_vl.cpp
+++ b/ABACUS.develop/source/src_pw/VL_in_pw.cpp
@@ -1,7 +1,7 @@
 #include "global.h"
 #include "tools.h"
 #include "../src_pw/myfunc.h"
-#include "pseudopot_cell_vl.h"
+#include "VL_in_pw.h"
 #include "../src_global/math_integral.h"
 
 pseudopot_cell_vl::pseudopot_cell_vl()
diff --git a/ABACUS.develop/source/src_pw/pseudopot_cell_vl.h b/ABACUS.develop/source/src_pw/VL_in_pw.h
similarity index 87%
rename from ABACUS.develop/source/src_pw/pseudopot_cell_vl.h
rename to ABACUS.develop/source/src_pw/VL_in_pw.h
index 7541e041d0..873df18f62 100644
--- a/ABACUS.develop/source/src_pw/pseudopot_cell_vl.h
+++ b/ABACUS.develop/source/src_pw/VL_in_pw.h
@@ -1,5 +1,5 @@
-#ifndef pseudopot_cell_vl_H
-#define pseudopot_cell_vl_H
+#ifndef VL_IN_PW_H 
+#define VL_IN_PW_H 
 
 #include "tools.h"
 
@@ -29,4 +29,4 @@ class pseudopot_cell_vl
 
 };
 
-#endif // PSEUDOPOT_CELL_VL_H
+#endif // VL_IN_PW 
diff --git a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp b/ABACUS.develop/source/src_pw/VNL_in_pw.cpp
similarity index 98%
rename from ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp
rename to ABACUS.develop/source/src_pw/VNL_in_pw.cpp
index a492bb00e5..e419dca1db 100644
--- a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp
+++ b/ABACUS.develop/source/src_pw/VNL_in_pw.cpp
@@ -1,9 +1,5 @@
-//==========================================================
-// AUTHOR : Lixin He, mohan
-// DATE : 2008-11-08
-//==========================================================
 #include "global.h"
-#include "pseudopot_cell_vnl.h"
+#include "VNL_in_pw.h"
 #include "tools.h"
 #include "wavefunc.h"
 #include "../src_lcao/ORB_gen_tables.h"
diff --git a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h b/ABACUS.develop/source/src_pw/VNL_in_pw.h
similarity index 95%
rename from ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h
rename to ABACUS.develop/source/src_pw/VNL_in_pw.h
index d9cd3c82d0..ebca083bc3 100644
--- a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h
+++ b/ABACUS.develop/source/src_pw/VNL_in_pw.h
@@ -1,8 +1,8 @@
-#ifndef PSEUDOPOT_CELL_VNL_H
-#define PSEUDOPOT_CELL_VNL_H
+#ifndef VNL_IN_PW_H 
+#define VNL_IN_PW_H 
 
 #include "tools.h"
-#include "pseudopot_cell_vl.h"
+#include "VL_in_pw.h"
 #include "../src_lcao/ORB_gen_tables.h"
 #include "wavefunc_in_pw.h"
 #include "unitcell_pseudo.h"
@@ -92,4 +92,4 @@ class pseudopot_cell_vnl: public pseudopot_cell_vl
 	ORB_gaunt_table MGT;
 };
 
-#endif // PSEUDOPOT_CELL_VNL_H
+#endif // VNL_IN_PW 
diff --git a/ABACUS.develop/source/src_pw/energy.cpp b/ABACUS.develop/source/src_pw/energy.cpp
index 8e77ff7f89..5eeae37758 100644
--- a/ABACUS.develop/source/src_pw/energy.cpp
+++ b/ABACUS.develop/source/src_pw/energy.cpp
@@ -1,7 +1,3 @@
-//==========================================================
-// AUTHOR : Lixin He, mohan
-// DATE : 2008-11-21
-//==========================================================
 #include "tools.h"
 #include "global.h"
 #include "energy.h"
diff --git a/ABACUS.develop/source/src_pw/global.cpp b/ABACUS.develop/source/src_pw/global.cpp
index 7c9a3c1871..96f7c28c10 100644
--- a/ABACUS.develop/source/src_pw/global.cpp
+++ b/ABACUS.develop/source/src_pw/global.cpp
@@ -1,9 +1,3 @@
-//==========================================================
-// Author: Lixin He,mohan
-// DATE : 2008-11-6
-// LAST MODIFIED : 2021-01-31 by Mohan
-//==========================================================
-
 #include "global.h"
 //----------------------------------------------------------
 // init "GLOBAL CLASS" object
diff --git a/ABACUS.develop/source/src_pw/global.h b/ABACUS.develop/source/src_pw/global.h
index c59edca018..d7dd83e170 100644
--- a/ABACUS.develop/source/src_pw/global.h
+++ b/ABACUS.develop/source/src_pw/global.h
@@ -6,7 +6,7 @@
 #include "../src_global/global_function.h"
 #include "pw_basis.h"
 #include "energy.h"
-#include "pseudopot_cell_vnl.h"
+#include "VNL_in_pw.h"
 #include "charge_broyden.h"
 #include "potential.h"
 #include "xc_type.h"
diff --git a/ABACUS.develop/source/src_pw/pseudo_vl.cpp b/ABACUS.develop/source/src_pw/pseudo_vl.cpp
deleted file mode 100644
index 5885eab99d..0000000000
--- a/ABACUS.develop/source/src_pw/pseudo_vl.cpp
+++ /dev/null
@@ -1,39 +0,0 @@
-/* pseudo_vl.cpp */
-
-#include "pseudo_vl.h"
-#include "global.h"
-
-pseudo_vl::pseudo_vl()
-{
-	vloc_at =  new double[1];
-}
-
-pseudo_vl::~pseudo_vl()
-{
-	delete[] vloc_at;
-}
-
-//---------------------------------------------------------------------
-void pseudo_vl::set_pseudo_vl(const Pseudopot_upf &upf)
-{
-//	cout << "\n === set_pseudo_vl() === ";
-	this->set_pseudo_at(upf);
-
-	delete[] vloc_at;
-	vloc_at = new double[mesh];
-	assert(vloc_at != 0);
-
-	for (int i = 0;i < mesh;i++)
-	{
-		vloc_at[i] = upf.vloc[i];
-	}
-//	cout << "\n End set_pseudo_vl() " << endl;
-} 
-
-void pseudo_vl::print_pseudo_vl(ofstream &ofs)
-{
-	ofs << "\n pseudo_vl:";
-	print_pseudo_at(ofs);
-	out.printr1_d(ofs, "vloc_at : ", vloc_at, mesh);
-	ofs << "\n ----------------------------------- ";
-}
diff --git a/ABACUS.develop/source/src_pw/pseudo_vl.h b/ABACUS.develop/source/src_pw/pseudo_vl.h
deleted file mode 100644
index 6bf761f881..0000000000
--- a/ABACUS.develop/source/src_pw/pseudo_vl.h
+++ /dev/null
@@ -1,21 +0,0 @@
-#ifndef PSEUDO_VL_H
-#define PSEUDO_VL_H
-
-#include "pseudopot_upf.h"
-#include "pseudo_atom.h"
-
-class pseudo_vl: public pseudo_atom
-{
-public:
-
-
-	double *vloc_at;	// [mesh], local potential( = pseudopot_upf.vloc )
-
-	pseudo_vl();
-	~pseudo_vl();
-
-	void set_pseudo_vl(const Pseudopot_upf &upf);
-	void print_pseudo_vl(ofstream &ofs);
-};
-
-#endif // PSEUDO_VL_H
diff --git a/ABACUS.develop/source/src_pw/pseudopot_upf.cpp b/ABACUS.develop/source/src_pw/pseudopot_upf.cpp
index 55db8ffd93..a857f4b8f4 100644
--- a/ABACUS.develop/source/src_pw/pseudopot_upf.cpp
+++ b/ABACUS.develop/source/src_pw/pseudopot_upf.cpp
@@ -7,7 +7,6 @@
 #include <sstream>
 #include <cstring> // Peize Lin fix bug about strcpy 2016-08-02
 
-int Number[2];
 
 using namespace std;
 
@@ -144,412 +143,6 @@ string Pseudopot_upf::trimend(string &in_str)
 
 
 
-//----------------------------------------------------------
-// This code is used to read in vwr pseudopotential format,
-// Now we only use LDA, so if PBE or other functionals are used, 
-// one needs to change the following code. The vwr format 
-// needs we to generate NL projectors by ourself. 
-// One way to check this is correct
-// is to use opium to generate .ncpp pseudopotential first, 
-// which contains the same informaiton in vwr, then we had
-// both UPF format from upftools in Quantum Espresso and
-// we can write a short code to transform ncpp to vwr.
-// Then compare the two results.
-// mohan 2013-05-25
-//-----------------------------------------------------------
-int Pseudopot_upf::read_pseudo_vwr(ifstream &ifs)
-{
-	ofs_running << " -------------------------------------------------" << endl;
-	cout << " READ IN VWR TYPE PSEUDOPOTENTIALS." << endl;
-	ofs_running << " Read in vwr type pseudopotentials " << endl;
-
-
-	// --------------------------------------
-	// (1) read in data 
-	// --------------------------------------
-	this->dft[0]="SLA";
-	this->dft[1]="PZ";
-	this->dft[2]="NOGX";
-	this->dft[3]="NOGC";
-	this->pp_type="NC";
-	this->tvanp=false;
-	ofs_running << " Always use PZ-LDA by now." << endl;
-	
-
-	// (1) read in mesh
-	string value;
-	int length=0;
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	mesh = std::atoi( value.c_str() );
-	//the mesh should be odd, which is forced in Simpson integration
-	if(mesh%2==0) 
-	{
-		mesh=mesh-1;
-		ofs_running << " Mesh number - 1, we need odd number, \n this may affect some polar atomic orbitals." << endl;
-	}
-	ofs_running << setw(15) << "MESH" << setw(15) << mesh << endl;
-	// (2) read in nlcc: nonlinear core correction
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	nlcc = std::atoi( value.c_str() );
-	ofs_running << setw(15) << "NLCC" << setw(15) << nlcc << endl;
-	// (3) iatom : index for atom 
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	psd = value;
-	ofs_running << setw(15) << "ATOM" << setw(15) << psd << endl;
-	// (4) valence electron number
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	zp = std::atoi( value.c_str() );
-	ofs_running << setw(15) << "Z(VALENCE)" << setw(15) << zp << endl;
-	// (5) spd_loc, which local pseudopotential should I choose
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	spd_loc = std::atoi( value.c_str() );
-	ofs_running << setw(15) << "LOC(spd)" << setw(15) << spd_loc << endl;
-	// (6) read in the occupations
-	double* tmp_oc = new double[3];
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	tmp_oc[0]= std::atoi( value.c_str() );
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	tmp_oc[1]= std::atoi( value.c_str() );
-	ifs >> value; length = value.find(","); value.erase(length,1);
-	tmp_oc[2]= std::atoi( value.c_str() );
-	ofs_running << setw(15) << "OCCUPATION" << setw(15) << tmp_oc[0] 
-	<< setw(15) << tmp_oc[1] << setw(15) << tmp_oc[2] << endl;
-	// (7) spin orbital
-	ifs >> has_so;
-
-
-	// label to count the projector or atomic wave functions	
-	getline(ifs,value);
-	int iref_s, iref_p, iref_d;
-	ifs >> iref_s >> iref_p >> iref_d;
-	ofs_running << setw(15) << "Vnl_USED" << setw(15) << iref_s 
-	<< setw(15) << iref_p << setw(15) << iref_d << endl;
-	if(spd_loc==1) iref_s=0;
-	else if(spd_loc==2) iref_p=0;
-	else if(spd_loc==3) iref_d=0;
-	ifs >> iTB_s >> iTB_p >> iTB_d;
-	ofs_running << setw(15) << "Orb_USED" << setw(15) << iTB_s 
-	<< setw(15) << iTB_p << setw(15) << iTB_d << endl;
-	
-	
-	// calculate the number of wave functions
-	this->nwfc = 0;
-	if(iTB_s) ++nwfc;
-	if(iTB_p) ++nwfc;
-	if(iTB_d) ++nwfc;
-	ofs_running << setw(15) << "NWFC" << setw(15) << nwfc << endl;
-	// allocate occupation number array for wave functions
-	delete[] this->oc;
-	delete[] els;
-	this->oc = new double[nwfc];
-	els = new string[nwfc];
-	// set the value of occupations
-	delete[] lchi;
-	lchi = new int[nwfc];
-	int iwfc=0;
-	if(iTB_s){oc[iwfc]=tmp_oc[0];lchi[iwfc]=0;els[iwfc]="S";++iwfc;}
-	if(iTB_p){oc[iwfc]=tmp_oc[1];lchi[iwfc]=1;els[iwfc]="P";++iwfc;}
-	if(iTB_d){oc[iwfc]=tmp_oc[2];lchi[iwfc]=2;els[iwfc]="D";++iwfc;}
-	delete[] tmp_oc;
-	getline(ifs,value);
-
-
-	// global variables that will be used
-	// in other classes.
-	delete[] r;
-	delete[] rab;
-	delete[] vloc;
-	this->r = new double[mesh];
-	this->rab = new double[mesh];
-	this->vloc = new double[mesh];
-	ZEROS(r,mesh);
-	ZEROS(rab,mesh);
-	ZEROS(vloc,mesh);
-	delete[] rho_at;
-	delete[] rho_atc;
-	rho_at = new double[mesh];
-	rho_atc = new double[mesh];
-	ZEROS(rho_at, mesh);
-	ZEROS(rho_atc, mesh);
-	// local variables in this function
-	this->vs = new double[mesh];
-	this->vp = new double[mesh];
-	this->vd = new double[mesh];
-	this->ws = new double[mesh];
-	this->wp = new double[mesh];
-	this->wd = new double[mesh];
-	ZEROS(vs,mesh);
-	ZEROS(vp,mesh);
-	ZEROS(vd,mesh);
-	ZEROS(ws,mesh);
-	ZEROS(wp,mesh);
-	ZEROS(wd,mesh);
-	string line;
-	if(spd_loc>0 && nlcc==0)
-	{
-		for(int ir=0; ir<mesh; ++ir)
-		{
-			// it's an interesting question whether
-			// ws[ir] has 1/sqrt(4pi)
-			ifs >> r[ir] >> vs[ir] >> vp[ir] >> vd[ir] 
-				>> ws[ir] >> wp[ir] >> wd[ir];
-			getline(ifs, line);
-		}
-	}
-	else if(spd_loc==0 && nlcc==0)
-	{
-		for(int ir=0; ir<mesh; ++ir)
-		{
-			ifs >> r[ir] >> vs[ir] >> vp[ir] >> vd[ir] 
-				>> ws[ir] >> wp[ir] >> wd[ir] >> vloc[ir];
-			getline(ifs, line);
-		}
-	}
-	else if(spd_loc>0 && nlcc==1)
-	{
-		for(int ir=0; ir<mesh; ++ir)
-		{
-			ifs >> r[ir] >> vs[ir] >> vp[ir] >> vd[ir] 
-				>> ws[ir] >> wp[ir] >> wd[ir] >> rho_atc[ir];
-			getline(ifs, line);
-		}
-	}
-	else if(spd_loc==0 && nlcc==1)
-	{
-		for(int ir=0; ir<mesh; ++ir)
-		{
-			ifs >> r[ir] >> vs[ir] >> vp[ir] >> vd[ir] 
-				>> ws[ir] >> wp[ir] >> wd[ir] >> vloc[ir] >> rho_atc[ir];
-			getline(ifs, line);
-		}
-	}
-	// Hartree to Rydberg
-	for(int ir=0; ir<mesh; ++ir)
-	{
-		vs[ir] *= 2.0;
-		vp[ir] *= 2.0;
-		vd[ir] *= 2.0;
-		vloc[ir] *= 2.0;
-	}
-
-
-	// --------------------------------------
-	// (2) check unit
-	// --------------------------------------
-	// calculate rab;
-	// rab may not be accurate enough
-	rab[0] = r[0];
-	for(int ir=1; ir<mesh-1; ++ir)
-	{
-		rab[ir]=(r[ir+1]-r[ir-1])/2.0;
-	}
-	// check unit of vs, vp, vd
-	double units = 0.0;
-	double unitp = 0.0;
-	double unitd = 0.0;
-	for(int ir=1; ir<mesh-1; ++ir)
-	{
-		double dr = (r[ir+1]-r[ir-1])/2.0;
-		units += ws[ir] * ws[ir] * r[ir] * r[ir] * dr;
-		unitp += wp[ir] * wp[ir] * r[ir] * r[ir] * dr;
-		unitd += wd[ir] * wd[ir] * r[ir] * r[ir] * dr;
-	}
-	ofs_running << setw(15) << "WFC_UNIT" << setw(15) << units 
-	<< setw(15) << unitp << setw(15) << unitd << endl; 
-
-
-	// because only the rank=0 procesor read the pseudopotential
-	// information, in order to make all the processors to stop
-	// the job, we need to return the error information first.
-	// we need to choose a threshold for the deviation of the
-	// norm of pseudo atomic orbitals, I set 0.2
-	// mohan 2013-06-28
-	if( abs(units-1.0) > 0.2 && (iTB_s==1 || iref_s==1)) {return 3;}
-	if( abs(unitp-1.0) > 0.2 && (iTB_p==1 || iref_p==1)) {return 3;}
-	if( abs(unitd-1.0) > 0.2 && (iTB_d==1 || iref_d==1)) {return 3;}
-
-
-	// calculate the phi*r*sqrt(4pi)
-	this->chi.create(nwfc,mesh);
-	for(int ir=0; ir<mesh; ++ir)
-	{
-		int iwfc=0;
-		if(iTB_s==1){chi(iwfc,ir) = ws[ir]*r[ir];++iwfc;}
-		if(iTB_p==1){chi(iwfc,ir) = wp[ir]*r[ir];++iwfc;}
-		if(iTB_d==1){chi(iwfc,ir) = wd[ir]*r[ir];++iwfc;}
-	}
-	// rho atom 
-	for(int ir=0; ir<mesh; ++ir)
-	{
-		for(int iwfc=0; iwfc<nwfc; ++iwfc)
-		{
-			rho_at[ir] += oc[iwfc]*chi(iwfc,ir)*chi(iwfc,ir);
-		}
-	}
-
-
-
-	// --------------------------------------
-	// (4) local pseudopotential 
-	// --------------------------------------
-	if(spd_loc==0) for(int ir=0; ir<mesh; ++ir)
-	{
-			// do nothing	
-	}
-	else if(spd_loc==1) for(int ir=0; ir<mesh; ++ir) vloc[ir] = vs[ir]; 
-	else if(spd_loc==2) for(int ir=0; ir<mesh; ++ir) vloc[ir] = vp[ir];
-	else if(spd_loc==3) for(int ir=0; ir<mesh; ++ir) vloc[ir] = vd[ir];
-	else if(spd_loc==12) for(int ir=0; ir<mesh; ++ir) vloc[ir] = (vs[ir]+vp[ir])/2.0;
-	else if(spd_loc==13) for(int ir=0; ir<mesh; ++ir) vloc[ir] = (vs[ir]+vd[ir])/2.0;
-	else if(spd_loc==23) for(int ir=0; ir<mesh; ++ir) vloc[ir] = (vp[ir]+vd[ir])/2.0;
-
-
-	// --------------------------------------
-	// (5) setup nonlocal pseudopotentials
-	// --------------------------------------
-	// for non-local pseudopotentials.
-	if(iref_d==1) lmax=2;
-	else if(iref_p==1) lmax=1;
-	else if(iref_s==1) lmax=0;
-	else
-	{
-		cout << "\n !!! READ THIS FIRST !!!" << endl;
-		cout << " Could not decide which is the max angular momentum from .vwr pseudopotential file." << endl;
-		cout << " No reference states in .vwr pseudopotential file." << endl;
-		cout << " That's incorrect, please check the refenrece states in .vwr file.";
-		cout << "\n !!! READ THIS FIRST !!!" << endl;
-		return 3;	
-	}
-	// no projectors now
-	this->nbeta = 0;
-	if(iref_s==1) ++nbeta; // add one s projector
-	if(iref_p==1) ++nbeta; // add one p projector
-	if(iref_d==1) ++nbeta; // add one p projector
-	ofs_running << setw(15) << "NPROJ" << setw(15) << nbeta << endl;
-	this->nd = this->nbeta;
-	ofs_running << setw(15) << "N-Dij" << setw(15) << nd << endl;
-	// calculate the angular momentum for each beta
-	delete[] lll;
-	lll = new int[nbeta]; 
-	int icount=0;
-	if(iref_s==1) {lll[icount]=0; ++icount;}// s projector
-	if(iref_p==1) {lll[icount]=1; ++icount;}// p projector
-	if(iref_d==1) {lll[icount]=2; ++icount;}// p projector
-	for(int i=0; i<nbeta; ++i)
-	{
-		ofs_running << " lll[" << i << "]=" << lll[i] << endl;
-	}
-	// kkbeta(nbeta): number of mesh points for projector i (must be .le.mesh )
-	delete[] kkbeta;
-	kkbeta = new int[nbeta];
-	for(int ib=0; ib<nbeta; ++ib)
-	{
-		kkbeta[ib] = mesh;
-	}
-	// nonlocal projector
-	beta.create(nbeta,mesh);
-	// coefficients
-	dion.create(nbeta,nbeta);
-
-
-	// --------------------------------------
-	// (6) generate nonlocal pseudopotentials
-	// --------------------------------------
-	// tmp function to evaluate < beta | delta_v | beta>
-	double* func = new double[mesh];
-	// tmp value (vs, vp or vd)
-	double* vl = new double[mesh];
-	// tmp wave function (ws, wp or wd with r)
-	double* wlr = new double[mesh];
-	ZEROS(func, mesh);
-	ZEROS(vl, mesh);
-	ZEROS(wlr, mesh);
-	double rcut = 5.0/1.03;
-	ofs_running << setw(15) << "RCUT_NL" << setw(15) << rcut << endl;
-	for(int ib=0; ib<nbeta; ++ib)
-	{
-		double coef = 0.0;
-		const int lnow = lll[ib];
-		if(lnow==0) for(int ir=0; ir<mesh; ++ir){vl[ir]=vs[ir]; wlr[ir]=ws[ir]*r[ir];}
-		else if(lnow==1) for(int ir=0; ir<mesh; ++ir){vl[ir]=vp[ir]; wlr[ir]=wp[ir]*r[ir];}
-		else if(lnow==2) for(int ir=0; ir<mesh; ++ir){vl[ir]=vd[ir]; wlr[ir]=wd[ir]*r[ir];}
-		// for non-local projectors
-		// note that < phi | dV | phi > integration must have 4pi,
-		// this 4pi is also needed in < phi | phi > = 1 integration.
-		// However, this phi has sqrt(sphi) already because I
-		// found < phi | phi > = 1 directly.
-		ofs_running << " Projector index = " << ib+1 << ", L = " << lnow << endl;
-		for(int ir=2; ir<mesh-1; ++ir)
-		{
-			// p nl
-			beta(ib,ir)=(vl[ir]-vloc[ir])*wlr[ir];
-			if(r[ir]<rcut)
-			{
-				coef=coef+(vl[ir]-vloc[ir])*wlr[ir]*wlr[ir]*(r[ir+1]-r[ir-1])/2.0;
-			}
-		}
-
-
-// In pw they did this:
-//		dion(ib,ib)=1.0/coef;
-
-        if(coef<0.0) dion(ib,ib) = -1.0;
-        if(coef>=0.0) dion(ib,ib) = 1.0;
-		//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-		// suppose wave function have sqrt(4pi) already
-		//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-        coef=1.0/sqrt(abs(coef));
-		ofs_running << setw(25) << "1/sqrt(<phi|deltaV|phi>)" << setw(15) << coef << endl;
-		for(int ir=0; ir<mesh; ++ir)
-		{
-			beta(ib,ir) *= coef;
-			// --------- FOR TEST ---------
-			if(ib>2)
-			{
-//				beta(ib,ir) *= 0.0; // for test, disable Non-local
-			}
-			// --------- FOR TEST ---------
-		}
-		
-	}
-
-	// print out the projector.
-	/*
-	ofs_running << " Nonlocal projector : " << endl;
-	for(int ir=0; ir<mesh; ++ir)
-	{
-		if(r[ir]<rcut)
-		{
-			ofs_running << " " << setw(15) << r[ir]; 
-			for(int ib=0; ib<nbeta; ++ib)
-			{
-				ofs_running << setw(25) << beta(ib,ir);
-			}
-			ofs_running << endl;
-		}
-	}
-	*/
-	ofs_running << " -------------------------------------------------" << endl;
-
-
-	// --------------------------------------
-	// (5) clean up 
-	// --------------------------------------
-	delete[] func;
-	delete[] vl;
-	delete[] wlr;
-
-	delete[] vs;
-	delete[] vp;
-	delete[] vd;
-
-	delete[] ws;
-	delete[] wp;
-	delete[] wd;
-
-	return 0;
-}
-
 //  read pseudopot_upf potential "upf" in the Unified
 //  Pseudopot_upfpotential Format
 int Pseudopot_upf::read_pseudo_upf(ifstream &ifs)
@@ -1009,794 +602,6 @@ void Pseudopot_upf::print_pseudo_upf(ofstream &ofs)
 }
 
 
-// added by zhangwenshuai
-int Pseudopot_upf::read_pseudo_upf201(ifstream &ifs)
-{
-    string dummy, word;
-    int i, j, ir, ONCVPSP;
-
-	while (ifs.good())
-	{
-		ifs >> dummy;
-		// We start from PP_Header
-		if(dummy=="<PP_HEADER")
-		{
-			// Read header
-			READ_VALUE(ifs, word);   // generated
-			READ_VALUE(ifs, word);   // author
-			READ_VALUE(ifs, word);   // date
-			READ_VALUE(ifs, word);   // comment
-
-			READ_VALUE(ifs, word);   // element
-			stringstream wdsstream(word);
-			getline(wdsstream,this->psd,'"'); 
-			getline(wdsstream,this->psd,'"'); 
-
-			ifs >> word;   // pseudo_type
-			if(word == "pseudo_type=\"")
-			{
-				ifs >> word;
-				get_char(word);
-				this->pp_type = word.substr(0,Number[0]);
-			}
-			else
-			{
-				get_char(word);
-				this->pp_type = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
-			}
-
-			if(pp_type!="NC") 
-			{
-				WARNING_QUIT("Pseudopot_upf::read_pseudo_header","unknown pseudo type");
-			}
-
-			READ_VALUE(ifs, word);   // relativistic
-			READ_VALUE(ifs, word);   // is_ultrasoft
-			if ( word.find("\"T\"") < word.length() ) // zws add 20160108
-			{
-				cout << "\n WARNING: ULTRASOFT PSEUDOPOTENTIAL IS NOT SUPPORTED !!! \n" << endl;
-			}
-			READ_VALUE(ifs, word);   // is_paw
-			if ( word.find("\"T\"") < word.length() )
-			{
-				cout << "\n WARNING: PAW PSEUDOPOTENTIAL IS NOT SUPPORTED !!! \n" << endl;
-			}
-
-			READ_VALUE(ifs, word);   // is_coulomb
-			ifs >> word;   // has_so
-			string so;
-
-			if(word == "has_so=\"")
-			{
-				ifs >> word;
-				get_char(word);
-				so = word.substr(0,Number[0]);
-			}
-			else
-			{
-				get_char(word);
-				so = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
-			}
-
-			if (so == "T")
-			{
-				this->has_so = true;
-			}
-			else
-			{
-				this->has_so = false;
-			}
-
-			READ_VALUE(ifs, word);   // has_wfc
-			READ_VALUE(ifs, word);   // has_gipaw
-
-			string nlc;
-			//char p[13] = "paw_as_gipaw";
-			ifs >> word;             // paw_as_gipaw?
-			//cout << "word.substr(0,30) = " << word.substr(0,30) << "."<< endl;
-			if( word.substr(0,13) == "paw_as_gipaw" )
-			{
-				ONCVPSP = 0;
-				ifs >> word;     // core_correction
-				if(word == "core_correction=\"")
-				{
-					ifs >> word;
-					get_char(word);
-					nlc = word.substr(0,Number[0]);
-				}
-				else
-				{
-					get_char(word);
-					nlc = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
-				}
-
-			}
-			else
-			{
-				ONCVPSP = 1; // Generated using ONCVPSP code by D. R. Hamann, SG15 DOJO
-				if(word == "core_correction=\"")
-				{
-					ifs >> word;
-					get_char(word);
-					nlc = word.substr(0,Number[0]);
-				}
-				else
-				{
-					get_char(word);
-					nlc = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
-				}
-			}
-
-			//cout << "nlc = " << nlc << endl;
-
-			if (nlc == "T")
-			{
-				this->nlcc = true;
-			}
-			else
-			{
-				this->nlcc = false;
-			}
-
-			READ_VALUE(ifs, word);   // functional
-			//cout << "word = " << word << endl;
-			//                        this->dft[0]="SLA";
-			//                        this->dft[1]="PZ";
-			//                        this->dft[2]="NOGX";
-			//                        this->dft[3]="NOGC";
-
-			string funcstr;  //{zws 01-06-16
-			wdsstream.str("");
-			wdsstream.clear();
-			wdsstream << word;
-			for ( int idft = 0; idft < 2; idft++)
-			{
-				getline(wdsstream,funcstr,'"');
-			}
-			wdsstream.str("");
-			wdsstream.clear();
-			wdsstream << funcstr;
-
-			for( int idft = 0; idft < 4; idft++ )
-			{
-				getline(wdsstream,dft[idft],'-');
-			}
-
-			do 
-			{
-				getline(ifs, word);
-				//cout << "word       = " << word << endl;
-				word.erase(0,word.find_first_not_of(" ") );
-				word.erase(word.find_last_not_of(" ")+1 );
-				//word = trim(word);
-				//cout << "trim(word) = " << word << endl;
-				get_char(word);
-				//cout << " Number = " << Number[0] << ", " << Number[1] << endl;
-				//cout << word.substr(0,Number[0])  << "__" << word.substr(Number[0]+1, Number[1]-Number[0]-1) << endl;
-				dummy = word.substr(0,Number[0]) ;
-				//cout << " dummy = " << dummy << endl;
-				if( dummy == "z_valence=" ) 
-				{
-					this->zp = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-				else if ( dummy == "total_psenergy=" ) 
-				{
-					this->etotps = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-				else if ( dummy == "rho_cutoff=" )
-				{
-				}
-				else if ( dummy == "wfc_cutoff=" ) 
-				{
-				}
-				else if ( dummy == "l_max=" ) 
-				{
-					this->lmax = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-				else if ( dummy == "mesh_size=" ) 
-				{
-					this->mesh = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-				else if ( dummy == "number_of_wfc=" ) 
-				{
-					this->nwfc = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-				else if ( dummy == "number_of_proj=" ) 
-				{
-					this->nbeta = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-				//break;
-
-			}while( word.substr(word.length()-1, 1) !=">" ); 
-			//cout << "word.substr(word.length()-1, 1)=" <<  word.substr(word.length()-1, 1)  << endl;
-			//exit(0);
-
-
-			//ifs >> word;   // zp
-			////cout << "word = " << word << endl;
-			//{
-			//     if(word == "z_valence=\"")
-			//     {
-			//        ifs >> word;
-			//        get_char(word);
-			//        this->zp = atoi(word.substr(0,Number[0]).c_str());
-			//     }
-			//     else
-			//     {
-			//        get_char(word);
-			//        this->zp = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-			//     }
-			//     //cout << "zp = " << this->zp << endl;
-			//}
-
-			//ifs >> word;   // total_psenergy
-			//{
-			//     if(word == "total_psenergy=\"")
-			//     {
-			//        ifs >> word;
-			//        get_char(word);
-			//        this->etotps = atof(word.substr(0,Number[0]).c_str());
-			//     }
-			//     else
-			//     {
-			//        get_char(word);
-			//        this->etotps = atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-			//     }
-			//     //cout << "etotps = " << this->etotps << endl;
-			//}
-			////cout << " word (total_psenergy) = " << word << endl;
-
-
-			//if(ONCVPSP == 0)    //zws modify 20160108
-			//{
-			//	READ_VALUE(ifs, word);   // wfc_cutoff
-			//	//cout << "word = " << word << endl;
-			//}
-			//READ_VALUE(ifs, word); // rho_cutoff
-			//cout << "word (cutoff) = " << word << endl;
-
-
-			//ifs >> word;             // lmax
-			////cout << "word (lmax) = " << word << endl;
-			//{
-			//        if(word == "l_max=\"")
-			//        {
-			//             ifs >> word;
-			//             get_char(word);
-			//             this->lmax = atoi(word.substr(0,Number[0]).c_str());
-			//        }
-			//        else
-			//        {
-			//             get_char(word);
-			//             this->lmax = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-			//        }
-
-			//}
-
-			////cout << "lmax = " << this->lmax << endl;
-
-			//if(ONCVPSP == 0)
-			//{
-			//   READ_VALUE(ifs, word);   // l_max_rho
-			//}
-
-			//READ_VALUE(ifs, word);   // l_local
-
-			//ifs >> word;   // mesh_size
-			////cout << "word (mesh) = " << word << endl;
-			//{
-			//     if(word == "mesh_size=\"")
-			//     {
-			//             ifs >> word;
-			//             get_char(word);
-			//             this->mesh = atoi(word.substr(0,Number[0]).c_str());
-			//     }
-			//     else
-			//     {
-			//             get_char(word);
-			//             this->mesh = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-			//     }
-			//     //cout << "mesh = " << this->mesh << endl;
-			//}
-
-
-
-			//ifs >> word;  // number_of_wfc
-			////cout << "word = " << word << endl;
-			//{
-			//     if(word == "number_of_wfc=\"")
-			//     {
-			//             ifs >> word;
-			//             get_char(word);
-			//             this->nwfc = atoi(word.substr(0,Number[0]).c_str());
-
-			//     }
-			//     else
-			//     {
-			//             get_char(word);
-			//             this->nwfc = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-			//     }
-			//     //cout << "nwfc = " << this->nwfc << endl;
-			//}
-			//     
-			//ifs >> word;   // number_of_proj
-			////cout << "word = " << word << endl;
-			//{
-			//     if(word == "number_of_proj=\"")
-			//     {
-			//             ifs >> word;
-			//             get_char(word);
-			//             this->nbeta = atoi(word.substr(0,Number[0]).c_str());
-
-			//     }
-			//     else
-			//     {
-			//             get_char(word);
-			//             this->nbeta = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-			//     }
-			//     //cout << "nbeta = " << this->nbeta << endl;
-			//}
-
-
-			// READ Mesh
-			if(ONCVPSP == 0)
-			{
-				SCAN_BEGIN(ifs, "<PP_MESH");
-			}
-			else
-			{
-				SCAN_BEGIN(ifs, "<PP_MESH>");
-			}
-
-			assert(mesh>0);
-			if(ONCVPSP == 0)
-			{
-				ifs >> word;             // dx
-				ifs >> word;             // mesh
-				ifs >> word;             // xmin
-				ifs >> word;             // rmax
-				ifs >> word;             // zmesh
-			}
-
-			SCAN_BEGIN(ifs, "<PP_R"); 
-			READ_VALUE(ifs, word);    // type  size  columns
-
-			//                        double  rmesh0 = 1;    //{zws add160108 delete160328
-			//                        int 	nmeshdel = 0;
-			//                        ifs >> rmesh0;
-			//                        if ( abs(rmesh0) < 1.0e-15 )
-			//                        {
-			//                            mesh       -= 1;
-			//                            nmeshdel   += 1;
-			//                        }
-			//                        cout << " mesh =" << mesh << endl;
-			//                    	if (mesh%2 == 0)
-			//                    	{
-			//                    	    mesh     -= 1;
-			//                    	    nmeshdel += 1;
-			//                    	}    //}zws add 20160108
-			//                    	cout << " nmeshdel =" << nmeshdel << endl;
-
-
-			delete[] r;
-			delete[] rab;
-			this->r = new double[mesh];
-			this->rab = new double[mesh];
-			ZEROS(r,mesh);
-			ZEROS(rab,mesh);
-
-
-			//                        if (nmeshdel == 0)    //{zws add160108 delete160328
-			//                        {
-			//                            this->r[0] = rmesh0;
-			//                            for (ir = 1;ir < mesh;ir++)
-			//                            {
-			//                                ifs >> this->r[ir];
-			//                            }
-			//                        }
-			//                        else
-			//                        {
-			//                            for ( int idel=0; idel < nmeshdel-1; idel++)
-			//                        	{
-			//                            	cout << "skip " << nmeshdel << "grid point(s) in PP mesh" << endl;
-			//                        	    double	tmpdel;
-			//                        	    ifs >> tmpdel;
-			//                        	}
-			//                            for (ir = 0;ir < mesh;ir++)
-			//                            {
-			//                                 ifs >> this->r[ir];
-			//                            }
-			//                        }    //}zws 20160108
-			for (ir = 0;ir < mesh;ir++)
-			{
-				ifs >> this->r[ir];
-			}
-			SCAN_END(ifs, "</PP_R>");
-
-			SCAN_BEGIN(ifs, "<PP_RAB");
-			READ_VALUE(ifs, word);    // type size columns
-
-			//                        for ( int idel=0; idel < nmeshdel; idel++)    //{zws add 20160108
-			//                    	{
-			//                    	    double	tmpdel;
-			//                    	    ifs >> tmpdel;
-			//                    	}    //}zws add 20160108
-			for (ir = 0;ir < mesh;ir++)
-			{
-				ifs >> this->rab[ir];
-			}
-			SCAN_END(ifs, "</PP_RAB>");
-			SCAN_END(ifs, "</PP_MESH>");
-
-			// READ NLCC
-			if (this->nlcc)
-			{
-				SCAN_BEGIN(ifs, "<PP_NLCC");
-				READ_VALUE(ifs, word);    // type size columns
-
-				assert(mesh>0);
-				delete[] rho_atc;
-				this->rho_atc = new double[mesh];
-				ZEROS(rho_atc, mesh);
-
-				for (ir = 0;ir < mesh;ir++)
-				{
-					ifs >> this->rho_atc[ir];
-				}
-				SCAN_END(ifs, "</PP_NLCC>");
-
-			}
-
-			// READ VLOCAL
-			SCAN_BEGIN(ifs, "<PP_LOCAL");
-			READ_VALUE(ifs, word);    // type size columns
-
-			assert(mesh>0);
-			delete[] vloc;
-			this->vloc = new double[mesh];
-			ZEROS(vloc, mesh);
-
-			for (ir = 0;ir < mesh;ir++)
-			{
-				ifs >> this->vloc[ir];
-			}
-
-			SCAN_END(ifs, "</PP_LOCAL>");
-
-			// READ NONLOCAL
-			SCAN_BEGIN(ifs, "<PP_NONLOCAL>");
-
-			delete[] kkbeta;
-			delete[] lll;
-			this->kkbeta = new int[nbeta];
-			this->lll = new int[nbeta];
-			this->beta.create(nbeta , mesh);
-			this->dion.create(nbeta , nbeta);
-
-			for(i=0;i<nbeta;i++)
-			{
-				ifs >> word;  //number
-				ifs >> word;  //type
-				if(word == "type=\"")
-				{
-					ifs >> word;
-				}
-				ifs >> word;  //size
-				if(word == "size=\"")
-				{
-					ifs >> word;
-				}
-
-				ifs >> word;  //columns
-				if(word == "columns=\"")
-				{
-					ifs >> word;
-				}
-
-				ifs >> word;  //index
-				{
-					if(word == "index=\"")
-					{
-						ifs >> word;
-						get_char(word);
-						//idum = atoi(word.substr(0,Number[0]).c_str());
-					}
-					else
-					{
-						get_char(word);
-						//idum = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-					}
-					//cout << "idum = " << idum << endl;
-				}
-
-				if(ONCVPSP == 0)
-				{
-					ifs >> word;  //label
-				}
-
-				ifs >> word;  //angular_momentum
-				if(word == "angular_momentum=\"")
-				{
-					ifs >> word;
-					get_char(word);
-					this->lll[i] = atoi(word.substr(0,Number[0]).c_str());
-
-				}
-				else
-				{
-					get_char(word);
-					this->lll[i] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-
-				ifs >> word;  //cutoff_radius_index
-				if(word == "cutoff_radius_index=\"")
-				{
-					ifs >> word;
-					get_char(word);
-					this->kkbeta[i] = atoi(word.substr(0,Number[0]).c_str());
-
-				}
-				else
-				{
-					get_char(word);
-					this->kkbeta[i] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-				//cout << "kkbeta[i] = " << this->kkbeta[i] << endl;
-
-				if(ONCVPSP ==0) 
-				{
-					ifs >> word;  //cutoff_radius
-					ifs >> word;  //ultrasoft_cutoff_radius
-				}
-				else
-				{
-					READ_VALUE(ifs, word); // cutoff_radius
-				}
-
-				for (ir=0;ir<mesh;ir++)
-				{
-					ifs >> this->beta(i, ir);
-
-				}
-
-				ifs >> word;  //number
-
-			}
-
-			// READ DIJ
-			SCAN_BEGIN(ifs, "<PP_DIJ");
-			READ_VALUE(ifs, word);  // type size columns
-
-			this->nd = nbeta * nbeta;
-			for(i=0;i<nbeta;i++)
-			{
-				for(j=0;j<nbeta;j++)
-				{
-					ifs >> dion(i,j);
-					if ( i != j  && dion(i,j) != 0.0 )
-					{
-						cout << " error: for i != j, Dij of Pseudopotential must be 0.0 " << endl;
-						exit(1);
-					}
-				}
-			}
-			SCAN_END(ifs, "</PP_DIJ>");
-
-			SCAN_END(ifs, "</PP_NONLOCAL>");
-
-			// READ PSWFC
-			SCAN_BEGIN(ifs, "<PP_PSWFC>");
-
-			delete[] els;
-			delete[] lchi;
-			delete[] oc;
-			this->els = new string[nwfc];
-			this->lchi = new int[nwfc];
-			this->oc = new double[nwfc];
-			ZEROS(lchi, nwfc); // angular momentum of each orbital
-			ZEROS(oc, nwfc);//occupation of each orbital
-
-			this->chi.create(this->nwfc, this->mesh);
-			for (i=0;i<nwfc;i++)
-			{
-				ifs >> word;  // number
-				ifs >> word;  // type
-				ifs >> word;  // size 
-				{
-					if(word == "size=\"")
-					{
-						ifs >> word;
-						word = "\"" + word ;
-					}
-				}
-				ifs >> word;  // columns
-				ifs >> word;  // index
-				ifs >> word;  // occupation
-				{
-					if(word == "occupation=\"")
-					{
-						ifs >> word;
-						word = "\"" + word ;
-					}
-					get_char(word);
-					oc[i] = atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				}
-
-				ifs >> word;  // pseudo_energy
-				if(word == "pseudo_energy=\"")
-				{
-					ifs >> word;
-					word = "\"" + word ;
-				}
-				get_char(word);
-
-				ifs >> word;  // label
-				if(word == "label=\"")
-				{
-					ifs >> word;
-					word = "\"" + word ;
-				}
-				get_char(word);
-				els[i] = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
-
-				ifs >> word;  // l
-				if(word == "l=\"")
-				{
-					ifs >> word;
-					word = "\"" + word ;
-				}
-				get_char(word);
-				lchi[i] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-				//cout << " lchi[i] = " << lchi[i] << endl;
-
-				ifs >> word; // >
-				if ( word !=  ">" )
-				{
-					cout << " error: bad end while reading CHI" << i <<  " of PSWFC" << endl;
-					exit(1);
-				}
-
-				for (ir = 0;ir < mesh;ir++)
-				{
-					ifs >> this->chi(i, ir);
-				}
-				ifs >> word;  // number
-			}
-
-			SCAN_END(ifs, "</PP_PSWFC>");
-
-			// READ RHOATOM
-			SCAN_BEGIN(ifs, "<PP_RHOATOM");
-			READ_VALUE(ifs, word); // type size columns
-
-			delete[] rho_at;
-			this->rho_at = new double[mesh];
-			ZEROS(rho_at, mesh);
-
-			for (ir = 0;ir < mesh;ir++)
-			{
-				ifs >> this->rho_at[ir];
-			}
-			SCAN_END(ifs, "</PP_RHOATOM>");
-
-			SCAN_BEGIN(ifs, "<PP_SPIN_ORB>");
-			//added by zhengdy-soc
-			delete[] this->jchi;
-			delete[] this->jjj;
-			delete[] this->nn;
-			this->jchi = new double [nwfc];
-			this->jjj = new double [nbeta];
-			this->nn = new int [nwfc];
-			ZEROS(jchi,nwfc);
-			ZEROS(jjj,nbeta);
-			ZEROS(nn,nwfc);
-
-			for(int round=0;round<2;round++)
-			{
-				ifs>>word;
-				if(word=="<PP_RELBETA.1")
-				{
-					for(int nb = 0;nb<nbeta;nb++)
-					{
-						if(nb!=0) ifs>>word; //RELBETA
-						ifs>>word;           //index
-						ifs>>word;           //lll
-						get_char(word);
-						this->lll[nb] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-						ifs>>word;           //jjj
-						get_char(word);
-						this->jjj[nb] = atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-					}
-				}
-				else if(word=="<PP_RELWFC.1")
-				{
-					if(round==0)
-					{
-						for(int nw = 0;nw<nwfc;nw++)
-						{
-							if(nw!=0) ifs>>word;     //RELWFC
-							ifs>>word;               //index
-							ifs>>word;               //els
-							get_char(word);
-							this->els[nw]= word.substr(Number[0]+1,(Number[1]-Number[0]-1));
-							ifs>>word;               //nn
-							get_char(word);
-							this->nn[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-							ifs>>word;               //lchi
-							get_char(word);
-							this->lchi[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-							ifs>>word;               //jchi
-							get_char(word);
-							this->jchi[nw]= atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-							ifs>>word;               //oc
-							get_char(word);
-							this->oc[nw]= atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-
-						}
-					}
-					else
-					{
-						for(int nw = 0;nw<nwfc;nw++)
-						{
-							if(nw!=0) ifs>>word;//RELWFC
-							ifs>>word;          //index
-							ifs>>word;          //lchi
-							get_char(word);
-							this->lchi[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-							ifs>>word;          //jchi
-							get_char(word);
-							this->jchi[nw]= atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-							ifs>>word;
-							get_char(word);
-							this->nn[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
-
-						}
-					}
-				}
-				else if(round==0)
-				{
-					this->has_so = 0;
-					//	cout<<"ignore SPIN_ORB part!"<<endl;
-					break;
-				}
-			}
-			SCAN_END(ifs, "</PP_SPIN_ORB>");
-
-			if (mesh%2 == 0)
-			{
-				mesh -= 1;
-			}
-
-			SCAN_END(ifs, "</UPF>");
-			break;
-		}
-	}
-	return 0;
-}
-
-
-void Pseudopot_upf::get_char( string ss)
-{
-    int i, q;
-    //char b[1]; //LiuXh 20171109
-    char b='\"'; //LiuXh 20171109
-    q =0;
-    //strcpy(b,"\""); //LiuXh 20171109
-
-    for(i=0;i<200;i++)
-    {
-        //if(ss[i]== b[0]) //LiuXh 20171109
-        if(ss[i]== b) //LiuXh 20171109
-        {
-           Number[q] = i;
-           q++;
-        }
-
-    }
-
-    return;
-}
-
 int Pseudopot_upf::average_p()
 {
 	int error = 0;
diff --git a/ABACUS.develop/source/src_pw/pseudopot_upf.h b/ABACUS.develop/source/src_pw/pseudopot_upf.h
index 819e13c118..87ad394661 100644
--- a/ABACUS.develop/source/src_pw/pseudopot_upf.h
+++ b/ABACUS.develop/source/src_pw/pseudopot_upf.h
@@ -59,16 +59,12 @@ class Pseudopot_upf
 	double *jjj;            // jjj(nbeta) j=l+1/2 or l-1/2 of beta
 
 	int  nd; 			// nl_5 // Number of nonzero Dij
-//	int  nqf;			// nl_6 // Number of expansion coefficients for q_{ij} (may be zero)
-//	int  nqlc;			// nl_7 // 2 * lmax  + 1
-//	matrix qqq;			// nl_9  // qqq(nbeta,nbeta):Q_{ij} = \int q_{ij}(r) dr
-//	realArray qfunc;	// nl_10 // qfunc(mesh,nbeta,nbeta):q_{ij}(r) for r > rinner(i)
-//	realArray qfcoef;	// nl_11 // qfcoef(nqf,0:2*lmax,nbeta,nbeta):expansion coefficients of q_{ij}(r) for r < rinner(i)
-						         //(not read if nqf=0)
 	
 	// the followings are for the vwr format
 	int spd_loc; 
-	int iTB_s, iTB_p, iTB_d;
+	int iTB_s;
+	int iTB_p;
+	int iTB_d;
 	double* vs;// local pseudopotential for s, unit is Hartree 
 	double* vp;// local pseudopotential for p 
 	double* vd;// local pseudopotential for d
@@ -92,7 +88,7 @@ class Pseudopot_upf
 
 	int read_pseudo_upf(ifstream &ifs);
 	int read_pseudo_vwr(ifstream &ifs);
-        int read_pseudo_upf201(ifstream &ifs);
+	int read_pseudo_upf201(ifstream &ifs);
 	void read_pseudo_header(ifstream &ifs);
 	void read_pseudo_mesh(ifstream &ifs);
 	void read_pseudo_nlcc(ifstream &ifs);
@@ -102,10 +98,10 @@ class Pseudopot_upf
 	void read_pseudo_rhoatom(ifstream &ifs);
 	void read_pseudo_addinfo(ifstream &ifs);
 	void read_pseudo_so(ifstream &ifs);
-        //string get_string( char ss[]);
-        //int get_int( char ss[]);
-        //double get_double( char ss[]);
-        void get_char( string ss);
+	//string get_string( char ss[]);
+	//int get_int( char ss[]);
+	//double get_double( char ss[]);
+	void get_char( string ss);
 };
 
 #endif //pseudopot_upf class
diff --git a/ABACUS.develop/source/src_pw/pseudopot_upf201.cpp b/ABACUS.develop/source/src_pw/pseudopot_upf201.cpp
new file mode 100644
index 0000000000..8474e60371
--- /dev/null
+++ b/ABACUS.develop/source/src_pw/pseudopot_upf201.cpp
@@ -0,0 +1,791 @@
+#include "pseudopot_upf.h"
+
+int Number[2];
+
+// added by zhangwenshuai
+int Pseudopot_upf::read_pseudo_upf201(ifstream &ifs)
+{
+    string dummy, word;
+    int i, j, ir, ONCVPSP;
+
+	while (ifs.good())
+	{
+		ifs >> dummy;
+		// We start from PP_Header
+		if(dummy=="<PP_HEADER")
+		{
+			// Read header
+			READ_VALUE(ifs, word);   // generated
+			READ_VALUE(ifs, word);   // author
+			READ_VALUE(ifs, word);   // date
+			READ_VALUE(ifs, word);   // comment
+
+			READ_VALUE(ifs, word);   // element
+			stringstream wdsstream(word);
+			getline(wdsstream,this->psd,'"'); 
+			getline(wdsstream,this->psd,'"'); 
+
+			ifs >> word;   // pseudo_type
+			if(word == "pseudo_type=\"")
+			{
+				ifs >> word;
+				get_char(word);
+				this->pp_type = word.substr(0,Number[0]);
+			}
+			else
+			{
+				get_char(word);
+				this->pp_type = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
+			}
+
+			if(pp_type!="NC") 
+			{
+				WARNING_QUIT("Pseudopot_upf::read_pseudo_header","unknown pseudo type");
+			}
+
+			READ_VALUE(ifs, word);   // relativistic
+			READ_VALUE(ifs, word);   // is_ultrasoft
+			if ( word.find("\"T\"") < word.length() ) // zws add 20160108
+			{
+				cout << "\n WARNING: ULTRASOFT PSEUDOPOTENTIAL IS NOT SUPPORTED !!! \n" << endl;
+			}
+			READ_VALUE(ifs, word);   // is_paw
+			if ( word.find("\"T\"") < word.length() )
+			{
+				cout << "\n WARNING: PAW PSEUDOPOTENTIAL IS NOT SUPPORTED !!! \n" << endl;
+			}
+
+			READ_VALUE(ifs, word);   // is_coulomb
+			ifs >> word;   // has_so
+			string so;
+
+			if(word == "has_so=\"")
+			{
+				ifs >> word;
+				get_char(word);
+				so = word.substr(0,Number[0]);
+			}
+			else
+			{
+				get_char(word);
+				so = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
+			}
+
+			if (so == "T")
+			{
+				this->has_so = true;
+			}
+			else
+			{
+				this->has_so = false;
+			}
+
+			READ_VALUE(ifs, word);   // has_wfc
+			READ_VALUE(ifs, word);   // has_gipaw
+
+			string nlc;
+			//char p[13] = "paw_as_gipaw";
+			ifs >> word;             // paw_as_gipaw?
+			//cout << "word.substr(0,30) = " << word.substr(0,30) << "."<< endl;
+			if( word.substr(0,13) == "paw_as_gipaw" )
+			{
+				ONCVPSP = 0;
+				ifs >> word;     // core_correction
+				if(word == "core_correction=\"")
+				{
+					ifs >> word;
+					get_char(word);
+					nlc = word.substr(0,Number[0]);
+				}
+				else
+				{
+					get_char(word);
+					nlc = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
+				}
+
+			}
+			else
+			{
+				ONCVPSP = 1; // Generated using ONCVPSP code by D. R. Hamann, SG15 DOJO
+				if(word == "core_correction=\"")
+				{
+					ifs >> word;
+					get_char(word);
+					nlc = word.substr(0,Number[0]);
+				}
+				else
+				{
+					get_char(word);
+					nlc = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
+				}
+			}
+
+			//cout << "nlc = " << nlc << endl;
+
+			if (nlc == "T")
+			{
+				this->nlcc = true;
+			}
+			else
+			{
+				this->nlcc = false;
+			}
+
+			READ_VALUE(ifs, word);   // functional
+			//cout << "word = " << word << endl;
+			//                        this->dft[0]="SLA";
+			//                        this->dft[1]="PZ";
+			//                        this->dft[2]="NOGX";
+			//                        this->dft[3]="NOGC";
+
+			string funcstr;  //{zws 01-06-16
+			wdsstream.str("");
+			wdsstream.clear();
+			wdsstream << word;
+			for ( int idft = 0; idft < 2; idft++)
+			{
+				getline(wdsstream,funcstr,'"');
+			}
+			wdsstream.str("");
+			wdsstream.clear();
+			wdsstream << funcstr;
+
+			for( int idft = 0; idft < 4; idft++ )
+			{
+				getline(wdsstream,dft[idft],'-');
+			}
+
+			do 
+			{
+				getline(ifs, word);
+				//cout << "word       = " << word << endl;
+				word.erase(0,word.find_first_not_of(" ") );
+				word.erase(word.find_last_not_of(" ")+1 );
+				//word = trim(word);
+				//cout << "trim(word) = " << word << endl;
+				get_char(word);
+				//cout << " Number = " << Number[0] << ", " << Number[1] << endl;
+				//cout << word.substr(0,Number[0])  << "__" << word.substr(Number[0]+1, Number[1]-Number[0]-1) << endl;
+				dummy = word.substr(0,Number[0]) ;
+				//cout << " dummy = " << dummy << endl;
+				if( dummy == "z_valence=" ) 
+				{
+					this->zp = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+				else if ( dummy == "total_psenergy=" ) 
+				{
+					this->etotps = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+				else if ( dummy == "rho_cutoff=" )
+				{
+				}
+				else if ( dummy == "wfc_cutoff=" ) 
+				{
+				}
+				else if ( dummy == "l_max=" ) 
+				{
+					this->lmax = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+				else if ( dummy == "mesh_size=" ) 
+				{
+					this->mesh = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+				else if ( dummy == "number_of_wfc=" ) 
+				{
+					this->nwfc = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+				else if ( dummy == "number_of_proj=" ) 
+				{
+					this->nbeta = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+				//break;
+
+			}while( word.substr(word.length()-1, 1) !=">" ); 
+			//cout << "word.substr(word.length()-1, 1)=" <<  word.substr(word.length()-1, 1)  << endl;
+			//exit(0);
+
+
+			//ifs >> word;   // zp
+			////cout << "word = " << word << endl;
+			//{
+			//     if(word == "z_valence=\"")
+			//     {
+			//        ifs >> word;
+			//        get_char(word);
+			//        this->zp = atoi(word.substr(0,Number[0]).c_str());
+			//     }
+			//     else
+			//     {
+			//        get_char(word);
+			//        this->zp = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+			//     }
+			//     //cout << "zp = " << this->zp << endl;
+			//}
+
+			//ifs >> word;   // total_psenergy
+			//{
+			//     if(word == "total_psenergy=\"")
+			//     {
+			//        ifs >> word;
+			//        get_char(word);
+			//        this->etotps = atof(word.substr(0,Number[0]).c_str());
+			//     }
+			//     else
+			//     {
+			//        get_char(word);
+			//        this->etotps = atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+			//     }
+			//     //cout << "etotps = " << this->etotps << endl;
+			//}
+			////cout << " word (total_psenergy) = " << word << endl;
+
+
+			//if(ONCVPSP == 0)    //zws modify 20160108
+			//{
+			//	READ_VALUE(ifs, word);   // wfc_cutoff
+			//	//cout << "word = " << word << endl;
+			//}
+			//READ_VALUE(ifs, word); // rho_cutoff
+			//cout << "word (cutoff) = " << word << endl;
+
+
+			//ifs >> word;             // lmax
+			////cout << "word (lmax) = " << word << endl;
+			//{
+			//        if(word == "l_max=\"")
+			//        {
+			//             ifs >> word;
+			//             get_char(word);
+			//             this->lmax = atoi(word.substr(0,Number[0]).c_str());
+			//        }
+			//        else
+			//        {
+			//             get_char(word);
+			//             this->lmax = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+			//        }
+
+			//}
+
+			////cout << "lmax = " << this->lmax << endl;
+
+			//if(ONCVPSP == 0)
+			//{
+			//   READ_VALUE(ifs, word);   // l_max_rho
+			//}
+
+			//READ_VALUE(ifs, word);   // l_local
+
+			//ifs >> word;   // mesh_size
+			////cout << "word (mesh) = " << word << endl;
+			//{
+			//     if(word == "mesh_size=\"")
+			//     {
+			//             ifs >> word;
+			//             get_char(word);
+			//             this->mesh = atoi(word.substr(0,Number[0]).c_str());
+			//     }
+			//     else
+			//     {
+			//             get_char(word);
+			//             this->mesh = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+			//     }
+			//     //cout << "mesh = " << this->mesh << endl;
+			//}
+
+
+
+			//ifs >> word;  // number_of_wfc
+			////cout << "word = " << word << endl;
+			//{
+			//     if(word == "number_of_wfc=\"")
+			//     {
+			//             ifs >> word;
+			//             get_char(word);
+			//             this->nwfc = atoi(word.substr(0,Number[0]).c_str());
+
+			//     }
+			//     else
+			//     {
+			//             get_char(word);
+			//             this->nwfc = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+			//     }
+			//     //cout << "nwfc = " << this->nwfc << endl;
+			//}
+			//     
+			//ifs >> word;   // number_of_proj
+			////cout << "word = " << word << endl;
+			//{
+			//     if(word == "number_of_proj=\"")
+			//     {
+			//             ifs >> word;
+			//             get_char(word);
+			//             this->nbeta = atoi(word.substr(0,Number[0]).c_str());
+
+			//     }
+			//     else
+			//     {
+			//             get_char(word);
+			//             this->nbeta = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+			//     }
+			//     //cout << "nbeta = " << this->nbeta << endl;
+			//}
+
+
+			// READ Mesh
+			if(ONCVPSP == 0)
+			{
+				SCAN_BEGIN(ifs, "<PP_MESH");
+			}
+			else
+			{
+				SCAN_BEGIN(ifs, "<PP_MESH>");
+			}
+
+			assert(mesh>0);
+			if(ONCVPSP == 0)
+			{
+				ifs >> word;             // dx
+				ifs >> word;             // mesh
+				ifs >> word;             // xmin
+				ifs >> word;             // rmax
+				ifs >> word;             // zmesh
+			}
+
+			SCAN_BEGIN(ifs, "<PP_R"); 
+			READ_VALUE(ifs, word);    // type  size  columns
+
+			//                        double  rmesh0 = 1;    //{zws add160108 delete160328
+			//                        int 	nmeshdel = 0;
+			//                        ifs >> rmesh0;
+			//                        if ( abs(rmesh0) < 1.0e-15 )
+			//                        {
+			//                            mesh       -= 1;
+			//                            nmeshdel   += 1;
+			//                        }
+			//                        cout << " mesh =" << mesh << endl;
+			//                    	if (mesh%2 == 0)
+			//                    	{
+			//                    	    mesh     -= 1;
+			//                    	    nmeshdel += 1;
+			//                    	}    //}zws add 20160108
+			//                    	cout << " nmeshdel =" << nmeshdel << endl;
+
+
+			delete[] r;
+			delete[] rab;
+			this->r = new double[mesh];
+			this->rab = new double[mesh];
+			ZEROS(r,mesh);
+			ZEROS(rab,mesh);
+
+
+			//                        if (nmeshdel == 0)    //{zws add160108 delete160328
+			//                        {
+			//                            this->r[0] = rmesh0;
+			//                            for (ir = 1;ir < mesh;ir++)
+			//                            {
+			//                                ifs >> this->r[ir];
+			//                            }
+			//                        }
+			//                        else
+			//                        {
+			//                            for ( int idel=0; idel < nmeshdel-1; idel++)
+			//                        	{
+			//                            	cout << "skip " << nmeshdel << "grid point(s) in PP mesh" << endl;
+			//                        	    double	tmpdel;
+			//                        	    ifs >> tmpdel;
+			//                        	}
+			//                            for (ir = 0;ir < mesh;ir++)
+			//                            {
+			//                                 ifs >> this->r[ir];
+			//                            }
+			//                        }    //}zws 20160108
+			for (ir = 0;ir < mesh;ir++)
+			{
+				ifs >> this->r[ir];
+			}
+			SCAN_END(ifs, "</PP_R>");
+
+			SCAN_BEGIN(ifs, "<PP_RAB");
+			READ_VALUE(ifs, word);    // type size columns
+
+			//                        for ( int idel=0; idel < nmeshdel; idel++)    //{zws add 20160108
+			//                    	{
+			//                    	    double	tmpdel;
+			//                    	    ifs >> tmpdel;
+			//                    	}    //}zws add 20160108
+			for (ir = 0;ir < mesh;ir++)
+			{
+				ifs >> this->rab[ir];
+			}
+			SCAN_END(ifs, "</PP_RAB>");
+			SCAN_END(ifs, "</PP_MESH>");
+
+			// READ NLCC
+			if (this->nlcc)
+			{
+				SCAN_BEGIN(ifs, "<PP_NLCC");
+				READ_VALUE(ifs, word);    // type size columns
+
+				assert(mesh>0);
+				delete[] rho_atc;
+				this->rho_atc = new double[mesh];
+				ZEROS(rho_atc, mesh);
+
+				for (ir = 0;ir < mesh;ir++)
+				{
+					ifs >> this->rho_atc[ir];
+				}
+				SCAN_END(ifs, "</PP_NLCC>");
+
+			}
+
+			// READ VLOCAL
+			SCAN_BEGIN(ifs, "<PP_LOCAL");
+			READ_VALUE(ifs, word);    // type size columns
+
+			assert(mesh>0);
+			delete[] vloc;
+			this->vloc = new double[mesh];
+			ZEROS(vloc, mesh);
+
+			for (ir = 0;ir < mesh;ir++)
+			{
+				ifs >> this->vloc[ir];
+			}
+
+			SCAN_END(ifs, "</PP_LOCAL>");
+
+			// READ NONLOCAL
+			SCAN_BEGIN(ifs, "<PP_NONLOCAL>");
+
+			delete[] kkbeta;
+			delete[] lll;
+			this->kkbeta = new int[nbeta];
+			this->lll = new int[nbeta];
+			this->beta.create(nbeta , mesh);
+			this->dion.create(nbeta , nbeta);
+
+			for(i=0;i<nbeta;i++)
+			{
+				ifs >> word;  //number
+				ifs >> word;  //type
+				if(word == "type=\"")
+				{
+					ifs >> word;
+				}
+				ifs >> word;  //size
+				if(word == "size=\"")
+				{
+					ifs >> word;
+				}
+
+				ifs >> word;  //columns
+				if(word == "columns=\"")
+				{
+					ifs >> word;
+				}
+
+				ifs >> word;  //index
+				{
+					if(word == "index=\"")
+					{
+						ifs >> word;
+						get_char(word);
+						//idum = atoi(word.substr(0,Number[0]).c_str());
+					}
+					else
+					{
+						get_char(word);
+						//idum = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+					}
+					//cout << "idum = " << idum << endl;
+				}
+
+				if(ONCVPSP == 0)
+				{
+					ifs >> word;  //label
+				}
+
+				ifs >> word;  //angular_momentum
+				if(word == "angular_momentum=\"")
+				{
+					ifs >> word;
+					get_char(word);
+					this->lll[i] = atoi(word.substr(0,Number[0]).c_str());
+
+				}
+				else
+				{
+					get_char(word);
+					this->lll[i] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+
+				ifs >> word;  //cutoff_radius_index
+				if(word == "cutoff_radius_index=\"")
+				{
+					ifs >> word;
+					get_char(word);
+					this->kkbeta[i] = atoi(word.substr(0,Number[0]).c_str());
+
+				}
+				else
+				{
+					get_char(word);
+					this->kkbeta[i] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+				//cout << "kkbeta[i] = " << this->kkbeta[i] << endl;
+
+				if(ONCVPSP ==0) 
+				{
+					ifs >> word;  //cutoff_radius
+					ifs >> word;  //ultrasoft_cutoff_radius
+				}
+				else
+				{
+					READ_VALUE(ifs, word); // cutoff_radius
+				}
+
+				for (ir=0;ir<mesh;ir++)
+				{
+					ifs >> this->beta(i, ir);
+
+				}
+
+				ifs >> word;  //number
+
+			}
+
+			// READ DIJ
+			SCAN_BEGIN(ifs, "<PP_DIJ");
+			READ_VALUE(ifs, word);  // type size columns
+
+			this->nd = nbeta * nbeta;
+			for(i=0;i<nbeta;i++)
+			{
+				for(j=0;j<nbeta;j++)
+				{
+					ifs >> dion(i,j);
+					if ( i != j  && dion(i,j) != 0.0 )
+					{
+						cout << " error: for i != j, Dij of Pseudopotential must be 0.0 " << endl;
+						exit(1);
+					}
+				}
+			}
+			SCAN_END(ifs, "</PP_DIJ>");
+
+			SCAN_END(ifs, "</PP_NONLOCAL>");
+
+			// READ PSWFC
+			SCAN_BEGIN(ifs, "<PP_PSWFC>");
+
+			delete[] els;
+			delete[] lchi;
+			delete[] oc;
+			this->els = new string[nwfc];
+			this->lchi = new int[nwfc];
+			this->oc = new double[nwfc];
+			ZEROS(lchi, nwfc); // angular momentum of each orbital
+			ZEROS(oc, nwfc);//occupation of each orbital
+
+			this->chi.create(this->nwfc, this->mesh);
+			for (i=0;i<nwfc;i++)
+			{
+				ifs >> word;  // number
+				ifs >> word;  // type
+				ifs >> word;  // size 
+				{
+					if(word == "size=\"")
+					{
+						ifs >> word;
+						word = "\"" + word ;
+					}
+				}
+				ifs >> word;  // columns
+				ifs >> word;  // index
+				ifs >> word;  // occupation
+				{
+					if(word == "occupation=\"")
+					{
+						ifs >> word;
+						word = "\"" + word ;
+					}
+					get_char(word);
+					oc[i] = atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				}
+
+				ifs >> word;  // pseudo_energy
+				if(word == "pseudo_energy=\"")
+				{
+					ifs >> word;
+					word = "\"" + word ;
+				}
+				get_char(word);
+
+				ifs >> word;  // label
+				if(word == "label=\"")
+				{
+					ifs >> word;
+					word = "\"" + word ;
+				}
+				get_char(word);
+				els[i] = word.substr(Number[0]+1,(Number[1]-Number[0]-1));
+
+				ifs >> word;  // l
+				if(word == "l=\"")
+				{
+					ifs >> word;
+					word = "\"" + word ;
+				}
+				get_char(word);
+				lchi[i] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+				//cout << " lchi[i] = " << lchi[i] << endl;
+
+				ifs >> word; // >
+				if ( word !=  ">" )
+				{
+					cout << " error: bad end while reading CHI" << i <<  " of PSWFC" << endl;
+					exit(1);
+				}
+
+				for (ir = 0;ir < mesh;ir++)
+				{
+					ifs >> this->chi(i, ir);
+				}
+				ifs >> word;  // number
+			}
+
+			SCAN_END(ifs, "</PP_PSWFC>");
+
+			// READ RHOATOM
+			SCAN_BEGIN(ifs, "<PP_RHOATOM");
+			READ_VALUE(ifs, word); // type size columns
+
+			delete[] rho_at;
+			this->rho_at = new double[mesh];
+			ZEROS(rho_at, mesh);
+
+			for (ir = 0;ir < mesh;ir++)
+			{
+				ifs >> this->rho_at[ir];
+			}
+			SCAN_END(ifs, "</PP_RHOATOM>");
+
+			SCAN_BEGIN(ifs, "<PP_SPIN_ORB>");
+			//added by zhengdy-soc
+			delete[] this->jchi;
+			delete[] this->jjj;
+			delete[] this->nn;
+			this->jchi = new double [nwfc];
+			this->jjj = new double [nbeta];
+			this->nn = new int [nwfc];
+			ZEROS(jchi,nwfc);
+			ZEROS(jjj,nbeta);
+			ZEROS(nn,nwfc);
+
+			for(int round=0;round<2;round++)
+			{
+				ifs>>word;
+				if(word=="<PP_RELBETA.1")
+				{
+					for(int nb = 0;nb<nbeta;nb++)
+					{
+						if(nb!=0) ifs>>word; //RELBETA
+						ifs>>word;           //index
+						ifs>>word;           //lll
+						get_char(word);
+						this->lll[nb] = atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+						ifs>>word;           //jjj
+						get_char(word);
+						this->jjj[nb] = atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+					}
+				}
+				else if(word=="<PP_RELWFC.1")
+				{
+					if(round==0)
+					{
+						for(int nw = 0;nw<nwfc;nw++)
+						{
+							if(nw!=0) ifs>>word;     //RELWFC
+							ifs>>word;               //index
+							ifs>>word;               //els
+							get_char(word);
+							this->els[nw]= word.substr(Number[0]+1,(Number[1]-Number[0]-1));
+							ifs>>word;               //nn
+							get_char(word);
+							this->nn[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+							ifs>>word;               //lchi
+							get_char(word);
+							this->lchi[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+							ifs>>word;               //jchi
+							get_char(word);
+							this->jchi[nw]= atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+							ifs>>word;               //oc
+							get_char(word);
+							this->oc[nw]= atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+
+						}
+					}
+					else
+					{
+						for(int nw = 0;nw<nwfc;nw++)
+						{
+							if(nw!=0) ifs>>word;//RELWFC
+							ifs>>word;          //index
+							ifs>>word;          //lchi
+							get_char(word);
+							this->lchi[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+							ifs>>word;          //jchi
+							get_char(word);
+							this->jchi[nw]= atof(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+							ifs>>word;
+							get_char(word);
+							this->nn[nw]= atoi(word.substr(Number[0]+1,(Number[1]-Number[0]-1)).c_str());
+
+						}
+					}
+				}
+				else if(round==0)
+				{
+					this->has_so = 0;
+					//	cout<<"ignore SPIN_ORB part!"<<endl;
+					break;
+				}
+			}
+			SCAN_END(ifs, "</PP_SPIN_ORB>");
+
+			if (mesh%2 == 0)
+			{
+				mesh -= 1;
+			}
+
+			SCAN_END(ifs, "</UPF>");
+			break;
+		}
+	}
+	return 0;
+}
+
+
+void Pseudopot_upf::get_char( string ss)
+{
+    int i, q;
+    //char b[1]; //LiuXh 20171109
+    char b='\"'; //LiuXh 20171109
+    q =0;
+    //strcpy(b,"\""); //LiuXh 20171109
+
+    for(i=0;i<200;i++)
+    {
+        //if(ss[i]== b[0]) //LiuXh 20171109
+        if(ss[i]== b) //LiuXh 20171109
+        {
+           Number[q] = i;
+           q++;
+        }
+
+    }
+
+    return;
+}
diff --git a/ABACUS.develop/source/src_pw/unitcell.cpp b/ABACUS.develop/source/src_pw/unitcell.cpp
index 5a47691423..518cc9dbab 100644
--- a/ABACUS.develop/source/src_pw/unitcell.cpp
+++ b/ABACUS.develop/source/src_pw/unitcell.cpp
@@ -1,7 +1,3 @@
-//==========================================================
-// AUTHOR : Lixin He , mohan
-// DATE : 2008-11-08
-//==========================================================
 #include <cstdlib>
 #ifdef __MPI
 #include <mpi.h>
diff --git a/ABACUS.develop/source/src_pw/unitcell_pseudo.cpp b/ABACUS.develop/source/src_pw/unitcell_pseudo.cpp
index 70a9a2188f..e5eb560552 100644
--- a/ABACUS.develop/source/src_pw/unitcell_pseudo.cpp
+++ b/ABACUS.develop/source/src_pw/unitcell_pseudo.cpp
@@ -177,7 +177,7 @@ void UnitCell_pseudo::setup_cell(
 	ofs_running << "\n\n\n\n";
 
 
-	this->read_pseudopot(s_pseudopot_dir);
+	this->read_cell_pseudopots(s_pseudopot_dir);
 	
 	if(MY_RANK == 0) 
 	{
diff --git a/ABACUS.develop/source/src_pw/unitcell_pseudo.h b/ABACUS.develop/source/src_pw/unitcell_pseudo.h
index 0c46ca679d..19b5585c59 100644
--- a/ABACUS.develop/source/src_pw/unitcell_pseudo.h
+++ b/ABACUS.develop/source/src_pw/unitcell_pseudo.h
@@ -36,11 +36,10 @@ class UnitCell_pseudo : public UnitCell
 	void check_dtau(void);
     void setup_cell_after_vc(const string &s_pseudopot_dir, const string &fn, ofstream &log); //LiuXh add 20180515
 
-private: // member variables
 	bool set_atom_flag;//added on 2009-3-8 by mohan
 
-private: // member functions
-	void read_pseudopot(const string &fn); // read in pseudopotential from files for each type of atom
+	// read in pseudopotential from files for each type of atom
+	void read_cell_pseudopots(const string &fn);
 
 	//================================================================
 	// cal_natomwfc : calculate total number of atomic wavefunctions
diff --git a/ABACUS.develop/source/src_pw/xc_type.cpp b/ABACUS.develop/source/src_pw/xc_type.cpp
index 0299cc1877..46d2463162 100644
--- a/ABACUS.develop/source/src_pw/xc_type.cpp
+++ b/ABACUS.develop/source/src_pw/xc_type.cpp
@@ -1,7 +1,3 @@
-//==========================================================
-// AUTHOR : Lixin He ,mohan
-// DATE : 2008-11-08
-//==========================================================
 #include "xc_type.h"
 #include "src_global/global_function.h"
 #include "src_pw/global.h"