Update output filenames for PCHG

source/module_io/ctrl_output_fp.cpp

@@ -0,0 +1,133 @@
+#include "module_io/ctrl_output_fp.h" // use ctrl_output_fp() 
+
+namespace ModuleIO
+{
+
+template <typename TK, typename TR>
+void ctrl_output_fp(UnitCell& ucell, 
+		elecstate::ElecStateLCAO<TK>* pelec, 
+		const int istep)
+{
+    ModuleBase::TITLE("ModuleIO", "ctrl_output_fp");
+    ModuleBase::timer::tick("ModuleIO", "ctrl_output_fp");
+
+    const bool out_app_flag = PARAM.inp.out_app_flag;
+    const bool gamma_only = PARAM.globalv.gamma_only_local;
+    const int nspin = PARAM.inp.nspin;
+    const std::string global_out_dir = PARAM.globalv.global_out_dir;
+
+	// 1) 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 + "/chgs" + std::to_string(is + 1) + ".cube";
+			ModuleIO::write_vdata_palgrid(Pgrid,
+					this->chr.rho_save[is],
+					is,
+					PARAM.inp.nspin,
+					istep,
+					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 + "/taus" + std::to_string(is + 1) + ".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));
+			}
+		}
+	}
+
+
+	// 2) 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 + "/pots" + std::to_string(is + 1) + ".cube";
+
+			ModuleIO::write_vdata_palgrid(Pgrid,
+					this->pelec->pot->get_effective_v(is),
+					is,
+					PARAM.inp.nspin,
+					istep,
+					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 + "/pot_es.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);
+	}
+
+
+	// 3) 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]);
+	}
+
+    ModuleBase::timer::tick("ModuleIO", "ctrl_output_fp");
+}
+
+} // End ModuleIO
+
+
+// For gamma only
+template void ModuleIO::ctrl_output_lcao<double, double>(UnitCell& ucell, 
+		const int istep);
+
+// For multiple k-points
+template void ModuleIO::ctrl_output_lcao<std::complex<double>, double>(UnitCell& ucell, 
+		const int istep);
+
+template void ModuleIO::ctrl_output_lcao<std::complex<double>, std::complex<double>>(UnitCell& ucell, 
+		const int istep);
+

source/module_io/ctrl_output_fp.h

@@ -0,0 +1,11 @@
+#ifndef CTRL_OUTPUT_FP_H 
+#define CTRL_OUTPUT_FP_H 
+
+namespace ModuleIO
+{
+	template <typename TK, typename TR>
+		void ctrl_output_fp(UnitCell& ucell, 
+				elecstate::ElecStateLCAO<TK>* pelec, 
+				const int istep);
+}
+#endif

source/module_io/get_pchg_lcao.cpp

@@ -55,7 +55,7 @@ void IState_Charge::begin(Gint_Gamma& gg,
 {
     ModuleBase::TITLE("IState_Charge", "begin");
 
-    std::cout << " Calculate |psi(i)|^2 for selected bands (band-decomposed charge densities, gamma only)."
+    std::cout << " Calculate |psi(i)|^2 for selected electronic states (gamma only)."
               << std::endl;
 
     // Determine the mode based on the input parameters
@@ -77,8 +77,8 @@ void IState_Charge::begin(Gint_Gamma& gg,
     // if ucell is even, it's also correct.
     // +1.0e-8 in case like (2.999999999+1)/2
     const int fermi_band = static_cast<int>((nelec + 1) / 2 + 1.0e-8);
-    std::cout << " number of electrons = " << nelec << std::endl;
-    std::cout << " number of occupied bands = " << fermi_band << std::endl;
+    GlobalV::ofs_running << " number of electrons = " << nelec << std::endl;
+    GlobalV::ofs_running << " number of occupied bands = " << fermi_band << std::endl;
 
     // Set this->bands_picked_ according to the mode
     select_bands(nbands_istate, out_pchg, nbands, nelec, mode, fermi_band);
@@ -101,7 +101,7 @@ void IState_Charge::begin(Gint_Gamma& gg,
                 ModuleBase::GlobalFunc::ZEROS(rho[is], rhopw_nrxx);
             }
 
-            std::cout << " Performing grid integral over real space grid for band " << ib + 1 << "..." << std::endl;
+            //std::cout << " Performing grid integral over real space grid for band " << ib + 1 << "..." << std::endl;
 
             DM.init_DMR(GridD_in, ucell_in);
             DM.cal_DMR();
@@ -120,21 +120,22 @@ void IState_Charge::begin(Gint_Gamma& gg,
                 ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is].data(), rhopw_nrxx); // Copy data
             }
 
-            std::cout << " Writing cube files...";
 
             for (int is = 0; is < nspin; ++is)
             {
                 // ssc should be inside the inner loop to reset the string stream each time
                 std::stringstream ssc;
-                ssc << global_out_dir << "BAND" << ib + 1 << "_GAMMA" << "_SPIN" << is + 1 << "_CHG.cube";
+                ssc << global_out_dir << "pchgs" << is + 1 << "i" << ib + 1 << ".cube";
+
+                GlobalV::ofs_running << " Writing cube file " << ssc.str() << std::endl;
 
                 // Use a const vector to store efermi for all spins, replace the original implementation:
                 // const double ef_tmp = pelec->eferm.get_efval(is);
                 double ef_spin = ef_all_spin[is];
                 ModuleIO::write_vdata_palgrid(pgrid, rho_save[is].data(), is, nspin, 0, ssc.str(), ef_spin, ucell_in);
             }
 
-            std::cout << " Complete!" << std::endl;
+            //std::cout << " Complete!" << std::endl;
         }
     }
 
