Update cell-relaxation output format (#7752)

* fix a bug when init_wfc=nao in pw basis for nspin=4

* update cell-relax output format, including rename the output STRU files

* remove useless words

* update documents

* update STRU output, the final one is named STRU_FINAL

* update out_stru command

* update

* update out_stru options (3 now)

* update fix bug in unitcell_test.cpp

* remove PARAM in relax_driver.cpp

* fix docs

* fix docs

* update docs

---------

Co-authored-by: abacus_fixer <mohanchen@pku.eud.cn>

Author: mohanchen

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
 

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.

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

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.

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.
 

source/source_cell/print_cell.cpp

@@ -80,7 +80,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,
@@ -90,12 +91,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; it<ucell.ntype; it++)
         { 
             str += FmtCore::format("%s %8.4f %s %s\n", 
@@ -119,12 +126,12 @@ namespace unitcell
             str += "\nNUMERICAL_DESCRIPTOR\n" + ucell.descriptor_file + "\n"; 
         }
         // LATTICE_CONSTANT
-        str += "\nLATTICE_CONSTANT\n" + FmtCore::format("%-.10f\n", ucell.lat0);
+        str += "\nLATTICE_CONSTANT\n" + FmtCore::format("%-.10f", ucell.lat0) + "  # in Bohr\n";
         // LATTICE_VECTORS
-        str += "\nLATTICE_VECTORS\n";
-        str += FmtCore::format("%20.10f%20.10f%20.10f\n", latvec.e11, latvec.e12, latvec.e13);
-        str += FmtCore::format("%20.10f%20.10f%20.10f\n", latvec.e21, latvec.e22, latvec.e23);
-        str += FmtCore::format("%20.10f%20.10f%20.10f\n", latvec.e31, latvec.e32, latvec.e33);
+        str += "\nLATTICE_VECTORS  # in units of lat0\n";
+        str += FmtCore::format("%24.16f%24.16f%24.16f\n", latvec.e11, latvec.e12, latvec.e13);
+        str += FmtCore::format("%24.16f%24.16f%24.16f\n", latvec.e21, latvec.e22, latvec.e23);
+        str += FmtCore::format("%24.16f%24.16f%24.16f\n", latvec.e31, latvec.e32, latvec.e33);
         // ATOMIC_POSITIONS
         str += "\nATOMIC_POSITIONS\n";
         const std::string scale = direct? "Direct": "Cartesian";

source/source_cell/print_cell.h

@@ -34,6 +34,7 @@ namespace unitcell
      * @param atoms Atom list [in]
      * @param latvec lattice parameter vector [in]
      * @param fn STRU file name [in]
+     * @param header comment lines written at the top of the file [in]
      * @param nspin number of spin channels [in]
      * @param direct true for direct coords, false for cartesian coords [in]
      * @param vel true for printing velocities [in]
@@ -46,6 +47,7 @@ namespace unitcell
                          const Atom*     atoms,
                          const ModuleBase::Matrix3& latvec,
                          const std::string& fn,
+                         const std::string& header,
                          const int& nspin = 1,
                          const bool& direct = false,
                          const bool& vel = false,

source/source_cell/test/unitcell_test.cpp

@@ -803,7 +803,7 @@ TEST_F(UcellTest, PrintSTRU)
      *
      */
     unitcell::print_stru_file(*ucell,ucell->atoms,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<char>(ifs)), std::istreambuf_iterator<char>());
@@ -813,9 +813,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"));
@@ -835,7 +835,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<char>(ifs)), std::istreambuf_iterator<char>()};
     EXPECT_THAT(str, testing::HasSubstr("ATOMIC_SPECIES"));
@@ -844,9 +844,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"));
@@ -877,7 +877,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<char>(ifs)), std::istreambuf_iterator<char>()};
     EXPECT_THAT(str, testing::HasSubstr("ATOMIC_SPECIES"));
@@ -891,9 +891,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"));

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;

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]};