Feature: support segment split in kline mode in KPT file and out_band band output precision control, 8 as default

docs/advanced/input_files/input-main.md

@@ -1494,8 +1494,8 @@ These variables are used to control the output of properties.
 
 ### out_band
 
-- **Type**: Boolean
-- **Description**: Whether to output the band structure (in eV). For more information, refer to the [band.md](../elec_properties/band.md)
+- **Type**: Boolean Integer(optional)
+- **Description**: Whether to output the band structure (in eV), optionally output precision can be set by a second parameter, default is 8. For more information, refer to the [band.md](../elec_properties/band.md)
 - **Default**: False
 
 ### out_proj_band

source/module_cell/klist.cpp

@@ -362,6 +362,10 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 
 			//recalculate nkstot.
 			nkstot = 0;
+            /* ISSUE#3482: to distinguish different kline segments */
+            std::vector<int> kpt_segids;
+            kl_segids.clear(); kl_segids.shrink_to_fit();
+            int kpt_segid = 0;
 			for(int iks=0; iks<nks_special; iks++)
 			{
 				ifk >> ksx[iks];
@@ -371,6 +375,9 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 				//std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl;
 				assert(nkl[iks] >= 0);
 				nkstot += nkl[iks];
+                /* ISSUE#3482: to distinguish different kline segments */
+                if((nkl[iks] == 1)&&(iks!=(nks_special-1))) kpt_segid++;
+                kpt_segids.push_back(kpt_segid);
 			}
 			assert( nkl[nks_special-1] == 1);
 
@@ -389,6 +396,7 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 					kvec_c[count].x = ksx[iks-1] + is*dx;
 					kvec_c[count].y = ksy[iks-1] + is*dy;
 					kvec_c[count].z = ksz[iks-1] + is*dz;
+                    kl_segids.push_back(kpt_segids[iks-1]); /* ISSUE#3482: to distinguish different kline segments */
 					++count;
 				}
 			}
@@ -397,15 +405,14 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 			kvec_c[count].x = ksx[nks_special-1];
 			kvec_c[count].y = ksy[nks_special-1];
 			kvec_c[count].z = ksz[nks_special-1];
+            kl_segids.push_back(kpt_segids[nks_special-1]); /* ISSUE#3482: to distinguish different kline segments */
 			++count;
 
 			//std::cout << " count = " << count << std::endl;
-			assert (count == nkstot );
-
-			for(int ik=0; ik<nkstot; ik++)
-			{
-				wk[ik] = 1.0;
-			}
+			assert(count == nkstot);
+            assert(kl_segids.size() == nkstot); /* ISSUE#3482: to distinguish different kline segments */
+
+            std::for_each(wk.begin(), wk.end(), [](double& d){d = 1.0;});
 
             this->kc_done = true;
 
@@ -439,15 +446,22 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 
 			//recalculate nkstot.
 			nkstot = 0;
+            /* ISSUE#3482: to distinguish different kline segments */
+            std::vector<int> kpt_segids;
+            kl_segids.clear(); kl_segids.shrink_to_fit();
+            int kpt_segid = 0;
 			for(int iks=0; iks<nks_special; iks++)
 			{
 				ifk >> ksx[iks];
 				ifk >> ksy[iks];
 				ifk >> ksz[iks];
-				ModuleBase::GlobalFunc::READ_VALUE( ifk, nkl[iks] );
+				ModuleBase::GlobalFunc::READ_VALUE( ifk, nkl[iks] ); /* so ifk is ifstream for kpoint, then nkl is number of kpoints on line */
 				//std::cout << " nkl[" << iks << "]=" << nkl[iks] << std::endl;
 				assert(nkl[iks] >= 0);
 				nkstot += nkl[iks];
+                /* ISSUE#3482: to distinguish different kline segments */
+                if((nkl[iks] == 1)&&(iks!=(nks_special-1))) kpt_segid++;
+                kpt_segids.push_back(kpt_segid);
 			}
 			assert( nkl[nks_special-1] == 1);
 
@@ -466,6 +480,7 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 					kvec_d[count].x = ksx[iks-1] + is*dx;
 					kvec_d[count].y = ksy[iks-1] + is*dy;
 					kvec_d[count].z = ksz[iks-1] + is*dz;
+                    kl_segids.push_back(kpt_segids[iks-1]); /* ISSUE#3482: to distinguish different kline segments */
 					++count;
 				}
 			}
@@ -474,18 +489,16 @@ bool K_Vectors::read_kpoints(const std::string &fn)
 			kvec_d[count].x = ksx[nks_special-1];
 			kvec_d[count].y = ksy[nks_special-1];
 			kvec_d[count].z = ksz[nks_special-1];
