Skip to content

Commit 59d1891

Browse files
mohanchenmohanchen
committed
solve conflicts
2 parents 540d845 + fc7bc7c commit 59d1891

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

46 files changed

+557
-504
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -643,6 +643,7 @@ OBJS_LCAO=evolve_elec.o\
643643
LCAO_init_basis.o\
644644
setup_exx.o\
645645
setup_deepks.o\
646+
setup_dm.o\
646647
rho_tau_lcao.o\
647648
center2_orb.o\
648649
center2_orb-orb11.o\

‎source/source_esolver/esolver_dm2rho.cpp‎

Lines changed: 3 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -46,30 +46,22 @@ void ESolver_DM2rho<TK, TR>::runner(UnitCell& ucell, const int istep)
4646

4747
ESolver_KS_LCAO<TK, TR>::before_scf(ucell, istep);
4848

49-
auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
50-
51-
if(!estate)
52-
{
53-
ModuleBase::WARNING_QUIT("ESolver_DM2rho::after_scf","pelec does not exist");
54-
}
55-
5649
// file name of DM
5750
std::string zipname = "output_DM0.npz";
58-
elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
5951

6052
// read DM from file
61-
ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(1)));
53+
ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(1)));
6254

6355
// if nspin=2, need extra reading
6456
if (PARAM.inp.nspin == 2)
6557
{
6658
zipname = "output_DM1.npz";
67-
ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(dm->get_DMR_pointer(2)));
59+
ModuleIO::read_mat_npz(&(this->pv), ucell, zipname, *(this->dmat.dm->get_DMR_pointer(2)));
6860
}
6961

7062
// it's dangerous to design psiToRho function like this, mohan note 20251024
7163
// this->pelec->psiToRho(*this->psi);
72-
LCAO_domain::dm2rho(estate->DM->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
64+
LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, &this->chr);
7365

7466
int nspin0 = PARAM.inp.nspin == 2 ? 2 : 1;
7567

‎source/source_esolver/esolver_double_xc.cpp‎

Lines changed: 27 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,7 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
8484
}
8585

8686
// 6) initialize the density matrix
87-
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin);
87+
this->dmat_base.allocate_dm(&this->kv, &this->pv, PARAM.inp.nspin);
8888

8989
// 10) inititlize the charge density
9090
this->chr_base.allocate(PARAM.inp.nspin);
@@ -138,8 +138,6 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
138138
}
139139
if (this->p_hamilt_base == nullptr)
140140
{
141-
elecstate::DensityMatrix<TK, double>* DM = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM();
142-
143141
this->p_hamilt_base = new hamilt::HamiltLCAO<TK, TR>(
144142
ucell,
145143
this->gd,
@@ -148,7 +146,7 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
148146
this->kv,
149147
this->two_center_bundle_,
150148
this->orb_,
151-
DM,
149+
this->dmat_base.dm,
152150
this->deepks,
153151
istep,
154152
this->exx_nao);
@@ -159,13 +157,11 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
159157
XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
160158

161159
// DMR should be same size with Hamiltonian(R)
162-
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)
163-
->get_DM()
164-
->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
160+
this->dmat_base.dm->init_DMR(*(dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt_base)->getHR()));
165161

166162
if (istep > 0)
167163
{
168-
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec_base)->get_DM()->cal_DMR();
164+
this->dmat_base.dm->cal_DMR();
169165
}
170166

171167
ModuleBase::timer::tick("ESolver_DoubleXC", "before_scf");
@@ -226,23 +222,23 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
226222
std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&this->deepks.ld, [](LCAO_Deepks<TK>*) {});
227223
LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
228224

