diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 6200ccc834..665c9f9da6 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -732,7 +732,7 @@ If only one value is set (such as `kspacing 0.5`), then kspacing values of a/b/c [back to top](#full-list-of-input-keywords) -## Variables related to input files +## Input files These variables are used to control parameters related to input files. @@ -1656,16 +1656,16 @@ These variables are used to control the geometry relaxation. [back to top](#full-list-of-input-keywords) -## Variables related to output information +## Output information These variables are used to control the output of properties. ### out_freq_ion - **Type**: Integer -- **Description**: After self-consistent-field calculations, control the interval of ionic movements for printing properties. These properties cover charge density, local potential, electrostatic potential, Hamiltonian matrix, overlap matrix, density matrix, Mulliken population analysis and so on. +- **Description**: Control the interval to print information every few ion steps. These properties cover charge density, local potential, electrostatic potential, Hamiltonian matrix, overlap matrix, density matrix, Mulliken population analysis and so on. - **Default**: 0 -- **Note**: If you want to use out_freq_elec, please set out_freq_ion to 1, otherwise out_freq_elec is useless +- **Note**: The integer indicates to print information every 'out_freq_ion' ion steps. ### out_freq_elec @@ -1679,11 +1679,11 @@ These variables are used to control the output of properties. - **Description**: The first integer controls whether to output the charge density on real space grids: - 1: Output the charge density (in Bohr^-3) on real space grids into the density files in the folder `OUT.${suffix}`. The files are named as: - - nspin = 1: `chgs1.cube`; + - nspin = 1: `chg.cube`; - nspin = 2: `chgs1.cube`, and `chgs2.cube`; - nspin = 4: `chgs1.cube`, `chgs2.cube`, `chgs3.cube`, and `chgs4.cube`; Note that by using the Meta-GGA functional, additional files containing the kinetic energy density will be output with the following names: - - nspin = 1: `taus1.cube`; + - nspin = 1: `tau.cube`; - nspin = 2: `taus1.cube`, and `taus2.cube`; - nspin = 4: `taus1.cube`, `taus2.cube`, `taus3.cube`, and `taus4.cube`; - 2: On top of 1, also output the initial charge density files with a suffix name as '_ini', such as `taus1_ini.cube`, etc. diff --git a/source/Makefile.Objects b/source/Makefile.Objects index ed356b4329..db5585dd3a 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -579,6 +579,7 @@ OBJS_IO=input_conv.o\ output_log.o\ output_mat_sparse.o\ ctrl_output_lcao.o\ + ctrl_output_fp.o\ para_json.o\ abacusjson.o\ general_info.o\ diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 46f00cf284..9d8326d2a3 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -8,20 +8,15 @@ #include "source_hamilt/module_vdw/vdw.h" #include "source_pw/module_pwdft/global.h" #include "source_io/cif_io.h" -#include "source_io/cube_io.h" +#include "source_io/cube_io.h" // use write_vdata_palgrid #include "source_io/json_output/init_info.h" #include "source_io/json_output/output_info.h" #include "source_io/output_log.h" #include "source_io/print_info.h" #include "source_io/rhog_io.h" -#include "source_io/write_elecstat_pot.h" -#include "source_io/write_elf.h" #include "source_io/module_parameter/parameter.h" #include "source_cell/k_vector_utils.h" - -#ifdef USE_LIBXC -#include "source_io/write_libxc_r.h" -#endif +#include "source_io/ctrl_output_fp.h" namespace ModuleESolver { @@ -162,169 +157,10 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso // 3) update delta_rho for charge extrapolation CE.update_delta_rho(ucell, &(this->chr), &(this->sf)); + // 4) print out charge density, potential, elf, etc. + ModuleIO::ctrl_output_fp(ucell, this->pelec, this->pw_big, this->pw_rhod, + this->chr, this->solvent, this->Pgrid, istep); - // print out the 'g' index when istep_in!=-1 - int istep_in = -1; - if (PARAM.inp.out_freq_ion>0) // default value of out_freq_ion is 0 - { - if (istep % PARAM.inp.out_freq_ion == 0) - { - istep_in = istep; - } - } - - std::string geom_block; - if(istep_in==-1) - { - // do nothing - } - else if(istep_in>=0) - { - geom_block = "g" + std::to_string(istep + 1); - } - - // 4) write charge density - if (PARAM.inp.out_chg[0] > 0) - { - for (int is = 0; is < PARAM.inp.nspin; ++is) - { - this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]); - - std::string fn =PARAM.globalv.global_out_dir + "chg"; - - std::string spin_block; - if(PARAM.inp.nspin == 2 || PARAM.inp.nspin == 4) - { - spin_block= "s" + std::to_string(is + 1); - } - else if(PARAM.inp.nspin == 1) - { - // do nothing - } - - fn += spin_block + geom_block + ".cube"; - - ModuleIO::write_vdata_palgrid(Pgrid, - this->chr.rho_save[is], - is, - PARAM.inp.nspin, - istep_in, - fn, - this->pelec->eferm.get_efval(is), - &(ucell), - PARAM.inp.out_chg[1], - 1); - - if (XC_Functional::get_ked_flag()) - { - fn = PARAM.globalv.global_out_dir + "tau"; - - fn += spin_block + geom_block + ".cube"; - - ModuleIO::write_vdata_palgrid(Pgrid, - this->chr.kin_r_save[is], - is, - PARAM.inp.nspin, - istep, - fn, - this->pelec->eferm.get_efval(is), - &(ucell)); - } - } - } - - // 5) write potential - if (PARAM.inp.out_pot == 1 || PARAM.inp.out_pot == 3) - { - for (int is = 0; is < PARAM.inp.nspin; is++) - { - std::string fn =PARAM.globalv.global_out_dir + "pot"; - - std::string spin_block; - if(PARAM.inp.nspin == 2 || PARAM.inp.nspin == 4) - { - spin_block= "s" + std::to_string(is + 1); - } - else if(PARAM.inp.nspin == 1) - { - // do nothing - } - - fn += spin_block + geom_block + ".cube"; - - ModuleIO::write_vdata_palgrid(Pgrid, - this->pelec->pot->get_effective_v(is), - is, - PARAM.inp.nspin, - istep_in, - fn, - 0.0, // efermi - &(ucell), - 3, // precision - 0); // out_fermi - } - } - else if (PARAM.inp.out_pot == 2) - { - std::string fn =PARAM.globalv.global_out_dir + "potes"; - fn += geom_block + ".cube"; - - ModuleIO::write_elecstat_pot( -#ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, -#endif - fn, - istep, - this->pw_rhod, - &this->chr, - &(ucell), - this->pelec->pot->get_fixed_v(), - this->solvent); - } - - // 6) write ELF - if (PARAM.inp.out_elf[0] > 0) - { - this->chr.cal_elf = true; - Symmetry_rho srho; - for (int is = 0; is < PARAM.inp.nspin; is++) - { - srho.begin(is, this->chr, this->pw_rhod, ucell.symm); - } - - std::string out_dir =PARAM.globalv.global_out_dir; - ModuleIO::write_elf( -#ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, -#endif - out_dir, - istep, - PARAM.inp.nspin, - this->chr.rho, - this->chr.kin_r, - this->pw_rhod, - this->Pgrid, - &(ucell), - PARAM.inp.out_elf[1]); - } - -#ifdef USE_LIBXC - // 7) write xc(r) - if(PARAM.inp.out_xc_r[0]>=0) - { - ModuleIO::write_libxc_r( - PARAM.inp.out_xc_r[0], - XC_Functional::get_func_id(), - this->pw_rhod->nrxx, // number of real-space grid - ucell.omega, // volume of cell - ucell.tpiba, - this->chr, - *this->pw_big, - *this->pw_rhod); - } -#endif } void ESolver_FP::before_scf(UnitCell& ucell, const int istep) @@ -402,19 +238,32 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) //---------------------------------------------------------- elecstate::cal_ux(ucell); + + //---------------------------------------------------------- //! output the initial charge density //---------------------------------------------------------- + const int nspin = PARAM.inp.nspin; if (PARAM.inp.out_chg[0] == 2) { - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < nspin; is++) { std::stringstream ss; - ss << PARAM.globalv.global_out_dir << "chgs" << is + 1 << "_ini.cube"; + ss << PARAM.globalv.global_out_dir << "chg"; + + if(nspin==1) + { + ss << "ini.cube"; + } + else if(nspin==2 || nspin==4) + { + ss << "s" << is + 1 << "ini.cube"; + } + ModuleIO::write_vdata_palgrid(this->Pgrid, this->chr.rho[is], is, - PARAM.inp.nspin, + nspin, istep, ss.str(), this->pelec->eferm.ef, @@ -427,14 +276,24 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) //---------------------------------------------------------- if (PARAM.inp.out_pot == 3) { - for (int is = 0; is < PARAM.inp.nspin; is++) + for (int is = 0; is < nspin; is++) { std::stringstream ss; - ss << PARAM.globalv.global_out_dir << "pots" << is + 1 << "_ini.cube"; + ss << PARAM.globalv.global_out_dir << "pot"; + + if(nspin==1) + { + ss << "ini.cube"; + } + else if(nspin==2 || nspin==4) + { + ss << "s" << is + 1 << "ini.cube"; + } + ModuleIO::write_vdata_palgrid(this->Pgrid, this->pelec->pot->get_effective_v(is), is, - PARAM.inp.nspin, + nspin, istep, ss.str(), 0.0, // efermi diff --git a/source/source_esolver/esolver_fp.h b/source/source_esolver/esolver_fp.h index fda47313bc..7d4f34736e 100644 --- a/source/source_esolver/esolver_fp.h +++ b/source/source_esolver/esolver_fp.h @@ -7,26 +7,13 @@ #include #endif -//! plane wave basis -#include "source_basis/module_pw/pw_basis.h" - -//! symmetry analysis -#include "source_cell/module_symmetry/symmetry.h" - -//! electronic states -#include "source_estate/elecstate.h" - -//! charge extrapolation -#include "source_estate/module_charge/charge_extra.h" - -//! solvation model -#include "source_hamilt/module_surchem/surchem.h" - -//! local pseudopotential -#include "source_pw/module_pwdft/VL_in_pw.h" - -//! structure factor related to plane wave basis -#include "source_pw/module_pwdft/structure_factor.h" +#include "source_basis/module_pw/pw_basis.h" // plane wave basis +#include "source_cell/module_symmetry/symmetry.h" // symmetry analysis +#include "source_estate/elecstate.h" // electronic states +#include "source_estate/module_charge/charge_extra.h" // charge extrapolation +#include "source_hamilt/module_surchem/surchem.h" // solvation model +#include "source_pw/module_pwdft/VL_in_pw.h" // local pseudopotential +#include "source_pw/module_pwdft/structure_factor.h" // structure factor #include @@ -44,10 +31,8 @@ namespace ModuleESolver class ESolver_FP: public ESolver { public: - //! Constructor ESolver_FP(); - //! Deconstructor virtual ~ESolver_FP(); //! Initialize of the first-principels energy solver @@ -56,13 +41,10 @@ class ESolver_FP: public ESolver virtual void after_all_runners(UnitCell& ucell) override; protected: - //! Something to do before SCF iterations. virtual void before_scf(UnitCell& ucell, const int istep); - //! Something to do after SCF iterations when SCF is converged or comes to the max iter step. virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver); - //! Something to do after hamilt2rho function in each iter loop. virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver); //! ------------------------------------------------------------------------------ diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index 1ece91eb62..45f26babcb 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -11,9 +11,7 @@ #include "source_pw/module_pwdft/global.h" #include "source_io/print_info.h" #include "source_estate/cal_ux.h" -//-----force------------------- #include "source_pw/module_pwdft/forces.h" -//-----stress------------------ #include "source_pw/module_ofdft/of_stress_pw.h" namespace ModuleESolver diff --git a/source/source_estate/elecstate.h b/source/source_estate/elecstate.h index bdf5f9cc77..7b1ff8fbcc 100644 --- a/source/source_estate/elecstate.h +++ b/source/source_estate/elecstate.h @@ -51,8 +51,6 @@ class ElecState { return; } - // virtual void updateRhoK(const psi::Psi> &psi) = 0; - // virtual void updateRhoK(const psi::Psi &psi)=0 virtual void cal_tau(const psi::Psi>& psi) { return; @@ -77,12 +75,9 @@ class ElecState return; } - - // use occupied weights from INPUT and skip calculate_weights // mohan updated on 2024-06-08 - // if nupdown is not 0(TWO_EFERMI case), // nelec_spin will be fixed and weights will be constrained void init_nelec_spin(); diff --git a/source/source_io/CMakeLists.txt b/source/source_io/CMakeLists.txt index 6c7b5bf349..ffa7598cd9 100644 --- a/source/source_io/CMakeLists.txt +++ b/source/source_io/CMakeLists.txt @@ -1,5 +1,6 @@ list(APPEND objects input_conv.cpp + ctrl_output_fp.cpp bessel_basis.cpp cal_test.cpp cal_dos.cpp diff --git a/source/source_io/ctrl_output_fp.cpp b/source/source_io/ctrl_output_fp.cpp index 477612ed27..0cde6e6011 100644 --- a/source/source_io/ctrl_output_fp.cpp +++ b/source/source_io/ctrl_output_fp.cpp @@ -1,14 +1,26 @@ #include "source_io/ctrl_output_fp.h" // use ctrl_output_fp() +#include "source_estate/module_charge/symmetry_rho.h" // use Symmetry_rho +#include "source_io/write_elecstat_pot.h" // use write_elecstat_pot +#include "source_io/write_elf.h" +#include "cube_io.h" // use write_vdata_palgrid +#include "source_hamilt/module_xc/xc_functional.h" // use XC_Functional + +#ifdef USE_LIBXC +#include "source_io/write_libxc_r.h" +#endif namespace ModuleIO { -template void ctrl_output_fp(UnitCell& ucell, - elecstate::ElecStateLCAO* pelec, + elecstate::ElecState* pelec, + ModulePW::PW_Basis_Big* pw_big, + ModulePW::PW_Basis* pw_rhod, + Charge &chr, + surchem &solvent, + Parallel_Grid ¶_grid, const int istep) { -/* ModuleBase::TITLE("ModuleIO", "ctrl_output_fp"); ModuleBase::timer::tick("ModuleIO", "ctrl_output_fp"); @@ -17,121 +29,172 @@ void ctrl_output_fp(UnitCell& ucell, const int nspin = PARAM.inp.nspin; const std::string global_out_dir = PARAM.globalv.global_out_dir; - // 1) write charge density - if (PARAM.inp.out_chg[0] > 0) - { - for (int is = 0; is < PARAM.inp.nspin; is++) - { - this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]); - std::string fn =PARAM.globalv.global_out_dir + "/chgs" + std::to_string(is + 1) + ".cube"; - ModuleIO::write_vdata_palgrid(Pgrid, - this->chr.rho_save[is], - is, - PARAM.inp.nspin, - istep, - fn, - this->pelec->eferm.get_efval(is), - &(ucell), - PARAM.inp.out_chg[1], - 1); - - if (XC_Functional::get_ked_flag()) - { - fn =PARAM.globalv.global_out_dir + "/taus" + std::to_string(is + 1) + ".cube"; - ModuleIO::write_vdata_palgrid(Pgrid, - this->chr.kin_r_save[is], - is, - PARAM.inp.nspin, - istep, - fn, - this->pelec->eferm.get_efval(is), - &(ucell)); - } - } - } - - - // 2) write potential - if (PARAM.inp.out_pot == 1 || PARAM.inp.out_pot == 3) - { - for (int is = 0; is < PARAM.inp.nspin; is++) - { - std::string fn =PARAM.globalv.global_out_dir + "/pots" + std::to_string(is + 1) + ".cube"; - - ModuleIO::write_vdata_palgrid(Pgrid, - this->pelec->pot->get_effective_v(is), - is, - PARAM.inp.nspin, - istep, - fn, - 0.0, // efermi - &(ucell), - 3, // precision - 0); // out_fermi - } - } - else if (PARAM.inp.out_pot == 2) - { - std::string fn =PARAM.globalv.global_out_dir + "/pot_es.cube"; - ModuleIO::write_elecstat_pot( + + // print out the 'g' index when istep_in!=-1 + int istep_in = -1; + if (PARAM.inp.out_freq_ion>0) // default value of out_freq_ion is 0 + { + if (istep % PARAM.inp.out_freq_ion == 0) + { + istep_in = istep; + } + } + + std::string geom_block; + if(istep_in==-1) + { + // do nothing + } + else if(istep_in>=0) + { + geom_block = "g" + std::to_string(istep + 1); + } + + + // 4) write charge density + if (PARAM.inp.out_chg[0] > 0) + { + for (int is = 0; is < nspin; ++is) + { + pw_rhod->real2recip(chr.rho_save[is], chr.rhog_save[is]); + + std::string fn =PARAM.globalv.global_out_dir + "chg"; + + std::string spin_block; + if(nspin == 2 || nspin == 4) + { + spin_block= "s" + std::to_string(is + 1); + } + else if(nspin == 1) + { + // do nothing + } + + fn += spin_block + geom_block + ".cube"; + + ModuleIO::write_vdata_palgrid(para_grid, + chr.rho_save[is], + is, + nspin, + istep_in, + fn, + pelec->eferm.get_efval(is), + &(ucell), + PARAM.inp.out_chg[1], + 1); + + if (XC_Functional::get_ked_flag()) + { + fn = PARAM.globalv.global_out_dir + "tau"; + + fn += spin_block + geom_block + ".cube"; + + ModuleIO::write_vdata_palgrid(para_grid, + chr.kin_r_save[is], + is, + nspin, + istep, + fn, + pelec->eferm.get_efval(is), + &(ucell)); + } + } + } + + // 5) write potential + if (PARAM.inp.out_pot == 1 || PARAM.inp.out_pot == 3) + { + for (int is = 0; is < nspin; is++) + { + std::string fn =PARAM.globalv.global_out_dir + "pot"; + + std::string spin_block; + if(nspin == 2 || nspin == 4) + { + spin_block= "s" + std::to_string(is + 1); + } + else if(nspin == 1) + { + // do nothing + } + + fn += spin_block + geom_block + ".cube"; + + ModuleIO::write_vdata_palgrid(para_grid, + pelec->pot->get_effective_v(is), + is, + nspin, + istep_in, + fn, + 0.0, // efermi + &(ucell), + 3, // precision + 0); // out_fermi + } + } + else if (PARAM.inp.out_pot == 2) + { + std::string fn =PARAM.globalv.global_out_dir + "potes"; + fn += geom_block + ".cube"; + + ModuleIO::write_elecstat_pot( #ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, + pw_big->bz, + pw_big->nbz, #endif - fn, - istep, - this->pw_rhod, - &this->chr, - &(ucell), - this->pelec->pot->get_fixed_v(), - this->solvent); - } - - - // 3) write ELF - if (PARAM.inp.out_elf[0] > 0) - { - this->chr.cal_elf = true; - Symmetry_rho srho; - for (int is = 0; is < PARAM.inp.nspin; is++) - { - srho.begin(is, this->chr, this->pw_rhod, ucell.symm); - } - - std::string out_dir =PARAM.globalv.global_out_dir; - ModuleIO::write_elf( + fn, + istep, + pw_rhod, + &chr, + &(ucell), + pelec->pot->get_fixed_v(), + solvent); + } + + // 6) write ELF + if (PARAM.inp.out_elf[0] > 0) + { + chr.cal_elf = true; + Symmetry_rho srho; + for (int is = 0; is < nspin; is++) + { + srho.begin(is, chr, pw_rhod, ucell.symm); + } + + std::string out_dir =PARAM.globalv.global_out_dir; + ModuleIO::write_elf( #ifdef __MPI - this->pw_big->bz, - this->pw_big->nbz, + pw_big->bz, + pw_big->nbz, +#endif + out_dir, + istep, + nspin, + chr.rho, + chr.kin_r, + pw_rhod, + para_grid, + &(ucell), + PARAM.inp.out_elf[1]); + } + +#ifdef USE_LIBXC + // 7) write xc(r) + if(PARAM.inp.out_xc_r[0]>=0) + { + ModuleIO::write_libxc_r( + PARAM.inp.out_xc_r[0], + XC_Functional::get_func_id(), + pw_rhod->nrxx, // number of real-space grid + ucell.omega, // volume of cell + ucell.tpiba, + chr, + *pw_big, + *pw_rhod); + } #endif - out_dir, - istep, - PARAM.inp.nspin, - this->chr.rho, - this->chr.kin_r, - this->pw_rhod, - this->Pgrid, - &(ucell), - PARAM.inp.out_elf[1]); - } ModuleBase::timer::tick("ModuleIO", "ctrl_output_fp"); - -*/ - } } // End ModuleIO - - -// For gamma only -template void ModuleIO::ctrl_output_lcao(UnitCell& ucell, - const int istep); - -// For multiple k-points -template void ModuleIO::ctrl_output_lcao, double>(UnitCell& ucell, - const int istep); - -template void ModuleIO::ctrl_output_lcao, std::complex>(UnitCell& ucell, - const int istep); - diff --git a/source/source_io/ctrl_output_fp.h b/source/source_io/ctrl_output_fp.h index 98a4ad1f6d..12ff1e12b0 100644 --- a/source/source_io/ctrl_output_fp.h +++ b/source/source_io/ctrl_output_fp.h @@ -1,11 +1,19 @@ #ifndef CTRL_OUTPUT_FP_H #define CTRL_OUTPUT_FP_H +#include "source_estate/elecstate_lcao.h" + namespace ModuleIO { - template - void ctrl_output_fp(UnitCell& ucell, - elecstate::ElecStateLCAO* pelec, - const int istep); + + void ctrl_output_fp(UnitCell& ucell, + elecstate::ElecState* pelec, + ModulePW::PW_Basis_Big* pw_big, + ModulePW::PW_Basis* pw_rhod, + Charge &chr, + surchem &solvent, + Parallel_Grid ¶_grid, + const int istep); + } #endif