Merge pull request #1097 from dyzheng/develop

Fix: berry phase calculation in LCAO base with nscf calculation

Author: dyzheng

source/input_conv.cpp

@@ -444,6 +444,10 @@ void Input_Conv::Convert(void)
     hsolver::HSolverLCAO::out_mat_hs = INPUT.out_mat_hs;
     hsolver::HSolverLCAO::out_mat_hsR = INPUT.out_mat_hs2; // LiuXh add 2019-07-16
     elecstate::ElecStateLCAO::out_wfc_lcao = INPUT.out_wfc_lcao;
+    if(INPUT.calculation == "nscf" && !INPUT.towannier90 && !INPUT.berry_phase)
+    {
+        elecstate::ElecStateLCAO::need_psi_grid = false;
+    }
 #endif
 
     GlobalC::en.dos_emin_ev = INPUT.dos_emin_ev;

source/module_elecstate/elecstate_lcao.cpp

@@ -8,6 +8,7 @@
 namespace elecstate
 {
 int ElecStateLCAO::out_wfc_lcao = 0;
+bool ElecStateLCAO::need_psi_grid = 1;
 
 // multi-k case
 void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
@@ -38,25 +39,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
         for (int ik = 0; ik < psi.get_nk(); ik++)
         {
             psi.fix_k(ik);
-            this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), this->lowf->wfc_k_grid[ik], ik, this->ekb, this->wg);
-            //added by zhengdy-soc, rearrange the wfc_k_grid from [up,down,up,down...] to [up,up...down,down...],
-            if(GlobalV::NSPIN==4)
-            {
-                int row = GlobalC::GridT.lgd;
-                std::vector<std::complex<double>> tmp(row);
-                for(int ib=0; ib<GlobalV::NBANDS; ib++)
-                {
-                    for(int iw=0; iw<row / GlobalV::NPOL; iw++)
-                    {
-                        tmp[iw] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL];
-                        tmp[iw + row / GlobalV::NPOL] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL + 1];
-                    }
-                    for(int iw=0; iw<row; iw++)
-                    {
-                        this->lowf->wfc_k_grid[ik][ib][iw] = tmp[iw];
-                    }
-                }
-            }
+            this->print_psi(psi);
         }
     }
 
@@ -112,8 +95,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
             if (GlobalV::KS_SOLVER == "genelpa" || GlobalV::KS_SOLVER == "scalapack_gvx")
             {
                 psi.fix_k(ik);
-                double** wfc_grid = nullptr; // output but not do "2d-to-grid" conversion
-                this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi.get_pointer(), wfc_grid, this->ekb, this->wg);
+                this->print_psi(psi);
             }
             //this->loc->dm2dToGrid(this->loc->dm_gamma[ik], this->loc->DM[ik]); // transform dm_gamma[is].c to this->loc->DM[is]
             this->loc->cal_dk_gamma_from_2D_pub();
