diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index b332f1fc80..f5b72d8125 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -129,6 +129,7 @@ - [out\_mul](#out_mul) - [out\_freq\_elec](#out_freq_elec) - [out\_chg](#out_chg) + - [out\_chg](#out_xc_r) - [out\_pot](#out_pot) - [out\_dm](#out_dm) - [out\_dm1](#out_dm1) @@ -1550,6 +1551,25 @@ These variables are used to control the output of properties. In molecular dynamics calculations, the output frequency is controlled by [out_interval](#out_interval). - **Default**: 0 3 +### out_xc_r + +- **Type**: Integer \[Integer\](optional) +- **Description**: + The first integer controls whether to output the exchange-correlation (in Bohr^-3) on real space grids using Libxc to folder `OUT.${suffix}/xc_r`: + - 0: rho, amag, sigma, exc + - 1: vrho, vsigma + - 2: v2rho2, v2rhosigma, v2sigma2 + - 3: v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3 + - 4: v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4 + The meaning of the files is presented in [Libxc](https://libxc.gitlab.io/manual/libxc-5.1.x/) + + The second integer controls the precision of the charge density output, if not given, will use `3` as default. + + --- + The circle order of the charge density on real space grids is: x is the outer loop, then y and finally z (z is moving fastest). + +- **Default**: -1 3 + ### out_pot - **Type**: Integer diff --git a/source/Makefile.Objects b/source/Makefile.Objects index d99287243d..950e978d41 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -507,6 +507,7 @@ OBJS_IO=input_conv.o\ write_dipole.o\ td_current_io.o\ write_wfc_r.o\ + write_libxc_r.o\ output_log.o\ output_mat_sparse.o\ para_json.o\ diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index 7bb5c83f26..0c1d76d533 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -11,6 +11,10 @@ #include "module_io/cif_io.h" #include "module_elecstate/module_charge/symmetry_rho.h" +#ifdef USE_LIBXC +#include "module_io/write_libxc_r.h" +#endif + namespace ModuleESolver { @@ -280,6 +284,23 @@ void ESolver_FP::after_scf(const int istep) &(GlobalC::ucell), PARAM.inp.out_elf[1]); } + + // 6) write xc(r) + if (PARAM.inp.out_xc_r[0]>=0) + { +#ifdef USE_LIBXC + ModuleIO::write_libxc_r( + PARAM.inp.out_xc_r[0], + XC_Functional::get_func_id(), + this->pw_rhod->nrxx, // number of real-space grid + GlobalC::ucell.omega, // volume of cell + GlobalC::ucell.tpiba, + &this->chr, + this); +#else + throw std::invalid_argument("out_xc_r must compile with libxc.\nSee "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); +#endif + } } } diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index 4014fb0c95..b641998f6a 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -27,6 +27,7 @@ list(APPEND objects write_dipole.cpp td_current_io.cpp write_wfc_r.cpp + write_libxc_r.cpp output_log.cpp para_json.cpp parse_args.cpp diff --git a/source/module_io/cube_io.h b/source/module_io/cube_io.h index a6a296b5dd..e316e0e0bb 100644 --- a/source/module_io/cube_io.h +++ b/source/module_io/cube_io.h @@ -77,22 +77,30 @@ void read_cube_core_mismatch( const int ny_read, const int nz_read); +#ifdef __MPI // when MPI: -// write data[ixy*nplane+iz] to file as order (ixy,iz) -// when serial: -// write data[iz*nxy+ixy] to file as order (ixy,iz) +// write data[ixy*nplane+iz*nld] to file as order (ixy,iz) void write_cube_core( std::ofstream &ofs_cube, -#ifdef __MPI const int bz, const int nbz, const int nplane, const int startz_current, -#endif const double*const data, const int nxy, const int nz, + const int nld, const int n_data_newline); +#else +// when serial: +// write data[iz*nxy+ixy] to file as order (ixy,iz) +void write_cube_core( + std::ofstream &ofs_cube, + const double*const data, + const int nxy, + const int nz, + const int n_data_newline); +#endif /** * @brief The trilinear interpolation method diff --git a/source/module_io/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp index 7709f60213..9678f24008 100644 --- a/source/module_io/read_input_item_output.cpp +++ b/source/module_io/read_input_item_output.cpp @@ -80,6 +80,20 @@ void ReadInput::item_output() read_sync_bool(input.out_wfc_r); this->add_item(item); } + { + Input_Item item("out_xc_r"); + item.annotation = "if >=0, output the derivatives of exchange correlation in realspace, second parameter controls the precision"; + item.read_value = [](const Input_Item& item, Parameter& para) { + size_t count = item.get_size(); + std::vector out_xc_r(count); // create a placeholder vector + std::transform(item.str_values.begin(), item.str_values.end(), out_xc_r.begin(), [](std::string s) { return std::stoi(s); }); + // assign non-negative values to para.input.out_xc_r + std::copy(out_xc_r.begin(), out_xc_r.end(), para.input.out_xc_r.begin()); + std::cout<add_item(item); + } { Input_Item item("printe"); item.annotation = "Print out energy for each band for every printe steps"; diff --git a/source/module_io/write_cube.cpp b/source/module_io/write_cube.cpp index 2573740218..56579697bd 100644 --- a/source/module_io/write_cube.cpp +++ b/source/module_io/write_cube.cpp @@ -127,10 +127,12 @@ void ModuleIO::write_cube( ofs_cube << std::scientific; } + constexpr int n_data_newline = 6; #ifdef __MPI - ModuleIO::write_cube_core(ofs_cube, bz, nbz, nplane, startz_current, data, nx*ny, nz, 6); + constexpr int nld = 1; + ModuleIO::write_cube_core(ofs_cube, bz, nbz, nplane, startz_current, data, nx*ny, nz, nld, n_data_newline); #else - ModuleIO::write_cube_core(ofs_cube, data, nx*ny, nz, 6); + ModuleIO::write_cube_core(ofs_cube, data, nx*ny, nz, n_data_newline); #endif if (my_rank == 0) @@ -141,28 +143,25 @@ void ModuleIO::write_cube( /// for cube file ofs_cube.close(); } - - return; } +#ifdef __MPI + void ModuleIO::write_cube_core( std::ofstream &ofs_cube, -#ifdef __MPI const int bz, const int nbz, const int nplane, const int startz_current, -#endif const double*const data, const int nxy, const int nz, + const int nld, const int n_data_newline) { ModuleBase::TITLE("ModuleIO", "write_cube_core"); -#ifdef __MPI - const int my_rank = GlobalV::MY_RANK; const int my_pool = GlobalV::MY_POOL; const int rank_in_pool = GlobalV::RANK_IN_POOL; @@ -177,7 +176,6 @@ void ModuleIO::write_cube_core( // num_z: how many planes on processor 'ip' std::vector num_z(nproc_in_pool, 0); - for (int iz = 0; iz < nbz; iz++) { const int ip = iz % nproc_in_pool; @@ -231,7 +229,7 @@ void ModuleIO::write_cube_core( // mohan change to rho_save on 2012-02-10 // because this can make our next restart calculation lead // to the same scf_thr as the one saved. - zpiece[ixy] = data[ixy * nplane + iz - startz_current]; + zpiece[ixy] = data[ixy * nplane + (iz - startz_current) * nld]; } } // case 2: > first part rho: send the rho to @@ -240,7 +238,7 @@ void ModuleIO::write_cube_core( { for (int ixy = 0; ixy < nxy; ixy++) { - zpiece[ixy] = data[ixy * nplane + iz - startz_current]; + zpiece[ixy] = data[ixy * nplane + (iz - startz_current) * nld]; } MPI_Send(zpiece.data(), nxy, MPI_DOUBLE, 0, tag, POOL_WORLD); } @@ -282,7 +280,18 @@ void ModuleIO::write_cube_core( /// for cube file } MPI_Barrier(MPI_COMM_WORLD); -#else +} + +#else // #ifdef __MPI + +void ModuleIO::write_cube_core( + std::ofstream &ofs_cube, + const double*const data, + const int nxy, + const int nz, + const int n_data_newline) +{ + ModuleBase::TITLE("ModuleIO", "write_cube_core"); for (int ixy = 0; ixy < nxy; ixy++) { for (int iz = 0; iz < nz; iz++) @@ -296,5 +305,6 @@ void ModuleIO::write_cube_core( } ofs_cube << "\n"; } -#endif } + +#endif // #ifdef __MPI diff --git a/source/module_io/write_libxc_r.cpp b/source/module_io/write_libxc_r.cpp new file mode 100644 index 0000000000..f260b82239 --- /dev/null +++ b/source/module_io/write_libxc_r.cpp @@ -0,0 +1,285 @@ +//====================== +// AUTHOR : Peize Lin +// DATE : 2024-09-12 +//====================== + +#ifdef USE_LIBXC + +#include "write_libxc_r.h" +#include "module_hamilt_general/module_xc/xc_functional.h" +#include "module_hamilt_general/module_xc/xc_functional_libxc.h" +#include "module_elecstate/module_charge/charge.h" +#include "module_esolver/esolver_fp.h" +#include "module_io/cube_io.h" +#include "module_base/global_variable.h" +#include "module_parameter/parameter.h" +#include "module_base/timer.h" +#include "module_base/tool_title.h" + +#include +#include +#include +#include +#include + +void ModuleIO::write_libxc_r( + const int order, + const std::vector &func_id, + const int &nrxx, // number of real-space grid + const double &omega, // volume of cell + const double tpiba, + const Charge*const chr, + const ModuleESolver::ESolver_FP*const esolver) +{ + ModuleBase::TITLE("ModuleIO","write_libxc_r"); + ModuleBase::timer::tick("ModuleIO","write_libxc_r"); + + const int nspin = + (PARAM.inp.nspin == 1 || ( PARAM.inp.nspin ==4 && !PARAM.globalv.domag && !PARAM.globalv.domag_z)) + ? 1 : 2; + + //---------------------------------------------------------- + // xc_func_type is defined in Libxc package + // to understand the usage of xc_func_type, + // use can check on website, for example: + // https://www.tddft.org/programs/libxc/manual/libxc-5.1.x/ + //---------------------------------------------------------- + + std::vector funcs = XC_Functional_Libxc::init_func( func_id, (1==nspin) ? XC_UNPOLARIZED:XC_POLARIZED ); + + const bool is_gga = [&funcs]() + { + for( xc_func_type &func : funcs ) + { + switch( func.info->family ) + { + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + return true; + } + } + return false; + }(); + + // converting rho + std::vector rho; + std::vector amag; + if(1==nspin || 2==PARAM.inp.nspin) + { + rho = XC_Functional_Libxc::convert_rho(nspin, nrxx, chr); + } + else + { + std::tuple,std::vector> rho_amag = XC_Functional_Libxc::convert_rho_amag_nspin4(nspin, nrxx, chr); + rho = std::get<0>(std::move(rho_amag)); + amag = std::get<1>(std::move(rho_amag)); + } + + std::vector sigma; + if(is_gga) + { + const std::vector>> gdr = XC_Functional_Libxc::cal_gdr(nspin, nrxx, rho, tpiba, chr); + sigma = XC_Functional_Libxc::convert_sigma(gdr); + } + + std::vector exc; + std::vector vrho; + std::vector vsigma; + std::vector v2rho2; + std::vector v2rhosigma; + std::vector v2sigma2; + std::vector v3rho3; + std::vector v3rho2sigma; + std::vector v3rhosigma2; + std::vector v3sigma3; + std::vector v4rho4; + std::vector v4rho3sigma; + std::vector v4rho2sigma2; + std::vector v4rhosigma3; + std::vector v4sigma4; + // attention: order 4321 don't break + switch( order ) + { + case 4: v4rho4.resize( nrxx * ((1==nspin)?1:5) ); + case 3: v3rho3.resize( nrxx * ((1==nspin)?1:4) ); + case 2: v2rho2.resize( nrxx * ((1==nspin)?1:3) ); + case 1: vrho .resize( nrxx * nspin ); + case 0: exc .resize( nrxx ); + break; + default: throw std::domain_error("order ="+std::to_string(order) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + if(is_gga) + { + switch( order ) + { + case 4: v4rho3sigma .resize( nrxx * ((1==nspin)?1:12) ); + v4rho2sigma2.resize( nrxx * ((1==nspin)?1:15) ); + v4rhosigma3 .resize( nrxx * ((1==nspin)?1:20) ); + v4sigma4 .resize( nrxx * ((1==nspin)?1:15) ); + case 3: v3rho2sigma .resize( nrxx * ((1==nspin)?1:9) ); + v3rhosigma2 .resize( nrxx * ((1==nspin)?1:12) ); + v3sigma3 .resize( nrxx * ((1==nspin)?1:10) ); + case 2: v2rhosigma .resize( nrxx * ((1==nspin)?1:6) ); + v2sigma2 .resize( nrxx * ((1==nspin)?1:6) ); + case 1: vsigma .resize( nrxx * ((1==nspin)?1:3) ); + case 0: break; + default: throw std::domain_error("order ="+std::to_string(order) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + } + + 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 std::vector sgn = XC_Functional_Libxc::cal_sgn(rho_threshold, grho_threshold, func, nspin, nrxx, rho, sigma); + + // call libxc function + // attention: order 432 don't break + switch( func.info->family ) + { + case XC_FAMILY_LDA: + { + switch( order ) + { + case 4: xc_lda_lxc ( &func, nrxx, rho.data(), v4rho4.data() ); + case 3: xc_lda_kxc ( &func, nrxx, rho.data(), v3rho3.data() ); + case 2: xc_lda_fxc ( &func, nrxx, rho.data(), v2rho2.data() ); + case 1: xc_lda_exc_vxc( &func, nrxx, rho.data(), exc.data(), vrho.data() ); + break; + case 0: xc_lda_exc ( &func, nrxx, rho.data(), exc.data() ); + break; + default: throw std::domain_error("order ="+std::to_string(order) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + break; + } + case XC_FAMILY_GGA: + case XC_FAMILY_HYB_GGA: + { + switch( order ) + { + case 4: xc_gga_lxc ( &func, nrxx, rho.data(), sigma.data(), v4rho4.data(), v4rho3sigma.data(), v4rho2sigma2.data(), v4rhosigma3.data(), v4sigma4.data() ); + case 3: xc_gga_kxc ( &func, nrxx, rho.data(), sigma.data(), v3rho3.data(), v3rho2sigma.data(), v3rhosigma2.data(), v3sigma3.data() ); + case 2: xc_gga_fxc ( &func, nrxx, rho.data(), sigma.data(), v2rho2.data(), v2rhosigma.data(), v2sigma2.data() ); + case 1: xc_gga_exc_vxc( &func, nrxx, rho.data(), sigma.data(), exc.data(), vrho.data(), vsigma.data() ); + break; + case 0: xc_gga_exc ( &func, nrxx, rho.data(), sigma.data(), exc.data() ); + break; + default: throw std::domain_error("order ="+std::to_string(order) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + break; + } + default: + { + throw std::domain_error("func.info->family ="+std::to_string(func.info->family) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + } // end switch( func.info->family ) + } // end for( xc_func_type &func : funcs ) + + auto write_data = [&esolver]( + const std::string data_name, + const std::vector &data, + const int number_spin) + { + for(int is=0; ispw_big->bz, + esolver->pw_big->nbz, + esolver->pw_rhod->nplane * number_spin, + esolver->pw_rhod->startz_current, + data.data()+is, + esolver->pw_rhod->nx * esolver->pw_rhod->ny, + esolver->pw_rhod->nz, + number_spin, + esolver->pw_rhod->nz); + #else + if(nspin!=1) + { throw std::invalid_argument("nspin="+std::to_string(nspin)+" is invalid for ModuleIO::write_cube_core without MPI. see "+std:string(__FILE__)+" line "+std::to_string(__LINE__)); } + ModuleIO::write_cube_core( + ofs, + data.data()+is, + esolver->pw_rhod->nx * esolver->pw_rhod->ny, + esolver->pw_rhod->nz, + esolver->pw_rhod->nz); + #endif + } + }; + + write_data( "rho", rho, nspin ); + + if(1!=nspin && 2!=PARAM.inp.nspin) + write_data( "amag", amag, 1 ); + + if(is_gga) + write_data( "sigma", sigma, (1==nspin)?1:3 ); + + switch( order ) + { + case 4: write_data( "v4rho4", v4rho4, (1==nspin)?1:5 ); + case 3: write_data( "v3rho3", v3rho3, (1==nspin)?1:4 ); + case 2: write_data( "v2rho2", v2rho2, (1==nspin)?1:3 ); + case 1: write_data( "vrho" , vrho , nspin ); + case 0: write_data( "exc" , exc , 1 ); + break; + default: throw std::domain_error("order ="+std::to_string(order) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + if(is_gga) + { + switch( order ) + { + case 4: write_data( "v4rho3sigma" , v4rho3sigma , (1==nspin)?1:12 ); + write_data( "v4rho2sigma2", v4rho2sigma2, (1==nspin)?1:15 ); + write_data( "v4rhosigma3" , v4rhosigma3 , (1==nspin)?1:20 ); + write_data( "v4sigma4" , v4sigma4 , (1==nspin)?1:15 ); + case 3: write_data( "v3rho2sigma" , v3rho2sigma , (1==nspin)?1:9 ); + write_data( "v3rhosigma2" , v3rhosigma2 , (1==nspin)?1:12 ); + write_data( "v3sigma3" , v3sigma3 , (1==nspin)?1:10 ); + case 2: write_data( "v2rhosigma" , v2rhosigma , (1==nspin)?1:6 ); + write_data( "v2sigma2" , v2sigma2 , (1==nspin)?1:6 ); + case 1: write_data( "vsigma" , vsigma , (1==nspin)?1:3 ); + case 0: break; + default: throw std::domain_error("order ="+std::to_string(order) + +" unfinished in "+std::string(__FILE__)+" line "+std::to_string(__LINE__)); + break; + } + } + + XC_Functional_Libxc::finish_func(funcs); + + ModuleBase::timer::tick("ModuleIO","write_libxc_r"); +} + +#endif // USE_LIBXC \ No newline at end of file diff --git a/source/module_io/write_libxc_r.h b/source/module_io/write_libxc_r.h new file mode 100644 index 0000000000..13085da022 --- /dev/null +++ b/source/module_io/write_libxc_r.h @@ -0,0 +1,30 @@ +//====================== +// AUTHOR : Peize Lin +// DATE : 2024-09-12 +//====================== + +#ifndef WRITE_LIBXC_R_H +#define WRITE_LIBXC_R_H + +#ifdef USE_LIBXC + +#include + + class Charge; + namespace ModuleESolver{ class ESolver_FP; } + +namespace ModuleIO +{ + extern void write_libxc_r( + const int order, + const std::vector &func_id, + const int &nrxx, // number of real-space grid + const double &omega, // volume of cell + const double tpiba, + const Charge* const chr, + const ModuleESolver::ESolver_FP*const esolver); +} + +#endif // USE_LIBXC + +#endif // WRITE_LIBXC_R_H \ No newline at end of file diff --git a/source/module_parameter/input_parameter.h b/source/module_parameter/input_parameter.h index f870dd7ba0..cc2b2f4fac 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -314,6 +314,7 @@ struct Input_para int out_freq_ion = 0; ///< the frequency ( >= 0 ) of ionic step to output charge density; ///< 0: output only when ion steps are finished std::vector out_chg = {0, 3}; ///< output charge density. 0: no; 1: yes + std::vector out_xc_r = {-1, 3}; ///< output xc(r). -1: no; >=0: output the order of xc(r) int out_pot = 0; ///< yes or no int out_wfc_pw = 0; ///< 0: no; 1: txt; 2: dat bool out_wfc_r = false; ///< 0: no; 1: yes