+            kl_segids.push_back(kpt_segids[nks_special-1]); /* ISSUE#3482: to distinguish different kline segments */
 			++count;
 
 			//std::cout << " count = " << count << std::endl;
-			assert (count == nkstot );
+			assert(count == nkstot );
+            assert(kl_segids.size() == nkstot); /* ISSUE#3482: to distinguish different kline segments */
 
-			for(int ik=0; ik<nkstot; ik++)
-			{
-				wk[ik] = 1.0;
-			}
+			std::for_each(wk.begin(), wk.end(), [](double& d){d = 1.0;});
 
             this->kd_done = true;
-
 		}
 
         else
@@ -1122,6 +1135,9 @@ void K_Vectors::mpi_k(void)
 
     Parallel_Common::bcast_int(nmp, 3);
 
+    kl_segids.resize(nkstot);
+    Parallel_Common::bcast_int(kl_segids.data(), nkstot);
+
     Parallel_Common::bcast_double(koffset, 3);
 
     this->nks = GlobalC::Pkpoints.nks_pool[GlobalV::MY_POOL];
@@ -1352,6 +1368,8 @@ void K_Vectors::mpi_k_after_vc(void)
     Parallel_Common::bcast_int(nspin);
     Parallel_Common::bcast_int(nkstot);
     Parallel_Common::bcast_int(nmp, 3);
+    kl_segids.resize(nkstot);
+    Parallel_Common::bcast_int(kl_segids.data(), nkstot);
     Parallel_Common::bcast_double(koffset, 3);
 
     this->nks = GlobalC::Pkpoints.nks_pool[GlobalV::MY_POOL];

source/module_cell/klist.h

@@ -29,6 +29,7 @@ class K_Vectors
     int nkstot_full;    /// number of k points in full k mesh
 
     int nmp[3];						// Number of Monhorst-Pack
+    std::vector<int> kl_segids;	// index of kline segment
 
     K_Vectors();
     ~K_Vectors();

source/module_esolver/esolver_ks_lcao.cpp

@@ -309,7 +309,7 @@ namespace ModuleESolver
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
     GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
 
-    if (INPUT.out_dos != 0 || INPUT.out_band != 0 || INPUT.out_proj_band != 0)
+    if (INPUT.out_dos != 0 || INPUT.out_band[0] != 0 || INPUT.out_proj_band != 0)
     {
         GlobalV::ofs_running << "\n\n\n\n";
         GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
@@ -331,7 +331,7 @@ namespace ModuleESolver
 
     int nspin0 = (GlobalV::NSPIN == 2) ? 2 : 1;
 
-    if (INPUT.out_band) // pengfei 2014-10-13
+    if (INPUT.out_band[0]) // pengfei 2014-10-13
     {
         int nks = 0;
         if (nspin0 == 1)
@@ -348,7 +348,15 @@ namespace ModuleESolver
             std::stringstream ss2;
             ss2 << GlobalV::global_out_dir << "BANDS_" << is + 1 << ".dat";
             GlobalV::ofs_running << "\n Output bands in file: " << ss2.str() << std::endl;
-            ModuleIO::nscf_band(is, ss2.str(), nks, GlobalV::NBANDS, 0.0, this->pelec->ekb, this->kv, &(GlobalC::Pkpoints));
+            ModuleIO::nscf_band(is, 
+                                ss2.str(), 
+                                nks, 
+                                GlobalV::NBANDS, 
+                                0.0, 
+                                INPUT.out_band[1],
+                                this->pelec->ekb, 
+                                this->kv, 
+                                &(GlobalC::Pkpoints));
         }
     } // out_band
 

source/module_esolver/esolver_ks_pw.cpp

@@ -959,7 +959,7 @@ void ESolver_KS_PW<T, Device>::postprocess()
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
     GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
 
-    if (INPUT.out_dos != 0 || INPUT.out_band != 0)
+    if (INPUT.out_dos != 0 || INPUT.out_band[0] != 0)
     {
         GlobalV::ofs_running << "\n\n\n\n";
         GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
@@ -1001,7 +1001,7 @@ void ESolver_KS_PW<T, Device>::postprocess()
         }
     }
 
-    if (INPUT.out_band) // pengfei 2014-10-13
+    if (INPUT.out_band[0]) // pengfei 2014-10-13
     {
         int nks = 0;
         if (nspin0 == 1)
@@ -1022,6 +1022,7 @@ void ESolver_KS_PW<T, Device>::postprocess()
                                 nks,
                                 GlobalV::NBANDS,
                                 0.0,
+                                INPUT.out_band[1],
                                 this->pelec->ekb,
                                 this->kv,
                                 &(GlobalC::Pkpoints));

