diff --git a/ABACUS.develop/source/Makefile.Objects b/ABACUS.develop/source/Makefile.Objects index 10825db039..caef758e51 100644 --- a/ABACUS.develop/source/Makefile.Objects +++ b/ABACUS.develop/source/Makefile.Objects @@ -30,11 +30,7 @@ dftd3_subroutine.o\ vdwd3.o\ vdwd3_parameters.o\ pseudopot_upf.o\ -pseudo_h.o \ -pseudo_atom.o \ -pseudo_vl.o \ pseudo_nc.o \ -pseudo_us.o \ pseudopot_cell_vl.o \ pseudopot_cell_vnl.o \ atom_pseudo.o\ diff --git a/ABACUS.develop/source/src_io/read_pseudopot.cpp b/ABACUS.develop/source/src_io/read_pseudopot.cpp index 645fdc3781..8ddb8c1b47 100644 --- a/ABACUS.develop/source/src_io/read_pseudopot.cpp +++ b/ABACUS.develop/source/src_io/read_pseudopot.cpp @@ -10,22 +10,19 @@ void UnitCell_pseudo::read_pseudopot(const string &pp_dir) { TITLE("UnitCell_pseudo","read_pseudopot"); -//---------------------------------------------------------- -// EXPLAIN : setup reading log for pseudopot_upf -//---------------------------------------------------------- + // setup reading log for pseudopot_upf stringstream ss; ss << global_out_dir << "atom_pseudo.log"; -//---------------------------------------------------------- -// EXPLAIN : Read in the atomic pseudo potential -//---------------------------------------------------------- + // Read in the atomic pseudo potentials string pp_address; for (int i = 0;i < ntype;i++) { Pseudopot_upf upf; // mohan update 2010-09-12 - int error = 0, error_ap = 0; + int error = 0; + int error_ap = 0; if(MY_RANK==0) { @@ -48,7 +45,10 @@ void UnitCell_pseudo::read_pseudopot(const string &pp_dir) Parallel_Common::bcast_int(error_ap); #endif - if(error_ap) WARNING_QUIT("UnitCell_pseudo::read_pseudopot","error when average the pseudopotential."); + if(error_ap) + { + WARNING_QUIT("UnitCell_pseudo::read_pseudopot","error when average the pseudopotential."); + } if(error==1) { @@ -64,7 +64,6 @@ void UnitCell_pseudo::read_pseudopot(const string &pp_dir) { WARNING_QUIT("read_pseudopot","Check the reference states in pseudopotential .vwr file.\n Also the norm of the read in pseudo wave functions\n explicitly please check S, P and D channels.\n If the norm of the wave function is \n unreasonable large (should be near 1.0), ABACUS would quit. \n The solution is to turn off the wave functions \n and the corresponding non-local projectors together\n in .vwr pseudopotential file."); } -// OUT(ofs_running,"PP_ERRROR",error); //xiaohui add 2015-03-24 #ifdef __MPI @@ -79,7 +78,7 @@ void UnitCell_pseudo::read_pseudopot(const string &pp_dir) if(MY_RANK==0) { // upf.print_pseudo_upf( ofs ); - atoms[i].set_pseudo_us( upf ); + atoms[i].set_pseudo_nc( upf ); ofs_running << "\n Read in pseudopotential file is " << pseudo_fn[i] << endl; OUT(ofs_running,"pseudopotential type",atoms[i].pp_type); diff --git a/ABACUS.develop/source/src_pdiag/pdiag_double.cpp b/ABACUS.develop/source/src_pdiag/pdiag_double.cpp index df36a8d0ec..4000479693 100644 --- a/ABACUS.develop/source/src_pdiag/pdiag_double.cpp +++ b/ABACUS.develop/source/src_pdiag/pdiag_double.cpp @@ -24,7 +24,16 @@ extern "C" #include "src_external/src_test/test_function.h" -inline int cart2blacs(MPI_Comm comm_2D, int nprows, int npcols, int N, int nblk, int lld, int *desc, int &mpi_comm_rows, int &mpi_comm_cols) +inline int cart2blacs( + MPI_Comm comm_2D, + int nprows, + int npcols, + int N, + int nblk, + int lld, + int *desc, + int &mpi_comm_rows, + int &mpi_comm_cols) { #ifdef __MPI int my_blacs_ctxt; @@ -53,10 +62,18 @@ inline int cart2blacs(MPI_Comm comm_2D, int nprows, int npcols, int N, int nblk, #endif } -inline int q2ZLOC_WFC(int pos, int naroc[2], int nb, - int dim0, int dim1, int iprow, int ipcol, - int loc_size, - double* work, double* ZLOC, double** WFC) +inline int q2ZLOC_WFC( + int pos, + int naroc[2], + int nb, + int dim0, + int dim1, + int iprow, + int ipcol, + int loc_size, + double* work, + double* ZLOC, + double** WFC) { //OUT(ofs_running,"start q2ZLOC_WFC"); for(int j=0; j* work, complex** WFC) +inline int q2WFC_complex( + int naroc[2], + int nb, + int dim0, + int dim1, + int iprow, + int ipcol, + complex* work, + complex** WFC) { for(int j=0; j* work, complex** WFC, complex** WFCAUG) +inline int q2WFC_WFCAUG_complex( + int naroc[2], + int nb, + int dim0, + int dim1, + int iprow, + int ipcol, + complex* work, + complex** WFC, + complex** WFCAUG) { for(int j=0; j* work, complex** WFC, complex** CTOT) +inline int q2WFC_CTOT_complex( + int myid, + int naroc[2], + int nb, + int dim0, + int dim1, + int iprow, + int ipcol, + complex* work, + complex** WFC, + complex** CTOT) { for(int j=0; j* work, complex** WFC, complex** WFCAUG, complex** CTOT) +inline int q2WFC_WFCAUG_CTOT_complex( + int myid, + int naroc[2], + int nb, + int dim0, + int dim1, + int iprow, + int ipcol, + complex* work, + complex** WFC, + complex** WFCAUG, + complex** CTOT) { for(int j=0; j **wfc, ComplexMatrix &wfc_2d, - complex* ch_mat, complex* cs_mat, double *ekb) +void Pdiag_Double::diago_complex_begin( + const int &ik, + complex **wfc, + ComplexMatrix &wfc_2d, + complex* ch_mat, + complex* cs_mat, + double *ekb) { #ifdef TEST_DIAG { @@ -804,12 +888,18 @@ void Pdiag_Double::diago_complex_begin(const int &ik, complex **wfc, Com if(std::norm(m[index])>1E-10) { if(std::imag(m[index])>1E-10) + { ofs< **wfc, Com if(std::norm(m[index])>1E-10) { if(std::imag(m[index])>1E-10) + { ofs< **wfc, Com &NLOCAL, h_tmp.c, &one, &one, desc, s_tmp.c, &one, &one, desc, NULL, NULL, &il, &iu, &abstol, &M, &NZ, ekb, &orfac, wfc_2d.c, &one, &one, desc, - work.data(), &lwork, rwork.data(), &lrwork, iwork.data(), &liwork, ifail.data(), iclustr.data(), gap.data(), &info); + work.data(), &lwork, rwork.data(), &lrwork, + iwork.data(), &liwork, ifail.data(), iclustr.data(), gap.data(), &info); + ofs_running<<"lwork="< **wfc, Com #ifdef __MPI -void Pdiag_Double::readin(const string &fa, const string &fb, const int &nlocal_tot, double *eigen, double *eigvr) +void Pdiag_Double::readin( + const string &fa, + const string &fb, + const int &nlocal_tot, + double *eigen, + double *eigvr) { TITLE("Pdiag_Double","readin"); @@ -1186,4 +1289,3 @@ void Pdiag_Double::readin(const string &fa, const string &fb, const int &nlocal_ delete[] Z; } #endif - diff --git a/ABACUS.develop/source/src_pw/atom_pseudo.cpp b/ABACUS.develop/source/src_pw/atom_pseudo.cpp index 64258bb5c6..b7c8458232 100644 --- a/ABACUS.develop/source/src_pw/atom_pseudo.cpp +++ b/ABACUS.develop/source/src_pw/atom_pseudo.cpp @@ -4,7 +4,6 @@ Atom_pseudo::Atom_pseudo() { - if(test_atom) TITLE("atom_pseudo","Constructor"); mbl = new Vector3[1]; pseudo_fn = "not_init"; mass = 0.0; @@ -28,7 +27,7 @@ void Atom_pseudo::print_atom(ofstream &ofs) #ifdef __MPI void Atom_pseudo::bcast_atom_pseudo(const int &na) { - if(test_pp) TITLE("Atom_pseudo","bcast_atom_pseudo"); + TITLE("Atom_pseudo","bcast_atom_pseudo"); Parallel_Common::bcast_double( mass ); Parallel_Common::bcast_string( pseudo_fn ); @@ -48,7 +47,7 @@ void Atom_pseudo::bcast_atom_pseudo(const int &na) void Atom_pseudo::bcast_atom_pseudo2(void) { - if(test_pp) TITLE("Atom_pseudo","bcast_atom_pseudo2"); + TITLE("Atom_pseudo","bcast_atom_pseudo2"); // == pseudo_h == //int Parallel_Common::bcast_int( lmax ); @@ -142,57 +141,14 @@ void Atom_pseudo::bcast_atom_pseudo2(void) betar.create(nr,nc); dion.create(nbeta, nbeta); } + + // below two 'bcast_double' lines of codes seem to have bugs, + // on some computers, the code will stuck here for ever. + // mohan note 2021-04-28 Parallel_Common::bcast_double( dion.c , nbeta * nbeta); Parallel_Common::bcast_double( betar.c, nr * nc ); // == end of psesudo_nc == -// == pseudo_us == -/* - Parallel_Common::bcast_int( nqf ); - Parallel_Common::bcast_int( nqlc ); - if(MY_RANK!=0) - { - rinner = new double[ nqlc ]; - } - Parallel_Common::bcast_double( rinner, nqlc); - - if(nbeta > 0) - { - if(MY_RANK!=0) - { - qqq.create(nbeta, nbeta); - } - Parallel_Common::bcast_double( qqq.c, nbeta*nbeta); - } - - if(mesh>0 && nbeta>0) - { - if(MY_RANK!=0) - { - qfunc.create(nbeta, nbeta, mesh); - } - Parallel_Common::bcast_double(qfunc.ptr, qfunc.getSize() ); - } - - int b[4]; - if(MY_RANK==0) - { - b[0] = qfcoef.getBound1(); - b[1] = qfcoef.getBound2(); - b[2] = qfcoef.getBound3(); - b[3] = qfcoef.getBound4(); - } - Parallel_Common::bcast_int( b, 4); - - if(MY_RANK!=0) - { - qfcoef.create(b[0],b[1],b[2],b[3]); - } - - Parallel_Common::bcast_double(qfcoef.ptr, qfcoef.getSize() ); -*/ -// == end of pseudo_us == - return; } diff --git a/ABACUS.develop/source/src_pw/atom_pseudo.h b/ABACUS.develop/source/src_pw/atom_pseudo.h index f66cc4f323..91b2f15d57 100644 --- a/ABACUS.develop/source/src_pw/atom_pseudo.h +++ b/ABACUS.develop/source/src_pw/atom_pseudo.h @@ -1,29 +1,21 @@ -//========================================================== -// AUTHOR : mohan -// DATE : 2008-11-08 -//========================================================== #ifndef ATOM_PSEUDO_H #define ATOM_PSEUDO_H #include "tools.h" -#include "pseudo_us.h" +#include "pseudo_nc.h" using namespace std; -//========================================================== -// CLASS : -// NAME : Atom_pseudo -//========================================================== -class Atom_pseudo : public pseudo_us +class Atom_pseudo : public pseudo_nc { public: Atom_pseudo(); ~Atom_pseudo(); - Vector3 *mbl; //If this atom can move - string pseudo_fn;// File name of pseudopotentia + Vector3 *mbl; // whether the atoms can move or not + string pseudo_fn; // File name of pseudopotentials double mass; // the mass of atom - bool flag_empty_element = false; // whether is the empty element for bsse. Peize Lin add 2021.04.07 + bool flag_empty_element = false; // whether is the empty element for bsse. Peize Lin add 2021.04.07 protected: @@ -31,7 +23,7 @@ class Atom_pseudo : public pseudo_us #ifdef __MPI void bcast_atom_pseudo(const int &na); - void bcast_atom_pseudo2(void); + void bcast_atom_pseudo2(void); // for upf201 #endif }; diff --git a/ABACUS.develop/source/src_pw/hamilt_pw.cpp b/ABACUS.develop/source/src_pw/hamilt_pw.cpp index 2a1181ace7..9e43540936 100644 --- a/ABACUS.develop/source/src_pw/hamilt_pw.cpp +++ b/ABACUS.develop/source/src_pw/hamilt_pw.cpp @@ -654,62 +654,62 @@ void Hamilt_PW::add_vuspsi(complex *hpsi_in,const complex *becp, // is calculated before. if(NSPIN!=4) { - for (int it=0; it becp1,becp2; - const int Nprojs = ucell.atoms[it].nh; - for (int ia=0; ia becp1,becp2; + const int Nprojs = ucell.atoms[it].nh; + for (int ia=0; ia rcut) - { - msh = ir + 1; - br = true; - break; - // goto 5; - } // endif - } // enddo - - if (br) - { - // force msh to be odd for simpson integration - msh = 2 * (int)((msh + 1) / 2) - 1; // 5 - } - else - { - msh = mesh ; - } -// cout << " msh=" << msh << endl; -} // end subroutine set_pseudo - -void pseudo_atom::print_pseudo_at(ofstream &ofs) -{ - print_pseudo_h(ofs); - ofs << "\n pseudo_atom : "; - ofs << "\n msh " << msh; -// ofs << "\n nchi " << nchi; - out.printr1_d(ofs, " r : ", r, mesh); - out.printr1_d(ofs, " rab : ", rab, mesh); - out.printr1_d(ofs, " rho_atc : ", rho_atc, mesh); - out.printr1_d(ofs, " rho_at : ", rho_at, mesh); - out.printr1_d(ofs," jchi : ", jchi, nchi); - out.printrm(ofs, " chi : ", chi); - ofs << "\n ----------------------"; -} diff --git a/ABACUS.develop/source/src_pw/pseudo_atom.h b/ABACUS.develop/source/src_pw/pseudo_atom.h deleted file mode 100644 index 470dfc910f..0000000000 --- a/ABACUS.develop/source/src_pw/pseudo_atom.h +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef PSEUDO_ATOM_H -#define PSEUDO_ATOM_H - -#include "tools.h" -#include "pseudopot_upf.h" -#include "pseudo_h.h" - -class pseudo_atom: public pseudo_h -{ - public: - - // - double *r; // radial logaritmic mesh, r[0:mesh-1] - double *rab; // derivative of the radial mesh, rab[0:mesh-1] - - // - double *rho_atc; // radial core charge density, rho_atc[0:mesh-1] - - // - double *rho_at; // radial atomic charge density, rho_at[0:mesh-1] - - // - matrix chi; // radial atomic orbitals, chi(nchi, mesh) - - //other - int msh; // the point at rcut - double rcut; // cut-off radius - - //member functions - pseudo_atom(); - ~pseudo_atom(); - - void set_pseudo_at(const Pseudopot_upf &upf); - void print_pseudo_at(ofstream &ofs); -}; - -#endif diff --git a/ABACUS.develop/source/src_pw/pseudo_h.cpp b/ABACUS.develop/source/src_pw/pseudo_h.cpp deleted file mode 100644 index 885c063829..0000000000 --- a/ABACUS.develop/source/src_pw/pseudo_h.cpp +++ /dev/null @@ -1,164 +0,0 @@ -//========================================================== -// AUTHOR : Lixin He,mohan -// DATE : 2008-11-11 -//========================================================== -#include "xc_type.h" -#include "global.h" -#include "pseudo_h.h" -#include "tools.h" - -pseudo_h::pseudo_h() -{ - els = new string[1]; - lchi = new int[1]; - oc = new double[1]; - jjj = new double[1]; - jchi = new double[1]; - nn = new int[1]; - has_so = false; - zv = 0; -} - -pseudo_h::~pseudo_h() -{ - delete[] els; - delete[] lchi; - delete[] oc; - delete[] jjj; - delete[] jchi; - delete[] nn; -} - -//--------------------------------------------------------------------- -void pseudo_h::set_pseudo_h(const Pseudopot_upf &upf) -{ - // TITLE("pseudo_h::set_pseudo_h"); - // set "is"-th pseudopotential using the Unified Pseudopotential Format - - int i=0; - this->nv = upf.nv;// UPF file version number - this->psd = upf.psd; - this->pp_type = upf.pp_type; - - this->tvanp = upf.tvanp;// if USPP - this->nlcc = upf.nlcc;// Non linear core corrections( bool ?) - - for(int i=0; i<4; i++) - { - this->dft[i] = upf.dft[i]; - } - - this->zv = upf.zp; - this->etotps = upf.etotps; - this->ecutwfc = upf.ecutwfc; - this->ecutrho = upf.ecutrho; - - this->lmax = upf.lmax; - this->mesh = upf.mesh; - - // mohan update 2021-02-22 - // max number of points in the atomic radial mesh - int ndmx = 2000; - if (this->mesh > ndmx) - { - cout << "\n set_pseudo_h, too many grid points,"; - } - - this->nchi = upf.nwfc; - this->nbeta = upf.nbeta; - - delete[] els; - this->els = new string[nchi]; - assert(this->els != 0); - - delete[] lchi; - this->lchi = new int[this->nchi]; - assert(this->lchi != 0); - - delete[] oc; - this->oc = new double[nchi]; - assert(this->oc != 0); - - for (i = 0;i < nchi;i++) - { - this->els[i] = upf.els[i]; - this->lchi[i] = upf.lchi[i]; - this->oc[i] = upf.oc[i]; - } - - delete[] jjj; - this->jjj = new double[nbeta]; - assert(this->jjj != 0); - - delete[] nn; - this->nn = new int[nchi]; - assert(this->nn != 0); - - delete[] jchi; - this->jchi = new double[nchi]; - assert(this->jchi != 0); - - this->has_so = upf.has_so;//added by zhengdy-soc - if (this->has_so) - { - for (i = 0;i < nchi;i++){ - this->nn[i] = upf.nn[i]; - this->jchi[i] = upf.jchi[i]; - } - for (i = 0;i < upf.nbeta;i++) - { - this->jjj[i] = upf.jjj [i]; - } - } - else - { - for (i = 0;i < nchi;i++){ - this->nn[i] = 0; - this->jchi[i] = 0; - } - for (i = 0;i < upf.nbeta;i++) - { - this->jjj[i] = 0; - } - } - - /* - jchi = new double[nwfc]; - if (upf.has_so){ - for(i=0;i -#include -#include - -using namespace std; - -#include "pseudopot_upf.h" - -class pseudo_h -{ -public: - - pseudo_h(); - ~pseudo_h(); - - // - bool has_so; // if .true. includes spin-orbit - int nv; // UPF file version number - string psd; // Element label - string pp_type; // Pseudo type ( NC or US ) - bool tvanp; // .true. if Ultrasoft - bool nlcc; // Non linear core corrections(bool) - string dft[4]; // Exch-Corr type - int zv; // z valence - double etotps; // total energy - double ecutwfc; // suggested cut-off for wfc - double ecutrho; // suggested cut-off for rho - int lmax; // maximum angular momentum component - int mesh; // number of point in the radial mesh - int nchi; // nwfc,number of wavefunctions - int nbeta; // number of projectors - string *els; // els[nchi] - int *lchi; // lchi[nchi] - double *oc; // oc[nchi] - - double *jjj; // total angual momentum, jjj[nbeta] - double *jchi; //jchi(nwfc), added by zhengdy-soc - int *nn; - - // member functions - void set_pseudo_h(const Pseudopot_upf &upf); - void print_pseudo_h(ofstream &ofs); -}; - -#endif // PSEUDOH_H diff --git a/ABACUS.develop/source/src_pw/pseudo_nc.cpp b/ABACUS.develop/source/src_pw/pseudo_nc.cpp index 1c338d6130..36a8787693 100644 --- a/ABACUS.develop/source/src_pw/pseudo_nc.cpp +++ b/ABACUS.develop/source/src_pw/pseudo_nc.cpp @@ -1,40 +1,79 @@ -/* pseudo_nc.cpp */ #include "pseudo_nc.h" #include "global.h" pseudo_nc::pseudo_nc() { -// cout << "\n === pseudo_nc === "; +// pseudo_h + els = new string[1]; + lchi = new int[1]; + oc = new double[1]; + jjj = new double[1]; + jchi = new double[1]; + nn = new int[1]; + has_so = false; + zv = 0; + +// pseudo local parts + vloc_at = new double[1]; + +// pseudo_atom + r = new double[1]; + rab = new double[1]; + rho_at = new double[1]; + rho_atc = new double[1]; + +// pseudo_nc lll = new int[1]; } pseudo_nc::~pseudo_nc() { -// cout << "\n exit pseudo_nc "; +// pseudo_h + delete[] els; + delete[] lchi; + delete[] oc; + delete[] jjj; + delete[] jchi; + delete[] nn; + +// pseudo local parts + delete[] vloc_at; + +// pseudo_atom + delete[] r; + delete[] rab; + delete[] rho_at; + delete[] rho_atc; + +// pseudo_nc delete[] lll; } + + //--------------------------------------------------------------------- void pseudo_nc::set_pseudo_nc(const Pseudopot_upf &upf) { - //----------------------------------------------------------------- - - int i, nb; + TITLE("pseudo_nc","set_pseudo_nc"); + // call subroutines + this->set_pseudo_h(upf); + this->set_pseudo_atom(upf); this->set_pseudo_vl(upf); + delete[] lll; lll = new int[nbeta]; assert(lll != 0); - for (i = 0;i < nbeta;i++) + for (int i = 0;i < nbeta;i++) { lll[i] = upf.lll[i]; } kkbeta = 0; - for (nb = 0;nb < nbeta;nb++) + for (int nb = 0;nb < nbeta;nb++) { kkbeta = (upf.kkbeta[nb] > kkbeta) ? upf.kkbeta[nb] : kkbeta; } @@ -52,7 +91,7 @@ void pseudo_nc::set_pseudo_nc(const Pseudopot_upf &upf) nh = 0; - for (nb = 0; nb < nbeta;nb++) + for (int nb = 0; nb < nbeta;nb++) { nh += 2 * lll [nb] + 1; } @@ -73,3 +112,284 @@ void pseudo_nc::print_pseudo_nc(ofstream &ofs) ofs << "\n ----------------------"; } + +void pseudo_nc::set_pseudo_h(const Pseudopot_upf &upf) +{ + TITLE("pseudo_nc::set_pseudo_h"); + // set pseudopotential for each atom type + // by using the Unified Pseudopotential Format + + this->nv = upf.nv;// UPF file version number + this->psd = upf.psd; + this->pp_type = upf.pp_type; + + this->tvanp = upf.tvanp;// if USPP + this->nlcc = upf.nlcc;// Non linear core corrections( bool ?) + + for(int i=0; i<4; i++) + { + this->dft[i] = upf.dft[i]; + } + + this->zv = upf.zp; + this->etotps = upf.etotps; + this->ecutwfc = upf.ecutwfc; + this->ecutrho = upf.ecutrho; + + this->lmax = upf.lmax; + this->mesh = upf.mesh; + + // mohan update 2021-02-22 + // max number of points in the atomic radial mesh + int ndmx = 2000; + if (this->mesh > ndmx) + { + cout << "\n set_pseudo_h, too many grid points,"; + } + + this->nchi = upf.nwfc; + this->nbeta = upf.nbeta; + + delete[] els; + this->els = new string[nchi]; + assert(this->els != 0); + + delete[] lchi; + this->lchi = new int[this->nchi]; + assert(this->lchi != 0); + + delete[] oc; + this->oc = new double[nchi]; + assert(this->oc != 0); + + for (int i = 0;i < nchi;i++) + { + this->els[i] = upf.els[i]; + this->lchi[i] = upf.lchi[i]; + this->oc[i] = upf.oc[i]; + } + + delete[] jjj; + this->jjj = new double[nbeta]; + assert(this->jjj != 0); + + delete[] nn; + this->nn = new int[nchi]; + assert(this->nn != 0); + + delete[] jchi; + this->jchi = new double[nchi]; + assert(this->jchi != 0); + + this->has_so = upf.has_so;//added by zhengdy-soc + + if (this->has_so) + { + for(int i=0;i < nchi;i++) + { + this->nn[i] = upf.nn[i]; + this->jchi[i] = upf.jchi[i]; + } + for(int i=0;i < upf.nbeta;i++) + { + this->jjj[i] = upf.jjj [i]; + } + } + else + { + for (int i=0; inn[i] = 0; + this->jchi[i] = 0; + } + for (int i=0; ijjj[i] = 0; + } + } + + /* + jchi = new double[nwfc]; + if (upf.has_so){ + for(i=0;ircut = 15.0;//(a.u.); + +// remember to update here if you need it. +// rcut = 25.0; + + OUT(ofs_running,"PAO radial cut off (Bohr)",rcut); + if(rcut <= 0.0) + { + WARNING_QUIT("pseudo_atom::set_pseudo_atom","PAO rcut<=0.0"); + } + + chi.create(nchi, mesh); + + delete[] r; + r = new double[mesh]; + assert(r != 0); + ZEROS(r, mesh); + + delete[] rab; + rab = new double[mesh]; + assert(rab != 0); + ZEROS(rab, mesh); + + delete[] rho_at; + rho_at = new double[mesh]; + assert(rho_at != 0); + ZEROS(rho_at,mesh); + + delete[] rho_atc; + rho_atc = new double[mesh]; + assert(rho_atc != 0); + ZEROS(rho_atc, mesh); + + for (int i = 0;i < nchi;i++) + { + for (int j = 0; j < mesh; j++) + { + chi(i, j) = upf.chi(i, j); + } + } + + for (int i = 0;i < mesh;i++) + { + r [i] = upf.r [i]; + rab[i] = upf.rab[i]; + rho_at [i] = upf.rho_at [i]; + } + + // nlcc: non-linear core corrections + // rho_atc: core atomic charge + if (nlcc) + { + for (int i = 0;i < mesh;i++) + { + rho_atc[i] = upf.rho_atc[i]; + } + } + else + { + for (int i = 0;i < upf.mesh;i++) + { + rho_atc[i] = 0.0; + } + } // end if + + + bool br = false; + + this->msh = 0; + + for (int ir = 0;ir < mesh;ir++) + { + if (r [ir] > rcut) + { + msh = ir + 1; + br = true; + break; + } + } + + if (br) + { + // force msh to be odd for simpson integration + msh = 2 * (int)((msh + 1) / 2) - 1; // 5 + } + else + { + msh = mesh ; + } + + return; +} // end subroutine set_pseudo + + + +void pseudo_nc::set_pseudo_vl(const Pseudopot_upf &upf) +{ + TITLE("pseudo_nc","set_pseudo_vl"); + + assert(mesh>0);//mohan add 2021-05-01 + + 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]; + } + + return; +} + + +void pseudo_nc::print_pseudo_atom(ofstream &ofs) +{ + print_pseudo_h(ofs); + ofs << "\n pseudo_atom : "; + ofs << "\n msh " << msh; +// ofs << "\n nchi " << nchi; + out.printr1_d(ofs, " r : ", r, mesh); + out.printr1_d(ofs, " rab : ", rab, mesh); + out.printr1_d(ofs, " rho_atc : ", rho_atc, mesh); + out.printr1_d(ofs, " rho_at : ", rho_at, mesh); + out.printr1_d(ofs," jchi : ", jchi, nchi); + out.printrm(ofs, " chi : ", chi); + ofs << "\n ----------------------"; +} + + +void pseudo_nc::print_pseudo_vl(ofstream &ofs) +{ + ofs << "\n pseudo_vl:"; + print_pseudo_atom(ofs); + out.printr1_d(ofs, "vloc_at : ", vloc_at, mesh); + ofs << "\n ----------------------------------- "; +} + +void pseudo_nc::print_pseudo_h(ofstream &ofs) +{ + ofs << "\n pseudo_info :"; + ofs << "\n nv " << nv; + ofs << "\n psd " << psd; + ofs << "\n pp_type " << pp_type; + ofs << "\n tvanp " << tvanp; + ofs << "\n nlcc " << nlcc; + ofs << "\n dft " << dft; + ofs << "\n zv " << zv; + ofs << "\n etotps " << etotps; + ofs << "\n ecutwfc " << ecutwfc; + ofs << "\n ecutrho " << ecutrho; + ofs << "\n lmax " << lmax; + ofs << "\n mesh " << mesh; + ofs << "\n nchi " << nchi; + ofs << "\n nbeta " << nbeta; +// out.printr1_d(ofs," els: ", els, nchi); + out.printr1_d(ofs, " lchi: ", lchi, nchi); + out.printr1_d(ofs, " oc: ", oc, nchi); + ofs << "\n ----------------------"; +} + diff --git a/ABACUS.develop/source/src_pw/pseudo_nc.h b/ABACUS.develop/source/src_pw/pseudo_nc.h index 80b8f4f61d..4f1869f05b 100644 --- a/ABACUS.develop/source/src_pw/pseudo_nc.h +++ b/ABACUS.develop/source/src_pw/pseudo_nc.h @@ -3,27 +3,82 @@ #include "tools.h" #include "pseudopot_upf.h" -#include "pseudo_vl.h" -class pseudo_nc: public pseudo_vl +//----------------------------------------- +// read in norm conserving pseudopotentials +// mohan update 2021-05-01 +//----------------------------------------- +class pseudo_nc { public: pseudo_nc(); ~pseudo_nc(); + // + bool has_so; // if .true. includes spin-orbit + int nv; // UPF file version number + string psd; // Element label + string pp_type; // Pseudo type ( NC or US ) + bool tvanp; // .true. if Ultrasoft + bool nlcc; // Non linear core corrections(bool) + string dft[4]; // Exch-Corr type + int zv; // z valence + double etotps; // total energy + double ecutwfc; // suggested cut-off for wfc + double ecutrho; // suggested cut-off for rho + int lmax; // maximum angular momentum component + int mesh; // number of point in the radial mesh + int nchi; // nwfc,number of wavefunctions + int nbeta; // number of projectors + string *els; // els[nchi] + int *lchi; // lchi[nchi] + double *oc; // oc[nchi] + + double *jjj; // total angual momentum, jjj[nbeta] + double *jchi; //jchi(nwfc), added by zhengdy-soc + int *nn; + + // Local pseudopotentials + double *vloc_at; // [mesh], local potential( = pseudopot_upf.vloc ) + + // + double *r; // radial logaritmic mesh, r[0:mesh-1] + double *rab; // derivative of the radial mesh, rab[0:mesh-1] + + // + double *rho_atc; // radial core charge density, rho_atc[0:mesh-1] + + // + double *rho_at; // radial atomic charge density, rho_at[0:mesh-1] + + // + matrix chi; // radial atomic orbitals, chi(nchi, mesh) + + //other + int msh; // number of points up to rcut + double rcut; // cut-off radius + + // int *lll; // lll(nbeta), angular momentum of the beta function int kkbeta; // kkbeta(nbeta), point where the beta are zero // matrix dion; // dion(nbeta,nbeta) - matrix betar; //(nbeta, mesh), radial beta_{mu} functions + matrix betar; // (nbeta, mesh), radial beta_{mu} functions // other int nh; // number of beta functions per atomic type + void set_pseudo_h(const Pseudopot_upf &upf); + void set_pseudo_atom(const Pseudopot_upf &upf); + void set_pseudo_vl(const Pseudopot_upf &upf); void set_pseudo_nc(const Pseudopot_upf &upf); + + void print_pseudo_h(ofstream &ofs); + void print_pseudo_atom(ofstream &ofs); + void print_pseudo_vl(ofstream &ofs); void print_pseudo_nc(ofstream &ofs); }; diff --git a/ABACUS.develop/source/src_pw/pseudo_us.cpp b/ABACUS.develop/source/src_pw/pseudo_us.cpp deleted file mode 100644 index 2bdabf1615..0000000000 --- a/ABACUS.develop/source/src_pw/pseudo_us.cpp +++ /dev/null @@ -1,73 +0,0 @@ -#include "../src_io/output.h" -#include "pseudo_us.h" - -pseudo_us::pseudo_us() -{ -// cout << "\n === pseudo_us === "; - rinner = new double[1]; -} - -pseudo_us::~pseudo_us() -{ -// cout << "\n exit pseudo_us "; - delete[] rinner; -} - -//--------------------------------------------------------------------- -void pseudo_us::set_pseudo_us(const Pseudopot_upf &upf) -{ -// cout << "\n === set_pseudo_us() === "; - //int i; - - this->set_pseudo_nc(upf); - -/* nqlc = upf.nqlc; - nqf = upf.nqf; - - delete[] rinner; - rinner = new double[nqlc]; - - for (i = 0;i < nqlc;i++) - { - this->rinner[i] = upf.rinner[i]; - } - - if (nbeta > 0) - { - this->qqq.create(nbeta, nbeta); - this->qqq = upf.qqq; - } - - if (mesh > 0 && nbeta > 0) - { - this->qfunc.create(nbeta, nbeta, mesh); - this->qfunc = upf.qfunc; - } - - this->qfcoef.create(nbeta, nbeta, upf.qfcoef.getBound3(), upf.qfcoef.getBound4() ); - this->qfcoef = upf.qfcoef; -*/ - -// if (lspinorb && !upf.has_so) // see USE spin_orb -// cout << "\n upf_to_internal,At least one non s.o. pseudo,"; //-1); -// -// lspinorb=lspinorb && upf.has_so; -// cout << "\n End set_pseudo_us() " << endl; -} // end subroutine set_pseudo_upf - -void pseudo_us::print_pseudo_us(ofstream &ofs) -{ - ofs << "\n\n internal format for pseudopotential "; - output out; - print_pseudo_nc(ofs); - ofs << "\n pseudo_us : "; - ofs << "\n nqf " << nqf - << "\n nqlc " << nqlc; - out.printr1_d(ofs, " rinner : ", rinner, nbeta); - out.printrm(ofs, " qqq : ", qqq); - out.printr3_d(ofs, " qfunc : ", qfunc); - out.printr4_d(ofs, " qfcoef : ", qfcoef); - ofs << "\n ----------------------\n"; -} - - diff --git a/ABACUS.develop/source/src_pw/pseudo_us.h b/ABACUS.develop/source/src_pw/pseudo_us.h deleted file mode 100644 index 3ee83b0d9d..0000000000 --- a/ABACUS.develop/source/src_pw/pseudo_us.h +++ /dev/null @@ -1,35 +0,0 @@ -#ifndef PSEUDO_US_H -#define PSEUDO_US_H - -#include "tools.h" -#include "pseudopot_upf.h" -#include "pseudo_nc.h" - -class pseudo_us: public pseudo_nc -{ - public: - - pseudo_us(); - ~pseudo_us(); - - // - int nqf; - - // - double *rinner; // rinner(0:2*lmax) - matrix qqq; // qqq(nbeta,nbeta) - realArray qfunc; // qfunc(mesh,nbeta,nbeta) - - // - realArray qfcoef; // qfcoef(nqf,0:2*lmax,nbeta,nbeta) - - // other ? - int nqlc; // number of angular momenta in Q - - void set_pseudo_us(const Pseudopot_upf &upf); - - void print_pseudo_us(ofstream &ofs); - -}; - -#endif //PSEUDO_US_H diff --git a/ABACUS.develop/source/src_pw/pseudo_vl.h b/ABACUS.develop/source/src_pw/pseudo_vl.h index 458cf14084..6bf761f881 100644 --- a/ABACUS.develop/source/src_pw/pseudo_vl.h +++ b/ABACUS.develop/source/src_pw/pseudo_vl.h @@ -8,6 +8,7 @@ class pseudo_vl: public pseudo_atom { public: + double *vloc_at; // [mesh], local potential( = pseudopot_upf.vloc ) pseudo_vl(); diff --git a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp b/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp index 67af7a0a6f..a492bb00e5 100644 --- a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp +++ b/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp @@ -68,15 +68,15 @@ void pseudopot_cell_vnl::init(const int ntype, const bool allocate_vkb) if( this->nhm > 0 ) { - indv.create(ntype, this->nhm); - nhtol.create(ntype, this->nhm); - nhtolm.create(ntype, this->nhm); - nhtoj.create(ntype, this->nhm); - deeq.create(NSPIN, ucell.nat, this->nhm, this->nhm); - deeq_nc.create(NSPIN, ucell.nat, this->nhm, this->nhm); - dvan.create(ntype, this->nhm, this->nhm); - dvan_so.create(NSPIN, ntype, this->nhm, this->nhm); - becsum.create(NSPIN, ucell.nat, this->nhm * (this->nhm + 1) / 2); + this->indv.create(ntype, this->nhm); + this->nhtol.create(ntype, this->nhm); + this->nhtolm.create(ntype, this->nhm); + this->nhtoj.create(ntype, this->nhm); + this->deeq.create(NSPIN, ucell.nat, this->nhm, this->nhm); + this->deeq_nc.create(NSPIN, ucell.nat, this->nhm, this->nhm); + this->dvan.create(ntype, this->nhm, this->nhm); + this->dvan_so.create(NSPIN, ntype, this->nhm, this->nhm); + this->becsum.create(NSPIN, ucell.nat, this->nhm * (this->nhm + 1) / 2); } else { @@ -127,18 +127,6 @@ void pseudopot_cell_vnl::init(const int ntype, const bool allocate_vkb) this->tab_at.create(ntype, nchix_nc, NQX); } - if(test_pp > 1) - { - cout - << "\n ntype = " << ntype - << "\n lmaxkb = " << this->lmaxkb - << "\n nhm = " << this->nhm - << "\n nkb = " << this->nkb - << "\n nbrx = " << nbrx - << "\n nchix = " << nchix - << endl; - } - timer::tick("ppcell_vnl","init",'C'); return; } diff --git a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h b/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h index a5636811b0..d9cd3c82d0 100644 --- a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h +++ b/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.h @@ -32,11 +32,9 @@ class pseudopot_cell_vnl: public pseudopot_cell_vl double cell_factor; //LiuXh add 20180619 -// NAME : nkb(total number of beta functions, with struct.fact.) - int nkb; // be called in hm.hpw.init + int nkb; // total number of beta functions considering all atoms -// NAME : lmaxkb(max angular momentum,(see pseudo_h)) - int lmaxkb; + int lmaxkb; // max angular momentum for non-local projectors void init_vnl(UnitCell_pseudo &cell); @@ -60,6 +58,7 @@ class pseudopot_cell_vnl: public pseudopot_cell_vl int calculate_nqx(const double &ecutwfc,const double &dq); int nhm; + int lmaxq; matrix indv; // indes linking atomic beta's to beta's in the solid @@ -78,13 +77,16 @@ class pseudopot_cell_vnl: public pseudopot_cell_vl ComplexArray deeq_nc; //(:,:,:,:), the spin-orbit case realArray becsum; //(:,:,:,:), \sum_i f(i) //used in charge + ComplexMatrix vkb; // all beta functions in reciprocal space complex ***vkb1_alpha; complex ***vkb_alpha; - + // other variables complex Cal_C(int alpha, int lu, int mu, int L, int M); + double CG(int l1, int m1, int l2, int m2, int L, int M); + void print_vnl(ofstream &ofs); ORB_gaunt_table MGT; diff --git a/ABACUS.develop/source/src_pw/pseudopot_upf.cpp b/ABACUS.develop/source/src_pw/pseudopot_upf.cpp index 001aa65b29..55db8ffd93 100644 --- a/ABACUS.develop/source/src_pw/pseudopot_upf.cpp +++ b/ABACUS.develop/source/src_pw/pseudopot_upf.cpp @@ -1,7 +1,3 @@ -//========================================================== -// Author: Lixin He,mohan -// DATE : 2009-02-26 -//========================================================== #include "../src_io/output.h" #include "pseudopot_upf.h" #include @@ -65,13 +61,12 @@ int Pseudopot_upf::init_pseudo_reader(const string &fn) // First check if this pseudo-potential has spin-orbit information ifstream ifs(fn.c_str(), ios::in); - if (!ifs) + // can't find the file. + if (!ifs) { - // can't find the file. return 1; } - //cout << "global_pseudo_type =" << global_pseudo_type << endl; if(global_pseudo_type=="auto") //zws { set_pseudo_type(fn); @@ -81,20 +76,17 @@ int Pseudopot_upf::init_pseudo_reader(const string &fn) if(global_pseudo_type=="upf") { int info = read_pseudo_upf(ifs); - //this->print_pseudo_upf(ofs_running); return info; } // read in the .vwr type of pseudopotentials else if(global_pseudo_type=="vwr") { int info = read_pseudo_vwr(ifs); - // this->print_pseudo_upf(ofs_running); return info; } else if(global_pseudo_type=="upf201") { int info = read_pseudo_upf201(ifs); - //this->print_pseudo_upf(ofs_running); return info; } @@ -108,22 +100,16 @@ int Pseudopot_upf::init_pseudo_reader(const string &fn) int Pseudopot_upf::set_pseudo_type(const string &fn) //zws add { ifstream pptype_ifs(fn.c_str(), ios::in); - string dummy, strversion; - //int ierr = 0; + string dummy; + string strversion; if (pptype_ifs.good()) { getline(pptype_ifs,dummy); - //cout << "dummy : " << dummy << endl; - //if ( dummy == "" ) - //{ - // global_pseudo_type = "upf"; - //} stringstream wdsstream(dummy); getline(wdsstream,strversion,'"'); getline(wdsstream,strversion,'"'); - //cout << "strversi:" << strversion << endl; if ( trim(strversion) == "2.0.1" ) { @@ -142,7 +128,9 @@ string& Pseudopot_upf::trim(string &in_str) static const string deltri = " \t" ; // delete tab or space string::size_type position = in_str.find_first_of(deltri, position); if (position == string::npos) + { return in_str; + } return trim(in_str.erase(position, 1) ); } @@ -150,20 +138,18 @@ string Pseudopot_upf::trimend(string &in_str) { const string &deltri =" \t" ; string::size_type position = in_str.find_last_not_of(deltri)+1; - //cout << "--" << position << "--" << endl; string tmpstr=in_str.erase(position); - //cout << "--" << tmpstr << "--" << endl; return tmpstr.erase(0,tmpstr.find_first_not_of(deltri)); -} //}zws +} //zws //---------------------------------------------------------- // This code is used to read in vwr pseudopotential format, -// which is used in PEtot mostly. Now I 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 +// 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 @@ -501,11 +487,10 @@ int Pseudopot_upf::read_pseudo_vwr(ifstream &ifs) { coef=coef+(vl[ir]-vloc[ir])*wlr[ir]*wlr[ir]*(r[ir+1]-r[ir-1])/2.0; } - // ofs_running << setw(15) << r[ir] << setw(25) << beta(ib,ir) << endl; } -// In quantum espresso they did this: +// In pw they did this: // dion(ib,ib)=1.0/coef; if(coef<0.0) dion(ib,ib) = -1.0; @@ -673,6 +658,7 @@ int Pseudopot_upf::read_pseudo_upf(ifstream &ifs) return 0; }// end subroutine read_pseudopot_upf + void Pseudopot_upf::read_pseudo_header(ifstream &ifs) { READ_VALUE(ifs, this->nv);// Version number @@ -855,18 +841,12 @@ void Pseudopot_upf::read_pseudo_nl(ifstream &ifs) if (nbeta == 0) { -// this->nqf = 0; -// this->nqlc = 0; delete[] kkbeta; delete[] lll; this->kkbeta = new int[1]; this->lll = new int[1]; this->beta.create(1, 1); this->dion.create(1, 1); -// this->rinner = new double[lmax]; -// this->qqq.create(1, 1); -// this->qfunc.create(1, 1, mesh); -// this->qfcoef.create(1, 1, 1, lmax); return; } else @@ -921,93 +901,9 @@ void Pseudopot_upf::read_pseudo_nl(ifstream &ifs) if (tvanp) { WARNING_QUIT("read_pseudo_nl","Ultra Soft Pseudopotential not available yet."); - /* - SCAN_BEGIN(ifs, "", false); - ifs >> this->nqf;//nl_6 - int xx = (this->nqf > 1) ? this->nqf : 1; - ifs.ignore(75, '\n'); - this->nqlc = 2 * lmax + 1; - this->rinner = new double[nqlc](); - - this->qqq.create(nbeta , nbeta); - this->qqq.zero_out(); - - this->qfunc.create(nbeta , nbeta , mesh); - this->qfcoef.create(nbeta , nbeta , nqlc , xx); - - if (nqf != 0) - { - SCAN_BEGIN(ifs, "", false); - - for (i = 0;i < nqlc;i++) - { - ifs >> idum >> rinner[i]; - } - - SCAN_END(ifs, ""); - } - - for (int nb = 0;nb < nbeta;nb++) - { - for (int mb = nb;mb < nbeta;mb++) - { - ifs >> idum >> idum >> ldum; - ifs.ignore(75, '\n'); - - if (ldum != lll[mb]) - { - cerr << " read_pseudo_nl inconsistent angular momentum for Q_ij \n "; - } - - ifs >> x; - - this->qqq(nb, mb) = x; - ifs.ignore(75, '\n'); - - // "Q_int" - // qqq.assign(mb , nb , x); - this->qqq(mb, nb) = x; - - for (n = 0;n < mesh;n++) - { - ifs >> x; - this->qfunc(mb , nb , n) = x ; - } - - for (n = 0;n < mesh;n++) - { - x = this->qfunc(mb,nb,n); // bad idea! - this->qfunc(nb,mb,n) = x; - } - - // QFCOEF - if (nqf > 0) - { - SCAN_BEGIN(ifs, "", false); - for (lp = 0;lp < nqlc;lp++) - { - for (i = 1;i < nqf;i++) - { - ifs >> x; - this->qfcoef(nb , mb , lp , i) = x; - } - } - SCAN_END(ifs, ""); - } - } - } - SCAN_END(ifs, ""); - */ } else // not tvanp { -// this->nqf = 1;//nl_6 -// this->nqlc = 2 * lmax + 1;//nl_7 -// this->rinner = new double[nqlc]();//nl_8 -// this->qqq.create(nbeta , nbeta);//nl_9 -// this->qqq.zero_out(); -// this->qfunc.create(nbeta , nbeta , mesh);//nl_10 -// this->qfcoef.create(nbeta , nbeta , nqlc , nqf);//nl_11 } } return; @@ -1079,46 +975,6 @@ void Pseudopot_upf::read_pseudo_so(ifstream &ifs) return; } -/* -// -// This routine reads from the new UPF file, -// and the total angual momentum jjj of the beta and jchi of the -// wave-functions. -void Pseudopot_upf::read_pseudo_addinfo(ifstream &ifs) -{ - for (int i = 0;i < 6;i++) - { - ifs.ignore(75, '\n'); - } - - int nb = 0; - - this->nn = new int[nwfc](); - this->rcut = new double[nwfc](); - this->rcutus = new double[nwfc](); - this->epseu = new double[nwfc](); - this->jchi = new double[nwfc](); - this->jjj = new double[nbeta](); - - char s[10]; - - for (nb = 0; nb < nwfc;nb++) - { - ifs >> s >> this->nn[nb] >> this->lchi[nb] >> this->jchi[nb] >> this->oc[nb]; - this->els[nb] = s; - } - - - for (nb = 0; nb < nbeta;nb++) - { - ifs >> lll[nb] >> jjj[nb]; - } - - ifs >> xmin >> rmax >> zmesh >> dx; - -// cerr<< " read_pseudo_addinfo Reading pseudo file \n"; -} -*/ void Pseudopot_upf::print_pseudo_upf(ofstream &ofs) { @@ -1146,64 +1002,6 @@ void Pseudopot_upf::print_pseudo_upf(ofstream &ofs) ofs << " iw=" << i << " els=" << els[i] << " lchi=" << lchi[i] << " oc=" << oc[i] << endl; } - // print PP_MESH - /* - ofs.setf(ios::scientific, ios::floatfield); - - for(int i=0; i<100; ++i) - { - ofs << setw(25) << r[i] << setw(25) << rab[i] << setw(25) << vloc[i] << setw(25) << rho_at[i] << endl; - } - - ofs << "chi:" << endl; - for(int ir=0; ir<100; ++ir) - { - for(int i=0; ihas_so && LSPINORB) {error++; cout<<"warning_quit! no soc upf used for lspinorb calculation, error!"<has_so || LSPINORB) return error; + if(!this->has_so && LSPINORB) + { + error++; + cout<<"warning_quit! no soc upf used for lspinorb calculation, error!"<has_so || LSPINORB) + { + return error; + } + int new_nbeta = 0; //calculate the new nbeta for(int nb=0; nb< this->nbeta; nb++) { new_nbeta++; if(this->lll[nb] != 0 && abs(this->jjj[nb] - this->lll[nb] - 0.5) < 1e-6) //two J = l +- 0.5 average to one + { new_nbeta--; + } } + this->nbeta = new_nbeta; matrix dion_new; dion_new.create(this->nbeta, this->nbeta); @@ -2111,8 +1835,6 @@ int Pseudopot_upf::average_p() int ind=0, ind1=0; if(l != 0) { -// cout<jjj[old_nbeta] <<" "<< this->lll[old_nbeta]<<" "<jjj[old_nbeta] - this->lll[old_nbeta] + 0.5)<jjj[old_nbeta+1] <<" "<< this->lll[old_nbeta+1]<<" "<jjj[old_nbeta+1] - this->lll[old_nbeta+1] + 0.5)<jjj[old_nbeta] - this->lll[old_nbeta] + 0.5) < 1e-6) { if(abs(this->jjj[old_nbeta+1]-this->lll[old_nbeta+1]-0.5)>1e-6) @@ -2141,17 +1863,24 @@ int Pseudopot_upf::average_p() this->dion(nb, nb) = vion1; old_nbeta++; } - else{ + else + { for(int ir = 0; irmesh;ir++) this->beta(nb, ir) = this->beta(old_nbeta, ir); this->dion(nb, nb) = this->dion(old_nbeta, old_nbeta); } this->lll[nb] = this->lll[old_nbeta]; //reset the lll index, ignore jjj index } + //store the old dion and then recreate dion for(int i=0;inbeta; i++) + { for(int j=0;jnbeta;j++) + { dion_new(i,j) = this->dion(i,j); + } + } + this->dion = dion_new; // this->dion.create(this->nbeta, this->nbeta); // for(int i=0;inbeta; i++) @@ -2163,8 +1892,11 @@ int Pseudopot_upf::average_p() { new_nwfc++; if(this->lchi[nb] != 0 && abs(this->jchi[nb] - this->lchi[nb] - 0.5)<1e-6) + { new_nwfc--; + } } + this->nwfc = new_nwfc; int old_nwfc=0; for(int nb=0; nbnwfc; nb++) diff --git a/ABACUS.develop/source/src_pw/pseudopot_upf.h b/ABACUS.develop/source/src_pw/pseudopot_upf.h index 0050cdb340..819e13c118 100644 --- a/ABACUS.develop/source/src_pw/pseudopot_upf.h +++ b/ABACUS.develop/source/src_pw/pseudopot_upf.h @@ -61,8 +61,6 @@ class Pseudopot_upf 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 -// double *rinner; // nl_8 // rinner(0:2*lmax) : for r < rinner(i) Q functions are pseudized - //(not read if nqf=0) // 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) diff --git a/ABACUS.develop/source/src_pw/stress.cpp b/ABACUS.develop/source/src_pw/stress.cpp index 4c5bc8c9d6..080066d194 100644 --- a/ABACUS.develop/source/src_pw/stress.cpp +++ b/ABACUS.develop/source/src_pw/stress.cpp @@ -11,7 +11,7 @@ #include "H_Hartree_pw.h" #include "H_XC_pw.h" -void Stress::cal_stress() +void Stress::cal_stress(void) { TITLE("Stress","cal_stress"); timer::tick("Stress","cal_stress",'E'); @@ -30,6 +30,7 @@ void Stress::cal_stress() sigmaxcc[i][j] = 0.0; } } + //kinetic contribution stres_knl(); @@ -231,23 +232,33 @@ void Stress::stres_knl(void) TITLE("Stress","stres_nkl"); double *kfac; - double **gk; - gk=new double* [3]; - double tbsp,gk2,arg; - int ik,l,m,i,j,ibnd,is; - int npw; + double **gk=new double*[3]; + double tbsp=0.0; + double gk2=0.0; + double arg=0.0; + int ik=0; + int l=0; + int m=0; + int i=0; + int j=0; + int ibnd=0; + int is=0; + int npw=0; tbsp=2.0/sqrt(PI); int npwx=0; - int qtot = 0; + int qtot=0; for(int ik=0; ik *Porter = UFFT.porter; ZEROS( Porter, pw.nrxx ); @@ -591,21 +595,6 @@ void Stress::stres_loc(void) } pw.FFT_chg.FFT3D(Porter, -1); - /* complex *psic = new complex [pw.nrxx]; - double *psic0 = new double[pw.nrxx]; - ZEROS( psic0, pw.nrxx); - for(int is=0; is(psic0[ir], 0.0); - } - } - - pw.FFT_chg.FFT3D(psic, -1) ; - */ - if(INPUT.gamma_only==1) fact=2.0; else fact=1.0;