11#include " esolver_ks_lcao.h"
22
3- // #include "source_base/formatter.h"
43#include " source_base/global_variable.h"
54#include " source_base/tool_title.h"
65#include " source_estate/elecstate_tools.h"
76
87#include " source_estate/module_dm/cal_dm_psi.h"
98#include " source_lcao/module_deltaspin/spin_constrain.h"
109#include " source_lcao/module_dftu/dftu.h"
11- // #include "source_io/berryphase.h"
1210#include " source_io/cube_io.h"
13- // #include "source_io/io_npz.h"
14- // #include "source_io/output_dmk.h"
1511#include " source_io/output_log.h"
16- // #include "source_io/output_mat_sparse.h"
17- // #include "source_io/output_mulliken.h"
18- // #include "source_io/output_sk.h"
1912#include " source_io/read_wfc_nao.h"
20- // #include "source_io/to_qo.h"
21- // #include "source_io/to_wannier90_lcao.h"
22- // #include "source_io/to_wannier90_lcao_in_pw.h"
23- // #include "source_io/write_HS.h"
2413#include " source_io/write_elecstat_pot.h"
2514#include " source_io/module_parameter/parameter.h"
2615
6049#include " source_lcao/module_gint/temp_gint/gint_info.h"
6150
6251#include " source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006
52+ #include " source_estate/module_dm/setup_dm.h" // setup dm from electronic wave functions
6353#include " source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao()
6454#include " source_io/ctrl_iter_lcao.h" // use ctrl_iter_lcao()
6555
66-
6756namespace ModuleESolver
6857{
6958
@@ -164,12 +153,8 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
164153 if (inp.init_wfc == " file" && inp.esolver_type != " tddft" )
165154 {
166155 if (!ModuleIO::read_wfc_nao (PARAM.globalv .global_readin_dir ,
167- this ->pv ,
168- *(this ->psi ),
169- this ->pelec ,
170- this ->pelec ->klist ->ik2iktot ,
171- this ->pelec ->klist ->get_nkstot (),
172- inp.nspin ))
156+ this ->pv , *(this ->psi ), this ->pelec , this ->pelec ->klist ->ik2iktot ,
157+ this ->pelec ->klist ->get_nkstot (), inp.nspin ))
173158 {
174159 ModuleBase::WARNING_QUIT (" ESolver_KS_LCAO" , " read electronic wave functions failed" );
175160 }
@@ -409,60 +394,35 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
409394 // call iter_init() of ESolver_KS
410395 ESolver_KS<TK>::iter_init (ucell, istep, iter);
411396
412- elecstate::DensityMatrix<TK, double >* dm
413- = dynamic_cast <const elecstate::ElecStateLCAO<TK>*>(this ->pelec )->get_DM ();
397+ // cast pointers
398+
399+ auto * estate = dynamic_cast <elecstate::ElecStateLCAO<TK>*>(this ->pelec );
400+
401+ if (!estate)
402+ {
403+ ModuleBase::WARNING_QUIT (" ESolver_KS_LCAO::iter_init" ," pelec does not exist" );
404+ }
405+
406+ elecstate::DensityMatrix<TK, double >* dm = estate->get_DM ();
414407
415408 module_charge::chgmixing_ks_lcao (iter, this ->p_chgmix , dm->get_DMR_pointer (1 )->get_nnr (), PARAM.inp );
416409
417410 // mohan update 2012-06-05
418- this -> pelec -> f_en .deband_harris = this -> pelec ->cal_delta_eband (ucell);
411+ estate-> f_en .deband_harris = estate ->cal_delta_eband (ucell);
419412
420- // first need to calculate the weight according to
421- // electrons number.
422413 if (istep == 0 && PARAM.inp .init_wfc == " file" )
423- {
424- int exx_two_level_step = 0 ;
414+ {
415+ int exx_two_level_step = 0 ;
425416#ifdef __EXX
426- if (GlobalC::exx_info.info_global .cal_exx )
427- {
428- // the following steps are only needed in the first outer exx loop
429- exx_two_level_step
430- = GlobalC::exx_info.info_ri .real_number ? this ->exd ->two_level_step : this ->exc ->two_level_step ;
431- }
417+ if (GlobalC::exx_info.info_global .cal_exx )
418+ {
419+ // the following steps are only needed in the first outer exx loop
420+ exx_two_level_step
421+ = GlobalC::exx_info.info_ri .real_number ? this ->exd ->two_level_step : this ->exc ->two_level_step ;
422+ }
432423#endif
433- if (iter == 1 && exx_two_level_step == 0 )
434- {
435- std::cout << " WAVEFUN -> CHARGE " << std::endl;
436-
437- // calculate the density matrix using read in wave functions
438- // and then calculate the charge density on grid.
439-
440- this ->pelec ->skip_weights = true ;
441- elecstate::calculate_weights (this ->pelec ->ekb ,
442- this ->pelec ->wg ,
443- this ->pelec ->klist ,
444- this ->pelec ->eferm ,
445- this ->pelec ->f_en ,
446- this ->pelec ->nelec_spin ,
447- this ->pelec ->skip_weights );
448-
449- auto _pelec = dynamic_cast <elecstate::ElecStateLCAO<TK>*>(this ->pelec );
450- elecstate::calEBand (_pelec->ekb , _pelec->wg , _pelec->f_en );
451- elecstate::cal_dm_psi (_pelec->DM ->get_paraV_pointer (), _pelec->wg , *this ->psi , *(_pelec->DM ));
452- _pelec->DM ->cal_DMR ();
453-
454- this ->pelec ->psiToRho (*this ->psi );
455- this ->pelec ->skip_weights = false ;
456-
457- elecstate::cal_ux (ucell);
458-
459- // ! update the potentials by using new electron charge density
460- this ->pelec ->pot ->update_from_charge (&this ->chr , &ucell);
461-
462- // ! compute the correction energy for metals
463- this ->pelec ->f_en .descf = this ->pelec ->cal_delta_escf ();
464- }
465- }
424+ elecstate::setup_dm<TK>(ucell, estate, this ->psi , this ->chr , iter, exx_two_level_step);
425+ }
466426
467427#ifdef __EXX
468428 // calculate exact-exchange
@@ -523,7 +483,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
523483{
524484 ModuleBase::TITLE (" ESolver_KS_LCAO" , " hamilt2rho_single" );
525485
526- // i1 ) reset energy
486+ // 1 ) reset energy
527487 this ->pelec ->f_en .eband = 0.0 ;
528488 this ->pelec ->f_en .demet = 0.0 ;
529489 bool skip_charge = PARAM.inp .calculation == " nscf" ? true : false ;
@@ -618,7 +578,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
618578
619579 const std::vector<std::vector<TK>>& dm_vec = estate->get_DM ()->get_DMK_vector ();
620580
621-
622581 // 1) calculate the local occupation number matrix and energy correction in DFT+U
623582 if (PARAM.inp .dft_plus_u )
624583 {
@@ -628,12 +587,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
628587 {
629588 if (GlobalC::dftu.omc != 2 )
630589 {
631- ModuleDFTU::dftu_cal_occup_m (iter,
632- ucell,
633- dm_vec,
634- this ->kv ,
635- this ->p_chgmix ->get_mixing_beta (),
636- hamilt_lcao);
590+ ModuleDFTU::dftu_cal_occup_m (iter, ucell, dm_vec, this ->kv ,
591+ this ->p_chgmix ->get_mixing_beta (), hamilt_lcao);
637592 }
638593 GlobalC::dftu.cal_energy_correction (ucell, istep);
639594 }
0 commit comments