forked from abacusmodeling/abacus-develop
-
Notifications
You must be signed in to change notification settings - Fork 134
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
7 changed files
with
318 additions
and
16 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
292 changes: 292 additions & 0 deletions
292
source/module_hamilt_general/module_xc/xc_functional_libxc_tools.cpp
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,292 @@ | ||
| #ifdef USE_LIBXC | ||
|
|
||
| #include "xc_functional_libxc.h" | ||
| #include "xc_functional.h" | ||
| #include "module_elecstate/module_charge/charge.h" | ||
|
|
||
| // converting rho (abacus=>libxc) | ||
| std::vector<double> XC_Functional_Libxc::convert_rho( | ||
| const int nspin, | ||
| const std::size_t nrxx, | ||
| const Charge* const chr) | ||
| { | ||
| std::vector<double> rho(nrxx*nspin); | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for collapse(2) schedule(static, 1024) | ||
| #endif | ||
| for( int is=0; is<nspin; ++is ) | ||
| for( int ir=0; ir<nrxx; ++ir ) | ||
| rho[ir*nspin+is] = chr->rho[is][ir] + 1.0/nspin*chr->rho_core[ir]; | ||
| return rho; | ||
| } | ||
|
|
||
| // converting rho (abacus=>libxc) | ||
| std::tuple<std::vector<double>, std::vector<double>> | ||
| XC_Functional_Libxc::convert_rho_amag_nspin4( | ||
| const int nspin, | ||
| const std::size_t nrxx, | ||
| const Charge* const chr) | ||
| { | ||
| assert(GlobalV::NSPIN==4); | ||
| std::vector<double> rho(nrxx*nspin); | ||
| std::vector<double> amag(nrxx); | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for | ||
| #endif | ||
| for( int ir=0; ir<nrxx; ++ir ) | ||
| { | ||
| const double arhox = std::abs( chr->rho[0][ir] + chr->rho_core[ir] ); | ||
| amag[ir] = std::sqrt( std::pow(chr->rho[1][ir],2) | ||
| + std::pow(chr->rho[2][ir],2) | ||
| + std::pow(chr->rho[3][ir],2) ); | ||
| const double amag_clip = (amag[ir]<arhox) ? amag[ir] : arhox; | ||
| rho[ir*nspin+0] = (arhox + amag_clip) / 2.0; | ||
| rho[ir*nspin+1] = (arhox - amag_clip) / 2.0; | ||
| } | ||
| return std::make_tuple(std::move(rho), std::move(amag)); | ||
| } | ||
|
|
||
| // calculating grho | ||
| std::vector<std::vector<ModuleBase::Vector3<double>>> | ||
| XC_Functional_Libxc::cal_gdr( | ||
| const int nspin, | ||
| const std::vector<double> &rho, | ||
| const double tpiba, | ||
| const Charge* const chr) | ||
| { | ||
| std::vector<std::vector<ModuleBase::Vector3<double>>> gdr(nspin); | ||
| const std::size_t nrxx = rho.size(); | ||
| for( int is=0; is!=nspin; ++is ) | ||
| { | ||
| std::vector<double> rhor(nrxx); | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for schedule(static, 1024) | ||
| #endif | ||
| for(std::size_t ir=0; ir<nrxx; ++ir) | ||
| rhor[ir] = rho[ir*nspin+is]; | ||
| //------------------------------------------ | ||
| // initialize the charge density array in reciprocal space | ||
| // bring electron charge density from real space to reciprocal space | ||
| //------------------------------------------ | ||
| std::vector<std::complex<double>> rhog(chr->rhopw->npw); | ||
| chr->rhopw->real2recip(rhor.data(), rhog.data()); | ||
|
|
||
| //------------------------------------------- | ||
| // compute the gradient of charge density and | ||
| // store the gradient in gdr[is] | ||
| //------------------------------------------- | ||
| gdr[is].resize(nrxx); | ||
| XC_Functional::grad_rho(rhog.data(), gdr[is].data(), chr->rhopw, tpiba); | ||
| } // end for(is) | ||
| return gdr; | ||
| } | ||
|
|
||
| // converting grho (abacus=>libxc) | ||
| std::vector<double> XC_Functional_Libxc::convert_sigma( | ||
| const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr) | ||
| { | ||
| const std::size_t nspin = gdr.size(); | ||
| assert(nspin>0); | ||
| const std::size_t nrxx = gdr[0].size(); | ||
| for(std::size_t is=1; is<nspin; ++is) | ||
| assert(nrxx==gdr[is].size()); | ||
|
|
||
| std::vector<double> sigma( nrxx * ((1==nspin)?1:3) ); | ||
| if( 1==nspin ) | ||
| { | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for schedule(static, 1024) | ||
| #endif | ||
| for( std::size_t ir=0; ir<nrxx; ++ir ) | ||
| sigma[ir] = gdr[0][ir]*gdr[0][ir]; | ||
| } | ||
| else | ||
| { | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for schedule(static, 256) | ||
| #endif | ||
| for( std::size_t ir=0; ir<nrxx; ++ir ) | ||
| { | ||
| sigma[ir*3] = gdr[0][ir]*gdr[0][ir]; | ||
| sigma[ir*3+1] = gdr[0][ir]*gdr[1][ir]; | ||
| sigma[ir*3+2] = gdr[1][ir]*gdr[1][ir]; | ||
| } | ||
| } | ||
| return sigma; | ||
| } | ||
|
|
||
| // sgn for threshold mask | ||
| std::vector<double> XC_Functional_Libxc::cal_sgn( | ||
| const double rho_threshold, | ||
| const double grho_threshold, | ||
| const xc_func_type &func, | ||
| const int nspin, | ||
| const std::size_t nrxx, | ||
| const std::vector<double> &rho, | ||
| const std::vector<double> &sigma) | ||
| { | ||
| std::vector<double> sgn(nrxx*nspin, 1.0); | ||
| // in the case of GGA correlation for polarized case, | ||
| // a cutoff for grho is required to ensure that libxc gives reasonable results | ||
| if(nspin==2 && func.info->family != XC_FAMILY_LDA && func.info->kind==XC_CORRELATION) | ||
| { | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for schedule(static, 512) | ||
| #endif | ||
| for( int ir=0; ir<nrxx; ++ir ) | ||
| { | ||
| if ( rho[ir*2]<rho_threshold || std::sqrt(std::abs(sigma[ir*3]))<grho_threshold ) | ||
| sgn[ir*2] = 0.0; | ||
| if ( rho[ir*2+1]<rho_threshold || std::sqrt(std::abs(sigma[ir*3+2]))<grho_threshold ) | ||
| sgn[ir*2+1] = 0.0; | ||
| } | ||
| } | ||
| return sgn; | ||
| } | ||
|
|
||
| // converting etxc from exc (libxc=>abacus) | ||
| double XC_Functional_Libxc::convert_etxc( | ||
| const int nspin, | ||
| const std::size_t nrxx, | ||
| const std::vector<double> &sgn, | ||
| const std::vector<double> &rho, | ||
| std::vector<double> exc) | ||
| { | ||
| double etxc = 0.0; | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256) | ||
| #endif | ||
| for( int is=0; is<nspin; ++is ) | ||
| for( int ir=0; ir<nrxx; ++ir ) | ||
| etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is]; | ||
| return etxc; | ||
| } | ||
|
|
||
| // converting etxc from exc (libxc=>abacus) | ||
| double XC_Functional_Libxc::convert_vtxc_v( | ||
| const xc_func_type &func, | ||
| const int nspin, | ||
| const std::size_t nrxx, | ||
| const std::vector<double> &sgn, | ||
| const std::vector<double> &rho, | ||
| const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr, | ||
| const std::vector<double> &vrho, | ||
| const std::vector<double> &vsigma, | ||
| const double tpiba, | ||
| const Charge* const chr, | ||
| ModuleBase::matrix &v) | ||
| { | ||
| assert(v.nr==nspin); | ||
| assert(v.nc==nrxx); | ||
|
|
||
| double vtxc = 0.0; | ||
|
|
||
| #ifdef _OPENMP | ||
| #pragma omp parallel for collapse(2) reduction(+:vtxc) schedule(static, 256) | ||
| #endif | ||
| for( int is=0; is<nspin; ++is ) | ||
| { | ||
| for( std::size_t ir=0; ir<nrxx; ++ir ) | ||
| { | ||
| const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is] * sgn[ir*nspin+is]; | ||
| v(is,ir) += v_tmp; | ||
| vtxc += v_tmp * rho[ir*nspin+is]; | ||
| } | ||
| } | ||
|
|
||
| if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA) | ||
| { | ||
| const std::vector<std::vector<double>> dh = XC_Functional_Libxc::cal_dh(nspin, nrxx, sgn, gdr, vsigma, tpiba, chr); | ||
|
|
||
| double rvtxc = 0.0; | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for collapse(2) reduction(+:rvtxc) schedule(static, 256) | ||
| #endif | ||
| for( int is=0; is<nspin; ++is ) | ||
| { | ||
| for( std::size_t ir=0; ir<nrxx; ++ir ) | ||
| { | ||
| rvtxc += dh[is][ir] * rho[ir*nspin+is]; | ||
| v(is,ir) -= dh[is][ir]; | ||
| } | ||
| } | ||
|
|
||
| vtxc -= rvtxc; | ||
| } // end if(func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)) | ||
|
|
||
| return vtxc; | ||
| } | ||
|
|
||
|
|
||
| // dh for gga v | ||
| std::vector<std::vector<double>> XC_Functional_Libxc::cal_dh( | ||
| const int nspin, | ||
| const std::size_t nrxx, | ||
| const std::vector<double> &sgn, | ||
| const std::vector<std::vector<ModuleBase::Vector3<double>>> &gdr, | ||
| const std::vector<double> &vsigma, | ||
| const double tpiba, | ||
| const Charge* const chr) | ||
| { | ||
| std::vector<std::vector<ModuleBase::Vector3<double>>> h( | ||
| nspin, | ||
| std::vector<ModuleBase::Vector3<double>>(nrxx) ); | ||
| if( 1==nspin ) | ||
| { | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for schedule(static, 1024) | ||
| #endif | ||
| for( std::size_t ir=0; ir<nrxx; ++ir ) | ||
| h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir]; | ||
| } | ||
| else | ||
| { | ||
| #ifdef _OPENMP | ||
| #pragma omp parallel for schedule(static, 1024) | ||
| #endif | ||
| for( std::size_t ir=0; ir< nrxx; ++ir ) | ||
| { | ||
| h[0][ir] = 2.0 * (gdr[0][ir] * vsigma[ir*3 ] * sgn[ir*2 ] * 2.0 | ||
| + gdr[1][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); | ||
| h[1][ir] = 2.0 * (gdr[1][ir] * vsigma[ir*3+2] * sgn[ir*2+1] * 2.0 | ||
| + gdr[0][ir] * vsigma[ir*3+1] * sgn[ir*2] * sgn[ir*2+1]); | ||
| } | ||
| } | ||
|
|
||
| // define two dimensional array dh [ nspin, nrxx ] | ||
| std::vector<std::vector<double>> dh(nspin, std::vector<double>(nrxx)); | ||
| for( int is=0; is!=nspin; ++is ) | ||
| XC_Functional::grad_dot( h[is].data(), dh[is].data(), chr->rhopw, tpiba); | ||
|
|
||
| return dh; | ||
| } | ||
|
|
||
|
|
||
| // convert v for NSPIN=4 | ||
| ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4( | ||
| const std::size_t nrxx, | ||
| const Charge* const chr, | ||
| const std::vector<double> &amag, | ||
| const ModuleBase::matrix &v) | ||
| { | ||
| assert(GlobalV::NSPIN==4); | ||
| constexpr double vanishing_charge = 1.0e-10; | ||
| ModuleBase::matrix v_nspin4(GlobalV::NSPIN, nrxx); | ||
| for( int ir=0; ir<nrxx; ++ir ) | ||
| v_nspin4(0,ir) = 0.5 * (v(0,ir)+v(1,ir)); | ||
| if(GlobalV::DOMAG || GlobalV::DOMAG_Z) | ||
| { | ||
| for( int ir=0; ir<nrxx; ++ir ) | ||
| { | ||
| if ( amag[ir] > vanishing_charge ) | ||
| { | ||
| const double vs = 0.5 * (v(0,ir)-v(1,ir)); | ||
| for(int ipol=1; ipol<GlobalV::NSPIN; ++ipol) | ||
| v_nspin4(ipol,ir) = vs * chr->rho[ipol][ir] / amag[ir]; | ||
| } | ||
| } | ||
| } | ||
| return v_nspin4; | ||
| } | ||
|
|
||
| #endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.