Merge branch 'develop' into develop

Author: lyb9812

docs/advanced/input_files/input-main.md

@@ -732,7 +732,7 @@ If only one value is set (such as `kspacing 0.5`), then kspacing values of a/b/c
 
 [back to top](#full-list-of-input-keywords)
 
-## Variables related to input files
+## Input files
 
 These variables are used to control parameters related to input files.
 
@@ -1656,16 +1656,16 @@ These variables are used to control the geometry relaxation.
 
 [back to top](#full-list-of-input-keywords)
 
-## Variables related to output information
+## Output information
 
 These variables are used to control the output of properties.
 
 ### out_freq_ion
 
 - **Type**: Integer
-- **Description**: After self-consistent-field calculations, control the interval of ionic movements for printing properties. These properties cover charge density, local potential, electrostatic potential, Hamiltonian matrix, overlap matrix, density matrix, Mulliken population analysis and so on.
+- **Description**: Control the interval to print information every few ion steps. These properties cover charge density, local potential, electrostatic potential, Hamiltonian matrix, overlap matrix, density matrix, Mulliken population analysis and so on.
 - **Default**: 0
-- **Note**: If you want to use out_freq_elec, please set out_freq_ion to 1, otherwise out_freq_elec is useless
+- **Note**: The integer indicates to print information every 'out_freq_ion' ion steps. 
 
 ### out_freq_elec
 
@@ -1679,11 +1679,11 @@ These variables are used to control the output of properties.
 - **Description**:
   The first integer controls whether to output the charge density on real space grids:
   - 1: Output the charge density (in Bohr^-3) on real space grids into the density files in the folder `OUT.${suffix}`. The files are named as:
-    - nspin = 1: `chgs1.cube`;
+    - nspin = 1: `chg.cube`;
     - nspin = 2: `chgs1.cube`, and `chgs2.cube`;
     - nspin = 4: `chgs1.cube`, `chgs2.cube`, `chgs3.cube`, and `chgs4.cube`;
   Note that by using the Meta-GGA functional, additional files containing the kinetic energy density will be output with the following names:
-    - nspin = 1: `taus1.cube`;
+    - nspin = 1: `tau.cube`;
     - nspin = 2: `taus1.cube`, and `taus2.cube`;
     - nspin = 4: `taus1.cube`, `taus2.cube`, `taus3.cube`, and `taus4.cube`;
   - 2: On top of 1, also output the initial charge density files with a suffix name as '_ini', such as `taus1_ini.cube`, etc.

source/Makefile.Objects

@@ -580,6 +580,7 @@ OBJS_IO=input_conv.o\
     output_log.o\
     output_mat_sparse.o\
     ctrl_output_lcao.o\
+    ctrl_output_fp.o\
     para_json.o\
     abacusjson.o\
     general_info.o\

source/source_esolver/esolver_fp.cpp

@@ -8,20 +8,15 @@
 #include "source_hamilt/module_vdw/vdw.h"
 #include "source_pw/module_pwdft/global.h"
 #include "source_io/cif_io.h"
-#include "source_io/cube_io.h"
+#include "source_io/cube_io.h" // use write_vdata_palgrid
 #include "source_io/json_output/init_info.h"
 #include "source_io/json_output/output_info.h"
 #include "source_io/output_log.h"
 #include "source_io/print_info.h"
 #include "source_io/rhog_io.h"
-#include "source_io/write_elecstat_pot.h"
-#include "source_io/write_elf.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_cell/k_vector_utils.h"
-
-#ifdef USE_LIBXC
-#include "source_io/write_libxc_r.h"
-#endif
+#include "source_io/ctrl_output_fp.h"
 
 namespace ModuleESolver
 {
@@ -162,169 +157,10 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
     // 3) update delta_rho for charge extrapolation
     CE.update_delta_rho(ucell, &(this->chr), &(this->sf));
 
+    // 4) print out charge density, potential, elf, etc.
+	ModuleIO::ctrl_output_fp(ucell, this->pelec, this->pw_big, this->pw_rhod, 
+			this->chr, this->solvent, this->Pgrid, istep); 
 
-    // print out the 'g' index when istep_in!=-1
-    int istep_in = -1;
-	if (PARAM.inp.out_freq_ion>0) // default value of out_freq_ion is 0
-	{
-		if (istep % PARAM.inp.out_freq_ion == 0)
-		{
-			istep_in = istep;
-		}
-	}
-
-    std::string geom_block;
-    if(istep_in==-1)
-	{
-		// do nothing
-	}
-    else if(istep_in>=0)
-	{
-		geom_block = "g" + std::to_string(istep + 1);
-	}
-
-	// 4) write charge density
-	if (PARAM.inp.out_chg[0] > 0)
-	{
-		for (int is = 0; is < PARAM.inp.nspin; ++is)
-		{
-			this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
-
-			std::string fn =PARAM.globalv.global_out_dir + "chg";
-
-            std::string spin_block;
-			if(PARAM.inp.nspin == 2 || PARAM.inp.nspin == 4)
-			{
-				spin_block= "s" + std::to_string(is + 1);
-			}
-			else if(PARAM.inp.nspin == 1)
-			{
-				// do nothing
-			}
-
-			fn += spin_block + geom_block + ".cube";
-
-			ModuleIO::write_vdata_palgrid(Pgrid,
-					this->chr.rho_save[is],
-					is,
-					PARAM.inp.nspin,
-					istep_in,
-					fn,
-					this->pelec->eferm.get_efval(is),
-					&(ucell),
-					PARAM.inp.out_chg[1],
-					1);
-
-			if (XC_Functional::get_ked_flag())
-			{
-				fn = PARAM.globalv.global_out_dir + "tau";
-
-                fn += spin_block + geom_block + ".cube";
-
-				ModuleIO::write_vdata_palgrid(Pgrid,
-						this->chr.kin_r_save[is],
-						is,
-						PARAM.inp.nspin,
-						istep,
-						fn,
-						this->pelec->eferm.get_efval(is),
-						&(ucell));
-			}
-		}
-	}
-
-	// 5) write potential
-	if (PARAM.inp.out_pot == 1 || PARAM.inp.out_pot == 3)
-	{
-		for (int is = 0; is < PARAM.inp.nspin; is++)
-		{
-			std::string fn =PARAM.globalv.global_out_dir + "pot";
-
-            std::string spin_block;
-			if(PARAM.inp.nspin == 2 || PARAM.inp.nspin == 4)
-			{
-				spin_block= "s" + std::to_string(is + 1);
-			}
-			else if(PARAM.inp.nspin == 1)
-			{
-				// do nothing
-			}
-
-			fn += spin_block + geom_block + ".cube";
-
-			ModuleIO::write_vdata_palgrid(Pgrid,
-					this->pelec->pot->get_effective_v(is),
-					is,
-					PARAM.inp.nspin,
-					istep_in,
-					fn,
-					0.0, // efermi
-					&(ucell),
-					3,  // precision
-					0); // out_fermi
-		}
-	}
-	else if (PARAM.inp.out_pot == 2)
-	{
-		std::string fn =PARAM.globalv.global_out_dir + "potes";
-        fn += geom_block + ".cube";
-
-		ModuleIO::write_elecstat_pot(
-#ifdef __MPI
-				this->pw_big->bz,
-				this->pw_big->nbz,
-#endif
-				fn,
-				istep,
-				this->pw_rhod,
-				&this->chr,
-				&(ucell),
-				this->pelec->pot->get_fixed_v(),
-				this->solvent);
-	}
-
-	// 6) write ELF
-	if (PARAM.inp.out_elf[0] > 0)
-	{
-		this->chr.cal_elf = true;
-		Symmetry_rho srho;
-		for (int is = 0; is < PARAM.inp.nspin; is++)
-		{
-			srho.begin(is, this->chr, this->pw_rhod, ucell.symm);
-		}
-
-		std::string out_dir =PARAM.globalv.global_out_dir;
-		ModuleIO::write_elf(
-#ifdef __MPI
-				this->pw_big->bz,
-				this->pw_big->nbz,
-#endif
-				out_dir,
-				istep,
-				PARAM.inp.nspin,
-				this->chr.rho,
-				this->chr.kin_r,
-				this->pw_rhod,
-				this->Pgrid,
-				&(ucell),
-				PARAM.inp.out_elf[1]);
-	}
-
-#ifdef USE_LIBXC
-	// 7) write xc(r)
-	if(PARAM.inp.out_xc_r[0]>=0)
-	{
-		ModuleIO::write_libxc_r(
-				PARAM.inp.out_xc_r[0],
-				XC_Functional::get_func_id(),
-				this->pw_rhod->nrxx, // number of real-space grid
-				ucell.omega, // volume of cell
-				ucell.tpiba,
-				this->chr,
-				*this->pw_big,
-				*this->pw_rhod);
-	}
-#endif
 }
 
 void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
@@ -402,19 +238,32 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     //----------------------------------------------------------
     elecstate::cal_ux(ucell);
 
+
+
     //----------------------------------------------------------
     //! output the initial charge density
     //----------------------------------------------------------
+    const int nspin = PARAM.inp.nspin;
     if (PARAM.inp.out_chg[0] == 2)
     {
-        for (int is = 0; is < PARAM.inp.nspin; is++)
+        for (int is = 0; is < nspin; is++)
         {
             std::stringstream ss;
-            ss << PARAM.globalv.global_out_dir << "chgs" << is + 1 << "_ini.cube";
+            ss << PARAM.globalv.global_out_dir << "chg";
+
+			if(nspin==1)
+			{
+				ss << "ini.cube";
+			}
+			else if(nspin==2 || nspin==4)
+			{
+				ss << "s" << is + 1 << "ini.cube";
+			}
+
             ModuleIO::write_vdata_palgrid(this->Pgrid,
                                           this->chr.rho[is],
                                           is,
-                                          PARAM.inp.nspin,
+                                          nspin,
                                           istep,
                                           ss.str(),
                                           this->pelec->eferm.ef,
@@ -427,14 +276,24 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     //----------------------------------------------------------
     if (PARAM.inp.out_pot == 3)
     {
-        for (int is = 0; is < PARAM.inp.nspin; is++)
+        for (int is = 0; is < nspin; is++)
         {
             std::stringstream ss;
-            ss << PARAM.globalv.global_out_dir << "pots" << is + 1 << "_ini.cube";
+            ss << PARAM.globalv.global_out_dir << "pot";
+
+			if(nspin==1)
+			{
+				ss << "ini.cube";
+			}
+			else if(nspin==2 || nspin==4)
+			{
+				ss << "s" << is + 1 << "ini.cube";
+			}
+
             ModuleIO::write_vdata_palgrid(this->Pgrid,
                                           this->pelec->pot->get_effective_v(is),
                                           is,
-                                          PARAM.inp.nspin,
+                                          nspin,
                                           istep,
                                           ss.str(),
                                           0.0, // efermi

source/source_esolver/esolver_fp.h

@@ -7,26 +7,13 @@
 #include <chrono>
 #endif
 
-//! plane wave basis
-#include "source_basis/module_pw/pw_basis.h"
-
-//! symmetry analysis
-#include "source_cell/module_symmetry/symmetry.h"
-
-//! electronic states
-#include "source_estate/elecstate.h"
-
-//! charge extrapolation
-#include "source_estate/module_charge/charge_extra.h"
-
-//! solvation model
-#include "source_hamilt/module_surchem/surchem.h"
-
-//! local pseudopotential
-#include "source_pw/module_pwdft/VL_in_pw.h"
-
-//! structure factor related to plane wave basis
-#include "source_pw/module_pwdft/structure_factor.h"
+#include "source_basis/module_pw/pw_basis.h" // plane wave basis
+#include "source_cell/module_symmetry/symmetry.h" // symmetry analysis
+#include "source_estate/elecstate.h" // electronic states
+#include "source_estate/module_charge/charge_extra.h" // charge extrapolation
+#include "source_hamilt/module_surchem/surchem.h" // solvation model
+#include "source_pw/module_pwdft/VL_in_pw.h" // local pseudopotential
+#include "source_pw/module_pwdft/structure_factor.h" // structure factor
 
 #include <fstream>
 
@@ -44,10 +31,8 @@ namespace ModuleESolver
 class ESolver_FP: public ESolver
 {
   public:
-    //! Constructor
     ESolver_FP();
 
-    //! Deconstructor
     virtual ~ESolver_FP();
 
     //! Initialize of the first-principels energy solver
@@ -56,13 +41,10 @@ class ESolver_FP: public ESolver
     virtual void after_all_runners(UnitCell& ucell) override;
 
   protected:
-    //! Something to do before SCF iterations.
     virtual void before_scf(UnitCell& ucell, const int istep);
 
-    //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
     virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver);
 
-    //! Something to do after hamilt2rho function in each iter loop.
     virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver);
 
     //! ------------------------------------------------------------------------------

source/source_esolver/esolver_of.cpp

@@ -11,9 +11,7 @@
 #include "source_pw/module_pwdft/global.h"
 #include "source_io/print_info.h"
 #include "source_estate/cal_ux.h"
-//-----force-------------------
 #include "source_pw/module_pwdft/forces.h"
-//-----stress------------------
 #include "source_pw/module_ofdft/of_stress_pw.h"
 
 namespace ModuleESolver