Skip to content

Commit 86e0291

Browse files
maki49maki49
committed
extract edm
1 parent efb0990 commit 86e0291

File tree

8 files changed

+125
-84
lines changed

8 files changed

+125
-84
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -560,6 +560,7 @@ OBJS_LCAO=evolve_elec.o\
560560
FORCE_gamma.o\
561561
FORCE_k.o\
562562
stress_tools.o\
563+
edm.o\
563564
fedm_gamma.o\
564565
fedm_k.o\
565566
ftvnl_dphi_gamma.o\

‎source/module_hamilt_lcao/hamilt_lcaodft/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@ if(ENABLE_LCAO)
1919
FORCE_gamma.cpp
2020
FORCE_k.cpp
2121
stress_tools.cpp
22+
edm.cpp
2223
fedm_gamma.cpp
2324
fedm_k.cpp
2425
ftvnl_dphi_gamma.cpp

‎source/module_hamilt_lcao/hamilt_lcaodft/FORCE.h‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -94,7 +94,7 @@ class Force_LCAO
9494
const bool isstress,
9595
ForceStressArrays& fsr,
9696
const UnitCell& ucell,
97-
const elecstate::DensityMatrix<T, double>* dm,
97+
const elecstate::DensityMatrix<T, double>& dm,
9898
const psi::Psi<T>* psi,
9999
const Parallel_Orbitals& pv,
100100
const elecstate::ElecState* pelec,
@@ -136,6 +136,16 @@ class Force_LCAO
136136
typename TGint<T>::type& gint,
137137
ModuleBase::matrix& fvl_dphi,
138138
ModuleBase::matrix& svl_dphi);
139+
140+
elecstate::DensityMatrix<T, double> cal_edm(const elecstate::ElecState* pelec,
141+
const psi::Psi<T>& psi,
142+
const elecstate::DensityMatrix<T, double>& dm,
143+
const K_Vectors& kv,
144+
const Parallel_Orbitals& pv,
145+
const int& nspin,
146+
const int& nbands,
147+
const UnitCell& ucell,
148+
Record_adj& ra) const;
139149
};
140150

141151
#endif

‎source/module_hamilt_lcao/hamilt_lcaodft/FORCE_gamma.cpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -215,7 +215,9 @@ void Force_LCAO<double>::ftable(const bool isforce,
215215
this->allocate(pv, fsr, two_center_bundle, orb);
216216

217217
// calculate the force related to 'energy density matrix'.
218-
this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap);
218+
this->cal_fedm(isforce, isstress, fsr, ucell,
219+
this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra),
220+
psi, pv, pelec, foverlap, soverlap);
219221

220222
this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi);
221223

‎source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -312,7 +312,9 @@ void Force_LCAO<std::complex<double>>::ftable(const bool isforce,
312312

313313
// calculate the energy density matrix
314314
// and the force related to overlap matrix and energy density matrix.
315-
this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap, kv, ra);
315+
this->cal_fedm(isforce, isstress, fsr, ucell,
316+
this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra),
317+
psi, pv, pelec, foverlap, soverlap, kv, ra);
316318

317319
this->cal_ftvnl_dphi(dm, pv, ucell, fsr, isforce, isstress, ftvnl_dphi, stvnl_dphi, ra);
318320

