Skip to content

Commit f7e3770

Browse files
LKFEIYILKFEIYI
authored
parabolic
1 parent f5c5fb4 commit f7e3770

6 files changed

Lines changed: 276 additions & 2 deletions

File tree

‎source/source_estate/module_pot/H_Hartree_pw.cpp‎

Lines changed: 36 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
#include "source_base/timer.h"
66
#include "source_base/parallel_reduce.h"
77
#include "source_hamilt/module_poisson/mt_poisson.h"
8+
#include "source_hamilt/module_poisson/parabolic_correction.h"
89

910
namespace elecstate
1011
{
@@ -48,6 +49,7 @@ ModuleBase::matrix H_Hartree_pw::v_hartree(const UnitCell &cell,
4849
const int ig0 = rho_basis->ig_gge0;
4950

5051
bool use_mt = (PARAM.inp.dim_corr == "mt");
52+
bool use_parabolic = (PARAM.inp.dim_corr == "parabolic");
5153
double tpiba = cell.tpiba;
5254
double tpiba2 = cell.tpiba2;
5355
int dir = PARAM.inp.dim_corr_dir;; // 0=x, 1=y, 2=z
@@ -114,7 +116,40 @@ ModuleBase::matrix H_Hartree_pw::v_hartree(const UnitCell &cell,
114116
for (int ir = 0; ir < rho_basis->nrxx; ir++)
115117
v(is, ir) = Porter[ir].real();
116118
}
117-
119+
if (do_parabolic)
120+
if (use_parabolic)
121+
{
122+
ParabolicCorrection pc;
123+
// int dir = PARAM.inp.dim_corr_dir;
124+
125+
// 1. 修正 Spin 1 (Up) 或 Total
126+
// 注意：这里传入完整的 rho 和 nspin，模块内部会自动处理密度求和
127+
pc.apply_correction(
128+
GlobalC::unitcell,
129+
rho_basis,
130+
&v(0, 0), // 修改第一列势场
131+
rho, // 传入完整的 rho 指针 (const double* const*)
132+
nspin, // 传入 nspin
133+
GlobalV::nelec,
134+
dir
135+
);
136+
137+
// 2. 如果是 Spin 2 (Down)，必须施加同样的修正
138+
// 因为 Hartree 势对 Up 和 Down 电子是一样的
139+
if (nspin == 2) {
140+
pc.apply_correction(
141+
GlobalC::unitcell,
142+
rho_basis,
143+
&v(1, 0), // 修改第二列势场
144+
rho, // 依然传入同样的 rho，计算出的偶极矩是一样的
145+
nspin,
146+
GlobalV::nelec,
147+
dir
148+
);
149+
}
150+
}
151+
// 对于 nspin=4，通常只处理 v(0, :)，因为它是 H_scalar
152+
}
118153
ModuleBase::timer::tick("H_Hartree_pw", "v_hartree");
119154
return v;
120155
} // end subroutine v_h

‎source/source_hamilt/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ list(APPEND objects
88
module_ewald/H_Ewald_pw.cpp
99
module_ewald/dnrm2.cpp
1010
module_poisson/mt_poisson.cpp
11+
module_poisson/parabolic_correction.cpp
1112
)
1213

1314
add_library(

‎source/source_hamilt/module_ewald/H_Ewald_pw.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -172,7 +172,7 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
172172
ewaldg += fact * std::abs(rhon) * std::abs(rhon)
173173
* exp(- rho_basis->gg[ig] * cell.tpiba2 / alpha / 4.0 ) / rho_basis->gg[ig] / cell.tpiba2;
174174
*/
175-
}
175+
176176

177177
ewaldg += fact * std::abs(rhon) * std::abs(rhon)
178178
* exp(- rho_basis->gg[ig] * cell.tpiba2 / alpha / 4.0 ) / rho_basis->gg[ig] / cell.tpiba2;

‎source/source_hamilt/module_poisson/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33

44
set( sources
55
mt_poisson.cpp
6+
parabolic.cpp
67
)
78

89
add_library(hamilt_poisson OBJECT ${sources})

‎source/source_hamilt/module_poisson/parabolic_correction.cpp‎

