From f53310c6131db00e866ec650dd01041f1b557753 Mon Sep 17 00:00:00 2001 From: linpz Date: Thu, 6 May 2021 11:44:02 +0800 Subject: [PATCH 1/2] 1. refactor H_Hartree_pw::v_hartree() --- ABACUS.develop/source/src_pw/H_Hartree_pw.cpp | 78 ++++++------------- ABACUS.develop/source/src_pw/H_Hartree_pw.h | 8 +- ABACUS.develop/source/src_pw/potential.cpp | 2 +- 3 files changed, 27 insertions(+), 61 deletions(-) diff --git a/ABACUS.develop/source/src_pw/H_Hartree_pw.cpp b/ABACUS.develop/source/src_pw/H_Hartree_pw.cpp index 84175bccf4..d74e21d199 100644 --- a/ABACUS.develop/source/src_pw/H_Hartree_pw.cpp +++ b/ABACUS.develop/source/src_pw/H_Hartree_pw.cpp @@ -5,41 +5,29 @@ double H_Hartree_pw::hartree_energy=0.0; //-------------------------------------------------------------------- // Transform charge density to hartree potential. //-------------------------------------------------------------------- -void H_Hartree_pw::v_hartree( +matrix H_Hartree_pw::v_hartree( const UnitCell &cell, PW_Basis &pwb, - const Use_FFT &ufft, - const int &nspin, - matrix &v, - double** rho) + const int &nspin, + const double*const*const rho) { TITLE("H_Hartree_pw","v_hartree"); timer::tick("H_Hartree_pw","v_hartree"); - complex *Porter = ufft.porter; - // Hartree potential VH(r) from n(r) - ZEROS( Porter, pwb.nrxx ); - int nspin0 = 1; - if(nspin==2) nspin0 = nspin; + vector> Porter(pwb.nrxx); + const int nspin0 = (nspin==2) ? 2 : 1; for(int is=0; is( rho[is][ir], 0.0 ); - } - } - //============================= // bring rho (aux) to G space //============================= - pwb.FFT_chg.FFT3D(Porter, -1); + pwb.FFT_chg.FFT3D(Porter.data(), -1); //double charge; //if (pwb.gstart == 1) - //{ // charge = cell.omega * Porter[pwb.ig2fftc[0]].real(); - //} //OUT(ofs_running, "v_h charge", charge); //======================================================= @@ -48,9 +36,7 @@ void H_Hartree_pw::v_hartree( double ehart = 0.0; - complex *vh_g = new complex[pwb.ngmc]; - ZEROS(vh_g, pwb.ngmc); - + vector> vh_g(pwb.ngmc); for (int ig = pwb.gstart; ig Date: Thu, 6 May 2021 18:20:20 +0800 Subject: [PATCH 2/2] 1. refactor Potential_Libxc --- ABACUS.develop/source/src_pw/potential.cpp | 7 +- .../source/src_pw/potential_libxc.cpp | 165 +++++++----------- .../source/src_pw/potential_libxc.h | 9 +- 3 files changed, 75 insertions(+), 106 deletions(-) diff --git a/ABACUS.develop/source/src_pw/potential.cpp b/ABACUS.develop/source/src_pw/potential.cpp index 16edd8e0a2..770d110756 100644 --- a/ABACUS.develop/source/src_pw/potential.cpp +++ b/ABACUS.develop/source/src_pw/potential.cpp @@ -292,7 +292,12 @@ void Potential::v_of_rho //---------------------------------------------------------- #ifdef USE_LIBXC - Potential_Libxc::v_xc(rho_in, H_XC_pw::etxc, H_XC_pw::vtxc, v_in); + H_XC_pw::etxc = 0.0; + H_XC_pw::vtxc = 0.0; + const auto etxc_vtxc_v = Potential_Libxc::v_xc(rho_in, CHR.rho_core); + H_XC_pw::etxc += std::get<0>(etxc_vtxc_v); + H_XC_pw::vtxc += std::get<1>(etxc_vtxc_v); + v_in += std::get<2>(etxc_vtxc_v); #else H_XC_pw::v_xc(pw.nrxx, pw.ncxyz, ucell.omega, rho_in, CHR.rho_core, v_in); #endif diff --git a/ABACUS.develop/source/src_pw/potential_libxc.cpp b/ABACUS.develop/source/src_pw/potential_libxc.cpp index 371f3e45ee..dc0cd45b5e 100644 --- a/ABACUS.develop/source/src_pw/potential_libxc.cpp +++ b/ABACUS.develop/source/src_pw/potential_libxc.cpp @@ -13,18 +13,16 @@ #include "src_global/global_function.h" #include "src_pw/xc_gga_pw.h" -void Potential_Libxc::v_xc( - const double * const * const rho_in, // CHR.rho_core may be needed in future - double &etxc, - double &vtxc, - matrix &v) +std::tuple Potential_Libxc::v_xc( + const double * const * const rho_in, + const double * const rho_core_in) { TITLE("Potential_Libxc","v_xc"); timer::tick("Potential_Libxc","v_xc"); - etxc = 0.0; - vtxc = 0.0; - v.zero_out(); + double etxc = 0.0; + double vtxc = 0.0; + matrix v(NSPIN,pw.nrxx); //---------------------------------------------------------- // xc_func_type is defined in Libxc package @@ -36,34 +34,47 @@ void Potential_Libxc::v_xc( // the type of input_tmp is automatically set to 'tuple' // [rho, sigma, gdr] = cal_input( funcs, rho_in ); - const auto input_tmp = cal_input( funcs, rho_in ); - - // obtain the electron charge density rho from 'get' function (first variable) + const auto input_tmp = cal_input( funcs, rho_in, rho_core_in ); const vector &rho = std::get<0>(input_tmp); - - // obtain the sigma from 'get' function (second variable) const vector &sigma = std::get<1>(input_tmp); for( xc_func_type &func : funcs ) { + // jiyy add for threshold + constexpr double rho_threshold = 1E-6; + constexpr double grho_threshold = 1E-10; + xc_func_set_dens_threshold(&func, rho_threshold); + // sgn for threshold mask + const vector sgn = [&]() -> vector + { + vector sgn( pw.nrxx * nspin0(), 1.0); + if(nspin0()==2 && func.info->family != XC_FAMILY_LDA && func.info->kind==XC_CORRELATION) + { + for( size_t ir=0; ir!=pw.nrxx; ++ir ) + { + if ( rho[ir*2] exc ( pw.nrxx ); vector vrho ( pw.nrxx * nspin0() ); vector vsigma( pw.nrxx * ((1==nspin0())?1:3) ); // cal etxc from rho, exc - auto process_exc = [&](vector &sgn) + auto process_exc = [&]() { for( size_t is=0; is!=nspin0(); ++is ) - { for( size_t ir=0; ir!=pw.nrxx; ++ir ) - { etxc += e2 * exc[ir] * rho[ir*nspin0()+is] * sgn[ir*nspin0()+is]; - } - } }; // cal vtx, v from rho_in, vrho - auto process_vrho = [&](vector &sgn) + auto process_vrho = [&]() { if(nspin0()==1 || NSPIN==2) { @@ -107,7 +118,7 @@ void Potential_Libxc::v_xc( }; // cal vtxc, v from rho_in, rho, gdr, vsigma - auto process_vsigma = [&](vector &sgn) + auto process_vsigma = [&]() { const std::vector>> &gdr = std::get<2>(input_tmp); @@ -123,38 +134,27 @@ void Potential_Libxc::v_xc( { for( size_t ir=0; ir!=pw.nrxx; ++ir ) { - h[0][ir] = e2 * (gdr[0][ir] * vsigma[ir*3 ] * 2.0 - * sgn[ir*2 ] + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); - h[1][ir] = e2 * (gdr[1][ir] * vsigma[ir*3+2] * 2.0 - * sgn[ir*2+1] + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + h[0][ir] = e2 * (gdr[0][ir] * vsigma[ir*3 ] * 2.0 * sgn[ir*2 ] + + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); + h[1][ir] = e2 * (gdr[1][ir] * vsigma[ir*3+2] * 2.0 * sgn[ir*2+1] + + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); } } // define two dimensional array dh [ nspin, pw.nrxx ] vector> dh(nspin0(), vector(pw.nrxx)); for( size_t is=0; is!=nspin0(); ++is ) - { GGA_PW::grad_dot( VECTOR_TO_PTR(h[is]), VECTOR_TO_PTR(dh[is]) ); - } - for( size_t is=0; is!=nspin0(); ++is ) - { for( size_t ir=0; ir!=pw.nrxx; ++ir ) - { vtxc -= dh[is][ir] * rho[ir*nspin0()+is]; - } - } if(nspin0()==1 || NSPIN==2) { for( size_t is=0; is!=nspin0(); ++is ) - { for( size_t ir=0; ir!=pw.nrxx; ++ir ) - { v(is,ir) -= dh[is][ir]; - } - } } else // may need updates for SOC { @@ -163,40 +163,17 @@ void Potential_Libxc::v_xc( { v(0,ir) -= 0.5 * (dh[0][ir] + dh[1][ir]); const double amag = sqrt( pow(rho_in[1][ir],2) + pow(rho_in[2][ir],2) + pow(rho_in[3][ir],2) ); - const double neg = (soc.lsign && rho_in[1][ir]*soc.ux[0] - +rho_in[2][ir]*soc.ux[1]+rho_in[3][ir]*soc.ux[2]<=0) ? -1 : 1; + const double neg = (soc.lsign + && rho_in[1][ir]*soc.ux[0]+rho_in[2][ir]*soc.ux[1]+rho_in[3][ir]*soc.ux[2]<=0) + ? -1 : 1; if(amag > vanishing_charge) { for(int i=1;i<4;i++) - { v(i,ir) -= neg * 0.5 * (dh[0][ir]-dh[1][ir]) * rho_in[i][ir] / amag; - } } } } }; - - // jiyy add for threshold - constexpr double rho_threshold = 1E-6; - constexpr double grho_threshold = 1E-10; - xc_func_set_dens_threshold(&func, rho_threshold); - - // LPZ: Explain what is sgn - vector sgn( pw.nrxx * nspin0(), 1.0); - if(nspin0()==2 && func.info->family != XC_FAMILY_LDA && func.info->kind==XC_CORRELATION) - { - for( size_t ir=0; ir!=pw.nrxx; ++ir ) - { - if ( rho[ir*2]family" includes LDA, GGA, Hybrid functional (HYB). @@ -205,20 +182,19 @@ void Potential_Libxc::v_xc( { case XC_FAMILY_LDA: // call Libxc function: xc_lda_exc_vxc - xc_lda_exc_vxc( &func, pw.nrxx, - VECTOR_TO_PTR(rho), VECTOR_TO_PTR(exc), VECTOR_TO_PTR(vrho) ); - process_exc(sgn); - process_vrho(sgn); + xc_lda_exc_vxc( &func, pw.nrxx, rho.data(), + exc.data(), vrho.data() ); + process_exc(); + process_vrho(); break; case XC_FAMILY_GGA: case XC_FAMILY_HYB_GGA: // call Libxc function: xc_gga_exc_vxc - xc_gga_exc_vxc( &func, pw.nrxx, - VECTOR_TO_PTR(rho), VECTOR_TO_PTR(sigma), VECTOR_TO_PTR(exc), - VECTOR_TO_PTR(vrho), VECTOR_TO_PTR(vsigma) ); - process_exc(sgn); - process_vrho(sgn); - process_vsigma(sgn); + xc_gga_exc_vxc( &func, pw.nrxx, rho.data(), sigma.data(), + exc.data(), vrho.data(), vsigma.data() ); + process_exc(); + process_vrho(); + process_vsigma(); break; default: throw domain_error("func.info->family ="+TO_STRING(func.info->family) @@ -229,7 +205,6 @@ void Potential_Libxc::v_xc( xc_func_end(&func); } - //------------------------------------------------- // for MPI, reduce the exchange-correlation energy //------------------------------------------------- @@ -240,7 +215,7 @@ void Potential_Libxc::v_xc( vtxc *= ucell.omega / pw.ncxyz; timer::tick("Potential_Libxc","v_xc"); - return; + return std::make_tuple( etxc, vtxc, v ); } @@ -306,8 +281,8 @@ std::vector Potential_Libxc::init_func() } else { - throw domain_error("iexch="+TO_STRING(xcf.iexch_now)+", igcx=" - +TO_STRING(xcf.igcx_now)+" unfinished in "+TO_STRING(__FILE__)+" line "+TO_STRING(__LINE__)); + throw domain_error("iexch="+TO_STRING(xcf.iexch_now)+", igcx="+TO_STRING(xcf.igcx_now) + +" unfinished in "+TO_STRING(__FILE__)+" line "+TO_STRING(__LINE__)); } //-------------------------------------- @@ -323,8 +298,8 @@ std::vector Potential_Libxc::init_func() } else { - throw domain_error("icorr="+TO_STRING(xcf.icorr_now)+", igcc=" - +TO_STRING(xcf.igcc_now)+" unfinished in "+TO_STRING(__FILE__)+" line "+TO_STRING(__LINE__)); + throw domain_error("icorr="+TO_STRING(xcf.icorr_now)+", igcc="+TO_STRING(xcf.igcc_now) + +" unfinished in "+TO_STRING(__FILE__)+" line "+TO_STRING(__LINE__)); } return funcs; @@ -338,13 +313,14 @@ std::tuple< std::vector, std::vector>> > Potential_Libxc::cal_input( const std::vector &funcs, - const double * const * const rho_in ) + const double * const * const rho_in, + const double * const rho_core_in ) { // ..., ↑_{i},↓_{i}, ↑_{i+1},↓_{i+1}, ... std::vector rho; bool finished_rho = false; - // here we assume CHR.rho_core equally exists in different spins, may need double check + // here we assume rho_core equally exists in different spins, may need double check auto cal_rho = [&]() { if(!finished_rho) @@ -353,12 +329,8 @@ Potential_Libxc::cal_input( if(nspin0()==1 || NSPIN==2) { for( size_t is=0; is!=nspin0(); ++is ) - { for( size_t ir=0; ir!=pw.nrxx; ++ir ) - { - rho[ir*nspin0()+is] = rho_in[is][ir] + 1.0/nspin0()*CHR.rho_core[ir]; - } - } + rho[ir*nspin0()+is] = rho_in[is][ir] + 1.0/nspin0()*rho_core_in[ir]; } else // may need updates for SOC { @@ -370,10 +342,11 @@ Potential_Libxc::cal_input( { const double amag = sqrt( pow(rho_in[1][ir],2) + pow(rho_in[2][ir],2) + pow(rho_in[3][ir],2) ); const double neg = (soc.lsign && rho_in[1][ir]*soc.ux[0] - +rho_in[2][ir]*soc.ux[1] - +rho_in[3][ir]*soc.ux[2]<=0) ? -1 : 1; - rho[ir*2] = 0.5 * (rho_in[0][ir] + neg * amag) + 0.5 * CHR.rho_core[ir]; - rho[ir*2+1] = 0.5 * (rho_in[0][ir] - neg * amag) + 0.5 * CHR.rho_core[ir]; + +rho_in[2][ir]*soc.ux[1] + +rho_in[3][ir]*soc.ux[2]<=0) + ? -1 : 1; + rho[ir*2] = 0.5 * (rho_in[0][ir] + neg * amag) + 0.5 * rho_core_in[ir]; + rho[ir*2+1] = 0.5 * (rho_in[0][ir] - neg * amag) + 0.5 * rho_core_in[ir]; } } } @@ -392,20 +365,13 @@ Potential_Libxc::cal_input( { vector rhor(pw.nrxx); for(int ir=0; ir> rhog(pw.ngmc); - - //------------------------------------------- - // bring electron charge density from real - // space to reciprocal space - //------------------------------------------- CHR.set_rhog(rhor.data(), rhog.data()); //------------------------------------------- @@ -462,9 +428,8 @@ Potential_Libxc::cal_input( cal_sigma(); break; default: - throw domain_error("func.info->family =" - +TO_STRING(func.info->family)+" unfinished in " - +TO_STRING(__FILE__)+" line "+TO_STRING(__LINE__)); + throw domain_error("func.info->family ="+TO_STRING(func.info->family) + +" unfinished in "+TO_STRING(__FILE__)+" line "+TO_STRING(__LINE__)); break; } } diff --git a/ABACUS.develop/source/src_pw/potential_libxc.h b/ABACUS.develop/source/src_pw/potential_libxc.h index 54c3352523..0328546af9 100644 --- a/ABACUS.develop/source/src_pw/potential_libxc.h +++ b/ABACUS.develop/source/src_pw/potential_libxc.h @@ -25,11 +25,9 @@ class Potential_Libxc // evaluate the exchange-correlation (XC) energy // by using the input charge density rho_in //------------------------------------------------ - static void v_xc( + static std::tuple v_xc( const double * const * const rho_in, - double &etxc, - double &vtxc, - matrix &v); + const double * const rho_core_in); private: @@ -50,7 +48,8 @@ class Potential_Libxc std::vector>> > cal_input( const std::vector &funcs, - const double * const * const rho_in ); + const double * const * const rho_in, + const double * const rho_core_in ); //---------------------------- // decide the value of spin