Skip to content

Commit 3e2e7b4

Browse files
LKFEIYILKFEIYI
authored
mt_correction
1 parent 24ff385 commit 3e2e7b4

7 files changed

Lines changed: 189 additions & 1 deletion

File tree

‎source/source_estate/module_pot/H_Hartree_pw.cpp‎

Lines changed: 32 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,16 +45,47 @@ ModuleBase::matrix H_Hartree_pw::v_hartree(const UnitCell &cell,
4545

4646
std::vector<std::complex<double>> vh_g(rho_basis->npw);
4747
const int ig0 = rho_basis->ig_gge0;
48+
49+
bool use_mt = (PARAM.inp.dim_corr == "mt");
50+
double tpiba = cell.tpiba;
51+
double tpiba2 = cell.tpiba2;
52+
int dir = INPUT.mt_special_dimension; // 0=x, 1=y, 2=z
53+
double L = 0.0;
54+
if (dir == 0) L = cell.a1.norm() * cell.lat0; // X方向, 注意乘以 lat0 (Bohr)
55+
else if (dir == 1) L = cell.a2.norm() * cell.lat0; // Y方向
56+
else if (dir == 2) L = cell.a3.norm() * cell.lat0; // Z方向
57+
4858
#ifdef _OPENMP
4959
#pragma omp parallel for reduction(+:ehart)
5060
#endif
5161
for (int ig = 0; ig < rho_basis->npw; ig++)
5262
{
63+
double g2 = tpiba2 * rho_basis->gg[ig];
64+
5365
if (ig == ig0)
5466
{
67+
if (use_mt)
68+
{
69+
double screen_val0 = MTPoisson::get_screen_val_g0(L, PARAM.inp.mt_type);
70+
double fac0 = ModuleBase::e2 * screen_val0;
71+
double rho_g0_sq = (conj(Porter[ig]) * Porter[ig]).real();
72+
ehart += rho_g0_sq * fac0;
73+
vh_g[ig] = fac0 * Porter[ig];
74+
}
75+
else{
76+
vh_g[ig] = std::complex<double>(0.0, 0.0);
77+
}
5578
continue; // skip G=0
5679
}
57-
const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / (cell.tpiba2 * rho_basis->gg[ig]);
80+
double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / g2;
81+
if (use_mt)
82+
{
83+
ModuleBase::Vector3<double> g_vec = rho_basis->gcar[ig] * tpiba;
84+
85+
double screen_val = MTPoisson::get_screen_val(g2, g_vec, L, 0, PARAM.inp.mt_type, dir); //alpha is not used now
86+
87+
fac += ModuleBase::e2 * screen_val;
88+
}
5889
ehart += (conj(Porter[ig]) * Porter[ig]).real() * fac;
5990
vh_g[ig] = fac * Porter[ig];
6091

‎source/source_hamilt/CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,13 @@
11
add_subdirectory(module_vdw)
22
add_subdirectory(module_surchem)
33
add_subdirectory(module_xc)
4+
add_subdirectory(module_poisson)
45

56
list(APPEND objects
67
operator.cpp
78
module_ewald/H_Ewald_pw.cpp
89
module_ewald/dnrm2.cpp
10+
module_poisson/mt_poisson.cpp
911
)
1012

1113
add_library(

‎source/source_hamilt/module_poisson/CMakeLists.txt‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
# source_hamilt/module_poisson/CMakeLists.txt
2+
3+
4+
set( sources
5+
mt_poisson.cpp
6+
)
7+
8+
add_library(hamilt_poisson OBJECT ${sources})
9+
10+
11+
if(ENABLE_COVERAGE)
12+
add_coverage(hamilt_poisson)
13+
endif()

‎source/source_hamilt/module_poisson/mt_poisson.cpp‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
#include "mt_poisson.h"
2+
#include "module_base/constants.h" // 获取 PI, FOUR_PI
3+
#include <iostream>
4+
5+
// 如果需要报错功能
6+
#include "module_base/global_function.h"
7+
8+
double MTPoisson::get_screen_val_g0(double L, std::string type) {
9+
// 对应 CP2K 逻辑: IF (grid%have_g0) screen_function%array(1) = pi*zlength*zlength/2.0_dp
10+
if (type == "2D") {
11+
return ModuleBase::PI * L * L / 2.0;
12+
}
13+
// 0D 暂时略过，因为需要实空间积分或者 alpha 处理
14+
return 0.0;
15+
}
16+
17+
double MTPoisson::get_screen_val(double g2,
18+
const ModuleBase::Vector3<double>& g_vec,
19+
double L,
20+
double alpha,
21+
std::string type,
22+
int direction)
23+
{
24+
// 安全检查：G^2 非常小视为 G=0，返回0 (外部应该处理了 G=0 的情况)
25+
if (g2 < 1e-10) return 0.0;
26+
27+
if (type == "2D") {
28+
// 1. 确定非周期方向的分量 (G_special)
29+
double g_special = 0.0;
30+
switch (direction) {
31+
case 0: g_special = g_vec.x; break; // X 方向真空
32+
case 1: g_special = g_vec.y; break; // Y 方向真空 (默认)
33+
case 2: g_special = g_vec.z; break; // Z 方向真空
34+
default:
35+
// 默认回落到 Y
36+
g_special = g_vec.y;
37+
break;
38+
}
39+
40+
// 2. 调用核心公式
41+
return calculate_mt2d(g2, g_special, L);
42+
}
43+
else if (type == "0D") {
44+
// 0D 逻辑 (尚未完全移植，需要 erfc 等特殊处理)
45+
return 0.0;
46+
}
47+
48+
return 0.0;
49+
}
50+
51+
double MTPoisson::calculate_mt2d(double g2, double g_special, double L) {
52+
// 公式来源: CP2K mt_util.F (Martyna-Tuckerman, 1999)
53+
// V_screen(G) = - (4pi/G^2) * cos(G_z * L / 2) * exp(-G_xy * L / 2)
54+
// 这里 G_z 泛指非周期方向分量(g_special)，G_xy 泛指周期平面分量
55+
56+
// 1. 计算标准库仑项 (注意符号，这里只计算系数)
57+
// 最终公式是: fac = (4pi/G^2) + V_screen
58+
// 所以 V_screen = (4pi/G^2) * [ -cos(...) * exp(...) ]
59+
double coulomb_term = ModuleBase::FOUR_PI / g2;
60+
61+
// 2. 计算周期性平面内的 G 分量模长
62+
// G^2 = G_special^2 + G_planar^2 => G_planar = sqrt(G^2 - G_special^2)
63+
double g_planar = std::sqrt(std::abs(g2 - g_special * g_special));
64+
65+
// 3. 计算指数项和余弦项
66+
double half_L = L / 2.0;
67+
double cos_term = std::cos(g_special * half_L);
68+
double exp_term = std::exp(-g_planar * half_L);
69+
70+
// 4. 组合
71+
// 注意负号
72+
return -coulomb_term * cos_term * exp_term;
73+
}

‎source/source_hamilt/module_poisson/mt_poisson.h‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
#ifndef MT_POISSON_H
2+
#define MT_POISSON_H
3+
4+
#include "module_base/vector3.h"
5+
#include "module_base/constants.h"
6+
#include <string>
7+
#include <cmath>
8+
9+
// 假设 PARAM 或者 INPUT 结构体已经包含了 mt_slab_size 等参数
10+
// 这里我们尽量让接口保持纯净，只传递数值，不直接依赖全局变量
11+
// 这样方便单元测试
12+
13+
class MTPoisson {
14+
public:
15+
/**
16+
* @brief 计算 G=0 处的修正值 (Gamma点)
17+
* * @param L 真空层厚度 (slab size), 单位: Bohr
18+
* @param type MT类型 ("2D", "0D" 等)
19+
* @return double 修正势的值 (Hartree)
20+
*/
21+
static double get_screen_val_g0(double L, std::string type);
22+
23+
/**
24+
* @brief 计算 G != 0 处的屏蔽势 V_screen(G)
25+
* * @param g2 G向量模的平方 (单位: Bohr^-2)
26+
* @param g_vec G向量 (单位: Bohr^-1)
27+
* @param L 真空层厚度 (slab size), 单位: Bohr
28+
* @param alpha 平滑参数 (用于0D/1D), 单位: Bohr^-1
29+
* @param type MT类型 ("2D", "0D" 等)
30+
* @param direction 非周期方向 (0=X, 1=Y, 2=Z), 默认为 1 (Y)
31+
* @return double 修正值，需与 4pi/G^2 叠加
32+
*/
33+
static double get_screen_val(double g2,
34+
const ModuleBase::Vector3<double>& g_vec,
35+
double L,
36+
double alpha,
37+
std::string type,
38+
int direction = 1);
39+
40+
private:
41+
// 具体的 MT2D 实现逻辑
42+
static double calculate_mt2d(double g2, double g_special, double L);
43+
};
44+
45+
#endif // MT_POISSON_H

‎source/source_io/module_parameter/input_parameter.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -697,5 +697,10 @@ struct Input_para
697697
// ============== #Parameters (25.uncommon hardware) =================
698698
int dsp_count = 4; /// the count of dsp hardwares in one node
699699

700+
// ============== #Parameters (26.uncommon hardware) =================
701+
std::string dim_corr = 'mt';
702+
int dim_corr_dir = 1; // x 0, y 1, z 3 abacus could have the best performance when x/y
703+
std::string mt_type = '2D';
704+
700705
};
701706
#endif

‎source/source_io/read_input_item_model.cpp‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,25 @@ void ReadInput::item_model()
130130
this->add_item(item);
131131
}
132132

133+
// MT correction
134+
{
135+
Input_Item item("dim_corr");
136+
item.annotation = "how to correct dimession";
137+
read_sync_double(input.dim_corr);
138+
this->add_item(item);
139+
}
140+
{
141+
Input_Item item("dim_corr_dir");
142+
item.annotation = "the direction for correction";
143+
read_sync_double(input.mt_special_dimension);
144+
this->add_item(item);
145+
}
146+
{
147+
Input_Item item("mt_type");
148+
item.annotation = "MT correction type 2D/0D 0D is not implemented";
149+
read_sync_double(input.mt_type);
150+
this->add_item(item);
151+
}
133152
// vdW Correction
134153
{
135154
Input_Item item("vdw_method");

0 commit comments

Comments
 (0)