229-
deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
230-
this->kv.get_nks(),
231-
ucell.nat,
232-
PARAM.globalv.nlocal,
233-
this->pelec->ekb,
234-
this->kv.kvec_d,
235-
ucell,
236-
this->orb_,
237-
this->gd,
238-
&(this->pv),
239-
*(this->psi),
240-
dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM(),
241-
p_ham_deepks,
242-
iter,
243-
conv_esolver,
244-
GlobalV::MY_RANK,
245-
GlobalV::ofs_running);
225+
deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
226+
this->kv.get_nks(),
227+
ucell.nat,
228+
PARAM.globalv.nlocal,
229+
this->pelec->ekb,
230+
this->kv.kvec_d,
231+
ucell,
232+
this->orb_,
233+
this->gd,
234+
&(this->pv),
235+
*(this->psi),
236+
this->dmat_base.dm,
237+
p_ham_deepks,
238+
iter,
239+
conv_esolver,
240+
GlobalV::MY_RANK,
241+
GlobalV::ofs_running);
246242
#endif
247243

248244
// restore to density after charge mixing
@@ -352,9 +348,12 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
352348
auto _pes_lcao = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
353349
for (int ik = 0; ik < nks; ik++)
354350
{
355-
_pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
351+
// mohan update 2025-11-03
352+
this->dmat_base.dm->set_DMK_pointer(ik, this->dmat.dm->get_DMK_pointer(ik));
353+
// _pes_lcao_base->get_DM()->set_DMK_pointer(ik, _pes_lcao->get_DM()->get_DMK_pointer(ik));
356354
}
357-
_pes_lcao_base->get_DM()->cal_DMR();
355+
this->dmat_base.dm->cal_DMR();
356+
// _pes_lcao_base->get_DM()->cal_DMR();
358357
_pes_lcao_base->ekb = _pes_lcao->ekb;
359358
_pes_lcao_base->wg = _pes_lcao->wg;
360359
}
@@ -386,6 +385,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
386385
this->gd,
387386
this->pv,
388387
this->pelec_base,
388+
this->dmat_base, // mohan add 2025-11-03
389389
this->psi,
390390
this->two_center_bundle_,
391391
this->orb_,

‎source/source_esolver/esolver_double_xc.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,9 @@ class ESolver_DoubleXC : public ESolver_KS_LCAO<TK, TR>
3232
//! Electronic states
3333
elecstate::ElecState* pelec_base = nullptr;
3434

35+
//! Density Matrix, mohan add 2025-11-03
36+
LCAO_domain::Setup_DM<TK> dmat_base;
37+
3538
//! Electorn charge density
3639
Charge chr_base;
3740
};

‎source/source_esolver/esolver_ks_lcao.cpp‎

Lines changed: 23 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
9191

9292

9393
// 7) init DMK, but DMR is constructed in before_scf()
94-
dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), inp.nspin);
94+
this->dmat.allocate_dm(&this->kv, &this->pv, inp.nspin);
9595

9696
// 8) init exact exchange calculations
9797
this->exx_nao.before_runner(ucell, this->kv, this->orb_, this->pv, inp);
@@ -153,12 +153,6 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
153153
//! 1) call before_scf() of ESolver_KS.
154154
ESolver_KS<TK>::before_scf(ucell, istep);
155155

156-
auto* estate = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec);
157-
if(!estate)
158-
{
159-
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","pelec does not exist");
160-
}
161-
162156
//! 2) find search radius
163157
double search_radius = atom_arrange::set_sr_NL(GlobalV::ofs_running,
164158
PARAM.inp.out_level, orb_.get_rcutmax_Phi(), ucell.infoNL.get_rcutmax_Beta(),
@@ -193,11 +187,9 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
193187
}
194188
if (this->p_hamilt == nullptr)
195189
{
196-
elecstate::DensityMatrix<TK, double>* DM = estate->get_DM();
197-
198190
this->p_hamilt = new hamilt::HamiltLCAO<TK, TR>(
199191
ucell, this->gd, &this->pv, this->pelec->pot, this->kv,
200-
two_center_bundle_, orb_, DM, this->deepks, istep, exx_nao);
192+
two_center_bundle_, orb_, this->dmat.dm, this->deepks, istep, exx_nao);
201193
}
202194

