Skip to content

Commit 8996e78

Browse files
JGHan7JGHan7
committed
modified according to Mr. Chen's suggestion
1 parent e432d86 commit 8996e78

File tree

6 files changed

+87
-41
lines changed

6 files changed

+87
-41
lines changed

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 4 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -1102,19 +1102,10 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
11021102
}
11031103

11041104
// 7) rdmft, added by jghan, 2024-10-25
1105-
if ( PARAM.inp.rdmft == true )
1105+
if ( PARAM.inp.rdmft == true && get_init_value_rdmft )
11061106
{
1107-
ModuleBase::TITLE("RDMFT", "E & Egradient");
1108-
ModuleBase::timer::tick("RDMFT", "E & Egradient");
1109-
1110-
if( get_init_value_rdmft )
1111-
{
11121107
ModuleBase::matrix occ_number_ks(this->pelec->wg);
1113-
for(int ik=0; ik < occ_number_ks.nr; ++ik)
1114-
{
1115-
for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik];
1116-
}
1117-
1108+
for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik]; }
11181109
this->rdmft_solver.update_elec(occ_number_ks, *(this->psi));
11191110

11201111
//initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0.
@@ -1123,11 +1114,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int& iter)
11231114
dE_dWfc.zero_out();
11241115

11251116
double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc);
1126-
1127-
ModuleBase::timer::tick("RDMFT", "E & Egradient");
1128-
11291117
// break;
1130-
}
11311118
}
11321119

11331120
}
@@ -1225,15 +1212,10 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
12251212
#endif
12261213

12271214
/******** test RDMFT *********/
1228-
1229-
// rdmft, added by jghan, 2024-10-17
1230-
if ( PARAM.inp.rdmft == true )
1215+
if ( PARAM.inp.rdmft == true ) // rdmft, added by jghan, 2024-10-17
12311216
{
12321217
ModuleBase::matrix occ_number_ks(this->pelec->wg);
1233-
for(int ik=0; ik < occ_number_ks.nr; ++ik)
1234-
{
1235-
for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik];
1236-
}
1218+
for(int ik=0; ik < occ_number_ks.nr; ++ik) { for(int inb=0; inb < occ_number_ks.nc; ++inb) occ_number_ks(ik, inb) /= this->kv.wk[ik]; }
12371219
this->rdmft_solver.update_elec(occ_number_ks, *(this->psi));
12381220

12391221
//initialize the gradients of Etotal on occupation numbers and wfc, and set all elements to 0.
@@ -1243,7 +1225,6 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
12431225

12441226
double Etotal_RDMFT = this->rdmft_solver.run(dE_dOccNum, dE_dWfc);
12451227
}
1246-
12471228
/******** test RDMFT *********/
12481229

12491230

