diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index e76826128c..c6bab7b036 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -1782,7 +1782,7 @@ ### out_freq_ion - **Type**: Integer -- **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, and Mulliken population analysis is printed every n ionic steps. +- **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. > Note: In RT-TDDFT calculations, this parameter is inactive; output frequency is instead controlled by out_freq_td. - **Default**: 0 @@ -1969,9 +1969,13 @@ ### out_stru -- **Type**: Boolean -- **Description**: Whether to output structure files per ionic step in geometry relaxation calculations into OUT.{istep}_D, where ${istep} is the ionic step. -- **Default**: False +- **Type**: Integer +- **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. + - 0: No structure files are output. + - 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. + - 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. + > Note: For backward compatibility, true/false (case insensitive) are accepted and converted to 1/0. +- **Default**: 1 ### out_level diff --git a/docs/advanced/input_files/stru.md b/docs/advanced/input_files/stru.md index a30064351a..37ae64d871 100644 --- a/docs/advanced/input_files/stru.md +++ b/docs/advanced/input_files/stru.md @@ -116,10 +116,10 @@ For general usage requirements, the APNSv1.0 pseudopotential and orbital set is ### LATTICE_CONSTANT - The lattice constant of the system in unit of Bohr. + The lattice constant of the system in unit of Bohr. In output structure files (e.g., `STRU`, `STRU1`, `STRU2`), a trailing comment `# in Bohr` is appended to the value line. ### LATTICE_VECTORS - The lattice vectors of the unit cell. It is a 3by3 matrix written in 3 lines. Please note that *the lattice vectors given here are scaled by the lattice constant*. This section must be removed if the type Bravais lattice is specified using the input parameter `latname`. (See [input parameters](input-main.md#latname).) + The lattice vectors of the unit cell. It is a 3by3 matrix written in 3 lines. Please note that *the lattice vectors given here are scaled by the lattice constant*. This section must be removed if the type Bravais lattice is specified using the input parameter `latname`. (See [input parameters](input-main.md#latname).) In output structure files, a trailing comment `# in units of lat0` is appended to the section header. ### LATTICE_PARAMETERS This section is only relevant when `latname` (see [input parameters](input-main.md#latname)) is used to specify the Bravais lattice type. The example above is a fcc lattice, where no additional information except the lattice constant is required to determine the geometry of the lattice. diff --git a/docs/parameters.yaml b/docs/parameters.yaml index c21f046f02..c9ffc651ed 100644 --- a/docs/parameters.yaml +++ b/docs/parameters.yaml @@ -2816,7 +2816,7 @@ parameters: category: Output information type: Integer 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, and Mulliken population analysis is printed every n ionic steps. + 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. [NOTE] In RT-TDDFT calculations, this parameter is inactive; output frequency is instead controlled by out_freq_td. default_value: "0" @@ -3024,10 +3024,14 @@ parameters: availability: "" - name: out_stru category: Output information - type: Boolean + type: Integer description: | - Whether to output structure files per ionic step in geometry relaxation calculations into OUT.{istep}_D, where ${istep} is the ionic step. - default_value: "False" + 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. + - 0: No structure files are output. + - 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. + - 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. + [NOTE] For backward compatibility, true/false (case insensitive) are accepted and converted to 1/0. + default_value: "1" unit: "" availability: "" - name: out_level diff --git a/docs/quick_start/hands_on.md b/docs/quick_start/hands_on.md index fe0111ea4c..ab280b233c 100644 --- a/docs/quick_start/hands_on.md +++ b/docs/quick_start/hands_on.md @@ -211,7 +211,7 @@ stress_thr 5 # the threshold of the stress convergence, in unit of kBar relax_nmax 100 # the maximal number of ionic iteration steps out_stru 1 ``` -Use the same `KPT`, `STRU`, pseudopotential, and orbital files as in the above SCF-LCAO example. The final optimized structure can be found in `STRU_NOW.cif` and `OUT.MgO/running_cell-relax.log`. +Use the same `KPT`, `STRU`, pseudopotential, and orbital files as in the above SCF-LCAO example. The final optimized structure can be found in `STRU_FINAL` and `OUT.MgO/running_cell-relax.log`. ### A quick PW example @@ -232,4 +232,4 @@ relax_nmax 100 # the maximal number of ionic iteration steps out_stru 1 ``` -Use the same `KPT`, `STRU`, and pseudopotential files as in the above SCF-PW examples. The final optimized structure can be found in `STRU_NOW.cif` and `STRU_ION_D` with different format. +Use the same `KPT`, `STRU`, and pseudopotential files as in the above SCF-PW examples. The final optimized structure can be found in `STRU_FINAL` and `STRU` with different format. diff --git a/examples/17_relax/README b/examples/17_relax/README index 2afedf336e..051032e96b 100644 --- a/examples/17_relax/README +++ b/examples/17_relax/README @@ -19,5 +19,5 @@ set `relax_method` to `cg`(default value) (3)`relax_nmax`: number of ion iteration steps; `force_thr_ev`: force threshold, unit: eV/Angstrom; \ `stress_ev`: stress threshold, unit: kBar. -(4)Output file OUT.ABACUS/STRU_NOW.cif contains the optimized atom positions. +(4)Output file OUT.ABACUS/STRU_FINAL contains the optimized atom positions. diff --git a/source/source_cell/print_cell.cpp b/source/source_cell/print_cell.cpp index 0c71b8bdbb..a83ec9a5e4 100644 --- a/source/source_cell/print_cell.cpp +++ b/source/source_cell/print_cell.cpp @@ -79,7 +79,8 @@ namespace unitcell void print_stru_file(const UnitCell& ucell, const Atom* atoms, const ModuleBase::Matrix3& latvec, - const std::string& fn, + const std::string& fn, + const std::string& header, const int& nspin, const bool& direct, const bool& vel, @@ -89,12 +90,18 @@ namespace unitcell const int& iproc) { ModuleBase::TITLE("UnitCell","print_stru_file"); - if (iproc != 0) + if (iproc != 0) { return; // old: if(GlobalV::MY_RANK != 0) return; } + // optional header comments + std::string str; + if (!header.empty()) + { + str = header; + } // ATOMIC_SPECIES - std::string str = "ATOMIC_SPECIES\n"; + str += "ATOMIC_SPECIES\n"; for(int it=0; itatoms,ucell->latvec, - fn, 1, false, false, false, false, false, 0); + fn, "", 1, false, false, false, false, false, 0); std::ifstream ifs; ifs.open("C1H2_STRU"); std::string str((std::istreambuf_iterator(ifs)), std::istreambuf_iterator()); @@ -829,9 +829,9 @@ TEST_F(UcellTest, PrintSTRU) EXPECT_THAT(str, testing::HasSubstr("LATTICE_CONSTANT")); EXPECT_THAT(str, testing::HasSubstr("1.8897261255")); EXPECT_THAT(str, testing::HasSubstr("LATTICE_VECTORS")); - EXPECT_THAT(str, testing::HasSubstr("10.0000000000 0.0000000000 0.0000000000")); - EXPECT_THAT(str, testing::HasSubstr(" 0.0000000000 10.0000000000 0.0000000000")); - EXPECT_THAT(str, testing::HasSubstr(" 0.0000000000 0.0000000000 10.0000000000")); + EXPECT_THAT(str, testing::HasSubstr("10.0000000000000000 0.0000000000000000 0.0000000000000000")); + EXPECT_THAT(str, testing::HasSubstr("0.0000000000000000 10.0000000000000000 0.0000000000000000")); + EXPECT_THAT(str, testing::HasSubstr("0.0000000000000000 0.0000000000000000 10.0000000000000000")); EXPECT_THAT(str, testing::HasSubstr("ATOMIC_POSITIONS")); EXPECT_THAT(str, testing::HasSubstr("Cartesian")); EXPECT_THAT(str, testing::HasSubstr("C #label")); @@ -851,7 +851,7 @@ TEST_F(UcellTest, PrintSTRU) * */ unitcell::print_stru_file(*ucell,ucell->atoms,ucell->latvec, - fn, 2, true, true, false, false, false, 0); + fn, "", 2, true, true, false, false, false, 0); ifs.open("C1H2_STRU"); str = {(std::istreambuf_iterator(ifs)), std::istreambuf_iterator()}; EXPECT_THAT(str, testing::HasSubstr("ATOMIC_SPECIES")); @@ -860,9 +860,9 @@ TEST_F(UcellTest, PrintSTRU) EXPECT_THAT(str, testing::HasSubstr("LATTICE_CONSTANT")); EXPECT_THAT(str, testing::HasSubstr("1.8897261255")); EXPECT_THAT(str, testing::HasSubstr("LATTICE_VECTORS")); - EXPECT_THAT(str, testing::HasSubstr("10.0000000000 0.0000000000 0.0000000000")); - EXPECT_THAT(str, testing::HasSubstr(" 0.0000000000 10.0000000000 0.0000000000")); - EXPECT_THAT(str, testing::HasSubstr(" 0.0000000000 0.0000000000 10.0000000000")); + EXPECT_THAT(str, testing::HasSubstr("10.0000000000000000 0.0000000000000000 0.0000000000000000")); + EXPECT_THAT(str, testing::HasSubstr("0.0000000000000000 10.0000000000000000 0.0000000000000000")); + EXPECT_THAT(str, testing::HasSubstr("0.0000000000000000 0.0000000000000000 10.0000000000000000")); EXPECT_THAT(str, testing::HasSubstr("ATOMIC_POSITIONS")); EXPECT_THAT(str, testing::HasSubstr("Direct")); EXPECT_THAT(str, testing::HasSubstr("C #label")); @@ -893,7 +893,7 @@ TEST_F(UcellTest, PrintSTRU) ucell->atom_mulliken = {{-1, 0.5}, {-1, 0.4}, {-1, 0.3}}; // first index is iat, the second is components, starts seems from 1 unitcell::print_stru_file(*ucell,ucell->atoms,ucell->latvec, - fn, 2, true, false, true, true, true, 0); + fn, "", 2, true, false, true, true, true, 0); ifs.open("C1H2_STRU"); str = {(std::istreambuf_iterator(ifs)), std::istreambuf_iterator()}; EXPECT_THAT(str, testing::HasSubstr("ATOMIC_SPECIES")); @@ -907,9 +907,9 @@ TEST_F(UcellTest, PrintSTRU) EXPECT_THAT(str, testing::HasSubstr("LATTICE_CONSTANT")); EXPECT_THAT(str, testing::HasSubstr("1.8897261255")); EXPECT_THAT(str, testing::HasSubstr("LATTICE_VECTORS")); - EXPECT_THAT(str, testing::HasSubstr("10.0000000000 0.0000000000 0.0000000000")); - EXPECT_THAT(str, testing::HasSubstr(" 0.0000000000 10.0000000000 0.0000000000")); - EXPECT_THAT(str, testing::HasSubstr(" 0.0000000000 0.0000000000 10.0000000000")); + EXPECT_THAT(str, testing::HasSubstr("10.0000000000000000 0.0000000000000000 0.0000000000000000")); + EXPECT_THAT(str, testing::HasSubstr("0.0000000000000000 10.0000000000000000 0.0000000000000000")); + EXPECT_THAT(str, testing::HasSubstr("0.0000000000000000 0.0000000000000000 10.0000000000000000")); EXPECT_THAT(str, testing::HasSubstr("ATOMIC_POSITIONS")); EXPECT_THAT(str, testing::HasSubstr("Direct")); EXPECT_THAT(str, testing::HasSubstr("C #label")); diff --git a/source/source_esolver/esolver_dp.cpp b/source/source_esolver/esolver_dp.cpp index 2c880f343e..9e3aded558 100644 --- a/source/source_esolver/esolver_dp.cpp +++ b/source/source_esolver/esolver_dp.cpp @@ -20,7 +20,6 @@ #include "esolver_dp.h" #include "source_base/parallel_common.h" #include "source_base/timer.h" -#include "source_io/module_output/cif_io.h" #include "source_io/module_output/output_log.h" #include "source_io/module_parameter/parameter.h" @@ -39,11 +38,6 @@ void ESolver_DP::before_all_runners(BaseCell& basecell, const Input_para& inp) dp_force.create(ucell.nat, 3); dp_virial.create(3, 3); - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", - ucell, - "# Generated by ABACUS ModuleIO::CifParser", - "data_?"); - atype.resize(ucell.nat); rescaling = inp.mdp.dp_rescaling; diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 754471313d..fc5e24b316 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -6,7 +6,6 @@ #include "source_estate/param_update.h" #include "source_hamilt/module_ewald/H_Ewald_pw.h" #include "source_hamilt/module_vdw/vdw.h" -#include "source_io/module_output/cif_io.h" #include "source_io/module_output/output_log.h" #include "source_io/module_output/print_info.h" #include "source_io/module_chgpot/rhog_io.h" @@ -73,14 +72,10 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp) //! 3) setup structure factors this->sf.set(this->pw_rhod, inp.nbspline); - //! 4) write geometry file - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", - ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?"); - - //! 5) init charge extrapolation + //! 4) init charge extrapolation this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap); - //! 6) symmetry analysis should be performed every time the cell is changed + //! 5) symmetry analysis should be performed every time the cell is changed if (ModuleSymmetry::Symmetry::symm_flag == 1) { const int cal_symm_repr[2] = {PARAM.inp.cal_symm_repr[0], PARAM.inp.cal_symm_repr[1]}; diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp index a8e8f40282..30b8ff78f5 100644 --- a/source/source_esolver/esolver_lj.cpp +++ b/source/source_esolver/esolver_lj.cpp @@ -3,7 +3,6 @@ #include "source_cell/module_neighbor/sltk_atom_arrange.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_io/module_output/output_log.h" -#include "source_io/module_output/cif_io.h" #include "source_cell/module_neighlist/neighbor_types.h" #include "source_cell/module_neighlist/neighbor_search.h" #include "source_base/global_variable.h" @@ -51,11 +50,6 @@ void ESolver_LJ::before_all_runners(BaseCell& cell, const Input_para& inp) lj_force.create(ucell.nat, 3); lj_virial.create(3, 3); - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", - ucell, - "# Generated by ABACUS ModuleIO::CifParser", - "data_?"); - // determine the maximum rcut and lj_rcut rcut_search_radius(ucell.ntype, inp.mdp.lj_rcut); diff --git a/source/source_esolver/esolver_nep.cpp b/source/source_esolver/esolver_nep.cpp index 1cb1aec63d..2adf610721 100644 --- a/source/source_esolver/esolver_nep.cpp +++ b/source/source_esolver/esolver_nep.cpp @@ -18,7 +18,6 @@ #include "esolver_nep.h" #include "source_base/parallel_common.h" #include "source_base/timer.h" -#include "source_io/module_output/cif_io.h" #include "source_io/module_output/output_log.h" #include "source_io/module_parameter/parameter.h" @@ -40,11 +39,6 @@ void ESolver_NEP::before_all_runners(BaseCell& basecell, const Input_para& inp) _f.resize(3 * ucell.nat); _v.resize(9 * ucell.nat); - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", - ucell, - "# Generated by ABACUS ModuleIO::CifParser", - "data_?"); - #ifdef __NEP /// determine the type map from STRU to NEP model type_map(ucell); diff --git a/source/source_estate/elecstate_print.cpp b/source/source_estate/elecstate_print.cpp index e7d9c9ee8e..51a7765463 100644 --- a/source/source_estate/elecstate_print.cpp +++ b/source/source_estate/elecstate_print.cpp @@ -386,9 +386,9 @@ void print_etot(const Magnetism& magnet, : PARAM.inp.ks_solver; elecstate::print_scf_iterinfo(iter_label, iter, - 6, + 4, mag, - 10, + 9, elec.f_en.etot * ModuleBase::Ry_to_eV, elec.f_en.etot_delta * ModuleBase::Ry_to_eV, 16, diff --git a/source/source_hsolver/diago_iter_assist.cpp b/source/source_hsolver/diago_iter_assist.cpp index c68dd4e5af..b547febefe 100644 --- a/source/source_hsolver/diago_iter_assist.cpp +++ b/source/source_hsolver/diago_iter_assist.cpp @@ -223,7 +223,9 @@ void DiagoIterAssist::diag_subspace_init(hamilt::Hamilt* p if (base_device::get_device_type(ctx) == base_device::GpuDevice) { - psi::Psi psi_temp(1, 1, psi_nc, dmin, true); + /// NOTE: current_nbasis must be npw (WITHOUT npol) for Nonlocal::act's + /// gemm K (vkb only has npw rows). See CPU branch comment above. + psi::Psi psi_temp(1, 1, psi_nc, evc.get_current_nbas(), true); T* ppsi = psi_temp.get_pointer(); // hpsi and spsi share the temp space @@ -270,7 +272,11 @@ void DiagoIterAssist::diag_subspace_init(hamilt::Hamilt* p } else if (base_device::get_device_type(ctx) == base_device::CpuDevice) { - psi::Psi psi_temp(1, nstart, psi_nc, dmin, true); + /// NOTE: the 4th arg (current_nbasis) must be npw (WITHOUT npol), + /// NOT dmin (= nbasis = npol*npwx in SOC). Nonlocal::act uses + /// psi_temp.get_current_nbas() as gemm K, but vkb only has npw rows. + /// dmin (still = nbasis) is kept for hcc/scc gemm K which needs npol. + psi::Psi psi_temp(1, nstart, psi_nc, evc.get_current_nbas(), true); T* ppsi = psi_temp.get_pointer(); syncmem_complex_op()(ppsi, psi, psi_temp.size()); @@ -295,7 +301,6 @@ void DiagoIterAssist::diag_subspace_init(hamilt::Hamilt* p delmem_complex_op()(temp); add_to_hcc(hcc, nstart); - } if (GlobalV::NPROC_IN_POOL > 1) diff --git a/source/source_io/module_output/print_info.cpp b/source/source_io/module_output/print_info.cpp index afc49757bf..7cf32256ee 100644 --- a/source/source_io/module_output/print_info.cpp +++ b/source/source_io/module_output/print_info.cpp @@ -376,11 +376,11 @@ void print_screen(const int& stress_step, const int& force_step, const int& iste else if(PARAM.inp.calculation=="cell-relax") { std::cout << " RELAX STEP: " << unsigned(istep); - std::cout << " (CELL_CHANGE# " << unsigned(stress_step); - std::cout << " IONS_CHANGE# " << unsigned(force_step) << ")" << std::endl; + std::cout << " (CELL# " << unsigned(stress_step); + std::cout << " IONS# " << unsigned(force_step) << ")" << std::endl; GlobalV::ofs_running << " RELAX STEP: " << unsigned(istep); - GlobalV::ofs_running << " (CELL_CHANGE# " << unsigned(stress_step); - GlobalV::ofs_running << " IONS_CHANGE# " << unsigned(force_step) << ")" << std::endl; + GlobalV::ofs_running << " (CELL# " << unsigned(stress_step); + GlobalV::ofs_running << " IONS# " << unsigned(force_step) << ")" << std::endl; } } diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index 8b7dae8f45..45f8d65b49 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -377,7 +377,8 @@ struct Input_para std::vector aims_nbasis = {}; ///< the number of basis functions for each atom type used in FHI-aims (for benchmark) // ============== #Parameters (11.Output) =========================== - bool out_stru = false; ///< outut stru file each ion step + int out_stru = 1; ///< output stru file each ion step + ///< 0: no output, 1: STRU format, 2: CIF format int out_freq_elec = 0; ///< the frequency of electronic iter to output charge and wavefunction int out_freq_ion = 0; ///< the frequency ( >= 0 ) of ionic step to output charge density; ///< 0: output only when ion steps are finished diff --git a/source/source_io/module_parameter/read_input_item_output.cpp b/source/source_io/module_parameter/read_input_item_output.cpp index 6274eca87d..93e6ed7be4 100644 --- a/source/source_io/module_parameter/read_input_item_output.cpp +++ b/source/source_io/module_parameter/read_input_item_output.cpp @@ -1,3 +1,4 @@ +#include "source_base/formatter.h" #include "source_base/global_function.h" #include "source_base/tool_quit.h" #include "read_input.h" @@ -14,7 +15,7 @@ void ReadInput::item_output() 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, and Mulliken population analysis is printed every n ionic steps." + 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 = ""; @@ -436,21 +437,67 @@ Also controled by out_freq_ion and out_app_flag. } { Input_Item item("out_stru"); - item.annotation = "output the structure files after each ion step"; + item.annotation = "output the structure files per ion step"; item.category = "Output information"; - item.type = "Boolean"; - item.description = "Whether to output structure files per ionic step in geometry relaxation calculations into OUT.{istep}_D, where ${istep} is the ionic step."; - item.default_value = "False"; + 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.availability = ""; + 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 offlist = {"nscf", "get_s", "get_pchg", "get_wf"}; if (std::find(offlist.begin(), offlist.end(), para.input.calculation) != offlist.end()) { - para.input.out_stru = false; + 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."); } }; - read_sync_bool(input.out_stru); + sync_int(input.out_stru); this->add_item(item); } { diff --git a/source/source_io/test/print_info_test.cpp b/source/source_io/test/print_info_test.cpp index a40d4afa0a..b02ff66acd 100644 --- a/source/source_io/test/print_info_test.cpp +++ b/source/source_io/test/print_info_test.cpp @@ -185,8 +185,8 @@ TEST_F(PrintInfoTest, PrintScreen) ModuleIO::print_screen(stress_step, force_step, istep); output = testing::internal::GetCapturedStdout(); EXPECT_THAT(output,testing::HasSubstr("RELAX STEP")); - EXPECT_THAT(output,testing::HasSubstr("CELL_CHANGE#")); - EXPECT_THAT(output,testing::HasSubstr("IONS_CHANGE#")); + EXPECT_THAT(output,testing::HasSubstr("CELL#")); + EXPECT_THAT(output,testing::HasSubstr("IONS#")); } } } diff --git a/source/source_io/test_serial/read_input_item_test.cpp b/source/source_io/test_serial/read_input_item_test.cpp index 9844e4b4ef..c59e96e484 100644 --- a/source/source_io/test_serial/read_input_item_test.cpp +++ b/source/source_io/test_serial/read_input_item_test.cpp @@ -780,9 +780,9 @@ TEST_F(InputTest, Item_test) { // out_stru auto it = find_label("out_stru", readinput.input_lists); param.input.calculation = "get_wf"; - param.input.out_stru = true; + param.input.out_stru = 1; it->second.reset_value(it->second, param); - EXPECT_EQ(param.input.out_stru, false); + EXPECT_EQ(param.input.out_stru, 0); } { // cal_stress auto it = find_label("cal_stress", readinput.input_lists); @@ -2146,3 +2146,120 @@ TEST_F(InputTest, Item_test_out_mat_vec) EXPECT_EQ(param.input.out_mat_xc2[1], 9); } } + +TEST_F(InputTest, OutStru) +{ + ModuleIO::ReadInput readinput(0); + readinput.check_ntype_flag = false; + Parameter param; + auto it = find_label("out_stru", readinput.input_lists); + ASSERT_NE(it, readinput.input_lists.end()); + + // --- Valid numeric values --- + { + it->second.str_values = {"0"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 0); + } + { + it->second.str_values = {"1"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 1); + } + { + it->second.str_values = {"2"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 2); + } + + // --- Backward-compatible boolean aliases (true -> 1, false -> 0) --- + { + it->second.str_values = {"true"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 1); + } + { + it->second.str_values = {"TRUE"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 1); + } + { + it->second.str_values = {".true."}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 1); + } + { + it->second.str_values = {"Yes"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 1); + } + { + it->second.str_values = {"false"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 0); + } + { + it->second.str_values = {"FALSE"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 0); + } + { + it->second.str_values = {".false."}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 0); + } + { + it->second.str_values = {"No"}; + it->second.read_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 0); + } + + // --- Valid value check_value passes --- + { + for (const int v : {0, 1, 2}) + { + param.input.out_stru = v; + // Expect no exit / no crash; check_value is a void function that only + // calls WARNING_QUIT on bad input. + it->second.check_value(it->second, param); + } + } + + // --- reset_value: calculation in offlist forces out_stru to 0 --- + { + param.input.calculation = "get_wf"; + param.input.out_stru = 1; + it->second.reset_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 0); + + param.input.calculation = "nscf"; + param.input.out_stru = 2; + it->second.reset_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 0); + + // Non-offlist calculation preserves value + param.input.calculation = "cell-relax"; + param.input.out_stru = 1; + it->second.reset_value(it->second, param); + EXPECT_EQ(param.input.out_stru, 1); + } + + // --- Invalid integer values -> WARNING_QUIT via check_value --- + { + for (const std::string& s : {"3", "-1", "-2", "4", "10"}) + { + it->second.str_values = {s}; + it->second.read_value(it->second, param); + EXPECT_EXIT(it->second.check_value(it->second, param), ::testing::ExitedWithCode(1), ""); + } + } + + // --- Non-numeric / malformed inputs -> WARNING_QUIT via read_value --- + { + for (const std::string& s : {"abc", "2.5", "2abc", "-1abc", "xyz", ""}) + { + it->second.str_values = {s}; + EXPECT_EXIT(it->second.read_value(it->second, param), ::testing::ExitedWithCode(1), ""); + } + } +} diff --git a/source/source_md/run_md.cpp b/source/source_md/run_md.cpp index ef28e5c897..b09e64482c 100644 --- a/source/source_md/run_md.cpp +++ b/source/source_md/run_md.cpp @@ -115,10 +115,11 @@ void md_line(UnitCell& unit_in, ModuleESolver::ESolver* p_esolver, const Paramet unitcell::print_stru_file(unit_in, unit_in.atoms, unit_in.latvec, - file.str(), - PARAM.inp.nspin, + file.str(), + "", + PARAM.inp.nspin, false, // Cartesian coordinates - PARAM.inp.calculation == "md", + PARAM.inp.calculation == "md", PARAM.inp.out_mul, need_orb, PARAM.globalv.deepks_setorb, diff --git a/source/source_psi/psi_init_nao.cpp b/source/source_psi/psi_init_nao.cpp index 916d87295a..e514ae063d 100644 --- a/source/source_psi/psi_init_nao.cpp +++ b/source/source_psi/psi_init_nao.cpp @@ -303,7 +303,6 @@ void psi_init_nao::init_psig(T* psig, const int& ik) { /* FOR EACH SPIN CHANNEL */ for (int is_N = 0; is_N < 2; is_N++) // rotate base - // for(int is_N = 0; is_N < 1; is_N++) { if (L == 0 && is_N == 1) { diff --git a/source/source_relax/relax_driver.cpp b/source/source_relax/relax_driver.cpp index 9c85fa2a31..5bbc00c7da 100644 --- a/source/source_relax/relax_driver.cpp +++ b/source/source_relax/relax_driver.cpp @@ -1,5 +1,7 @@ #include "relax_driver.h" +#include "source_base/formatter.h" #include "source_base/global_file.h" +#include "source_base/version.h" #include "source_io/module_output/cif_io.h" #include "source_io/module_json/output_info.h" #include "source_io/module_output/output_log.h" @@ -8,6 +10,8 @@ #include "source_io/module_parameter/parameter.h" #include "source_cell/print_cell.h" +#include + void Relax_Driver::relax_driver( ModuleESolver::ESolver* p_esolver, UnitCell& ucell, @@ -27,14 +31,16 @@ void Relax_Driver::relax_driver( // Main iteration loop for relaxation calculations // For scf/nscf calculations, relax_step returns true immediately, // so the loop exits after one iteration + double etot = 0.0; + ModuleBase::matrix stress(3, 3); + while (steps[0] < inp.relax_nmax) { ModuleBase::matrix force(ucell.nat, 3); - ModuleBase::matrix stress(3, 3); - double etot = 0.0; this->iter_info(steps, inp); this->esolve(steps[0], p_esolver, ucell, inp, force, stress, etot); + this->stru_out(steps[0], ucell, inp, etot, stress); bool converged = this->relax_step(steps, p_esolver, ucell, inp, force, stress, etot, ofs_running); this->json_out(p_esolver, ucell, inp, force, stress); @@ -53,7 +59,7 @@ void Relax_Driver::relax_driver( ++steps[0]; } - this->final_out(steps[0], ucell, inp); + this->final_out(steps[0], ucell, inp, etot, stress); ModuleBase::timer::end("Relax_Driver", "relax_driver"); return; @@ -147,57 +153,108 @@ bool Relax_Driver::relax_step(std::vector& steps, stress, steps[1], steps[2], ofs_running); } - this->stru_out(steps[0], ucell, inp); - ModuleIO::output_after_relax(converged, p_esolver->conv_esolver, ofs_running); return converged; } -void Relax_Driver::stru_out(const int istep, UnitCell& ucell, const Input_para& inp) +void Relax_Driver::stru_out(const int istep, UnitCell& ucell, const Input_para& inp, const double etot, const ModuleBase::matrix& stress) { + // Guard: only output structure files for relaxation calculations + if (inp.calculation != "relax" && inp.calculation != "cell-relax") + { + return; + } + + // out_stru: -1 no output, 0 final only, 1 STRU format, 2 CIF format + // For -1 and 0, no per-step structure output + if (inp.out_stru <= 0) + { + return; + } + + // cache global parameters to reduce repeated PARAM access + const std::string& out_dir = PARAM.globalv.global_out_dir; + const bool deepks_setorb = PARAM.globalv.deepks_setorb; + + // Build header comment with version, timestamp, energy and stress + std::time_t now = std::time(nullptr); + char time_buf[64]; + std::strftime(time_buf, sizeof(time_buf), "%Y-%m-%d %H:%M:%S", std::localtime(&now)); + std::string header = FmtCore::format("# ABACUS version: %s\n# Written at %s\n# RELAX STEP %d, Energy: %.8f eV\n", + VERSION, + time_buf, + istep + 1, + etot * ModuleBase::Ry_to_eV); + // stress in kbar: Ry/Bohr^3 -> kbar, 3 rows + const double stress_transform = ModuleBase::RYDBERG_SI + / (ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI) + * 1.0e-8; + for (int i = 0; i < 3; i++) + { + header += FmtCore::format("# Stress (kbar): %.6f %.6f %.6f\n", + stress(i, 0) * stress_transform, + stress(i, 1) * stress_transform, + stress(i, 2) * stress_transform); + } + bool need_orb = inp.basis_type == "pw"; need_orb = need_orb && inp.init_wfc.substr(0, 3) == "nao"; need_orb = need_orb || inp.basis_type == "lcao"; need_orb = need_orb || inp.basis_type == "lcao_in_pw"; - std::stringstream ss, ss1; - ss << PARAM.globalv.global_out_dir << "STRU_ION_D"; + const bool freq_ok = (inp.out_freq_ion > 0 && istep % inp.out_freq_ion == 0); - unitcell::print_stru_file(ucell, - ucell.atoms, - ucell.latvec, - ss.str(), - inp.nspin, - true, - inp.calculation == "md", - inp.out_mul, - need_orb, - PARAM.globalv.deepks_setorb, - GlobalV::MY_RANK); + // STRU_NOW: overwrite each step (for out_stru 1 and 2) + if (inp.out_stru == 1) + { + unitcell::print_stru_file(ucell, + ucell.atoms, + ucell.latvec, + out_dir + "STRU_NOW", + header, + inp.nspin, + true, + inp.calculation == "md", + inp.out_mul, + need_orb, + deepks_setorb, + GlobalV::MY_RANK); + } + else if (inp.out_stru == 2) + { + ModuleIO::CifParser::write(out_dir + "STRU_NOW.cif", + ucell, + header, + "data_?", + GlobalV::MY_RANK); + } - if (inp.out_stru) + // Numbered files per out_freq_ion (for out_stru 1 and 2 only) + if (freq_ok) { - if (inp.out_freq_ion == 0 || istep % inp.out_freq_ion == 0) + if (inp.out_stru == 1) { - ss1 << PARAM.globalv.global_out_dir << "STRU_ION"; - ss1 << istep+1 << "_D"; unitcell::print_stru_file(ucell, ucell.atoms, ucell.latvec, - ss1.str(), + out_dir + "STRU" + std::to_string(istep + 1), + header, inp.nspin, true, inp.calculation == "md", inp.out_mul, need_orb, - PARAM.globalv.deepks_setorb, + deepks_setorb, GlobalV::MY_RANK); - - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU_NOW.cif", + } + else if (inp.out_stru == 2) + { + ModuleIO::CifParser::write(out_dir + "STRU" + std::to_string(istep + 1) + ".cif", ucell, - "# Generated by ABACUS ModuleIO::CifParser", - "data_?"); + header, + "data_?", + GlobalV::MY_RANK); } } } @@ -213,17 +270,70 @@ void Relax_Driver::json_out(ModuleESolver::ESolver* p_esolver, UnitCell& ucell, #endif } -void Relax_Driver::final_out(const int istep, UnitCell& ucell, const Input_para& inp) +void Relax_Driver::final_out(const int istep, UnitCell& ucell, const Input_para& inp, const double etot, const ModuleBase::matrix& stress) { if (inp.calculation != "relax" && inp.calculation != "cell-relax") { return; } - ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU_FINAL.cif", - ucell, - "# Generated by ABACUS ModuleIO::CifParser", - "data_?"); + // out_stru: 0 no output, 1 STRU format, 2 CIF format + // 1: write STRU_FINAL; 2: write STRU_FINAL.cif + if (inp.out_stru == 1 || inp.out_stru == 2) + { + // cache global parameters to reduce repeated PARAM access + const std::string& out_dir = PARAM.globalv.global_out_dir; + const bool deepks_setorb = PARAM.globalv.deepks_setorb; + + // Build header comment for STRU_FINAL + std::time_t now = std::time(nullptr); + char time_buf[64]; + std::strftime(time_buf, sizeof(time_buf), "%Y-%m-%d %H:%M:%S", std::localtime(&now)); + std::string header = FmtCore::format("# ABACUS version: %s\n# Written at %s\n# RELAX STEP %d (FINAL), Energy: %.8f eV\n", + VERSION, + time_buf, + istep + 1, + etot * ModuleBase::Ry_to_eV); + const double stress_transform = ModuleBase::RYDBERG_SI + / (ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI * ModuleBase::BOHR_RADIUS_SI) + * 1.0e-8; + for (int i = 0; i < 3; i++) + { + header += FmtCore::format("# Stress (kbar): %.6f %.6f %.6f\n", + stress(i, 0) * stress_transform, + stress(i, 1) * stress_transform, + stress(i, 2) * stress_transform); + } + + if (inp.out_stru == 1) + { + bool need_orb = inp.basis_type == "pw"; + need_orb = need_orb && inp.init_wfc.substr(0, 3) == "nao"; + need_orb = need_orb || inp.basis_type == "lcao"; + need_orb = need_orb || inp.basis_type == "lcao_in_pw"; + + unitcell::print_stru_file(ucell, + ucell.atoms, + ucell.latvec, + out_dir + "STRU_FINAL", + header, + inp.nspin, + true, + inp.calculation == "md", + inp.out_mul, + need_orb, + deepks_setorb, + GlobalV::MY_RANK); + } + else if (inp.out_stru == 2) + { + ModuleIO::CifParser::write(out_dir + "STRU_FINAL.cif", + ucell, + header, + "data_?", + GlobalV::MY_RANK); + } + } if (istep == inp.relax_nmax) { diff --git a/source/source_relax/relax_driver.h b/source/source_relax/relax_driver.h index 8348a2d4f7..8636c98fc7 100644 --- a/source/source_relax/relax_driver.h +++ b/source/source_relax/relax_driver.h @@ -94,13 +94,15 @@ class Relax_Driver const double etot, std::ofstream& ofs_running); /** - * @brief Output structure files after relaxation step. + * @brief Output structure files before relaxation move. * * @param istep Current iteration step. * @param ucell Reference to the unit cell. * @param inp Input parameters for the calculation. + * @param etot Total energy in Ry corresponding to this structure. + * @param stress Stress matrix (3x3) in Ry/Bohr^3 corresponding to this structure. */ - void stru_out(const int istep, UnitCell& ucell, const Input_para& inp); + void stru_out(const int istep, UnitCell& ucell, const Input_para& inp, const double etot, const ModuleBase::matrix& stress); /** * @brief Output JSON format results. @@ -120,8 +122,10 @@ class Relax_Driver * @param istep Final iteration step. * @param ucell Reference to the unit cell. * @param inp Input parameters for the calculation. + * @param etot Total energy of the final step. + * @param stress Stress tensor of the final step. */ - void final_out(const int istep, UnitCell& ucell, const Input_para& inp); + void final_out(const int istep, UnitCell& ucell, const Input_para& inp, const double etot, const ModuleBase::matrix& stress); }; #endif diff --git a/source/source_relax/relax_sync.cpp b/source/source_relax/relax_sync.cpp index 212d883354..463f68a578 100644 --- a/source/source_relax/relax_sync.cpp +++ b/source/source_relax/relax_sync.cpp @@ -25,6 +25,7 @@ void Relax::init_relax(const int nat_in) srp_srp = 100000; etot = 0; etot_p = 0; + omega_p = 0.0; force_thr_eva = PARAM.inp.force_thr * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A; // convert to eV/Angstrom fac_force = PARAM.inp.relax_scale_force * 0.1; @@ -66,6 +67,7 @@ bool Relax::relax_step(UnitCell& ucell, if (istep == 0) { etot_p = etot; + omega_p = ucell.omega * pow(ModuleBase::BOHR_TO_A, 3); } bool relax_done = this->setup_gradient(ucell, force, stress, ofs_running); @@ -158,6 +160,16 @@ bool Relax::setup_gradient(const UnitCell& ucell, const ModuleBase::matrix& forc } if (PARAM.inp.out_level == "ie") { + if (if_cell_moves) + { + const double omega_ang = ucell.omega * pow(ModuleBase::BOHR_TO_A, 3); + const double omega_diff = omega_ang - omega_p; + const double omega_ratio = (std::abs(omega_p) > 0.0) ? omega_diff / omega_p * 100.0 : 0.0; + std::cout << " CELL VOLUME (Angstroms^3) : " << omega_ang << std::endl; + std::cout << " VOLUME DIFF (Angstroms^3) : " << omega_diff << std::endl; + std::cout << " VOLUME RATIO (%) : " << omega_ratio << std::endl; + omega_p = omega_ang; + } std::cout << " ETOT DIFF (eV) : " << etot - etot_p << std::endl; std::cout << " LARGEST GRAD (eV/Angstrom) : " << max_grad << std::endl; etot_p = etot; diff --git a/source/source_relax/relax_sync.h b/source/source_relax/relax_sync.h index c022d195b3..9438c6f499 100644 --- a/source/source_relax/relax_sync.h +++ b/source/source_relax/relax_sync.h @@ -92,6 +92,8 @@ class Relax double dmoveh = 0.0; double etot = 0.0; double etot_p = 0.0; + /// previous cell volume in Angstrom^3, used to print volume diff during cell-relax + double omega_p = 0.0; double force_thr_eva = 0.0; bool brent_done = false; // if brent line search is finished