diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index bfcd443bea..2c83ed749d 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -134,6 +134,7 @@ - [Output Variables](#variables-related-to-output-information) - [out\_freq\_elec](#out_freq_elec) - [out\_chg](#out_chg) + - [out\_xc\_r](#out_xc_r) - [out\_pot](#out_pot) - [out\_dm](#out_dmk) - [out\_dm1](#out_dmr) @@ -1652,6 +1653,25 @@ These variables are used to control the output of properties. - **Default**: 0 3 - **Note**: In the 3.10-LTS version, the file names are SPIN1_CHG.cube and SPIN1_CHG_INI.cube, etc. +### 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}`: + - 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 d65c057c9c..41bea656b7 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -539,6 +539,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\ ctrl_output_lcao.o\ diff --git a/source/module_io/CMakeLists.txt b/source/module_io/CMakeLists.txt index 70286ccbe7..a41c1b5179 100644 --- a/source/module_io/CMakeLists.txt +++ b/source/module_io/CMakeLists.txt @@ -31,6 +31,7 @@ list(APPEND objects write_mlkedf_descriptors.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/read_input_item_output.cpp b/source/module_io/read_input_item_output.cpp index c0a342b6c9..3d2f15dc64 100644 --- a/source/module_io/read_input_item_output.cpp +++ b/source/module_io/read_input_item_output.cpp @@ -514,6 +514,27 @@ void ReadInput::item_output() add_intvec_bcast(input.out_wfc_re_im, para.input.out_wfc_re_im.size(), 0); 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()); + }; + item.check_value = [](const Input_Item& item, const Parameter& para) { + if (para.input.out_xc_r[0] >= 0) + { +#ifndef USE_LIBXC + ModuleBase::WARNING_QUIT("ReadInput", "INPUT out_xc_r is only aviailable with Libxc"); +#endif + } + }; + sync_intvec(input.out_xc_r, 2, -1); + this->add_item(item); + } { Input_Item item("if_separate_k"); item.annotation = "specify whether to write the partial charge densities for all k-points to individual files " diff --git a/source/module_io/write_libxc_r.cpp b/source/module_io/write_libxc_r.cpp new file mode 100644 index 0000000000..076a8d02c9 --- /dev/null +++ b/source/module_io/write_libxc_r.cpp @@ -0,0 +1,448 @@ +//====================== +// 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_basis/module_pw/pw_basis_big.h" +#include "module_basis/module_pw/pw_basis.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 &chr, + const ModulePW::PW_Basis_Big &pw_big, + const ModulePW::PW_Basis &pw_rhod) +{ + 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 = [&pw_big, &pw_rhod]( + const std::string data_name, + const std::vector &data, + const int number_spin) + { + for(int is=0; is data_cube(nxyz, 0.0); + + // 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; + num_z[ip] += bz; + } + + // start_z: start position of z in + // processor ip. + std::vector start_z(nproc_in_pool, 0); + for (int ip = 1; ip < nproc_in_pool; ip++) + { + start_z[ip] = start_z[ip - 1] + num_z[ip - 1]; + } + + // which_ip: found iz belongs to which ip. + std::vector which_ip(nz, 0); + for (int iz = 0; iz < nz; iz++) + { + for (int ip = 0; ip < nproc_in_pool; ip++) + { + if (iz >= start_z[nproc_in_pool - 1]) + { + which_ip[iz] = nproc_in_pool - 1; + break; + } + else if (iz >= start_z[ip] && iz < start_z[ip + 1]) + { + which_ip[iz] = ip; + break; + } + } + } + + int count = 0; + std::vector zpiece(nxy, 0.0); + + // save the rho one z by one z. + for (int iz = 0; iz < nz; iz++) + { + zpiece.assign(nxy, 0.0); + + // tag must be different for different iz. + const int tag = iz; + MPI_Status ierror; + + // case 1: the first part of rho in processor 0. + if (which_ip[iz] == 0 && rank_in_pool == 0) + { + for (int ixy = 0; ixy < nxy; ixy++) + { + // 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) * nld]; + } + } + // case 2: > first part rho: send the rho to + // processor 0. + else if (which_ip[iz] == rank_in_pool) + { + for (int ixy = 0; ixy < nxy; ixy++) + { + zpiece[ixy] = data[ixy * nplane + (iz - startz_current) * nld]; + } + MPI_Send(zpiece.data(), nxy, MPI_DOUBLE, 0, tag, POOL_WORLD); + } + + // case 2: > first part rho: processor 0 receive the rho + // from other processors + else if (rank_in_pool == 0) + { + MPI_Recv(zpiece.data(), nxy, MPI_DOUBLE, which_ip[iz], tag, POOL_WORLD, &ierror); + } + + if (my_rank == 0) + { + /// for cube file + for (int ixy = 0; ixy < nxy; ixy++) + { + data_cube[ixy * nz + iz] = zpiece[ixy]; + } + /// for cube file + } + } // end iz + + // for cube file + if (my_rank == 0) + { + for (int ixy = 0; ixy < nxy; ixy++) + { + for (int iz = 0; iz < nz; iz++) + { + ofs_cube << " " << data_cube[ixy * nz + iz]; + if ((iz % n_data_newline == n_data_newline-1) && (iz != nz - 1)) + { + ofs_cube << "\n"; + } + } + ofs_cube << "\n"; + } + } + /// for cube file + } + MPI_Barrier(MPI_COMM_WORLD); +} + +#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++) + { + ofs_cube << " " << data[iz * nxy + ixy]; + // ++count_cube; + if ((iz % n_data_newline == n_data_newline-1) && (iz != nz - 1)) + { + ofs_cube << "\n"; + } + } + ofs_cube << "\n"; + } +} + +#endif // #ifdef __MPI + +#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..ad82e68b57 --- /dev/null +++ b/source/module_io/write_libxc_r.h @@ -0,0 +1,54 @@ +//====================== +// AUTHOR : Peize Lin +// DATE : 2024-09-12 +//====================== + +#ifndef WRITE_LIBXC_R_H +#define WRITE_LIBXC_R_H + +#ifdef USE_LIBXC + +#include +#include + + class Charge; + namespace ModulePW{ class PW_Basis_Big; } + namespace ModulePW{ class PW_Basis; } + +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 &chr, + const ModulePW::PW_Basis_Big &pw_big, + const ModulePW::PW_Basis &pw_rhod); + + #ifdef __MPI + extern void write_cube_core( + std::ofstream &ofs_cube, + const int bz, + const int nbz, + const int nplane, + const int startz_current, + const double*const data, + const int nxy, + const int nz, + const int nld, + const int n_data_newline); + #else + extern void write_cube_core( + std::ofstream &ofs_cube, + const double*const data, + const int nxy, + const int nz, + const int n_data_newline); + #endif +} + +#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 fd3e29aade..548e7da575 100644 --- a/source/module_parameter/input_parameter.h +++ b/source/module_parameter/input_parameter.h @@ -361,6 +361,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 int printe = 0; ///< Print out energy for each band for every printe step, default is scf_nmax diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 9e7e62b60b..5aebb6724e 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -19,6 +19,10 @@ #include "module_parameter/parameter.h" #include "module_cell/k_vector_utils.h" +#ifdef USE_LIBXC +#include "module_io/write_libxc_r.h" +#endif + namespace ModuleESolver { @@ -85,7 +89,7 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp) #ifdef __MPI this->pw_rho->initmpi(GlobalV::NPROC_IN_POOL, GlobalV::RANK_IN_POOL, POOL_WORLD); #endif - if (this->classname == "ESolver_OF" || PARAM.inp.of_ml_gene_data == 1) + if (this->classname == "ESolver_OF" || PARAM.inp.of_ml_gene_data == 1) { this->pw_rho->setfullpw(inp.of_full_pw, inp.of_full_pw_dim); } @@ -255,6 +259,22 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso &(ucell), PARAM.inp.out_elf[1]); } + +#ifdef USE_LIBXC + // 7) write xc(r) + if(PARAM.inp.out_xc_r[0]>=0) + { + ModuleIO::write_libxc_r( + PARAM.inp.out_xc_r[0], + XC_Functional::get_func_id(), + this->pw_rhod->nrxx, // number of real-space grid + ucell.omega, // volume of cell + ucell.tpiba, + this->chr, + *this->pw_big, + *this->pw_rhod); + } +#endif } }