203195
// 9) for each ionic step, the overlap <phi|alpha> must be rebuilt
@@ -210,7 +202,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
210202
spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
211203
sc.init_sc(PARAM.inp.sc_thr, PARAM.inp.nsc, PARAM.inp.nsc_min, PARAM.inp.alpha_trial,
212204
PARAM.inp.sccut, PARAM.inp.sc_drop_thr, ucell, &(this->pv),
213-
PARAM.inp.nspin, this->kv, this->p_hamilt, this->psi, this->pelec);
205+
PARAM.inp.nspin, this->kv, this->p_hamilt, this->psi, this->dmat.dm, this->pelec);
214206
}
215207

216208
// 11) set xc type before the first cal of xc in pelec->init_scf, Peize Lin add 2016-12-03
@@ -226,7 +218,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
226218
{
227219
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_scf","p_hamilt does not exist");
228220
}
229-
estate->get_DM()->init_DMR(*hamilt_lcao->getHR());
221+
this->dmat.dm->init_DMR(*hamilt_lcao->getHR());
230222

231223
#ifdef __MLALGO
232224
// 14) initialize DMR of DeePKS
@@ -238,7 +230,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
238230
// 2. DMK in DensityMatrix is empty (istep == 0), then DMR is initialized by zeros
239231
if (istep > 0)
240232
{
241-
estate->get_DM()->cal_DMR();
233+
this->dmat.dm->cal_DMR();
242234
}
243235

244236
// 16) the electron charge density should be symmetrized,
@@ -282,7 +274,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
282274

283275
fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress,
284276
PARAM.inp.test_force, PARAM.inp.test_stress,
285-
this->gd, this->pv, this->pelec, this->psi,
277+
this->gd, this->pv, this->pelec, this->dmat, this->psi,
286278
two_center_bundle_, orb_, force, this->scs,
287279
this->locpp, this->sf, this->kv,
288280
this->pw_rho, this->solvent, this->deepks,
@@ -341,7 +333,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
341333
}
342334

343335
ModuleIO::ctrl_runner_lcao<TK, TR>(ucell,
344-
PARAM.inp, this->kv, estate, this->pv, this->Pgrid,
336+
PARAM.inp, this->kv, estate, this->dmat, this->pv, this->Pgrid,
345337
this->gd, this->psi, this->chr, hamilt_lcao,
346338
this->two_center_bundle_,
347339
this->orb_, this->pw_rho, this->pw_rhod,
@@ -365,9 +357,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
365357
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_init","pelec does not exist");
366358
}
367359

368-
elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
369-
370-
module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp);
360+
module_charge::chgmixing_ks_lcao(iter, this->p_chgmix, this->dmat.dm->get_DMR_pointer(1)->get_nnr(), PARAM.inp);
371361

372362
// mohan update 2012-06-05
373363
estate->f_en.deband_harris = estate->cal_delta_eband(ucell);
@@ -384,7 +374,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
384374
this->exx_nao.exd->two_level_step : this->exx_nao.exc->two_level_step;
385375
}
386376
#endif
387-
elecstate::init_dm<TK>(ucell, estate, this->psi, this->chr, iter, exx_two_level_step);
377+
elecstate::init_dm<TK>(ucell, estate, this->dmat, this->psi, this->chr, iter, exx_two_level_step);
388378
}
389379

390380
#ifdef __EXX
@@ -393,11 +383,11 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
393383
{
394384
if (GlobalC::exx_info.info_ri.real_number)
395385
{
396-
this->exx_nao.exd->exx_eachiterinit(istep, ucell, *dm, this->kv, iter);
386+
this->exx_nao.exd->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, iter);
397387
}
398388
else
399389
{
400-
this->exx_nao.exc->exx_eachiterinit(istep, ucell, *dm, this->kv, iter);
390+
this->exx_nao.exc->exx_eachiterinit(istep, ucell, *this->dmat.dm, this->kv, iter);
401391
}
402392
}
403393
#endif
@@ -406,7 +396,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
406396
{
407397
if (istep != 0 || iter != 1)
408398
{
409-
GlobalC::dftu.set_dmr(dm);
399+
GlobalC::dftu.set_dmr(this->dmat.dm);
410400
}
411401
// Calculate U and J if Yukawa potential is used
412402
GlobalC::dftu.cal_slater_UJ(ucell, this->chr.rho, this->pw_rho->nrxx);
@@ -432,7 +422,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
432422
// save density matrix DMR for mixing
433423
if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0)
434424
{
435-
dm->save_DMR();
425+
this->dmat.dm->save_DMR();
436426
}
437427
}
438428

