Merge deepmodeling/SDFT

Author: Qianruipku

docs/examples/band-struc.md

@@ -27,6 +27,7 @@ pw_diag_thr 1.0e-7
 #Parameters (File)
 init_chg file
 out_band 1
+out_proj_band 1
 
 #Parameters (Smearing)
 smearing_method gaussian
@@ -58,4 +59,35 @@ points.
 Run the program, and you will see a file named BANDS_1.dat in the output directory. Plot it
 to get energy band structure.
 
-[back to top](#band-structure)
+If "out_proj_band" set 1, it will also produce the projected band structure in a file called PBAND_1 in xml format.
+
+The PBAND_1 file starts with number of atomic orbitals in the system, the text contents of element <band structure> is the same as data in the BANDS_1.dat file, such as:
+```
+<pband>
+<nspin>1</nspin>
+<norbitals>153</norbitals>
+<band_structure nkpoints="96" nbands="50" units="eV">
+...
+
+```
+
+The rest of the files arranged in sections, each section with a header such as below:
+    
+```
+<orbital
+ index="                                       1"
+ atom_index="                                       1"
+ species="Si"
+ l="                                       0"
+ m="                                       0"
+ z="                                       1"
+>
+<data>
+...
+</data>
+
+```
+
+The shape of text contents of element <data> is (Number of k-points, Number of bands) 
+ 
+[back to top](#band-structure)

docs/input-main.md

@@ -8,7 +8,7 @@
 
 - [System variables](#system-variables)
 
-    [suffix](#suffix) | [ntype](#ntype) | [calculation](#calculation) | [symmetry](#symmetry) | [kpar](#kpar) | [latname](#latname) | [init_wfc](#init_wfc) | [init_chg](#init_chg) | [init_vel](#init_vel) | [nelec](#nelec) | [tot_magnetization](#tot-magnetization) | [dft_functional](#dft-functional) | [pseudo_type](#pseudo-type) |  [pseudo_rcut](#pseudo-rcut) | [pseudo_mesh](#pseudo_mesh) | [mem_saver](#mem-saver) | [diago_proc](#diago_proc) | [nbspline](#nbspline)
+    [suffix](#suffix) | [ntype](#ntype) | [calculation](#calculation) | [symmetry](#symmetry) | [kpar](#kpar) | [bndpar](#bndpar) | [latname](#latname) | [init_wfc](#init_wfc) | [init_chg](#init_chg) | [init_vel](#init_vel) | [nelec](#nelec) | [tot_magnetization](#tot-magnetization) | [dft_functional](#dft-functional) | [pseudo_type](#pseudo-type) |  [pseudo_rcut](#pseudo-rcut) | [pseudo_mesh](#pseudo_mesh) | [mem_saver](#mem-saver) | [diago_proc](#diago_proc) | [nbspline](#nbspline)
 
 - [Variables related to input files](#variables-related-to-input-files)
 
@@ -24,15 +24,15 @@
 
 - [Electronic structure](#electronic-structure)
 
-    [basis_type](#basis-type) | [ks_solver](#ks-solver) | [nbands](#nbands) | [nbands_istate](#nbands-istate) | [nspin](#nspin) | [occupations](#occupations) | [smearing_method](#smearing_method) | [smearing_sigma](#smearing_sigma) | [mixing_type](#mixing-type) | [mixing_beta](#mixing-beta) | [mixing_ndim](#mixing-ndim) | [mixing_gg0](#mixing-gg0) | [gamma_only](#gamma-only) | [printe](#printe) | [scf_nmax](#scf_nmax) | [scf_thr](#scf_thr) | [chg_extrap](#chg_extrap)
+    [basis_type](#basis-type) | [ks_solver](#ks-solver) | [nbands](#nbands) | [nbands_sto](#nbands_sto) | [nbands_istate](#nbands-istate) | [nspin](#nspin) | [occupations](#occupations) | [smearing_method](#smearing_method) | [smearing_sigma](#smearing_sigma) | [mixing_type](#mixing-type) | [mixing_beta](#mixing-beta) | [mixing_ndim](#mixing-ndim) | [mixing_gg0](#mixing-gg0) | [gamma_only](#gamma-only) | [printe](#printe) | [scf_nmax](#scf_nmax) | [scf_thr](#scf_thr) | [chg_extrap](#chg_extrap) | [nche_sto](#nche_sto) | [emin_sto](#emin_sto) | [emax_sto](#emax_sto) | [seed_sto](#seed_sto)
 
 - [Geometry relaxation](#geometry-relaxation)
 
     [relax_nmax](#relax_nmax) | [relax_method](#relax_method) | [relax_cg_thr](#relax_cg_thr) | [relax_bfgs_w1](#bfgs-w1) | [relax_bfgs_w2](#bfgs-w2) | [relax_bfgs_rmax](#relax_bfgs_rmax) | [relax_bfgs_rmin](#relax_bfgs_rmin) | [relax_bfgs_init](#relax_bfgs_init) | [cal_force](#cal_force) | [force_thr](#force-thr) | [force_thr_ev](#force-thr-ev) | [cal_stress](#cal_stress) | [stress_thr](#stress-thr) | [press1, press2, press3](#press) | [fixed_axes](#fixed-axes) | [cell_factor](#cell-factor)
 
 - [Variables related to output information](#variables-related-to-output-information)
 
-    [out_force](#out_force) | [out_mul](#out_mul) | [out_freq_elec](#out_freq_elec) | [out_freq_ion](#out_freq_ion) | [out_chg](#out_chg) | [out_pot](#out_pot) | [out_dm](#out-dm) | [out_wfc_pw](#out_wfc_pw) | [out_wfc_r](#out_wfc_r) | [out_wfc_lcao](#out_wfc_lcao) | [out_dos](#out-dos) | [out_band](#out-band) | [out_stru](#out-stru) | [out_level](#out_level) | [out_alllog](#out-alllog) | [out_mat_hs](#out_mat_hs) | [out_mat_r](#out_mat_r) | [out_mat_hs2](#out_mat_hs2) | [out_element_info](#out-element-info) | [restart_save](#restart_save) | [restart_load](#restart_load)
+    [out_force](#out_force) | [out_mul](#out_mul) | [out_freq_elec](#out_freq_elec) | [out_freq_ion](#out_freq_ion) | [out_chg](#out_chg) | [out_pot](#out_pot) | [out_dm](#out-dm) | [out_wfc_pw](#out_wfc_pw) | [out_wfc_r](#out_wfc_r) | [out_wfc_lcao](#out_wfc_lcao) | [out_dos](#out-dos) | [out_band](#out-band) | [out_proj_band](#out-proj-band) | [out_stru](#out-stru) | [out_level](#out_level) | [out_alllog](#out-alllog) | [out_mat_hs](#out_mat_hs) | [out_mat_r](#out_mat_r) | [out_mat_hs2](#out_mat_hs2) | [out_element_info](#out-element-info) | [restart_save](#restart_save) | [restart_load](#restart_load)
 
 - [Density of states](#density-of-states)
 
@@ -150,6 +150,12 @@ This part of variables are used to control general system parameters.
 - **Description**: devide all processors into kpar groups, and k points will be distributed among each group. The value taken should be less than or equal to the number of k points as well as the number of MPI threads.
 - **Default**: 1
 
+#### bndpar
+
+- **Type**: Integer
+- **Description**: devide all processors into bndpar groups, and bands (only stochastic orbitals now) will be distributed among each group. It should be larger than 0.
+- **Default**: 1
+
 #### latname
 
 - **Type**: String
@@ -453,6 +459,14 @@ calculations.
   - nspin=1: 1.2\*occupied_bands, occupied_bands+10)
   - nspin=2: max(1.2\*nelec, nelec+20)
 
+#### nbands_sto
+
+- **Type**: Integer
+- **Description**: 
+  - nbands_sto>0: Number of stochastic orbitals to calculate in stochastic DFT (SDFT) or mix stochastic-deterministic (MDFT). More bands obtain more precise results or smaller stochastic errors ($ \propto 1/\sqrt{N_{\chi}}$); 
+  - nbands_sto=0: Complete basis will be used to replace stochastic orbitals with the Chebyshev method (CT) and it will get the results the same as KSDFT without stochastic errors.
+- **Default**: 0
+
 #### nbands_istate
 
 - **Type**: Integer
@@ -551,6 +565,32 @@ calculations.
   - second-order: second-order extrapolation
 - **Default**:atomic
 
+#### nche_sto
+
+- **Type**: Integer
+- **Description**: Chebyshev expansion orders for SDFT, MDFT, CT methods.
+- **Default**:5
+
+#### emin_sto
+
+- **Type**: Real
+- **Description**: Trial energy to guess the lower bound of eigen energies of the Hamitonian Operator $\hat{H}$. The unit is Ry.
+- **Default**:0.0
+
+#### emax_sto
+
+- **Type**: Real
+- **Description**: Trial energy to guess the upper bound of eigen energies of the Hamitonian Operator $\hat{H}$. The unit is Ry.
+- **Default**:0.0
+
+#### seed_sto
+
+- **Type**: Integer
+- **Description**: The random seed to generate stochastic orbitals.
+  - seed_sto>=0: Stochastic orbitals have the form of $\exp(i2\pi\theta(G))$, where $\theta$ is a uniform distribution in $(0,1)$. If seed_sto=0, the seed is decided by time(NULL).
+  - seed_sto<=-1: Stochastic orbitals have the form of $\pm1$ with the equal probability. If seed_sto=-1, the seed is decided by time(NULL).
+- **Default**:0
+
 ### Geometry relaxation
 
 This part of variables are used to control the geometry relaxation.
@@ -735,6 +775,12 @@ This part of variables are used to control the output of properties.
 - **Description**: Controls whether to output the band structure. For mroe information, refer to the [worked example](examples/band-struc.md)
 - **Default**: 0
 
+#### out_proj_band
+
+- **Type**: Integer
+- **Description**: Controls whether to output the projected band structure. For mroe information, refer to the [worked example](examples/band-struc.md)
+- **Default**: 0
+
 #### out_stru
 
 - **Type**: Boolean

source/input.cpp

@@ -138,7 +138,7 @@ void Input::Default(void)
     nbands_sto = 0;
     nbands_istate = 5;
     pw_seed = 1;
-    nche_sto = 0;
+    nche_sto = 5;
     seed_sto = 0;
     bndpar = 1;
     kpar = 1;
@@ -268,6 +268,7 @@ void Input::Default(void)
     out_wfc_r = 0;
     out_dos = 0;
     out_band = 0;
+    out_proj_band = 0;
     out_mat_hs = 0;
     out_mat_hs2 = 0; // LiuXh add 2019-07-15
     out_mat_r = 0; // jingan add 2019-8-14
@@ -995,6 +996,10 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, out_band);
         }
+        else if (strcmp("out_proj_band", word) == 0)
+        {
+            read_value(ifs, out_proj_band);
+        }
 
         else if (strcmp("out_mat_hs", word) == 0)
         {
@@ -1942,6 +1947,7 @@ void Input::Bcast()
     Parallel_Common::bcast_int(out_wfc_r);
     Parallel_Common::bcast_int(out_dos);
     Parallel_Common::bcast_int(out_band);
+    Parallel_Common::bcast_int(out_proj_band);
     Parallel_Common::bcast_int(out_mat_hs);
     Parallel_Common::bcast_int(out_mat_hs2); // LiuXh add 2019-07-15
     Parallel_Common::bcast_int(out_mat_r); // jingan add 2019-8-14
@@ -2228,6 +2234,7 @@ void Input::Check(void)
         out_stru = 0;
         out_dos = 0;
         out_band = 0;
+        out_proj_band = 0;
         cal_force = 0;
         init_wfc = "file";
         init_chg = "atomic"; // useless,
@@ -2249,6 +2256,7 @@ void Input::Check(void)
         out_stru = 0;
         out_dos = 0;
         out_band = 0;
+        out_proj_band = 0;
         cal_force = 0;
         init_wfc = "file";
         init_chg = "atomic";

source/input.h

@@ -206,6 +206,7 @@ class Input
     int out_wfc_r; // 0: no; 1: yes
     int out_dos; // dos calculation. mohan add 20090909
     int out_band; // band calculation pengfei 2014-10-13
+    int out_proj_band; // projected band structure calculation jiyy add 2022-05-11
     int out_mat_hs; // output H matrix and S matrix in local basis.
     int out_mat_hs2; // LiuXh add 2019-07-16, output H(R) matrix and S(R) matrix in local basis.
     int out_mat_r; // jingan add 2019-8-14, output r(R) matrix.

source/input_conv.cpp

@@ -435,6 +435,7 @@ void Input_Conv::Convert(void)
     GlobalC::wf.out_wfc_r = INPUT.out_wfc_r;
     GlobalC::en.out_dos = INPUT.out_dos;
     GlobalC::en.out_band = INPUT.out_band;
+    GlobalC::en.out_proj_band = INPUT.out_proj_band;
 #ifdef __LCAO
     Local_Orbital_Charge::out_dm = INPUT.out_dm;
     Pdiag_Double::out_mat_hs = INPUT.out_mat_hs;

source/module_esolver/esolver_ks_lcao.cpp

@@ -143,6 +143,19 @@ namespace ModuleESolver
         int ion_step = 0;
         GlobalC::pot.init_pot(ion_step, GlobalC::pw.strucFac);
 
+        //------------------init Basis_lcao----------------------
+        // Init Basis should be put outside of Ensolver.
+        // * reading the localized orbitals/projectors
+        // * construct the interpolation tables.
+        this->Init_Basis_lcao(this->orb_con, inp, ucell);
+        //------------------init Basis_lcao----------------------
+
+        //------------------init Hamilt_lcao----------------------
+        // * allocate H and S matrices according to computational resources
+        // * set the 'trace' between local H/S and global H/S
+        this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, orb_con.ParaV);
+        //------------------init Hamilt_lcao----------------------
+
 #ifdef __MPI  
         // PLEASE simplify the Exx_Global interface
         // mohan add 2021-03-25
@@ -175,19 +188,7 @@ namespace ModuleESolver
 
         if (INPUT.dft_plus_dmft) GlobalC::dmft.init(INPUT, ucell);
 
-        //------------------init Basis_lcao----------------------
-        // Init Basis should be put outside of Ensolver.
-        // * reading the localized orbitals/projectors
-        // * construct the interpolation tables.
-        this->Init_Basis_lcao(this->orb_con, inp, ucell);
-        //------------------init Basis_lcao----------------------
-
-        //------------------init Hamilt_lcao----------------------
-        // * allocate H and S matrices according to computational resources
-        // * set the 'trace' between local H/S and global H/S
-        this->LM.divide_HS_in_frag(GlobalV::GAMMA_ONLY_LOCAL, orb_con.ParaV);
-        //------------------init Hamilt_lcao----------------------
-      //init Psi
+        //init Psi
         if (GlobalV::GAMMA_ONLY_LOCAL)
             this->LOWF.wfc_gamma.resize(GlobalV::NSPIN);
         else

source/module_orbital/ORB_gen_tables.cpp

@@ -140,7 +140,18 @@ void ORB_gen_tables::snap_psibeta_half(
 
 	//find number of projectors on atom R0
 	const int nproj = infoNL_.nproj[T0];
-	assert(nproj>0); // mohan add 2021-04-25
+	if(nproj==0)
+	{
+		if(calc_deri)
+		{
+			nlm.resize(4);
+		}
+		else
+		{
+			nlm.resize(1);
+		}
+		return;	
+	}
 
 	std::vector<bool> calproj;
 	calproj.resize(nproj);