read_inp_out.cpp

#include "source_base/formatter.h"
#include "source_base/global_function.h"
#include "source_base/tool_quit.h"
#include "read_input.h"
#include "read_input_tool.h"
namespace ModuleIO
{
void ReadInput::item_output()
{
    // NOTE: The order of add_item() calls below determines the parameter order
    // in the generated documentation (docs/advanced/input_files/input-main.md).
    // Please preserve this ordering when adding new parameters.
    {
        Input_Item item("out_freq_ion");
        item.annotation = "print information every few ionic steps";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Controls the output interval in ionic steps. When set to a positive integer, information such as charge density, local potential, electrostatic potential, Hamiltonian matrix, overlap matrix, density matrix, Mulliken population analysis, and structure files (STRU{istep} or STRU{istep}.cif, when out_stru is 1 or 2) is printed every n ionic steps."
                          "\n\n[NOTE] In RT-TDDFT calculations, this parameter is inactive; output frequency is instead controlled by out_freq_td.";
        item.default_value = "0";
        item.unit = "";
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.out_freq_ion <= 0)
            {
                para.input.out_freq_ion = 0; // 0 means no output of info
            }
        };
        read_sync_int(input.out_freq_ion);
        this->add_item(item);
    }
    {
        Input_Item item("out_freq_td");
        item.annotation = "print information every few completed electronic iterations in RT-TDDFT";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Controls the output interval in completed electronic evolution steps during RT-TDDFT calculations. When set to a positive integer n, detailed information (see out_freq_ion) is printed every n electron time-evolution steps (i.e., every STEP OF ELECTRON EVOLVE). For example, if you wish to output information once per ionic step, you should set out_freq_td equal to estep_per_md, since one ionic step corresponds to estep_per_md electronic evolution steps."
                          "\n\n[NOTE] This parameter is only active in RT-TDDFT mode (esolver_type = tddft). It has no effect in ground-state calculations.";
        item.default_value = "0";
        item.unit = "";
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.out_freq_td <= 0)
            {
                para.input.out_freq_td = 0; // 0 means no output of info
            }
        };
        read_sync_int(input.out_freq_td);
        this->add_item(item);
    }
    {
        Input_Item item("out_freq_elec");
        item.annotation = "print information every few electronic steps";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Output the charge density (only binary format, controlled by out_chg), wavefunction (controlled by out_wfc_pw) per out_freq_elec electronic iterations. Note that they are always output when converged or reach the maximum iterations scf_nmax.";
        item.default_value = "scf_nmax";
        item.unit = "";
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.out_freq_elec <= 0)
            {
                para.input.out_freq_elec = para.input.scf_nmax;
            }
        };
        read_sync_int(input.out_freq_elec);
        this->add_item(item);
    }
    {
        Input_Item item("out_chg");
        item.annotation = "> 0 output charge density for selected electron steps"
                          ", second parameter controls the precision, default is 3.";
        item.category = "Output information";
        item.type = R"(Integer \[Integer\](optional))";
                item.description = R"(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: `chg.cube`;
        - nspin = 2: `chgs1.cube`, and `chgs2.cube`;
        - nspin = 4: `chgs1.cube`, `chgs2.cube`, `chgs3.cube`, and `chgs4.cube`;
        - When using the Meta-GGA functional, additional files containing the kinetic energy density are also output:
            - 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. The files are named as:
        - out_freq_ion = 0:
            - nspin = 1: `chg_ini.cube`;
            - nspin = 2: `chgs1_ini.cube` and `chgs2_ini.cube`;
            - nspin = 4: `chgs1_ini.cube`, `chgs2_ini.cube`, `chgs3_ini.cube`, and `chgs4_ini.cube`;
            - output at every step (overwrite same file)
        - out_freq_ion > 0:
            - nspin = 1: `chgg{geom_step}_ini.cube` (e.g., `chgg1_ini.cube`);
            - nspin = 2: `chgs1g{geom_step}_ini.cube` and `chgs2g{geom_step}_ini.cube`;
            - nspin = 4: `chgs1g{geom_step}_ini.cube`, `chgs2g{geom_step}_ini.cube`, `chgs3g{geom_step}_ini.cube`, and `chgs4g{geom_step}_ini.cube`.
            - output every out_freq_ion steps
        Here, {geom_step} denotes the geometry step index, starting from 1 (geom_step = istep + 1).
    - -1: Disable the charge density auto-back-up file `{suffix}-CHARGE-DENSITY.restart`, useful for large systems.

The second integer controls the precision of the charge density output. If not given, `3` is used as default. For restarting from this file and other high-precision calculations, `10` is recommended.

In molecular dynamics simulations, the output frequency is controlled by out_freq_ion.

[NOTE] In the 3.10-LTS version, the file names are SPIN1_CHG.cube and SPIN1_CHG_INI.cube, etc.)";
        item.default_value = "0 3";
        item.unit = "";
			item.read_value = [](const Input_Item& item, Parameter& para) {
				const size_t count = item.get_size();
				if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_chg needs at least 1 value");
				para.input.out_chg[0] = std::stoi(item.str_values[0]);
            para.input.out_chg[1] = 3;
			if (count >= 2) try { para.input.out_chg[1] = std::stoi(item.str_values[1]); }
			catch (const std::invalid_argument&) { /* do nothing */ }
			catch (const std::out_of_range&) {/* do nothing */}
		};
        // reset value in some special case
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            para.input.out_chg[0] = (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
                                        ? 1
                                        : para.input.out_chg[0];
        };
        sync_intvec(input.out_chg, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_pot");
        item.annotation = "output real space potential (with precision 8)";
        item.category = "Output information";
        item.type = R"(Integer \[Integer\](optional))";
        item.description = R"(* 1: Output the total local potential (i.e., local pseudopotential + Hartree potential + XC potential + external electric field (if exists) + dipole correction potential (if exists) + ...) on real space grids (in Ry) into files in the folder OUT.{suffix}. The files are named as:
 * nspin = 1: pots1.cube;
 * nspin = 2: pots1.cube and pots2.cube;
 * nspin = 4: pots1.cube, pots2.cube, pots3.cube, and pots4.cube
* 2: Output the electrostatic potential on real space grids into OUT.{suffix}/pot_es.cube. The Python script named tools/02_postprocessing/average_pot/aveElecStatPot.py can be used to calculate the average electrostatic potential along the z-axis and outputs it into ElecStaticPot_AVE. Please note that the total local potential refers to the local component of the self-consistent potential, excluding the non-local pseudopotential. The distinction between the local potential and the electrostatic potential is as follows: local potential = electrostatic potential + XC potential.
* 3: Apart from 1, also output the total local potential of the initial charge density. The files are named as:
 * out_freq_ion = 0:
   * nspin = 1: `pot_ini.cube`;
   * nspin = 2: `pots1_ini.cube` and `pots2_ini.cube`;
   * nspin = 4: `pots1_ini.cube`, `pots2_ini.cube`, `pots3_ini.cube`, and `pots4_ini.cube`;
   * output at every step (overwrite same file)
 * out_freq_ion > 0:
   * nspin = 1: `potg{geom_step}_ini.cube` (e.g., `potg1_ini.cube`);
   * nspin = 2: `pots1g{geom_step}_ini.cube` and `pots2g{geom_step}_ini.cube`;
   * nspin = 4: `pots1g{geom_step}_ini.cube`, `pots2g{geom_step}_ini.cube`, `pots3g{geom_step}_ini.cube`, and `pots4g{geom_step}_ini.cube`.
   * output every out_freq_ion steps
 Here, {geom_step} denotes the geometry step index, starting from 1 (geom_step = istep + 1).

The optional second integer controls the output precision. If not provided, the default precision is 8.

In molecular dynamics calculations, the output frequency is controlled by out_freq_ion.

[NOTE] In the 3.10-LTS version, the file names are SPIN1_POT.cube and SPIN1_POT_INI.cube, etc.)";
        item.default_value = "0";
        item.unit = "";
        item.read_value = [](const Input_Item& item, Parameter& para) {
                    const size_t count = item.get_size();
                    if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_pot needs at least 1 value");
                    para.input.out_pot[0] = std::stoi(item.str_values[0]);
                    para.input.out_pot[1] = 8;
                    if (count >= 2) try { para.input.out_pot[1] = std::stoi(item.str_values[1]); }
                    catch (const std::invalid_argument&) { /* do nothing */ }
                    catch (const std::out_of_range&) {/* do nothing */}
            };

        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
            {
                para.input.out_pot[0] = 0;
            }
        };
        sync_intvec(input.out_pot, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_dmk");
        item.annotation = ">0 output density matrix DM(k) for each k-point";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = R"(Whether to output the density matrix for each k-point into files in the folder OUT.${suffix}. For current develop versions, out_dmk writes *_nao.txt files and includes a g{istep} index in the file name:
    * For gamma only case:
     * nspin = 1 and 4: dmg1_nao.txt;
     * nspin = 2: dms1g1_nao.txt and dms2g1_nao.txt for the two spin channels.
    * For multi-k points case:
     * nspin = 1 and 4: dmk1g1_nao.txt, dmk2g1_nao.txt, ...;
     * nspin = 2: dmk1s1g1_nao.txt... and dmk1s2g1_nao.txt... for the two spin channels.

    Here, g{istep} denotes the geometry/step index in the output file name.

    [NOTE] Version difference (develop vs 3.10-LTS):
    * In develop, out_dmk supports both gamma-only and multi-k-point density-matrix output.
    * In 3.10-LTS, the corresponding keyword is out_dm, and the output files are SPIN1_DM and SPIN2_DM, etc.)";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao");
			item.read_value = [](const Input_Item& item, Parameter& para) {
				const size_t count = item.get_size();
				if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_dmk needs at least 1 value");
				para.input.out_dmk[0] = assume_as_boolean(item.str_values[0]);
            para.input.out_dmk[1] = 8;
			if (count >= 2) try { para.input.out_dmk[1] = std::stoi(item.str_values[1]); }
			catch (const std::invalid_argument&) { /* do nothing */ }
			catch (const std::out_of_range&) {/* do nothing */}
			};
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
            {
                para.input.out_dmk[0] = 0;
            }
        };
        sync_intvec(input.out_dmk, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_dmr");
	    item.annotation = "output density matrix DM(R) with respect to lattice vector R (with precision 8)";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = R"(Whether to output the density matrix with Bravias lattice vector R index into files in the folder OUT.${suffix}. The files are named as dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in out_hsr. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps:
* nspin = 1: dmrs1_nao.csr;
* nspin = 2: dmrs1_nao.csr and dmrs2_nao.csr for the two spin channels.

[NOTE] In the 3.10-LTS version, the parameter is named out_dm1, and the file names are data-DMR-sparse_SPIN0.csr and data-DMR-sparse_SPIN1.csr, etc.)";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao and gamma_only==0");
        item.read_value = [](const Input_Item& item, Parameter& para) {
		    const size_t count = item.get_size();
		    if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_dmr needs at least 1 value");
		    para.input.out_dmr[0] = assume_as_boolean(item.str_values[0]);
		    para.input.out_dmr[1] = 8;
		    if (count >= 2) try { para.input.out_dmr[1] = std::stoi(item.str_values[1]); }
		    catch (const std::invalid_argument&) { /* do nothing */ }
		    catch (const std::out_of_range&) {/* do nothing */}
	    };
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
            {
                para.input.out_dmr[0] = 0;
            }
        };
	    item.check_value = [](const Input_Item& item, const Parameter& para) {
		    if (para.sys.gamma_only_local == true && para.input.out_dmr[0])
		    {
			    ModuleBase::WARNING_QUIT("ReadInput", "out_dmr is only valid for multi-k calculation");
		    }
	    };

	    sync_intvec(input.out_dmr, 2, 0);
	    this->add_item(item);
    }
    {
        Input_Item item("out_wfc_pw");
        item.annotation = "output wave functions";
        item.category = "Output information";
        item.type = "Integer";
        item.description = R"(Whether to output the electronic wavefunction coefficients into files and store them in the folder OUT.${suffix}. The files are named as wf{k}{k-point index}{s}{spin index}{g}{geometry index}{e}{electronic iteration index}{_pw} + {".txt"/".dat"}. Here, the s index refers to spin but the label will not show up for non-spin-polarized calculations, where s1 means spin up channel while s2 means spin down channel, and s4 refers to spinor wave functions that contains both spin channels with spin-orbital coupling or noncollinear calculations enabled. For scf or nscf calculations, g index will not appear, but the g index appears for geometry relaxation and molecular dynamics, where one can use the out_freq_ion command to control. To print out the electroinc wave functions every few SCF iterations, use the out_freq_elec command and the e index will appear in the file name.
* 0: no output
* 1: (txt format)
 * non-gamma-only with nspin=1: wfk1_pw.txt, wfk2_pw.txt, ...;
 * non-gamma-only with nspin=2: wfk1s1_pw.txt, wfk1s2_pw.txt, wfk2s1_pw.txt, wfk2s2_pw.txt, ...;
 * non-gamma-only with nspin=4: wfk1s4_pw.txt, wfk2s4_pw.txt, ...;
* 2: (binary format)
 * non-gamma-only with nspin=1: wfk1_pw.dat, wfk2_pw.dat, ...;
 * non-gamma-only with nspin=2: wfk1s1_pw.dat, wfk1s2_pw.dat, wfk2s1_pw.dat, wfk2s2_pw.dat, ...;
 * non-gamma-only with nspin=4: wfk1s4_pw.dat, wfk2s4_pw.dat, ...;

[NOTE] In the 3.10-LTS version, the file names are WAVEFUNC1.dat, WAVEFUNC2.dat, etc.)";
        item.default_value = "0";
        item.unit = "";
        item.set_availability("basis_type==pw or (basis_type==lcao and calculation==get_wf)");
        read_sync_int(input.out_wfc_pw);
        this->add_item(item);
    }
    {
        Input_Item item("out_wfc_lcao");
        item.annotation = "ouput LCAO wave functions, 0, no output 1: text, 2: binary";
        item.category = "Output information";
        item.type = "Integer";
        item.description = R"(Whether to output the electronic wavefunction coefficients into files and store them in the folder OUT.${suffix}. The files are named as wf{s}{spin index}{k(optional)}{k-point index}{g(optional)}{geometry index1}{_nao} + {".txt"/".dat"}. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and 's12' refer to spinor wave functions that contains both spin channels with spin-orbital coupling or noncollinear calculations enabled. In addition, if 'gamma_only' is set to 0, then the optinoal k-point sampling index appears with the k-point index attached to the electronic wave function file names. Finally, if out_app_flag is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if out_app_flag is set to true, the wave functions accumulate during ionic steps. If the out_app_flag is set to false, a new folder named WFC will be created, and the wave function files will be saved into it.
* 0: no output
* 1: (txt format)
 * gamma-only: wfs1_nao.txt or wfs2_nao.txt, ...;
 * non-gamma-only: wfs1k1_nao.txt or wfs1k2_nao.txt, ...;
* 2: (binary format)
 * gamma-only: wfs1_nao.dat or wfs2_nao.dat, ...;
 * non-gamma-only: wfs1k1_nao.dat or wfs1k2_nao.dat, ....

The corresponding sequence of the orbitals can be seen in Basis Set.

Also controled by out_freq_ion and out_app_flag.

[NOTE] In the 3.10-LTS version, the file names are WFC_NAO_GAMMA1_ION1.txt and WFC_NAO_K1_ION1.txt, etc.)";
        item.default_value = "0";
        item.unit = "";
        item.set_availability("basis_type==lcao");
        read_sync_int(input.out_wfc_lcao);
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.qo_switch)
            {
                para.input.out_wfc_lcao = 1;
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_wfc_lcao < 0 || para.input.out_wfc_lcao > 2)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_wfc_lcao should be 0, 1, or 2");
            }
            if (para.input.basis_type != "lcao" && para.input.out_wfc_lcao != 0)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_wfc_lcao is only available for basis_type = lcao");
            }
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_dos");
        item.annotation = "output energy and dos";
        item.category = "Output information";
        item.type = "Integer";
        item.description = R"(Whether to output the density of states (DOS). For more information, refer to the dos.md.
