Reformat running log, delete 'printe' command in INPUT and fix the output eigenvalue information issue (#6326)

* update output files

* keep updating the output information in running log

* update output formats

* change the file name of write_istate_info write_eig_occ

* update output of band information

* update output of eigenvalues, delete useless functions

* adjust output formats in running log

* close screen output of stress

* update documents, delete printe variable and use out_band command

* update examples

* update lbfgs output of coordinates

* update CMakeLists.txt in module_io

* update write_eig_occ_test.cpp

* fix some tests

* fix io tests

* update bessel_basis warning information

* fix extracting the keep information ibz

* fix tests in cell

* fix some output formats

* update some INPUTs

* update outputtest

* update pband

* update DOS info

* fix two tests

* update dH and T(R) matrix output information

* update elecstate_print_test

* update file name of dH/dR

* update examples of dhrxs1

* update test

Author: mohanchen

docs/advanced/input_files/input-main.md

@@ -86,7 +86,6 @@
     - [mixing\_tau](#mixing_tau)
     - [mixing\_dftu](#mixing_dftu)
     - [gamma\_only](#gamma_only)
-    - [printe](#printe)
     - [scf\_nmax](#scf_nmax)
     - [scf\_thr](#scf_thr)
     - [scf\_ene\_thr](#scf_ene_thr)
@@ -1253,12 +1252,6 @@ Note: In new angle mixing, you should set `mixing_beta_mag >> mixing_beta`. The
 
 - **Default**: 0
 
-### printe
-
-- **Type**: Integer
-- **Description**: Print out energy for each band for every printe step
-- **Default**: `scf_nmax`
-
 ### scf_nmax
 
 - **Type**: Integer
@@ -1772,7 +1765,7 @@ These variables are used to control the output of properties.
 ### out_band
 
 - **Type**: Boolean \[Integer\](optional)
-- **Description**: Whether to output the eigenvalues of Hamiltonian matrix (in eV), optionally output precision can be set by a second parameter, default is 8. The output file names are:
+- **Description**: 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: `eigs1.txt`;
     - nspin = 2: `eigs1.txt` and `eigs2.txt`;
     - For more information, refer to the [band.md](../elec_properties/band.md)
@@ -1864,7 +1857,7 @@ These variables are used to control the output of properties.
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis (not gamma-only algorithm)
-- **Description**: Generate files containing the kinetic energy matrix $T(R)$. The format will be the same as the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `trs1.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
+- **Description**: Generate files containing the kinetic energy matrix $T(R)$. The format will be the same as the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `trs1_nao.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
 - **Default**: False
 - **Unit**: Ry
 - **Note**: In the 3.10-LTS version, the file name is data-TR-sparse_SPIN0.csr. 
@@ -1873,7 +1866,7 @@ These variables are used to control the output of properties.
 
 - **Type**: Boolean
 - **Availability**: Numerical atomic orbital basis (not gamma-only algorithm)
-- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `dhrxs1.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
+- **Description**: Whether to print files containing the derivatives of the Hamiltonian matrix. The format will be the same as the Hamiltonian matrix $H(R)$ and overlap matrix $S(R)$ as mentioned in [out_mat_hs2](#out_mat_hs2). The name of the files will be `dhrxs1_nao.csr`, `dhrys1_nao.csr`, `dhrzs1_nao.csr` and so on. Also controled by [out_interval](#out_interval) and [out_app_flag](#out_app_flag).
 - **Default**: False
 - **Unit**: Ry/Bohr
 - **Note**: In the 3.10-LTS version, the file name is data-dHRx-sparse_SPIN0.csr and so on.

examples/relax/lcao_output/INPUT

@@ -3,9 +3,9 @@ INPUT_PARAMETERS
 suffix          autotest
 
 nbands			8
-calculation     cell-relax
+calculation     relax #cell-relax
 ecutwfc			10
-scf_nmax		20
+scf_nmax		2
 
 basis_type		lcao
 relax_nmax      5
@@ -20,6 +20,7 @@ mixing_type     broyden
 mixing_beta     0.7
 
 relax_new       0
+relax_method    bfgs
 
 pseudo_dir      ../../../tests/PP_ORB	
 orbital_dir	    ../../../tests/PP_ORB

examples/relax/lcao_output/STRU

@@ -15,7 +15,7 @@ LATTICE_VECTORS
 ATOMIC_POSITIONS
 Cartesian 		#Cartesian(Unit is LATTICE_CONSTANT)
 Si 			#Name of element	
-0.0			#Magnetic for this element.
+4.0			#Magnetic for this element.
 2			#Number of atoms
-0.00 0.00 0.00 0 0 0	#x,y,z, move_x, move_y, move_z
-0.25 0.25 0.25 0 0 0
+0.00 0.00 0.00 1 1 1	#x,y,z, move_x, move_y, move_z
+0.27 0.25 0.25 1 1 1

examples/relax/pw_output/INPUT

@@ -3,7 +3,7 @@ INPUT_PARAMETERS
 suffix          autotest
 
 nbands			8
-calculation     cell-relax
+calculation     relax
 ecutwfc			10
 scf_nmax		20
 
@@ -20,6 +20,7 @@ mixing_type     broyden
 mixing_beta     0.7
 
 relax_new       0
+relax_method    lbfgs 
 
 pseudo_dir      ../../../tests/PP_ORB	
 orbital_dir	    ../../../tests/PP_ORB
@@ -30,8 +31,12 @@ out_chg         1
 out_pot         1
 out_wfc_pw      1
 out_dos         1
-out_band        1
 out_stru        1
 out_app_flag    0
 
 out_interval    1
+kpar            2
+symmetry        -1
+
+#out_freq_elec   2
+#out_band        1

examples/relax/pw_output/STRU

@@ -18,4 +18,4 @@ Si 			#Name of element
 0.0			#Magnetic for this element.
 2			#Number of atoms
 0.00 0.00 0.00 0 0 0	#x,y,z, move_x, move_y, move_z
-0.25 0.25 0.25 0 0 0
+0.27 0.243 0.233 1 1 1

source/Makefile.Objects

@@ -571,7 +571,7 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_orb_info.o\
       write_dos_lcao.o\
       write_proj_band_lcao.o\
-      write_istate_info.o\
+      write_eig_occ.o\
       get_pchg_lcao.o\
       get_wf_lcao.o\
       io_dmk.o\

source/module_elecstate/cal_nelec_nband.cpp

@@ -7,22 +7,22 @@ namespace elecstate {
 void cal_nelec(const Atom* atoms, const int& ntype, double& nelec)
 {
     ModuleBase::TITLE("UnitCell", "cal_nelec");
-    GlobalV::ofs_running << "\n SETUP THE ELECTRONS NUMBER" << std::endl;
+    //GlobalV::ofs_running << "\n Setup number of electrons" << std::endl;
 
     if (nelec == 0)
     {
         for (int it = 0; it < ntype; it++)
         {
             std::stringstream ss1, ss2;
-            ss1 << "electron number of element " << atoms[it].label;
+            ss1 << "Electron number of element " << atoms[it].label;
             const double nelec_it = atoms[it].ncpp.zv * atoms[it].na;
             nelec += nelec_it;
-            ss2 << "total electron number of element " << atoms[it].label;
+            ss2 << "Total electron number of element " << atoms[it].label;
 
             ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss1.str(), atoms[it].ncpp.zv);
             ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, ss2.str(), nelec_it);
         }
-        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "AUTOSET number of electrons: ", nelec);
+        ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Autoset the number of electrons", nelec);
     }
     if (PARAM.inp.nelec_delta != 0)
     {
@@ -42,20 +42,22 @@ void cal_nbands(const int& nelec, const int& nlocal, const std::vector<double>&
     {
         return;
     }
+
     //=======================================
     // calculate number of bands (setup.f90)
     //=======================================
     double occupied_bands = static_cast<double>(nelec / ModuleBase::DEGSPIN);
-    if (PARAM.inp.lspinorb == 1) {
-        occupied_bands = static_cast<double>(nelec);
-    }
+	if (PARAM.inp.lspinorb == 1) 
+	{
+		occupied_bands = static_cast<double>(nelec);
+	}
 
     if ((occupied_bands - std::floor(occupied_bands)) > 0.0)
     {
         occupied_bands = std::floor(occupied_bands) + 1.0; // mohan fix 2012-04-16
     }
 
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "occupied bands", occupied_bands);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Occupied electronic states", occupied_bands);
 
     if (nbands == 0)
     {
@@ -89,8 +91,6 @@ void cal_nbands(const int& nelec, const int& nlocal, const std::vector<double>&
         }
         ModuleBase::GlobalFunc::AUTO_SET("NBANDS", nbands);
     }
-    // else if ( PARAM.inp.calculation=="scf" || PARAM.inp.calculation=="md" || PARAM.inp.calculation=="relax") //pengfei
-    // 2014-10-13
     else
     {
         if (nbands < occupied_bands) {
@@ -112,7 +112,6 @@ void cal_nbands(const int& nelec, const int& nlocal, const std::vector<double>&
     }
 
     // mohan add 2010-09-04
-    // std::cout << "nbands(this-> = " <<nbands <<std::endl;
     if (nbands == occupied_bands)
     {
         if (PARAM.inp.smearing_method != "fixed")
@@ -127,16 +126,16 @@ void cal_nbands(const int& nelec, const int& nlocal, const std::vector<double>&
     {
         if (nbands > nlocal)
         {
-            ModuleBase::WARNING_QUIT("ElecState::cal_nbandsc", "NLOCAL < NBANDS");
+            ModuleBase::WARNING_QUIT("ElecState::cal_nbands", 
+            "Number of basis (NLOCAL) < Number of electronic states (NBANDS)");
         }
         else
         {
-            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "NLOCAL", nlocal);
-            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "NBANDS", nbands);
+            ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Number of basis (NLOCAL)", nlocal);
         }
     }
 
-    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "NBANDS", nbands);
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Number of electronic states (NBANDS)", nbands);
 }
 
-}
+}

source/module_elecstate/cal_wfc.cpp

@@ -41,8 +41,8 @@ namespace elecstate
             }
         }
 
-        log << " " << std::setw(40) << "NLOCAL"
-            << " = " << nlocal_tmp << std::endl;
+//        log << " " << std::setw(40) << "NLOCAL"
+  //          << " = " << nlocal_tmp << std::endl;
         //========================================================
         // (4) set index for itia2iat, itiaiw2iwt
         //========================================================
@@ -140,33 +140,38 @@ namespace elecstate
     void cal_natomwfc(std::ofstream& log,int& natomwfc,const int ntype,const Atom* atoms) 
     {
         natomwfc = 0;
-        for (int it = 0; it < ntype; it++) {
-            //============================
-            // Use pseudo-atomic orbitals
-            //============================
-            int tmp = 0;
-            for (int l = 0; l < atoms[it].ncpp.nchi; l++) {
-                if (atoms[it].ncpp.oc[l] >= 0) {
-                    if (PARAM.inp.nspin == 4) {
-                        if (atoms[it].ncpp.has_so) {
-                            tmp += 2 * atoms[it].ncpp.lchi[l];
-                            if (fabs(atoms[it].ncpp.jchi[l] - atoms[it].ncpp.lchi[l] - 0.5)< 1e-6) 
-                                {
-                                    tmp += 2;
-                                }
-                        } else {
-                            tmp += 2 * (2 * atoms[it].ncpp.lchi[l] + 1);
-                        }
-                    } else {
-                        tmp += 2 * atoms[it].ncpp.lchi[l] + 1;
-                    }
-                }
-            }
-            natomwfc += tmp * atoms[it].na;
-        }
-        ModuleBase::GlobalFunc::OUT(log,
-                                    "initial pseudo atomic orbital number",
-                                    natomwfc);
-        return;
+		for (int it = 0; it < ntype; it++) 
+		{
+			//============================
+			// Use pseudo-atomic orbitals
+			//============================
+			int tmp = 0;
+			for (int l = 0; l < atoms[it].ncpp.nchi; l++) 
+			{
+				if (atoms[it].ncpp.oc[l] >= 0) 
+				{
+					if (PARAM.inp.nspin == 4) 
+					{
+						if (atoms[it].ncpp.has_so) 
+						{
+							tmp += 2 * atoms[it].ncpp.lchi[l];
+							if (fabs(atoms[it].ncpp.jchi[l] - atoms[it].ncpp.lchi[l] - 0.5)< 1e-6) 
+							{
+								tmp += 2;
+							}
+						} else 
+						{
+							tmp += 2 * (2 * atoms[it].ncpp.lchi[l] + 1);
+						}
+					} else 
+					{
+						tmp += 2 * atoms[it].ncpp.lchi[l] + 1;
+					}
+				}
+			}
+			natomwfc += tmp * atoms[it].na;
+		}
+		ModuleBase::GlobalFunc::OUT(log, "Number of pseudo atomic orbitals", natomwfc);
+		return;
     }
-}
+}

source/module_elecstate/elecstate_print.cpp

@@ -149,56 +149,6 @@ void print_scf_iterinfo(const std::string& ks_solver,
     std::cout << buf;
 }
 
-/// @brief function for printing eigenvalues : ekb
-/// @param ik: index of kpoints
-/// @param printe: print energy every 'printe' electron iteration.
-/// @param iter: index of iterations
-void print_band(const ModuleBase::matrix& ekb,
-                const ModuleBase::matrix& wg,
-                const K_Vectors* klist,
-                const int& ik,
-                const int& printe,
-                const int& iter,
-                std::ofstream &ofs)
-{
-    const double largest_eig = 1.0e10;
-
-    // check the band energy.
-    bool wrong = false;
-    for (int ib = 0; ib < PARAM.globalv.nbands_l; ++ib)
-    {
-        if (std::abs(ekb(ik, ib)) > largest_eig)
-        {
-            GlobalV::ofs_warning << " ik=" << ik + 1 << " ib=" << ib + 1 << " " << ekb(ik, ib) << " Ry" << std::endl;
-            wrong = true;
-        }
-    }
-    if (wrong)
-    {
-        ModuleBase::WARNING_QUIT("print_eigenvalue", "Eigenvalues are too large!");
-    }
-
-    if (GlobalV::MY_RANK == 0)
-    {
-        if (printe > 0 && ((iter + 1) % printe == 0))
-        {
-            ofs << std::setprecision(6);
-            ofs << " Energy (eV) & Occupations for spin=" << klist->isk[ik] + 1
-                                 << " k-point=" << ik + 1 << std::endl;
-            ofs << std::setiosflags(std::ios::showpoint);
-            for (int ib = 0; ib < PARAM.globalv.nbands_l; ib++)
-            {
-                ofs << " " << std::setw(6) << ib + 1 << std::setw(15)
-                                     << ekb(ik, ib) * ModuleBase::Ry_to_eV;
-                // for the first electron iteration, we don't have the energy
-                // spectrum, so we can't get the occupations.
-                ofs << std::setw(15) << wg(ik, ib);
-                ofs << std::endl;
-            }
-        }
-    }
-    return;
-}
 
 /// @brief print total free energy and other energies
 /// @param ucell: unit cell
@@ -216,7 +166,6 @@ void print_etot(const Magnetism& magnet,
                 const double& scf_thr,
                 const double& scf_thr_kin,
                 const double& duration,
-                const int printe,
                 const double& pw_diag_thr,
                 const double& avg_iter,
                 const bool print)
@@ -228,7 +177,8 @@ void print_etot(const Magnetism& magnet,
 
     GlobalV::ofs_running << std::setprecision(12);
     GlobalV::ofs_running << std::setiosflags(std::ios::right);
-    GlobalV::ofs_running << " Electron density deviation is " << scf_thr << std::endl;
+
+    ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"Electron density deviation",scf_thr);
 
     if (PARAM.inp.basis_type == "pw")
     {
@@ -239,8 +189,8 @@ void print_etot(const Magnetism& magnet,
     std::vector<double> energies_Ry;
     std::vector<double> energies_eV;
 
-    if (printe > 0 && ((iter + 1) % printe == 0 || converged || iter == PARAM.inp.scf_nmax))
-    {
+	if( (iter % PARAM.inp.out_freq_elec == 0) || converged || iter == PARAM.inp.scf_nmax )
+	{
         int n_order = std::max(0, Occupy::gaussian_type);
         titles.push_back("E_KohnSham");
         energies_Ry.push_back(elec.f_en.etot);
@@ -350,7 +300,7 @@ void print_etot(const Magnetism& magnet,
     FmtTable table(/*titles=*/{"Energy", "Rydberg", "eV"},
                    /*nrows=*/titles.size(),
                    /*formats=*/{"%-14s", "%20.10f", "%20.10f"}, 
-                   /*indents=*/0,
+                   /*indents=*/1,
                    /*align=*/{/*value*/FmtTable::Align::LEFT, /*title*/FmtTable::Align::CENTER});
     // print out the titles
     table << titles << energies_Ry << energies_eV;