Skip to content

Commit feeabe1

Browse files
QianruipkuQianruipku
committed
Merge branch 'SDFT' of https://github.com/deepmodeling/abacus-develop into SDFT
2 parents f99da37 + 9c412e6 commit feeabe1

File tree

88 files changed

+4917
-607
lines changed

Some content is hidden

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

88 files changed

+4917
-607
lines changed

‎source/Makefile‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@ VPATH=./src_global\
1414
:./module_cell\
1515
:./module_base\
1616
:./module_md\
17+
:./module_surchem\
1718
:./module_deepks\
1819
:./module_xc\
1920
:./module_esolver\
@@ -69,6 +70,7 @@ FP_OBJS_0=$(OBJS_FIRST_PRINCIPLES)\
6970
$(OBJS_PARALLEL)\
7071
$(OBJS_TOOLS) \
7172
$(OBJS_COMMON)\
73+
$(OBJS_SURCHEM)\
7274
main.o\
7375

7476
FP_OBJS=$(patsubst %.o, ${OBJ_DIR}/%.o, ${FP_OBJS_0})

‎source/Makefile.Objects‎

Lines changed: 9 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -63,15 +63,6 @@ sto_iter.o\
6363
sto_hchi.o\
6464
sto_che.o\
6565

66-
OBJS_SURCHEM=cal_epsilon.o\
67-
cal_pseudo.o\
68-
cal_totn.o\
69-
cal_vcav.o\
70-
cal_vel.o\
71-
corrected_energy.o\
72-
H_correction_pw.o\
73-
minimize_cg.o
74-
7566
OBJS_TOOLS=complexarray.o\
7667
complexmatrix.o \
7768
realarray.o \
@@ -261,7 +252,6 @@ OBJS_FIRST_PRINCIPLES=$(OBJS_MAIN)\
261252
$(OBJS_PW)\
262253
$(OBJS_LCAO)\
263254
$(OBJS_ESOLVER)\
264-
$(OBJS_SURCHEM)\
265255
$(OBJS_XC)\
266256
charge.o \
267257
charge_mixing.o \
@@ -378,6 +368,15 @@ pzt2s.o\
378368
pdtrsm.o\
379369
pzhtrsm.o\
380370

371+
OBJS_SURCHEM=cal_epsilon.o\
372+
cal_pseudo.o\
373+
cal_totn.o\
374+
cal_vcav.o\
375+
cal_vel.o\
376+
corrected_energy.o\
377+
H_correction_pw.o\
378+
minimize_cg.o\
379+
381380
PDIAG_MR_0=dcopy.o\
382381
dlae2.o\
383382
dlaebz.o\

‎source/module_elecstate/CMakeLists.txt‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,8 @@
11
add_library(
22
elecstate
33
OBJECT
4+
elecstate.cpp
45
elecstate_pw.cpp
6+
elecstate_lcao.cpp
7+
dm2d_to_grid.cpp
58
)