Lines changed: 217 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,217 @@
1+
#include "parabolic_correction.h"
2+
#include "source_base/constants.h" // ModuleBase::e2, PI, FOUR_PI
3+
#include "source_base/parallel_reduce.h" // Parallel_Reduce
4+
#include "source_basis/module_pw/pw_basis.h"
5+
#include "source_cell/unitcell.h"
6+
7+
#include <cmath>
8+
#include <iostream>
9+
#include <algorithm> // for std::sort
10+
11+
// ---------------------------------------------------------
12+
// 1. 自动寻找真空层中心 (逻辑完全借鉴 Efield::autoset)
13+
// ---------------------------------------------------------
14+
double ParabolicCorrection::find_vacuum_center(const UnitCell& cell, int dir)
15+
{
16+
std::vector<double> pos;
17+
for (int it = 0; it < cell.ntype; ++it)
18+
{
19+
for (int ia = 0; ia < cell.atoms[it].na; ++ia)
20+
{
21+
double p = cell.atoms[it].taud[ia][dir];
22+
p -= std::floor(p); // 归一化到 [0, 1)
23+
pos.push_back(p);
24+
}
25+
}
26+
27+
std::sort(pos.begin(), pos.end());
28+
29+
double max_gap = 0.0;
30+
double center = 0.5;
31+
32+
// 内部间隙
33+
for (size_t i = 1; i < pos.size(); i++)
34+
{
35+
double gap = pos[i] - pos[i - 1];
36+
if (gap > max_gap)
37+
{
38+
max_gap = gap;
39+
center = (pos[i] + pos[i - 1]) / 2.0;
40+
}
41+
}
42+
43+
// 跨边界间隙
44+
if (!pos.empty()) {
45+
double tail_gap = pos[0] + 1.0 - pos.back();
46+
if (tail_gap > max_gap)
47+
{
48+
max_gap = tail_gap;
49+
center = (pos[0] + pos.back() + 1.0) / 2.0;
50+
if (center >= 1.0) center -= 1.0;
51+
}
52+
}
53+
54+
return center;
55+
}
56+
57+
ParabolicCorrection::ParabolicCorrection() {}
58+
ParabolicCorrection::~ParabolicCorrection() {}
59+
60+
// ---------------------------------------------------------
61+
// 2. 应用校正的主函数
62+
// ---------------------------------------------------------
63+
void ParabolicCorrection::apply_correction(const UnitCell& cell,
64+
const ModulePW::PW_Basis* rho_basis,
65+
double* v_hartree,
66+
const double* const* rho_elec, // 【修正1】指针的指针
67+
int nspin, // 【修正2】添加 nspin 参数
68+
double nelec,
69+
int dir)
70+
{
71+
// 基础几何参数
72+
double omega = cell.omega;
73+
double lat_vec = 0.0;
74+
double area = 0.0;
75+
76+
if (dir == 0) {
77+
lat_vec = cell.a1.norm() * cell.lat0;
78+
} else if (dir == 1) {
79+
lat_vec = cell.a2.norm() * cell.lat0;
80+
} else { // dir == 2
81+
lat_vec = cell.a3.norm() * cell.lat0;
82+
}
83+
area = omega / lat_vec;
84+
85+
// 1. 自动寻找真空层中心
86+
double vacuum_center = find_vacuum_center(cell, dir);
87+
88+
// 2. 定义 Slab 的几何中心 (真空中心的对面)
89+
double slab_center = vacuum_center + 0.5;
90+
if(slab_center >= 1.0) slab_center -= 1.0;
91+
92+
// 3. 计算净电荷 (Q_ion - Q_elec)
93+
double net_charge = calc_net_charge(cell, nelec);
94+
95+
// 4. 计算总偶极矩 (传递 slab_center 作为参考原点)
96+
double total_dipole = calc_total_dipole(cell, rho_basis, rho_elec, nspin, dir, slab_center);
97+
98+
// 5. 构造 1D 修正势
99+
// 系数 factor = 4pi / Area * e^2
100+
double factor = (ModuleBase::FOUR_PI / area) * ModuleBase::e2;
101+
102+
int nrxx = rho_basis->nrxx;
103+
104+
for (int ir = 0; ir < nrxx; ++ir)
105+
{
106+
int i = ir / (rho_basis->ny * rho_basis->nplane);
107+
int j = ir / rho_basis->nplane - i * rho_basis->ny;
108+
int k = ir % rho_basis->nplane + rho_basis->startz_current;
109+
110+
double coord_frac = 0.0;
111+
if (dir == 0) coord_frac = (double)i / rho_basis->nx;
112+
else if (dir == 1) coord_frac = (double)j / rho_basis->ny;
113+
else coord_frac = (double)k / rho_basis->nz;
114+
115+
// 计算到 Slab 中心的物理距离 (考虑 PBC)
116+
double dist_frac = coord_frac - slab_center;
117+
if (dist_frac > 0.5) dist_frac -= 1.0;
118+
if (dist_frac < -0.5) dist_frac += 1.0;
119+
120+
double dist_bohr = dist_frac * lat_vec;
121+
122+
// 核心抛物线修正公式
123+
double v_corr = factor * ( -0.5 * net_charge * dist_bohr * dist_bohr
124+
+ total_dipole * dist_bohr );
125+
126+
v_hartree[ir] += v_corr;
127+
}
128+
}
129+
130+
double ParabolicCorrection::calc_net_charge(const UnitCell& cell, double nelec)
131+
{
132+
double ion_charge = 0.0;
133+
for(int it=0; it<cell.ntype; ++it) {
134+
ion_charge += cell.atoms[it].na * cell.atoms[it].ncpp.zv;
135+
}
136+
return ion_charge - nelec;
137+
}
138+
139+
double ParabolicCorrection::calc_total_dipole(const UnitCell& cell,
140+
const ModulePW::PW_Basis* rho_basis,
141+
const double* const* rho_elec, // 【修正】类型匹配
142+
int nspin, // 【修正】传递 nspin
143+
int dir,
144+
double center)
145+
{
146+
double lat_vec = 0.0;
147+
if(dir==0) lat_vec = cell.a1.norm() * cell.lat0;
148+
else if(dir==1) lat_vec = cell.a2.norm() * cell.lat0;
149+
else lat_vec = cell.a3.norm() * cell.lat0;
150+
151+
double d_ion = calc_ion_dipole(cell, dir, center) * lat_vec;
152+
double d_elec = calc_elec_dipole(rho_basis, rho_elec, nspin, dir, center, cell.omega) * lat_vec;
153+
154+
return d_ion - d_elec;
155+
}
156+
157+
double ParabolicCorrection::calc_ion_dipole(const UnitCell& cell, int dir, double center)
158+
{
159+
double d = 0.0;
160+
for(int it=0; it<cell.ntype; ++it) {
161+
for(int ia=0; ia<cell.atoms[it].na; ++ia) {
162+
double pos = cell.atoms[it].taud[ia][dir];
163+
double dist = pos - center;
164+
if(dist > 0.5) dist -= 1.0;
165+
if(dist < -0.5) dist += 1.0;
166+
167+
d += cell.atoms[it].ncpp.zv * dist;
168+
}
169+
}
170+
return d;
171+
}
172+
173+
// ---------------------------------------------------------
174+
// 3. 计算电子偶极矩 (包含多自旋求和)
175+
// ---------------------------------------------------------
176+
double ParabolicCorrection::calc_elec_dipole(const ModulePW::PW_Basis* rho_basis,
177+
const double* const* rho_elec, // 【修正】指针的指针
178+
int nspin, // 【修正】接收 nspin
179+
int dir,
180+
double center,
181+
double omega)
182+
{
183+
double d = 0.0;
184+
int nrxx = rho_basis->nrxx;
185+
186+
// 【修正】Runtime 下，rho[0]=Up, rho[1]=Down，需要求和
187+
int n_components = (nspin == 2) ? 2 : 1;
188+
189+
for (int ir = 0; ir < nrxx; ++ir)
190+
{
191+
int i = ir / (rho_basis->ny * rho_basis->nplane);
192+
int j = ir / rho_basis->nplane - i * rho_basis->ny;
193+
int k = ir % rho_basis->nplane + rho_basis->startz_current;
194+
195+
double pos = 0.0;
196+
if(dir==0) pos = (double)i / rho_basis->nx;
197+
else if(dir==1) pos = (double)j / rho_basis->ny;
198+
else pos = (double)k / rho_basis->nz;
199+
200+
double dist = pos - center;
201+
if(dist > 0.5) dist -= 1.0;
202+
if(dist < -0.5) dist += 1.0;
203+
204+
double rho_val = 0.0;
205+
for(int is=0; is<n_components; ++is) {
206+
rho_val += rho_elec[is][ir]; // 【修正】正确访问二维数组
207+
}
208+
209+
d += rho_val * dist;
210+
} // 【修正】补回了丢失的括号
211+
212+
Parallel_Reduce::reduce_pool(d);
213+
214+
d *= (omega / rho_basis->nxyz);
215+
216+
return d;
217+
}

‎source/source_hamilt/module_poisson/parabolic_correction.h‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
// 修改前：const double* rho_elec
2+
// 修改后：const double* const* rho_elec, int nspin
3+
void apply_correction(const UnitCell& cell,
4+
const ModulePW::PW_Basis* rho_basis,
5+
double* v_hartree,
6+
const double* const* rho_elec, // 传入指针的指针
7+
int nspin, // 传入自旋数
8+
double nelec,
9+
int correction_dir);
10+
11+
private:
12+
// ... 其他函数 ...
13+
14+
// 修改后
15+
double calc_elec_dipole(const ModulePW::PW_Basis* rho_basis,
16+
const double* const* rho_elec, // 指针的指针
17+
int nspin, // 自旋
18+
int dir,
19+
double center_coord,
20+
double omega);

0 commit comments

Comments
 (0)