@@ -470,7 +460,8 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
470460
if (!skip_solve)
471461
{
472462
hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
473-
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, this->chr, PARAM.inp.nspin, skip_charge);
463+
hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, *this->dmat.dm,
464+
this->chr, PARAM.inp.nspin, skip_charge);
474465
}
475466

476467
// 4) EXX
@@ -518,7 +509,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
518509
ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::iter_finish","p_hamilt does not exist");
519510
}
520511

521-
const std::vector<std::vector<TK>>& dm_vec = estate->get_DM()->get_DMK_vector();
512+
const std::vector<std::vector<TK>>& dm_vec = this->dmat.dm->get_DMK_vector();
522513

523514
// 1) calculate the local occupation number matrix and energy correction in DFT+U
524515
if (PARAM.inp.dft_plus_u)
@@ -559,8 +550,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
559550
{
560551
if (PARAM.inp.mixing_restart > 0 && this->p_chgmix->mixing_restart_count > 0 && PARAM.inp.mixing_dmr)
561552
{
562-
elecstate::DensityMatrix<TK, double>* dm = estate->get_DM();
563-
this->p_chgmix->mix_dmr(dm);
553+
this->p_chgmix->mix_dmr(this->dmat.dm);
564554
}
565555
}
566556

@@ -571,7 +561,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
571561
}
572562

573563
// control the output related to the finished iteration
574-
ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, PARAM.inp, this->kv, estate,
564+
ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, PARAM.inp, this->kv, estate, *this->dmat.dm,
575565
this->pv, this->gd, this->psi, this->chr, this->p_chgmix,
576566
hamilt_lcao, this->orb_, this->deepks,
577567
this->exx_nao, iter, istep, conv_esolver, this->scf_ene_thr);
@@ -599,16 +589,19 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
599589

600590
if (PARAM.inp.out_elf[0] > 0)
601591
{
602-
LCAO_domain::dm2tau(estate->DM->get_DMR_vector(), PARAM.inp.nspin, estate->charge);
592+
LCAO_domain::dm2tau(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, this->pelec->charge);
603593
}
604594

605595
//! 1) call after_scf() of ESolver_KS
606596
ESolver_KS<TK>::after_scf(ucell, istep, conv_esolver);
607597

608598

609599
//! 2) output of lcao every few ionic steps
600+
601+
602+
610603
ModuleIO::ctrl_scf_lcao<TK, TR>(ucell,
611-
PARAM.inp, this->kv, estate, this->pv,
604+
PARAM.inp, this->kv, estate, this->dmat.dm, this->pv,
612605
this->gd, this->psi, hamilt_lcao,
613606
this->two_center_bundle_,
614607
this->orb_, this->pw_wfc, this->pw_rho,

‎source/source_esolver/esolver_ks_lcao.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include "source_lcao/setup_deepks.h" // for deepks, mohan add 20251008
1010
#include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
1111
#include "source_lcao/module_rdmft/rdmft.h" // rdmft
12+
#include "source_lcao/setup_dm.h" // mohan add 2025-10-30
1213

1314
#include <memory>
1415

@@ -71,6 +72,9 @@ class ESolver_KS_LCAO : public ESolver_KS<TK>
7172
//! NAO orbitals: two-center integrations
7273
TwoCenterBundle two_center_bundle_;
7374

75+
//! Add density matrix class, mohan add 2025-10-30
76+
LCAO_domain::Setup_DM<TK> dmat;
77+
7478
//! For RDMFT calculations, added by jghan, 2024-03-16
7579
rdmft::RDMFT<TK, TR> rdmft_solver;
7680

0 commit comments

Comments
 (0)