* 0: no output
* 1: output the density of states (DOS)
 * nspin=1 or 4: doss1g{geom}_{basis}.txt, where geom is the geometry index when cell changes or ions move while basis is either pw or nao.
 * nspin=2: doss1g{geom}_{basis}.txt and doss2g{geom}_{basis}.txt for two spin channles.
* 2: (LCAO) output the density of states (DOS) and the projected density of states (PDOS)
* 3: output the Fermi surface file (fermi.bxsf) in BXSF format that can be visualized by XCrySDen)";
        item.default_value = "0";
        item.unit = "";
        read_sync_int(input.out_dos);
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
            {
                para.input.out_dos = 0;
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_dos == 3 && para.input.symmetry == "1")
            {
                ModuleBase::WARNING_QUIT("ReadInput",
                                         "symmetry can't be used for out_dos==3(Fermi Surface "
                                         "Plotting) by now.");
            }
            if (para.input.basis_type == "pw" && para.input.out_dos == 3)
            {
                ModuleBase::WARNING_QUIT("ReadInput",
                                         "Fermi Surface Plotting not "
                                         "implemented for plane wave now.");
            }
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_ldos");
        item.annotation = "output mode of local density of states, second parameter controls the precision";
        item.category = "Output information";
        item.type = R"(Integer \[Integer\](optional))";
        item.description = R"(Whether to output the local density of states (LDOS), optionally output precision can be set by a second parameter, default is 3.
