abacus-mc/source/source_estate/module_pot/pot_xc_fdm.cpp at 02af32118737f97c6980bbbbcc15a0af2aa2bec0 · mohanchen/abacus-mc · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
//=======================
// AUTHOR : Peize Lin
// DATE :   2025-10-01
//=======================

#include "pot_xc_fdm.h"
#include "source_hamilt/module_xc/xc_functional.h"
#include "source_io/module_parameter/parameter.h"

namespace elecstate
{

PotXC_FDM::PotXC_FDM(
	const ModulePW::PW_Basis* rho_basis_in,
	const Charge*const chg_0_in,
	const UnitCell*const ucell)
	: chg_0(chg_0_in)
{
	this->rho_basis_ = rho_basis_in;
	this->dynamic_mode = true;
	this->fixed_mode = false;

	const double hybrid_alpha = XC_Functional::get_hybrid_alpha();
#ifdef __EXX
	const double hse_omega = XC_Functional::get_hse_omega();
#else
	const double hse_omega = 0.0;
#endif
	const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v_0
		= XC_Functional::v_xc(this->chg_0->nrxx, this->chg_0, ucell,
							  PARAM.inp.nspin,
							  PARAM.globalv.domag,
							  PARAM.globalv.domag_z,
							  hybrid_alpha,
							  hse_omega);
	this->v_xc_0 = std::get<2>(etxc_vtxc_v_0);
}

void PotXC_FDM::cal_v_eff(
	const Charge*const chg_1,
	const UnitCell*const ucell,
	ModuleBase::matrix& v_eff)
{
	ModuleBase::TITLE("PotXC_FDM", "cal_veff");
	ModuleBase::timer::start("PotXC_FDM", "cal_veff");

	assert(this->chg_0->nrxx == chg_1->nrxx);
	assert(this->chg_0->nspin == chg_1->nspin);

	Charge chg_01;
	chg_01.set_rhopw(chg_1->rhopw);
	chg_01.allocate(chg_1->nspin, chg_01.kin_density());

	for(int ir=0; ir<chg_01.nrxx; ++ir)
	{
		for(int is=0; is<chg_01.nspin; ++is)
			{ chg_01.rho[is][ir] = chg_0->rho[is][ir] + chg_1->rho[is][ir]; }
		chg_01.rho_core[ir] = chg_0->rho_core[ir] + chg_1->rho_core[ir];
	}

	const double hybrid_alpha = XC_Functional::get_hybrid_alpha();
#ifdef __EXX
	const double hse_omega = XC_Functional::get_hse_omega();
#else
	const double hse_omega = 0.0;
#endif
	const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v_01
		= XC_Functional::v_xc(chg_01.nrxx, &chg_01, ucell,
							  PARAM.inp.nspin,
							  PARAM.globalv.domag,
							  PARAM.globalv.domag_z,
							  hybrid_alpha,
							  hse_omega);
	const ModuleBase::matrix &v_xc_01 = std::get<2>(etxc_vtxc_v_01);

	v_eff += v_xc_01 - this->v_xc_0;

	ModuleBase::timer::end("PotXC_FDM", "cal_veff");
}

} // namespace elecstate