@@ -251,7 +252,7 @@ void IState_Charge::begin(Gint_k& gk,
                     {
                         // ssc should be inside the inner loop to reset the string stream each time
                         std::stringstream ssc;
-                        ssc << global_out_dir << "BAND" << ib + 1 << "_K" << ik + 1 << "_SPIN" << is + 1 << "_CHG.cube";
+                        ssc << global_out_dir << "pchgs" << is + 1 << "k" << ik+1 << "i" << ib + 1 << ".cube";
 
                         double ef_spin = ef_all_spin[is];
                         ModuleIO::write_vdata_palgrid(pgrid,
@@ -292,7 +293,7 @@ void IState_Charge::begin(Gint_k& gk,
                 }
 
                 // Symmetrize the charge density, otherwise the results are incorrect if the symmetry is on
-                std::cout << " Symmetrizing band-decomposed charge density..." << std::endl;
+                // std::cout << " Symmetrizing band-decomposed charge density..." << std::endl;
                 Symmetry_rho srho;
                 for (int is = 0; is < nspin; ++is)
                 {
@@ -310,7 +311,7 @@ void IState_Charge::begin(Gint_k& gk,
                 {
                     // ssc should be inside the inner loop to reset the string stream each time
                     std::stringstream ssc;
-                    ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is + 1 << "_CHG.cube";
+                    ssc << global_out_dir << "pchgs" << is + 1 << "i" << ib + 1 << ".cube";
 
                     double ef_spin = ef_all_spin[is];
                     ModuleIO::write_vdata_palgrid(pgrid,

source/module_io/get_pchg_pw.h

@@ -135,8 +135,8 @@ void get_pchg_pw(const std::vector<int>& out_pchg,
                     }
 
                     std::stringstream ssc;
-                    ssc << global_out_dir << "BAND" << ib + 1 << "_K" << k_number << "_SPIN" << spin_index + 1
-                        << "_CHG.cube";
+                    ssc << global_out_dir << "pchgs" << spin_index + 1 << "k" << k_number << "i" << ib + 1
+                        << ".cube";
 
                     ModuleIO::write_vdata_palgrid(pgrid,
                                                   rho_band[spin_index].data(),
@@ -182,7 +182,7 @@ void get_pchg_pw(const std::vector<int>& out_pchg,
 #endif
 
                 // Symmetrize the charge density, otherwise the results are incorrect if the symmetry is on
-                std::cout << " Symmetrizing band-decomposed charge density..." << std::endl;
+                // std::cout << " Symmetrizing band-decomposed charge density..." << std::endl;
                 Symmetry_rho srho;
                 for (int is = 0; is < nspin; ++is)
                 {
@@ -208,7 +208,7 @@ void get_pchg_pw(const std::vector<int>& out_pchg,
                 for (int is = 0; is < nspin; ++is)
                 {
                     std::stringstream ssc;
-                    ssc << global_out_dir << "BAND" << ib + 1 << "_SPIN" << is + 1 << "_CHG.cube";
+                    ssc << global_out_dir << "pchgs" << is + 1 << "i" << ib + 1 << ".cube";
 
                     ModuleIO::write_vdata_palgrid(pgrid, rho_band[is].data(), is, nspin, 0, ssc.str(), 0.0, ucell);
                 }

tests/01_PW/088_PW_OP_nscf/INPUT → tests/01_PW/054_PW_OP_nscf/INPUT

@@ -5,23 +5,23 @@ calculation     nscf
 
 nbands			8
 symmetry		1
-pseudo_dir	../../PP_ORB
 out_pot			1
 
 #Parameters (2.Iteration)
 ecutwfc			5
-scf_thr				1e-9
-scf_nmax			100
+scf_thr			1e-9
+scf_nmax		100
 
 #Parameters (3.Basis)
 basis_type		pw
 
 #Parameters (4.Smearing)
-smearing_method		gauss
-smearing_sigma			0.002
+smearing_method	gauss
+smearing_sigma	0.002
 
 #Parameters (5.Mixing)
 mixing_type		broyden
 mixing_beta		0.7
 
-read_file_dir ./
+read_file_dir   ./
+pseudo_dir	    ../../PP_ORB

tests/01_PW/085_PW_get_pchg/result.ref

@@ -0,0 +1,10 @@
+etotref -197.1405644417785
+etotperatomref -98.5702822209
+pchgs1i1.cube 2
+pchgs1i2.cube 2
+pchgs1i3.cube 2
+pchgs1i4.cube 2
+pointgroupref T_d
+spacegroupref O_h
+nksibzref 1
+totaltimeref 0.36