@@ -158,11 +140,36 @@ void ElecStateLCAO::print_psi(const psi::Psi<double>& psi_in)
 }
 void ElecStateLCAO::print_psi(const psi::Psi<std::complex<double>>& psi_in)
 {
-    if(!ElecStateLCAO::out_wfc_lcao) return;
+    if(!ElecStateLCAO::out_wfc_lcao && !ElecStateLCAO::need_psi_grid) return;
 
     // output but not do "2d-to-grid" conversion
-    std::complex<double>** wfc_grid = nullptr; 
-    this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi_in.get_pointer(), wfc_grid, psi_in.get_current_k(), this->ekb, this->wg);
+    std::complex<double>** wfc_grid = nullptr;
+    int ik = psi_in.get_current_k();
+    if(ElecStateLCAO::need_psi_grid)
+    {
+        wfc_grid = this->lowf->wfc_k_grid[ik];
+    }
+    this->lowf->wfc_2d_to_grid(ElecStateLCAO::out_wfc_lcao, psi_in.get_pointer(), wfc_grid, ik, this->ekb, this->wg);
+
+    //added by zhengdy-soc, rearrange the wfc_k_grid from [up,down,up,down...] to [up,up...down,down...],
+    if(ElecStateLCAO::need_psi_grid && GlobalV::NSPIN==4)
+    {
+        int row = GlobalC::GridT.lgd;
+        std::vector<std::complex<double>> tmp(row);
+        for(int ib=0; ib<GlobalV::NBANDS; ib++)
+        {
+            for(int iw=0; iw<row / GlobalV::NPOL; iw++)
+            {
+                tmp[iw] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL];
+                tmp[iw + row / GlobalV::NPOL] = this->lowf->wfc_k_grid[ik][ib][iw * GlobalV::NPOL + 1];
+            }
+            for(int iw=0; iw<row; iw++)
+            {
+                this->lowf->wfc_k_grid[ik][ib][iw] = tmp[iw];
+            }
+        }
+    }
+
     return;
 }
 

source/module_elecstate/elecstate_lcao.h

@@ -41,6 +41,7 @@ class ElecStateLCAO : public ElecState
     virtual void print_psi(const psi::Psi<std::complex<double>>& psi_in)override;
 
     static int out_wfc_lcao;
+    static bool need_psi_grid;
 
   private:
     // calculate electronic charge density on grid points or density matrix in real space

source/module_elecstate/test/updaterhok_pw_test.cpp

@@ -12,7 +12,8 @@
 #include "module_elecstate/elecstate.h"
 
 #include "updaterhok_pw_test.h"
-
+#include "src_io/berryphase.h"
+bool berryphase::berry_phase_flag;
 using::testing::AtLeast;
 using::testing::Assign;
 

source/module_esolver/esolver_ks_lcao_elec.cpp

@@ -395,7 +395,7 @@ namespace ModuleESolver
         if (berryphase::berry_phase_flag && ModuleSymmetry::Symmetry::symm_flag == 0)
         {
             berryphase bp(this->LOWF);
-            bp.Macroscopic_polarization(nullptr);
+            bp.Macroscopic_polarization(this->psi);
         }
 
         return;

source/src_io/berryphase.cpp