‎source/module_hamilt_lcao/hamilt_lcaodft/edm.cpp‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
#include "FORCE.h"
2+
#include "module_elecstate/module_dm/cal_dm_psi.h"
3+
#include "module_base/memory.h"
4+
template<>
5+
elecstate::DensityMatrix<double, double> Force_LCAO<double>::cal_edm(const elecstate::ElecState* pelec,
6+
const psi::Psi<double>& psi,
7+
const elecstate::DensityMatrix<double, double>& dm,
8+
const K_Vectors& kv,
9+
const Parallel_Orbitals& pv,
10+
const int& nspin,
11+
const int& nbands,
12+
const UnitCell& ucell,
13+
Record_adj& ra) const
14+
{
15+
ModuleBase::matrix wg_ekb;
16+
wg_ekb.create(nspin, nbands);
17+
18+
for(int is=0; is<nspin; is++)
19+
{
20+
for(int ib=0; ib<nbands; ib++)
21+
{
22+
wg_ekb(is,ib) = pelec->wg(is,ib) * pelec->ekb(is, ib);
23+
}
24+
}
25+
26+
// construct a DensityMatrix for Gamma-Only
27+
elecstate::DensityMatrix<double, double> edm(&pv, nspin);
28+
29+
#ifdef __PEXSI
30+
if (GlobalV::KS_SOLVER == "pexsi")
31+
{
32+
auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
33+
for (int ik = 0; ik < nspin; ik++)
34+
{
35+
edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]);
36+
}
37+
38+
}
39+
else
40+
#endif
41+
{
42+
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm);
43+
}
44+
return edm;
45+
}
46+
47+
template<>
48+
elecstate::DensityMatrix<std::complex<double>, double> Force_LCAO<std::complex<double>>::cal_edm(const elecstate::ElecState* pelec,
49+
const psi::Psi<std::complex<double>>& psi,
50+
const elecstate::DensityMatrix<std::complex<double>, double>& dm,
51+
const K_Vectors& kv,
52+
const Parallel_Orbitals& pv,
53+
const int& nspin,
54+
const int& nbands,
55+
const UnitCell& ucell,
56+
Record_adj& ra) const
57+
{
58+
59+
// construct a DensityMatrix object
60+
elecstate::DensityMatrix<std::complex<double>, double> edm(&kv, &pv, nspin);
61+
62+
//--------------------------------------------
63+
// calculate the energy density matrix here.
64+
//--------------------------------------------
65+
66+
ModuleBase::matrix wg_ekb;
67+
wg_ekb.create(kv.get_nks(), nbands);
68+
ModuleBase::Memory::record("Force::wg_ekb", sizeof(double) * kv.get_nks() * nbands);
69+
#ifdef _OPENMP
70+
#pragma omp parallel for collapse(2) schedule(static, 1024)
71+
#endif
72+
for (int ik = 0; ik < kv.get_nks(); ik++)
73+
{
74+
for (int ib = 0; ib < nbands; ib++)
75+
{
76+
wg_ekb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib);
77+
}
78+
}
79+
80+
// use the original formula (Hamiltonian matrix) to calculate energy density matrix
81+
if (dm.EDMK.size())
82+
{
83+
#ifdef _OPENMP
84+
#pragma omp parallel for schedule(static)
85+
#endif
86+
for (int ik = 0; ik < kv.get_nks(); ++ik)
87+
{
88+
edm.set_DMK_pointer(ik, dm.EDMK[ik].c);
89+
}
90+
}
91+
else
92+
{
93+
// cal_dm_psi
94+
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi, edm);
95+
}
96+
97+
// cal_dm_2d
98+
edm.init_DMR(ra, &ucell);
99+
edm.cal_DMR();
100+
edm.sum_DMR_spin();
101+
return edm;
102+
}

‎source/module_hamilt_lcao/hamilt_lcaodft/fedm_gamma.cpp‎

