55
66namespace elecstate
77{
8+ int ElecStateLCAO::out_wfc_lcao = 0 ;
89
910// for gamma_only(double case) and multi-k(complex<double> case)
1011template <typename T> void ElecStateLCAO::cal_dm (const ModuleBase::matrix& wg,
@@ -53,21 +54,37 @@ template<typename T> void ElecStateLCAO::cal_dm(const ModuleBase::matrix& wg,
5354 return ;
5455}
5556
57+ // multi-k case
5658void ElecStateLCAO::psiToRho (const psi::Psi<std::complex <double >>& psi)
5759{
5860 ModuleBase::TITLE (" ElecStateLCAO" , " psiToRho" );
5961 ModuleBase::timer::tick (" ElecStateLCAO" , " psiToRho" );
6062
63+ ModuleBase::GlobalFunc::NOTE (" Calculate the density matrix." );
64+
65+ // this part for calculating dm_k in 2d-block format, not used for charge now
6166 psi::Psi<std::complex <double >> dm_k_2d (psi.get_nk (), psi.get_nbasis (), psi.get_nbasis ());
6267
63- ModuleBase::GlobalFunc::NOTE (" Calculate the density matrix." );
64- // this->loc->cal_dk_k(GlobalC::GridT);
6568 if (GlobalV::KS_SOLVER == " genelpa" || GlobalV::KS_SOLVER == " scalapack_gvx"
6669 || GlobalV::KS_SOLVER == " lapack" ) // Peize Lin test 2019-05-15
6770 {
6871 this ->cal_dm (this ->wg , psi, dm_k_2d);
6972 }
7073
74+ // this part for steps:
75+ // 1. psi_k transform from 2d-block to grid format
76+ // 2. psi_k_grid -> DM_R
77+ // 3. DM_R -> rho(r)
78+ if (GlobalV::KS_SOLVER == " genelpa" || GlobalV::KS_SOLVER == " scalapack_gvx" )
79+ {
80+ for (int ik = 0 ; ik<psi.get_nk (); ik++)
81+ {
82+ psi.fix_k (ik);
83+ this ->lowf ->wfc_2d_to_grid (ElecStateLCAO::out_wfc_lcao, psi.get_pointer (), this ->lowf ->wfc_k_grid [ik], ik);
84+ }
85+ }
86+
87+ this ->loc ->cal_dk_k (GlobalC::GridT);
7188 for (int is = 0 ; is < GlobalV::NSPIN; is++)
7289 {
7390 ModuleBase::GlobalFunc::ZEROS (this ->charge ->rho [is], this ->charge ->nrxx ); // mohan 2009-11-10
@@ -78,7 +95,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<std::complex<double>>& psi)
7895 // ------------------------------------------------------------
7996
8097 ModuleBase::GlobalFunc::NOTE (" Calculate the charge on real space grid!" );
81- // uhm. GK.cal_rho_k(this->loc->DM_R);
98+ this -> uhm -> GK .cal_rho_k (this ->loc ->DM_R );
8299
83100 this ->charge ->renormalize_rho ();
84101
@@ -91,12 +108,12 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
91108{
92109 ModuleBase::TITLE (" ElecStateLCAO" , " psiToRho" );
93110 ModuleBase::timer::tick (" ElecStateLCAO" , " psiToRho" );
111+
112+ this ->calculate_weights ();
94113
95114 if (GlobalV::KS_SOLVER == " genelpa" || GlobalV::KS_SOLVER == " scalapack_gvx"
96115 || GlobalV::KS_SOLVER == " lapack" )
97116 {
98- // LiuXh modify 2021-09-06, clear memory, cal_dk_gamma() not used for genelpa solver.
99- // density matrix has already been calculated.
100117 ModuleBase::timer::tick (" ElecStateLCAO" , " cal_dm_2d" );
101118
102119 psi::Psi<double > dm_gamma_2d (psi.get_nk (), psi.get_nbasis (), psi.get_nbasis ());
@@ -105,7 +122,17 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
105122
106123 ModuleBase::timer::tick (" ElecStateLCAO" , " cal_dm_2d" );
107124
108- // this->loc->cal_dk_gamma_from_2D(); // transform dm_gamma[is].c to this->loc->DM[is]
125+ for (int ik = 0 ; ik< psi.get_nk (); ++ik)
126+ {
127+ // for gamma_only case, no convertion occured, just for print.
128+ if (GlobalV::KS_SOLVER == " genelpa" || GlobalV::KS_SOLVER == " scalapack_gvx" )
129+ {
130+ psi.fix_k (ik);
131+ double ** wfc_grid = nullptr ; // output but not do "2d-to-grid" conversion
132+ this ->lowf ->wfc_2d_to_grid (ElecStateLCAO::out_wfc_lcao, psi.get_pointer (), wfc_grid);
133+ }
134+ this ->loc ->dm2dToGrid (dm_gamma_2d, this ->loc ->DM [ik]); // transform dm_gamma[is].c to this->loc->DM[is]
135+ }
109136 }
110137
111138 for (int is = 0 ; is < GlobalV::NSPIN; is++)
@@ -117,7 +144,7 @@ void ElecStateLCAO::psiToRho(const psi::Psi<double>& psi)
117144 // calculate the charge density on real space grid.
118145 // ------------------------------------------------------------
119146 ModuleBase::GlobalFunc::NOTE (" Calculate the charge on real space grid!" );
120- // uhm. GG.cal_rho(this->loc->DM);
147+ this -> uhm -> GG .cal_rho (this ->loc ->DM );
121148
122149 this ->charge ->renormalize_rho ();
123150
0 commit comments