source/module_io/input.cpp

@@ -336,7 +336,7 @@ void Input::Default(void)
     out_wfc_pw = 0;
     out_wfc_r = 0;
     out_dos = 0;
-    out_band = 0;
+    out_band = {0, 8};
     out_proj_band = 0;
     out_mat_hs = {0, 8};
     out_mat_xc = 0;
@@ -1378,13 +1378,13 @@ bool Input::Read(const std::string& fn)
         }
         else if (strcmp("out_band", word) == 0)
         {
-            read_bool(ifs, out_band);
+            read_value2stdvector(ifs, out_band);
+            if(out_band.size() == 1) out_band.push_back(8);
         }
         else if (strcmp("out_proj_band", word) == 0)
         {
             read_bool(ifs, out_proj_band);
         }
-
         else if (strcmp("out_mat_hs", word) == 0)
         {
             read_value2stdvector(ifs, out_mat_hs);
@@ -2826,7 +2826,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
         this->relax_nmax = 1;
         out_stru = 0;
         out_dos = 0;
-        out_band = 0;
+        out_band[0] = 0;
         out_proj_band = 0;
         cal_force = 0;
         init_wfc = "file";
@@ -2843,7 +2843,7 @@ void Input::Default_2(void) // jiyy add 2019-08-04
         this->relax_nmax = 1;
         out_stru = 0;
         out_dos = 0;
-        out_band = 0;
+        out_band[0] = 0;
         out_proj_band = 0;
         cal_force = 0;
         init_wfc = "file";
@@ -3325,7 +3325,8 @@ void Input::Bcast()
     Parallel_Common::bcast_int(out_wfc_pw);
     Parallel_Common::bcast_bool(out_wfc_r);
     Parallel_Common::bcast_int(out_dos);
-    Parallel_Common::bcast_bool(out_band);
+    if(GlobalV::MY_RANK != 0) out_band.resize(2); /* If this line is absent, will cause segmentation fault in io_input_test_para */
+    Parallel_Common::bcast_int(out_band.data(), 2);
     Parallel_Common::bcast_bool(out_proj_band);
     if(GlobalV::MY_RANK != 0) out_mat_hs.resize(2); /* If this line is absent, will cause segmentation fault in io_input_test_para */
     Parallel_Common::bcast_int(out_mat_hs.data(), 2);

source/module_io/input.h

@@ -263,7 +263,7 @@ class Input
     int out_wfc_pw; // 0: no; 1: txt; 2: dat
     bool out_wfc_r; // 0: no; 1: yes
     int out_dos; // dos calculation. mohan add 20090909
-    bool out_band; // band calculation pengfei 2014-10-13
+    std::vector<int> out_band; // band calculation pengfei 2014-10-13
     bool out_proj_band; // projected band structure calculation jiyy add 2022-05-11
     std::vector<int> out_mat_hs; // output H matrix and S matrix in local basis.
     bool out_mat_xc; // output exchange-correlation matrix in KS-orbital representation.
@@ -667,7 +667,15 @@ class Input
     template <typename T>
     typename std::enable_if<std::is_same<T, double>::value, T>::type cast_string(const std::string& str) { return std::stod(str); }
     template <typename T>
-    typename std::enable_if<std::is_same<T, int>::value, T>::type cast_string(const std::string& str) { return std::stoi(str); }
+    typename std::enable_if<std::is_same<T, int>::value, T>::type cast_string(const std::string& str)
+    {
+        if (str == "true" || str == "1")
+            return 1;
+        else if (str == "false" || str == "0")
+            return 0;
+        else
+            return std::stoi(str);
+    }
     template <typename T>
     typename std::enable_if<std::is_same<T, bool>::value, T>::type cast_string(const std::string& str) { return (str == "true" || str == "1"); }
     template <typename T>

source/module_io/nscf_band.cpp

@@ -3,13 +3,15 @@
 #include "module_base/global_variable.h"
 #include "module_base/timer.h"
 #include "module_base/tool_title.h"
+#include "module_base/formatter_physfmt.h"
 
 void ModuleIO::nscf_band(
 	const int &is,
 	const std::string &out_band_dir, 
 	const int &nks, 
 	const int &nband,
 	const double &fermie,
+	const int &precision,
 	const ModuleBase::matrix& ekb,
 	const K_Vectors& kv,
 	const Parallel_Kpoints* Pkpoints)
@@ -33,23 +35,28 @@ void ModuleIO::nscf_band(
 		if (ik>0)
 		{
 			auto delta=kv.kvec_c[ik]-kv.kvec_c[ik-1];
-			klength[ik] = klength[ik-1] + delta.norm();
+			klength[ik] = klength[ik-1];
+			klength[ik] += (kv.kl_segids[ik] == kv.kl_segids[ik-1]) ? delta.norm() : 0.0;
 		}
+		/* first find if present kpoint in present pool */
 		if ( GlobalV::MY_POOL == Pkpoints->whichpool[ik] )
 		{
+			/* then get the local kpoint index, which starts definitly from 0 */
 			const int ik_now = ik - Pkpoints->startk_pool[GlobalV::MY_POOL];
+			/* if present kpoint corresponds the spin of the present one */
 			if( kv.isk[ik_now+is*nks] == is )
 			{ 
 				if ( GlobalV::RANK_IN_POOL == 0)
 				{
-					std::ofstream ofs(out_band_dir.c_str(),std::ios::app);
-					ofs << std::setprecision(8);
-					//start from 1
-					ofs << ik+1;
-					ofs << " " << klength[ik] << " ";
+					formatter::PhysicalFmt physfmt; // create a physical formatter temporarily
+					std::ofstream ofs(out_band_dir.c_str(), std::ios::app);
+					physfmt.adjust_formatter_flexible(4, 0, false); // for integer
+					ofs << physfmt.get_p_formatter()->format(ik+1);
+					physfmt.adjust_formatter_flexible(precision, 4.0/double(precision), false); // for decimal
+					ofs << physfmt.get_p_formatter()->format(klength[ik]);
 					for(int ib = 0; ib < nband; ib++)
 					{
-						ofs << " " << (ekb(ik_now+is*nks, ib)-fermie) * ModuleBase::Ry_to_eV;
+						ofs << physfmt.get_p_formatter()->format((ekb(ik_now+is*nks, ib)-fermie) * ModuleBase::Ry_to_eV);
 					}
 					ofs << std::endl;
 					ofs.close();	
@@ -83,18 +90,30 @@ void ModuleIO::nscf_band(
 #else
 //	std::cout<<"\n nband = "<<nband<<std::endl;
 //	std::cout<<out_band_dir<<std::endl;
-
+	formatter::PhysicalFmt physfmt; // create a physical formatter temporarily
+	std::vector<double> klength;
+	klength.resize(nks);
+	klength[0] = 0.0;
 	std::ofstream ofs(out_band_dir.c_str());
 	for(int ik=0;ik<nks;ik++)
 	{
+		if (ik>0)
+		{
+			auto delta=kv.kvec_c[ik]-kv.kvec_c[ik-1];
+			klength[ik] = klength[ik-1];
+			klength[ik] += (kv.kl_segids[ik] == kv.kl_segids[ik-1]) ? delta.norm() : 0.0;
+		}
 		if( kv.isk[ik] == is)
 		{
-			ofs<<std::setw(12)<<ik + 1;
+			physfmt.adjust_formatter_flexible(4, 0, false); // for integer
+			ofs << physfmt.get_p_formatter()->format(ik+1);
+			physfmt.adjust_formatter_flexible(precision, 4.0/double(precision), false); // for decimal
+			ofs << physfmt.get_p_formatter()->format(klength[ik]); // add klength, in accordance with the MPI version
 			for(int ibnd = 0; ibnd < nband; ibnd++)
 			{
-				ofs <<std::setw(15) << (ekb(ik, ibnd)-fermie) * ModuleBase::Ry_to_eV;
+				ofs << physfmt.get_p_formatter()->format((ekb(ik, ibnd)-fermie) * ModuleBase::Ry_to_eV);
 			}
-			ofs<<std::endl;
+			ofs << std::endl;
 		}
 	}
 	ofs.close();

source/module_io/nscf_band.h

@@ -12,6 +12,7 @@ namespace ModuleIO
 		const int &nks, 
 		const int &nband, 
 		const double &fermie,
+		const int &precision,
 		const ModuleBase::matrix &ekb,
 		const K_Vectors& kv,
 		const Parallel_Kpoints* Pkpoints);

source/module_io/parameter_pool.cpp

@@ -917,7 +917,7 @@ bool input_parameters_set(std::map<std::string, InputParameter> input_parameters
     }
     else if (input_parameters.count("out_band") != 0)
     {
-        INPUT.out_band = *static_cast<bool*>(input_parameters["out_band"].get());
+        INPUT.out_band = *static_cast<std::vector<int>*>(input_parameters["out_band"].get());
     }
     else if (input_parameters.count("out_proj_band") != 0)
     {