‎source/module_rdmft/rdmft.cpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -428,6 +428,9 @@ void RDMFT<TK, TR>::cal_Energy(const int cal_type)
428428
template <typename TK, typename TR>
429429
double RDMFT<TK, TR>::run(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>& E_gradient_wfc)
430430
{
431+
ModuleBase::TITLE("RDMFT", "E & Egradient");
432+
ModuleBase::timer::tick("RDMFT", "E & Egradient");
433+
431434
// this->cal_V_hartree();
432435
// this->cal_V_XC();
433436
this->cal_Hk_Hpsi();
@@ -439,8 +442,9 @@ double RDMFT<TK, TR>::run(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>& E
439442

440443
TK* pwfc = &occNum_HamiltWfc(0, 0, 0);
441444
TK* pwfc_out = &E_gradient_wfc(0, 0, 0);
442-
for(int i=0; i<wfc.size(); ++i) { pwfc_out[i] = pwfc[i];
443-
}
445+
for(int i=0; i<wfc.size(); ++i) { pwfc_out[i] = pwfc[i]; }
446+
447+
ModuleBase::timer::tick("RDMFT", "E & Egradient");
444448
// return E_RDMFT[3];
445449
return Etotal;
446450
}

‎source/module_rdmft/rdmft.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -199,8 +199,11 @@ class RDMFT
199199
// update occ_number for optimization algorithms that depend on Hamilton
200200
void update_occNumber(const ModuleBase::matrix& occ_number_in);
201201

202+
// update occ_number for optimization algorithms that depend on Hamilton
203+
void update_wg(const ModuleBase::matrix& wg_in);
204+
202205
// do all calculation after update occNum&wfc, get Etotal and the gradient of energy with respect to the occNum&wfc
203-
double run(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>&E_gradient_wfc);
206+
double run(ModuleBase::matrix& E_gradient_occNum, psi::Psi<TK>& E_gradient_wfc);
204207

205208

206209

‎source/module_rdmft/rdmft_tools.cpp‎

Lines changed: 25 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -110,8 +110,12 @@ void psiDotPsi<double>(const Parallel_Orbitals* ParaV, const Parallel_2D& para_E
110110

111111
// occNum_wfcHwfc = occNum*wfcHwfc + occNum_wfcHwfc
112112
// When symbol = 0, 1, 2, 3, 4, occNum = occNum, 0.5*occNum, g(occNum), 0.5*g(occNum), d_g(occNum)/d_occNum respectively. Default symbol=0.
113-
void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc, ModuleBase::matrix& occNum_wfcHwfc,
114-
int symbol, const std::string XC_func_rdmft, const double alpha)
113+
void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number,
114+
const ModuleBase::matrix& wfcHwfc,
115+
ModuleBase::matrix& occNum_wfcHwfc,
116+
int symbol,
117+
const std::string XC_func_rdmft,
118+
const double alpha)
115119
{
116120
for(int ir=0; ir<occ_number.nr; ++ ir)
117121
{
@@ -122,8 +126,15 @@ void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number, const ModuleBase::
122126

123127

124128
// for the gradient of Etotal with respect to occupation numbers
125-
void add_occNum(const K_Vectors& kv, const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
126-
const ModuleBase::matrix& wfcHwfc_dft_XC_in, const ModuleBase::matrix& wfcHwfc_exx_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha)
129+
void add_occNum(const K_Vectors& kv,
130+
const ModuleBase::matrix& occ_number,
131+
const ModuleBase::matrix& wfcHwfc_TV_in,
132+
const ModuleBase::matrix& wfcHwfc_hartree_in,
133+
const ModuleBase::matrix& wfcHwfc_dft_XC_in,
134+
const ModuleBase::matrix& wfcHwfc_exx_XC_in,
135+
ModuleBase::matrix& occNum_wfcHwfc,
136+
const std::string XC_func_rdmft,
137+
const double alpha)
127138
{
128139
occNum_wfcHwfc.zero_out();
129140
occNum_Mul_wfcHwfc(occ_number, wfcHwfc_exx_XC_in, occNum_wfcHwfc, 4, XC_func_rdmft, alpha);
@@ -141,8 +152,14 @@ void add_occNum(const K_Vectors& kv, const ModuleBase::matrix& occ_number, const
141152

142153

143154
// do wk*g(occNum)*wfcHwfc and add for TV, hartree, XC. This function just use once, so it can be replace and delete
144-
void add_wfcHwfc(const ModuleBase::matrix& wg, const ModuleBase::matrix& wk_fun_occNum, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
145-
const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha)
155+
void add_wfcHwfc(const ModuleBase::matrix& wg,
156+
const ModuleBase::matrix& wk_fun_occNum,
157+
const ModuleBase::matrix& wfcHwfc_TV_in,
158+
const ModuleBase::matrix& wfcHwfc_hartree_in,
159+
const ModuleBase::matrix& wfcHwfc_XC_in,
160+
ModuleBase::matrix& occNum_wfcHwfc,
161+
const std::string XC_func_rdmft,
162+
const double alpha)
146163
{
147164
occNum_wfcHwfc.zero_out();
148165
occNum_Mul_wfcHwfc(wg, wfcHwfc_TV_in, occNum_wfcHwfc);
@@ -216,7 +233,8 @@ void Veff_rdmft<TK, TR>::initialize_HR(const UnitCell* ucell_in,
216233
for (int iat1 = 0; iat1 < ucell_in->nat; iat1++)
217234
{
218235
auto tau1 = ucell_in->get_tau(iat1);
219-
int T1, I1;
236+
int T1 = 0;
237+
int I1 = 0;
220238
ucell_in->iat2iait(iat1, &I1, &T1);
221239
AdjacentAtomInfo adjs;
222240
GridD->Find_atom(*ucell_in, tau1, T1, I1, &adjs);

‎source/module_rdmft/rdmft_tools.h‎

Lines changed: 33 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -221,8 +221,16 @@ void occNum_MulPsi(const Parallel_Orbitals* ParaV, const ModuleBase::matrix& occ
221221

222222
// add psi with eta and g(eta)
223223
template <typename TK>
224-
void add_psi(const Parallel_Orbitals* ParaV, const K_Vectors* kv, const ModuleBase::matrix& occ_number, psi::Psi<TK>& psi_TV, psi::Psi<TK>& psi_hartree,
225-
psi::Psi<TK>& psi_dft_XC, psi::Psi<TK>& psi_exx_XC, psi::Psi<TK>& occNum_Hpsi, const std::string XC_func_rdmft = "hf", const double alpha = 1.0)
224+
void add_psi(const Parallel_Orbitals* ParaV,
225+
const K_Vectors* kv,
226+
const ModuleBase::matrix& occ_number,
227+
psi::Psi<TK>& psi_TV,
228+
psi::Psi<TK>& psi_hartree,
229+
psi::Psi<TK>& psi_dft_XC,
230+
psi::Psi<TK>& psi_exx_XC,
231+
psi::Psi<TK>& occNum_Hpsi,
232+
const std::string XC_func_rdmft = "hf",
233+
const double alpha = 1.0)
226234
{
227235
const int nk = psi_TV.get_nk();
228236
const int nbn_local = psi_TV.get_nbands();
@@ -255,21 +263,38 @@ void add_psi(const Parallel_Orbitals* ParaV, const K_Vectors* kv, const ModuleBa
255263

256264
// occNum_wfcHwfc = occNum*wfcHwfc + occNum_wfcHwfc
257265
// When symbol = 0, 1, 2, 3, 4, occNum = occNum, 0.5*occNum, g(occNum), 0.5*g(occNum), d_g(occNum)/d_occNum respectively. Default symbol=0.
258-
void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc, ModuleBase::matrix& occNum_wfcHwfc,
259-
int symbol = 0, const std::string XC_func_rdmft = "hf", const double alpha = 1.0);
266+
void occNum_Mul_wfcHwfc(const ModuleBase::matrix& occ_number,
267+
const ModuleBase::matrix& wfcHwfc,
268+
ModuleBase::matrix& occNum_wfcHwfc,
269+
int symbol = 0,
270+
const std::string XC_func_rdmft = "hf",
271+
const double alpha = 1.0);
260272

261273

262274
// Default symbol = 0 for the gradient of Etotal with respect to occupancy
263275
// symbol = 1 for the relevant calculation of Etotal
264-
void add_occNum(const K_Vectors& kv, const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
265-
const ModuleBase::matrix& wfcHwfc_dft_XC_in, const ModuleBase::matrix& wfcHwfc_exx_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft = "hf", const double alpha = 1.0);
276+
void add_occNum(const K_Vectors& kv,
277+
const ModuleBase::matrix& occ_number,
278+
const ModuleBase::matrix& wfcHwfc_TV_in,
279+
const ModuleBase::matrix& wfcHwfc_hartree_in,
280+
const ModuleBase::matrix& wfcHwfc_dft_XC_in,
281+
const ModuleBase::matrix& wfcHwfc_exx_XC_in,
282+
ModuleBase::matrix& occNum_wfcHwfc,
283+
const std::string XC_func_rdmft = "hf",
284+
const double alpha = 1.0);
266285

267286

268287
// // do wk*g(occNum)*wfcHwfc and add for TV, hartree, XC. This function just use once, so it can be replace and delete
269288
// void add_wfcHwfc(const std::vector<double>& wk_in, const ModuleBase::matrix& occ_number, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
270289
// const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha);
271-
void add_wfcHwfc(const ModuleBase::matrix& wg, const ModuleBase::matrix& wk_fun_occNum, const ModuleBase::matrix& wfcHwfc_TV_in, const ModuleBase::matrix& wfcHwfc_hartree_in,
272-
const ModuleBase::matrix& wfcHwfc_XC_in, ModuleBase::matrix& occNum_wfcHwfc, const std::string XC_func_rdmft, const double alpha);
290+
void add_wfcHwfc(const ModuleBase::matrix& wg,
291+
const ModuleBase::matrix& wk_fun_occNum,
292+
const ModuleBase::matrix& wfcHwfc_TV_in,
293+
const ModuleBase::matrix& wfcHwfc_hartree_in,
294+
const ModuleBase::matrix& wfcHwfc_XC_in,
295+
ModuleBase::matrix& occNum_wfcHwfc,
296+
const std::string XC_func_rdmft,
297+
const double alpha);
273298

274299

275300
//give certain occNum_wfcHwfc, get the corresponding energy

‎source/module_rdmft/update_state_rdmft.cpp‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -188,6 +188,21 @@ void RDMFT<TK, TR>::update_occNumber(const ModuleBase::matrix& occ_number_in)
188188
}
189189

190190

191+
template <typename TK, typename TR>
192+
void RDMFT<TK, TR>::update_wg(const ModuleBase::matrix& wg_in)
193+
{
194+
wg = (wg_in);
195+
occ_number = (wg);
196+
for(int ik=0; ik < wg.nr; ++ik)
197+
{
198+
for(int inb=0; inb < wg.nc; ++inb)
199+
{
200+
occ_number(ik, inb) /= kv->wk[ik];
201+
wk_fun_occNum(ik, inb) = kv->wk[ik] * occNum_func(occ_number(ik, inb), 2, XC_func_rdmft, alpha_power);
202+
}
203+
}
204+
}
205+
191206

192207
template class RDMFT<double, double>;
193208
template class RDMFT<std::complex<double>, double>;

0 commit comments

Comments
 (0)