‎source/module_elecstate/dm2d_to_grid.cpp‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
#include "module_base/timer.h"
2+
#include "module_psi/psi.h"
3+
#include "src_lcao/local_orbital_charge.h"
4+
5+
/// transformation from 2d block density matrix to grid points for one k point
6+
void Local_Orbital_Charge::dm2dToGrid(const psi::Psi<double>& dm2d, double** dm_grid)
7+
{
8+
ModuleBase::timer::tick("Local_Orbital_Charge", "dm_2dTOgrid");
9+
10+
#ifdef __MPI
11+
// put data from dm_gamma[ik] to sender index
12+
int nNONZERO = 0;
13+
for (int i = 0; i < this->sender_size; ++i)
14+
{
15+
const int idx = this->sender_2D_index[i];
16+
const int icol = idx % GlobalV::NLOCAL;
17+
const int irow = (idx - icol) / GlobalV::NLOCAL;
18+
// sender_buffer[i]=wfc_dm_2d.dm_gamma(ik, irow, icol);
19+
this->sender_buffer[i] = dm2d(icol, irow); // sender_buffer is clomun major,
20+
// so the row and column index should be switched
21+
if (this->sender_buffer[i] != 0)
22+
++nNONZERO;
23+
}
24+
25+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "number of non-zero elements in sender_buffer", nNONZERO);
26+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "sender_size", this->sender_size);
27+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "last sender_buffer", this->sender_buffer[this->sender_size - 1]);
28+
29+
// transform data via MPI_Alltoallv
30+
MPI_Alltoallv(this->sender_buffer,
31+
this->sender_size_process,
32+
this->sender_displacement_process,
33+
MPI_DOUBLE,
34+
this->receiver_buffer,
35+
this->receiver_size_process,
36+
this->receiver_displacement_process,
37+
MPI_DOUBLE,
38+
this->ParaV->comm_2D);
39+
40+
// put data from receiver buffer to dm_grid[ik]
41+
nNONZERO = 0;
42+
for (int i = 0; i < this->receiver_size; ++i)
43+
{
44+
const int idx = this->receiver_local_index[i];
45+
const int icol = idx % this->lgd_now;
46+
const int irow = (idx - icol) / this->lgd_now;
47+
dm_grid[irow][icol] = this->receiver_buffer[i];
48+
if (this->receiver_buffer[i] != 0)
49+
++nNONZERO;
50+
}
51+
52+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "number of non-zero elements in receiver_buffer", nNONZERO);
53+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "receiver_size", this->receiver_size);
54+
ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "last receiver_buffer", receiver_buffer[this->receiver_size - 1]);
55+
#else
56+
for (int irow = 0; irow < dm2d.get_nbasis(); ++irow)
57+
{
58+
for (int icol = 0; icol < dm2d.get_nbands(); ++icol)
59+
{
60+
dm_grid[irow][icol] = dm2d(icol, irow);
61+
}
62+
}
63+
64+
#endif
65+
66+
ModuleBase::timer::tick("Local_Orbital_Charge", "dm_2dTOgrid");
67+
return;
68+
}

‎source/module_elecstate/elecstate.cpp‎

