Update out_dmk and out_dmr (#6481)

* update SIAB codes

* add README.md in tools

* update SIAB package

* small fix for the file

* add geometry index for outputting density matrix in k space

* update output of dmk

* fix bug in io_dmk_test

* update output context of dmk

* rename io_dmk to write_dmk

* update dmr and dmk parameters, add precision control of the two parameters

* update sparse matrix format

* update write_dmk

* update write dmk

* update documents for out_dmk

* fix bug about dmk and dmr

* update input-main.md

* fix a bug in out_mat_hs and out_chg that cannot read notes

* fix bug when reading the second parameter of dmk, dmr, out_chg, out_mat_hs, out_mat_tk, out_mat_l

* update tests of input parameters

* update out_chg command

* fix exx error due to nrxx>0

* fix bug

* update makefile

* update assert of nrxx, which can be 0 in some cases

Author: mohanchen

docs/advanced/input_files/input-main.md

@@ -1699,21 +1699,21 @@ These variables are used to control the output of properties.
 
 ### out_dmk
 
-- **Type**: Boolean
+- **Type**: Boolean \[Integer\](optional) 
 - **Availability**: Numerical atomic orbital basis
 - **Description**: Whether to output the density matrix for each k-point into files in the folder `OUT.${suffix}`. The files are named as:
   - For gamma only case:
-    - nspin = 1: `dms1_nao.csr`;
+    - nspin = 1 and 4: `dm_nao.csr`;
     - nspin = 2: `dms1_nao.csr` and `dms2_nao.csr` for the two spin channels. 
   - For multi-k points case:
-    - nspin = 1: `dms1k1_nao.csr`, `dms1k2_nao.csr`, ...;
-    - nspin = 2: `dms1k1_nao.csr`... and `dms2k1_nao.csr`... for the two spin channels. 
+    - nspin = 1 and 4: `dmk1_nao.csr`, `dmk2_nao.csr`, ...;
+    - nspin = 2: `dmk1s1_nao.csr`... and `dmk1s2_nao.csr`... for the two spin channels. 
 - **Default**: False
 - **Note**: In the 3.10-LTS version, the parameter is named `out_dm` and the file names are SPIN1_DM and SPIN2_DM, etc.
 
 ### out_dmr
 
-- **Type**: Boolean
+- **Type**: Boolean \[Integer\](optional)
 - **Availability**: Numerical atomic orbital basis (multi-k points)
 - **Description**: Whether to output the density matrix with Bravias lattice vector R index into files in the folder `OUT.${suffix}`. The files are named as `dmr{s}{spin index}{g}{geometry index}{_nao} + {".csr"}`. Here, 's' refers to spin, where s1 means spin up channel while s2 means spin down channel, and the sparse matrix format 'csr' is mentioned in [out_mat_hs2](#out_mat_hs2). Finally, if [out_app_flag](#out_app_flag) is set to false, the file name contains the optional 'g' index for each ionic step that may have different geometries, and if [out_app_flag](#out_app_flag) is set to true, the density matrix with respect to Bravias lattice vector R accumulates during ionic steps:
   - nspin = 1: `dmrs1_nao.csr`;

examples/relax/lcao_output/INPUT

@@ -5,10 +5,10 @@ suffix          autotest
 nbands			8
 calculation     relax #cell-relax
 ecutwfc			10
-scf_nmax		2
+scf_nmax		20
 
 basis_type		lcao
-relax_nmax      5
+relax_nmax      1
 
 cal_stress      1
 stress_thr      1e-6
@@ -25,16 +25,19 @@ relax_method    bfgs
 pseudo_dir      ../../../tests/PP_ORB	
 orbital_dir	    ../../../tests/PP_ORB
 
-nspin           4
+nspin           1
 
-out_chg         1
-out_pot         1
-out_dmk         1
-out_dmr         1
+out_mat_hs      0 // note 
+out_mat_tk      1 // note
+out_mat_l       1 // note
+out_chg         0 8 // note 
+out_pot         0
+out_dmk         0 // note 
+out_dmr         0 // note
 #out_wfc_lcao    1   // OK
-out_dos         1
-out_band        1
-out_stru        1
+out_dos         0
+out_band        0
+out_stru        0
 out_app_flag    0
 
 out_interval    1

source/Makefile.Objects

@@ -607,7 +607,7 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_eig_occ.o\
       get_pchg_lcao.o\
       get_wf_lcao.o\
-      io_dmk.o\
+      write_dmk.o\
       unk_overlap_lcao.o\
       read_wfc_nao.o\
       read_wfc_lcao.o\

source/source_base/tool_quit.cpp

@@ -91,21 +91,21 @@ void WARNING_QUIT(const std::string &file,const std::string &description,int ret
     std::cout << " https://github.com/deepmodeling/abacus-develop/issues" << std::endl;
 
 #else
-		std::cout << " " << std::endl;
-		std::cout << " ---------------------------------------------------------" << std::endl;
-		std::cout << "                         !NOTICE!                         " << std::endl;
-		std::cout << " ---------------------------------------------------------" << std::endl;
-		std::cout << " " << std::endl;
-		std::cout << " " << description << std::endl;
-		std::cout << " CHECK IN FILE : " << PARAM.globalv.global_out_dir << "warning.log" << std::endl;
-		std::cout << " " << std::endl;
-		std::cout << " For detailed manual of ABACUS, please see the website" << std::endl;
-		std::cout << " https://abacus.deepmodeling.com" << std::endl;
-		std::cout << " For any questions, propose issues on the website" << std::endl;
-		std::cout << " https://github.com/deepmodeling/abacus-develop/issues" << std::endl;
-		std::cout << " ---------------------------------------------------------" << std::endl;
-		std::cout << "                         !NOTICE!                         " << std::endl;
-		std::cout << " ---------------------------------------------------------" << std::endl;
+	std::cout << " " << std::endl;
+	std::cout << " ---------------------------------------------------------" << std::endl;
+	std::cout << "                         !NOTICE!                         " << std::endl;
+	std::cout << " ---------------------------------------------------------" << std::endl;
+	std::cout << " " << std::endl;
+	std::cout << " " << description << std::endl;
+	std::cout << " CHECK IN FILE : " << PARAM.globalv.global_out_dir << "warning.log" << std::endl;
+	std::cout << " " << std::endl;
+	std::cout << " For detailed manual of ABACUS, please see the website" << std::endl;
+	std::cout << " https://abacus.deepmodeling.com" << std::endl;
+	std::cout << " For any questions, propose issues on the website" << std::endl;
+	std::cout << " https://github.com/deepmodeling/abacus-develop/issues" << std::endl;
+	std::cout << " ---------------------------------------------------------" << std::endl;
+	std::cout << "                         !NOTICE!                         " << std::endl;
+	std::cout << " ---------------------------------------------------------" << std::endl;
 
 
 	GlobalV::ofs_running << " ---------------------------------------------------------" << std::endl;

source/source_esolver/esolver_ks_lcao.cpp

@@ -13,8 +13,7 @@
 #include "source_io/berryphase.h"
 #include "source_io/cal_ldos.h"
 #include "source_io/cube_io.h"
-#include "source_io/io_dmk.h"
-#include "source_io/io_npz.h"
+//#include "source_io/io_npz.h"
 #include "source_io/output_dmk.h"
 #include "source_io/output_log.h"
 #include "source_io/output_mat_sparse.h"

source/source_hamilt/module_xc/xc_functional_libxc_tools.cpp

@@ -16,8 +16,12 @@ std::vector<double> XC_Functional_Libxc::convert_rho(
 	#pragma omp parallel for collapse(2) schedule(static, 1024)
 	#endif
 	for( int is=0; is<nspin; ++is )
+	{
 		for( int ir=0; ir<nrxx; ++ir )
+		{
 			rho[ir*nspin+is] = chr->rho[is][ir] + 1.0/nspin*chr->rho_core[ir];
+		}
+	}
 	return rho;
 }
 
@@ -64,7 +68,9 @@ XC_Functional_Libxc::cal_gdr(
 		#pragma omp parallel for schedule(static, 1024)
 		#endif
 		for(std::size_t ir=0; ir<nrxx; ++ir)
+		{
 			rhor[ir] = rho[ir*nspin+is];
+		}
 		//------------------------------------------
 		// initialize the charge density array in reciprocal space
 		// bring electron charge density from real space to reciprocal space
@@ -90,7 +96,9 @@ std::vector<double> XC_Functional_Libxc::convert_sigma(
 	assert(nspin>0);
 	const std::size_t nrxx = gdr[0].size();
 	for(std::size_t is=1; is<nspin; ++is)
+	{
 		assert(nrxx==gdr[is].size());
+	}
 
 	std::vector<double> sigma( nrxx * ((1==nspin)?1:3) );
 	if( 1==nspin )
@@ -126,6 +134,7 @@ std::vector<double> XC_Functional_Libxc::cal_sgn(
 	const std::vector<double> &rho,
 	const std::vector<double> &sigma)
 {
+    //assert(nrxx>0); // adding this once will cause error in examples
 	std::vector<double> sgn(nrxx*nspin, 1.0);
 	// in the case of GGA correlation for polarized case,
 	// a cutoff for grho is required to ensure that libxc gives reasonable results
@@ -137,9 +146,13 @@ std::vector<double> XC_Functional_Libxc::cal_sgn(
 		for( int ir=0; ir<nrxx; ++ir )
 		{
 			if ( rho[ir*2]<rho_threshold || std::sqrt(std::abs(sigma[ir*3]))<grho_threshold )
+			{
 				sgn[ir*2] = 0.0;
+			}
 			if ( rho[ir*2+1]<rho_threshold || std::sqrt(std::abs(sigma[ir*3+2]))<grho_threshold )
+			{
 				sgn[ir*2+1] = 0.0;
+			}
 		}
 	}
 	return sgn;
@@ -158,8 +171,12 @@ double XC_Functional_Libxc::convert_etxc(
 	#pragma omp parallel for collapse(2) reduction(+:etxc) schedule(static, 256)
 	#endif
 	for( int is=0; is<nspin; ++is )
+	{
 		for( int ir=0; ir<nrxx; ++ir )
+		{
 			etxc += ModuleBase::e2 * exc[ir] * rho[ir*nspin+is] * sgn[ir*nspin+is];
+		}
+	}
 	return etxc;
 }
 
@@ -176,6 +193,7 @@ std::pair<double,ModuleBase::matrix> XC_Functional_Libxc::convert_vtxc_v(
 	const double tpiba,
 	const Charge* const chr)
 {
+    // assert(nrxx>0); // will cause error
 	double vtxc = 0.0;
 	ModuleBase::matrix v(nspin, nrxx);
 
@@ -186,9 +204,10 @@ std::pair<double,ModuleBase::matrix> XC_Functional_Libxc::convert_vtxc_v(
 	{
 		for( std::size_t ir=0; ir<nrxx; ++ir )
 		{
-			const double v_tmp = ModuleBase::e2 * vrho[ir*nspin+is] * sgn[ir*nspin+is];
+            const std::size_t index = ir*nspin+is;
+			const double v_tmp = ModuleBase::e2 * vrho[index] * sgn[index];
 			v(is,ir) += v_tmp;
-			vtxc += v_tmp * rho[ir*nspin+is];
+			vtxc += v_tmp * rho[index];
 		}
 	}
 
@@ -226,16 +245,20 @@ std::vector<std::vector<double>> XC_Functional_Libxc::cal_dh(
 	const double tpiba,
 	const Charge* const chr)
 {
+    //assert(nrxx>0); // this line will cause bug
 	std::vector<std::vector<ModuleBase::Vector3<double>>> h(
 		nspin,
 		std::vector<ModuleBase::Vector3<double>>(nrxx) );
-	if( 1==nspin )
+
+	if( nspin==1 )
 	{
 		#ifdef _OPENMP
 		#pragma omp parallel for schedule(static, 1024)
 		#endif
 		for( std::size_t ir=0; ir<nrxx; ++ir )
+		{
 			h[0][ir] = 2.0 * gdr[0][ir] * vsigma[ir] * 2.0 * sgn[ir];
+		}
 	}
 	else
 	{
@@ -254,7 +277,9 @@ std::vector<std::vector<double>> XC_Functional_Libxc::cal_dh(
 	// define two dimensional array dh [ nspin, nrxx ]
 	std::vector<std::vector<double>> dh(nspin, std::vector<double>(nrxx));
 	for( int is=0; is!=nspin; ++is )
+	{
 		XC_Functional::grad_dot( h[is].data(), dh[is].data(), chr->rhopw, tpiba);
+	}
 
 	return dh;
 }
@@ -267,11 +292,14 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(
 	const std::vector<double> &amag,
 	const ModuleBase::matrix &v)
 {
+    //assert(nrxx>0);
 	assert(PARAM.inp.nspin==4);
 	constexpr double vanishing_charge = 1.0e-10;
 	ModuleBase::matrix v_nspin4(PARAM.inp.nspin, nrxx);
 	for( int ir=0; ir<nrxx; ++ir )
+	{
 		v_nspin4(0,ir) = 0.5 * (v(0,ir)+v(1,ir));
+	}
 	if(PARAM.globalv.domag || PARAM.globalv.domag_z)
 	{
 		for( int ir=0; ir<nrxx; ++ir )
@@ -280,11 +308,13 @@ ModuleBase::matrix XC_Functional_Libxc::convert_v_nspin4(
 			{
 				const double vs = 0.5 * (v(0,ir)-v(1,ir));
 				for(int ipol=1; ipol<PARAM.inp.nspin; ++ipol)
+				{
 					v_nspin4(ipol,ir) = vs * chr->rho[ipol][ir] / amag[ir];
+				}
 			}
 		}
 	}
 	return v_nspin4;
 }
 
-#endif
+#endif

source/source_io/CMakeLists.txt

@@ -60,7 +60,7 @@ if(ENABLE_LCAO)
       read_wfc_nao.cpp
       read_wfc_lcao.cpp
       write_wfc_nao.cpp
-      io_dmk.cpp
+      write_dmk.cpp
       write_dmr.cpp
       sparse_matrix.cpp
       file_reader.cpp

source/source_io/ctrl_output_lcao.cpp

@@ -8,7 +8,7 @@
 // functions
 #include "source_io/write_dos_lcao.h" // use ModuleIO::write_dos_lcao() 
 #include "source_io/write_dmr.h" // use ModuleIO::write_dmr() 
-#include "source_io/io_dmk.h" // use ModuleIO::write_dmk()
+#include "source_io/write_dmk.h" // use ModuleIO::write_dmk()
 #include "source_io/write_HS.h" // use ModuleIO::write_hsk()
 #include "source_io/write_wfc_nao.h" // use ModuleIO::write_wfc_nao() 
 #include "source_io/output_mat_sparse.h" // use ModuleIO::output_mat_sparse() 
@@ -92,26 +92,30 @@ void ctrl_output_lcao(UnitCell& ucell,
 	//------------------------------------------------------------------
 	//! 1) Output density matrix DM(R)
 	//------------------------------------------------------------------
-    if(PARAM.inp.out_dmr)
+    if(PARAM.inp.out_dmr[0])
 	{
 		const auto& dmr_vector = pelec->get_DM()->get_DMR_vector();
-		ModuleIO::write_dmr(dmr_vector, pv,	out_app_flag,
+
+        const int precision = PARAM.inp.out_dmr[1];
+
+		ModuleIO::write_dmr(dmr_vector, precision, pv, out_app_flag,
 				ucell.get_iat2iwt(), ucell.nat, istep);
 	}
 
 	//------------------------------------------------------------------
 	//! 2) Output density matrix DM(k)
 	//------------------------------------------------------------------
-	if (PARAM.inp.out_dmk)
+	if (PARAM.inp.out_dmk[0])
 	{
 		std::vector<double> efermis(nspin == 2 ? 2 : 1);
 		for (int ispin = 0; ispin < efermis.size(); ispin++)
 		{
 			efermis[ispin] = pelec->eferm.get_efval(ispin);
 		}
-		const int precision = 3;
+		const int precision = PARAM.inp.out_dmk[1];
+
 		ModuleIO::write_dmk(pelec->get_DM()->get_DMK_vector(),
-				precision, efermis, &(ucell), pv);
+				precision, efermis, &(ucell), pv, istep);
 	}
 
     //------------------------------------------------------------------

source/source_io/module_parameter/input_parameter.h

@@ -374,8 +374,8 @@ struct Input_para
     bool out_mul = false;                 ///< qifeng add 2019-9-10
     bool out_proj_band = false;           ///< projected band structure calculation jiyy add 2022-05-11
     std::string out_level = "ie";         ///< control the output information.
-    bool out_dmk = false;                 ///< output density matrix DM(k)
-    bool out_dmr = false;                 ///< output density matrix DM(R)
+    std::vector<int> out_dmr = {0, 8};    ///< output density matrix in real space DM(R)
+    std::vector<int> out_dmk = {0, 8};    ///< output density matrix in reciprocal space DM(k)   
     bool out_bandgap = false;             ///< QO added for bandgap printing
     std::vector<int> out_mat_hs = {0, 8}; ///< output H matrix and S matrix in local basis.
     std::vector<int> out_mat_tk = {0, 8}; ///< output T(k) matrix in local basis.

source/source_io/nscf_band.cpp

@@ -21,6 +21,8 @@ void ModuleIO::nscf_band(
     ModuleBase::TITLE("ModuleIO","nscf_band");
     ModuleBase::timer::tick("ModuleIO", "nscf_band");
 
+    assert(precision>0);
+
 	GlobalV::ofs_running << "\n";
 	GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
 	GlobalV::ofs_running << " |                                                                    |" << std::endl;