From 1d7e79f82bd5b4d997e287f7cf6b7634c99faf75 Mon Sep 17 00:00:00 2001 From: Wenfei Li <38569667+wenfei-li@users.noreply.github.com> Date: Thu, 19 Oct 2023 20:44:56 +0800 Subject: [PATCH] Feature : update energy calculation of PAW (#3077) * Feature : parallelization of PAW * update total energy for paw * update test result * fix bug * update libpaw --------- Co-authored-by: wenfei-li --- deps/libpaw_interface | 2 +- source/module_cell/module_paw/paw_cell.h | 6 ++ .../module_paw/paw_cell_libpaw.cpp | 26 ++++++-- source/module_cell/module_paw/paw_sphbes.cpp | 3 +- source/module_elecstate/elecstate_energy.cpp | 66 ++++++++++++++----- source/module_elecstate/fp_energy.cpp | 43 +++++++++++- tests/integrate/101_PW_15_paw/result.ref | 7 +- 7 files changed, 120 insertions(+), 33 deletions(-) diff --git a/deps/libpaw_interface b/deps/libpaw_interface index 1fe93a0a90..bdf13fb3a8 160000 --- a/deps/libpaw_interface +++ b/deps/libpaw_interface @@ -1 +1 @@ -Subproject commit 1fe93a0a9048bad5585d74715f31867a527775ee +Subproject commit bdf13fb3a837e71c5a9ecccd2caf71aeaa3578ad diff --git a/source/module_cell/module_paw/paw_cell.h b/source/module_cell/module_paw/paw_cell.h index af20a81d9a..34f60381c2 100644 --- a/source/module_cell/module_paw/paw_cell.h +++ b/source/module_cell/module_paw/paw_cell.h @@ -219,12 +219,16 @@ class Paw_Cell int get_libpaw_ixc() {return ixc;} int get_libpaw_xclevel() {return xclevel;} int get_nspin() {return nspden;} + + double get_epawdc() {return epawdc * 2.0;} int get_nfft() {return nfft;} int get_nrxx() {return nx*ny*num_z[GlobalV::RANK_IN_POOL];} int get_val(const int it) {return paw_element_list[it].get_zval();} int get_zat(const int it) {return paw_element_list[it].get_zat();} + double calculate_ecore(); + private: // Info to be passed to libpaw_interface: // 1. ecut, ecutpaw : kinetic energy cutoff of the planewave basis set @@ -242,6 +246,7 @@ class Paw_Cell double ecut, ecutpaw; std::vector rprimd, gprimd, gmet; + std::vector epsatm; double ucvol; std::vector ngfft, ngfftdg; int nfft; @@ -251,6 +256,7 @@ class Paw_Cell char* filename_list; int xclevel, ixc; int nspden, nsppol; + double epawdc; // Part IV. Calling Fortran subroutines from libpaw_interface public: diff --git a/source/module_cell/module_paw/paw_cell_libpaw.cpp b/source/module_cell/module_paw/paw_cell_libpaw.cpp index 21a693d3da..aae1bce182 100644 --- a/source/module_cell/module_paw/paw_cell_libpaw.cpp +++ b/source/module_cell/module_paw/paw_cell_libpaw.cpp @@ -180,6 +180,7 @@ void Paw_Cell::set_libpaw_atom(const int natom_in, const int ntypat_in, const in typat.resize(natom); xred.resize(3 * natom); + epsatm.resize(ntypat); for(int iat = 0; iat < natom; iat ++) { @@ -250,7 +251,7 @@ extern "C" { void prepare_libpaw_(double&,double&,double*,double*,double*,double&,int*,int*, // ecut ecutpaw gmet rprimd gprimd ucvol ngfft ngfftdg - int&,int&,int*,double*,int&,int&,char*,int&,int&); + int&,int&,int*,double*,int&,int&,char*,int&,int&,double*); // natom ntypat typat xred ixc xclevel filename_list nspden nsppol void get_vloc_ncoret_(int*, int&,int&, int&, double*,double*,double*,double&,double*,double*,double*); @@ -262,8 +263,8 @@ extern "C" void get_nhat_(int&, int&, double*, int*, int&, int&, double*, double*, double&, double*, double*); // natom,ntypat,xred, ngfft,nfft,nspden,gprimd, rprimd, ucvol, nhat, nhatgr - void calculate_dij_(int&, int&, int&, int&, int&, int&, double*, double&, double*, double*, double*); - // natom,ntypat,ixc, xclevel,nfft, nspden,xred, ucvol, gprimd, vks, vxc + void calculate_dij_(int&, int&, int&, int&, int&, int&, double*, double&, double*, double*, double*, double&); + // natom,ntypat,ixc, xclevel,nfft, nspden,xred, ucvol, gprimd, vks, vxc, epawdc void get_dij_(int&, int&, int&, double*); // iatom,size_dij,nspden,dij @@ -278,7 +279,7 @@ void Paw_Cell::prepare_paw() { prepare_libpaw_(ecut, ecutpaw, gmet.data(), rprimd.data(), gprimd.data(), ucvol, ngfft.data(), ngfftdg.data(), natom, ntypat, typat.data(), xred.data(), - ixc, xclevel, filename_list, nspden, nsppol); + ixc, xclevel, filename_list, nspden, nsppol, epsatm.data()); } void Paw_Cell::get_vloc_ncoret(double* vloc, double* ncoret) @@ -459,7 +460,7 @@ void Paw_Cell::calculate_dij(double* vks, double* vxc) if(GlobalV::RANK_IN_POOL == 0) { - calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree); + calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree,epawdc); } if(GlobalV::RANK_IN_POOL == 0) @@ -489,7 +490,7 @@ void Paw_Cell::calculate_dij(double* vks, double* vxc) } } } - calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree); + calculate_dij_(natom,ntypat,ixc,xclevel,nfft,nspden,xred.data(),ucvol,gprimd.data(),vks_hartree,vxc_hartree,epawdc); delete[] vks_hartree; delete[] vxc_hartree; #endif @@ -645,4 +646,17 @@ void Paw_Cell::set_sij() delete[] sij_libpaw; delete[] sij; } +} + +double Paw_Cell::calculate_ecore() +{ + double charge = 0.0; + double esum = 0.0; + for(int ia = 0; ia < natom; ia ++) + { + const int it = atom_type[ia]; + esum += epsatm[it]; + charge += paw_element_list[it].get_zval(); + } + return esum * charge / ucvol * 2.0; } \ No newline at end of file diff --git a/source/module_cell/module_paw/paw_sphbes.cpp b/source/module_cell/module_paw/paw_sphbes.cpp index c928ac451c..71360ba34e 100644 --- a/source/module_cell/module_paw/paw_sphbes.cpp +++ b/source/module_cell/module_paw/paw_sphbes.cpp @@ -218,7 +218,7 @@ void Paw_Element::spherical_bessel_function(const int l, const double xx, // search for case("r=a*(exp(d*i)-1)") in m_pawpsp.F90, for example double Paw_Element::simpson_integration(std::vector & f) const { - int ir=nr; + int ir=nr-1; while (std::abs(f[ir]) < 1e-20) { ir=ir-1; @@ -275,6 +275,7 @@ void Paw_Element::prepare_simpson_integration(const double r_for_intg, int & mes meshsz = ir; } + if(meshsz == nr) meshsz -= 1; //std::cout << "meshsz : " << meshsz << std::endl; //get simp_fact diff --git a/source/module_elecstate/elecstate_energy.cpp b/source/module_elecstate/elecstate_energy.cpp index 32c3b30f6c..8ff695bd2c 100644 --- a/source/module_elecstate/elecstate_energy.cpp +++ b/source/module_elecstate/elecstate_energy.cpp @@ -4,6 +4,10 @@ #include "elecstate_getters.h" #include "module_base/global_variable.h" #include "module_base/parallel_reduce.h" +#ifdef USE_PAW +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#endif namespace elecstate { @@ -98,41 +102,67 @@ double ElecState::cal_delta_eband() const const double* v_eff = this->pot->get_effective_v(0); const double* v_fixed = this->pot->get_fixed_v(); const double* v_ofk = nullptr; - if (get_xc_func_type() == 3 || get_xc_func_type() == 5) + +#ifdef USE_PAW + if(GlobalV::use_paw) { - v_ofk = this->pot->get_effective_vofk(0); + ModuleBase::matrix v_xc; + const std::tuple etxc_vtxc_v + = XC_Functional::v_xc(this->charge->nrxx, this->charge, &GlobalC::ucell); + v_xc = std::get<2>(etxc_vtxc_v); + + for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) + { + deband_aux -= this->charge->rho[0][ir] * v_xc(0,ir); + } + if (GlobalV::NSPIN == 2) + { + for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) + { + deband_aux -= this->charge->rho[1][ir] * v_xc(1,ir); + } + } } +#endif - for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) + if(!GlobalV::use_paw) { - deband_aux -= this->charge->rho[0][ir] * (v_eff[ir] - v_fixed[ir]); if (get_xc_func_type() == 3 || get_xc_func_type() == 5) { - deband_aux -= this->charge->kin_r[0][ir] * v_ofk[ir]; + v_ofk = this->pot->get_effective_vofk(0); } - } - if (GlobalV::NSPIN == 2) - { - v_eff = this->pot->get_effective_v(1); - v_ofk = this->pot->get_effective_vofk(1); for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { - deband_aux -= this->charge->rho[1][ir] * (v_eff[ir] - v_fixed[ir]); + deband_aux -= this->charge->rho[0][ir] * (v_eff[ir] - v_fixed[ir]); if (get_xc_func_type() == 3 || get_xc_func_type() == 5) { - deband_aux -= this->charge->kin_r[1][ir] * v_ofk[ir]; + deband_aux -= this->charge->kin_r[0][ir] * v_ofk[ir]; } } - } - else if (GlobalV::NSPIN == 4) - { - for (int is = 1; is < 4; is++) + + if (GlobalV::NSPIN == 2) { - v_eff = this->pot->get_effective_v(is); + v_eff = this->pot->get_effective_v(1); + v_ofk = this->pot->get_effective_vofk(1); for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) { - deband_aux -= this->charge->rho[is][ir] * v_eff[ir]; + deband_aux -= this->charge->rho[1][ir] * (v_eff[ir] - v_fixed[ir]); + if (get_xc_func_type() == 3 || get_xc_func_type() == 5) + { + deband_aux -= this->charge->kin_r[1][ir] * v_ofk[ir]; + } + } + } + else if (GlobalV::NSPIN == 4) + { + for (int is = 1; is < 4; is++) + { + v_eff = this->pot->get_effective_v(is); + for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++) + { + deband_aux -= this->charge->rho[is][ir] * v_eff[ir]; + } } } } diff --git a/source/module_elecstate/fp_energy.cpp b/source/module_elecstate/fp_energy.cpp index a6bbf82841..d6934ed99a 100644 --- a/source/module_elecstate/fp_energy.cpp +++ b/source/module_elecstate/fp_energy.cpp @@ -1,4 +1,8 @@ #include "fp_energy.h" +#include "module_base/global_variable.h" +#ifdef USE_PAW +#include "module_cell/module_paw/paw_cell.h" +#endif #include #include @@ -10,16 +14,49 @@ namespace elecstate /// @brief calculate etot double fenergy::calculate_etot() { - etot = eband + deband + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield + gatefield - + evdw + esol_el + esol_cav + edftu + edeepks_scf; + if(GlobalV::use_paw) + { + etot = eband + deband + etxc + ewald_energy - hartree_energy + demet + descf + exx + efield + gatefield + + evdw + esol_el + esol_cav + edftu + edeepks_scf; + } + else + { + etot = eband + deband + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield + gatefield + + evdw + esol_el + esol_cav + edftu + edeepks_scf; + } + +#ifdef USE_PAW + if(GlobalV::use_paw) + { + double ecore = GlobalC::paw_cell.calculate_ecore(); + double epawdc = GlobalC::paw_cell.get_epawdc(); + etot += ( ecore + epawdc ); + } +#endif return etot; } /// @brief calculate etot_harris double fenergy::calculate_harris() { - etot_harris = eband + deband_harris + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield + if(GlobalV::use_paw) + { + etot_harris = eband + deband_harris + etxc + ewald_energy - hartree_energy + demet + descf + exx + efield + + gatefield + evdw + esol_el + esol_cav + edftu + edeepks_scf; + } + else + { + etot_harris = eband + deband_harris + (etxc - etxcc) + ewald_energy + hartree_energy + demet + descf + exx + efield + gatefield + evdw + esol_el + esol_cav + edftu + edeepks_scf; + } +#ifdef USE_PAW + if(GlobalV::use_paw) + { + double ecore = GlobalC::paw_cell.calculate_ecore(); + double epawdc = GlobalC::paw_cell.get_epawdc(); + etot_harris += ( ecore + epawdc ); + } +#endif return etot_harris; } diff --git a/tests/integrate/101_PW_15_paw/result.ref b/tests/integrate/101_PW_15_paw/result.ref index 5a12da9f60..a3228cabc9 100644 --- a/tests/integrate/101_PW_15_paw/result.ref +++ b/tests/integrate/101_PW_15_paw/result.ref @@ -1,7 +1,6 @@ -etotref -261.3678632375260 -etotperatomref -130.6839316188 +etotref -200.9340475810059 +etotperatomref -100.4670237905 pointgroupref T_d spacegroupref O_h nksibzref 1 -totaltimeref 29.2697 -1.07 +totaltimeref 1.17