Lines changed: 252 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,252 @@
1+
#include "elecstate.h"
2+
3+
#include "module_base/global_variable.h"
4+
#include "module_base/tool_title.h"
5+
#include "src_parallel/parallel_reduce.h"
6+
#include "src_pw/global.h"
7+
#include "src_pw/occupy.h"
8+
9+
namespace elecstate
10+
{
11+
12+
const double* ElecState::getRho(int spin) const
13+
{
14+
// hamilt::MatrixBlock<double> temp{&(this->charge->rho[spin][0]), 1, this->charge->nrxx}; //
15+
// this->chr->get_nspin(), this->chr->get_nrxx()};
16+
return &(this->charge->rho[spin][0]);
17+
}
18+
19+
void ElecState::calculate_weights(void)
20+
{
21+
ModuleBase::TITLE("ElecState", "calculate_weights");
22+
23+
// for test
24+
// std::cout << " gaussian_broadening = " << use_gaussian_broadening << std::endl;
25+
// std::cout << " tetrahedron_method = " << use_tetrahedron_method << std::endl;
26+
// std::cout << " fixed_occupations = " << fixed_occupations << std::endl;
27+
double** ekb_tmp = new double*[this->ekb.nr];
28+
for (int i = 0; i < this->ekb.nr; ++i)
29+
{
30+
ekb_tmp[i] = &(this->ekb(i, 0));
31+
}
32+
int nbands = this->ekb.nc;
33+
34+
if (GlobalV::KS_SOLVER == "selinv")
35+
{
36+
GlobalV::ofs_running << " Could not calculate occupation." << std::endl;
37+
return;
38+
}
39+
40+
if (!Occupy::use_gaussian_broadening && !Occupy::use_tetrahedron_method && !Occupy::fixed_occupations)
41+
{
42+
if (GlobalV::TWO_EFERMI)
43+
{
44+
Occupy::iweights(GlobalC::kv.nks,
45+
GlobalC::kv.wk,
46+
nbands,
47+
GlobalC::ucell.magnet.get_nelup(),
48+
ekb_tmp,
49+
GlobalC::en.ef_up,
50+
this->wg,
51+
0,
52+
GlobalC::kv.isk);
53+
Occupy::iweights(GlobalC::kv.nks,
54+
GlobalC::kv.wk,
55+
nbands,
56+
GlobalC::ucell.magnet.get_neldw(),
57+
ekb_tmp,
58+
GlobalC::en.ef_dw,
59+
this->wg,
60+
1,
61+
GlobalC::kv.isk);
62+
// ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16
63+
}
64+
else
65+
{
66+
// -1 means don't need to consider spin.
67+
Occupy::iweights(GlobalC::kv.nks,
68+
GlobalC::kv.wk,
69+
nbands,
70+
GlobalC::CHR.nelec,
71+
ekb_tmp,
72+
this->ef,
73+
this->wg,
74+
-1,
75+
GlobalC::kv.isk);
76+
}
77+
}
78+
else if (Occupy::use_tetrahedron_method)
79+
{
80+
ModuleBase::WARNING_QUIT("calculate_weights", "not implemented yet,coming soon!");
81+
// if(my_rank == 0)
82+
// {
83+
// tweights(GlobalC::kv.nkstot, nspin, nbands, GlobalC::CHR.nelec, ntetra,tetra, GlobalC::wf.et,
84+
// this->ef, this->wg);
85+
// }
86+
}
87+
else if (Occupy::use_gaussian_broadening)
88+
{
89+
if (GlobalV::TWO_EFERMI)
90+
{
91+
double demet_up = 0.0;
92+
double demet_dw = 0.0;
93+
Occupy::gweights(GlobalC::kv.nks,
94+
GlobalC::kv.wk,
95+
nbands,
96+
GlobalC::ucell.magnet.get_nelup(),
97+
Occupy::gaussian_parameter,
98+
Occupy::gaussian_type,
99+
ekb_tmp,
100+
GlobalC::en.ef_up,
101+
demet_up,
102+
this->wg,
103+
0,
104+
GlobalC::kv.isk);
105+
Occupy::gweights(GlobalC::kv.nks,
106+
GlobalC::kv.wk,
107+
nbands,
108+
GlobalC::ucell.magnet.get_neldw(),
109+
Occupy::gaussian_parameter,
110+
Occupy::gaussian_type,
111+
ekb_tmp,
112+
GlobalC::en.ef_dw,
113+
demet_dw,
114+
this->wg,
115+
1,
116+
GlobalC::kv.isk);
117+
GlobalC::en.demet = demet_up + demet_dw;
118+
}
119+
else
120+
{
121+
// -1 means is no related to spin.
122+
Occupy::gweights(GlobalC::kv.nks,
123+
GlobalC::kv.wk,
124+
nbands,
125+
GlobalC::CHR.nelec,
126+
Occupy::gaussian_parameter,
127+
Occupy::gaussian_type,
128+
ekb_tmp,
129+
this->ef,
130+
GlobalC::en.demet,
131+
this->wg,
132+
-1,
133+
GlobalC::kv.isk);
134+
}
135+
136+
// qianrui fix a bug on 2021-7-21
137+
Parallel_Reduce::reduce_double_allpool(GlobalC::en.demet);
138+
}
139+
else if (Occupy::fixed_occupations)
140+
{
141+
// fix occupations need nelup and neldw.
142+
// mohan add 2011-04-03
143+
this->ef = -1.0e+20;
144+
for (int ik = 0; ik < GlobalC::kv.nks; ik++)
145+
{
146+
for (int ibnd = 0; ibnd < nbands; ibnd++)
147+
{
148+
if (this->wg(ik, ibnd) > 0.0)
149+
{
150+
this->ef = std::max(this->ef, ekb_tmp[ik][ibnd]);
151+
}
152+
}
153+
}
154+
}
155+
156+
if (GlobalV::TWO_EFERMI)
157+
{
158+
Parallel_Reduce::gather_max_double_all(GlobalC::en.ef_up);
159+
Parallel_Reduce::gather_max_double_all(GlobalC::en.ef_dw);
160+
}
161+
else
162+
{
163+
double ebotom = ekb_tmp[0][0];
164+
double etop = ekb_tmp[0][0];
165+
for (int ik = 0; ik < GlobalC::kv.nks; ik++)
166+
{
167+
for (int ib = 0; ib < nbands; ib++)
168+
{
169+
ebotom = min(ebotom, ekb_tmp[ik][ib]);
170+
etop = max(etop, ekb_tmp[ik][ib]);
171+
}
172+
}
173+
174+
// parallel
175+
Parallel_Reduce::gather_max_double_all(this->ef);
176+
Parallel_Reduce::gather_max_double_all(etop);
177+
Parallel_Reduce::gather_min_double_all(ebotom);
178+
179+
// not parallel yet!
180+
// OUT(GlobalV::ofs_running,"Top Energy (eV)", etop * ModuleBase::Ry_to_eV);
181+
// OUT(GlobalV::ofs_running,"Fermi Energy (eV)", this->ef * ModuleBase::Ry_to_eV);
182+
// OUT(GlobalV::ofs_running,"Bottom Energy (eV)", ebotom * ModuleBase::Ry_to_eV);
183+
// OUT(GlobalV::ofs_running,"Range Energy (eV)", etop-ebotom * ModuleBase::Ry_to_eV);
184+
}
185+
186+
delete[] ekb_tmp;
187+
188+
return;
189+
}
190+
191+
double ElecState::eBandK(const int& ik)
192+
{
193+
ModuleBase::TITLE("ElecStatePW", "eBandK");
194+
double eband_k = 0.0;
195+
for (int ibnd = 0; ibnd < this->ekb.nc; ibnd++)
196+
{
197+
eband_k += this->ekb(ik, ibnd) * this->wg(ik, ibnd);
198+
}
199+
return eband_k;
200+
}
201+
202+
void ElecState::print_band(const int& ik, const int& printe, const int& iter)
203+
{
204+
// check the band energy.
205+
bool wrong = false;
206+
int nbands = this->ekb.nc;
207+
for (int ib = 0; ib < nbands; ++ib)
208+
{
209+
if (abs(this->ekb(ik, ib)) > 1.0e10)
210+
{
211+
GlobalV::ofs_warning << " ik=" << ik + 1 << " ib=" << ib + 1 << " " << this->ekb(ik, ib) << " Ry"
212+
<< std::endl;
213+
wrong = true;
214+
}
215+
}
216+
if (wrong)
217+
{
218+
ModuleBase::WARNING_QUIT("Threshold_Elec::print_eigenvalue", "Eigenvalues are too large!");
219+
}
220+
221+
if (GlobalV::MY_RANK == 0)
222+
{
223+
// if( GlobalV::DIAGO_TYPE == "selinv" ) xiaohui modify 2013-09-02
224+
if (GlobalV::KS_SOLVER == "selinv") // xiaohui add 2013-09-02
225+
{
226+
GlobalV::ofs_running << " No eigenvalues are available for selected inversion methods." << std::endl;
227+
}
228+
else
229+
{
230+
if (printe > 0 && ((iter + 1) % printe == 0))
231+
{
232+
// NEW_PART("ENERGY BANDS (Rydberg), (eV)");
233+
GlobalV::ofs_running << std::setprecision(6);
234+
GlobalV::ofs_running << " Energy (eV) & Occupations for spin=" << GlobalV::CURRENT_SPIN + 1
235+
<< " K-point=" << ik + 1 << std::endl;
236+
GlobalV::ofs_running << std::setiosflags(ios::showpoint);
237+
for (int ib = 0; ib < nbands; ib++)
238+
{
239+
GlobalV::ofs_running << " " << std::setw(6) << ib + 1 << std::setw(15)
240+
<< this->ekb(ik, ib) * ModuleBase::Ry_to_eV;
241+
// for the first electron iteration, we don't have the energy
242+
// spectrum, so we can't get the occupations.
243+
GlobalV::ofs_running << std::setw(15) << this->wg(ik, ib);
244+
GlobalV::ofs_running << std::endl;
245+
}
246+
}
247+
}
248+
}
249+
return;
250+
}
251+
252+
} // namespace elecstate

0 commit comments

Comments
 (0)