diff --git a/source/module_cell/atom_pseudo.cpp b/source/module_cell/atom_pseudo.cpp
index 863cbbbb3a..7f15bce71e 100644
--- a/source/module_cell/atom_pseudo.cpp
+++ b/source/module_cell/atom_pseudo.cpp
@@ -2,25 +2,10 @@
 
 Atom_pseudo::Atom_pseudo()
 {
-    for (int is = 0; is < 4; is++) {
-        this->index1_soc[is] = nullptr;
-}
-    for (int is = 0; is < 4; is++) {
-        this->index2_soc[is] = nullptr;
-}
 }
 
 Atom_pseudo::~Atom_pseudo()
 {
-    for (int is = 0; is < 4; is++)
-    {
-        if (this->index1_soc[is] != nullptr) {
-            delete[] this->index1_soc[is];
-}
-        if (this->index2_soc[is] != nullptr) {
-            delete[] this->index2_soc[is];
-}
-    }
 }
 
 // mohan add 2021-05-07
@@ -42,10 +27,8 @@ void Atom_pseudo::set_d_so(ModuleBase::ComplexMatrix& d_so_in,
     for (int is = 0; is < spin_dimension; is++)
     {
         this->non_zero_count_soc[is] = 0;
-        delete[] this->index1_soc[is];
-        this->index1_soc[is] = new int[nproj_soc * nproj_soc];
-        delete[] this->index2_soc[is];
-        this->index2_soc[is] = new int[nproj_soc * nproj_soc];
+        this->index1_soc[is] = std::vector<int>(nproj_soc * nproj_soc, 0);
+        this->index2_soc[is] = std::vector<int>(nproj_soc * nproj_soc, 0);
     }
 
     if (!has_so)
@@ -189,26 +172,20 @@ void Atom_pseudo::bcast_atom_pseudo()
 
     if (GlobalV::MY_RANK != 0)
     {
-        delete[] jjj;
-        delete[] els;
-        delete[] lchi;
-        delete[] oc;
-        delete[] jchi;
-        delete[] nn;
-        jjj = new double[nbeta];
-        els = new std::string[nchi];
-        lchi = new int[nchi];
-        oc = new double[nchi];
-        jchi = new double[nchi];
-        nn = new int[nchi];
+        jjj = std::vector<double>(nbeta, 0.0);
+        els = std::vector<std::string>(nchi, "");
+        lchi = std::vector<int>(nchi, 0);
+        oc = std::vector<double>(nchi, 0.0);
+        jchi = std::vector<double>(nchi, 0.0);
+        nn = std::vector<int>(nchi, 0);
     }
 
-    Parallel_Common::bcast_double(jjj, nbeta);
-    Parallel_Common::bcast_string(els, nchi);
-    Parallel_Common::bcast_int(lchi, nchi);
-    Parallel_Common::bcast_double(oc, nchi);
-    Parallel_Common::bcast_double(jchi, nchi);
-    Parallel_Common::bcast_int(nn, nchi);
+    Parallel_Common::bcast_double(jjj.data(), nbeta);
+    Parallel_Common::bcast_string(els.data(), nchi);
+    Parallel_Common::bcast_int(lchi.data(), nchi);
+    Parallel_Common::bcast_double(oc.data(), nchi);
+    Parallel_Common::bcast_double(jchi.data(), nchi);
+    Parallel_Common::bcast_int(nn.data(), nchi);
     // == end of pseudo_h
 
     // == pseudo_atom ==
@@ -217,31 +194,26 @@ void Atom_pseudo::bcast_atom_pseudo()
     if (GlobalV::MY_RANK != 0)
     {
         assert(mesh != 0);
-        delete[] r;
-        delete[] rab;
-        delete[] rho_atc;
-        delete[] rho_at;
-        r = new double[mesh];
-        rab = new double[mesh];
-        rho_atc = new double[mesh];
-        rho_at = new double[mesh];
+        r = std::vector<double>(mesh, 0.0);
+        rab = std::vector<double>(mesh, 0.0);
+        rho_atc = std::vector<double>(mesh, 0.0);
+        rho_at = std::vector<double>(mesh, 0.0);
         chi.create(nchi, mesh);
     }
 
-    Parallel_Common::bcast_double(r, mesh);
-    Parallel_Common::bcast_double(rab, mesh);
-    Parallel_Common::bcast_double(rho_atc, mesh);
-    Parallel_Common::bcast_double(rho_at, mesh);
+    Parallel_Common::bcast_double(r.data(), mesh);
+    Parallel_Common::bcast_double(rab.data(), mesh);
+    Parallel_Common::bcast_double(rho_atc.data(), mesh);
+    Parallel_Common::bcast_double(rho_at.data(), mesh);
     Parallel_Common::bcast_double(chi.c, nchi * mesh);
     // == end of pseudo_atom ==
 
     // == pseudo_vl ==
     if (GlobalV::MY_RANK != 0)
     {
-        delete[] vloc_at;
-        vloc_at = new double[mesh];
+        vloc_at = std::vector<double>(mesh, 0.0);
     }
-    Parallel_Common::bcast_double(vloc_at, mesh);
+    Parallel_Common::bcast_double(vloc_at.data(), mesh);
     // == end of pseudo_vl ==
 
     // == pseudo ==
@@ -251,10 +223,9 @@ void Atom_pseudo::bcast_atom_pseudo()
 
     if (GlobalV::MY_RANK != 0)
     {
-        delete[] lll;
-        lll = new int[nbeta];
+        lll = std::vector<int>(nbeta, 0);
     }
-    Parallel_Common::bcast_int(lll, nbeta);
+    Parallel_Common::bcast_int(lll.data(), nbeta);
     Parallel_Common::bcast_int(kkbeta);
     Parallel_Common::bcast_int(nh);
 
diff --git a/source/module_cell/atom_pseudo.h b/source/module_cell/atom_pseudo.h
index def72ac6d9..c02687c78f 100644
--- a/source/module_cell/atom_pseudo.h
+++ b/source/module_cell/atom_pseudo.h
@@ -1,12 +1,11 @@
 #ifndef ATOM_PSEUDO_H
 #define ATOM_PSEUDO_H
 
-#include "../module_base/global_function.h"
-#include "../module_base/global_variable.h"
-#include "../module_base/vector3.h"
-#include "../module_io/output.h"
-#include "../module_base/complexarray.h"
-#include "../module_base/complexmatrix.h"
+#include "module_base/global_variable.h"
+#include "module_base/vector3.h"
+#include "module_io/output.h"
+#include "module_base/complexarray.h"
+#include "module_base/complexmatrix.h"
 #include "pseudo.h"
 
 
@@ -22,9 +21,9 @@ class Atom_pseudo : public pseudo
 	ModuleBase::matrix d_real; //(:,:), non-spin-orbit case
 	int nproj;
 	int nproj_soc; // dimension of D_ij^so
-	int non_zero_count_soc[4];
-	int *index1_soc[4];
-	int *index2_soc[4];
+	std::vector<int> non_zero_count_soc = {0, 0, 0, 0};
+	std::vector<std::vector<int>> index1_soc = {{}, {}, {}, {}};
+	std::vector<std::vector<int>> index2_soc = {{}, {}, {}, {}};
 
 	void set_d_so( // mohan add 2021-05-07
 		ModuleBase::ComplexMatrix &d_so_in,
diff --git a/source/module_cell/pseudo.cpp b/source/module_cell/pseudo.cpp
index b1793628a5..7743e45447 100644
--- a/source/module_cell/pseudo.cpp
+++ b/source/module_cell/pseudo.cpp
@@ -7,18 +7,6 @@ pseudo::pseudo()
 
 pseudo::~pseudo()
 {
-    delete[] els;
-    delete[] lchi;
-    delete[] oc;
-    delete[] jjj;
-    delete[] jchi;
-    delete[] nn;
-    delete[] vloc_at;
-    delete[] r;
-    delete[] rab;
-    delete[] rho_at;
-    delete[] rho_atc;
-    delete[] lll;
 }
 
 void pseudo::print_pseudo(std::ofstream& ofs)
@@ -27,7 +15,7 @@ void pseudo::print_pseudo(std::ofstream& ofs)
 	ofs << "\n pseudo : ";
 	ofs << "\n kkbeta	" << kkbeta;
 	ofs << "\n nh  " << nh;
-	output::printr1_d(ofs, " lll : ", lll, nbeta);
+	output::printr1_d(ofs, " lll : ", lll.data(), nbeta);
 	output::printrm(ofs, " betar : ", betar);
 	output::printrm(ofs, " dion : ", dion);
 	ofs << "\n ----------------------";
@@ -39,11 +27,11 @@ void pseudo::print_pseudo_atom(std::ofstream &ofs)
 	ofs << "\n pseudo_atom : ";
 	ofs << "\n msh	" << msh;
 //	ofs	<< "\n nchi	" << nchi;
-	output::printr1_d(ofs, " r : ", r, mesh);
-	output::printr1_d(ofs, " rab : ", rab, mesh);
-	output::printr1_d(ofs, " rho_atc : ", rho_atc, mesh);
-	output::printr1_d(ofs, " rho_at : ", rho_at, mesh);
-	output::printr1_d(ofs," jchi : ", jchi, nchi);
+	output::printr1_d(ofs, " r : ", r.data(), mesh);
+	output::printr1_d(ofs, " rab : ", rab.data(), mesh);
+	output::printr1_d(ofs, " rho_atc : ", rho_atc.data(), mesh);
+	output::printr1_d(ofs, " rho_at : ", rho_at.data(), mesh);
+	output::printr1_d(ofs," jchi : ", jchi.data(), nchi);
 	output::printrm(ofs, " chi : ", chi);
 	ofs << "\n ----------------------";
 }
@@ -53,7 +41,7 @@ void pseudo::print_pseudo_vl(std::ofstream &ofs)
 {
 	ofs << "\n pseudo_vl:";
 	print_pseudo_atom(ofs);
-	output::printr1_d(ofs, "vloc_at : ", vloc_at, mesh);
+	output::printr1_d(ofs, "vloc_at : ", vloc_at.data(), mesh);
 	ofs << "\n ----------------------------------- ";
 }
 
@@ -75,8 +63,8 @@ void pseudo::print_pseudo_h(std::ofstream &ofs)
     ofs << "\n nchi " << nchi;
     ofs << "\n nbeta    " << nbeta;
 //  out.printr1_d(ofs," els: ", els, nchi);
-    output::printr1_d(ofs, " lchi: ", lchi, nchi);
-    output::printr1_d(ofs, " oc: ", oc, nchi);
+    output::printr1_d(ofs, " lchi: ", lchi.data(), nchi);
+    output::printr1_d(ofs, " oc: ", oc.data(), nchi);
     ofs << "\n ----------------------";
 }
 
diff --git a/source/module_cell/pseudo.h b/source/module_cell/pseudo.h
index 8bebf278da..22af587557 100644
--- a/source/module_cell/pseudo.h
+++ b/source/module_cell/pseudo.h
@@ -33,26 +33,26 @@ class pseudo
     int nqlc = 0;         // number of angular momenta in Q
     int kkbeta = 0;       // kkbeta, point where the beta are zero
 
-    std::string* els = nullptr; // els[nchi]
-    int* lchi = nullptr;        // lchi[nchi]
-    double* oc = nullptr;       // oc[nchi]
+    std::vector<std::string> els = {}; // els[nchi]
+    std::vector<int> lchi = {};        // lchi[nchi]
+    std::vector<double> oc = {};       // oc[nchi]
 
