Merge branch 'develop' into mdp

docs/advanced/input_files/input-main.md

@@ -147,6 +147,7 @@
     - [out\_mat\_t](#out_mat_t)
     - [out\_mat\_dh](#out_mat_dh)
     - [out\_mat\_xc](#out_mat_xc)
+    - [out\_eband\_terms](#out_eband_terms)
     - [out\_hr\_npz/out\_dm\_npz](#out_hr_npzout_dm_npz)
     - [dm\_to\_rho](#dm_to_rho)
     - [out\_app\_flag](#out_app_flag)
@@ -1680,9 +1681,16 @@ These variables are used to control the output of properties.
 ### out_mat_xc
 
 - **Type**: Boolean
-- **Availability**: Numerical atomic orbital basis
+- **Availability**: Numerical atomic orbital (NAO) and NAO-in-PW basis
 - **Description**: Whether to print the upper triangular part of the exchange-correlation matrices in **Kohn-Sham orbital representation** (unit: Ry): $\braket{\psi_i|V_\text{xc}^\text{(semi-)local}+V_\text{exx}+V_\text{DFTU}|\psi_j}$ for each k point into files in the directory `OUT.${suffix}`, which is useful for the subsequent GW calculation. (Note that currently DeePKS term is not included. ) The files are named `k-$k-Vxc`, the meaning of `$k`corresponding to k point and spin  is same as [hs_matrix.md](../elec_properties/hs_matrix.md#out_mat_hs).
-The band (KS orbital) energy for each (k-point, spin, band) will be printed in the file `OUT.${suffix}/vxc_out`. If EXX is calculated, the local and EXX part of band energy will also be printed in `OUT.${suffix}/vxc_local_out`and `OUT.${suffix}/vxc_exx_out`, respectively. All the `vxc*_out` files contains 3 integers (nk, nspin, nband) followed by nk\*nspin\*nband lines of energy Hartree and eV.
+The band (KS orbital) energy for each (k-point, spin, band) will be printed in the file `OUT.${suffix}/vxc_out.dat`. If EXX is calculated, the local and EXX part of band energy will also be printed in `OUT.${suffix}/vxc_local_out.dat`and `OUT.${suffix}/vxc_exx_out.dat`, respectively. All the `vxc*_out.dat` files contains 3 integers (nk, nspin, nband) followed by nk\*nspin\*nband lines of energy Hartree and eV.
+- **Default**: False
+
+### out_eband_terms
+
+- **Type**: Boolean
+- **Availability**: Numerical atomic orbital basis
+- **Description**: Whether to print the band energy terms separately in the file `OUT.${suffix}/${term}_out.dat`. The terms include the kinetic, pseudopotential (local + nonlocal), Hartree and exchange-correlation (including exact exchange if calculated).
 - **Default**: False
 
 ### out_hr_npz/out_dm_npz

source/Makefile.Objects

@@ -185,16 +185,16 @@ OBJS_CELL=atom_pseudo.o\
     print_cif.o\
 
 OBJS_DEEPKS=LCAO_deepks.o\
-        LCAO_deepks_fgamma.o\
-        LCAO_deepks_fk.o\
+        deepks_fgamma.o\
+        deepks_fk.o\
         LCAO_deepks_odelta.o\
         LCAO_deepks_io.o\
         LCAO_deepks_mpi.o\
         LCAO_deepks_pdm.o\
         LCAO_deepks_psialpha.o\
         LCAO_deepks_torch.o\
         LCAO_deepks_vdelta.o\
-        LCAO_deepks_hmat.o\
+        deepks_hmat.o\
         LCAO_deepks_interface.o\
         orbital_precalc.o\
         orbital_precalc_k.o\
@@ -240,12 +240,15 @@ OBJS_ESOLVER=esolver.o\
     esolver_of.o\
     esolver_of_tool.o\
     esolver_of_interface.o\
+    print_funcs.o\
 
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_elec.o\
       esolver_ks_lcao_tddft.o\
-      esolver_ks_lcao_tmpfunc.o\
+      dpks_cal_e_delta_band.o\
+      dftu_cal_occup_m.o\
       io_npz.o\
+      set_matrix_grid.o\
 
 OBJS_GINT=gint.o\
       gint_gamma_env.o\
@@ -517,8 +520,8 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\
       write_proj_band_lcao.o\
       write_istate_info.o\
       nscf_fermi_surf.o\
-      istate_charge.o\
-      istate_envelope.o\
+      get_pchg.o\
+      get_wf.o\
       io_dmk.o\
       unk_overlap_lcao.o\
       read_wfc_nao.o\

source/module_base/scalapack_connector.h

@@ -180,7 +180,21 @@ class ScalapackConnector
 		pzgeadd_(&transa, &m, &n, &alpha, a, &ia, &ja, desca, &beta, c, &ic, &jc, descc);
 	}
 
-	static inline
+    static inline
+        void gemm(
+            const char transa, const char transb,
+            const int M, const int N, const int K,
+            const double alpha,
+            const double* A, const int IA, const int JA, const int* DESCA,
+            const double* B, const int IB, const int JB, const int* DESCB,
+            const double beta,
+            double* C, const int IC, const int JC, const int* DESCC)
+    {
+        pdgemm_(&transa, &transb, &M, &N, &K, &alpha, A, &IA, &JA, DESCA,
+            B, &IB, &JB, DESCB, &beta, C, &IC, &JC, DESCC);
+    }
+
+    static inline
 	void gemm(
 		const char transa, const char transb,
 		const int M, const int N, const int K,

source/module_elecstate/module_charge/symmetry_rho.cpp

@@ -1,14 +1,13 @@
 #include "symmetry_rho.h"
+
 #include "module_hamilt_general/module_xc/xc_functional.h"
 
 Symmetry_rho::Symmetry_rho()
 {
-
 }
 
 Symmetry_rho::~Symmetry_rho()
 {
-
 }
 
 void Symmetry_rho::begin(const int& spin_now,
@@ -17,96 +16,137 @@ void Symmetry_rho::begin(const int& spin_now,
                          Parallel_Grid& Pgrid,
                          ModuleSymmetry::Symmetry& symm) const
 {
-	assert(spin_now < 4);//added by zhengdy-soc
+    assert(spin_now < 4); // added by zhengdy-soc
+
+    if (ModuleSymmetry::Symmetry::symm_flag != 1) {
+        return;
+}
+    // both parallel and serial
+    // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
+    // {
+    // 	psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm);
+    // 	if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now],
+    // rho_basis,Pgrid,symm);
+    // }
+    // else	//space group, do rho_symm in reciprocal space
+    {
+        rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]);
+        psymmg(CHR.rhog[spin_now], rho_basis, Pgrid, symm); // need to modify
+        rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]);
+
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        {
+            // Use std::vector to manage kin_g instead of raw pointer
+            std::vector<std::complex<double>> kin_g(CHR.ngmc);
+            rho_basis->real2recip(CHR.kin_r[spin_now], kin_g.data());
+            psymmg(kin_g.data(), rho_basis, Pgrid, symm);
+            rho_basis->recip2real(kin_g.data(), CHR.kin_r[spin_now]);
+        }
+    }
+    return;
+}
+
+void Symmetry_rho::begin(const int& spin_now,
+                         double** rho,
+                         std::complex<double>** rhog,
+                         int ngmc,
+                         double** kin_r,
+                         const ModulePW::PW_Basis* rho_basis,
+                         Parallel_Grid& Pgrid,
+                         ModuleSymmetry::Symmetry& symm) const
+{
+    assert(spin_now < 4); // added by zhengdy-soc
 
-	if(ModuleSymmetry::Symmetry::symm_flag != 1) return;
-	// both parallel and serial
+    if (ModuleSymmetry::Symmetry::symm_flag != 1)
+    {
+        return;
+    }
+    // both parallel and serial
     // if(symm.nrot==symm.nrotk) //pure point-group, do rho_symm in real space
     // {
     // 	psymm(CHR.rho[spin_now], rho_basis, Pgrid, symm);
-    // 	if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now], rho_basis,Pgrid,symm);
+    // 	if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) psymm(CHR.kin_r[spin_now],
+    // rho_basis,Pgrid,symm);
     // }
     // else	//space group, do rho_symm in reciprocal space
-	{
-		rho_basis->real2recip(CHR.rho[spin_now], CHR.rhog[spin_now]);
-		psymmg(CHR.rhog[spin_now], rho_basis, Pgrid, symm);	//need to modify
-		rho_basis->recip2real(CHR.rhog[spin_now], CHR.rho[spin_now]);
+    {
+        rho_basis->real2recip(rho[spin_now], rhog[spin_now]);
+        psymmg(rhog[spin_now], rho_basis, Pgrid, symm);
+        rho_basis->recip2real(rhog[spin_now], rho[spin_now]);
 
-		if(XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5) 
-		{
-			std::complex<double>* kin_g = new std::complex<double>[CHR.ngmc];
-			rho_basis->real2recip(CHR.kin_r[spin_now], kin_g);
-			psymmg(kin_g,  rho_basis,Pgrid,symm);
-			rho_basis->recip2real(kin_g, CHR.kin_r[spin_now]);
-			delete[] kin_g;
-		}	
-	}
-	return;
+        if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
+        {
+            // Use std::vector to manage kin_g instead of raw pointer
+            std::vector<std::complex<double>> kin_g(ngmc);
+            rho_basis->real2recip(kin_r[spin_now], kin_g.data());
+            psymmg(kin_g.data(), rho_basis, Pgrid, symm);
+            rho_basis->recip2real(kin_g.data(), kin_r[spin_now]);
+        }
+    }
+    return;
 }
 
-void Symmetry_rho::psymm(double* rho_part, const ModulePW::PW_Basis *rho_basis, Parallel_Grid &Pgrid, ModuleSymmetry::Symmetry &symm) const
+void Symmetry_rho::psymm(double* rho_part,
+                         const ModulePW::PW_Basis* rho_basis,
+                         Parallel_Grid& Pgrid,
+                         ModuleSymmetry::Symmetry& symm) const
 {
 #ifdef __MPI
-	// (1) reduce all rho from the first pool.
-	double* rhotot;
-	if(GlobalV::MY_RANK == 0)
-	{
-		rhotot = new double[rho_basis->nxyz];
-		ModuleBase::GlobalFunc::ZEROS(rhotot, rho_basis->nxyz);
-	}
-	Pgrid.reduce_to_fullrho(rhotot, rho_part);
+    // (1) reduce all rho from the first pool.
+    std::vector<double> rhotot;
+    if (GlobalV::MY_RANK == 0)
+    {
+        rhotot.resize(rho_basis->nxyz);
+        ModuleBase::GlobalFunc::ZEROS(rhotot.data(), rho_basis->nxyz);
+    }
+    Pgrid.reduce_to_fullrho(rhotot.data(), rho_part);
 
-	// (2)
-	if(GlobalV::MY_RANK==0)
-	{
-		symm.rho_symmetry(rhotot, rho_basis->nx, rho_basis->ny, rho_basis->nz);
+    // (2)
+    if (GlobalV::MY_RANK == 0)
+    {
+        symm.rho_symmetry(rhotot.data(), rho_basis->nx, rho_basis->ny, rho_basis->nz);
 #else
-		symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz);
+    symm.rho_symmetry(rho_part, rho_basis->nx, rho_basis->ny, rho_basis->nz);
 #endif
-		/*
-		int count = 0;
-		GlobalV::ofs_running << scientific;
-		for(int iz=0; iz<rho_basis->nz; iz++)
-		{
-			GlobalV::ofs_running << "\n iz=" << iz;
-			for(int iy=0; iy<rho_basis->ny; iy++)
-			{
-				for(int ix=0; ix<rho_basis->nx; ix++)
-				{
-					if(count%5==0) GlobalV::ofs_running << "\n";
-					++count;
-					GlobalV::ofs_running << " " << rhotot[ix*rho_basis->ny*rho_basis->nz+iy*rho_basis->nz+iz];
-				}
-			}
-		}
-		*/
-	#ifdef __MPI
-	}
-
-	// (3)
-	const int ncxy = rho_basis->nx * rho_basis->ny;
-	double* zpiece = new double[ncxy];
-	for(int iz=0; iz<rho_basis->nz; iz++)
-	{
-		//GlobalV::ofs_running << "\n iz=" << iz;
-		ModuleBase::GlobalFunc::ZEROS(zpiece, ncxy);
-		if(GlobalV::MY_RANK==0)
-		{
-			for(int ix=0; ix<rho_basis->nx; ix++)
-			{
-				for(int iy=0; iy<rho_basis->ny; iy++)
-				{
-					const int ir = ix * rho_basis->ny + iy;
-					zpiece[ir] = rhotot[ix * rho_basis->ny * rho_basis->nz + iy * rho_basis->nz + iz];
-					//rho[ir*nczp+znow] = zpiece[ir];
-				}
-			}
-		}
-		Pgrid.zpiece_to_all(zpiece,iz, rho_part);
-	}
+        /*
+        int count = 0;
+        GlobalV::ofs_running << scientific;
+        for(int iz=0; iz<rho_basis->nz; iz++)
+        {
+            GlobalV::ofs_running << "\n iz=" << iz;
+            for(int iy=0; iy<rho_basis->ny; iy++)
+            {
+                for(int ix=0; ix<rho_basis->nx; ix++)
+                {
+                    if(count%5==0) GlobalV::ofs_running << "\n";
+                    ++count;
+                    GlobalV::ofs_running << " " << rhotot[ix*rho_basis->ny*rho_basis->nz+iy*rho_basis->nz+iz];
+                }
+            }
+        }
+        */
+#ifdef __MPI
+    }
 
-	if(GlobalV::MY_RANK==0)		delete[] rhotot;
-	delete[] zpiece;
+    // (3)
+    const int ncxy = rho_basis->nx * rho_basis->ny;
+    std::vector<double> zpiece(ncxy);
+    for (int iz = 0; iz < rho_basis->nz; iz++)
+    {
+        ModuleBase::GlobalFunc::ZEROS(zpiece.data(), ncxy);
+        if (GlobalV::MY_RANK == 0)
+        {
+            for (int ix = 0; ix < rho_basis->nx; ix++)
+            {
+                for (int iy = 0; iy < rho_basis->ny; iy++)
+                {
+                    const int ir = ix * rho_basis->ny + iy;
+                    zpiece[ir] = rhotot[ix * rho_basis->ny * rho_basis->nz + iy * rho_basis->nz + iz];
+                }
+            }
+        }
+        Pgrid.zpiece_to_all(zpiece.data(), iz, rho_part);
+    }
 #endif
-	return;
-}
+    return;
+}