Skip to content

Commit 6866e07

Browse files
LKFEIYILKFEIYI
authored
fix
1 parent 72a73c7 commit 6866e07

3 files changed

Lines changed: 170 additions & 20 deletions

File tree

‎source/source_estate/module_pot/H_Hartree_pw.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -155,10 +155,10 @@ if (use_parabolic)
155155
// 2. 将离子修正能加到 Hartree 能量中
156156
// 这样总能量就正确了
157157
H_Hartree_pw::hartree_energy += e_ion;
158-
// 3. (可选) 日志输出
159-
if (GlobalV::RANK_IN_POOL == 0) {
160-
std::cout << " Parabolic Correction Applied. E_ion_corr = " << e_ion << " Ry" << std::endl;
161-
}
158+
// // 3. (可选) 日志输出
159+
// if (GlobalV::RANK_IN_POOL == 0) {
160+
// std::cout << " Parabolic Correction Applied. E_ion_corr = " << e_ion << " Ry" << std::endl;
161+
// }
162162

163163
// 4. 处理 nspin=2 的第二列势场 (能量只需加一次，所以这里不加)
164164
if (nspin == 2) {

‎source/source_hamilt/module_poisson/parabolic_correction.cpp‎

Lines changed: 151 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -84,16 +84,26 @@ double ParabolicCorrection::apply_correction(const UnitCell& cell,
8484
// 2. 计算电荷与偶极
8585
double net_charge = -PARAM.inp.nelec_delta;
8686
double total_dipole = calc_total_dipole(cell, rho_basis, rho_elec, nspin, dir, slab_center);
87+
88+
double total_quadrupole = calc_total_quadrupole(cell, rho_basis, rho_elec, nspin, dir, slab_center);
89+
double factor_common = (ModuleBase::FOUR_PI / omega) * ModuleBase::e2; // 4pi/Omega * e2
90+
double factor_quad = factor_common * 0.5;
91+
double madelung_term = - (ModuleBase::PI / 3.0) * net_charge / lat_vec * ModuleBase::e2;
92+
double quadrupole_term = - factor_quad * total_quadrupole;
93+
94+
double v_const = madelung_term + quadrupole_term;
8795

8896
// 缓存偶极矩供后续力修正使用
8997
this->last_total_dipole_ = total_dipole;
9098

9199
// 3. 计算离子修正能 (这是你想要加入的！)
92-
double e_ion_corr = calc_energy_correction(cell, dir, net_charge, total_dipole, vacuum_center);
100+
double e_ion_corr = calc_energy_correction(cell, dir, net_charge, total_dipole, vacuum_center,v_const);
101+
double e_elec_corr = 0.0;
93102

94103
// 4. 构造 1D 修正势并叠加到 v_hartree
95-
double factor = (ModuleBase::FOUR_PI / omega) * ModuleBase::e2;
104+
// double factor = (ModuleBase::FOUR_PI / omega) * ModuleBase::e2;
96105
int nrxx = rho_basis->nrxx;
106+
int nspin_eff = (nspin == 2) ? 2 : 1;
97107

98108
#ifdef _OPENMP
99109
#pragma omp parallel for
@@ -117,23 +127,38 @@ double ParabolicCorrection::apply_correction(const UnitCell& cell,
117127
double dist_bohr = dist_frac * lat_vec;
118128

119129
// 核心抛物线修正公式
120-
double v_corr = factor * ( -0.5 * net_charge * dist_bohr * dist_bohr
130+
double v_corr = factor_common * ( -0.5 * net_charge * dist_bohr * dist_bohr
121131
+ total_dipole * dist_bohr );
122132

123133
v_hartree[ir] += v_corr;
134+
double rho_val = 0.0;
135+
for(int is=0; is<nspin_eff; ++is) rho_val += rho_elec[is][ir];
136+
e_elec_corr += rho_val * v_corr;
137+
}
138+
Parallel_Reduce::reduce_pool(e_elec_corr);
139+
e_elec_corr *= (cell.omega / rho_basis->nxyz);
140+
double total_correction_energy = e_ion_corr - e_elec_corr;
141+
142+
if (GlobalV::RANK_IN_POOL == 0) {
143+
std::cout << "DEBUG PARABOLIC:" << std::endl;
144+
std::cout << " nelec_delta (Input): " << PARAM.inp.nelec_delta << std::endl;
145+
std::cout << " net_charge (Used): " << net_charge << std::endl;
146+
std::cout << " Total Dipole: " << total_dipole << std::endl;
147+
std::cout << " Total Quadrupole: " << total_quadrupole << std::endl; // 建议打印四极矩
148+
std::cout << " Madelung Term: " << madelung_term << std::endl; // 建议打印各项贡献
149+
std::cout << " Quadrupole Term: " << quadrupole_term << std::endl;
150+
std::cout << " Constant Shift: " << v_const << std::endl;
151+
152+
std::cout << " ----------------------------------------" << std::endl;
153+
std::cout << " Electronic Energy Corr: " << e_elec_corr << " Ry" << std::endl;
154+
std::cout << " Ionic Energy Corr: " << e_ion_corr << " Ry" << std::endl;
155+
std::cout << " TOTAL Parabolic Corr: " << total_correction_energy << " Ry" << std::endl;
156+
std::cout << " ----------------------------------------" << std::endl;
124157
}
125-
if (GlobalV::RANK_IN_POOL == 0) {
126-
std::cout << "DEBUG PARABOLIC:" << std::endl;
127-
std::cout << " nelec_delta (Input): " << PARAM.inp.nelec_delta << std::endl;
128-
std::cout << " net_charge (Used): " << net_charge << std::endl;
129-
//std::cout << " Calculated Q (Ion-Elec): " << calc_net_charge(cell, GlobalV::nelec) << std::endl;
130-
std::cout << " Total Dipole: " << total_dipole << std::endl;
131-
std::cout << " Slab Center: " << slab_center << std::endl;
132-
std::cout << " Factor: " << factor << std::endl;
133-
}
134158

135159
// 返回离子修正能，方便外部加到 Total Energy
136-
return e_ion_corr;
160+
// double total_ion_charge = net_charge - (-PARAM.inp.nelec_delta);
161+
return total_correction_energy;
137162
}
138163

139164
// ---------------------------------------------------------
@@ -228,7 +253,8 @@ double ParabolicCorrection::calc_energy_correction(const UnitCell& cell,
228253
int dir,
229254
double net_charge,
230255
double total_dipole,
231-
double vacuum_center)
256+
double vacuum_center,
257+
double v_const)
232258
{
233259
double lat_vec = 0.0;
234260
if (dir == 0) lat_vec = cell.a1.norm() * cell.lat0;
@@ -255,7 +281,7 @@ double ParabolicCorrection::calc_energy_correction(const UnitCell& cell,
255281

256282
// V_corr at atom position
257283
double v_at_atom = factor * ( -0.5 * net_charge * dist_bohr * dist_bohr
258-
+ total_dipole * dist_bohr );
284+
+ total_dipole * dist_bohr ) + v_const;
259285

260286
e_corr += Z * v_at_atom;
261287
}
@@ -311,4 +337,114 @@ void ParabolicCorrection::calc_force_correction(const UnitCell& cell,
311337
iat++;
312338
}
313339
}
340+
}
341+
342+
double ParabolicCorrection::calc_total_quadrupole(const UnitCell& cell,
343+
const ModulePW::PW_Basis* rho_basis,
344+
const double* const* rho_elec,
345+
int nspin,
346+
int dir,
347+
double center)
348+
{
349+
double lat_vec = 0.0;
350+
if(dir==0) lat_vec = cell.a1.norm() * cell.lat0;
351+
else if(dir==1) lat_vec = cell.a2.norm() * cell.lat0;
352+
else lat_vec = cell.a3.norm() * cell.lat0;
353+
354+
// Q_ion
355+
double q_ion = calc_ion_quadrupole(cell, dir, center);
356+
// 注意：离子位置已经是物理长度（Bohr），所以算出来直接是 Bohr^2
357+
// 如果 calc_ion_dipole 用的是分数坐标差，这里要注意单位转换。
358+
// 看原代码 calc_ion_dipole 中 dist 是分数坐标差，最后乘了 lat_vec。
359+
// 这里我们需要 (dist_frac * lat_vec)^2
360+
361+
// Q_elec
362+
double q_elec_frac = calc_elec_quadrupole(rho_basis, rho_elec, nspin, dir, center, cell.omega);
363+
// 同样需要乘以 lat_vec^2
364+
double q_elec = q_elec_frac * (lat_vec * lat_vec);
365+
// 假设辅助函数内部处理单位，或者在这里统一处理
366+
// 让我们看辅助函数实现细节：
367+
return q_ion - q_elec;
368+
}
369+
370+
double ParabolicCorrection::calc_ion_quadrupole(const UnitCell& cell, int dir, double center)
371+
{
372+
double q = 0.0;
373+
374+
double lat_vec = 0.0;
375+
if(dir==0) lat_vec = cell.a1.norm() * cell.lat0;
376+
else if(dir==1) lat_vec = cell.a2.norm() * cell.lat0;
377+
else lat_vec = cell.a3.norm() * cell.lat0;
378+
379+
for(int it=0; it<cell.ntype; ++it) {
380+
for(int ia=0; ia<cell.atoms[it].na; ++ia) {
381+
double pos = cell.atoms[it].taud[ia][dir];
382+
double dist_frac = pos - center;
383+
384+
// 最小镜像处理
385+
if(dist_frac > 0.5) dist_frac -= 1.0;
386+
if(dist_frac < -0.5) dist_frac += 1.0;
387+
388+
double dist_bohr = dist_frac * lat_vec;
389+
390+
// sum Z * r^2
391+
q += cell.atoms[it].ncpp.zv * (dist_bohr * dist_bohr);
392+
}
393+
}
394+
return q;
395+
}
396+
397+
double ParabolicCorrection::calc_elec_quadrupole(const ModulePW::PW_Basis* rho_basis,
398+
const double* const* rho_elec,
399+
int nspin,
400+
int dir,
401+
double center,
402+
double omega)
403+
{
404+
double q = 0.0;
405+
int nrxx = rho_basis->nrxx;
406+
int n_components = (nspin == 2) ? 2 : 1;
407+
408+
// 获取晶格长度
409+
// 这里的 omega 参数虽然传进来了，但我们需要长度来转换坐标
410+
// 简单起见，我们在外部计算好 lat_vec 比较麻烦，
411+
// 不如在这里根据 dir 重新算一下或者假设输入保证一致。
412+
// 为了严谨，建议像 calc_elec_dipole 一样，只算分数坐标部分的积分，最后在外面乘长度平方。
413+
// 但 calc_elec_dipole 是最后乘了 lat_vec。
414+
// 我们这里直接在循环里乘好吧，或者模仿原结构。
415+
416+
// 为了复用代码结构，我们计算 sum [ rho * (dist_frac)^2 ]
417+
// 最后再乘以 (lat_vec^2) * (omega / nxyz)
418+
419+
for (int ir = 0; ir < nrxx; ++ir)
420+
{
421+
int i = ir / (rho_basis->ny * rho_basis->nplane);
422+
int j = ir / rho_basis->nplane - i * rho_basis->ny;
423+
int k = ir % rho_basis->nplane + rho_basis->startz_current;
424+
425+
double pos = 0.0;
426+
if(dir==0) pos = (double)i / rho_basis->nx;
427+
else if(dir==1) pos = (double)j / rho_basis->ny;
428+
else pos = (double)k / rho_basis->nz;
429+
430+
double dist = pos - center;
431+
if(dist > 0.5) dist -= 1.0;
432+
if(dist < -0.5) dist += 1.0;
433+
434+
double rho_val = 0.0;
435+
for(int is=0; is<n_components; ++is) {
436+
rho_val += rho_elec[is][ir];
437+
}
438+
439+
q += rho_val * (dist * dist);
440+
}
441+
442+
Parallel_Reduce::reduce_pool(q);
443+
444+
// 积分体积元 dV = omega / nxyz
445+
q *= (omega / rho_basis->nxyz);
446+
447+
// 此时 q 是 sum rho * (dist_frac)^2 * dV
448+
// 我们需要在外部乘以 lat_vec^2 才能变成 Bohr^2 单位
449+
return q;
314450
}

‎source/source_hamilt/module_poisson/parabolic_correction.h‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,12 @@ class ParabolicCorrection
3838

3939
// 获取偶极矩的接口，方便力修正使用
4040
double get_last_total_dipole() const { return last_total_dipole_; }
41+
double calc_total_quadrupole(const UnitCell& cell,
42+
const ModulePW::PW_Basis* rho_basis,
43+
const double* const* rho_elec,
44+
int nspin,
45+
int dir,
46+
double center);
4147

4248
private:
4349
double calc_total_dipole(const UnitCell& cell,
@@ -60,10 +66,18 @@ class ParabolicCorrection
6066
int dir,
6167
double net_charge,
6268
double total_dipole,
63-
double vacuum_center);
69+
double vacuum_center,
70+
double v_const);
6471

6572
// 缓存一些中间变量供力修正使用，避免重复计算
6673
double last_total_dipole_ = 0.0;
74+
double calc_ion_quadrupole(const UnitCell& cell, int dir, double center);
75+
double calc_elec_quadrupole(const ModulePW::PW_Basis* rho_basis,
76+
const double* const* rho_elec,
77+
int nspin,
78+
int dir,
79+
double center,
80+
double omega);
6781
};
6882

6983
#endif

0 commit comments

Comments
 (0)