* 0: no output
* 1: output the partial charge density for given bias (controlled by stm_bias) in cube file format, which can be used to plot scanning tunneling spectroscopys to mimick STM images using the Python script plot.py.
* 2: output LDOS along a line in real space (controlled by ldos_line). Parameters used to control DOS output are also valid for LDOS.
* 3: output both two LDOS modes above.)";
        item.default_value = "0";
        item.unit = "";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count != 1 && count != 2)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_ldos should have 1 or 2 values");
            }
            para.input.out_ldos[0] = std::stoi(item.str_values[0]);
            para.input.out_ldos[1] = (count == 2) ? std::stoi(item.str_values[1]) : 3;
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_ldos[0] < 0 || para.input.out_ldos[0] > 3)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_ldos should be 0, 1, 2 or 3");
            }
        };
        sync_intvec(input.out_ldos, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_band");
        item.annotation = "output energy and band structure (with precision 8)";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = R"(Whether to output the eigenvalues of the Hamiltonian matrix (in eV) into the running log during electronic iterations and into a file at the end of calculations. The former can be used with the 'out_freq_elec' parameter while the latter option allows the output precision to be set via a second parameter, with a default value of 8. The output file names are:
 * nspin = 1 or 4: eig.txt;
 * nspin = 2: eigs1.txt and eigs2.txt;
 * For more information, refer to the band.md)";
        item.default_value = "False";
        item.unit = "";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count != 1 && count != 2)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_band should have 1 or 2 values");
            }
            para.input.out_band[0] = assume_as_boolean(item.str_values[0]);
            para.input.out_band[1] = (count == 2) ? std::stoi(item.str_values[1]) : 8;
        };
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
            {
                para.input.out_band[0] = 0;
            }
        };
        sync_intvec(input.out_band, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_proj_band");
        item.annotation = "output projected band structure";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to output the projected band structure. For more information, refer to the band.md";
        item.default_value = "False";
        item.unit = "";
        read_sync_bool(input.out_proj_band);
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (para.input.calculation == "get_wf" || para.input.calculation == "get_pchg")
            {
                para.input.out_proj_band = false;
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.basis_type == "pw" && para.input.out_proj_band)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_proj_band is only for lcao");
            }
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_stru");
        item.annotation = "output the structure files per ion step";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Controls the output of structure files per ionic step in geometry relaxation calculations. The files are written to the OUT.{suffix}/ directory. Each file corresponds to the structure at RELAX STEP ${istep}, i.e., the structure for which that step's energy was computed (before the relax move), and includes a header comment with the ABACUS version, timestamp, energy, and stress tensor. When out_freq_ion is positive, the numbered files STRU{istep} (or STRU{istep}.cif) are written every out_freq_ion steps; when out_freq_ion is 0, no numbered files are output.\n"
                          "    - 0: No structure files are output.\n"
                          "    - 1: ABACUS STRU format files are output. The latest structure is written to STRU_NOW (overwritten each step), the numbered file STRU{istep} (e.g., STRU1, STRU2) is written every out_freq_ion steps (when out_freq_ion is positive), and the final converged structure is written to STRU_FINAL. No CIF files are output.\n"
                          "    - 2: CIF format files are output. The latest structure is written to STRU_NOW.cif (overwritten each step), the numbered file STRU{istep}.cif (e.g., STRU1.cif, STRU2.cif) is written every out_freq_ion steps (when out_freq_ion is positive), and the final converged structure is written to STRU_FINAL.cif. No non-CIF files are output.\n"
                          "[NOTE] For backward compatibility, true/false (case insensitive) are accepted and converted to 1/0.";
        item.default_value = "1";
        item.unit = "";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const std::string val = FmtCore::lower(item.str_values[0]);
            if (val == "true" || val == "t" || val == "yes" || val == "y" || val == "on" || val == ".true.")
            {
                para.input.out_stru = 1;
            }
            else if (val == "false" || val == "f" || val == "no" || val == "n" || val == "off" || val == ".false.")
            {
                para.input.out_stru = 0;
            }
            else
            {
                try
                {
                    size_t pos = 0;
                    const int parsed = std::stoi(item.str_values[0], &pos);
                    if (pos != item.str_values[0].size())
                    {
                        ModuleBase::WARNING_QUIT("ReadInput",
                            "out_stru must be one of 0, 1, 2. For backward compatibility, true/false are also accepted. Got: '" + item.str_values[0] + "'.");
                    }
                    para.input.out_stru = parsed;
                }
                catch (const std::invalid_argument&)
                {
                    ModuleBase::WARNING_QUIT("ReadInput",
                        "out_stru must be one of 0, 1, 2. For backward compatibility, true/false are also accepted. Got: '" + item.str_values[0] + "'.");
                }
                catch (const std::out_of_range&)
                {
                    ModuleBase::WARNING_QUIT("ReadInput",
                        "out_stru must be one of 0, 1, 2. For backward compatibility, true/false are also accepted. Got: '" + item.str_values[0] + "'.");
                }
            }
        };
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            const std::vector<std::string> offlist = {"nscf", "get_s", "get_pchg", "get_wf"};
            if (std::find(offlist.begin(), offlist.end(), para.input.calculation) != offlist.end())
            {
                para.input.out_stru = 0;
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_stru < 0 || para.input.out_stru > 2)
            {
                ModuleBase::WARNING_QUIT("ReadInput",
                    "out_stru must be one of 0, 1, 2. For backward compatibility, true/false are also accepted.");
            }
        };
        sync_int(input.out_stru);
        this->add_item(item);
    }
    {
        Input_Item item("out_level");
        item.annotation = "ie(for electrons); i(for ions);";
        item.category = "Output information";
        item.type = "String";
        item.description = R"(Control the output level of information in OUT.{calculation}.log.
* ie: electronic iteration level, which prints useful information for electronic iterations;
* i: geometry relaxation level, which prints some information for geometry relaxations additionally;
* m: molecular dynamics level, which does not print some information for simplicity.)";
        item.default_value = "ie";
        item.unit = "";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            para.input.out_level = strvalue;
            para.sys.out_md_control = true;
        };
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            if (!para.sys.out_md_control && para.input.calculation == "md")
            {
                para.input.out_level = "m"; // zhengdy add 2019-04-07
            }
        };
        sync_string(input.out_level);
        add_bool_bcast(sys.out_md_control);
        this->add_item(item);
    }
    {
        Input_Item item("out_hsk");
        item.annotation = "output H(k) and S(k) matrices in reciprocal space";
        item.category = "Output information";
        item.type = R"(Integer \[Integer\](optional))";
        item.description = R"(Output the upper triangular part of the Hamiltonian and overlap matrices in reciprocal space for each k-point into files in the directory OUT.${suffix}. The first integer selects the format:
* 0: disabled;
* 1: text output; the optional second integer controls precision and defaults to 8;
* 2: binary output in the native ABACUS .dat format;
* 3: NPZ output, which is not implemented for H(k)/S(k).

The output is also controlled by out_freq_ion and out_app_flag. For more information, refer to hs_matrix.md.
* Gamma-only, nspin = 1: hk_nao.txt for the Hamiltonian matrix and sk_nao.txt for the overlap matrix.
* Gamma-only, nspin = 2: hks1_nao.txt and hks2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk_nao.txt for the overlap matrix. Only one overlap matrix is written because it is identical for both spin channels.
* Gamma-only, nspin = 4: not available with the gamma-only algorithm.
* Multi-k, nspin = 1: hk1_nao.txt for the Hamiltonian matrix and sk1_nao.txt for the overlap matrix at the first k-point.
* Multi-k, nspin = 2: hk1s1_nao.txt and hk1s2_nao.txt for the two spin channels of the Hamiltonian matrix, and sk1_nao.txt for the overlap matrix at the first k-point. Only one overlap matrix is written because it is identical for both spin channels.
* Multi-k, nspin = 4: hk1s4_nao.txt for the spinor Hamiltonian matrix and sk1_nao.txt for the spinor overlap matrix at the first k-point.
For binary output, the same names use the .dat suffix. Each native binary record contains the matrix dimension as an int followed by the row-major upper triangle. Gamma-only elements are doubles; multi-k and spinor elements are pairs of doubles containing the real and imaginary parts. Native integer representation and byte order are used.
When out_app_flag is true, the first ionic step truncates the file and later steps append complete records.
When out_app_flag is false, g followed by the one-based ionic-step index is inserted before _nao, for example hk1s1g1_nao.txt.

[NOTE] In the 3.10-LTS version, the file names are data-0-H and data-0-S, etc.)";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1 || count > 2)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_hsk expects a format and optional precision");
            }
            try
            {
                para.input.out_hsk[0] = std::stoi(item.str_values[0]);
                para.input.out_hsk[1] = count == 2 ? std::stoi(item.str_values[1]) : 8;
            }
            catch (const std::exception&)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_hsk format and precision must be integers");
            }
            if (count == 2 && para.input.out_hsk[0] != 1)
            {
                ModuleBase::WARNING("ReadInput", "out_hsk precision is ignored unless format is 1");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            const int format = para.input.out_hsk[0];
            if (format < 0 || format > 3)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_hsk format must be 0, 1, 2, or 3");
            }
            if (format == 3)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_hsk NPZ output is not implemented");
            }
        };
        sync_intvec(input.out_hsk, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_hs");
        item.annotation = "legacy alias for text H(k) and S(k) output in reciprocal space";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Legacy alias for out_hsk 1, which outputs Hamiltonian and overlap matrices in reciprocal space for each k-point. The optional second integer controls text precision. If both out_hsk and out_mat_hs are present, out_hsk takes precedence.";
        item.default_value = "False 8";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_hs needs at least 1 value");
            para.input.out_mat_hs[0] = assume_as_boolean(item.str_values[0]);
            para.input.out_mat_hs[1] = 8;
            if (count >= 2) try { para.input.out_mat_hs[1] = std::stoi(item.str_values[1]); }
            catch (const std::invalid_argument&) { /* do nothing */ }
            catch (const std::out_of_range&) {/* do nothing */}
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_hsr");
        item.annotation = "output H(R) and S(R) matrices in real space";
        item.category = "Output information";
        item.type = R"(Integer \[Integer\](optional))";
        item.description = R"(Output Hamiltonian and overlap matrices in real space, indexed by the Bravais lattice vector R, in the directory OUT.${suffix}. The first integer selects the format:
* 0: disabled;
* 1: text CSR output; the optional second integer controls precision and defaults to 8;
* 2: native binary CSR output using .dat files;
* 3: NPZ output using hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz.

For multi-k calculations, the output contains the individual real-space blocks stored for the Bravais lattice vectors R. For gamma-only calculations, the internal real-space contributions are folded into a single R = (0, 0, 0) block. This folded result cannot recover the original R-resolved contributions or interpolate arbitrary k points. Terms added only while constructing H(k) are not guaranteed to be present.

For binary output, each file uses the same basename as text output with a .dat suffix. Every native record contains the zero-based ionic step, matrix dimension, and number of R blocks as ints. Each R block contains three int coordinates, an int nonzero count, native double values (real/imaginary double pairs for complex matrices), int column indices, and long long row pointers. Native integer representation and byte order are used. When out_app_flag is true, the first ionic step truncates the file and later steps append complete records.

[NOTE] In the 3.10-LTS version, the file names are data-HR-sparse_SPIN0.csr and data-SR-sparse_SPIN0.csr, etc.)";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1 || count > 2)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_hsr expects a format and optional precision");
            }
            try
            {
                para.input.out_hsr[0] = std::stoi(item.str_values[0]);
                para.input.out_hsr[1] = count == 2 ? std::stoi(item.str_values[1]) : 8;
            }
            catch (const std::exception&)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_hsr format and precision must be integers");
            }
            if (count == 2 && para.input.out_hsr[0] != 1)
            {
                ModuleBase::WARNING("ReadInput", "out_hsr precision is ignored unless format is 1");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            const int format = para.input.out_hsr[0];
            if (format < 0 || format > 3)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_hsr format must be 0, 1, 2, or 3");
            }
            if (format == 3)
            {
#ifndef __CNPY
                ModuleBase::WARNING_QUIT("ReadInput",
                                         "to write in npz format, please "
                                         "recompile with -DENABLE_CNPY=1");
#endif
            }
        };
        sync_intvec(input.out_hsr, 2, 0);
        add_bool_bcast(input.out_hsr_npz_compat);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_hs2");
        item.annotation = "legacy alias for text H(R) and S(R) output in real space";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Legacy alias for out_hsr 1, which outputs Hamiltonian and overlap matrices in real space indexed by the Bravais lattice vector R. The optional second integer controls text precision. If both out_hsr and out_mat_hs2 are present, out_hsr takes precedence.";
        item.default_value = "False 8";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_hs2 needs at least 1 value");
            para.input.out_mat_hs2[0] = assume_as_boolean(item.str_values[0]);
            para.input.out_mat_hs2[1] = 8;
            if (count >= 2) try { para.input.out_mat_hs2[1] = std::stoi(item.str_values[1]); }
            catch (const std::invalid_argument&) { /* do nothing */ }
            catch (const std::out_of_range&) {/* do nothing */}
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_tk");
        item.annotation = "output kinetic matrix of electrons T(k)";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Whether to print the upper triangular part of the kinetic matrices for each k-point into OUT.${suffix}/tks1ki_nao.txt, where i is the index of k points. One may optionally provide a second parameter to specify the precision."
                          "\n\n[NOTE] In the 3.10-LTS version, the file names are data-TR-sparse_SPIN0.csr, etc.";
        item.default_value = "False [8]";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao");
			item.read_value = [](const Input_Item& item, Parameter& para) {
				const size_t count = item.get_size();
				if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_tk needs at least 1 value");
				para.input.out_mat_tk[0] = assume_as_boolean(item.str_values[0]);
            para.input.out_mat_tk[1] = 8;
			if (count >= 2) try { para.input.out_mat_tk[1] = std::stoi(item.str_values[1]); }
			catch (const std::invalid_argument&) { /* do nothing */ }
			catch (const std::out_of_range&) {/* do nothing */}
        };
        sync_intvec(input.out_mat_tk, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_r");
        item.annotation = "output r(R) matrix";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Whether to print the matrix representation of the position matrix into files named rxrs1_nao.csr, ryrs1_nao.csr, rzrs1_nao.csr in the directory OUT.${suffix}. The optional second parameter controls text output precision. If calculation is set to get_s, the position matrix can be obtained without scf iterations. For more information, please refer to position_matrix.md."
                          "\n\n[NOTE] In the 3.10-LTS version, the file name is data-rR-sparse.csr.";
        item.default_value = "False 8";
        item.unit = "Bohr";
        item.set_availability("basis_type==lcao and gamma_only==0");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_r needs at least 1 value");
            try {
                para.input.out_mat_r[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_r[1] = 8;
                if (count >= 2) try { para.input.out_mat_r[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument& e) {
                    ModuleBase::WARNING("Input", "out_mat_r precision must be an integer, using default 8");
                }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_r enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if ((para.inp.out_mat_r[0] || para.inp.out_mat_t[0]
                 || para.inp.out_hr_npz || para.inp.out_dm_npz || para.inp.dm_to_rho)
                && para.sys.gamma_only_local)
            {
                ModuleBase::WARNING_QUIT("ReadInput",
                                         "output of r(R)/T(R), H(R)-only/DM(R) in NPZ format, "
                                         "or conversion from DM(R) to rho is not "
                                         "available for gamma only calculations");
            }
        };
        sync_intvec(input.out_mat_r, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_t");
        item.annotation = "output T(R) matrix";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Generate files containing the kinetic energy matrix. The optional second parameter controls text output precision. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be trs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag."
                          "\n\n[NOTE] In the 3.10-LTS version, the file name is data-TR-sparse_SPIN0.csr.";
        item.default_value = "False 8";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao and gamma_only==0");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_t needs at least 1 value");
            try {
                para.input.out_mat_t[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_t[1] = 8;
                if (count >= 2) try { para.input.out_mat_t[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument& e) {
                    ModuleBase::WARNING("Input", "out_mat_t precision must be an integer, using default 8");
                }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_t enable flag must be 0/1, using default 0");
            }
        };
        sync_intvec(input.out_mat_t, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_dh");
        item.annotation = "output Hamiltonian derivatives dH/dR matrices";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix and overlap matrix as mentioned in out_hsr. The name of the files will be dhrxs1_nao.csr, dhrys1_nao.csr, dhrzs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag."
                          "\n\nFormat: <enable> [precision] [iat1 iat2 ...]. The first value (0/1) enables/disables output. The second optional value sets the output precision (default: 8). Starting from the third value, 1-based atom indices can be listed to restrict output to derivatives with respect to those specific atoms only; if no atom indices are given, all atoms are written."
                          "\n\n[NOTE] In the 3.10-LTS version, the file name is data-dHRx-sparse_SPIN0.csr and so on.";
        item.default_value = "0 8";
        item.unit = "Ry/Bohr";
        item.set_availability("basis_type==lcao and gamma_only==0");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh needs at least 1 value");
            try {
                para.input.out_mat_dh[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_dh[1] = 8;
                if (count >= 2) try { para.input.out_mat_dh[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument& e) {
                    ModuleBase::WARNING("Input", "out_mat_dh precision must be an integer, using default 8");
                }
                for (size_t i = 2; i < count; ++i)
                    try { para.input.out_mat_dh.push_back(std::stoi(item.str_values[i]) - 1); }
                    catch (const std::invalid_argument&) {
                        ModuleBase::WARNING("Input", "out_mat_dh atom index must be an integer, skipping");
                    }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_dh enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_dh[0] && para.input.nspin == 4)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh is not available for nspin = 4");
            }
        };
        sync_intvec(input.out_mat_dh, para.input.out_mat_dh.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_dh_t");
        item.annotation = "output kinetic energy dH/dR (dT/dR) matrices";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the derivatives of the kinetic energy matrix dT/dR."
                          "\n\nSee out_mat_dh for format details (enable, precision, atom indices).";
        item.default_value = "0 8";
        item.unit = "Ry/Bohr";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_dh_t[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_dh_t[1] = 8;
                if (count >= 2) try { para.input.out_mat_dh_t[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument&) {
                    ModuleBase::WARNING("Input", "out_mat_dh_t precision must be an integer, using default 8");
                }
                for (size_t i = 2; i < count; ++i)
                    try { para.input.out_mat_dh_t.push_back(std::stoi(item.str_values[i]) - 1); }
                    catch (const std::invalid_argument&) {
                        ModuleBase::WARNING("Input", "out_mat_dh_t atom index must be an integer, skipping");
                    }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_dh_t enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_dh_t[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh_t is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_dh_t, para.input.out_mat_dh_t.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_dh_vl");
        item.annotation = "output local pseudopotential dH/dR (dV^L/dR) matrices";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the derivatives of the local pseudopotential matrix dV^L/dR."
                          "\n\nSee out_mat_dh for format details.";
        item.default_value = "0 8";
        item.unit = "Ry/Bohr";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_dh_vl[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_dh_vl[1] = 8;
                if (count >= 2) try { para.input.out_mat_dh_vl[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument&) {
                    ModuleBase::WARNING("Input", "out_mat_dh_vl precision must be an integer, using default 8");
                }
                for (size_t i = 2; i < count; ++i)
                    try { para.input.out_mat_dh_vl.push_back(std::stoi(item.str_values[i]) - 1); }
                    catch (const std::invalid_argument&) {
                        ModuleBase::WARNING("Input", "out_mat_dh_vl atom index must be an integer, skipping");
                    }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_dh_vl enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_dh_vl[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh_vl is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_dh_vl, para.input.out_mat_dh_vl.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_dh_vnl");
        item.annotation = "output nonlocal pseudopotential dH/dR (dV^NL/dR) matrices";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the derivatives of the nonlocal pseudopotential matrix dV^NL/dR."
                          "\n\nSee out_mat_dh for format details.";
        item.default_value = "0 8";
        item.unit = "Ry/Bohr";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_dh_vnl[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_dh_vnl[1] = 8;
                if (count >= 2) try { para.input.out_mat_dh_vnl[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument&) {
                    ModuleBase::WARNING("Input", "out_mat_dh_vnl precision must be an integer, using default 8");
                }
                for (size_t i = 2; i < count; ++i)
                    try { para.input.out_mat_dh_vnl.push_back(std::stoi(item.str_values[i]) - 1); }
                    catch (const std::invalid_argument&) {
                        ModuleBase::WARNING("Input", "out_mat_dh_vnl atom index must be an integer, skipping");
                    }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_dh_vnl enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_dh_vnl[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh_vnl is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_dh_vnl, para.input.out_mat_dh_vnl.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_dh_vh");
        item.annotation = "output Hartree dH/dR (dV^H/dR) matrices";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the derivatives of the Hartree matrix dV^H/dR."
                          "\n\nSee out_mat_dh for format details.";
        item.default_value = "0 8";
        item.unit = "Ry/Bohr";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_dh_vh[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_dh_vh[1] = 8;
                if (count >= 2) try { para.input.out_mat_dh_vh[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument&) {
                    ModuleBase::WARNING("Input", "out_mat_dh_vh precision must be an integer, using default 8");
                }
                for (size_t i = 2; i < count; ++i)
                    try { para.input.out_mat_dh_vh.push_back(std::stoi(item.str_values[i]) - 1); }
                    catch (const std::invalid_argument&) {
                        ModuleBase::WARNING("Input", "out_mat_dh_vh atom index must be an integer, skipping");
                    }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_dh_vh enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_dh_vh[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh_vh is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_dh_vh, para.input.out_mat_dh_vh.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_dh_vxc");
        item.annotation = "output exchange-correlation dH/dR (dV^XC/dR) matrices";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the derivatives of the XC matrix dV^XC/dR."
                          "\n\nSee out_mat_dh for format details.";
        item.default_value = "0 8";
        item.unit = "Ry/Bohr";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_dh_vxc[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_dh_vxc[1] = 8;
                if (count >= 2) try { para.input.out_mat_dh_vxc[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument&) {
                    ModuleBase::WARNING("Input", "out_mat_dh_vxc precision must be an integer, using default 8");
                }
                for (size_t i = 2; i < count; ++i)
                    try { para.input.out_mat_dh_vxc.push_back(std::stoi(item.str_values[i]) - 1); }
                    catch (const std::invalid_argument&) {
                        ModuleBase::WARNING("Input", "out_mat_dh_vxc atom index must be an integer, skipping");
                    }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_dh_vxc enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_dh_vxc[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh_vxc is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_dh_vxc, para.input.out_mat_dh_vxc.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_dh_exx");
        item.annotation = "output exact-exchange dH/dR (dV^EXX/dR) matrices";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the derivatives of the exact-exchange matrix dV^EXX/dR."
                          "\n\nSee out_mat_dh for format details.";
        item.default_value = "0 8";
        item.unit = "Ry/Bohr";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_dh_exx[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_dh_exx[1] = 8;
                if (count >= 2) try { para.input.out_mat_dh_exx[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument&) {
                    ModuleBase::WARNING("Input", "out_mat_dh_exx precision must be an integer, using default 8");
                }
                for (size_t i = 2; i < count; ++i)
                    try { para.input.out_mat_dh_exx.push_back(std::stoi(item.str_values[i]) - 1); }
                    catch (const std::invalid_argument&) {
                        ModuleBase::WARNING("Input", "out_mat_dh_exx atom index must be an integer, skipping");
                    }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_dh_exx enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_dh_exx[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_dh_exx is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_dh_exx, para.input.out_mat_dh_exx.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_h_t");
        item.annotation = "output kinetic energy T(R) matrix";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the kinetic energy matrix T(R) in CSR format."
                          "\n\nSee out_hsr for format details.";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_h_t[0] = assume_as_boolean(item.str_values[0]);
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_h_t enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_h_t[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_h_t is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_h_t, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_h_vnl");
        item.annotation = "output nonlocal pseudopotential Vnl(R) matrix";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the nonlocal pseudopotential matrix Vnl(R) in CSR format."
                          "\n\nSee out_hsr for format details.";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_h_vnl[0] = assume_as_boolean(item.str_values[0]);
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_h_vnl enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_h_vnl[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_h_vnl is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_h_vnl, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_h_vl");
        item.annotation = "output local pseudopotential Vl(R) matrix";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the local pseudopotential matrix Vl(R) in CSR format."
                          "\n\nSee out_hsr for format details.";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_h_vl[0] = assume_as_boolean(item.str_values[0]);
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_h_vl enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_h_vl[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_h_vl is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_h_vl, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_h_vh");
        item.annotation = "output Hartree Vh(R) matrix";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the Hartree matrix Vh(R) in CSR format."
                          "\n\nSee out_hsr for format details.";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_h_vh[0] = assume_as_boolean(item.str_values[0]);
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_h_vh enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_h_vh[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_h_vh is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_h_vh, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_h_vxc");
        item.annotation = "output exchange-correlation Vxc(R) matrix";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the XC matrix Vxc(R) in CSR format."
                          "\n\nSee out_hsr for format details.";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_h_vxc[0] = assume_as_boolean(item.str_values[0]);
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_h_vxc enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_h_vxc[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_h_vxc is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_h_vxc, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_h_exx");
        item.annotation = "output exact-exchange Vexx(R) matrix";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Whether to print files containing the exact-exchange matrix Vexx(R) in CSR format."
                          "\n\nSee out_hsr for format details.";
        item.default_value = "0 8";
        item.unit = "Ry";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            try {
                para.input.out_mat_h_exx[0] = assume_as_boolean(item.str_values[0]);
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_h_exx enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_h_exx[0] && para.input.nspin == 4)
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_h_exx is not available for nspin = 4");
        };
        sync_intvec(input.out_mat_h_exx, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_ds");
        item.annotation = "output of derivative of S(R) matrix";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Whether to print files containing the derivatives of the overlap matrix. The optional second parameter controls text output precision. The format will be the same as the overlap matrix as mentioned in out_mat_dh. The name of the files will be dsxrs1_nao.csr and so on. Also controled by out_freq_ion and out_app_flag. This feature can be used with calculation get_s."
                          "\n\n[NOTE] In the 3.10-LTS version, the file name is data-dSRx-sparse_SPIN0.csr and so on.";
        item.default_value = "False 8";
        item.unit = "Ry/Bohr";
        item.set_availability("basis_type==lcao and gamma_only==0");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_ds needs at least 1 value");
            try {
                para.input.out_mat_ds[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_ds[1] = 8;
                if (count >= 2) try { para.input.out_mat_ds[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument& e) {
                    ModuleBase::WARNING("Input", "out_mat_ds precision must be an integer, using default 8");
                }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_ds enable flag must be 0/1, using default 0");
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_mat_ds[0] && para.input.nspin == 4)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_mat_ds is not available for nspin = 4");
            }
        };
        sync_intvec(input.out_mat_ds, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_xc");
        item.annotation = "output exchange-correlation matrix in KS-orbital representation";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to print the upper triangular part of the exchange-correlation matrices in Kohn-Sham orbital representation: for each k point into files in the directory OUT.i_nao.txt, where {suffix}/vxc_out.dat. If EXX is calculated, the local and EXX part of band energy will also be printed in OUT.{suffix}/vxc_exx_out.dat, respectively. All the vxc_out.dat files contains 3 integers (nk, nspin, nband) followed by nk*nspin*nband lines of energy Hartree and eV."
                          "\n\n[NOTE] In the 3.10-LTS version, the file name is k-$k-Vxc and so on.";
        item.default_value = "False";
        item.unit = "Ry";
        item.set_availability("basis_type in [lcao, lcao_in_pw]");
        read_sync_bool(input.out_mat_xc);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_xc2");
        item.annotation = "output exchange-correlation matrix in NAO representation";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Whether to print the exchange-correlation matrices in numerical orbital representation: in CSR format in the directory OUT.${suffix}. The name of the files will be vxcrs1_nao.csr and so on."
                          "\n\n[NOTE] In the 3.10-LTS version, the file name is Vxc_R_spin$s and so on.";
        item.default_value = "False 8";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_xc2 needs at least 1 value");
            try {
                para.input.out_mat_xc2[0] = assume_as_boolean(item.str_values[0]);
                para.input.out_mat_xc2[1] = 8;
                if (count >= 2) try { para.input.out_mat_xc2[1] = std::stoi(item.str_values[1]); }
                catch (const std::invalid_argument& e) {
                    ModuleBase::WARNING("Input", "out_mat_xc2 precision must be an integer, using default 8");
                }
            }
            catch (const std::invalid_argument& e) {
                ModuleBase::WARNING("Input", "out_mat_xc2 enable flag must be 0/1, using default 0");
            }
        };
        sync_intvec(input.out_mat_xc2, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_mat_l");
        item.annotation = "output the expectation values of angular momentum operators";
        item.category = "Output information";
        item.type = R"(Boolean \[Integer\](optional))";
        item.description = "Whether to print the expectation value of the angular momentum operator , , and in the basis of the localized atomic orbitals. The files are named OUT.{suffix}_Lx.dat, OUT.{suffix}_Ly.dat, and OUT.{suffix}_Lz.dat. The second integer controls the precision of the output.";
        item.default_value = "False 8";
        item.unit = "";
        item.set_availability("basis_type==lcao");
			item.read_value = [](const Input_Item& item, Parameter& para) {
				const size_t count = item.get_size();
				if (count < 1) ModuleBase::WARNING_QUIT("ReadInput", "out_mat_l needs at least 1 value");
				para.input.out_mat_l[0] = assume_as_boolean(item.str_values[0]);
            para.input.out_mat_l[1] = 8;
			if (count >= 2) try { para.input.out_mat_l[1] = std::stoi(item.str_values[1]); }
			catch (const std::invalid_argument&) { /* do nothing */ }
			catch (const std::out_of_range&) {/* do nothing */}
        };
        sync_intvec(input.out_mat_l, 2, 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_xc_r");
        item.annotation = "if >=0, output the derivatives of exchange correlation in realspace, second parameter controls the precision";
        item.category = "Output information";
        item.type = R"(Integer \[Integer\](optional))";
        item.description = R"(The first integer controls whether to output the exchange-correlation (in Bohr^-3) on real space grids using Libxc to folder OUT.${suffix}:
* 0: rho, amag, sigma, exc
* 1: vrho, vsigma
* 2: v2rho2, v2rhosigma, v2sigma2
* 3: v3rho3, v3rho2sigma, v3rhosigma2, v3sigma3
* 4: v4rho4, v4rho3sigma, v4rho2sigma2, v4rhosigma3, v4sigma4 The meaning of the files is presented in Libxc

The second integer controls the precision of the charge density output, if not given, will use 3 as default.

The circle order of the charge density on real space grids is: x is the outer loop, then y and finally z (z is moving fastest).)";
        item.default_value = "-1 3";
        item.unit = "";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count >= 1)
            {
                para.input.out_xc_r[0] = std::stoi(item.str_values[0]);
            }
            if (count >= 2)
            {
                para.input.out_xc_r[1] = std::stoi(item.str_values[1]);
            }
        };
        // check value
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_xc_r[0] >= 0)
            {
#ifndef __LIBXC
                ModuleBase::WARNING_QUIT("ReadInput", "INPUT out_xc_r is only aviailable with Libxc");
#endif
            }
        };
        sync_intvec(input.out_xc_r, 2, -1);
        this->add_item(item);
    }
    {
        Input_Item item("out_eband_terms");
        item.annotation = "output the band energy terms separately";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to print the band energy terms separately in the file OUT.{term}_out.dat. The terms include the kinetic, pseudopotential (local + nonlocal), Hartree and exchange-correlation (including exact exchange if calculated).";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao");
        read_sync_bool(input.out_eband_terms);
        this->add_item(item);
    }
    {
        Input_Item item("out_hr_npz");
        item.annotation = "output H(R) matrix in npz format";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to print Hamiltonian matrices H(R) in NPZ format as hrs1_nao.npz and, for nspin = 2, hrs2_nao.npz. This feature does not work for gamma-only calculations.";
        item.default_value = "False";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao and gamma_only==0");
        read_sync_bool(input.out_hr_npz);
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_hr_npz)
            {
#ifndef __CNPY
                ModuleBase::WARNING_QUIT("ReadInput",
                                         "to write in npz format, please "
                                         "recompile with -DENABLE_CNPY=1");
#endif
            }
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_hsr_npz");
        item.annotation = "legacy alias for H(R) and S(R) NPZ output";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Legacy alias for out_hsr 3, writing hrs1_nao.npz, hrs2_nao.npz when needed, and sr_nao.npz. If both out_hsr and out_hsr_npz are present, out_hsr takes precedence. Gamma-only calculations write the folded R = (0, 0, 0) representation.";
        item.default_value = "False";
        item.unit = "Ry";
        item.set_availability("basis_type==lcao");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            para.input.out_hsr_npz = assume_as_boolean(item.str_values[0]);
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_hsr_npz)
            {
#ifndef __CNPY
                ModuleBase::WARNING_QUIT("ReadInput",
                                         "to write in npz format, please "
                                         "recompile with -DENABLE_CNPY=1");
#endif
            }
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_dm_npz");
        item.annotation = "output DM(R) matrix in npz format";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to print density matrices DM(R) in npz format. This feature does not work for gamma-only calculations.";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao and gamma_only==0");
        read_sync_bool(input.out_dm_npz);
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_dm_npz)
            {
#ifndef __CNPY
                ModuleBase::WARNING_QUIT("ReadInput",
                                         "to write in npz format, please "
                                         "recompile with -DENABLE_CNPY=1");
#endif
            }
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_mul");
        item.annotation = "mulliken charge or not";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to print the Mulliken population analysis result into OUT.${suffix}/mulliken.txt. In molecular dynamics calculations, the output frequency is controlled by out_freq_ion.";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao");
        read_sync_bool(input.out_mul);
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.basis_type == "pw" && para.input.out_mul)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "out_mul is only for lcao");
            }
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_app_flag");
        item.annotation = "whether output r(R), H(R), S(R), T(R), and dH(R) "
                          "matrices in an append manner during MD";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to output r(R), H(R), S(R), T(R), dH(R), dS(R), and wfc matrices in an append manner during molecular dynamics calculations. Check input parameters out_mat_r, out_hsr, out_mat_t, out_mat_dh, out_hsk and out_wfc_lcao for more information.";
        item.default_value = "true";
        item.unit = "";
        item.set_availability("basis_type==lcao and gamma_only==0");
        read_sync_bool(input.out_app_flag);
        this->add_item(item);
    }
    {
        Input_Item item("out_ndigits");
        item.annotation = "the length of decimal part of output data";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "Controls the length of decimal part of output data, such as charge density, Hamiltonian matrix, Overlap matrix and so on.";
        item.default_value = "8";
        item.unit = "";
        item.set_availability("basis_type==lcao and out_hsk==1");
        read_sync_int(input.out_ndigits);
        this->add_item(item);
    }
    {
        Input_Item item("out_element_info");
        item.annotation = "output (projected) wavefunction of each element";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Whether to print element information into files in the directory OUT.{element_label}, including pseudopotential and orbital information of the element (in atomic Ryberg units).";
        item.default_value = "False";
        item.unit = "";
        read_sync_bool(input.out_element_info);
        this->add_item(item);
    }
    {
        Input_Item item("restart_save");
        item.annotation = "print to disk every step for restart";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = R"(Whether to save charge density files per ionic step, which are used to restart calculations. According to the value of read_file_dir:
* auto: These files are saved in folder OUT.{read_file_dir}/restart/.

If EXX(exact exchange) is calculated (i.e. dft_fuctional==hse/hf/pbe0/scan0 or rpa==True), the Hexx(R) files for each processor will also be saved in the above folder, which can be read in EXX calculation with restart_load==True.)";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao");
        read_sync_bool(input.restart_save);
        this->add_item(item);
    }
    {
        Input_Item item("rpa");
        item.annotation = "true:generate output files used in rpa calculation; "
                          "false:(default)";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = "Generate output files used in rpa calculations."
                          "\n\n[NOTE] If symmetry is set to 1, additional files containing the necessary information for "
                          "exploiting symmetry in the subsequent rpa calculation will be output: "
                          "irreducible_sector.txt, symrot_k.txt and symrot_R.txt.";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao");
        read_sync_bool(input.rpa);
        this->add_item(item);
    }
    {
        Input_Item item("rpa_out_vel");
        item.annotation = "whether to output velocity matrix for librpa";
        item.category = "Output information";
        item.type = "Boolean";
        item.description = " Velocity matrix in KS basis (in unit of eV *Angstrom). Loop layer: spin -> k -> direction -> KS_basis1 -> KS_basis2.";
        item.default_value = "False";
        item.unit = "eV * A";
        item.set_availability("basis_type==lcao");
        read_sync_bool(input.rpa_out_vel);
        this->add_item(item);
    }
    {
        Input_Item item("rpa_outdir");
        item.annotation = "output directory for librpa";
        item.category = "Output information";
        item.type = "String";
        item.description = "The directory to save files for LibRPA.";
        item.default_value = "\"OUT.librpa\"";
        item.unit = "";
        item.set_availability("basis_type==lcao");
        read_sync_string(input.rpa_outdir);
        item.reset_value = [](const Input_Item& item, Parameter& para) {
            para.input.rpa_outdir = to_dir(para.input.rpa_outdir);
        };
        this->add_item(item);
    }
    {
        Input_Item item("out_pchg");
        item.annotation = "specify the bands to be calculated for the partial (band-decomposed) charge densities";
        item.category = "Output information";
        item.type = "String";
        item.description = R"(Selects electronic states for partial (band-decomposed) charge-density output using a space-separated string of `0`s and `1`s, where `1` selects a state and `0` skips it. Repetition follows the `ocp_set` syntax, for example `1 4*0 5*1 0`; the expanded list must not exceed `nbands`. Each output represents a complete one-particle state rather than its SCF occupation. The spin degeneracy is 2 for `nspin=1` and 1 for `nspin=2` or `nspin=4`. For `nspin=1`, `s1` contains the charge density. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down charge densities, respectively. For `nspin=4`, `s1`, `s2`, `s3`, and `s4` respectively contain $\rho_0$, $m_x$, $m_y$, and $m_z$. With `if_separate_k=true`, files are named `pchgi[state]s[component]k[kpoint].cube`; otherwise, the weighted k-point sum is named `pchgi[state]s[component].cube`.

[NOTE] Enabling symmetry may produce unintended partial charge densities because of reduced k-point weights and real-space symmetry operations. If the desired symmetry treatment is uncertain, set `symmetry = -1`. Use the same symmetry setting as in the SCF calculation.)";
        item.default_value = "none";
        item.unit = "";
        item.set_availability("basis_type==pw or (basis_type==lcao and calculation==get_pchg)");
        item.read_value
            = [](const Input_Item& item, Parameter& para) { parse_expression(item.str_values, para.input.out_pchg); };
        item.get_final_value = [](Input_Item& item, const Parameter& para) {
            if (item.is_read())
            {
                item.final_value.str(longstring(item.str_values));
            }
        };
        add_intvec_bcast(input.out_pchg, para.input.out_pchg.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_wfc_norm");
        item.annotation = "specify the bands to be calculated for the norm of wavefunctions";
        item.category = "Output information";
        item.type = "String";
        item.description = R"(Selects electronic states for real-space wavefunction-modulus output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction modulus. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunction moduli, respectively. For `nspin=4`, `s1` contains the total spinor modulus. Files are named `wfi[state]s[spin]k[kpoint].cube`.)";
        item.default_value = "none";
        item.unit = "";
        item.set_availability("basis_type==pw or (basis_type==lcao and calculation==get_wf)");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            parse_expression(item.str_values, para.input.out_wfc_norm);
        };
        item.get_final_value = [](Input_Item& item, const Parameter& para) {
            if (item.is_read())
            {
                item.final_value.str(longstring(item.str_values));
            }
        };
        add_intvec_bcast(input.out_wfc_norm, para.input.out_wfc_norm.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("out_wfc_re_im");
        item.annotation = "specify the bands to be calculated for the real and imaginary parts of wavefunctions";
        item.category = "Output information";
        item.type = "String";
        item.description = R"(Selects electronic states for real-space wavefunction real- and imaginary-part output using the selection syntax and complete-state normalization of `out_pchg`. For `nspin=1`, `s1` contains the wavefunction. For `nspin=2`, `s1` and `s2` contain the spin-up and spin-down wavefunctions, respectively. For `nspin=4`, `s1` and `s2` contain the upper and lower spinor components, respectively. Files are named `wfi[state]s[spin]k[kpoint][re/im].cube`.)";
        item.default_value = "none";
        item.unit = "";
        item.set_availability("basis_type==pw or (basis_type==lcao and calculation==get_wf)");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            parse_expression(item.str_values, para.input.out_wfc_re_im);
        };
        item.get_final_value = [](Input_Item& item, const Parameter& para) {
            if (item.is_read())
            {
                item.final_value.str(longstring(item.str_values));
            }
        };
        add_intvec_bcast(input.out_wfc_re_im, para.input.out_wfc_re_im.size(), 0);
        this->add_item(item);
    }
    {
        Input_Item item("if_separate_k");
        item.annotation = "specify whether to write the partial charge densities for all k-points to individual files "
                          "or merge them";
        item.category = "Output information";
        item.type = "Boolean";
        item.description
            = "Specifies whether to write partial charge densities for individual k-points or merge them.";
        item.default_value = "false";
        item.unit = "";
        item.set_availability("(basis_type==pw and out_pchg!=none) or (basis_type==lcao and calculation==get_pchg and gamma_only==0)");
        read_sync_bool(input.if_separate_k);
        this->add_item(item);
    }
    {
        Input_Item item("out_elf");
        item.annotation = "> 0 output electron localization function (ELF) for selected electron steps"
                          ", second parameter controls the precision, default is 3.";
        item.category = "Output information";
        item.type = R"(Integer \[Integer\](optional))";
        item.description = R"(Whether to output the electron localization function (ELF) in the folder `OUT.${suffix}`. The files are named as
* nspin = 1:
    * elftot.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$;
* nspin = 2:
    * elfs1.cube, elfs2.cube: ${\rm{ELF}}_\sigma = \frac{1}{1+\chi_\sigma^2}$, $\chi_\sigma = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i,\sigma}|^2} - \frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}$;
    * elftot.cube: ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i,\sigma}{f_i |\nabla\psi_{i,\sigma}|^2} - \sum_{\sigma}{\frac{|\nabla\rho_\sigma|^2}{8\rho_\sigma}}}{\sum_{\sigma}{\frac{3}{10}(6\pi^2)^{2/3}\rho_\sigma^{5/3}}}$;
* nspin = 4 (noncollinear):
    * elftot.cube: ELF for total charge density, ${\rm{ELF}} = \frac{1}{1+\chi^2}$, $\chi = \frac{\frac{1}{2}\sum_{i}{f_i |\nabla\psi_{i}|^2} - \frac{|\nabla\rho|^2}{8\rho}}{\frac{3}{10}(3\pi^2)^{2/3}\rho^{5/3}}$

When `out_freq_ion > 0`, a geometry step suffix `g{#}` is appended to the file names (e.g., `elftotg1.cube`, `elfs1g1.cube`).

The second integer controls the precision of the kinetic energy density output, if not given, will use 3 as default. For purpose restarting from this file and other high-precision involved calculation, recommend to use 10.

In molecular dynamics calculations, the output frequency is controlled by out_freq_ion.)";
        item.default_value = "0 3";
        item.unit = "";
        item.set_availability("esolver_type in [ksdft, ofdft]");
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count >= 1)
            {
                para.input.out_elf[0] = std::stoi(item.str_values[0]);
            }
            if (count >= 2)
            {
                para.input.out_elf[1] = std::stoi(item.str_values[1]);
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.out_elf[0] > 0 && para.input.esolver_type != "ksdft" && para.input.esolver_type != "ofdft" && para.input.esolver_type != "tddft")
            {
                ModuleBase::WARNING_QUIT("ReadInput", "ELF is only available for ksdft, ofdft and tddft");
            }
        };
        sync_intvec(input.out_elf, 2, 0);
        this->add_item(item);
    }
    {
        // refactored from the removal of wannier input file, ISSUE 6469
        Input_Item item("out_spillage");
        item.annotation = "output spillage of wavefunctions. This parameter only accepts 0 or 2.";
        item.category = "Output information";
        item.type = "Integer";
        item.description = "This output is only intentively needed by the ABACUS numerical atomic orbital generation workflow. This parameter is used to control whether to output the overlap integrals between truncated spherical Bessel functions (TSBFs) and plane-wave basis expanded wavefunctions (named as OVERLAP_Q), and between TSBFs (named as OVERLAP_Sq), also their first order derivatives. The output files are named starting with orb_matrix. A value of 2 would enable the output.";
        item.default_value = "0";
        item.unit = "";
        item.set_availability("esolver_type==ksdft and basis_type==pw");
        read_sync_int(input.out_spillage);
        this->add_item(item);
    }
    {
        Input_Item item("out_dipole");
        item.annotation = "output dipole or not";
        item.category = "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory";
        item.type = "Boolean";
        item.description = R"(Controls electric-dipole output. In RT-TDDFT, each enabled spin channel is written to OUT.{suffix}/dipole_s[spin].txt using a one-based spin number. Every row contains the one-based electronic-step index followed by the Cartesian electronic-dipole components $P_x$, $P_y$, and $P_z$ in atomic units. The running log additionally reports the electronic, ionic, and total dipoles and the norm of the total dipole.
* True: Output the electric dipole information.
* False: Do not output the electric dipole information.)";
        item.default_value = "False";
        item.unit = "";
        read_sync_bool(input.out_dipole);
        this->add_item(item);
    }
    {
        Input_Item item("out_current");
        item.annotation = "output current or not";
        item.category = "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory";
        item.type = "Integer";
        item.description = R"(Controls the current-density output method for LCAO RT-TDDFT. Output rows contain the one-based electronic-step index followed by $J_x$, $J_y$, and $J_z$ in atomic units.
* 0: Do not output current.
* 1: Explicitly construct the velocity operator from the momentum, vector-potential, and KB nonlocal-pseudopotential terms using two-center and spherical-grid integrals: $\hat{v}_{\alpha}=-\mathrm{i}\nabla_{\alpha}+A_{\alpha}(t)+\mathrm{i}\left[\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}},r_{\alpha}\right]$, where $\widetilde{V}_{\mathrm{NL}}^{\mathrm{KB}}=\mathrm{e}^{-\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}\hat{V}_{\mathrm{NL}}^{\mathrm{KB}}\mathrm{e}^{\mathrm{i}\boldsymbol{A}(t)\cdot\boldsymbol{r}}$. $\boldsymbol{A}(t)$ is nonzero only for the velocity gauge (td_stype=1); otherwise $\boldsymbol{A}(t)=0$. Other nonlocal Hamiltonian terms, such as EXX, are not included explicitly. The total current is written to OUT.{suffix}/current_tot.txt.
* 2: Use the full Hamiltonian to construct the generalized velocity matrix in a nonorthogonal NAO basis, $\widetilde{v}_{\alpha}=\partial_{\alpha}H+\mathrm{i}HS^{-1}\mathcal{R}_{\alpha}-\mathrm{i}\mathcal{R}_{\alpha}S^{-1}H-HS^{-1}\partial_{\alpha}S$. This includes all contributions available in the real-space Hamiltonian matrix when enabled. This method is more general but more expensive. The total current is written to OUT.{suffix}/current_tot_comm.txt.)";
        item.default_value = "0";
        item.unit = "";
        item.set_availability("basis_type==lcao and esolver_type==tddft");
        read_sync_int(input.out_current);
        this->add_item(item);
    }
    {
        Input_Item item("out_current_k");
        item.annotation = "output current for each k";
        item.category = "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory";
        item.type = "Boolean";
        item.description = R"(Controls whether LCAO RT-TDDFT current density is also resolved by spin and k-point. The total-current file is always written when out_current is 1 or 2.
* True: In addition to the total, out_current=1 writes OUT.{suffix}/current_s[spin]k[kpoint].txt; out_current=2 writes OUT.{suffix}/current_s[spin]k[kpoint]_comm.txt. Both use one-based spin and k-point numbers, with k-points numbered independently within each spin channel. Each row contains the one-based electronic-step index followed by $J_x$, $J_y$, and $J_z$ in atomic units.
* False: Output only current_tot.txt for out_current=1 or current_tot_comm.txt for out_current=2.)";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao and esolver_type==tddft and out_current>0");
        read_sync_bool(input.out_current_k);
        this->add_item(item);
    }
    {
        Input_Item item("out_efield");
        item.annotation = "output dipole or not";
        item.category = "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory";
        item.type = "Boolean";
        item.description = R"(Controls time-dependent electric-field output. For each configured field, OUT.{suffix}/efield_[index].txt contains two columns: physical time in fs and the field value in V/Angstrom. The one-based field index follows the occurrence order shared by td_ttype and td_vext_dire, so fields assigned to the same direction remain in separate files. At initialization, a fresh calculation with md_restart=False truncates the files corresponding to the currently configured fields, whereas a calculation with md_restart=True preserves them and appends new samples.
* True: Output electric-field values on active electronic steps.
* False: Do not output electric-field values.)";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("esolver_type==tddft and td_vext==true");
        read_sync_bool(input.out_efield);
        this->add_item(item);
    }
    {
        Input_Item item("out_vecpot");
        item.annotation = "output TDDFT vector potential or not";
        item.category = "RT-TDDFT: Real-Time Time-Dependent Density Functional Theory";
        item.type = "Boolean";
        item.description = R"(Controls Cartesian vector-potential output for LCAO RT-TDDFT. OUT.{suffix}/vector_pot.txt contains four columns: the one-based electronic-step index followed by $A_x$, $A_y$, and $A_z$ in atomic units. At initialization, a fresh calculation with md_restart=False truncates the file and writes a new header, whereas a calculation with md_restart=True preserves a nonempty existing file and appends new samples. If the restart output file is missing or empty, a new file with a header is created.
* True: Write vector-potential samples on electronic propagation steps.
* False: Do not output the vector potential.)";
        item.default_value = "False";
        item.unit = "";
        item.set_availability("basis_type==lcao and esolver_type==tddft");
        read_sync_bool(input.out_vecpot);
        this->add_item(item);
    }
    {
        // recover the functionality of test_symmetry by introducing a new keyword "out_symm_mat"
        // the "out_symm_mat" keyword will be a
        Input_Item item("cal_symm_repr");
        item.annotation = "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.";
        item.category = "System variables";
        item.type = R"(Integer \[Integer\](optional))";
        item.description = "Whether to print the matrix representation of symmetry operation to running log file. If the first value is given as 1, then all matrix representations will be printed. The second optional parameter controls the precision (number of digits) to print, default is 3, which is enough for a quick check.";
        item.default_value = "1 3";
        item.unit = "";
        item.read_value = [](const Input_Item& item, Parameter& para) {
            const size_t count = item.get_size();
            if (count >= 1)
            {
                para.input.cal_symm_repr[0] = std::stoi(item.str_values[0]);
            }
            if (count >= 2)
            {
                para.input.cal_symm_repr[1] = std::stoi(item.str_values[1]);
            }
        };
        item.check_value = [](const Input_Item& item, const Parameter& para) {
            if (para.input.cal_symm_repr[0] < 0 || para.input.cal_symm_repr[0] > 1)
            {
                ModuleBase::WARNING_QUIT("ReadInput", "cal_symm_repr should be 0 or 1");
            }
        };
        sync_intvec(input.cal_symm_repr, 2, 0);
        this->add_item(item);
    }
    {
        // refactored from the removal of wannier input file, ISSUE 6469
        Input_Item item("spillage_outdir");
        item.annotation = "output directory for spillage of wavefunctions.";
        item.category = "Input files";
        item.type = "String";
        item.description = "The directory to save the spillage files.";
        item.default_value = "\"./\"";
        item.unit = "";
        item.set_availability("basis_type==pw");
        read_sync_string(input.spillage_outdir);
        this->add_item(item);
    }
}
} // namespace ModuleIO