Skip to content

Commit 27869ba

Browse files
maki49maki49
committed
extract edm
1 parent 49414e6 commit 27869ba

File tree

8 files changed

+123
-164
lines changed

8 files changed

+123
-164
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
@@ -217,7 +217,9 @@ void Force_LCAO<double>::ftable(const bool isforce,
217217
this->allocate(pv, fsr, two_center_bundle, orb);
218218

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

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

‎source/module_hamilt_lcao/hamilt_lcaodft/FORCE_k.cpp‎

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

316316
// calculate the energy density matrix
317317
// and the force related to overlap matrix and energy density matrix.
318-
this->cal_fedm(isforce, isstress, fsr, ucell, dm, psi, pv, pelec, foverlap, soverlap, kv, ra);
318+
this->cal_fedm(isforce, isstress, fsr, ucell,
319+
this->cal_edm(pelec, *psi, *dm, *kv, pv, GlobalV::NSPIN, GlobalV::NBANDS, ucell, *ra),
320+
psi, pv, pelec, foverlap, soverlap, kv, ra);
319321

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

‎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 (PARAM.inp.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: 0 additions & 113 deletions
This file was deleted.

‎source/module_hamilt_lcao/hamilt_lcaodft/fedm_k.cpp‎

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

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

0 commit comments

Comments
 (0)