Skip to content

Commit 377b9aa

Browse files
maki49maki49
committed
fix the overwriting between xc kernel components
1 parent a6631a6 commit 377b9aa

File tree

5 files changed

+56
-91
lines changed

5 files changed

+56
-91
lines changed

‎source/module_lr/potentials/xc_kernel.cpp‎

Lines changed: 51 additions & 86 deletions
Original file line numberDiff line numberDiff line change
@@ -144,13 +144,13 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
144144
}
145145
// -----------------------------------------------------------------------------------
146146
//==================== XC Kernels (f_xc)=============================
147-
this->vrho_.resize(nspin * nrxx);
148-
this->v2rho2_.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11
147+
this->vrho_.resize(nspin * nrxx, 0.);
148+
this->v2rho2_.resize(((1 == nspin) ? 1 : 3) * nrxx, 0.);//(nrxx* ((1 == nspin) ? 1 : 3)): 00, 01, 11
149149
if (is_gga)
150150
{
151-
this->vsigma_.resize(((1 == nspin) ? 1 : 3) * nrxx);//(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12
152-
this->v2rhosigma_.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12
153-
this->v2sigma2_.resize(((1 == nspin) ? 1 : 6) * nrxx); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22
151+
this->vsigma_.resize(((1 == nspin) ? 1 : 3) * nrxx, 0.);//(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12
152+
this->v2rhosigma_.resize(((1 == nspin) ? 1 : 6) * nrxx, 0.); //(nrxx*): 2 for rho * 3 for sigma: 00, 01, 02, 10, 11, 12
153+
this->v2sigma2_.resize(((1 == nspin) ? 1 : 6) * nrxx, 0.); //(nrxx* ((1 == nspin) ? 1 : 6)): 00, 01, 02, 11, 12, 22
154154
}
155155
//MetaGGA ...
156156