-    double* jjj = nullptr;  // total angual momentum, jjj[nbeta]
-    double* jchi = nullptr; // jchi(nwfc), added by zhengdy-soc
-    int* nn = nullptr;
+    std::vector<double> jjj = {};  // total angual momentum, jjj[nbeta]
+    std::vector<double> jchi = {}; // jchi(nwfc), added by zhengdy-soc
+    std::vector<int> nn = {};
 
     // Local pseudopotentials
-    double* vloc_at = nullptr; // [mesh], local potential( = pseudopot_upf.vloc )
+    std::vector<double> vloc_at = {}; // [mesh], local potential( = pseudopot_upf.vloc )
 
     // <PP_MESH>
-    double* r = nullptr;   // radial logaritmic mesh, r[0:mesh-1]
-    double* rab = nullptr; // derivative of the radial mesh, rab[0:mesh-1]
+    std::vector<double> r = {};   // radial logaritmic mesh, r[0:mesh-1]
+    std::vector<double> rab = {}; // derivative of the radial mesh, rab[0:mesh-1]
 
     //<PP_NLCC>
-    double* rho_atc = nullptr; // radial core charge density, rho_atc[0:mesh-1]
+    std::vector<double> rho_atc = {}; // radial core charge density, rho_atc[0:mesh-1]
 
     //<PP_RHOATOM>
-    double* rho_at = nullptr; // radial atomic charge density, rho_at[0:mesh-1]
+    std::vector<double> rho_at = {}; // radial atomic charge density, rho_at[0:mesh-1]
 
     // <PP_PSWFC>
     ModuleBase::matrix chi; // radial atomic orbitals, chi(nchi, mesh)
@@ -62,7 +62,7 @@ class pseudo
     double rcut = 0.0; // cut-off radius
 
     // <PP_BETA>
-    int* lll = nullptr; // lll(nbeta), angular momentum of the beta function
+    std::vector<int> lll = {}; // lll(nbeta), angular momentum of the beta function
 
     // <PP_DIJ>
     ModuleBase::matrix dion;  // dion(nbeta,nbeta)
diff --git a/source/module_cell/read_pp.cpp b/source/module_cell/read_pp.cpp
index fa3b6ed83a..f314be2ec7 100644
--- a/source/module_cell/read_pp.cpp
+++ b/source/module_cell/read_pp.cpp
@@ -16,15 +16,6 @@ Pseudopot_upf::Pseudopot_upf()
 
 Pseudopot_upf::~Pseudopot_upf()
 {
-    delete[] kbeta;
-    delete[] els_beta;
-    delete[] nchi;
-    delete[] epseu;
-    delete[] rcut_chi;
-    delete[] rcutus_chi;
-    delete[] rinner;
-    delete[] rcut;
-    delete[] rcutus;
 }
 
 int Pseudopot_upf::init_pseudo_reader(const std::string &fn, std::string &type, Atom_pseudo& pp)
@@ -460,7 +451,7 @@ void Pseudopot_upf::set_upf_q(Atom_pseudo& pp)
                                 break;
                             }
                         }
-                        this->setqfnew(nqf, ilast, l, 2, &(qfcoef(nb, mb, l, 0)), pp.r, &(pp.qfuncl(l, nmb, 0)));
+                        this->setqfnew(nqf, ilast, l, 2, &(qfcoef(nb, mb, l, 0)), pp.r.data(), &(pp.qfuncl(l, nmb, 0)));
                     }
                 }
             }
diff --git a/source/module_cell/read_pp.h b/source/module_cell/read_pp.h
index 3f4160bdeb..956690b6b5 100644
--- a/source/module_cell/read_pp.h
+++ b/source/module_cell/read_pp.h
@@ -4,11 +4,8 @@
 #include <string>
 
 #include "atom_pseudo.h"
-#include "module_base/global_function.h"
-#include "module_base/global_variable.h"
 #include "module_base/matrix.h"
 #include "module_base/realarray.h"
-#include "module_io/output.h"
 
 class Pseudopot_upf
 {
@@ -43,19 +40,19 @@ class Pseudopot_upf
     // need 'new' and 'delete'
     bool coulomb_potential = false;  // coulomb potentail : z/r
     ModuleBase::matrix chi;          // chi(nwfc,mesh) atomic wavefcts
-    int* kbeta = nullptr;            // kbeta(nbeta):number of mesh points for projector i (must be .le. mesh )
-    std::string* els_beta = nullptr; // els_beta(nwfc):label for the beta
-    int* nchi = nullptr;             // nchi(nwfc) value of pseudo-n for wavefcts
-    double* epseu = nullptr;         // epseu(nwfc) pseudo one-particle energy
-    double* rcut_chi = nullptr;      // rcut_chi(nwfc) cutoff inner radius
-    double* rcutus_chi = nullptr;    // rcutus_chi(nwfc) ultrasoft outer radius
-    double* rinner = nullptr;        // rinner(2*lmax+1) r_L
+    std::vector<int> kbeta = {};            // kbeta(nbeta):number of mesh points for projector i (must be .le. mesh )
+    std::vector<std::string> els_beta = {}; // els_beta(nwfc):label for the beta
+    std::vector<int> nchi = {};             // nchi(nwfc) value of pseudo-n for wavefcts
+    std::vector<double> epseu = {};         // epseu(nwfc) pseudo one-particle energy
+    std::vector<double> rcut_chi = {};      // rcut_chi(nwfc) cutoff inner radius
+    std::vector<double> rcutus_chi = {};    // rcutus_chi(nwfc) ultrasoft outer radius
+    std::vector<double> rinner = {};        // rinner(2*lmax+1) r_L
     ModuleBase::matrix qfunc;        // qfunc(nbeta*(nbeta+1)/2,mesh) Q_{mu,nu}(|r|) function for |r|> r_L
     ModuleBase::realArray qfcoef;    // qfcoef(nbeta,nbeta,2*lmax+1,nqf) coefficients for Q for |r|<r_L
     // ModuleBase::matrix aewfc;        // wfc(nbeta,mesh) all-electron wfc
     // ModuleBase::matrix pswfc;        // wfc(nbeta,mesh) pseudo wfc
-    double* rcut = nullptr;          // cut-off radius(nbeta)
-    double* rcutus = nullptr;        // ultrasoft cut-off radius (nbeta)
+    std::vector<double> rcut = {};          // cut-off radius(nbeta)
+    std::vector<double> rcutus = {};        // ultrasoft cut-off radius (nbeta)
 
     int nd; // nl_5 // Number of nonzero Dij
 
diff --git a/source/module_cell/read_pp_blps.cpp b/source/module_cell/read_pp_blps.cpp
index f8425ba388..ca8a54c001 100644
--- a/source/module_cell/read_pp_blps.cpp
+++ b/source/module_cell/read_pp_blps.cpp
@@ -11,23 +11,14 @@ int Pseudopot_upf::read_pseudo_blps(std::ifstream &ifs, Atom_pseudo& pp)
 
     pp.nbeta = 0;
     pp.kkbeta = 0;
-    delete[] this->kbeta;
-    delete[] pp.lll;
-    this->kbeta = nullptr;
-    pp.lll = nullptr;
-    pp.betar.create(1, 1);
-    pp.dion.create(1, 1);
+    pp.lll = std::vector<int>(pp.nbeta, 0);
+    pp.betar.create(0, 0);
+    pp.dion.create(pp.nbeta, pp.nbeta);
 
     pp.nchi = 0;
-    delete[] pp.nn;
-    delete[] pp.jchi;
-    delete[] pp.jjj;
-    pp.nn = new int[pp.nchi];
-    pp.jchi = new double[pp.nchi];
-    pp.jjj = new double[pp.nchi];
-    ModuleBase::GlobalFunc::ZEROS(pp.nn, pp.nchi);
-    ModuleBase::GlobalFunc::ZEROS(pp.jchi, pp.nchi);
-    ModuleBase::GlobalFunc::ZEROS(pp.jjj, pp.nchi);
+    pp.nn = std::vector<int>(pp.nchi, 0);
+    pp.jchi = std::vector<double>(pp.nchi, 0.0);
+    pp.jjj = std::vector<double>(pp.nchi, 0.0);
 
     ifs >> pp.psd;
     // if(!SCAN_BEGIN(ifs,"BLPS")) WARNING_QUIT("read_pp_blps","Find no PP_HEADER");
@@ -92,15 +83,9 @@ int Pseudopot_upf::read_pseudo_blps(std::ifstream &ifs, Atom_pseudo& pp)
 
     assert(pp.mesh > 0);
 
-    delete[] pp.r;
-    delete[] pp.rab;
-    delete[] pp.vloc_at;
-    pp.r = new double[pp.mesh]; // Bohr
-    pp.rab = new double[pp.mesh];
-    pp.vloc_at = new double[pp.mesh]; // Hartree
-    ModuleBase::GlobalFunc::ZEROS(pp.r, pp.mesh);
-    ModuleBase::GlobalFunc::ZEROS(pp.rab, pp.mesh);
-    ModuleBase::GlobalFunc::ZEROS(pp.vloc_at, pp.mesh);
+    pp.r = std::vector<double>(pp.mesh, 0.0); // Bohr
+    pp.rab = std::vector<double>(pp.mesh, 0.0);
+    pp.vloc_at = std::vector<double>(pp.mesh, 0.0); // Hartree
     int num = 0;
     if (pspcod == 8)
     {
@@ -126,9 +111,7 @@ int Pseudopot_upf::read_pseudo_blps(std::ifstream &ifs, Atom_pseudo& pp)
     }
     pp.rab[pp.mesh - 1] = pp.r[pp.mesh - 1] - pp.r[pp.mesh - 2];
 
-    delete[] pp.rho_at;
-    pp.rho_at = new double[pp.mesh];
-    ModuleBase::GlobalFunc::ZEROS(pp.rho_at, pp.mesh);
+    pp.rho_at = std::vector<double>(pp.mesh, 0.0);
     double charge = zion/pp.r[pp.mesh - 1];
     for(int i = 0;i < pp.mesh; ++i)
     {
diff --git a/source/module_cell/read_pp_complete.cpp b/source/module_cell/read_pp_complete.cpp
index f9b79d8187..7bfdae210d 100644
--- a/source/module_cell/read_pp_complete.cpp
+++ b/source/module_cell/read_pp_complete.cpp
@@ -14,10 +14,9 @@ void Pseudopot_upf::complete_default(Atom_pseudo& pp)
 		return;
 	}
 	
-	if (pp.lll == nullptr)
+	if (pp.lll.empty())
 	{
-		pp.lll = new int[pp.nbeta];
-		assert(pp.lll != nullptr);
+		pp.lll = std::vector<int>(pp.nbeta, 0);
 	}
 
 	pp.nh = 0;
@@ -42,48 +41,42 @@ void Pseudopot_upf::complete_default_h(Atom_pseudo& pp)
 		std::cout << "\n complete_default_h, too many grid points,";
 	}
 
-	if (pp.els == nullptr)
+	if (pp.els.empty())
 	{
-		pp.els = new std::string[pp.nchi];
-		assert(pp.els != nullptr);
+		pp.els = std::vector<std::string>(pp.nchi, "");
 	}
 
-	if (pp.lchi == nullptr)
+	if (pp.lchi.empty())
 	{
-		pp.lchi = new int[pp.nchi];
-		assert(pp.lchi != nullptr);
+		pp.lchi = std::vector<int>(pp.nchi, 0);
 	}
 