@@ -325,7 +325,7 @@ double berryphase::stringPhase(int index_str, int nbands, const psi::Psi<std::co
 			if(GlobalV::NSPIN!=4)
 			{
 				//std::complex<double> my_det = lcao_method.det_berryphase(ik_1,ik_2,dk,nbands);
-				zeta = zeta * lcao_method.det_berryphase(ik_1,ik_2,dk,nbands, *this->lowf);
+				zeta = zeta * lcao_method.det_berryphase(ik_1,ik_2,dk,nbands, *this->lowf, psi_in);
 				// test by jingan
 				//GlobalV::ofs_running << "methon 1: det = " << my_det << std::endl;
 				// test by jingan
@@ -456,7 +456,7 @@ void berryphase::Macroscopic_polarization(const psi::Psi<std::complex<double>>*
 	GlobalV::ofs_running << "\n\n\n\n";
 	GlobalV::ofs_running << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>" << std::endl;
 	GlobalV::ofs_running << " |                                                                    |" << std::endl;
-	GlobalV::ofs_running << " | POLARIZATION GlobalV::CALCULATION:                                          |" << std::endl;
+	GlobalV::ofs_running << " | POLARIZATION CALCULATION:                                          |" << std::endl;
 	GlobalV::ofs_running << " |                  Modern Theory of Polarization                     |" << std::endl;
 	GlobalV::ofs_running << " | calculate the Macroscopic polarization of a crystalline insulator  |" << std::endl;
 	GlobalV::ofs_running << " | by using Berry Phase method.                                       |" << std::endl;

source/src_io/unk_overlap_lcao.cpp

@@ -723,7 +723,8 @@ void unkOverlap_lcao::prepare_midmatrix_pblas(const int ik_L, const int ik_R, co
 
 std::complex<double> unkOverlap_lcao::det_berryphase(const int ik_L, const int ik_R,
     const ModuleBase::Vector3<double> dk, const int occ_bands,
-    Local_Orbital_wfc &lowf)
+    Local_Orbital_wfc &lowf,
+	const psi::Psi<std::complex<double>>* psi_in)
 {
 	const std::complex<double> minus = std::complex<double>(-1.0,0.0);
 	std::complex<double> det = std::complex<double>(1.0,0.0);
@@ -743,14 +744,14 @@ std::complex<double> unkOverlap_lcao::det_berryphase(const int ik_L, const int i
 	int one = 1;
 #ifdef __MPI
 	pzgemm_(&transa,&transb,&occBands,&nlocal,&nlocal,&alpha[0],
-			lowf.wfc_k.at(ik_L).c,&one,&one,lowf.ParaV->desc,
+			&psi_in[0](ik_L, 0, 0), &one, &one, lowf.ParaV->desc,
 							  midmatrix,&one,&one,lowf.ParaV->desc,
 													   &beta[0],
 							   C_matrix,&one,&one,lowf.ParaV->desc);
 
 	pzgemm_(&transb,&transb,&occBands,&occBands,&nlocal,&alpha[0],
 								 C_matrix,&one,&one,lowf.ParaV->desc,
-			lowf.wfc_k.at(ik_R).c,&one,&one,lowf.ParaV->desc,
+			&psi_in[0](ik_R, 0, 0), &one, &one, lowf.ParaV->desc,
 														 &beta[0],
 							   out_matrix,&one,&one,lowf.ParaV->desc);	
 

source/src_io/unk_overlap_lcao.h

@@ -68,7 +68,8 @@ class unkOverlap_lcao
 	void prepare_midmatrix_pblas(const int ik_L, const int ik_R, const ModuleBase::Vector3<double> dk, std::complex<double> *&midmatrix, const Parallel_Orbitals &pv);
     std::complex<double> det_berryphase(const int ik_L, const int ik_R,
         const ModuleBase::Vector3<double> dk, const int occ_bands,
-        Local_Orbital_wfc &lowf);
+        Local_Orbital_wfc &lowf,
+		const psi::Psi<std::complex<double>>* psi_in);
 
 	void test(std::complex<double>*** wfc_k_grid);
 

source/src_lcao/test/gamma_rho_mock.h

@@ -3107,3 +3107,6 @@ void Pseudopot_upf::read_pseudo_upf201_rhoatom(std::ifstream &ifs)
 		ifs >> this->rho_at[ir];
 	}
 }
+
+#include "src_io/berryphase.h"
+bool berryphase::berry_phase_flag;

source/src_pw/klist.cpp

@@ -5,6 +5,7 @@
 #include "../src_parallel/parallel_reduce.h"
 #include "../src_parallel/parallel_common.h"
 #include "../module_base/memory.h"
+#include "src_io/berryphase.h"
 
 K_Vectors::K_Vectors()
 {	
@@ -80,11 +81,16 @@ void K_Vectors::set(
 	}
 
     // (2)
-    this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag);
-    if(ModuleSymmetry::Symmetry::symm_flag || is_mp)
+    //only berry phase need all kpoints including time-reversal symmetry!
+    //if symm_flag is not set, only time-reversal symmetry would be considered.
+    if(!berryphase::berry_phase_flag)
     {
-        this->update_use_ibz();
-        this->nks = this->nkstot = this->nkstot_ibz;
+        this->ibz_kpoint(symm, ModuleSymmetry::Symmetry::symm_flag);
+        if(ModuleSymmetry::Symmetry::symm_flag || is_mp)
+        {
+            this->update_use_ibz();
+            this->nks = this->nkstot = this->nkstot_ibz;
+        }
     }
 
     // (3)