Lines changed: 2 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ void Force_LCAO<double>::cal_fedm(
1616
const bool isstress,
1717
ForceStressArrays &fsr, // mohan add 2024-06-16
1818
const UnitCell& ucell,
19-
const elecstate::DensityMatrix<double, double>* dm,
19+
const elecstate::DensityMatrix<double, double>& dm,
2020
const psi::Psi<double>* psi,
2121
const Parallel_Orbitals& pv,
2222
const elecstate::ElecState *pelec,
@@ -32,37 +32,6 @@ void Force_LCAO<double>::cal_fedm(
3232
const int nbands = GlobalV::NBANDS;
3333
const int nlocal = GlobalV::NLOCAL;
3434

35-
ModuleBase::matrix wg_ekb;
36-
wg_ekb.create(nspin, nbands);
37-
38-
for(int is=0; is<nspin; is++)
39-
{
40-
for(int ib=0; ib<nbands; ib++)
41-
{
42-
wg_ekb(is,ib) = pelec->wg(is,ib) * pelec->ekb(is, ib);
43-
}
44-
}
45-
46-
// construct a DensityMatrix for Gamma-Only
47-
elecstate::DensityMatrix<double, double> edm(&pv, nspin);
48-
49-
#ifdef __PEXSI
50-
if (GlobalV::KS_SOLVER == "pexsi")
51-
{
52-
auto pes = dynamic_cast<const elecstate::ElecStateLCAO<double>*>(pelec);
53-
for (int ik = 0; ik < nspin; ik++)
54-
{
55-
edm.set_DMK_pointer(ik, pes->get_DM()->pexsi_EDM[ik]);
56-
}
57-
58-
}
59-
else
60-
#endif
61-
{
62-
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi[0], edm);
63-
}
64-
65-
6635
for(int i=0; i<nlocal; i++)
6736
{
6837
const int iat = ucell.iwt2iat[i];
@@ -78,7 +47,7 @@ void Force_LCAO<double>::cal_fedm(
7847
for(int is=0; is<nspin; ++is)
7948
{
8049
//sum += edm_gamma[is](nu, mu);
81-
sum += edm.get_DMK(is+1, 0, nu, mu);
50+
sum += dm.get_DMK(is+1, 0, nu, mu);
8251
}
8352
sum *= 2.0;
8453

‎source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp‎

Lines changed: 2 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,6 @@
66
#include "module_cell/module_neighbor/sltk_grid_driver.h"
77
#include "module_elecstate/cal_dm.h"
88
#include "module_elecstate/elecstate_lcao.h"
9-
#include "module_elecstate/module_dm/cal_dm_psi.h"
109
#include "module_hamilt_pw/hamilt_pwdft/global.h"
1110
#include "module_io/write_HS.h"
1211
#include "module_hamilt_lcao/hamilt_lcaodft/stress_tools.h"
@@ -29,7 +28,7 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
2928
const bool isstress,
3029
ForceStressArrays& fsr,
3130
const UnitCell& ucell,
32-
const elecstate::DensityMatrix<std::complex<double>, double>* dm,
31+
const elecstate::DensityMatrix<std::complex<double>, double>& dm,
3332
const psi::Psi<std::complex<double>>* psi,
3433
const Parallel_Orbitals& pv,
3534
const elecstate::ElecState* pelec,
@@ -42,51 +41,6 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
4241
ModuleBase::timer::tick("Force_LCAO", "cal_fedm");
4342

4443
const int nspin = GlobalV::NSPIN;
45-
const int nbands = GlobalV::NBANDS;
46-
47-
// construct a DensityMatrix object
48-
elecstate::DensityMatrix<std::complex<double>, double> edm(kv, &pv, nspin);
49-
50-
//--------------------------------------------
51-
// calculate the energy density matrix here.
52-
//--------------------------------------------
53-
54-
ModuleBase::matrix wg_ekb;
55-
wg_ekb.create(kv->get_nks(), nbands);
56-
ModuleBase::Memory::record("Force::wg_ekb", sizeof(double) * kv->get_nks() * nbands);
57-
#ifdef _OPENMP
58-
#pragma omp parallel for collapse(2) schedule(static, 1024)
59-
#endif
60-
for (int ik = 0; ik < kv->get_nks(); ik++)
61-
{
62-
for (int ib = 0; ib < nbands; ib++)
63-
{
64-
wg_ekb(ik, ib) = pelec->wg(ik, ib) * pelec->ekb(ik, ib);
65-
}
66-
}
67-
std::vector<ModuleBase::ComplexMatrix> edm_k(kv->get_nks());
68-
69-
// use the original formula (Hamiltonian matrix) to calculate energy density matrix
70-
if (dm->EDMK.size())
71-
{
72-
#ifdef _OPENMP
73-
#pragma omp parallel for schedule(static)
74-
#endif
75-
for (int ik = 0; ik < kv->get_nks(); ++ik)
76-
{
77-
edm.set_DMK_pointer(ik, dm->EDMK[ik].c);
78-
}
79-
}
80-
else
81-
{
82-
// cal_dm_psi
83-
elecstate::cal_dm_psi(edm.get_paraV_pointer(), wg_ekb, psi[0], edm);
84-
}
85-
86-
// cal_dm_2d
87-
edm.init_DMR(*ra, &ucell);
88-
edm.cal_DMR();
89-
edm.sum_DMR_spin();
9044
//--------------------------------------------
9145
// summation \sum_{i,j} E(i,j)*dS(i,j)
9246
// BEGIN CALCULATION OF FORCE OF EACH ATOM
@@ -146,7 +100,7 @@ void Force_LCAO<std::complex<double>>::cal_fedm(const bool isforce,
146100
double Ry = ra->info[iat][cb][1];
147101
double Rz = ra->info[iat][cb][2];
148102
// get BaseMatrix
149-
hamilt::BaseMatrix<double>* tmp_matrix = edm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
103+
hamilt::BaseMatrix<double>* tmp_matrix = dm.get_DMR_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz);
150104
if (tmp_matrix == nullptr)
151105
{
152106
continue;

0 commit comments

Comments
 (0)