-	if (pp.oc == nullptr)
+	if (pp.oc.empty())
 	{
-		pp.oc = new double[pp.nchi];
-		assert(pp.oc != nullptr);
+		pp.oc = std::vector<double>(pp.nchi, 0.0);
 	}
 
-	if (pp.jjj == nullptr) {
-		pp.jjj = new double[pp.nbeta];
-		assert(pp.jjj != nullptr);
-		assert(!pp.has_so);
+	if (pp.jjj.empty()) {
+		pp.jjj = std::vector<double>(pp.nbeta, 0.0);
+		assert(!pp.has_so or pp.nbeta == 0);
 		for (int i=0; i<pp.nbeta; i++)
 		{
 			pp.jjj[i]  = 0;
 		}
 	}
 
-	if (pp.nn == nullptr) {
-		pp.nn = new int[pp.nchi];
-		assert(pp.nn != nullptr);
-		assert(!pp.has_so);
+	if (pp.nn.empty()) {
+		pp.nn = std::vector<int>(pp.nchi, 0);
+		assert(!pp.has_so or pp.nchi == 0);
 		for (int i=0; i<pp.nchi; i++)
 		{
 			pp.nn[i] = 0;
 		}
 	}
 
-	if (pp.jchi == nullptr) {
-		pp.jchi = new double[pp.nchi];
-		assert(pp.jchi != nullptr);
-		assert(!pp.has_so);
+	if (pp.jchi.empty()) {
+		pp.jchi = std::vector<double>(pp.nchi, 0.0);
+		assert(!pp.has_so or pp.nchi == 0);
 		for (int i=0; i<pp.nchi; i++)
 		{
 			pp.jchi[i] = 0;
@@ -105,7 +98,7 @@ void Pseudopot_upf::complete_default_atom(Atom_pseudo& pp)
 	// remember to update here if you need it.
 	//	rcut = 25.0; 
 
-	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"PAO radial cut off (Bohr)",rcut);
+	ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"PAO radial cut off (Bohr)", pp.rcut);
 	if(pp.rcut <= 0.0)
 	{
 		ModuleBase::WARNING_QUIT("Pseudopot_upf::complete_default_atom","PAO rcut<=0.0");
@@ -113,29 +106,21 @@ void Pseudopot_upf::complete_default_atom(Atom_pseudo& pp)
 
 	// chi.create(nchi, mesh);
 
-	if (pp.r == nullptr) {
-		pp.r = new double[pp.mesh];
-		assert(pp.r != nullptr);
-		ModuleBase::GlobalFunc::ZEROS(pp.r, pp.mesh);
+	if (pp.r.empty()) {
+		pp.r = std::vector<double>(pp.mesh, 0.0);
 	}
 
-	if (pp.rab == nullptr) {
-		pp.rab = new double[pp.mesh];
-		assert(pp.rab != nullptr);
-		ModuleBase::GlobalFunc::ZEROS(pp.rab, pp.mesh);
+	if (pp.rab.empty()) {
+		pp.rab = std::vector<double>(pp.mesh, 0.0);
 	}
 
-	if (pp.rho_at == nullptr) {
-		pp.rho_at = new double[pp.mesh];
-		assert(pp.rho_at != nullptr);
-		ModuleBase::GlobalFunc::ZEROS(pp.rho_at, pp.mesh);
+	if (pp.rho_at.empty()) {
+		pp.rho_at = std::vector<double>(pp.mesh, 0.0);
 	}
 
-	if (pp.rho_atc == nullptr) {
-		pp.rho_atc = new double[pp.mesh];
-		assert(pp.rho_atc != nullptr);
-		assert(!pp.nlcc);
-		ModuleBase::GlobalFunc::ZEROS(pp.rho_atc, pp.mesh);
+	if (pp.rho_atc.empty()) {
+		pp.rho_atc = std::vector<double>(pp.mesh, 0.0);
+		assert(!pp.nlcc or pp.mesh == 0);
 	}
 
 	bool br = false;
@@ -171,9 +156,8 @@ void Pseudopot_upf::complete_default_vl(Atom_pseudo& pp)
 
 	assert(pp.mesh>0);//mohan add 2021-05-01
 
-	if (pp.vloc_at == nullptr) {
-		pp.vloc_at = new double[pp.mesh];
-		assert(pp.vloc_at != nullptr);
+	if (pp.vloc_at.empty()) {
+		pp.vloc_at = std::vector<double>(pp.mesh, 0.0);
 	}
 
 	return;
diff --git a/source/module_cell/read_pp_upf100.cpp b/source/module_cell/read_pp_upf100.cpp
index 8f4bfd8feb..743df13235 100644
--- a/source/module_cell/read_pp_upf100.cpp
+++ b/source/module_cell/read_pp_upf100.cpp
@@ -202,15 +202,9 @@ void Pseudopot_upf::read_pseudo_header(std::ifstream &ifs, Atom_pseudo& pp)
 	ifs.ignore(75, '\n');
 	ifs.ignore(75, '\n');
 
-	delete[] pp.els;
-	delete[] pp.lchi;
-	delete[] pp.oc;
-	pp.els = new std::string[pp.nchi];
-	pp.lchi = new int[pp.nchi];
-	pp.oc = new double[pp.nchi];
-
-	ModuleBase::GlobalFunc::ZEROS(pp.lchi, pp.nchi); // angular momentum of each orbital
-	ModuleBase::GlobalFunc::ZEROS(pp.oc, pp.nchi);//occupation of each orbital
+	pp.els = std::vector<std::string>(pp.nchi, "");
+	pp.lchi = std::vector<int>(pp.nchi, 0);
+	pp.oc = std::vector<double>(pp.nchi, 0.0);
 
 	for(int i=0;i<pp.nchi;i++)
 	{
@@ -229,12 +223,8 @@ void Pseudopot_upf::read_pseudo_mesh(std::ifstream &ifs, Atom_pseudo& pp)
 {
 	assert(pp.mesh>0);
 
-	delete[] pp.r;
-	delete[] pp.rab;
-	pp.r = new double[pp.mesh];
-	pp.rab = new double[pp.mesh];
-	ModuleBase::GlobalFunc::ZEROS(pp.r,pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(pp.rab,pp.mesh);
+	pp.r = std::vector<double>(pp.mesh, 0.0);
+	pp.rab = std::vector<double>(pp.mesh, 0.0);
 
 	int ir = 0;
 
@@ -261,9 +251,7 @@ void Pseudopot_upf::read_pseudo_mesh(std::ifstream &ifs, Atom_pseudo& pp)
 void Pseudopot_upf::read_pseudo_nlcc(std::ifstream &ifs, Atom_pseudo& pp)
 {
 	assert(pp.mesh>0);
-	delete[] pp.rho_atc;
-	pp.rho_atc = new double[pp.mesh];
-	ModuleBase::GlobalFunc::ZEROS(pp.rho_atc, pp.mesh);
+	pp.rho_atc = std::vector<double>(pp.mesh, 0.0);
 	for (int ir = 0;ir < pp.mesh;ir++)
 	{
 		ifs >> pp.rho_atc[ir];
@@ -274,9 +262,7 @@ void Pseudopot_upf::read_pseudo_nlcc(std::ifstream &ifs, Atom_pseudo& pp)
 void Pseudopot_upf::read_pseudo_local(std::ifstream &ifs, Atom_pseudo& pp)
 {
 	assert(pp.mesh>0);
-	delete[] pp.vloc_at;
-	pp.vloc_at = new double[pp.mesh];
-	ModuleBase::GlobalFunc::ZEROS(pp.vloc_at, pp.mesh);
+	pp.vloc_at = std::vector<double>(pp.mesh, 0.0);
 
 	for (int ir = 0;ir < pp.mesh;ir++)
 	{
@@ -304,10 +290,8 @@ void Pseudopot_upf::read_pseudo_nl(std::ifstream &ifs, Atom_pseudo& pp)
     }
     else
     {
-        delete[] this->kbeta;
-        delete[] pp.lll;
-        this->kbeta = new int[pp.nbeta];
-        pp.lll = new int[pp.nbeta];
+        this->kbeta = std::vector<int>(pp.nbeta, 0);
+        pp.lll = std::vector<int>(pp.nbeta, 0);
         pp.betar.create(pp.nbeta, pp.mesh);
         pp.dion.create(pp.nbeta, pp.nbeta);
         pp.kkbeta = 0;
@@ -354,8 +338,7 @@ void Pseudopot_upf::read_pseudo_nl(std::ifstream &ifs, Atom_pseudo& pp)
             // If nqf is not zero, Qij's inside rinner are computed using qfcoef's
             ModuleBase::GlobalFunc::READ_VALUE(ifs, this->nqf);
             pp.nqlc = 2 * pp.lmax + 1;
-            delete[] rinner;
-            this->rinner = new double[pp.nqlc];
+            this->rinner = std::vector<double>(pp.nqlc, 0.0);
             pp.qqq.create(pp.nbeta, pp.nbeta);
             if (q_with_l)
             {
@@ -368,7 +351,6 @@ void Pseudopot_upf::read_pseudo_nl(std::ifstream &ifs, Atom_pseudo& pp)
 
             if (nqf <= 0)
             {
-                ModuleBase::GlobalFunc::ZEROS(rinner, pp.nqlc);
                 this->qfcoef.create(1, 1, 1, 1);
             }
             else
@@ -463,9 +445,7 @@ void Pseudopot_upf::read_pseudo_pswfc(std::ifstream &ifs, Atom_pseudo& pp)
 
 void Pseudopot_upf::read_pseudo_rhoatom(std::ifstream &ifs, Atom_pseudo& pp)
 {
-	delete[] pp.rho_at;
-	pp.rho_at = new double[pp.mesh];
-	ModuleBase::GlobalFunc::ZEROS(pp.rho_at, pp.mesh);
+	pp.rho_at = std::vector<double>(pp.mesh, 0.0);
 	for (int ir = 0;ir < pp.mesh;ir++)
 	{
 		ifs >> pp.rho_at[ir];
@@ -477,15 +457,9 @@ void Pseudopot_upf::read_pseudo_so(std::ifstream &ifs, Atom_pseudo& pp)
 {
        //read soc info from upf, added by zhengdy-soc
        if(!pp.has_so) return;
-       delete[] pp.nn;
-       delete[] pp.jchi;
-       delete[] pp.jjj;
-       pp.nn = new int[pp.nchi];
-       pp.jchi = new double[pp.nchi];
-       pp.jjj = new double[pp.nbeta];
-       ModuleBase::GlobalFunc::ZEROS(pp.nn,pp.nchi);
-       ModuleBase::GlobalFunc::ZEROS(pp.jchi,pp.nchi);
-       ModuleBase::GlobalFunc::ZEROS(pp.jjj,pp.nbeta);
+       pp.nn = std::vector<int>(pp.nchi, 0);
+       pp.jchi = std::vector<double>(pp.nchi, 0.0);
+       pp.jjj = std::vector<double>(pp.nbeta, 0.0);
        //RELWFC
        for(int nw=0;nw< pp.nchi;nw++)
        {
diff --git a/source/module_cell/read_pp_upf201.cpp b/source/module_cell/read_pp_upf201.cpp
index b665557ffa..1d96077046 100644
--- a/source/module_cell/read_pp_upf201.cpp
+++ b/source/module_cell/read_pp_upf201.cpp
@@ -27,9 +27,7 @@ int Pseudopot_upf::read_pseudo_upf201(std::ifstream &ifs, Atom_pseudo& pp)
         {
             ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_NLCC>");
         }
-        delete[] pp.rho_atc;
-        pp.rho_atc = new double[pp.mesh];
-        ModuleBase::GlobalFunc::ZEROS(pp.rho_atc, pp.mesh);
+        pp.rho_atc = std::vector<double>(pp.mesh, 0.0);
         for (int ir = 0; ir < pp.mesh; ir++)
         {
             ifs >> pp.rho_atc[ir];
@@ -50,9 +48,7 @@ int Pseudopot_upf::read_pseudo_upf201(std::ifstream &ifs, Atom_pseudo& pp)
         {
             ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_LOCAL>");
         }
-        delete[] pp.vloc_at;
-        pp.vloc_at = new double[pp.mesh];
-        ModuleBase::GlobalFunc::ZEROS(pp.vloc_at, pp.mesh);
+        pp.vloc_at = std::vector<double>(pp.mesh, 0.0);
         for (int ir = 0; ir < pp.mesh; ir++)
         {
             ifs >> pp.vloc_at[ir];
@@ -89,9 +85,7 @@ int Pseudopot_upf::read_pseudo_upf201(std::ifstream &ifs, Atom_pseudo& pp)
     {
         ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_RHOATOM>");
     }
-    delete[] pp.rho_at;
-    pp.rho_at = new double[pp.mesh];
-    ModuleBase::GlobalFunc::ZEROS(pp.rho_at, pp.mesh);
+    pp.rho_at = std::vector<double>(pp.mesh, 0.0);
     for (int ir = 0; ir < pp.mesh; ir++)
     {
         ifs >> pp.rho_at[ir];
@@ -401,13 +395,9 @@ void Pseudopot_upf::read_pseudo_upf201_mesh(std::ifstream& ifs, Atom_pseudo& pp)
     {
         ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_R>");
     }
-    delete[] pp.r;
-    delete[] pp.rab;
     assert(pp.mesh > 0);
-    pp.r = new double[pp.mesh];
-    pp.rab = new double[pp.mesh];
-    ModuleBase::GlobalFunc::ZEROS(pp.r, pp.mesh);
-    ModuleBase::GlobalFunc::ZEROS(pp.rab, pp.mesh);
+    pp.r = std::vector<double>(pp.mesh, 0.0);
+    pp.rab = std::vector<double>(pp.mesh, 0.0);
     for (int ir = 0; ir < pp.mesh; ir++)
     {
         ifs >> pp.r[ir];
@@ -441,16 +431,11 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs, Atom_pseudo&
     std::string name[50];
     std::string val[50];
     int nparameter;
-    delete[] this->kbeta;
-    delete[] pp.lll;
-    delete[] els_beta;
-    delete[] rcut;
-    delete[] rcutus;
-    this->kbeta = new int[pp.nbeta];
-    pp.lll = new int[pp.nbeta];
-    this->els_beta = new std::string[pp.nbeta];
-    this->rcut = new double[pp.nbeta];
-    this->rcutus = new double[pp.nbeta];
+    this->kbeta = std::vector<int>(pp.nbeta);
+    pp.lll = std::vector<int>(pp.nbeta);
+    this->els_beta = std::vector<std::string>(pp.nbeta);
+    this->rcut = std::vector<double>(pp.nbeta, 0.0);
+    this->rcutus = std::vector<double>(pp.nbeta, 0.0);
     pp.betar.create(pp.nbeta, pp.mesh);
     pp.dion.create(pp.nbeta, pp.nbeta);
     for (int ib = 0; ib < pp.nbeta; ib++)
@@ -586,11 +571,9 @@ void Pseudopot_upf::read_pseudo_upf201_nonlocal(std::ifstream& ifs, Atom_pseudo&
         ModuleBase::GlobalFunc::SCAN_END(ifs, "</PP_Q>");
 
         // Read polinomial coefficients for Q_ij expansion at small radius
-        delete[] rinner;
-        this->rinner = new double[pp.nqlc];
+        this->rinner = std::vector<double>(pp.nqlc, 0.0);
         if (nqf <= 0)
         {
-            ModuleBase::GlobalFunc::ZEROS(rinner, pp.nqlc);
             this->qfcoef.create(1, 1, 1, 1);
         }
         else
@@ -689,22 +672,13 @@ void Pseudopot_upf::read_pseudo_upf201_pswfc(std::ifstream& ifs, Atom_pseudo& pp
     std::string val[50];
     int nparameter;
     ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_PSWFC>");
-    delete[] pp.els;
-    delete[] pp.lchi;
-    delete[] nchi;
-    delete[] pp.oc;
-    delete[] epseu;
-    delete[] rcut_chi;
-    delete[] rcutus_chi;
-    pp.els = new std::string[pp.nchi];
-    pp.lchi = new int[pp.nchi];
-    this->nchi = new int[pp.nchi];
-    pp.oc = new double[pp.nchi];
-    this->epseu = new double[pp.nchi];
-    this->rcut_chi = new double[pp.nchi];
-    this->rcutus_chi = new double[pp.nchi];
-    ModuleBase::GlobalFunc::ZEROS(pp.lchi, pp.nchi); // angular momentum of each orbital
-    ModuleBase::GlobalFunc::ZEROS(pp.oc, pp.nchi);   // occupation of each orbital
+    pp.els = std::vector<std::string>(pp.nchi, "");
+    pp.lchi = std::vector<int>(pp.nchi, 0);
+    this->nchi = std::vector<int>(pp.nchi, 0);
+    pp.oc = std::vector<double>(pp.nchi, 0.0);
+    this->epseu = std::vector<double>(pp.nchi, 0.0);
+    this->rcut_chi = std::vector<double>(pp.nchi, 0.0);
+    this->rcutus_chi = std::vector<double>(pp.nchi, 0.0);
     pp.chi.create(pp.nchi, pp.mesh);
     for (int iw = 0; iw < pp.nchi; iw++)
     {
@@ -832,15 +806,9 @@ void Pseudopot_upf::read_pseudo_upf201_so(std::ifstream& ifs, Atom_pseudo& pp)
     std::string val[50];
     int nparameter;
     ModuleBase::GlobalFunc::SCAN_BEGIN(ifs, "<PP_SPIN_ORB>");
-    delete[] pp.jchi;
-    delete[] pp.jjj;
-    delete[] pp.nn;
-    pp.jchi = new double[pp.nchi];
-    pp.jjj = new double[pp.nbeta];
-    pp.nn = new int[pp.nchi];
-    ModuleBase::GlobalFunc::ZEROS(pp.jchi, pp.nchi);
-    ModuleBase::GlobalFunc::ZEROS(pp.jjj, pp.nbeta);
-    ModuleBase::GlobalFunc::ZEROS(pp.nn, pp.nchi);
+    pp.jchi = std::vector<double>(pp.nchi, 0.0);
+    pp.jjj = std::vector<double>(pp.nbeta, 0.0);
+    pp.nn = std::vector<int>(pp.nchi, 0);
 
     for (int nw = 0; nw < pp.nchi; nw++)
     {
diff --git a/source/module_cell/read_pp_vwr.cpp b/source/module_cell/read_pp_vwr.cpp
index 15e6f4ae59..29a241557d 100644
--- a/source/module_cell/read_pp_vwr.cpp
+++ b/source/module_cell/read_pp_vwr.cpp
@@ -56,7 +56,7 @@ int Pseudopot_upf::read_pseudo_vwr(std::ifstream &ifs, Atom_pseudo& pp)
 	spd_loc = std::atoi( value.c_str() );
 	GlobalV::ofs_running << std::setw(15) << "LOC(spd)" << std::setw(15) << spd_loc << std::endl;
 	// (6) read in the occupations
-	double* tmp_oc = new double[3];
+	std::vector<double> tmp_oc(3, 0.0);
 	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);
@@ -90,51 +90,31 @@ int Pseudopot_upf::read_pseudo_vwr(std::ifstream &ifs, Atom_pseudo& pp)
 	if(iTB_d) ++pp.nchi;
 	GlobalV::ofs_running << std::setw(15) << "NWFC" << std::setw(15) << pp.nchi << std::endl;
 	// allocate occupation number array for wave functions
-	delete[] pp.oc;
-	delete[] pp.els;
-	pp.oc = new double[pp.nchi];
-	pp.els = new std::string[pp.nchi];
+	pp.oc = std::vector<double>(pp.nchi, 0.0);
+	pp.els = std::vector<std::string>(pp.nchi, "");
 	// set the value of occupations
-	delete[] pp.lchi;
-	pp.lchi = new int[pp.nchi];
+	pp.lchi = std::vector<int>(pp.nchi, 0);
 	int iwfc=0;
 	if(iTB_s){pp.oc[iwfc]=tmp_oc[0];pp.lchi[iwfc]=0;pp.els[iwfc]="S";++iwfc;}
 	if(iTB_p){pp.oc[iwfc]=tmp_oc[1];pp.lchi[iwfc]=1;pp.els[iwfc]="P";++iwfc;}
 	if(iTB_d){pp.oc[iwfc]=tmp_oc[2];pp.lchi[iwfc]=2;pp.els[iwfc]="D";++iwfc;}
-	delete[] tmp_oc;
 	getline(ifs,value);
 
 
 	// global variables that will be used
 	// in other classes.
-	delete[] pp.r;
-	delete[] pp.rab;
-	delete[] pp.vloc_at;
-	pp.r = new double[pp.mesh];
-	pp.rab = new double[pp.mesh];
-	pp.vloc_at = new double[pp.mesh];
-	ModuleBase::GlobalFunc::ZEROS(pp.r,pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(pp.rab,pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(pp.vloc_at,pp.mesh);
-	delete[] pp.rho_at;
-	delete[] pp.rho_atc;
-	pp.rho_at = new double[pp.mesh];
-	pp.rho_atc = new double[pp.mesh];
-	ModuleBase::GlobalFunc::ZEROS(pp.rho_at, pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(pp.rho_atc, pp.mesh);
+	pp.r = std::vector<double>(pp.mesh, 0.0);
+	pp.rab = std::vector<double>(pp.mesh, 0.0);
+	pp.vloc_at = std::vector<double>(pp.mesh, 0.0);
+	pp.rho_at = std::vector<double>(pp.mesh, 0.0);
+	pp.rho_atc = std::vector<double>(pp.mesh, 0.0);
 	// local variables in this function
-    double* vs = new double[pp.mesh]; // local pseudopotential for s, unit is Hartree
-    double* vp = new double[pp.mesh]; // local pseudopotential for p
-    double* vd = new double[pp.mesh]; // local pseudopotential for d
-    double* ws = new double[pp.mesh]; // wave function for s
-    double* wp = new double[pp.mesh]; // wave function for p
-    double* wd = new double[pp.mesh]; // wave function for d
-    ModuleBase::GlobalFunc::ZEROS(vs,pp.mesh);
-    ModuleBase::GlobalFunc::ZEROS(vp,pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(vd,pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(ws,pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(wp,pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(wd,pp.mesh);
+    std::vector<double> vs = std::vector<double>(pp.mesh, 0.0); // local pseudopotential for s, unit is Hartree
+    std::vector<double> vp = std::vector<double>(pp.mesh, 0.0); // local pseudopotential for p
+    std::vector<double> vd = std::vector<double>(pp.mesh, 0.0); // local pseudopotential for d
+    std::vector<double> ws = std::vector<double>(pp.mesh, 0.0); // wave function for s
+    std::vector<double> wp = std::vector<double>(pp.mesh, 0.0); // wave function for p
+    std::vector<double> wd = std::vector<double>(pp.mesh, 0.0); // wave function for d
 	std::string line;
 	if(spd_loc>0 && pp.nlcc==0)
 	{
@@ -280,8 +260,7 @@ int Pseudopot_upf::read_pseudo_vwr(std::ifstream &ifs, Atom_pseudo& pp)
 	this->nd = pp.nbeta;
 	GlobalV::ofs_running << std::setw(15) << "N-Dij" << std::setw(15) << nd << std::endl;
 	// calculate the angular momentum for each pp.betar
-	delete[] pp.lll;
-	pp.lll = new int[pp.nbeta]; 
+	pp.lll = std::vector<int>(pp.nbeta, 0); 
 	int icount=0;
 	if(iref_s==1) {pp.lll[icount]=0; ++icount;}// s projector
 	if(iref_p==1) {pp.lll[icount]=1; ++icount;}// p projector
@@ -291,8 +270,7 @@ int Pseudopot_upf::read_pseudo_vwr(std::ifstream &ifs, Atom_pseudo& pp)
 		GlobalV::ofs_running << " lll[" << i << "]=" << pp.lll[i] << std::endl;
 	}
     // this->kbeta(pp.nbeta): number of mesh points for projector i (must be .le.mesh )
-    delete[] this->kbeta;
-    this->kbeta = new int[pp.nbeta];
+    this->kbeta = std::vector<int>(pp.nbeta, 0);
     pp.kkbeta = 0;
     for (int ib = 0; ib < pp.nbeta; ++ib)
     {
@@ -309,14 +287,11 @@ int Pseudopot_upf::read_pseudo_vwr(std::ifstream &ifs, Atom_pseudo& pp)
 	// (6) generate nonlocal pseudopotentials
 	// --------------------------------------
 	// tmp function to evaluate < pp.betar | delta_v | pp.betar>
-	double* func = new double[pp.mesh];
+	std::vector<double> func = std::vector<double>(pp.mesh, 0.0);
 	// tmp value (vs, vp or vd)
-	double* vl = new double[pp.mesh];
+	std::vector<double> vl = std::vector<double>(pp.mesh, 0.0);
 	// tmp wave function (ws, wp or wd with r)
-	double* wlr = new double[pp.mesh];
-	ModuleBase::GlobalFunc::ZEROS(func, pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(vl, pp.mesh);
-	ModuleBase::GlobalFunc::ZEROS(wlr, pp.mesh);
+	std::vector<double> wlr = std::vector<double>(pp.mesh, 0.0);
 	double rcut = 5.0/1.03;
 	GlobalV::ofs_running << std::setw(15) << "RCUT_NL" << std::setw(15) << rcut << std::endl;
 	for(int ib=0; ib<pp.nbeta; ++ib)
@@ -383,22 +358,5 @@ int Pseudopot_upf::read_pseudo_vwr(std::ifstream &ifs, Atom_pseudo& pp)
 	}
 	*/
 	GlobalV::ofs_running << " -------------------------------------------------" << std::endl;
-
-
-	// --------------------------------------
-	// (5) clean up 
-	// --------------------------------------
-	delete[] func;
-	delete[] vl;
-	delete[] wlr;
-
-	delete[] vs;
-	delete[] vp;
-	delete[] vd;
-
-	delete[] ws;
-	delete[] wp;
-	delete[] wd;
-
 	return 0;
 }
diff --git a/source/module_cell/setup_nonlocal.cpp b/source/module_cell/setup_nonlocal.cpp
index ed82a03673..96c3ac33b0 100644
--- a/source/module_cell/setup_nonlocal.cpp
+++ b/source/module_cell/setup_nonlocal.cpp
@@ -134,8 +134,8 @@ void InfoNonlocal::Set_NonLocal(const int& it,
                                    it,       // type
                                    lnow,     // angular momentum L
                                    cut_mesh, // number of radial mesh
-                                   atom->ncpp.rab,
-                                   atom->ncpp.r, // radial mesh value (a.u.)
+                                   atom->ncpp.rab.data(),
+                                   atom->ncpp.r.data(), // radial mesh value (a.u.)
                                    beta_r,
                                    kmesh,
                                    dk,
diff --git a/source/module_cell/test/read_pp_test.cpp b/source/module_cell/test/read_pp_test.cpp
index 46d81c4785..81ce976dc3 100644
--- a/source/module_cell/test/read_pp_test.cpp
+++ b/source/module_cell/test/read_pp_test.cpp
@@ -76,7 +76,7 @@ TEST_F(ReadPPTest, ReadUPF100_Coulomb)
 	std::ifstream ifs;
 	ifs.open("./support/Te.pbe-coulomb.UPF");
 	read_pp->read_pseudo_upf(ifs, *upf);
-	EXPECT_EQ(upf->vloc_at, nullptr);
+	EXPECT_TRUE(upf->vloc_at.empty());
 	EXPECT_EQ(read_pp->coulomb_potential, true);
 	EXPECT_EQ(upf->tvanp, false);
 	EXPECT_EQ(upf->nbeta, 0);
@@ -287,7 +287,7 @@ TEST_F(ReadPPTest, ReadUPF201_Coulomb)
 	std::ifstream ifs;
 	ifs.open("./support/Al.pbe-coulomb.UPF");
 	read_pp->read_pseudo_upf201(ifs, *upf);
-	EXPECT_EQ(upf->vloc_at, nullptr);
+	EXPECT_TRUE(upf->vloc_at.empty());
 	EXPECT_EQ(read_pp->coulomb_potential, true);
 	EXPECT_EQ(upf->nbeta, 0);
 	EXPECT_EQ(upf->lmax, 0);
@@ -383,7 +383,7 @@ TEST_F(ReadPPTest, ReadUSPPUPF201)
     EXPECT_DOUBLE_EQ(upf->rab[892], 3.344685763390000e0);
     EXPECT_DOUBLE_EQ(upf->vloc_at[0], 3.456089057550000e0);
     EXPECT_DOUBLE_EQ(upf->vloc_at[892], -1.096266796840000e-1);
-    EXPECT_EQ(upf->rho_atc, nullptr);
+	EXPECT_TRUE(upf->rho_atc.empty());
     EXPECT_EQ(upf->lll[0], 0);
     EXPECT_EQ(read_pp->kbeta[0], 617);
     EXPECT_EQ(read_pp->els_beta[0], "2S");
@@ -419,9 +419,9 @@ TEST_F(ReadPPTest, ReadUSPPUPF201)
     EXPECT_DOUBLE_EQ(upf->chi(0, 892), 0.0);
     EXPECT_DOUBLE_EQ(upf->rho_at[0], 0.0);
     EXPECT_DOUBLE_EQ(upf->rho_at[892], 0.0);
-    EXPECT_EQ(upf->jchi, nullptr);
-    EXPECT_EQ(upf->jjj, nullptr);
-    EXPECT_EQ(upf->nn, nullptr);
+    EXPECT_TRUE(upf->jchi.empty());
+	EXPECT_TRUE(upf->jjj.empty());
+	EXPECT_TRUE(upf->nn.empty());
     ifs.close();
 }
 
@@ -634,8 +634,8 @@ TEST_F(ReadPPTest, SetEmptyElement)
 {
 	upf->mesh = 10;
 	upf->nbeta = 10;
-	upf->vloc_at = new double[upf->mesh];
-	upf->rho_at = new double[upf->mesh];
+	upf->vloc_at = std::vector<double>(upf->mesh, 0.0);
+	upf->rho_at = std::vector<double>(upf->mesh, 0.0);
 	upf->dion.create(upf->nbeta,upf->nbeta);
 	read_pp->set_empty_element(*upf);
 	for(int ir=0;ir<upf->mesh;++ir)
diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp
index e659c46213..400f13f624 100644
--- a/source/module_elecstate/module_charge/charge.cpp
+++ b/source/module_elecstate/module_charge/charge.cpp
@@ -372,8 +372,8 @@ void Charge::atomic_rho(const int spin_number_need,
 
                                 double charge = 0.0;
                                 ModuleBase::Integral::Simpson_Integral(atom->ncpp.msh,
-                                                                       atom->ncpp.rho_at,
-                                                                       atom->ncpp.rab,
+                                                                       atom->ncpp.rho_at.data(),
+                                                                       atom->ncpp.rab.data(),
                                                                        charge);
                                 ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning, "charge from rho_at", charge);
                                 assert(charge != 0.0
@@ -410,7 +410,7 @@ void Charge::atomic_rho(const int spin_number_need,
                             {
                                 rho1d[ir] = rhoatm[ir];
                             }
-                            ModuleBase::Integral::Simpson_Integral(mesh, rho1d.data(), atom->ncpp.rab, rho_lgl[0]);
+                            ModuleBase::Integral::Simpson_Integral(mesh, rho1d.data(), atom->ncpp.rab.data(), rho_lgl[0]);
                             gstart = 1;
                         }
                         if (GlobalV::test_charge > 0)
@@ -445,7 +445,7 @@ void Charge::atomic_rho(const int spin_number_need,
                                         rho1d[ir] = rhoatm[ir] * ModuleBase::libm::sin(gxx) / gxx;
                                     }
                                 }
-                                ModuleBase::Integral::Simpson_Integral(mesh, rho1d.data(), atom->ncpp.rab, rho_lgl[igg]);
+                                ModuleBase::Integral::Simpson_Integral(mesh, rho1d.data(), atom->ncpp.rab.data(), rho_lgl[igg]);
                             }
     #ifdef _OPENMP
     #pragma omp single
diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp
index fd95ca2ff5..239025a845 100644
--- a/source/module_elecstate/module_charge/charge_init.cpp
+++ b/source/module_elecstate/module_charge/charge_init.cpp
@@ -254,9 +254,9 @@ void Charge::set_rho_core(
             this->non_linear_core_correction(
                 GlobalC::ppcell.numeric,
                 GlobalC::ucell.atoms[it].ncpp.msh,
-                GlobalC::ucell.atoms[it].ncpp.r,
-                GlobalC::ucell.atoms[it].ncpp.rab,
-                GlobalC::ucell.atoms[it].ncpp.rho_atc,
+                GlobalC::ucell.atoms[it].ncpp.r.data(),
+                GlobalC::ucell.atoms[it].ncpp.rab.data(),
+                GlobalC::ucell.atoms[it].ncpp.rho_atc.data(),
                 rhocg);
 //----------------------------------------------------------
 // EXPLAIN : multiply by the structure factor and sum
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp
index 45d095a7b2..3401593a36 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_T_NL_cd.cpp
@@ -67,10 +67,10 @@ class TNLTest : public ::testing::Test
         ucell.atoms[0].ncpp.non_zero_count_soc[1] = 0;
         ucell.atoms[0].ncpp.non_zero_count_soc[2] = 0;
         ucell.atoms[0].ncpp.non_zero_count_soc[3] = 5;
-        ucell.atoms[0].ncpp.index1_soc[0] = new int[5];
-        ucell.atoms[0].ncpp.index2_soc[0] = new int[5];
-        ucell.atoms[0].ncpp.index1_soc[3] = new int[5];
-        ucell.atoms[0].ncpp.index2_soc[3] = new int[5];
+        ucell.atoms[0].ncpp.index1_soc[0] = std::vector<int>(5, 0);
+        ucell.atoms[0].ncpp.index2_soc[0] = std::vector<int>(5, 0);
+        ucell.atoms[0].ncpp.index1_soc[3] = std::vector<int>(5, 0);
+        ucell.atoms[0].ncpp.index2_soc[3] = std::vector<int>(5, 0);
         for (int i = 0; i < 5; ++i)
         {
             ucell.atoms[0].ncpp.d_real(i, i) = 1.0;
@@ -95,10 +95,6 @@ class TNLTest : public ::testing::Test
         delete[] ucell.atoms[0].iw2l;
         delete[] ucell.atoms[0].iw2m;
         delete[] ucell.atoms[0].iw2n;
-        delete[] ucell.atoms[0].ncpp.index1_soc[0];
-        delete[] ucell.atoms[0].ncpp.index2_soc[0];
-        delete[] ucell.atoms[0].ncpp.index1_soc[3];
-        delete[] ucell.atoms[0].ncpp.index2_soc[3];
         delete[] ucell.atoms;
         delete[] ucell.iat2it;
         delete[] ucell.iat2ia;
diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp
index beb5925915..d5b00598d5 100644
--- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp
+++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/test/test_nonlocalnew.cpp
@@ -65,10 +65,10 @@ class NonlocalNewTest : public ::testing::Test
         ucell.atoms[0].ncpp.non_zero_count_soc[1] = 0;
         ucell.atoms[0].ncpp.non_zero_count_soc[2] = 0;
         ucell.atoms[0].ncpp.non_zero_count_soc[3] = 5;
-        ucell.atoms[0].ncpp.index1_soc[0] = new int[5];
-        ucell.atoms[0].ncpp.index2_soc[0] = new int[5];
-        ucell.atoms[0].ncpp.index1_soc[3] = new int[5];
-        ucell.atoms[0].ncpp.index2_soc[3] = new int[5];
+        ucell.atoms[0].ncpp.index1_soc[0] = std::vector<int>(5, 0);
+        ucell.atoms[0].ncpp.index2_soc[0] = std::vector<int>(5, 0);
+        ucell.atoms[0].ncpp.index1_soc[3] = std::vector<int>(5, 0);
+        ucell.atoms[0].ncpp.index2_soc[3] = std::vector<int>(5, 0);
         for (int i = 0; i < 5; ++i)
         {
             ucell.atoms[0].ncpp.d_real(i, i) = 1.0;
@@ -93,10 +93,6 @@ class NonlocalNewTest : public ::testing::Test
         delete[] ucell.atoms[0].iw2l;
         delete[] ucell.atoms[0].iw2m;
         delete[] ucell.atoms[0].iw2n;
-        delete[] ucell.atoms[0].ncpp.index1_soc[0];
-        delete[] ucell.atoms[0].ncpp.index2_soc[0];
-        delete[] ucell.atoms[0].ncpp.index1_soc[3];
-        delete[] ucell.atoms[0].ncpp.index2_soc[3];
         delete[] ucell.atoms;
         delete[] ucell.iat2it;
         delete[] ucell.iat2ia;
diff --git a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp
index d484ea9b2d..e49512ccaa 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/VL_in_pw.cpp
@@ -47,9 +47,9 @@ void pseudopot_cell_vl::init_vloc(ModuleBase::matrix& vloc_in, const ModulePW::P
 		{
 			this->vloc_of_g(
 					atom->ncpp.msh, // after cutoff 
-					atom->ncpp.rab,
-		          	atom->ncpp.r, 
-					atom->ncpp.vloc_at, // local potential in real space radial form.  
+					atom->ncpp.rab.data(),
+		          	atom->ncpp.r.data(), 
+					atom->ncpp.vloc_at.data(), // local potential in real space radial form.  
 		          	this->zp[it],
 					vloc1d,
 					rho_basis);
diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp
index 37a008a9a3..def8276a8e 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_grad_pw.cpp
@@ -39,7 +39,7 @@ void pseudopot_cell_vnl::initgradq_vnl(const UnitCell &cell)
 			for (int iq=0; iq<GlobalV::NQX; iq++)  
 			{
 				const double q = iq * GlobalV::DQ;
-				ModuleBase::Sphbes::dSpherical_Bessel_dx(kkbeta, cell.atoms[it].ncpp.r, q, l, djl);
+				ModuleBase::Sphbes::dSpherical_Bessel_dx(kkbeta, cell.atoms[it].ncpp.r.data(), q, l, djl);
 
 				for (int ir = 0;ir < kkbeta;ir++)
 				{
@@ -47,7 +47,7 @@ void pseudopot_cell_vnl::initgradq_vnl(const UnitCell &cell)
 					          djl[ir] * pow(cell.atoms[it].ncpp.r[ir],2);
 				} 
 				double vqint;
-				ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab, vqint);
+				ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab.data(), vqint);
 				this->tab_dq(it, ib, iq) = vqint * pref;
 			} 
 		} 
diff --git a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp
index 32d94adbe8..04e2b62b98 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/VNL_in_pw.cpp
@@ -20,15 +20,17 @@ pseudopot_cell_vnl::pseudopot_cell_vnl()
 
 pseudopot_cell_vnl::~pseudopot_cell_vnl()
 {
-    if (GlobalV::use_paw)
+    if (GlobalV::use_paw) {
         return;
+}
     delete[] indv_ijkb0;
 }
 
 void pseudopot_cell_vnl::release_memory()
 {
-    if (this->nhm <= 0 || memory_released)
+    if (this->nhm <= 0 || memory_released) {
         return;
+}
     if (GlobalV::device_flag == "gpu")
     {
         if (GlobalV::precision_flag == "single")
@@ -84,8 +86,9 @@ void pseudopot_cell_vnl::init(const int ntype,
                               const ModulePW::PW_Basis_K* wfc_basis,
                               const bool allocate_vkb)
 {
-    if (GlobalV::use_paw)
+    if (GlobalV::use_paw) {
         return;
+}
 
     ModuleBase::TITLE("pseudopot_cell_vnl", "init");
     ModuleBase::timer::tick("ppcell_vnl", "init");
@@ -300,10 +303,12 @@ void pseudopot_cell_vnl::init(const int ntype,
 //----------------------------------------------------------
 void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in) const
 {
-    if (GlobalV::use_paw)
+    if (GlobalV::use_paw) {
         return;
-    if (GlobalV::test_pp)
+}
+    if (GlobalV::test_pp) {
         ModuleBase::TITLE("pseudopot_cell_vnl", "getvnl");
+}
     ModuleBase::timer::tick("pp_cell_vnl", "getvnl");
 
     if (lmaxkb < 0)
@@ -344,13 +349,15 @@ void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in
         const int nbeta = GlobalC::ucell.atoms[it].ncpp.nbeta;
         const int nh = GlobalC::ucell.atoms[it].ncpp.nh;
 
-        if (GlobalV::test_pp > 1)
+        if (GlobalV::test_pp > 1) {
             ModuleBase::GlobalFunc::OUT("nbeta", nbeta);
+}
 
         for (int nb = 0; nb < nbeta; nb++)
         {
-            if (GlobalV::test_pp > 1)
+            if (GlobalV::test_pp > 1) {
                 ModuleBase::GlobalFunc::OUT("ib", nb);
+}
             for (int ig = 0; ig < npw; ig++)
             {
                 const double gnorm = gk[ig].norm() * GlobalC::ucell.tpiba;
@@ -406,10 +413,12 @@ void pseudopot_cell_vnl::getvnl(const int& ik, ModuleBase::ComplexMatrix& vkb_in
 template <typename FPTYPE, typename Device>
 void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex<FPTYPE>* vkb_in) const
 {
-    if (GlobalV::use_paw)
+    if (GlobalV::use_paw) {
         return;
-    if (GlobalV::test_pp)
+}
+    if (GlobalV::test_pp) {
         ModuleBase::TITLE("pseudopot_cell_vnl", "getvnl");
+}
     ModuleBase::timer::tick("pp_cell_vnl", "getvnl");
 
     using cal_vnl_op = hamilt::cal_vnl_op<FPTYPE, Device>;
@@ -532,8 +541,9 @@ void pseudopot_cell_vnl::getvnl(Device* ctx, const int& ik, std::complex<FPTYPE>
 
 void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_basis)
 {
-    if (GlobalV::use_paw)
+    if (GlobalV::use_paw) {
         return;
+}
     ModuleBase::TITLE("pseudopot_cell_vnl", "init_vnl");
     ModuleBase::timer::tick("ppcell_vnl", "init_vnl");
 
@@ -555,8 +565,9 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_
         if (cell.atoms[it].ncpp.tvanp)
         {
             cell.atoms[it].ncpp.nqlc = std::min(cell.atoms[it].ncpp.nqlc, lmaxq);
-            if (cell.atoms[it].ncpp.nqlc < 0)
+            if (cell.atoms[it].ncpp.nqlc < 0) {
                 cell.atoms[it].ncpp.nqlc = 0;
+}
         }
     }
 
@@ -673,14 +684,15 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_
                             this->dvan_so(ijs, it, ip, ip2)
                                 = cell.atoms[it].ncpp.dion(ir, is) * soc.fcoef(it, is1, is2, ip, ip2);
                             ++ijs;
-                            if (ir != is)
+                            if (ir != is) {
                                 soc.fcoef(it, is1, is2, ip, ip2) = std::complex<double>(0.0, 0.0);
+}
                         }
                     }
                 }
             }
         }
-        else
+        else {
             for (int ip = 0; ip < Nprojectors; ip++)
             {
                 for (int ip2 = 0; ip2 < Nprojectors; ip2++)
@@ -701,6 +713,7 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_
                     }
                 }
             }
+}
     }
 
     // e) It computes the coefficients c_{LM}^{nm} which relates the
@@ -835,14 +848,14 @@ void pseudopot_cell_vnl::init_vnl(UnitCell& cell, const ModulePW::PW_Basis* rho_
             for (int iq = 0; iq < GlobalV::NQX; iq++)
             {
                 const double q = iq * GlobalV::DQ;
-                ModuleBase::Sphbes::Spherical_Bessel(kkbeta, cell.atoms[it].ncpp.r, q, l, jl);
+                ModuleBase::Sphbes::Spherical_Bessel(kkbeta, cell.atoms[it].ncpp.r.data(), q, l, jl);
 
                 for (int ir = 0; ir < kkbeta; ir++)
                 {
                     aux[ir] = cell.atoms[it].ncpp.betar(ib, ir) * jl[ir] * cell.atoms[it].ncpp.r[ir];
                 }
                 double vqint;
-                ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab, vqint);
+                ModuleBase::Integral::Simpson_Integral(kkbeta, aux, cell.atoms[it].ncpp.rab.data(), vqint);
                 this->tab(it, ib, iq) = vqint * pref;
             }
         }
@@ -915,7 +928,7 @@ void pseudopot_cell_vnl::compute_qrad(UnitCell& cell)
                 {
                     const double q = iq * GlobalV::DQ;
                     // here we compute the spherical bessel function for each q_i
-                    ModuleBase::Sphbes::Spherical_Bessel(kkbeta, upf->r, q, l, besr);
+                    ModuleBase::Sphbes::Spherical_Bessel(kkbeta, upf->r.data(), q, l, besr);
                     for (int nb = 0; nb < nbeta; nb++)
                     {
                         // the Q are symmetric with respect to indices nb and mb
@@ -931,7 +944,7 @@ void pseudopot_cell_vnl::compute_qrad(UnitCell& cell)
                                 }
                                 // then we integrate with all the Q functions
                                 double vqint;
-                                ModuleBase::Integral::Simpson_Integral(kkbeta, aux, upf->rab, vqint);
+                                ModuleBase::Integral::Simpson_Integral(kkbeta, aux, upf->rab.data(), vqint);
                                 qrad(it, l, ijv, iq) = vqint * pref;
                             }
                         }
@@ -1279,10 +1292,11 @@ double pseudopot_cell_vnl::CG(int l1, int m1, int l2, int m2, int L, int M) // p
 // }
 #endif
 
-void pseudopot_cell_vnl::init_vnl_alpha(void) // pengfei Li 2018-3-23
+void pseudopot_cell_vnl::init_vnl_alpha() // pengfei Li 2018-3-23
 {
-    if (GlobalV::test_pp)
+    if (GlobalV::test_pp) {
         ModuleBase::TITLE("pseudopot_cell_vnl", "init_vnl_alpha");
+}
     ModuleBase::timer::tick("ppcell_vnl", "init_vnl_alpha");
 
     for (int it = 0; it < GlobalC::ucell.ntype; it++)
@@ -1331,7 +1345,7 @@ void pseudopot_cell_vnl::init_vnl_alpha(void) // pengfei Li 2018-3-23
                 for (int iq = 0; iq < GlobalV::NQX; iq++)
                 {
                     const double q = iq * GlobalV::DQ;
-                    ModuleBase::Sphbes::Spherical_Bessel(kkbeta, GlobalC::ucell.atoms[it].ncpp.r, q, L, jl);
+                    ModuleBase::Sphbes::Spherical_Bessel(kkbeta, GlobalC::ucell.atoms[it].ncpp.r.data(), q, L, jl);
 
                     for (int ir = 0; ir < kkbeta; ir++)
                     {
@@ -1339,7 +1353,7 @@ void pseudopot_cell_vnl::init_vnl_alpha(void) // pengfei Li 2018-3-23
                                   * GlobalC::ucell.atoms[it].ncpp.r[ir] * GlobalC::ucell.atoms[it].ncpp.r[ir];
                     }
                     double vqint;
-                    ModuleBase::Integral::Simpson_Integral(kkbeta, aux, GlobalC::ucell.atoms[it].ncpp.rab, vqint);
+                    ModuleBase::Integral::Simpson_Integral(kkbeta, aux, GlobalC::ucell.atoms[it].ncpp.rab.data(), vqint);
                     this->tab_alpha(it, ib, L, iq) = vqint * pref;
                 }
             }
@@ -1362,8 +1376,9 @@ void pseudopot_cell_vnl::cal_effective_D(const ModuleBase::matrix& veff,
                                          const ModulePW::PW_Basis* rho_basis,
                                          UnitCell& cell)
 {
-    if (GlobalV::use_paw)
+    if (GlobalV::use_paw) {
         return;
+}
     ModuleBase::TITLE("pseudopot_cell_vnl", "cal_effective_D");
 
     /*
diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp
index 7c6db88fd9..c111a78938 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp
@@ -65,8 +65,9 @@ void Forces<FPTYPE, Device>::cal_force_cc(ModuleBase::matrix& forcecc,
     }
     else
     {
-        if (GlobalV::NSPIN == 4)
+        if (GlobalV::NSPIN == 4) {
             GlobalC::ucell.cal_ux();
+}
         const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &GlobalC::ucell);
 
         // etxc = std::get<0>(etxc_vtxc_v);
@@ -117,9 +118,9 @@ void Forces<FPTYPE, Device>::cal_force_cc(ModuleBase::matrix& forcecc,
             //                                 rhocg);
             this->deriv_drhoc(GlobalC::ppcell.numeric,
                               GlobalC::ucell.atoms[it].ncpp.msh,
-                              GlobalC::ucell.atoms[it].ncpp.r,
-                              GlobalC::ucell.atoms[it].ncpp.rab,
-                              GlobalC::ucell.atoms[it].ncpp.rho_atc,
+                              GlobalC::ucell.atoms[it].ncpp.r.data(),
+                              GlobalC::ucell.atoms[it].ncpp.rab.data(),
+                              GlobalC::ucell.atoms[it].ncpp.rho_atc.data(),
                               rhocg,
                               rho_basis,
                               1);
diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp
index 33b953e5ea..aef02eaa2d 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/forces_scc.cpp
@@ -77,8 +77,9 @@ void Forces<FPTYPE, Device>::cal_force_scc(ModuleBase::matrix& forcescc,
 
     int igg0 = 0;
     const int ig0 = rho_basis->ig_gge0;
-    if (rho_basis->gg_uniq[0] < 1.0e-8)
+    if (rho_basis->gg_uniq[0] < 1.0e-8) {
         igg0 = 1;
+}
 
     double fact = 2.0;
     for (int nt = 0; nt < ucell_in.ntype; nt++) {
@@ -86,9 +87,9 @@ void Forces<FPTYPE, Device>::cal_force_scc(ModuleBase::matrix& forcescc,
         const int mesh = ucell_in.atoms[nt].ncpp.msh;
         this->deriv_drhoc_scc(GlobalC::ppcell.numeric,
                             mesh,
-                            ucell_in.atoms[nt].ncpp.r,
-                            ucell_in.atoms[nt].ncpp.rab,
-                            ucell_in.atoms[nt].ncpp.rho_at,
+                            ucell_in.atoms[nt].ncpp.r.data(),
+                            ucell_in.atoms[nt].ncpp.rab.data(),
+                            ucell_in.atoms[nt].ncpp.rho_at.data(),
                             rhocgnt.data(),
                             rho_basis,
                             ucell_in);        
@@ -105,8 +106,9 @@ void Forces<FPTYPE, Device>::cal_force_scc(ModuleBase::matrix& forcescc,
 #pragma omp parallel for reduction(+ : force0) reduction(+ : force1) reduction(+ : force2)
 #endif
                     for (int ig = 0; ig < rho_basis->npw; ++ig) {
-                        if (ig == ig0)
+                        if (ig == ig0) {
                             continue;
+}
                         const ModuleBase::Vector3<double> gv
                             = rho_basis->gcar[ig];
                         const double rhocgntigg
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp
index 3b0a0a34d0..b92e7a6f3c 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp
@@ -104,9 +104,9 @@ void Stress_Func<FPTYPE, Device>::stress_cc(ModuleBase::matrix& sigma,
 			this->deriv_drhoc(
 				GlobalC::ppcell.numeric,
 				GlobalC::ucell.atoms[nt].ncpp.msh,
-				GlobalC::ucell.atoms[nt].ncpp.r,
-				GlobalC::ucell.atoms[nt].ncpp.rab,
-				GlobalC::ucell.atoms[nt].ncpp.rho_atc,
+				GlobalC::ucell.atoms[nt].ncpp.r.data(),
+				GlobalC::ucell.atoms[nt].ncpp.rab.data(),
+				GlobalC::ucell.atoms[nt].ncpp.rho_atc.data(),
 				rhocg,
 				rho_basis,
 				1);
@@ -129,9 +129,9 @@ void Stress_Func<FPTYPE, Device>::stress_cc(ModuleBase::matrix& sigma,
 			this->deriv_drhoc (
 				GlobalC::ppcell.numeric,
 				GlobalC::ucell.atoms[nt].ncpp.msh,
-				GlobalC::ucell.atoms[nt].ncpp.r,
-				GlobalC::ucell.atoms[nt].ncpp.rab,
-				GlobalC::ucell.atoms[nt].ncpp.rho_atc,
+				GlobalC::ucell.atoms[nt].ncpp.r.data(),
+				GlobalC::ucell.atoms[nt].ncpp.rab.data(),
+				GlobalC::ucell.atoms[nt].ncpp.rho_atc.data(),
 				rhocg,
 				rho_basis,
 				0);
diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp
index b54060135a..2b340e2dc1 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_loc.cpp
@@ -95,8 +95,8 @@ void Stress_Func<FPTYPE, Device>::stress_loc(ModuleBase::matrix& sigma,
 		//
 		// normal case: dvloc contains dV_loc(G)/dG
 		//
-			this->dvloc_of_g ( atom->ncpp.msh, atom->ncpp.rab, atom->ncpp.r,
-					atom->ncpp.vloc_at, atom->ncpp.zv, dvloc.data(), rho_basis);
+			this->dvloc_of_g ( atom->ncpp.msh, atom->ncpp.rab.data(), atom->ncpp.r.data(),
+					atom->ncpp.vloc_at.data(), atom->ncpp.zv, dvloc.data(), rho_basis);
 		//
 		}
 #ifndef _OPENMP
diff --git a/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp b/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp
index e08c12a9ec..64967bf64e 100644
--- a/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp
+++ b/source/module_hamilt_pw/hamilt_pwdft/wf_atomic.cpp
@@ -100,7 +100,7 @@ void WF_atomic::init_at_1(Structure_Factor *sf_in)
                 inner_part[ir] = atom->ncpp.chi(ic,ir) * atom->ncpp.chi(ic,ir);
             }
             double unit = 0.0;
-            ModuleBase::Integral::Simpson_Integral(nmesh, inner_part, atom->ncpp.rab, unit);
+            ModuleBase::Integral::Simpson_Integral(nmesh, inner_part, atom->ncpp.rab.data(), unit);
             delete[] inner_part;
 
             // liuyu add 2023-10-06
@@ -150,7 +150,7 @@ void WF_atomic::init_at_1(Structure_Factor *sf_in)
                         {
                             norm_beta[ik] = atom->ncpp.betar(ib, ik) * atom->ncpp.chi(ic, ik);
                         }
-                        ModuleBase::Integral::Simpson_Integral(kkbeta, norm_beta, atom->ncpp.rab, work[ib]);
+                        ModuleBase::Integral::Simpson_Integral(kkbeta, norm_beta, atom->ncpp.rab.data(), work[ib]);
                     }
                     else
                     {
@@ -188,14 +188,14 @@ void WF_atomic::init_at_1(Structure_Factor *sf_in)
                 for (int iq=startq; iq<GlobalV::NQX; iq++)
                 {
                     const double q = GlobalV::DQ * iq;
-                    ModuleBase::Sphbes::Spherical_Bessel(atom->ncpp.msh, atom->ncpp.r, q, l, aux);
+                    ModuleBase::Sphbes::Spherical_Bessel(atom->ncpp.msh, atom->ncpp.r.data(), q, l, aux);
                     for (int ir = 0;ir < atom->ncpp.msh;ir++)
                     {
                         vchi[ir] = atom->ncpp.chi(ic,ir) * aux[ir] * atom->ncpp.r[ir];
                     }
 
                     double vqint = 0.0;
-                    ModuleBase::Integral::Simpson_Integral(atom->ncpp.msh, vchi, atom->ncpp.rab, vqint);
+                    ModuleBase::Integral::Simpson_Integral(atom->ncpp.msh, vchi, atom->ncpp.rab.data(), vqint);
 
                     GlobalC::ppcell.tab_at(it, ic, iq) =  vqint * pref;
                     //				if( it == 0 && ic == 0 )
diff --git a/source/module_psi/psi_initializer_atomic.cpp b/source/module_psi/psi_initializer_atomic.cpp
index ac06497d56..b89fe3e62f 100644
--- a/source/module_psi/psi_initializer_atomic.cpp
+++ b/source/module_psi/psi_initializer_atomic.cpp
@@ -123,12 +123,12 @@ void psi_initializer_atomic<T, Device>::tabulate()
             int n_rgrid = (PARAM.inp.pseudo_mesh)?atom->ncpp.mesh:atom->ncpp.msh;
             std::vector<double> pswfcr(n_rgrid);
             for (int ir=0; ir<n_rgrid; ir++) pswfcr[ir] = atom->ncpp.chi(ic, ir);
-            normalize(n_rgrid, pswfcr, atom->ncpp.rab);
+            normalize(n_rgrid, pswfcr, atom->ncpp.rab.data());
             if (atom->ncpp.oc[ic] >= 0.0) // reasonable occupation number, but is it always true?
             {
                 const int l = atom->ncpp.lchi[ic];
                 std::vector<double> ovlp_pswfcjlq_q(GlobalV::NQX);
-                this->sbt.direct(l, atom->ncpp.msh, atom->ncpp.r, pswfcr.data(), GlobalV::NQX, qgrid.data(), ovlp_pswfcjlq_q.data(), 1);
+                this->sbt.direct(l, atom->ncpp.msh, atom->ncpp.r.data(), pswfcr.data(), GlobalV::NQX, qgrid.data(), ovlp_pswfcjlq_q.data(), 1);
                 for (int iq = 0; iq < GlobalV::NQX; iq++)
                 {
                     this->ovlp_pswfcjlq_(it, ic, iq) = pref * ovlp_pswfcjlq_q[iq];
diff --git a/source/module_psi/test/psi_initializer_unit_test.cpp b/source/module_psi/test/psi_initializer_unit_test.cpp
index 8ae213ab35..98288db9f4 100644
--- a/source/module_psi/test/psi_initializer_unit_test.cpp
+++ b/source/module_psi/test/psi_initializer_unit_test.cpp
@@ -83,7 +83,8 @@ std::complex<double>* Structure_Factor::get_sk(int ik, int it, int ia, ModulePW:
 {
     int npw = wfc_basis->npwk[ik];
     std::complex<double> *sk = new std::complex<double>[npw];
-    for(int ipw = 0; ipw < npw; ++ipw) sk[ipw] = std::complex<double>(0.0, 0.0);
+    for(int ipw = 0; ipw < npw; ++ipw) { sk[ipw] = std::complex<double>(0.0, 0.0);
+}
     return sk;
 }
 
@@ -137,7 +138,8 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_ucell->tpiba = 2.0 * M_PI / this->p_ucell->lat0;
             this->p_ucell->tpiba2 = this->p_ucell->tpiba * this->p_ucell->tpiba;
             // atom
-            if(this->p_ucell->atom_label != nullptr) delete[] this->p_ucell->atom_label;
+            if(this->p_ucell->atom_label != nullptr) { delete[] this->p_ucell->atom_label;
+}
             this->p_ucell->atom_label = new std::string[1];
             this->p_ucell->atom_label[0] = "Si";
             // atom properties
@@ -156,7 +158,8 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_ucell->atoms[0].taud[0] = {0.25, 0.25, 0.25};
             this->p_ucell->atoms[0].mbl[0] = {0, 0, 0};
             // atom pseudopotential
-            if(this->p_ucell->pseudo_fn != nullptr) delete[] this->p_ucell->pseudo_fn;
+            if(this->p_ucell->pseudo_fn != nullptr) { delete[] this->p_ucell->pseudo_fn;
+}
             this->p_ucell->pseudo_fn = new std::string[1];
             this->p_ucell->pseudo_fn[0] = "Si_NCSR_ONCVPSP_v0.5_dojo.upf";
             this->p_ucell->natomwfc = 4;
@@ -165,29 +168,34 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_ucell->atoms[0].ncpp.msh = 10;
             this->p_ucell->atoms[0].ncpp.lmax = 2;
             //if(this->p_ucell->atoms[0].ncpp.rab != nullptr) delete[] this->p_ucell->atoms[0].ncpp.rab;
-            this->p_ucell->atoms[0].ncpp.rab = new double[10];
-            for(int i = 0; i < 10; ++i) this->p_ucell->atoms[0].ncpp.rab[i] = 0.01;
+            this->p_ucell->atoms[0].ncpp.rab = std::vector<double>(10, 0.0);
+            for(int i = 0; i < 10; ++i) { this->p_ucell->atoms[0].ncpp.rab[i] = 0.01;
+}
             //if(this->p_ucell->atoms[0].ncpp.r != nullptr) delete[] this->p_ucell->atoms[0].ncpp.r;
-            this->p_ucell->atoms[0].ncpp.r = new double[10];
-            for(int i = 0; i < 10; ++i) this->p_ucell->atoms[0].ncpp.r[i] = 0.01*i;
+            this->p_ucell->atoms[0].ncpp.r = std::vector<double>(10, 0.0);
+            for(int i = 0; i < 10; ++i) { this->p_ucell->atoms[0].ncpp.r[i] = 0.01*i;
+}
             this->p_ucell->atoms[0].ncpp.chi.create(2, 10);
-            for(int i = 0; i < 2; ++i) for(int j = 0; j < 10; ++j) this->p_ucell->atoms[0].ncpp.chi(i, j) = 0.01;
+            for(int i = 0; i < 2; ++i) { for(int j = 0; j < 10; ++j) { this->p_ucell->atoms[0].ncpp.chi(i, j) = 0.01;
+}
+}
             //if(this->p_ucell->atoms[0].ncpp.lchi != nullptr) delete[] this->p_ucell->atoms[0].ncpp.lchi;
-            this->p_ucell->atoms[0].ncpp.lchi = new int[2];
+            this->p_ucell->atoms[0].ncpp.lchi = std::vector<int>(2, 0);
             this->p_ucell->atoms[0].ncpp.lchi[0] = 0;
             this->p_ucell->atoms[0].ncpp.lchi[1] = 1;
             this->p_ucell->lmax_ppwf = 1;
-            this->p_ucell->atoms[0].ncpp.oc = new double[2];
+            this->p_ucell->atoms[0].ncpp.oc = std::vector<double>(2, 0.0);
             this->p_ucell->atoms[0].ncpp.oc[0] = 1.0;
             this->p_ucell->atoms[0].ncpp.oc[1] = 1.0;
 
             this->p_ucell->atoms[0].ncpp.has_so = false;
-            this->p_ucell->atoms[0].ncpp.jchi = new double[2];
+            this->p_ucell->atoms[0].ncpp.jchi = std::vector<double>(2, 0.0);
             this->p_ucell->atoms[0].ncpp.jchi[0] = 0.5;
             this->p_ucell->atoms[0].ncpp.jchi[1] = 1.5;
             // atom numerical orbital
             this->p_ucell->lmax = 2;
-            if(this->p_ucell->orbital_fn != nullptr) delete[] this->p_ucell->orbital_fn;
+            if(this->p_ucell->orbital_fn != nullptr) { delete[] this->p_ucell->orbital_fn;
+}
             this->p_ucell->orbital_fn = new std::string[1];
             this->p_ucell->orbital_fn[0] = "Si_gga_8au_60Ry_2s2p1d.orb";
             this->p_ucell->atoms[0].nwl = 2;
@@ -201,29 +209,36 @@ class PsiIntializerUnitTest : public ::testing::Test {
             // can support function PW_Basis::getfftixy2is
             this->p_pw_wfc->nks = 1;
             this->p_pw_wfc->npwk_max = 1;
-            if(this->p_pw_wfc->npwk != nullptr) delete[] this->p_pw_wfc->npwk;
+            if(this->p_pw_wfc->npwk != nullptr) { delete[] this->p_pw_wfc->npwk;
+}
             this->p_pw_wfc->npwk = new int[1];
             this->p_pw_wfc->npwk[0] = 1;
             this->p_pw_wfc->fftnxy = 1;
             this->p_pw_wfc->fftnz = 1;
             this->p_pw_wfc->nst = 1;
             this->p_pw_wfc->nz = 1;
-            if(this->p_pw_wfc->is2fftixy != nullptr) delete[] this->p_pw_wfc->is2fftixy;
+            if(this->p_pw_wfc->is2fftixy != nullptr) { delete[] this->p_pw_wfc->is2fftixy;
+}
             this->p_pw_wfc->is2fftixy = new int[1];
             this->p_pw_wfc->is2fftixy[0] = 0;
-            if(this->p_pw_wfc->fftixy2ip != nullptr) delete[] this->p_pw_wfc->fftixy2ip;
+            if(this->p_pw_wfc->fftixy2ip != nullptr) { delete[] this->p_pw_wfc->fftixy2ip;
+}
             this->p_pw_wfc->fftixy2ip = new int[1];
             this->p_pw_wfc->fftixy2ip[0] = 0;
-            if(this->p_pw_wfc->igl2isz_k != nullptr) delete[] this->p_pw_wfc->igl2isz_k;
+            if(this->p_pw_wfc->igl2isz_k != nullptr) { delete[] this->p_pw_wfc->igl2isz_k;
+}
             this->p_pw_wfc->igl2isz_k = new int[1];
             this->p_pw_wfc->igl2isz_k[0] = 0;
-            if(this->p_pw_wfc->gcar != nullptr) delete[] this->p_pw_wfc->gcar;
+            if(this->p_pw_wfc->gcar != nullptr) { delete[] this->p_pw_wfc->gcar;
+}
             this->p_pw_wfc->gcar = new ModuleBase::Vector3<double>[1];
             this->p_pw_wfc->gcar[0] = {0.0, 0.0, 0.0};
-            if(this->p_pw_wfc->igl2isz_k != nullptr) delete[] this->p_pw_wfc->igl2isz_k;
+            if(this->p_pw_wfc->igl2isz_k != nullptr) { delete[] this->p_pw_wfc->igl2isz_k;
+}
             this->p_pw_wfc->igl2isz_k = new int[1];
             this->p_pw_wfc->igl2isz_k[0] = 0;
-            if(this->p_pw_wfc->gk2 != nullptr) delete[] this->p_pw_wfc->gk2;
+            if(this->p_pw_wfc->gk2 != nullptr) { delete[] this->p_pw_wfc->gk2;
+}
             this->p_pw_wfc->gk2 = new double[1];
             this->p_pw_wfc->gk2[0] = 0.0;
             this->p_pw_wfc->latvec.e11 = this->p_ucell->latvec.e11; this->p_pw_wfc->latvec.e12 = this->p_ucell->latvec.e12; this->p_pw_wfc->latvec.e13 = this->p_ucell->latvec.e13;
@@ -235,10 +250,12 @@ class PsiIntializerUnitTest : public ::testing::Test {
             this->p_pw_wfc->lat0 = this->p_ucell->lat0;
             this->p_pw_wfc->tpiba = 2.0 * M_PI / this->p_ucell->lat0;
             this->p_pw_wfc->tpiba2 = this->p_pw_wfc->tpiba * this->p_pw_wfc->tpiba;
-            if(this->p_pw_wfc->kvec_c != nullptr) delete[] this->p_pw_wfc->kvec_c;
+            if(this->p_pw_wfc->kvec_c != nullptr) { delete[] this->p_pw_wfc->kvec_c;
+}
             this->p_pw_wfc->kvec_c = new ModuleBase::Vector3<double>[1];
             this->p_pw_wfc->kvec_c[0] = {0.0, 0.0, 0.0};
-            if(this->p_pw_wfc->kvec_d != nullptr) delete[] this->p_pw_wfc->kvec_d;
+            if(this->p_pw_wfc->kvec_d != nullptr) { delete[] this->p_pw_wfc->kvec_d;
+}
             this->p_pw_wfc->kvec_d = new ModuleBase::Vector3<double>[1];
             this->p_pw_wfc->kvec_d[0] = {0.0, 0.0, 0.0};