@@ -163,102 +163,67 @@ void LR::KernelXC::f_xc_libxc(const int& nspin, const double& omega, const doubl
163163

164164
//cut off grho if not LDA (future subfunc)
165165

166+
// Libxc interfaces overwrite (instead of add onto) the output arrays, so we need temporary copies
167+
std::vector<double> vrho_tmp(this->vrho_.size());
168+
std::vector<double> v2rho2_tmp(this->v2rho2_.size());
169+
std::vector<double> vsigma_tmp(this->vsigma_.size());
170+
std::vector<double> v2rhosigma_tmp(this->v2rhosigma_.size());
171+
std::vector<double> v2sigma2_tmp(this->v2sigma2_.size());
166172
switch (func.info->family)
167173
{
168174
case XC_FAMILY_LDA:
169-
xc_lda_vxc(&func, nrxx, rho.data(), vrho_.data());
170-
xc_lda_fxc(&func, nrxx, rho.data(), v2rho2_.data());
175+
xc_lda_vxc(&func, nrxx, rho.data(), vrho_tmp.data());
176+
xc_lda_fxc(&func, nrxx, rho.data(), v2rho2_tmp.data());
171177
break;
172178
case XC_FAMILY_GGA:
173179
case XC_FAMILY_HYB_GGA:
174-
xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), vrho_.data(), vsigma_.data());
175-
xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2_.data(), v2rhosigma_.data(), v2sigma2_.data());
180+
xc_gga_vxc(&func, nrxx, rho.data(), sigma.data(), vrho_tmp.data(), vsigma_tmp.data());
181+
xc_gga_fxc(&func, nrxx, rho.data(), sigma.data(), v2rho2_tmp.data(), v2rhosigma_tmp.data(), v2sigma2_tmp.data());
176182
break;
177183
default:
178184
throw std::domain_error("func.info->family =" + std::to_string(func.info->family)
179185
+ " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
180186
break;
181187
}
182-
// some formulas for GGA
183-
if (func.info->family == XC_FAMILY_GGA || func.info->family == XC_FAMILY_HYB_GGA)
184-
{
185-
const std::vector<double>& v2r2 = this->v2rho2_;
186-
const std::vector<double>& v2rs = this->v2rhosigma_;
187-
const std::vector<double>& v2s2 = this->v2sigma2_;
188-
const std::vector<double>& vs = this->vsigma_;
189-
const double tpiba2 = tpiba * tpiba;
190-
191-
if (1 == nspin)
192-
{
193-
using V3 = ModuleBase::Vector3<double>;
194-
// 0. drho
195-
this->drho_gs_ = gradrho[0];
196-
// 1. $2f^{\rho\sigma}*\nabla\rho$
197-
this->v2rhosigma_2drho_.resize(nrxx);
198-
std::transform(gradrho[0].begin(), gradrho[0].end(), v2rs.begin(), this->v2rhosigma_2drho_.begin(),
199-
[](const V3& a, const V3& b) {return a * b * 2.; });
200-
// 2. $4f^{\sigma\sigma}*\nabla\rho$
201-
this->v2sigma2_4drho_.resize(nrxx);
202-
std::transform(sigma.begin(), sigma.end(), v2s2.begin(), v2sigma2_4drho_.begin(),
203-
[](const V3& a, const V3& b) {return a * b * 4.; });
204-
}
205-
// else if (2 == nspin) // wrong, to be fixed
206-
// {
207-
// // 1. $\nabla\cdot(f^{\rho\sigma}*\nabla\rho)$
208-
// std::vector<double> div_v2rhosigma_gdrho_r(3 * nrxx);
209-
// std::vector<ModuleBase::Vector3<double>> v2rhosigma_gdrho_r(3 * nrxx);
210-
// for (int ir = 0; ir < nrxx; ++ir)
211-
// {
212-
// v2rhosigma_gdrho_r[ir] = gradrho[0][ir] * v2rs.at(ir * 6) * 4.0
213-
// + gradrho[1][ir] * v2rs.at(ir * 6 + 1) * 2.0; //up-up
214-
// v2rhosigma_gdrho_r[nrxx + ir] = gradrho[0][ir] * (v2rs.at(ir * 6 + 3) * 2.0 + v2rs.at(ir * 6 + 1))
215-
// + gradrho[1][ir] * (v2rs.at(ir * 6 + 2) * 2.0 + v2rs.at(ir * 6 + 4)); //up-down
216-
// v2rhosigma_gdrho_r[2 * nrxx + ir] = gradrho[1][ir] * v2rs.at(ir * 6 + 5) * 4.0
217-
// + gradrho[0][ir] * v2rs.at(ir * 6 + 4) * 2.0; //down-down
218-
// }
219-
// for (int isig = 0;isig < 3;++isig)
220-
// XC_Functional::grad_dot(v2rhosigma_gdrho_r.data() + isig * nrxx, div_v2rhosigma_gdrho_r.data() + isig * nrxx, chg_gs.rhopw, tpiba);
221-
// // 2. $\nabla^2(f^{\sigma\sigma}*\sigma)$
222-
// std::vector<double> v2sigma2_sigma_r(3 * nrxx);
223-
// for (int ir = 0; ir < nrxx; ++ir)
224-
// {
225-
// v2sigma2_sigma_r[ir] = v2s2.at(ir * 6) * sigma[ir * 3] * 4.0
226-
// + v2s2.at(ir * 6 + 1) * sigma[ir * 3 + 1] * 4.0
227-
// + v2s2.at(ir * 6 + 3) * sigma[ir * 3 + 2]; //up-up
228-
// v2sigma2_sigma_r[nrxx + ir] = v2s2.at(ir * 6 + 1) * sigma[ir * 3] * 2.0
229-
// + v2s2.at(ir * 6 + 4) * sigma[ir * 3 + 2] * 2.0
230-
// + (v2s2.at(ir * 6 + 2) * 4.0 + v2s2.at(ir * 6 + 3)) * sigma[ir * 3 + 1]; //up-down
231-
// v2sigma2_sigma_r[2 * nrxx + ir] = v2s2.at(ir * 6 + 5) * sigma[ir * 3 + 2] * 4.0
232-
// + v2s2.at(ir * 6 + 4) * sigma[ir * 3 + 1] * 4.0
233-
// + v2s2.at(ir * 6 + 3) * sigma[ir * 3]; //down-down
234-
// }
235-
// for (int isig = 0;isig < 3;++isig)
236-
// LR_Util::lapl(v2sigma2_sigma_r.data() + isig * nrxx, v2sigma2_sigma_r.data() + isig * nrxx, *(chg_gs.rhopw), tpiba2);
237-
// // 3. $\nabla^2(v^\sigma)$
238-
// std::vector<double> lap_vsigma(3 * nrxx);
239-
// for (int ir = 0;ir < nrxx;++ir)
240-
// {
241-
// lap_vsigma[ir] = vs.at(ir * 3) * 2.0;
242-
// lap_vsigma[nrxx + ir] = vs.at(ir * 3 + 1);
243-
// lap_vsigma[2 * nrxx + ir] = vs.at(ir * 3 + 2) * 2.0;
244-
// }
245-
// for (int isig = 0;isig < 3;++isig)
246-
// LR_Util::lapl(lap_vsigma.data() + isig * nrxx, lap_vsigma.data() + isig * nrxx, *(chg_gs.rhopw), tpiba2);
247-
// // add to v2rho2_
248-
// BlasConnector::axpy(3 * nrxx, 1.0, v2r2.data(), 1, to_mul_rho_.data(), 1);
249-
// BlasConnector::axpy(3 * nrxx, -1.0, div_v2rhosigma_gdrho_r.data(), 1, to_mul_rho_.data(), 1);
250-
// BlasConnector::axpy(3 * nrxx, 1.0, v2sigma2_sigma_r.data(), 1, to_mul_rho_.data(), 1);
251-
// BlasConnector::axpy(3 * nrxx, 1.0, lap_vsigma.data(), 1, to_mul_rho_.data(), 1);
252-
// }
253-
else
254-
{
255-
throw std::domain_error("nspin =" + std::to_string(nspin)
256-
+ " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
257-
}
258-
}
188+
// add onto the total components
189+
std::transform(vrho_tmp.begin(), vrho_tmp.end(), this->vrho_.begin(), this->vrho_.begin(), std::plus<double>());
190+
std::transform(v2rho2_tmp.begin(), v2rho2_tmp.end(), this->v2rho2_.begin(), this->v2rho2_.begin(), std::plus<double>());
191+
std::transform(vsigma_tmp.begin(), vsigma_tmp.end(), this->vsigma_.begin(), this->vsigma_.begin(), std::plus<double>());
192+
std::transform(v2rhosigma_tmp.begin(), v2rhosigma_tmp.end(), this->v2rhosigma_.begin(), this->v2rhosigma_.begin(), std::plus<double>());
193+
std::transform(v2sigma2_tmp.begin(), v2sigma2_tmp.end(), this->v2sigma2_.begin(), this->v2sigma2_.begin(), std::plus<double>());
259194
} // end for( xc_func_type &func : funcs )
195+
260196
XC_Functional_Libxc::finish_func(funcs);
261197

