diff --git a/docs/advanced/elec_properties/band.md b/docs/advanced/elec_properties/band.md index 2a0f7256a4..c33f5fc829 100644 --- a/docs/advanced/elec_properties/band.md +++ b/docs/advanced/elec_properties/band.md @@ -1,5 +1,7 @@ # Extracting Band Structure +For plane-wave screened hybrid calculations, follow the [hybrid NSCF workflow](hybrid_nscf.md), which additionally requires the frozen SCF orbitals and source q mesh. + In ABACUS, in order to obtain the eigenvalues of Hamiltonian, or generally called band structure, examples can be found in [examples/band](https://github.com/deepmodeling/abacus-develop/tree/develop/examples/12_band). Similar to the [DOS case](https://abacus-rtd.readthedocs.io/en/latest/advanced/elec_properties/dos.html), one first needs to perform a ground-state energy calculation ***with one additional keyword "[out_chg](https://abacus-rtd.readthedocs.io/en/latest/advanced/input_files/input-main.html#out-chg)" in the INPUT file***: diff --git a/docs/advanced/elec_properties/hybrid_nscf.md b/docs/advanced/elec_properties/hybrid_nscf.md new file mode 100644 index 0000000000..36caf3c2a1 --- /dev/null +++ b/docs/advanced/elec_properties/hybrid_nscf.md @@ -0,0 +1,69 @@ +# PW screened hybrid NSCF + +Hybrid NSCF reads the charge density, binary `wfk*_pw.dat` wavefunctions and +weighted occupations in `eig_occ.txt` from a converged SCF calculation. It keeps +that source fixed while solving a new K-point list or additional target bands. +No new file format is required. + +## SCF + +Use a complete, uniformly weighted source mesh and: + +```text +calculation scf +suffix hybrid +basis_type pw +dft_functional HSE +smearing_method fixed +symmetry -1 +kpar 1 +bndpar 1 +exxace false +exx_gamma_extra false +out_chg 1 +out_wfc_pw 2 +out_freq_ion 0 +cal_force 0 +cal_stress 0 +``` + +## NSCF + +Keep the same structure, pseudopotentials and cutoffs. In a separate output +directory, retain the shared settings above, replace `KPT` with the target +points, and change: + +```text +calculation nscf +read_file_dir ../scf/OUT.hybrid +init_wfc random +nbands 20 +out_band 1 +out_wfc_pw 0 +out_chg 0 +``` + +Supports CPU/CUDA and `nspin 1/2`; use one MPI rank for CUDA. ACE, SOC, +unscreened exchange, forces and stress are unsupported in this workflow. +Set `symmetry -1`: input `0` is accepted but PW EXX resets it to `-1` at runtime. + +HSE defaults to `exx_singularity_correction limits`, as in LCAO. Set `gygi` to +retain the historical PW correction, which can produce jumps at source K points. +Use matching SCF/NSCF settings for consistent bands. Do not use Methfessel-Paxton +smearing (`mp`/`mp2`/`mp3`) for the source: it can write occupations outside +`[0, k-point weight]`, which hybrid NSCF rejects; the non-negative schemes +(`fixed`, `gauss`, `mv`, `fd`) are accepted. Configuration differences +or missing `INPUT.info` produce warnings; missing or incompatible numerical +source data produces an error. + +## Verification + +With MPI and LibXC enabled, substitute your configured executable/build paths: + +```bash +cmake --build build --target MODULE_ESOLVER_pw_hybrid_source +OMP_NUM_THREADS=1 ctest --test-dir build -R '^MODULE_ESOLVER_pw_hybrid_source$' +# Frozen-source band regressions: tests/01_PW/nscf_hse_same and nscf_hse_path +cd tests/01_PW +OMP_NUM_THREADS=1 bash ../integrate/Autotest.sh -a -r '^nscf_hse' +``` diff --git a/docs/advanced/elec_properties/index.rst b/docs/advanced/elec_properties/index.rst index d75e28e28c..b7efcfe731 100644 --- a/docs/advanced/elec_properties/index.rst +++ b/docs/advanced/elec_properties/index.rst @@ -4,6 +4,7 @@ Electronic Properties and Outputs .. toctree:: band + hybrid_nscf dos Mulliken potential diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 43a4072889..83f1f49a47 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -334,6 +334,8 @@ OBJS_ESOLVER=esolver_factory.o\ esolver_of_tool.o\ esolver_dfpt_pw.o\ pw_others.o\ + pw_hybrid_nscf.o\ + pw_hybrid_source.o\ OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\ esolver_ks_lcao_tddft.o\ diff --git a/source/source_esolver/CMakeLists.txt b/source/source_esolver/CMakeLists.txt index 95488f4be3..059b0540a0 100644 --- a/source/source_esolver/CMakeLists.txt +++ b/source/source_esolver/CMakeLists.txt @@ -3,6 +3,8 @@ list(APPEND objects esolver_ks.cpp esolver_fp.cpp esolver_ks_pw.cpp + pw_hybrid_nscf.cpp + pw_hybrid_source.cpp esolver_ks_pw_tddft.cpp esolver_ks_lcaopw.cpp esolver_sdft_pw.cpp diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 7355d60a8c..592f505d1c 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -1,4 +1,5 @@ #include "esolver_ks_pw.h" +#include "pw_hybrid_nscf.h" #include "source_estate/elecstate_pw.h" #include "source_estate/module_charge/chg_symm.h" @@ -82,6 +83,8 @@ void ESolver_KS_PW::before_all_runners(BaseCell& basecell, const Inpu ESolver_KS::before_all_runners(ucell, inp); + validate_hybrid_nscf(inp, this->general_exx_info_); + //! setup and allocation for pelec, potentials, etc. elecstate::setup_estate_pw(ucell, this->kv, @@ -193,6 +196,13 @@ void ESolver_KS_PW::before_scf(UnitCell& ucell, const int istep) //! Setup EXX helper for Hamiltonian and psi exx_helper->before_scf(this->p_hamilt, this->stp.template get_psi_t(), *this->inp_, this->general_exx_info_); + if (this->inp_->calculation == "nscf" && this->general_exx_info_.cal_exx) + { + const Input_para& inp = *this->inp_; + const std::string& readin_dir = PARAM.globalv.global_readin_dir; + ModuleESolver::prepare_exx_nscf(inp, ucell, this->pw_wfc, this->exx_helper, readin_dir); + } + ModuleBase::timer::end("ESolver_KS_PW", "before_scf"); } @@ -315,7 +325,7 @@ void ESolver_KS_PW::iter_finish(UnitCell& ucell, const int istep, int // Related to EXX bool cal_exx = general_exx_info_.cal_exx; double hybrid_alpha = general_exx_info_.hybrid_alpha; - if (cal_exx && !exx_helper->get_op_first_iter()) + if (cal_exx && this->inp_->calculation != "nscf" && !exx_helper->get_op_first_iter()) { this->pelec->set_exx(exx_helper->cal_exx_energy(this->stp.template get_psi_t()), cal_exx, hybrid_alpha); } diff --git a/source/source_esolver/pw_hybrid_nscf.cpp b/source/source_esolver/pw_hybrid_nscf.cpp new file mode 100644 index 0000000000..3279698484 --- /dev/null +++ b/source/source_esolver/pw_hybrid_nscf.cpp @@ -0,0 +1,350 @@ +#include "pw_hybrid_nscf.h" +#include "esolver_ks_pw.h" +#include "pw_hybrid_source.h" +#include "source_io/module_wf/read_wfc_pw.h" +#include "source_base/module_out/filename.h" +#include "source_hamilt/module_xc/general_exx_info.h" +#include "source_pw/module_pwdft/exx_helper.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace ModuleESolver +{ +namespace +{ +// Check the stored Miller mapping before delegating coefficient loading to the unchanged IO reader. +void validate_source_grid(const std::string& filename, const ModulePW::PW_Basis_K& basis, int iq) +{ + std::ifstream in(filename, std::ios::binary); + const auto header = read_exx_source_header(in); + in.seekg(160); + int record = 0; + in.read(reinterpret_cast(&record), sizeof(record)); + if (record != 12LL * header.npw) + { + throw std::runtime_error("Invalid EXX source Miller record"); + } + const int grid[3] = {basis.nx, basis.ny, basis.nz}; + std::set indices; + for (int ig = 0; ig < header.npw; ++ig) + { + int miller[3] = {}; + ModuleBase::Vector3 direct_g; + for (int axis = 0; axis < 3; ++axis) + { + in.read(reinterpret_cast(&miller[axis]), sizeof(int)); + if (!in || miller[axis] < 0 || miller[axis] >= grid[axis]) + { + throw std::runtime_error("EXX source Miller index is incompatible with FFT grid"); + } + direct_g[axis] = miller[axis] <= grid[axis] / 2 ? miller[axis] : miller[axis] - grid[axis]; + } + const auto gk = direct_g * basis.G + basis.kvec_c[iq]; + const long long index = (static_cast(miller[0]) * basis.ny + miller[1]) * basis.nz + miller[2]; + if (gk.norm2() > basis.gk_ecut * (1.0 + 1e-10) || !indices.insert(index).second) + { + throw std::runtime_error("EXX source Miller mapping is incompatible with FFT basis"); + } + } +} + +bool source_parameter_matches(const std::string& key, std::string source_value, std::string target_value) +{ + if (key == "dft_functional") + { + const auto lowercase = [](unsigned char value) { return std::tolower(value); }; + std::transform(source_value.begin(), source_value.end(), source_value.begin(), lowercase); + std::transform(target_value.begin(), target_value.end(), target_value.begin(), lowercase); + } + if (source_value == target_value) + { + return true; + } + std::istringstream source(source_value); + std::istringstream target(target_value); + double source_number = 0; + double target_number = 0; + if (!(source >> source_number) || !(target >> target_number)) + { + return false; + } + std::string remaining; + if ((source >> remaining) || (target >> remaining)) + { + return false; + } + // INPUT.info may print doubles with the standard six significant digits. + const double scale = std::max(std::abs(source_number), std::abs(target_number)); + return std::isfinite(source_number) && std::isfinite(target_number) + && std::abs(source_number - target_number) <= 5e-6 * scale; +} + +void warn_source_configuration(const Input_para& inp, const std::string& readin_dir) +{ + std::map expected; + expected["dft_functional"] = inp.dft_functional; + expected["exx_singularity_correction"] = inp.exx_singularity_correction; + const std::vector>> terms = { + {"exx_erfc_alpha", inp.exx_erfc_alpha}, + {"exx_erfc_omega", inp.exx_erfc_omega}, + {"exx_fock_alpha", inp.exx_fock_alpha}}; + for (const auto& term : terms) + { + std::ostringstream value; + for (const auto& number : term.second) + { + // Zero-valued Fock coefficients are equivalent to omitting the term. + if (term.first == "exx_fock_alpha" && std::stod(number) == 0.0) + { + continue; + } + value << number << ' '; + } + expected[term.first] = value.str(); + } + std::ostringstream cutoff; + cutoff << inp.ecutexx; + expected["ecutexx"] = cutoff.str(); + expected["exx_gamma_extra"] = inp.exx_gamma_extra ? "1" : "0"; + std::ifstream input(readin_dir + "INPUT.info"); + std::string line; + std::map saved; + while (std::getline(input, line)) + { + line = line.substr(0, line.find('#')); + std::istringstream record(line); + std::string key; + record >> key; + if (expected.count(key)) + { + std::getline(record, saved[key]); + } + } + for (const auto& item : expected) + { + if (!saved.count(item.first)) + { + ModuleBase::WARNING("ESolver_KS_PW", "EXX source configuration is incomplete or unavailable in INPUT.info; " + "continuing with existing wavefunctions and eig_occ.txt. Verify the source SCF settings."); + return; + } + std::istringstream original(saved[item.first]); + std::istringstream current(item.second); + std::string source_value; + std::string target_value; + bool mismatch = false; + while (original >> source_value) + { + // Mirror the target-side normalization for saved zero Fock coefficients. + if (item.first == "exx_fock_alpha" && std::stod(source_value) == 0.0) + { + continue; + } + if (!(current >> target_value)) + { + mismatch = true; + break; + } + if (!source_parameter_matches(item.first, source_value, target_value)) + { + mismatch = true; + } + } + if (current >> target_value) + { + mismatch = true; + } + if (mismatch) + { + const std::string message = "EXX source configuration differs for " + item.first + + "; continuing with the frozen SCF source and current NSCF exchange settings. " + "Use matching settings for consistent SCF/NSCF bands."; + ModuleBase::WARNING("ESolver_KS_PW", message); + } + } +} +} + +void validate_hybrid_nscf(const Input_para& inp, const General_Exx_Info& info) +{ + if (inp.calculation != "nscf" || !info.cal_exx) + { + return; + } + bool has_fock = false; + const auto fock = info.coulomb_param.find(Conv_Coulomb_Pot_K::Coulomb_Type::Fock); + if (fock != info.coulomb_param.end()) + { + // An explicit zero coefficient (`exx_fock_alpha 0`) retains a Fock entry but is + // equivalent to omitting the term; only a nonzero coefficient is unscreened exchange. + for (const auto& term : fock->second) + { + if (std::stod(term.at("alpha")) != 0.0) + { + has_fock = true; + break; + } + } + } + bool compatible_device = inp.device == "cpu"; +#ifdef __CUDA + compatible_device = compatible_device || inp.device == "gpu"; +#endif + const bool compatible_execution = compatible_device && inp.kpar == 1 && inp.bndpar == 1 && inp.nspin != 4; + const bool compatible_exchange = !inp.exxace && !inp.exx_gamma_extra && !has_fock; + const bool compatible_targets = inp.symmetry != "1" && inp.init_wfc != "file" && !inp.cal_force && !inp.cal_stress; + if (!compatible_execution || !compatible_exchange || !compatible_targets) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", + "Hybrid NSCF currently requires CPU or CUDA GPU, kpar=bndpar=1, nspin=1/2, symmetry=-1/0, " + "exxace=false, exx_gamma_extra=false, screened exchange, " + "cal_force=cal_stress=false and fresh target wavefunctions"); + } +} + +template +void prepare_exx_nscf(const Input_para& inp, + const UnitCell& ucell, + const ModulePW::PW_Basis_K* pw_wfc, + Exx_HelperBase* exx_helper, + const std::string& readin_dir) +{ + ModuleBase::timer::start("ESolver_KS_PW", "prepare_exx_nscf"); + std::unique_ptr source_points(new K_Vectors); + ModuleBase::matrix source_weights; + warn_source_configuration(inp, readin_dir); + const int spin_mult = inp.nspin == 2 ? 2 : 1; + const std::vector first_index(1, 0); + const std::string first_filename = ModuleIO::filename_output(readin_dir, "wf", "pw", 0, + first_index, inp.nspin, spin_mult, 2, false, false, -1); + std::vector headers; + try + { + std::ifstream first(first_filename, std::ios::binary); + const auto header = read_exx_source_header(first); + const int nks = header.nks; + if (nks % spin_mult != 0) + { + throw std::runtime_error("Invalid EXX source spin dimensions"); + } + source_points->set_nks(nks); + source_points->set_nkstot(nks); + const int nq = nks / spin_mult; + source_points->set_nkstot_nospin(nq); + source_points->set_spin_mult(spin_mult); + source_points->kvec_d.resize(nks); + source_points->kvec_c.resize(nks); + source_points->wk.resize(nks); + source_points->isk.resize(nks); + source_points->ngk.resize(nks); + source_points->ik2iktot.resize(nks); + for (int iq = 0; iq < nks; ++iq) + { + source_points->ik2iktot[iq] = iq; + } + const ModuleBase::Matrix3 inverse_reciprocal = ucell.G.Inverse(); + for (int iq = 0; iq < nks; ++iq) + { + const std::string filename = ModuleIO::filename_output(readin_dir, "wf", "pw", iq, + source_points->ik2iktot, inp.nspin, nks, 2, false, false, -1); + std::ifstream in(filename, std::ios::binary); + const auto metadata = read_exx_source_header(in); + if (metadata.ecutwfc != inp.ecutwfc) + { + throw std::runtime_error("EXX source wavefunction cutoff differs from the current basis"); + } + headers.push_back(metadata); + auto& q = source_points->kvec_c[iq]; + q.x = metadata.kvec_c[0]; + q.y = metadata.kvec_c[1]; + q.z = metadata.kvec_c[2]; + source_points->kvec_d[iq] = q * inverse_reciprocal; + source_points->wk[iq] = metadata.weight; + source_points->isk[iq] = iq / nq; + } + std::ifstream occupations(readin_dir + "eig_occ.txt"); + read_exx_source_occupations(occupations, headers, inp.nspin, source_weights); + } + catch (const std::exception& error) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", error.what()); + } + const int nqs = source_points->get_nks(); + int nq = source_points->get_nkstot_nospin(); + source_points->para_k.kinfo(nq, 1, 0, pw_wfc->poolrank, pw_wfc->poolnproc, inp.nspin); + std::unique_ptr source_basis( + new ModulePW::PW_Basis_K(inp.device, inp.precision)); +#ifdef __MPI + source_basis->initmpi(pw_wfc->poolnproc, pw_wfc->poolrank, pw_wfc->pool_world); +#endif + source_basis->initgrids(ucell.lat0, ucell.latvec, + pw_wfc->nx, pw_wfc->ny, pw_wfc->nz); + source_basis->initparameters(false, inp.ecutwfc, nqs, source_points->kvec_d.data()); + // Preserve the exact stored Cartesian coordinates expected by the existing reader. + for (int iq = 0; iq < nqs; ++iq) + { + source_basis->kvec_c[iq] = source_points->kvec_c[iq]; + } + source_basis->fft_bundle.initfftmode(inp.fft_mode); + source_basis->setuptransform(); + source_basis->collect_local_pw(inp.erf_ecut, inp.erf_height, inp.erf_sigma); + for (int iq = 0; iq < nqs; ++iq) + { + source_points->ngk[iq] = source_basis->npwk[iq]; + } + const int nbands = source_weights.nc; + const int nbasis = source_basis->npwk_max; + std::unique_ptr> source_psi( + new psi::Psi(nqs, nbands, nbasis, source_points->ngk, true)); + for (int iq = 0; iq < nqs; ++iq) + { + const std::string wfc_filename = ModuleIO::filename_output(readin_dir, "wf", "pw", iq, + source_points->ik2iktot, inp.nspin, nqs, 2, false, false, -1); + try + { + validate_source_grid(wfc_filename, *source_basis, iq); + } + catch (const std::exception& error) + { + ModuleBase::WARNING_QUIT("ESolver_KS_PW", error.what()); + } + ModuleBase::ComplexMatrix wfc(nbands, nbasis); + ModuleIO::read_wfc_pw(wfc_filename, source_basis.get(), pw_wfc->poolrank, + pw_wfc->poolnproc, nbands, 1, iq, iq, nqs, wfc); + source_psi->fix_k(iq); + std::vector coefficients(nbands * nbasis); + for (int ib = 0; ib < nbands; ++ib) + { + for (int ig = 0; ig < nbasis; ++ig) + { + coefficients[ib * nbasis + ig] = static_cast(wfc(ib, ig)); + } + } + base_device::memory::synchronize_memory_op()( + source_psi->get_pointer(), coefficients.data(), coefficients.size()); + } + auto* helper = static_cast*>(exx_helper); + helper->op_exx->set_source(std::move(source_basis), std::move(source_points), + *source_psi, std::move(source_weights)); + ModuleBase::timer::end("ESolver_KS_PW", "prepare_exx_nscf"); +} + + +template void prepare_exx_nscf, base_device::DEVICE_CPU>( + const Input_para&, const UnitCell&, const ModulePW::PW_Basis_K*, Exx_HelperBase*, const std::string&); +template void prepare_exx_nscf, base_device::DEVICE_CPU>( + const Input_para&, const UnitCell&, const ModulePW::PW_Basis_K*, Exx_HelperBase*, const std::string&); +#if defined(__CUDA) || defined(__ROCM) +template void prepare_exx_nscf, base_device::DEVICE_GPU>( + const Input_para&, const UnitCell&, const ModulePW::PW_Basis_K*, Exx_HelperBase*, const std::string&); +template void prepare_exx_nscf, base_device::DEVICE_GPU>( + const Input_para&, const UnitCell&, const ModulePW::PW_Basis_K*, Exx_HelperBase*, const std::string&); +#endif +} // namespace ModuleESolver diff --git a/source/source_esolver/pw_hybrid_nscf.h b/source/source_esolver/pw_hybrid_nscf.h new file mode 100644 index 0000000000..91805d29a4 --- /dev/null +++ b/source/source_esolver/pw_hybrid_nscf.h @@ -0,0 +1,27 @@ +#ifndef PW_HYBRID_NSCF_H +#define PW_HYBRID_NSCF_H + +#include + +struct Input_para; +struct General_Exx_Info; +class UnitCell; +class Exx_HelperBase; +namespace ModulePW +{ +class PW_Basis_K; +} +namespace ModuleESolver +{ +void validate_hybrid_nscf(const Input_para& inp, const General_Exx_Info& info); + +// Read the frozen SCF ensemble in readin_dir and install it into the EXX operator +// held by exx_helper, which takes ownership for its whole remaining lifetime. +template +void prepare_exx_nscf(const Input_para& inp, + const UnitCell& ucell, + const ModulePW::PW_Basis_K* pw_wfc, + Exx_HelperBase* exx_helper, + const std::string& readin_dir); +} +#endif diff --git a/source/source_esolver/pw_hybrid_source.cpp b/source/source_esolver/pw_hybrid_source.cpp new file mode 100644 index 0000000000..fefcff3109 --- /dev/null +++ b/source/source_esolver/pw_hybrid_source.cpp @@ -0,0 +1,168 @@ +#include "pw_hybrid_source.h" + +#include "source_base/matrix.h" + +#include +#include +#include +#include +#include +#include + +namespace +{ +template +void read_binary(std::istream& in, T& value) +{ + in.read(reinterpret_cast(&value), sizeof(T)); + if (!in) + { + throw std::runtime_error("Missing or truncated EXX source wavefunction"); + } +} +} + +ModuleESolver::ExxSourceHeader ModuleESolver::read_exx_source_header(std::istream& in) +{ + ExxSourceHeader header{}; + int record = 0; + read_binary(in, record); + if (record != 72) + { + throw std::runtime_error("Invalid EXX source wavefunction header"); + } + read_binary(in, header.ik); + read_binary(in, header.nks); + for (double& coordinate : header.kvec_c) + { + read_binary(in, coordinate); + if (!std::isfinite(coordinate)) + { + throw std::runtime_error("Invalid EXX source wavefunction k point"); + } + } + read_binary(in, header.weight); + read_binary(in, header.npw); + read_binary(in, header.nbands); + read_binary(in, header.ecutwfc); + read_binary(in, header.lat0); + read_binary(in, header.tpiba); + read_binary(in, record); + if (record != 72 || header.ik <= 0 || header.nks < header.ik || header.nks > 100000 + || header.nbands <= 0 || header.nbands > 100000 || header.npw <= 0 + || !std::isfinite(header.weight) || header.weight <= 0 + || !std::isfinite(header.ecutwfc) || header.ecutwfc <= 0 + || !std::isfinite(header.lat0) || header.lat0 <= 0 + || !std::isfinite(header.tpiba) || header.tpiba <= 0) + { + throw std::runtime_error("Invalid EXX source wavefunction dimensions or metadata"); + } + // Two 72-byte records, a Miller-index record, and complex band records. + const long long expected_bytes = 160LL + 8 + 12LL * header.npw + + header.nbands * (8LL + 16LL * header.npw); + in.seekg(0, std::ios::end); + if (!in || in.tellg() != expected_bytes) + { + throw std::runtime_error("Truncated or invalid EXX source wavefunction size"); + } + return header; +} + +void ModuleESolver::read_exx_source_occupations(std::istream& in, + const std::vector& headers, + int nspin, + ModuleBase::matrix& weights) +{ + const int spin_mult = nspin == 2 ? 2 : 1; + const int nks = headers.size(); + if ((nspin != 1 && nspin != 2) || nks == 0 || nks % spin_mult != 0) + { + throw std::runtime_error("Invalid EXX source spin dimensions"); + } + const int nq = nks / spin_mult; + const int nbands = headers.front().nbands; + if (static_cast(nks) * nbands > 100000000) + { + throw std::runtime_error("Invalid EXX source occupation dimensions"); + } + std::string line; + std::getline(in, line); + std::getline(in, line); + std::getline(in, line); + int saved_spin = 0; + if (std::sscanf(line.c_str(), " Spin number %d", &saved_spin) != 1 || saved_spin != nspin) + { + throw std::runtime_error("Missing or incompatible EXX source eig_occ.txt"); + } + weights.create(nks, nbands); + const double expected_weight = nspin == 1 ? 2.0 / nq : 1.0 / nq; + for (int iq = 0; iq < nks; ++iq) + { + const auto& header = headers[iq]; + do + { + if (!std::getline(in, line)) + { + throw std::runtime_error("Truncated EXX source eig_occ.txt"); + } + } while (line.find_first_not_of(" \t\r") == std::string::npos); + int spin = 0; + int index = 0; + int total = 0; + int npw = 0; + double q[3] = {}; + const int fields = std::sscanf(line.c_str(), + " spin=%d k-point=%d/%d Cartesian=%lf %lf %lf (%d plane wave)", + &spin, &index, &total, &q[0], &q[1], &q[2], &npw); + if (fields != 7 || spin != iq / nq + 1 || index != iq % nq + 1 || total != nq + || npw != header.npw || header.ik != iq + 1 || header.nks != nks + || header.nbands != nbands || std::abs(header.weight - expected_weight) > 1e-10) + { + throw std::runtime_error("EXX source requires matching eig_occ.txt and a complete uniform q mesh"); + } + for (int axis = 0; axis < 3; ++axis) + { + if (nspin == 2 && iq >= nq && header.kvec_c[axis] != headers[iq - nq].kvec_c[axis]) + { + throw std::runtime_error("EXX source spin channels must use the same q mesh"); + } + // eig_occ.txt prints Cartesian coordinates with at least eight significant digits. + const double tolerance = 5e-8 * std::max(1.0, std::abs(header.kvec_c[axis])); + if (!std::isfinite(q[axis]) || std::abs(q[axis] - header.kvec_c[axis]) > tolerance) + { + throw std::runtime_error("EXX source eig_occ.txt k point differs from wavefunction"); + } + } + for (int ib = 0; ib < nbands; ++ib) + { + std::getline(in, line); + std::istringstream band(line); + double energy = 0; + band >> index >> energy >> weights(iq, ib); + const bool parsed = static_cast(band); + std::string extra; + const bool trailing = static_cast(band >> extra); + if (!parsed || trailing || index != ib + 1 || !std::isfinite(energy) + || !std::isfinite(weights(iq, ib))) + { + throw std::runtime_error("Invalid EXX source occupation in eig_occ.txt"); + } + if (weights(iq, ib) < 0 || weights(iq, ib) > expected_weight + 1e-10) + { + // Methfessel-Paxton smearing can write occupations below zero or above + // the k-point weight, and the EXX operator silently drops wg < 1e-12; + // only non-negative occupation schemes are supported. + throw std::runtime_error( + "EXX source eig_occ.txt occupation is outside [0, k-point weight]; " + "Methfessel-Paxton smearing (mp/mp2/mp3) can produce such occupations, " + "which hybrid NSCF does not support. Rerun the SCF source with " + "smearing_method fixed, gauss, mv or fd"); + } + } + } + in >> std::ws; + if (!in.eof()) + { + throw std::runtime_error("Unexpected trailing EXX source eig_occ.txt data; use a single SCF snapshot"); + } +} diff --git a/source/source_esolver/pw_hybrid_source.h b/source/source_esolver/pw_hybrid_source.h new file mode 100644 index 0000000000..98a67cdcb6 --- /dev/null +++ b/source/source_esolver/pw_hybrid_source.h @@ -0,0 +1,37 @@ +#ifndef PW_HYBRID_SOURCE_H +#define PW_HYBRID_SOURCE_H + +#include +#include + +namespace ModuleBase +{ +class matrix; +} + +namespace ModuleESolver +{ +// Metadata already present in the existing binary PW wavefunction format. +struct ExxSourceHeader +{ + int ik; + int nks; + int npw; + int nbands; + double kvec_c[3]; + double weight; + double ecutwfc; + double lat0; + double tpiba; +}; + +ExxSourceHeader read_exx_source_header(std::istream& in); + +// Read the existing, k-weighted occupation column without changing its format. +void read_exx_source_occupations(std::istream& in, + const std::vector& headers, + int nspin, + ModuleBase::matrix& weights); + +} +#endif diff --git a/source/source_esolver/test/CMakeLists.txt b/source/source_esolver/test/CMakeLists.txt index 2b7a04a122..7cc6b83c49 100644 --- a/source/source_esolver/test/CMakeLists.txt +++ b/source/source_esolver/test/CMakeLists.txt @@ -27,3 +27,9 @@ AddTest( ../../source_cell/module_neighlist/neighbor_list.cpp ../../source_io/module_output/output_log.cpp ) + +AddTest( + TARGET MODULE_ESOLVER_pw_hybrid_source + LIBS base device + SOURCES test_pw_hybrid_source.cpp ../pw_hybrid_source.cpp +) diff --git a/source/source_esolver/test/test_pw_hybrid_source.cpp b/source/source_esolver/test/test_pw_hybrid_source.cpp new file mode 100644 index 0000000000..a0f03632aa --- /dev/null +++ b/source/source_esolver/test/test_pw_hybrid_source.cpp @@ -0,0 +1,147 @@ +#include "source_esolver/pw_hybrid_source.h" +#include "source_base/matrix.h" + +#include "gtest/gtest.h" + +#include +#include +#include +#include + +namespace +{ +std::vector source_headers(int nspin) +{ + const int nks = nspin == 2 ? 4 : 2; + std::vector headers(nks); + for (int iq = 0; iq < nks; ++iq) + { + auto& header = headers[iq]; + header.ik = iq + 1; + header.nks = nks; + header.npw = 100; + header.nbands = 3; + header.kvec_c[0] = 0.123456789123456 * (iq % 2); + header.weight = nspin == 1 ? 1.0 : 0.5; + } + return headers; +} + +std::string eig_occ(int nspin) +{ + const auto headers = source_headers(nspin); + std::ostringstream out; + out << "1 # ionic step\n Electronic state energy (eV) and occupations\n Spin number " << nspin << '\n'; + for (int iq = 0; iq < headers.size(); ++iq) + { + const auto& header = headers[iq]; + out << " spin=" << iq / 2 + 1 << " k-point=" << iq % 2 + 1 + << "/2 Cartesian=" << std::setprecision(8) << header.kvec_c[0] + << " 0 0 (100 plane wave)\n"; + out << "1 -1 " << header.weight << "\n2 0 " << header.weight * 0.125 << "\n3 1 0\n\n"; + } + return out.str(); +} +} + +TEST(PWHybridSource, ReadsExistingWeightedOccupationsForBothSpins) +{ + for (int nspin = 1; nspin <= 2; ++nspin) + { + const auto headers = source_headers(nspin); + std::istringstream in(eig_occ(nspin)); + ModuleBase::matrix weights; + ModuleESolver::read_exx_source_occupations(in, headers, nspin, weights); + EXPECT_EQ(weights.nr, headers.size()); + EXPECT_EQ(weights.nc, 3); + EXPECT_DOUBLE_EQ(weights(1, 1), headers[1].weight * 0.125); + } +} + +TEST(PWHybridSource, RejectsInvalidExistingOccupationFiles) +{ + const auto headers = source_headers(1); + ModuleBase::matrix weights; + const std::string saved = eig_occ(1); + std::istringstream truncated(saved.substr(0, saved.size() - 6)); + EXPECT_THROW(ModuleESolver::read_exx_source_occupations(truncated, headers, 1, weights), std::runtime_error); + std::istringstream trailing(saved + saved); + EXPECT_THROW(ModuleESolver::read_exx_source_occupations(trailing, headers, 1, weights), std::runtime_error); + std::istringstream wrong_spin(saved); + EXPECT_THROW(ModuleESolver::read_exx_source_occupations(wrong_spin, headers, 2, weights), std::runtime_error); + auto changed = headers; + changed[1].kvec_c[0] += 0.1; + std::istringstream wrong_coordinates(saved); + EXPECT_THROW(ModuleESolver::read_exx_source_occupations(wrong_coordinates, changed, 1, weights), std::runtime_error); + std::string invalid = saved; + invalid.replace(invalid.find("0.125"), 5, "-0.25"); + std::istringstream negative(invalid); + try + { + // A negative occupation is what a Methfessel-Paxton source can produce; the + // diagnostic must name the unsupported smearing rather than report invalid data. + ModuleESolver::read_exx_source_occupations(negative, headers, 1, weights); + FAIL() << "Expected std::runtime_error"; + } + catch (const std::runtime_error& error) + { + EXPECT_NE(std::string(error.what()).find("Methfessel-Paxton"), std::string::npos); + EXPECT_NE(std::string(error.what()).find("smearing_method fixed"), std::string::npos); + } +} + +TEST(PWHybridSource, RejectsDifferentMeshesInSpinChannels) +{ + auto headers = source_headers(2); + headers[2].kvec_c[0] = 0.2; + std::string saved = eig_occ(2); + const std::string original = "spin=2 k-point=1/2 Cartesian=0"; + const std::size_t location = saved.find(original); + ASSERT_NE(location, std::string::npos); + saved.replace(location, original.size(), "spin=2 k-point=1/2 Cartesian=0.2"); + std::istringstream in(saved); + ModuleBase::matrix weights; + EXPECT_THROW(ModuleESolver::read_exx_source_occupations(in, headers, 2, weights), std::runtime_error); +} + +TEST(PWHybridSource, ReadsExistingBinaryWavefunctionHeader) +{ + // Native binary record layout from write_wfc_pw; no new wavefunction format. + std::ostringstream out(std::ios::binary); + const int record = 72; + const int ik = 1; + const int nks = 2; + const int npw = 1; + const int nbands = 3; + const double coordinate = 0.123456789123456; + const double weight = 1; + const double cutoff = 10; + const double lattice = 5; + const double tpiba = 1.256637061435917; + out.write(reinterpret_cast(&record), sizeof(record)); + out.write(reinterpret_cast(&ik), sizeof(ik)); + out.write(reinterpret_cast(&nks), sizeof(nks)); + for (int axis = 0; axis < 3; ++axis) + { + out.write(reinterpret_cast(&coordinate), sizeof(coordinate)); + } + out.write(reinterpret_cast(&weight), sizeof(weight)); + out.write(reinterpret_cast(&npw), sizeof(npw)); + out.write(reinterpret_cast(&nbands), sizeof(nbands)); + out.write(reinterpret_cast(&cutoff), sizeof(cutoff)); + out.write(reinterpret_cast(&lattice), sizeof(lattice)); + out.write(reinterpret_cast(&tpiba), sizeof(tpiba)); + out.write(reinterpret_cast(&record), sizeof(record)); + std::string saved = out.str(); + saved.resize(160 + 8 + 12 * npw + nbands * (8 + 16 * npw), '\0'); + std::istringstream in(saved, std::ios::binary); + const auto header = ModuleESolver::read_exx_source_header(in); + EXPECT_EQ(header.nbands, 3); + EXPECT_DOUBLE_EQ(header.kvec_c[0], coordinate); + EXPECT_DOUBLE_EQ(header.weight, weight); + std::istringstream truncated(saved.substr(0, saved.size() - 1), std::ios::binary); + EXPECT_THROW(ModuleESolver::read_exx_source_header(truncated), std::runtime_error); + saved[0] = 0; + std::istringstream invalid(saved, std::ios::binary); + EXPECT_THROW(ModuleESolver::read_exx_source_header(invalid), std::runtime_error); +} diff --git a/source/source_hamilt/module_xc/exx_info.cpp b/source/source_hamilt/module_xc/exx_info.cpp index 7098c1dd7d..98896c4490 100644 --- a/source/source_hamilt/module_xc/exx_info.cpp +++ b/source/source_hamilt/module_xc/exx_info.cpp @@ -114,7 +114,8 @@ bool init_general_exx_info(General_Exx_Info& info, const Input_para& inp) { info.coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc] = {{ {"alpha", ModuleBase::GlobalFunc::TO_STRING(erfc_alpha[i])}, - {"omega", ModuleBase::GlobalFunc::TO_STRING(inp.exx_erfc_omega[i])} }}; + {"omega", ModuleBase::GlobalFunc::TO_STRING(inp.exx_erfc_omega[i])}, + {"singularity_correction", inp.exx_singularity_correction} }}; } } } diff --git a/source/source_io/module_parameter/read_inp_exx_dftu.cpp b/source/source_io/module_parameter/read_inp_exx_dftu.cpp index bbecf5cace..308c4df63e 100644 --- a/source/source_io/module_parameter/read_inp_exx_dftu.cpp +++ b/source/source_io/module_parameter/read_inp_exx_dftu.cpp @@ -466,14 +466,19 @@ void ReadInput::item_exx() { Input_Item item("exx_singularity_correction"); item.annotation = "set the scheme of Coulomb singularity correction"; - item.category = "Exact Exchange (LCAO)"; + item.category = "Exact Exchange (Common)"; item.type = "String"; - item.description = R"(* spencer: see Phys. Rev. B 77, 193110 (2008). -* revised_spencer: see Phys. Rev. Mater. 5, 013807 (2021). Set the scheme of Coulomb singularity correction.)"; + item.description = R"(Scheme for Coulomb singularity / zero-transfer treatment in exact exchange. +* limits: use the finite screened Coulomb-kernel limit, without an auxiliary correction. Default for HSE and cwp22 in both PW and LCAO. PW requires screened exchange and exx_gamma_extra=false; the latter defaults to false for this scheme. +* gygi: retain the historical PW Gygi-Baldereschi auxiliary-function correction. PW only; default for unscreened PW hybrids. For independent-k screened hybrid NSCF this can cause discontinuities at source k points. +* spencer: LCAO spherical truncation, Phys. Rev. B 77, 193110 (2008). +* revised_spencer: LCAO, Phys. Rev. Mater. 5, 013807 (2021). +* massidda / carrier: LCAO auxiliary singularity corrections. +Use the same scheme for consistent SCF and NSCF bands. A source configuration mismatch produces a warning and does not prevent reading the existing SCF files.)"; item.default_value = "default"; item.unit = ""; read_sync_string(input.exx_singularity_correction); - item.reset_value = [](const Input_Item& item, Parameter& para) { + item.reset_value = [this](const Input_Item& item, Parameter& para) { if (para.input.exx_singularity_correction == "default") { std::string& dft_functional = para.input.dft_functional; @@ -490,13 +495,58 @@ void ReadInput::item_exx() || dft_functional_lower == "lrc_wpbeh" || dft_functional_lower == "cam_pbeh") { - para.input.exx_singularity_correction = "spencer"; + if (para.input.basis_type == "pw") + { + para.input.exx_singularity_correction = "gygi"; + } + else + { + para.input.exx_singularity_correction = "spencer"; + } } else if (dft_functional_lower == "hse" || dft_functional_lower == "cwp22") { para.input.exx_singularity_correction = "limits"; } } + if (para.input.basis_type == "pw" && para.input.exx_singularity_correction == "limits") + { + const auto gamma_item = std::find_if(this->input_lists.begin(), this->input_lists.end(), + [](const std::pair& entry) { + return entry.first == "exx_gamma_extra"; + }); + if (gamma_item != this->input_lists.end() && !gamma_item->second.is_read()) + { + para.input.exx_gamma_extra = false; + } + } + }; + item.check_value = [](const Input_Item& item, const Parameter& para) { + const std::string& scheme = para.input.exx_singularity_correction; + if (para.input.basis_type == "pw" && scheme != "default") + { + if (scheme != "limits" && scheme != "gygi") + { + ModuleBase::WARNING_QUIT("ReadInput", "PW exx_singularity_correction must be limits or gygi"); + } + if (scheme == "limits") + { + bool has_fock = false; + for (const auto& alpha : para.input.exx_fock_alpha) + { + has_fock = has_fock || std::stod(alpha) != 0.0; + } + if (has_fock || para.input.exx_gamma_extra) + { + ModuleBase::WARNING_QUIT("ReadInput", + "PW limits requires screened exchange and exx_gamma_extra=false"); + } + } + } + else if (para.input.basis_type != "pw" && scheme == "gygi") + { + ModuleBase::WARNING_QUIT("ReadInput", "exx_singularity_correction gygi is only supported for PW"); + } }; this->add_item(item); } diff --git a/source/source_io/module_parameter/read_inp_other.cpp b/source/source_io/module_parameter/read_inp_other.cpp index 2a99ac6022..68f17c0bdf 100644 --- a/source/source_io/module_parameter/read_inp_other.cpp +++ b/source/source_io/module_parameter/read_inp_other.cpp @@ -886,8 +886,8 @@ When false (default), both the direction and magnitude of the magnetic moment ar item.annotation = "whether to perform gamma extrapolation in exxpw"; item.category = "Exact Exchange (PW)"; item.type = "Boolean"; - item.description = "Whether to use the gamma point extrapolation method to calculate the Fock exchange operator. See https://doi.org/10.1103/PhysRevB.79.205114 for details. Should be set to true most of the time."; - item.default_value = "True"; + item.description = "Whether to use the gamma point extrapolation method to calculate the Fock exchange operator. See https://doi.org/10.1103/PhysRevB.79.205114 for details. Defaults to false for PW exx_singularity_correction=limits (finite screened kernel), and true otherwise. Explicit true is incompatible with the PW limits scheme."; + item.default_value = "False for PW limits; True otherwise"; item.unit = ""; read_sync_bool(input.exx_gamma_extra); this->add_item(item); diff --git a/source/source_io/module_parameter/read_inp_sys.cpp b/source/source_io/module_parameter/read_inp_sys.cpp index a0ba7ff398..68b04e8005 100644 --- a/source/source_io/module_parameter/read_inp_sys.cpp +++ b/source/source_io/module_parameter/read_inp_sys.cpp @@ -151,7 +151,7 @@ void ReadInput::item_system() item.description = R"(Specify the type of calculation. * scf: perform self-consistent electronic structure calculations -* nscf: perform non-self-consistent electronic structure calculations. A charge density file is required +* nscf: perform non-self-consistent electronic structure calculations using a saved charge density. PW screened hybrids also require SCF binary wavefunctions and eig_occ.txt; see [hybrid NSCF](../elec_properties/hybrid_nscf.md). * relax: perform structure relaxation calculations, the relax_nmax parameter depicts the maximal number of ionic iterations * cell-relax: perform cell relaxation calculations * md: perform molecular dynamics simulations diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 6107585e73..27a93accc7 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -1787,6 +1787,7 @@ TEST_F(InputTest, Item_test2) } { // exx_singluarity_correction auto it = find_label("exx_singularity_correction", readinput.input_lists); + TestParameters::input(param).basis_type = "lcao"; TestParameters::input(param).exx_singularity_correction = "default"; TestParameters::input(param).dft_functional = "HF"; it->second.reset_value(it->second, param); @@ -1801,6 +1802,27 @@ TEST_F(InputTest, Item_test2) TestParameters::input(param).dft_functional = "HSE"; it->second.reset_value(it->second, param); EXPECT_EQ(TestParameters::input(param).exx_singularity_correction, "limits"); + + TestParameters::input(param).basis_type = "pw"; + TestParameters::input(param).dft_functional = "PBE0"; + TestParameters::input(param).exx_singularity_correction = "default"; + it->second.reset_value(it->second, param); + EXPECT_EQ(TestParameters::input(param).exx_singularity_correction, "gygi"); + + TestParameters::input(param).dft_functional = "HSE"; + TestParameters::input(param).exx_singularity_correction = "default"; + TestParameters::input(param).exx_gamma_extra = true; + it->second.reset_value(it->second, param); + EXPECT_EQ(TestParameters::input(param).exx_singularity_correction, "limits"); + EXPECT_FALSE(TestParameters::input(param).exx_gamma_extra); + + auto gamma = find_label("exx_gamma_extra", readinput.input_lists); + gamma->second.str_values = {"true"}; + TestParameters::input(param).exx_gamma_extra = true; + it->second.reset_value(it->second, param); + EXPECT_TRUE(TestParameters::input(param).exx_gamma_extra); + gamma->second.str_values.clear(); + TestParameters::input(param).exx_gamma_extra = false; } { // exx_ccp_rmesh_times auto it = find_label("exx_ccp_rmesh_times", readinput.input_lists); diff --git a/source/source_pw/module_pwdft/exx_helper.cpp b/source/source_pw/module_pwdft/exx_helper.cpp index b183fc0e07..79d102becd 100644 --- a/source/source_pw/module_pwdft/exx_helper.cpp +++ b/source/source_pw/module_pwdft/exx_helper.cpp @@ -12,7 +12,7 @@ template void Exx_Helper::init(const UnitCell& ucell, const Input_para& inp, const ModuleBase::matrix& wg, const General_Exx_Info& exx_info) { if (inp.calculation != "scf" && inp.calculation != "relax" - && inp.calculation != "cell-relax" && inp.calculation != "md") + && inp.calculation != "cell-relax" && inp.calculation != "md" && inp.calculation != "nscf") { return; } @@ -22,7 +22,7 @@ void Exx_Helper::init(const UnitCell& ucell, const Input_para& inp, c return; } - if (exx_info.separate_loop) + if (exx_info.separate_loop && inp.calculation != "nscf") { XC_Functional::set_xc_first_loop(ucell); this->set_firstiter(); @@ -36,7 +36,7 @@ void Exx_Helper::before_scf(void* p_hamilt, void* psi, const Input_pa { /// Return if not a valid calculation type if (inp.calculation != "scf" && inp.calculation != "relax" - && inp.calculation != "cell-relax" && inp.calculation != "md") + && inp.calculation != "cell-relax" && inp.calculation != "md" && inp.calculation != "nscf") { return; } @@ -60,6 +60,11 @@ bool Exx_Helper::iter_finish(void* p_elec, Charge* p_charge, void* ps UnitCell& ucell, const Input_para& inp, bool& conv_esolver, int iter) { + if (inp.calculation == "nscf") + { + return false; + } + /// Return if EXX is not enabled if (op_exx == nullptr) { diff --git a/source/source_pw/module_pwdft/op_pw_exx.cpp b/source/source_pw/module_pwdft/op_pw_exx.cpp index a6057a1e2c..b580cad1f6 100644 --- a/source/source_pw/module_pwdft/op_pw_exx.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include namespace hamilt @@ -53,6 +54,8 @@ OperatorEXXPW::OperatorEXXPW(const int* isk_in, { ModuleBase::WARNING_QUIT("OperatorEXXPW", "EXX Calculation does not support k-point parallelism when exxace is set to false"); } + source_basis_ = wfcpw_in; + source_points_ = kv_in; gamma_extrapolation = exx_info.gamma_extrapolation; bool is_mp = kv_in->get_is_mp(); #ifdef __MPI @@ -88,8 +91,12 @@ OperatorEXXPW::OperatorEXXPW(const int* isk_in, Real tpiba2 = tpiba * tpiba; // initialize rhopw_dev on the resolved EXX cutoff - rhopw_dev = new ModulePW::PW_Basis(wfcpw->get_device(), rhopw->get_precision()); - rhopw_dev->fft_bundle.setfft(wfcpw->get_device(), rhopw->get_precision()); + // Pair densities have the wavefunction precision, while the ordinary + // charge density may retain double precision in a single-precision run. + const std::string exchange_device = wfcpw->get_device(); + const std::string exchange_precision = wfcpw->get_precision(); + rhopw_dev = new ModulePW::PW_Basis(exchange_device, exchange_precision); + rhopw_dev->fft_bundle.setfft(exchange_device, exchange_precision); #ifdef __MPI rhopw_dev->initmpi(rhopw->poolnproc, rhopw->poolrank, rhopw->pool_world); #endif @@ -115,14 +122,20 @@ OperatorEXXPW::OperatorEXXPW(const int* isk_in, auto param_erfc = this->coulomb_param[Conv_Coulomb_Pot_K::Coulomb_Type::Erfc]; for (auto param: param_erfc) { - erfc_div.push_back(exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, - std::stod(param["omega"]), - kv, - wfcpw, - rhopw_dev, - tpiba, - gamma_extrapolation, - ucell->omega)); + double correction = 0.0; + if (param.at("singularity_correction") == "gygi") + { + const double omega = std::stod(param.at("omega")); + correction = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, + omega, + kv, + wfcpw, + rhopw_dev, + tpiba, + gamma_extrapolation, + ucell->omega); + } + erfc_div.push_back(correction); } } // end of constructor @@ -225,6 +238,58 @@ void OperatorEXXPW::act(const int nbands, } } +template +void OperatorEXXPW::set_source(std::unique_ptr source_basis, + std::unique_ptr source_points, + const psi::Psi& source_psi, + ModuleBase::matrix source_weights) +{ + if (exxace_ || gamma_extrapolation || exx_sg_init) + { + ModuleBase::WARNING_QUIT("OperatorEXXPW", "Independent EXX sources require direct exchange without gamma extrapolation, before FFT initialization"); + } + if (source_basis->nks != source_points->get_nks() || source_psi.get_nk() != source_basis->nks + || source_weights.nr != source_basis->nks || source_weights.nc != source_psi.get_nbands() + || source_psi.get_nbasis() != source_basis->npwk_max + || source_basis->nx != wfcpw->nx || source_basis->ny != wfcpw->ny || source_basis->nz != wfcpw->nz + || source_basis->nrxx != wfcpw->nrxx || source_basis->poolnproc != wfcpw->poolnproc) + { + ModuleBase::WARNING_QUIT("OperatorEXXPW", "Incompatible EXX source dimensions or FFT distribution"); + } + owned_source_basis_ = std::move(source_basis); + owned_source_points_ = std::move(source_points); + owned_source_weights_ = std::move(source_weights); + source_basis_ = owned_source_basis_.get(); + source_points_ = owned_source_points_.get(); + wg = &owned_source_weights_; + psi = source_psi; + // Cache the Gygi-Baldereschi auxiliary correction on the frozen source mesh, rather than + // the independent target path. Limits keeps this correction at zero. + const auto screened = coulomb_param.find(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc); + if (screened != coulomb_param.end()) + { + for (std::size_t i = 0; i < screened->second.size(); ++i) + { + const auto& param = screened->second[i]; + erfc_div[i] = 0.0; + if (param.at("singularity_correction") == "gygi") + { + const double omega = std::stod(param.at("omega")); + erfc_div[i] = exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type::Erfc, + omega, + source_points_, + source_basis_, + rhopw_dev, + tpiba, + gamma_extrapolation, + ucell->omega); + } + } + } + q_points.clear(); + first_iter = false; +} + template void OperatorEXXPW::act_op(const int nbands, const int nbasis, @@ -245,7 +310,7 @@ void OperatorEXXPW::act_op(const int nbands, auto q_points = get_q_points(this->ik); int nk_fac = nspin_ == 2 ? 2 : 1; - int nk = wfcpw->nks / nk_fac; + int nk = source_basis_->nks / nk_fac; const Real nqs = q_points.size(); maybe_setup_exx_grid(); @@ -254,7 +319,8 @@ void OperatorEXXPW::act_op(const int nbands, for (int iq: q_points) { - get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, this->ik, iq % nk, false, this->coulomb_param); + const int iq_nospin = iq % nk; + get_exx_potential(source_points_, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, this->ik, iq_nospin, false, this->coulomb_param); for (int m_iband = 0; m_iband < psi.get_nbands(); m_iband++) { // occupation of the source state (m, iq), not of the target k-point @@ -264,10 +330,18 @@ void OperatorEXXPW::act_op(const int nbands, continue; } - wfc_to_real_exx(get_pw(m_iband, iq), iq, nbasis); + const T* source_band = get_pw(m_iband, iq); + if (source_basis_ == wfcpw) + { + wfc_to_real_exx(source_band, iq, psi.get_nbasis()); + } + else + { + source_basis_->recip_to_real(ctx, source_band, psi_mq_real, iq); + } // full accumulation weight, hybrid_alpha included - const Real factor = this->hybrid_alpha * wg_mqb_real / kv->wk[iq] / nqs; + const Real factor = this->hybrid_alpha * wg_mqb_real / source_points_->wk[iq] / nqs; apply_fock_all_bands(nbands, nbasis, iq, factor, tmhpsi); } // end of m_iband @@ -390,6 +464,11 @@ void OperatorEXXPW::maybe_setup_exx_grid() const return; } #endif + if (source_basis_ != wfcpw) + { + exx_sg_init = true; + return; + } setup_exx_small_grid(); if (std::is_same::value && !exx_sg_ok) { @@ -402,6 +481,10 @@ void OperatorEXXPW::maybe_setup_exx_grid() const template bool OperatorEXXPW::exx_grid_active() const { + if (source_basis_ != wfcpw) + { + return false; + } #if !defined(__CUDA) if (!std::is_same::value) { @@ -872,7 +955,7 @@ std::vector OperatorEXXPW::get_q_points(const int ik) const // if () // downsampling { - for (int iq = 0; iq < wfcpw->nks; iq++) + for (int iq = 0; iq < source_basis_->nks; iq++) { if (nspin_ ==1 ) { @@ -880,9 +963,7 @@ std::vector OperatorEXXPW::get_q_points(const int ik) const } else if (nspin_ == 2) { - int nk_fac = 2; - int nk = wfcpw->nks / nk_fac; - if (iq / nk == ik / nk) + if (source_points_->isk[iq] == isk[ik]) { q_points_ik.push_back(iq); } @@ -895,7 +976,7 @@ std::vector OperatorEXXPW::get_q_points(const int ik) const } // else // { - // for (int iq = 0; iq < wfcpw->nks; iq++) + // for (int iq = 0; iq < source_basis_->nks; iq++) // { // kv-> // } @@ -932,9 +1013,16 @@ template template OperatorEXXPW::OperatorEXXPW(const OperatorEXXPW *op) { + if (op->owned_source_basis_ || op->owned_source_points_) + { + ModuleBase::WARNING_QUIT("OperatorEXXPW", + "Precision conversion after EXX source installation is unsupported"); + } // copy all the datas this->isk = op->isk; this->wfcpw = op->wfcpw; + this->source_basis_ = op->source_basis_; + this->source_points_ = op->source_points_; this->rhopw = op->rhopw; this->rhopw_dev = op->rhopw_dev; this->psi = op->psi; @@ -1030,7 +1118,8 @@ double OperatorEXXPW::cal_exx_energy_op(psi::Psi *ppsi_) c for (int iq: q_points_ik) { - get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, ik, iq % nk, false, this->coulomb_param); + const int iq_nospin = iq % nk; + get_exx_potential(kv, wfcpw, rhopw_dev, pot, tpiba, gamma_extrapolation, ucell->omega, ik, iq_nospin, false, this->coulomb_param); for (int m_iband = 0; m_iband < nb; m_iband++) { const double wg_iqb_real = (*wg)(iq, m_iband); diff --git a/source/source_pw/module_pwdft/op_pw_exx.h b/source/source_pw/module_pwdft/op_pw_exx.h index 6e0131ab83..2715a49716 100644 --- a/source/source_pw/module_pwdft/op_pw_exx.h +++ b/source/source_pw/module_pwdft/op_pw_exx.h @@ -59,6 +59,14 @@ class OperatorEXXPW : public OperatorPW void set_wg(const ModuleBase::matrix *wg_in) { wg = wg_in; } + // The source orbitals remain fixed while target states are diagonalized. + // Takes ownership of the frozen ensemble: the operator keeps observer pointers + // to these objects, so their lifetime must match the operator's own. + void set_source(std::unique_ptr source_basis, + std::unique_ptr source_points, + const psi::Psi& source_psi, + ModuleBase::matrix source_weights); + void construct_ace() const; bool first_iter = true; @@ -71,6 +79,13 @@ class OperatorEXXPW : public OperatorPW private: const int* isk = nullptr; const ModulePW::PW_Basis_K* wfcpw = nullptr; + const ModulePW::PW_Basis_K* source_basis_ = nullptr; + const K_Vectors* source_points_ = nullptr; + // Owned frozen SCF source for hybrid NSCF; null in ordinary SCF, where the + // observer pointers above alias wfcpw/kv instead. + std::unique_ptr owned_source_basis_; + std::unique_ptr owned_source_points_; + ModuleBase::matrix owned_source_weights_; const ModulePW::PW_Basis* rhopw = nullptr; ModulePW::PW_Basis* rhopw_dev = nullptr; // for device const UnitCell *ucell = nullptr; diff --git a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp index 6094868def..6ae2a6d2c3 100644 --- a/source/source_pw/module_pwdft/op_pw_exx_pot.cpp +++ b/source/source_pw/module_pwdft/op_pw_exx_pot.cpp @@ -137,7 +137,7 @@ void get_exx_potential(const K_Vectors* kv, } } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk_fac = kv->get_spin_mult(); const int nk = nks / nk_fac; Real gg = (k_c - q_c + rhopw_dev->gcar[ig]).norm2() * tpiba2; @@ -167,9 +167,8 @@ void get_exx_potential(const K_Vectors* kv, double erfc_omega = std::stod(param["omega"]); double erfc_omega2 = erfc_omega * erfc_omega; double alpha = std::stod(param["alpha"]); - // Like the Fock correction, this global q sum is initialized by - // the operator. Recomputing it here would require idle pools to - // join a collective after their local k-point bounds check. + // Initialized on the SCF mesh, or replaced with the frozen source + // correction by set_source. The shared limits scheme stores zero. const double exx_div = OperatorEXXPW, Device>::erfc_div[i]; const ModuleBase::Vector3 k_c = wfcpw->kvec_c[ik]; const ModuleBase::Vector3 k_d = wfcpw->kvec_d[ik]; @@ -205,7 +204,7 @@ void get_exx_potential(const K_Vectors* kv, } } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk_fac = kv->get_spin_mult(); const int nk = nks / nk_fac; // const int ig_kq = ik * nks * npw + iq * npw + ig; @@ -324,7 +323,7 @@ void get_exx_stress_potential(const K_Vectors* kv, } } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk_fac = kv->get_spin_mult(); const int nk = nks / nk_fac; // const int ig_kq = ik * nks * npw + iq * npw + ig; @@ -391,7 +390,7 @@ void get_exx_stress_potential(const K_Vectors* kv, } } - const int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk_fac = kv->get_spin_mult(); const int nk = nks / nk_fac; // const int ig_kq = ik * nks * npw + iq * npw + ig; @@ -442,7 +441,7 @@ double exx_divergence(Conv_Coulomb_Pot_K::Coulomb_Type coulomb_type, double nqs_half2 = 0.5 * kv->nmp[1]; double nqs_half3 = 0.5 * kv->nmp[2]; - int nk_fac = PARAM.inp.nspin == 2 ? 2 : 1; + const int nk_fac = kv->get_spin_mult(); // here we follow the exx_divergence subroutine in q-e (PW/src/exx_base.f90) double alpha = 10.0 / wfcpw->gk_ecut; diff --git a/tests/01_PW/CASES_CPU.txt b/tests/01_PW/CASES_CPU.txt index b0d6cade3e..b6c569b958 100644 --- a/tests/01_PW/CASES_CPU.txt +++ b/tests/01_PW/CASES_CPU.txt @@ -123,5 +123,7 @@ relax_dftu_spin2 relax_dftu_spin4 scf_deltaspin2 scf_deltaspin4 +nscf_hse_same +nscf_hse_path 221_PW_PCC0D_H3O 222_PW_PCC2D_H3O diff --git a/tests/01_PW/CASES_GPU.txt b/tests/01_PW/CASES_GPU.txt index fe156d9365..224c940ced 100644 --- a/tests/01_PW/CASES_GPU.txt +++ b/tests/01_PW/CASES_GPU.txt @@ -115,3 +115,5 @@ scf_brav_bco_smear_mv #relax_dftu_spin4 scf_deltaspin2 scf_deltaspin4 +nscf_hse_same +nscf_hse_path diff --git a/tests/01_PW/nscf_hse_path/INPUT b/tests/01_PW/nscf_hse_path/INPUT new file mode 100644 index 0000000000..907cc77a8b --- /dev/null +++ b/tests/01_PW/nscf_hse_path/INPUT @@ -0,0 +1,32 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation nscf +nbands 3 +symmetry -1 + +#Parameters (2.Iteration) +ecutwfc 10 +scf_thr 1e-9 +scf_nmax 100 +pw_diag_thr 1e-10 +pw_seed 1 + +#Parameters (3.Basis) +basis_type pw + +#Parameters (4.Smearing) +smearing_method fixed + +#Parameters (5.ExactExchange) +dft_functional HSE +exxace false +exx_gamma_extra false + +read_file_dir ../nscf_hse_same/ +pseudo_dir ../../PP_ORB +out_band 1 +out_chg 0 +out_wfc_pw 0 +cal_force 0 +cal_stress 0 diff --git a/tests/01_PW/nscf_hse_path/KPT b/tests/01_PW/nscf_hse_path/KPT new file mode 100644 index 0000000000..56f407f00f --- /dev/null +++ b/tests/01_PW/nscf_hse_path/KPT @@ -0,0 +1,6 @@ +K_POINTS +3 +Direct +0 0 0 1 +0.125 0 0 1 +0.25 0 0 1 diff --git a/tests/01_PW/nscf_hse_path/README b/tests/01_PW/nscf_hse_path/README new file mode 100644 index 0000000000..748334c0a6 --- /dev/null +++ b/tests/01_PW/nscf_hse_path/README @@ -0,0 +1,2 @@ +PW HSE NSCF on the frozen SCF source committed in ../nscf_hse_same, solved on +an independent target k-point list that shares only the Gamma point. diff --git a/tests/01_PW/nscf_hse_path/STRU b/tests/01_PW/nscf_hse_path/STRU new file mode 100644 index 0000000000..c35003c65b --- /dev/null +++ b/tests/01_PW/nscf_hse_path/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +H 1 H_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.889766 + +LATTICE_VECTORS +0.0 2.708337 2.708337 +2.708337 0.0 2.708337 +2.708337 2.708337 0.0 + +ATOMIC_POSITIONS +Direct + +H +0 +2 +0.125 0.125 0.125 0 0 0 +0.675 0.835 0.875 0 0 0 diff --git a/tests/01_PW/nscf_hse_path/band.txt.ref b/tests/01_PW/nscf_hse_path/band.txt.ref new file mode 100644 index 0000000000..6c586d1df5 --- /dev/null +++ b/tests/01_PW/nscf_hse_path/band.txt.ref @@ -0,0 +1,3 @@ + 1 0.00000000 -11.12767527 0.23237114 5.91382475 + 2 0.03997035 -10.53435615 -0.56032894 5.88787937 + 3 0.07994070 -9.30092037 -2.20463726 5.45376065 diff --git a/tests/01_PW/nscf_hse_path/result.ref b/tests/01_PW/nscf_hse_path/result.ref new file mode 100644 index 0000000000..d1062bfbc3 --- /dev/null +++ b/tests/01_PW/nscf_hse_path/result.ref @@ -0,0 +1,4 @@ +etotref -34.96111659019459 +etotperatomref -17.4805582951 +CompareBand_pass 0 +totaltimeref 0.09 diff --git a/tests/01_PW/nscf_hse_same/INPUT b/tests/01_PW/nscf_hse_same/INPUT new file mode 100644 index 0000000000..a81d44ca26 --- /dev/null +++ b/tests/01_PW/nscf_hse_same/INPUT @@ -0,0 +1,33 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation nscf +nbands 5 +symmetry -1 + +#Parameters (2.Iteration) +ecutwfc 10 +scf_thr 1e-9 +scf_nmax 100 +pw_diag_thr 1e-10 +pw_seed 1 + +#Parameters (3.Basis) +basis_type pw + +#Parameters (4.Smearing) +smearing_method fixed + +#Parameters (5.ExactExchange) +dft_functional HSE +exx_fock_alpha 0 +exxace false +exx_gamma_extra false + +read_file_dir ./ +pseudo_dir ../../PP_ORB +out_band 1 +out_chg 0 +out_wfc_pw 0 +cal_force 0 +cal_stress 0 diff --git a/tests/01_PW/nscf_hse_same/INPUT.info b/tests/01_PW/nscf_hse_same/INPUT.info new file mode 100644 index 0000000000..a6818d20fd --- /dev/null +++ b/tests/01_PW/nscf_hse_same/INPUT.info @@ -0,0 +1,612 @@ +INPUT_PARAMETERS +#Parameters (1.System) +suffix hybrid #the name of main output directory +ntype 1 #atom species number +cell_replica 1 1 1 #replicate the input structure along the three lattice vectors +calculation scf #scf; relax; md; cell-relax; nscf; get_s; get_wf; get_pchg; gen_bessel; gen_opt_abfs; test_memory; test_neighbour +socket_driver 0 #run as a socket client for external drivers using the i-PI protocol +esolver_type ksdft #the energy solver: ksdft, sdft, ofdft, tdofdft, tddft, lj, dp, ks-lr, lr, dfpt +symmetry -1 #the control of symmetry +symmetry_prec 1e-06 #accuracy for symmetry +symmetry_autoclose 1 #whether to close symmetry automatically when error occurs in symmetry analysis +cal_force 0 #if calculate the force at the end of the electronic iteration +kpar 1 #devide all processors into kpar groups and k points will be distributed among +bndpar 1 #divide each k-point pool into band-parallel groups +latname user_defined_lattice #the name of lattice name +assume_isolated none #isolated-system electrostatic correction +init_wfc atomic #atomic; atomic+random; random; nao; nao+random; file txt; file binary +init_chg atomic #start charge is from 'atomic' or file +init_vel 0 #read velocity from STRU or not +mem_saver 0 #Only for nscf calculations. if set to 1, then a memory saving technique will be used for many k point calculations. +cal_stress 0 #calculate the stress or not +diago_proc 1 #the number of procs used to do diagonalization +nbspline -1 #the order of B-spline basis +kspacing 0 0 0 #unit in 1/bohr, should be > 0, default is 0 which means read KPT file +koffset 0 0 0 #offset for kspacing-generated automatic k-point mesh +kmesh_type gamma #mesh type for kspacing-generated automatic k-point mesh +min_dist_coef 0.2 #factor related to the allowed minimum distance between two atoms +device cpu #the computing device for ABACUS +precision double #the computing precision for ABACUS +gint_precision double #the computing precision for LCAO grid integral +timer_enable_nvtx 0 #enable NVTX labeling for profiling or not +cell_factor 1.2 #used in the construction of the pseudopotential tables +dm_to_rho 0 #reads dmr in npz format and calculates electron density +chg_extrap atomic #atomic; first-order; second-order; dm:coefficients of SIA +ecutwfc 10 #energy cutoff for wave functions +ecutrho 40 #energy cutoff for charge density and potential +nx 0 #number of points along x axis for FFT grid +ny 0 #number of points along y axis for FFT grid +nz 0 #number of points along z axis for FFT grid +ndx 0 #number of points along x axis for FFT smooth grid +ndy 0 #number of points along y axis for FFT smooth grid +ndz 0 #number of points along z axis for FFT smooth grid +pw_seed 1 #random seed for initializing wave functions +diag_subspace 0 #method of subspace diagonalization in dav_subspace. 0:LaPack; 1:genelpa, 2:scalapack +erf_ecut 0 #the value of the constant energy cutoff +fft_mode 0 #mode of FFTW +erf_height 0 #the height of the energy step for reciprocal vectors +erf_sigma 0.1 #the width of the energy step for reciprocal vectors +stru_file STRU #the filename of file containing atom positions +kpoint_file KPT #the name of file containing k points +pseudo_dir /tmp/abacus-hybrid-nscf-pr/tests/PP_ORB/ #the directory containing pseudo files +orbital_dir #the directory containing orbital files +read_file_dir OUT.hybrid/ #directory of files for reading +restart_load 0 #restart from disk +basis_type pw #PW; LCAO in pw; LCAO + +#Parameters (2.Electronic structure) +ks_solver cg #cg; dav; lapack; genelpa; elpa; scalapack_gvx; cusolver +nbands 3 #number of bands +nelec 0 #input number of electrons +nelec_delta 0 #change in the number of total electrons +nupdown 0 #the difference number of electrons between spin-up and spin-down +dft_functional HSE #exchange correlation functional +xc_temperature 0 #temperature for finite temperature functionals +xc_exch_ext #customize Libxc exchange functional parameters +xc_corr_ext #customize Libxc correlation functional parameters +pseudo_rcut 15 #default #exchange correlation functional +pseudo_mesh 0 #0: use our own mesh to do radial renormalization; 1: use mesh as in QE +nspin 1 #1: single spin; 2: up and down spin; 4: noncollinear spin +gga_grad 0 #Noncollinear GGA gradient method: 0 original, 1 local axis, 2 regularized projected LCA +smearing_method fixed #type of smearing_method: gauss; fd; fixed; mp; mp2; mv +smearing_sigma 0.015 #energy range for smearing +mixing_type broyden #plain; pulay; broyden +mixing_beta 0.8 #mixing parameter: 0 means no new charge +mixing_beta_mag -10 #mixing parameter for magnetic density +mixing_ndim 8 #mixing dimension in pulay or broyden +mixing_restart 0 #threshold to restart mixing during SCF +mixing_dmr 0 #whether to mix real-space density matrix +mixing_gg0 1 #mixing parameter in kerker +mixing_gg0_mag 0 #mixing parameter in kerker +mixing_gg0_min 0.1 #the minimum kerker coefficient +mixing_angle -10 #angle mixing parameter for non-colinear calculations +mixing_tau 0 #whether to mix tau in mGGA calculation +mixing_dftu 0 #whether to mix occ_mat in DFT+U calculation +gamma_only 0 #Only for localized orbitals set and gamma point. If set to 1, a fast algorithm is used +scf_nmax 100 #number of electron iterations +scf_thr 1e-09 #charge density error +scf_ene_thr -1 #total energy error threshold +scf_thr_type 1 #type of the criterion of scf_thr, 1: reci drho for pw, 2: real drho for lcao +scf_os_stop 0 #whether to stop scf when oscillation is detected +scf_os_thr -0.01 #charge density threshold for oscillation +scf_os_ndim 8 #number of old iterations used for oscillation detection +sc_os_ndim 5 #number of old iterations used for oscillation detection, for Spin-Constrained DFT +lspinorb 0 #consider the spin-orbit interaction +noncolin 0 #using non-collinear-spin +soc_lambda 1 #The fraction of SOC based on scalar relativity (SR) of the pseudopotential +dfthalf_type 0 #DFT-1/2 type, 0:off; 1:shell DFT-1/2 +pw_diag_thr 1e-10 #threshold for eigenvalues is cg electron iterations +diago_smooth_ethr 0 #smooth ethr for iter methods +use_k_continuity 0 #whether to use k-point continuity for initializing wave functions +pw_diag_nmax 50 #max iteration number for cg +pw_diag_ndim 4 #dimension of workspace for Davidson diagonalization +diago_cg_prec 1 #diago_cg_prec + +#Parameters (3.LCAO) +nb2d 0 #matrix 2d division +lmaxmax 2 #maximum of l channels used +lcao_ecut 0 #energy cutoff for LCAO +lcao_dk 0.01 #delta k for 1D integration in LCAO +lcao_dr 0.01 #delta r for 1D integration in LCAO +lcao_rmax 30 #max R for 1D two-center integration table +search_radius -1 #input search radius (Bohr) +bx 1 #division of an element grid in FFT grid along x +by 1 #division of an element grid in FFT grid along y +bz 1 #division of an element grid in FFT grid along z +elpa_num_thread -1 #Number of threads need to use in elpa +num_stream 4 #the nstream in compute the LCAO with CUDA +bessel_nao_ecut 10.000000 #energy cutoff for spherical bessel functions(Ry) +bessel_nao_tolerence 1e-12 #tolerence for spherical bessel root +bessel_nao_rcut #radial cutoff for spherical bessel functions(a.u.) +bessel_nao_smooth 1 #spherical bessel smooth or not +bessel_nao_sigma 0.1 #spherical bessel smearing_sigma + +#Parameters (4.Relaxation) +relax_method cg 2 #cg; bfgs; sd; cg_bfgs; lbfgs +relax_scale_force 0.5 #controls the size of the first CG 2 step +relax_nmax 1 #number of ion iteration steps +relax_cg_thr 0.5 #threshold for switching from cg to bfgs, unit: eV/Angstrom +force_thr 0.001 #force threshold, unit: Ry/Bohr +force_thr_ev 0.0257112 #force threshold, unit: eV/Angstrom +force_zero_out 0 #force invalid threshold, unit: eV/Angstrom +relax_bfgs_w1 0.01 #wolfe condition 1 for bfgs +relax_bfgs_w2 0.5 #wolfe condition 2 for bfgs +relax_bfgs_rmax 0.2 #maximal trust radius, unit: Bohr +relax_bfgs_rmin 1e-05 #minimal trust radius, unit: Bohr +relax_bfgs_init 0.5 #initial trust radius, unit: Bohr +stress_thr 0.5 #stress threshold +press1 0 #target pressure, unit: KBar +press2 0 #target pressure, unit: KBar +press3 0 #target pressure, unit: KBar +fixed_axes None #which axes are fixed +fixed_ibrav 0 #whether to preseve lattice type during relaxation +fixed_atoms 0 #whether to preseve direct coordinates of atoms during relaxation +relax 0 #allow relaxation along the specific direction + +#Parameters (5.Molecular dynamics) +md_type nvt #choose ensemble +md_nstep 10 #md steps +md_dt 1 #time step +md_thermostat nhc #choose thermostat +md_tfirst -1 #temperature first +md_tlast -1 #temperature last +md_prec_level 0 #precision level for vc-md +md_restart 0 #whether restart +md_restartfreq 5 #The period to output MD restart information +md_dumpfreq 1 #The period to dump MD information +md_neighbor_skin 0 #extra MDCell neighbor-list radius +md_out_force 1 #output all atomic forces into running_md.log or not +dump_force 1 #output atomic forces into the file MD_dump or not +dump_vel 1 #output atomic velocities into the file MD_dump or not +dump_virial 1 #output lattice virial into the file MD_dump or not +md_seed -1 #random seed for MD +md_tfreq 0 #oscillation frequency, used to determine qmass of NHC +md_tchain 1 #number of Nose-Hoover chains +md_pmode iso #NPT ensemble mode: iso, aniso, tri +ref_cell_factor 1 #construct a reference cell bigger than the initial cell +md_pcouple xyz #whether couple different components: xyz, xy, yz, xz, none +md_pfirst -1 #initial target pressure +md_plast -1 #final target pressure +md_pfreq 0 #oscillation frequency, used to determine qmass of thermostats coupled with barostat +md_pchain 1 #num of thermostats coupled with barostat +lj_rule 2 #combination rules used to construct the parameter matrix for LJ potential +lj_eshift 0 #whether to use energy shift for LJ potential +lj_rcut #cutoff radius of LJ potential +lj_epsilon #the value of epsilon for LJ potential +lj_sigma #the value of sigma for LJ potential +pot_file graph.pb #the filename of potential files for CMD such as DP +dp_rescaling 1 #rescaling factor for dp potential +dp_fparam #the frame parameter for dp potential +dp_aparam #the atomic parameter for dp potential +msst_direction 2 #the direction of shock wave +msst_vel 0 #the velocity of shock wave +msst_vis 0 #artificial viscosity +msst_tscale 0.01 #reduction in initial temperature +msst_qmass 1 #mass of thermostat +md_damp 1 #damping parameter (time units) used to add force in Langevin method +md_csvr_tau 100 #CSVR thermostat characteristic time scale +md_tolerance 100 #tolerance for velocity rescaling (K) +md_nraise 1 #parameters used when md_type=nvt +cal_syns 0 8 #calculate asynchronous overlap matrix to output for Hefei-NAMD +dmax 0.01 #maximum displacement of all atoms in one step (bohr) +plumed 0 #use PLUMED for collective variables and enhanced sampling +plumed_file plumed.dat #input file of PLUMED + +#Parameters (6.orbital free density functional theory) +of_kinetic wt #kinetic energy functional, such as tf, vw, wt +of_method tn #optimization method used in OFDFT, including cg1, cg2, tn (default) +of_conv energy #the convergence criterion, potential, energy (default), or both +of_tole 1e-06 #tolerance of the energy change (in Ry) for determining the convergence, default=2e-6 Ry +of_tolp 1e-05 #tolerance of potential for determining the convergence, default=1e-5 in a.u. +of_tf_weight 1 #weight of TF KEDF +of_vw_weight 1 #weight of vW KEDF +of_wt_alpha 0.833333 #parameter alpha of WT KEDF +of_wt_beta 0.833333 #parameter beta of WT KEDF +of_extwt_kappa 0.832442 #parameter kappa of EXT-WT KEDF +of_wt_rho0 0 #the average density of system, used in WT KEDF, in Bohr^-3 +of_hold_rho0 0 #If set to 1, the rho0 will be fixed even if the volume of system has changed, it will be set to 1 automaticly if of_wt_rho0 is not zero +of_lkt_a 1.3 #parameter a of LKT KEDF +of_xwm_rho_ref 0 #The reference density of XWM KEDF +of_xwm_kappa 0 #The parameter kappa of XWM KEDF +of_read_kernel 0 #If set to 1, the kernel of WT KEDF will be filled from file of_kernel_file, not from formula. Only usable for WT KEDF +of_kernel_file WTkernel.txt #The name of WT kernel file. +of_full_pw 1 #If set to 1, ecut will be ignored when collect planewaves, so that all planewaves will be used +of_full_pw_dim 0 #If of_full_pw = true, dimention of FFT is testricted to be (0) either odd or even; (1) odd only; (2) even only +of_ml_gene_data 0 #Generate training data or not +of_ml_device cpu #Run NN on GPU or CPU +of_ml_feg 0 #The Free Electron Gas limit: 0: no, 3: yes +of_ml_nkernel 1 #Number of kernels +of_ml_kernel 1 #Type of kernel, 1 for wt, 2 for yukawa, and 3 for TKK +of_ml_kernel_scaling 1 #Scaling parameter of kernel, w(r-r') = scaling^3 * w(scaling (r-r')) +of_ml_yukawa_alpha 1 #Parameter alpha of yukawa kernel +of_ml_kernel_file none #The file of TKK +of_ml_gamma 0 #Descriptor: gamma = (rho / rho0)^(1/3) +of_ml_p 0 #Descriptor: p = |nabla rho|^2 / [2 (3 pi^2)^(1/3) rho^(4/3)]^2 +of_ml_q 0 #Descriptor: q = nabla^2 rho / [4 (3 pi^2)^(2/3) rho^(5/3)] +of_ml_tanhp 0 #Descriptor: tanhp = tanh(chi_p * p) +of_ml_tanhq 0 #Descriptor: tanhq = tanh(chi_q * q) +of_ml_chi_p 1 #Hyperparameter: tanhp = tanh(chi_p * p) +of_ml_chi_q 1 #Hyperparameter: tanhq = tanh(chi_q * q) +of_ml_gammanl 0 #Descriptor: gammanl = int{gamma(r') * w(r-r') dr'} +of_ml_pnl 0 #Descriptor: pnl = int{p(r') * w(r-r') dr'} +of_ml_qnl 0 #Descriptor: qnl = int{q(r') * w(r-r') dr'} +of_ml_xi 0 #Descriptor: xi = int{rho(r')^(1/3) * w(r-r') dr'} / rho^(1/3) +of_ml_tanhxi 0 #Descriptor: tanhxi = tanh(chi_xi * xi) +of_ml_tanhxi_nl 0 #Descriptor: tanhxi_nl = int{tanhxi(r') * w(r-r') dr'} +of_ml_tanh_pnl 0 #Descriptor: tanh_pnl = tanh(chi_pnl * pnl) +of_ml_tanh_qnl 0 #Descriptor: tanh_qnl = tanh(chi_qnl * qnl) +of_ml_tanhp_nl 0 #Descriptor: tanhp_nl = int{tanhp(r') * w(r-r') dr'} +of_ml_tanhq_nl 0 #Descriptor: tanhq_nl = int{tanhq(r') * w(r-r') dr'} +of_ml_chi_xi 1 #Hyperparameter: tanhpxi = tanh(chi_xi * xi) +of_ml_chi_pnl 1 #Hyperparameter: tanh_pnl = tanh(chi_pnl * pnl) +of_ml_chi_qnl 1 #Hyperparameter: tanh_qnl = tanh(chi_qnl * qnl) +of_ml_local_test 0 #Test: read in the density, and output the F and Pauli potential +ml_exx 0 #Use ML EXX or not + +#Parameters (7.Stochastic DFT) +method_sto 2 #1: slow and save memory, 2: fast and waste memory +nbands_sto 256 #number of stochstic orbitals +nche_sto 100 #Chebyshev expansion orders +emin_sto 0 #trial energy to guess the lower bound of eigen energies of the Hamitonian operator +emax_sto 0 #trial energy to guess the upper bound of eigen energies of the Hamitonian operator +seed_sto 0 #the random seed to generate stochastic orbitals +initsto_ecut 0 #maximum ecut to init stochastic bands +initsto_freq 0 #frequency to generate new stochastic orbitals when running md +npart_sto 1 #Reduce memory when calculating Stochastic DOS + +#Parameters (8.DeepKS) +deepks_out_labels 0 #>0 compute descriptor for deepks. 1 used during training, 2 used for label production +deepks_out_freq_elec 0 #>0 frequency of electronic iteration to output descriptors and labels for deepks. +deepks_out_base none #base functional for output files, with dft_functional as target functional +deepks_scf 0 #>0 add V_delta to Hamiltonian +deepks_equiv 0 #whether to use equivariant version of DeePKS +deepks_model None #file dir of traced pytorch model: 'model.ptg +bessel_descriptor_lmax 2 #lmax used in generating spherical bessel functions +bessel_descriptor_ecut 10.000000 #energy cutoff for spherical bessel functions(Ry) +bessel_descriptor_tolerence 1e-12 #tolerence for spherical bessel root +bessel_descriptor_rcut 6 #radial cutoff for spherical bessel functions(a.u.) +bessel_descriptor_smooth 1 #spherical bessel smooth or not +bessel_descriptor_sigma 0.1 #sphereical bessel smearing_sigma +deepks_bandgap 0 #>0 for bandgap label +deepks_band_range -1 0 #(int, int) range of bands for bandgap label +deepks_v_delta 0 #!=0 for v_delta/v_delta_R label. when output, 1 for vdpre, 2 for phialpha and grad_evdm (can save memory ) -1 for vdrpre, -2 for phialpha_r and gevdm (can save memory) +deepks_out_unittest 0 #if set 1, prints intermediate quantities that shall be used for making unit test +lin_solver gmres #linear solver for real-time propagation +lin_precond kinetic_recycle #right preconditioner for PW real-time propagation +lin_thr 0 #residual tolerance of the linear solver +lin_maxiter 500 #maximum iterations per linear solve +lin_gmres_restart 20 #Arnoldi dimension of restarted GMRES +td_cn_init 1 #CN subspace initial guess and residual reuse +lin_reconstruct 1 #explicit GMRES residual reconstruction +estep_per_md 1 #steps of force change + +#Parameters (9.rt-tddft) +td_dt 1 #time step for evolving wavefunction +td_edm 0 #the method to calculate the energy density matrix +td_print_eij -1 #print eij or not +td_propagator 0 #method of propagator +td_vext 0 #add extern potential or not +td_vext_dire 1 #extern potential direction +td_stype 0 #type of electric field in space domain +td_ttype 0 #type of electric field in time domain +td_tstart 1 # number of steps where electric field starts +td_tend 1000 #number of steps where electric field ends +td_lcut1 0.05 #cut1 of interval in length gauge +td_lcut2 0.95 #cut2 of interval in length gauge +td_gauss_freq 22.13 #frequency (freq) of Gauss type electric field +td_gauss_phase 0 #phase of Gauss type electric field +td_gauss_sigma 30 #sigma of Gauss type electric field +td_gauss_t0 100 #step number of time center (t0) of Gauss type electric field +td_gauss_amp 0.25 #amplitude of Gauss type electric field +td_trape_freq 1.6 #frequency of Trapezoid type electric field +td_trape_phase 0 #phase of Trapezoid type electric field +td_trape_t1 1875 #t1 of Trapezoid type electric field +td_trape_t2 5625 #t2 of Trapezoid type electric field +td_trape_t3 7500 #t3 of Trapezoid type electric field +td_trape_amp 2.74 #amplitude of Trapezoid type electric field +td_trigo_freq1 1.16466 #frequency 1 of Trigonometric type electric field +td_trigo_freq2 0.029116 #frequency 2 of Trigonometric type electric field +td_trigo_phase1 0 #phase 1 of Trigonometric type electric field +td_trigo_phase2 0 #phase 2 of Trigonometric type electric field +td_trigo_amp 2.74 #amplitude of Trigonometric type electric field +td_heavi_t0 100 #t0 of Heaviside type electric field +td_heavi_amp 1 #amplitude of Heaviside type electric field +td_supsine_amp 0.27 #carrier electric-field scale of the supersine pulse +td_supsine_freq 0.18737 #carrier frequency of the supersine pulse +td_supsine_phase 0 #carrier phase at the center of the supersine pulse +td_supsine_sigma 0.75 #shape parameter of the supersine envelope +td_supsine_tstart default #start boundary step of the supersine pulse +td_supsine_tend default #end boundary step of the supersine pulse +init_vecpot_file 0 #init vector potential through file or not +ocp 0 #change occupation or not +ocp_set #set occupation +of_cd 0 #add CD Potential or not +of_mcd_alpha 1 #parameter of modified CD Potential +xc_kernel LDA #exchange correlation (XC) kernel for LR-TDDFT +lr_init_xc_kernel default #The method to initalize the xc kernel +lr_solver dav #the eigensolver for LR-TDDFT +lr_thr 0.01 #convergence threshold of the LR-TDDFT eigensolver +nocc 3 #the number of occupied orbitals to form the 2-particle basis ( <= nelec/2) +nvirt 1 #the number of virtual orbitals to form the 2-particle basis (nocc + nvirt <= nbands) + +#Parameters (10.lr-tddft) +lr_nstates 1 #the number of 2-particle states to be solved +lr_unrestricted 0 #Whether to use unrestricted construction for LR-TDDFT +abs_wavelen_range 0 0 #the range of wavelength(nm) to output the absorption spectrum +out_wfc_lr 0 #whether to output the eigenvectors (excitation amplitudes) in the particle-hole basis +abs_gauge velocity #whether to use length or velocity gauge to calculate the absorption spectrum in LR-TDDFT +abs_broadening 0.01 #the broadening (eta) for LR-TDDFT absorption spectrum +bse_tda tda #whether Tamm-Dancoff Approximation is used (can be 'tda', 'full' or 'both') +bse_spin_types singlet triplet #which spin type is calculated (can be 'singlet', 'triplet', also for test 'rpa', 'ipa') +bse_mem_save 0 #whether to save memory by adding V and W to BSE matrix directly +bse_ri_hartree 1 #whether to use RI approximation for Hartree term in BSE +bse_use_fine_kgrid 0 #Fine k-grid mode for BSE +bse_q_approx_mode 0 #q->kpair mapping mode: 0=exact, 1=coarse q grid, 2=mixed, 3=truncate +bse_q_approx_threshold 0.1 #threshold radius (in unit of 2*pi/lat0) for exact q mapping in mode 2; in mode 3 pairs with larger |q| are dropped entirely +out_bse_ab 0 #whether to output the AB matrix to file +bse_continue 0 #which step to continue from previous BSE calculation +plot_istate 0 #which state of exciton to be ploted +exciton_plot_type average #exciton density plot type: average or conditional +exciton_plot_format cube #exciton plot format: average supports cube, slice, or both; conditional supports slice only +exciton_fixed_coordinate 0 0 0 0 0 0 #fixed hole and electron Cartesian coordinates (Bohr): hx hy hz ex ey ez +exciton_slice_plane ab #cross-section plane: ab, bc, or ca +exciton_slice_pos 0 #offset along perpendicular direction (Bohr) for slice +exciton_slice_npoints 200 #grid points per dimension for slice +exciton_slice_range -1 2 -1 2 #slice cell range: ustart uend vstart vend +dfpt_qmesh 1 1 1 #Monkhorst-Pack q mesh for DFPT, should be >= 1, and commensurate with the KPT mesh +dfpt_qfile #file containing the DFPT q points; empty means the dfpt_qmesh Monkhorst-Pack mesh +dfpt_compute_q0 0 #whether to compute epsilon_inf and Born effective charges at q = 0 +dfpt_loto 0 #whether to apply the LO-TO non-analytic correction at q = 0 +dfpt_conv_thr 1e-08 #convergence threshold of the DFPT first-order density +dfpt_max_iter 100 #max number of the DFPT first-order density mixing iterations +dfpt_mix_beta 0.4 #mixing coefficient of the DFPT first-order density +out_freq_ion 0 #print information every few ionic steps +out_freq_td 0 #print information every few completed electronic iterations in RT-TDDFT +out_freq_elec 100 #print information every few electronic steps +out_chg 1 3 #> 0 output charge density for selected electron steps, second parameter controls the precision, default is 3. +out_pot 0 8 #output real space potential (with precision 8) +out_dmk 0 8 #>0 output density matrix DM(k) for each k-point +out_dmr 0 8 #output density matrix DM(R) with respect to lattice vector R (with precision 8) +out_wfc_pw 2 #output wave functions +out_wfc_lcao 0 #ouput LCAO wave functions, 0, no output 1: text, 2: binary +out_dos 0 #output energy and dos +out_ldos 0 3 #output mode of local density of states, second parameter controls the precision +out_band 1 8 #output energy and band structure (with precision 8) +out_proj_band 0 #output projected band structure + +#Parameters (11.Output) +out_stru 0 #output the structure files per ion step +out_level ie #ie(for electrons); i(for ions); +out_hsk 0 8 #output H(k) and S(k) matrices in reciprocal space +out_hsr 0 8 #output H(R) and S(R) matrices in real space +out_mat_tk 0 8 #output kinetic matrix of electrons T(k) +out_mat_r 0 8 #output r(R) matrix +out_mat_t 0 8 #output T(R) matrix +out_mat_dh 0 8 #output Hamiltonian derivatives dH/dR matrices +out_mat_dh_t 0 8 #output kinetic energy dH/dR (dT/dR) matrices +out_mat_dh_vl 0 8 #output local pseudopotential dH/dR (dV^L/dR) matrices +out_mat_dh_vnl 0 8 #output nonlocal pseudopotential dH/dR (dV^NL/dR) matrices +out_mat_dh_vh 0 8 #output Hartree dH/dR (dV^H/dR) matrices +out_mat_dh_vxc 0 8 #output exchange-correlation dH/dR (dV^XC/dR) matrices +out_mat_h_t 0 8 #output kinetic energy T(R) matrix +out_mat_h_vnl 0 8 #output nonlocal pseudopotential Vnl(R) matrix +out_mat_h_vl 0 8 #output local pseudopotential Vl(R) matrix +out_mat_h_vh 0 8 #output Hartree Vh(R) matrix +out_mat_h_vxc 0 8 #output exchange-correlation Vxc(R) matrix +out_mat_h_exx 0 8 #output exact-exchange Vexx(R) matrix +out_mat_ds 0 8 #output of derivative of S(R) matrix +out_mat_xc 0 #output exchange-correlation matrix in KS-orbital representation +out_mat_xc2 0 8 #output exchange-correlation matrix in NAO representation +out_mat_l 0 8 #output the expectation values of angular momentum operators +out_xc_r -1 3 #if >=0, output the derivatives of exchange correlation in realspace, second parameter controls the precision +out_eband_terms 0 #output the band energy terms separately +out_hr_npz 0 #output H(R) matrix in npz format +out_dm_npz 0 #output DM(R) matrix in npz format +out_mul 0 #mulliken charge or not +out_app_flag 1 #whether output r(R), H(R), S(R), T(R), and dH(R) matrices in an append manner during MD +out_ndigits 8 #the length of decimal part of output data +out_element_info 0 #output (projected) wavefunction of each element +restart_save 0 #print to disk every step for restart +rpa 0 #true:generate output files used in rpa calculation; false:(default) +rpa_out_vel 0 #whether to output velocity matrix for librpa +rpa_outdir OUT.librpa/ #output directory for librpa +out_pchg #specify the bands to be calculated for the partial (band-decomposed) charge densities +out_wfc_norm #specify the bands to be calculated for the norm of wavefunctions +out_wfc_re_im #specify the bands to be calculated for the real and imaginary parts of wavefunctions +if_separate_k 0 #specify whether to write the partial charge densities for all k-points to individual files or merge them +out_elf 0 3 #> 0 output electron localization function (ELF) for selected electron steps, second parameter controls the precision, default is 3. +out_spillage 0 #output spillage of wavefunctions. This parameter only accepts 0 or 2. +out_dipole 0 #output dipole or not +out_current 0 #output current or not +out_current_k 0 #output current for each k +out_efield 0 #output dipole or not +out_vecpot 0 #output TDDFT vector potential or not +cal_symm_repr 0 3 #output matrix representation of symmetry operation into running log file > 0 output the matrix representation of symmetry operation , the second parameter controls the precision, default is 3. +spillage_outdir ./ #output directory for spillage of wavefunctions. +dos_edelta_ev 0.03 #delta energy for dos +dos_sigma 0.07 #gauss b coefficeinet(default=0.07) +dos_scale 0.01 #scale dos range by + +#Parameters (12.Postprocess) +dos_emin_ev -15 #minimal range for dos +dos_emax_ev 15 #maximal range for dos +dos_nche 100 #orders of Chebyshev expansions for dos +stm_bias 1 0.1 1 #bias voltage used to calculate ldos +ldos_line 0 0 0 0 0 1 100 #start and end point of the line (direct coordinates) and number of points +cal_cond 0 #calculate electronic conductivities +cond_che_thr 1e-08 #control the error of Chebyshev expansions for conductivities +cond_dw 0.1 #frequency interval for conductivities +cond_wcut 10 #cutoff frequency (omega) for conductivities +cond_dt 0.02 #t interval to integrate Onsager coefficiencies +cond_dtbatch 0 #exp(iH*dt*cond_dtbatch) is expanded with Chebyshev expansion +cond_smear 1 #Smearing method for conductivities +cond_fwhm 0.4 #FWHM for conductivities +cond_nonlocal 1 #Nonlocal effects for conductivities +cond_mgga_vel 1 #Meta-GGA velocity correction for conductivities +berry_phase 0 #calculate berry phase or not +gdir 3 #calculate the polarization in the direction of the lattice vector +towannier90 0 #use wannier90 code interface or not +nnkpfile seedname.nnkp #the wannier90 code nnkp file name +wannier_method 1 #different implementation methods under Lcao basis set +wannier_spin up #calculate spin in wannier90 code interface +out_wannier_mmn 1 #output .mmn file or not +out_wannier_amn 1 #output .amn file or not +out_wannier_eig 1 #output .eig file or not +out_wannier_unk 0 #output UNK. file or not +out_wannier_wvfn_formatted 1 #output UNK. file in text format or in binary format + +#Parameters (13.Model) +efield_flag 0 #add electric field +dip_cor_flag 0 #dipole correction +efield_dir 2 #the direction of the electric field or dipole correction +efield_pos_max -1 #position of the maximum of the saw-like potential along crystal axis efield_dir +efield_pos_dec -1 #zone in the unit cell where the saw-like potential decreases +efield_amp 0 #amplitude of the electric field +gate_flag 0 #compensating charge or not +zgate 0.5 #position of charged plate +block 0 #add a block potential or not +block_down 0.45 #low bound of the block +block_up 0.55 #high bound of the block +block_height 0.1 #height of the block +imp_sol 0 #calculate implicit solvation correction or not +eb_k 80 #the relative permittivity of the bulk solvent +tau 1.0798e-05 #the effective surface tension parameter +sigma_k 0.6 #the width of the diffuse cavity +nc_k 0.00037 #the cut-off charge density + +#Parameters (14.vdW Correction) +vdw_method none #the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj ; d4) +vdw_d4_xc default #functional name passed to DFT-D4 +vdw_d4_model d4 #DFT-D4 dispersion model +vdw_s6 default #scale parameter of d2/d3_0/d3_bj +vdw_s8 default #scale parameter of d3_0/d3_bj +vdw_a1 default #damping parameter of d3_0/d3_bj +vdw_a2 default #damping parameter of d3_0/d3_bj +vdw_d 20 #damping parameter of d2 +vdw_abc 0 #third-order term? +vdw_c6_file default #filename of C6 +vdw_c6_unit Jnm6/mol #unit of C6, Jnm6/mol or eVA6 +vdw_r0_file default #filename of R0 +vdw_r0_unit A #unit of R0, A or Bohr +vdw_cutoff_type radius #expression model of periodic structure, radius or period +vdw_cutoff_radius 0 #radius cutoff for periodic structure +vdw_radius_unit Bohr #unit of radius cutoff for periodic structure +vdw_cutoff_width2 0.05 #smooth cutoff width for two-body dispersion +vdw_cutoff_width3 0 #smooth cutoff width for three-body dispersion +vdw_cutoff_period 3 3 3 #periods of periodic structure +vdw_cn_thr 40 #radius cutoff for cn +vdw_cn_thr_unit Bohr #unit of cn_thr, Bohr or Angstrom + +#Parameters (15.exx) +exx_fock_alpha #fraction of Fock exchange 1/r in hybrid functionals +exx_erfc_alpha 0.25 #fraction of exchange erfc(wr)/r in hybrid functionals +exx_erfc_omega 0.11 #range-separation parameter erfc(wr)/r in hybrid functionals +exx_separate_loop 1 #if 1, a two-step method is employed, else it will start with a GGA-Loop, and then Hybrid-Loop +exx_hybrid_step 100 #the maximal electronic iteration number in the evaluation of Fock exchange +exx_mixing_beta 1 #mixing_beta for outer-loop when exx_separate_loop=1 +exx_fock_lambda #used to compensate for divergence points at G=0 in the evaluation of Fock exchange using lcao_in_pw method +exx_pca_threshold 0.0001 #threshold to screen on-site ABFs in exx +exx_c_threshold 0.0001 #threshold to screen C matrix in exx +exx_cs_inv_thr -1 #threshold to inverse Vq in abfs for generating Cs +exx_v_threshold 0.1 #threshold to screen C matrix in exx +exx_dm_threshold 0.0001 #threshold to screen density matrix in exx +exx_c_grad_threshold 0.0001 #threshold to screen nabla C matrix in exx +exx_v_grad_threshold 0.1 #threshold to screen nabla V matrix in exx +exx_c_grad_r_threshold 0.0001 #threshold to screen nabla C * R matrix in exx +exx_v_grad_r_threshold 0.1 #threshold to screen nabla V * R matrix in exx +exx_ccp_rmesh_times 1.5 #how many times larger the radial mesh required for calculating Columb potential is to that of atomic orbitals +exx_opt_orb_lmax 0 #the maximum l of the spherical Bessel functions for opt ABFs +exx_opt_orb_ecut 0 #the cut-off of plane wave expansion for opt ABFs +exx_opt_orb_tolerence 0 #the threshold when solving for the zeros of spherical Bessel functions for opt ABFs +exx_real_number 0 #exx calculated in real or complex +exx_singularity_correction limits #set the scheme of Coulomb singularity correction +rpa_ccp_rmesh_times 10 #how many times larger the radial mesh required for calculating Columb potential is to that of atomic orbitals +exx_symmetry_realspace 0 #whether to reduce real-space sector in Hexx calculation +out_ri_cv 0 #Whether to output the coefficient tensor C and ABFs-representation Coulomb matrix V +out_unshrinked_v 0 #whether to output the large Vq matrix in unshrinked auxiliary basis +exx_coul_moment 0 #whether to use moment method for Coulomb calculation +exx_rotate_abfs 0 #whether to rotate auxiliary basis for Coulomb calculation +exx_multip_moments_threshold 1e-10 #threshold to screen multipole moments in Coulomb calculation +shrink_abfs_pca_thr -1 #threshold to shrink auxiliary basis for GW/RPA +shrink_lu_inv_thr 1e-06 #threshold to get inverse of overlap matrix by LU decomposition in auxiliary basis representation + +#Parameters (16.dft+u) +dft_plus_u 0 #DFT+U correction method +dft_plus_dmft 0 #true:DFT+DMFT; false: standard DFT calcullation(default) +orbital_corr -1 #which correlated orbitals need corrected ; d:2 ,f:3, do not need correction:-1 +hubbard_u 0 #Hubbard Coulomb interaction parameter U(ev) +yukawa_potential 0 #default: false +yukawa_lambda -1 #default:0.0 +uramping -1 #increasing U values during SCF +omc 0 #the mode of occupation matrix control +onsite_radius 0 #radius of the sphere for onsite projection (Bohr) + +#Parameters (17.non-collinear spin-constrained DFT) +sc_mag_switch 0 #switch to control spin-constrained DFT +decay_grad_switch 0 #switch to control gradient break condition +sc_thr 1e-06 #Convergence criterion of spin-constrained iteration (RMS) in uB +nsc 100 #Maximal number of spin-constrained iteration +nsc_min 2 #Minimum number of spin-constrained iteration +alpha_trial 0.01 #Initial trial step size for lambda in eV/uB^2 +sccut 3 #Maximal step size for lambda in eV/uB +sc_drop_thr 0.001 #Convergence criterion ratio of lambda iteration in Spin-constrained DFT +sc_scf_thr 0.001 #Density error threshold for inner loop of spin-constrained SCF +sc_direction_only 0 #only optimize the direction of magnetization +sc_lambda_strategy bfgs #lambda update strategy for spin-constrained DFT +sc_scan_lambda_start 0 #start value for linear lambda scan (eV/uB) +sc_scan_lambda_end 1 #end value for linear lambda scan (eV/uB) +sc_scan_steps 20 #number of steps for linear lambda scan + +#Parameters (18.Quasiatomic Orbital analysis) +qo_switch 0 #switch to control quasiatomic orbital analysis +qo_basis szv #type of QO basis function: hydrogen: hydrogen-like basis, pswfc: read basis from pseudopotential +qo_strategy all #strategy to generate generate radial orbitals +qo_screening_coeff 0.1 #rescale the shape of radial orbitals +qo_thr 1e-06 #accuracy for evaluating cutoff radius of QO basis function + +#Parameters (19.PEXSI) +pexsi_npole 40 #Number of poles in expansion +pexsi_inertia 1 #Whether inertia counting is used at the very beginning of PEXSI process +pexsi_nmax 80 #Maximum number of PEXSI iterations after each inertia counting procedure +pexsi_comm 1 #Whether to construct PSelInv communication pattern +pexsi_storage 1 #Storage space used by the Selected Inversion algorithm for symmetric matrices +pexsi_ordering 0 #Ordering strategy for factorization and selected inversion +pexsi_row_ordering 1 #Row permutation strategy for factorization and selected inversion, 0: NoRowPerm, 1: LargeDiag +pexsi_nproc 1 #Number of processors for parmetis +pexsi_symm 1 #Matrix symmetry +pexsi_trans 0 #Whether to transpose +pexsi_method 1 #pole expansion method, 1: Cauchy Contour Integral, 2: Moussa optimized method +pexsi_nproc_pole 1 #Number of processes used by each pole +pexsi_temp 0.015 #Temperature, in the same unit as H +pexsi_gap 0 #Spectral gap +pexsi_delta_e 20 #An upper bound for the spectral radius of S^{-1} H +pexsi_mu_lower -10 #Initial guess of lower bound for mu +pexsi_mu_upper 10 #Initial guess of upper bound for mu +pexsi_mu 0 #Initial guess for mu (for the solver) +pexsi_mu_thr 0.05 #Stopping criterion in terms of the chemical potential for the inertia counting procedure +pexsi_mu_expand 0.3 #If the chemical potential is not in the initial interval, the interval is expanded by muInertiaExpansion +pexsi_mu_guard 0.2 #Safe guard criterion in terms of the chemical potential to reinvoke the inertia counting procedure +pexsi_elec_thr 0.001 #Stopping criterion of the PEXSI iteration in terms of the number of electrons compared to numElectronExact +pexsi_zero_thr 1e-10 #if the absolute value of matrix element is less than ZERO_Limit, it will be considered as 0 + +#Parameters (20.Test) +out_alllog 0 #output information for each processor, when parallel +nurse 0 #for coders +t_in_h 1 #calculate the kinetic energy or not +vl_in_h 1 #calculate the local potential or not +vnl_in_h 1 #calculate the nonlocal potential or not +vh_in_h 1 #calculate the hartree potential or not +vion_in_h 1 #calculate the local ionic potential or not +test_force 0 #test the force +test_stress 0 #test the stress +test_skip_ewald 0 #whether to skip ewald +ri_hartree_benchmark none #whether to use the RI approximation for the Hartree term in LR-TDDFT for benchmark (with FHI-aims/ABACUS read-in style) +aims_nbasis #the number of basis functions for each atom type used in FHI-aims (for benchmark) +rdmft 0 #whether to perform rdmft calculation, default is false +rdmft_power_alpha 0.656 #the alpha parameter of power-functional, g(occ_number) = occ_number^alpha used in exx-type functionals such as muller and power +exxace 0 #whether to perform ace calculation in exxpw +exx_gamma_extra 0 #whether to perform gamma extrapolation in exxpw +ecutexx 0 #energy cutoff for exx calculation, Ry +exx_batch_size 0 #band chunk width of the EXX batched path, 0 = all bands +exx_thr_type density #threshold type for exx outer loop, energy or density +exx_ene_thr 1e-05 #threshold for exx outer loop when exx_thr_type = energy diff --git a/tests/01_PW/nscf_hse_same/KPT b/tests/01_PW/nscf_hse_same/KPT new file mode 100644 index 0000000000..e769af7638 --- /dev/null +++ b/tests/01_PW/nscf_hse_same/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 1 1 0 0 0 diff --git a/tests/01_PW/nscf_hse_same/README b/tests/01_PW/nscf_hse_same/README new file mode 100644 index 0000000000..d93148add8 --- /dev/null +++ b/tests/01_PW/nscf_hse_same/README @@ -0,0 +1,3 @@ +PW HSE NSCF on the frozen SCF source committed in this directory (binary +wavefunctions, eig_occ.txt, charge restart). Same k mesh with more target +bands; an explicit zero exx_fock_alpha is accepted like an omitted Fock term. diff --git a/tests/01_PW/nscf_hse_same/STRU b/tests/01_PW/nscf_hse_same/STRU new file mode 100644 index 0000000000..c35003c65b --- /dev/null +++ b/tests/01_PW/nscf_hse_same/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +H 1 H_ONCV_PBE-1.0.upf + +LATTICE_CONSTANT +1.889766 + +LATTICE_VECTORS +0.0 2.708337 2.708337 +2.708337 0.0 2.708337 +2.708337 2.708337 0.0 + +ATOMIC_POSITIONS +Direct + +H +0 +2 +0.125 0.125 0.125 0 0 0 +0.675 0.835 0.875 0 0 0 diff --git a/tests/01_PW/nscf_hse_same/autotest-CHARGE-DENSITY.restart b/tests/01_PW/nscf_hse_same/autotest-CHARGE-DENSITY.restart new file mode 100644 index 0000000000..f7c30f6865 Binary files /dev/null and b/tests/01_PW/nscf_hse_same/autotest-CHARGE-DENSITY.restart differ diff --git a/tests/01_PW/nscf_hse_same/band.txt.ref b/tests/01_PW/nscf_hse_same/band.txt.ref new file mode 100644 index 0000000000..6648ef56ac --- /dev/null +++ b/tests/01_PW/nscf_hse_same/band.txt.ref @@ -0,0 +1,2 @@ + 1 0.00000000 -11.12767527 0.23237114 5.91382475 9.84140426 10.33312736 + 2 0.15988140 -8.18994369 -4.74422564 4.63092625 7.73162354 8.50832814 diff --git a/tests/01_PW/nscf_hse_same/eig_occ.txt b/tests/01_PW/nscf_hse_same/eig_occ.txt new file mode 100644 index 0000000000..cf5895f7dc --- /dev/null +++ b/tests/01_PW/nscf_hse_same/eig_occ.txt @@ -0,0 +1,13 @@ +1 # ionic step + Electronic state energy (eV) and occupations + Spin number 1 + spin=1 k-point=1/2 Cartesian=0.0000000 0.0000000 0.0000000 (137 plane wave) + 1 -11.12768112632454 1.000000000000000 + 2 0.2323374417247880 0.000000000000000 + 3 5.913816152912471 0.000000000000000 + + spin=1 k-point=2/2 Cartesian=-0.09230756733744730 0.09230756733744730 0.09230756733744730 (138 plane wave) + 1 -8.189957518901499 1.000000000000000 + 2 -4.744242172101769 0.000000000000000 + 3 4.630928303470863 0.000000000000000 + diff --git a/tests/01_PW/nscf_hse_same/result.ref b/tests/01_PW/nscf_hse_same/result.ref new file mode 100644 index 0000000000..d98de69975 --- /dev/null +++ b/tests/01_PW/nscf_hse_same/result.ref @@ -0,0 +1,4 @@ +etotref -33.63676769132638 +etotperatomref -16.8183838457 +CompareBand_pass 0 +totaltimeref 0.08 diff --git a/tests/01_PW/nscf_hse_same/wfk1_pw.dat b/tests/01_PW/nscf_hse_same/wfk1_pw.dat new file mode 100644 index 0000000000..da05022468 Binary files /dev/null and b/tests/01_PW/nscf_hse_same/wfk1_pw.dat differ diff --git a/tests/01_PW/nscf_hse_same/wfk2_pw.dat b/tests/01_PW/nscf_hse_same/wfk2_pw.dat new file mode 100644 index 0000000000..64e9328b42 Binary files /dev/null and b/tests/01_PW/nscf_hse_same/wfk2_pw.dat differ