198+
// some postprocess if there're GGA funcs in the list
199+
if (is_gga)
200+
{
201+
const std::vector<double>& v2r2 = this->v2rho2_;
202+
const std::vector<double>& v2rs = this->v2rhosigma_;
203+
const std::vector<double>& v2s2 = this->v2sigma2_;
204+
const std::vector<double>& vs = this->vsigma_;
205+
const double tpiba2 = tpiba * tpiba;
206+
207+
if (1 == nspin)
208+
{
209+
using V3 = ModuleBase::Vector3<double>;
210+
// 0. drho
211+
this->drho_gs_ = gradrho[0];
212+
// 1. $2f^{\rho\sigma}*\nabla\rho$
213+
this->v2rhosigma_2drho_.resize(nrxx);
214+
std::transform(gradrho[0].begin(), gradrho[0].end(), v2rs.begin(), this->v2rhosigma_2drho_.begin(),
215+
[](const V3& a, const V3& b) {return a * b * 2.; });
216+
// 2. $4f^{\sigma\sigma}*\nabla\rho$
217+
this->v2sigma2_4drho_.resize(nrxx);
218+
std::transform(sigma.begin(), sigma.end(), v2s2.begin(), v2sigma2_4drho_.begin(),
219+
[](const V3& a, const V3& b) {return a * b * 4.; });
220+
}
221+
else
222+
{
223+
throw std::domain_error("nspin =" + std::to_string(nspin)
224+
+ " unfinished in " + std::string(__FILE__) + " line " + std::to_string(__LINE__));
225+
}
226+
}
262227
if (1 == PARAM.inp.nspin || 2 == PARAM.inp.nspin) return;
263228
// else if (4 == PARAM.inp.nspin)
264229
else//NSPIN != 1,2,4 is not supported

‎tests/integrate/291_NO_KP_LR/result.ref‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
totexcitationenergyref 0.784274
1+
totexcitationenergyref 0.755653

‎tests/integrate/391_NO_GO_LR/result.ref‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
totexcitationenergyref 2.641190
1+
totexcitationenergyref 2.510666

‎tests/integrate/393_NO_GO_ULR_HF/INPUT‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -27,12 +27,12 @@ mixing_type pulay
2727
mixing_beta 0.4
2828
mixing_gg0 0
2929

30-
dft_functional lda
30+
dft_functional hf
3131
exx_real_number 1
3232
exx_ccp_rmesh_times 1
3333

3434
lr_nstates 3
35-
xc_kernel lda
35+
xc_kernel hf
3636
lr_solver lapack
3737
lr_thr 1e-2
3838
lr_unrestricted 1

‎tests/integrate/393_NO_GO_ULR_HF/result.ref‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1 +1 @@
1-
totexcitationenergyref 1.940703
1+
totexcitationenergyref 2.187535

0 commit comments

Comments
 (0)