ion_cor

Author: LKFEIYI

source/source_estate/module_pot/H_Hartree_pw.cpp

@@ -116,40 +116,39 @@ ModuleBase::matrix H_Hartree_pw::v_hartree(const UnitCell &cell,
             for (int ir = 0; ir < rho_basis->nrxx; ir++)
                 v(is, ir) = Porter[ir].real();
     }
-    if (do_parabolic) 
 if (use_parabolic) 
     {
         ParabolicCorrection pc;
-        // int dir = PARAM.inp.dim_corr_dir;
-
-        // 1. 修正 Spin 1 (Up) 或 Total
-        // 注意：这里传入完整的 rho 和 nspin，模块内部会自动处理密度求和
-        pc.apply_correction(
+        
+        // 1. 调用 apply_correction (同时获取离子修正能)
+        double e_ion = pc.apply_correction(
             GlobalC::unitcell,
             rho_basis,
-            &v(0, 0),       // 修改第一列势场
-            rho,            // 传入完整的 rho 指针 (const double* const*)
-            nspin,          // 传入 nspin
+            &v(0, 0),       // 势场指针
+            rho,            // 密度指针 (const double* const*)
+            nspin,          // 自旋
             GlobalV::nelec,
             dir
         );
 
-        // 2. 如果是 Spin 2 (Down)，必须施加同样的修正
-        // 因为 Hartree 势对 Up 和 Down 电子是一样的
+        // 2. 将离子修正能加到 Hartree 能量中
+        // 这样总能量就正确了
+        H_Hartree_pw::hartree_energy += e_ion;
+        
+        // 3. (可选) 日志输出
+        if (GlobalV::RANK_IN_POOL == 0) {
+            std::cout << " Parabolic Correction Applied. E_ion_corr = " << e_ion << " Ry" << std::endl;
+        }
+
+        // 4. 处理 nspin=2 的第二列势场 (能量只需加一次，所以这里不加)
         if (nspin == 2) {
              pc.apply_correction(
-                GlobalC::unitcell,
-                rho_basis,
-                &v(1, 0),   // 修改第二列势场
-                rho,        // 依然传入同样的 rho，计算出的偶极矩是一样的
-                nspin,
-                GlobalV::nelec,
-                dir
+                GlobalC::unitcell, rho_basis, &v(1, 0), 
+                rho, nspin, GlobalV::nelec, dir
             );
         }
     }
-        // 对于 nspin=4，通常只处理 v(0, :)，因为它是 H_scalar
-    }
+
     ModuleBase::timer::tick("H_Hartree_pw", "v_hartree");
     return v;
 } // end subroutine v_h

source/source_hamilt/module_ewald/H_Ewald_pw.cpp

@@ -168,14 +168,7 @@ double H_Ewald_pw::compute_ewald(const UnitCell& cell,
         // 修改累加公式：加上 term_mt
         ewaldg += fact * std::abs(rhon) * std::abs(rhon) * (term_std + term_mt);
 
-        /* 原代码参考：
-        ewaldg += fact * std::abs(rhon) * std::abs(rhon)
-                  * exp(- rho_basis->gg[ig] * cell.tpiba2 / alpha / 4.0 ) / rho_basis->gg[ig] / cell.tpiba2;
-        */
 
-
-        ewaldg += fact * std::abs(rhon) * std::abs(rhon)
-                  * exp(- rho_basis->gg[ig] * cell.tpiba2 / alpha / 4.0 ) / rho_basis->gg[ig] / cell.tpiba2;
     }
 
     ewaldg = ModuleBase::FOUR_PI / cell.omega * ewaldg;

source/source_hamilt/module_poisson/CMakeLists.txt

@@ -3,7 +3,7 @@
 
 set( sources
     mt_poisson.cpp
-    parabolic.cpp
+    parabolic_correction.cpp
 )
 
 add_library(hamilt_poisson OBJECT ${sources})

source/source_hamilt/module_poisson/parabolic_correction.cpp

@@ -1,15 +1,18 @@
 #include "parabolic_correction.h"
-#include "source_base/constants.h" // ModuleBase::e2, PI, FOUR_PI
-#include "source_base/parallel_reduce.h" // Parallel_Reduce
+#include "source_base/constants.h" 
+#include "source_base/parallel_reduce.h" 
 #include "source_basis/module_pw/pw_basis.h"
 #include "source_cell/unitcell.h"
 
 #include <cmath>
 #include <iostream>
-#include <algorithm> // for std::sort
+#include <algorithm> 
+
+ParabolicCorrection::ParabolicCorrection() {}
+ParabolicCorrection::~ParabolicCorrection() {}
 
 // ---------------------------------------------------------
-// 1. 自动寻找真空层中心 (逻辑完全借鉴 Efield::autoset)
+// 1. 自动寻找真空层中心
 // ---------------------------------------------------------
 double ParabolicCorrection::find_vacuum_center(const UnitCell& cell, int dir)
 {
@@ -19,17 +22,15 @@ double ParabolicCorrection::find_vacuum_center(const UnitCell& cell, int dir)
         for (int ia = 0; ia < cell.atoms[it].na; ++ia)
         {
             double p = cell.atoms[it].taud[ia][dir];
-            p -= std::floor(p); // 归一化到 [0, 1)
+            p -= std::floor(p); 
             pos.push_back(p);
         }
     }
-
     std::sort(pos.begin(), pos.end());
 
     double max_gap = 0.0;
     double center = 0.5;
 
-    // 内部间隙
     for (size_t i = 1; i < pos.size(); i++)
     {
         double gap = pos[i] - pos[i - 1];
@@ -40,7 +41,6 @@ double ParabolicCorrection::find_vacuum_center(const UnitCell& cell, int dir)
         }
     }
 
-    // 跨边界间隙
     if (!pos.empty()) {
         double tail_gap = pos[0] + 1.0 - pos.back();
         if (tail_gap > max_gap)
@@ -50,57 +50,53 @@ double ParabolicCorrection::find_vacuum_center(const UnitCell& cell, int dir)
             if (center >= 1.0) center -= 1.0;
         }
     }
-
     return center;
 }
 
-ParabolicCorrection::ParabolicCorrection() {}
-ParabolicCorrection::~ParabolicCorrection() {}
-
 // ---------------------------------------------------------
-// 2. 应用校正的主函数
+// 2. 应用校正的主函数 (集成离子能量计算)
 // ---------------------------------------------------------
-void ParabolicCorrection::apply_correction(const UnitCell& cell, 
-                                           const ModulePW::PW_Basis* rho_basis, 
-                                           double* v_hartree, 
-                                           const double* const* rho_elec, // 【修正1】指针的指针
-                                           int nspin,                     // 【修正2】添加 nspin 参数
-                                           double nelec,
-                                           int dir)
+double ParabolicCorrection::apply_correction(const UnitCell& cell, 
+                                             const ModulePW::PW_Basis* rho_basis, 
+                                             double* v_hartree, 
+                                             const double* const* rho_elec, 
+                                             int nspin,                     
+                                             double nelec,
+                                             int dir)
 {
     // 基础几何参数
     double omega = cell.omega;       
     double lat_vec = 0.0;            
     double area = 0.0;               
 
-    if (dir == 0) {
-        lat_vec = cell.a1.norm() * cell.lat0;
-    } else if (dir == 1) {
-        lat_vec = cell.a2.norm() * cell.lat0;
-    } else { // dir == 2
-        lat_vec = cell.a3.norm() * cell.lat0;
-    }
+    if (dir == 0) lat_vec = cell.a1.norm() * cell.lat0;
+    else if (dir == 1) lat_vec = cell.a2.norm() * cell.lat0;
+    else lat_vec = cell.a3.norm() * cell.lat0;
+    
     area = omega / lat_vec;
 
-    // 1. 自动寻找真空层中心
+    // 1. 寻找中心
     double vacuum_center = find_vacuum_center(cell, dir);
-
-    // 2. 定义 Slab 的几何中心 (真空中心的对面)
     double slab_center = vacuum_center + 0.5;
     if(slab_center >= 1.0) slab_center -= 1.0;
 
-    // 3. 计算净电荷 (Q_ion - Q_elec)
+    // 2. 计算电荷与偶极
     double net_charge = calc_net_charge(cell, nelec); 
-    
-    // 4. 计算总偶极矩 (传递 slab_center 作为参考原点)
     double total_dipole = calc_total_dipole(cell, rho_basis, rho_elec, nspin, dir, slab_center);
+    
+    // 缓存偶极矩供后续力修正使用
+    this->last_total_dipole_ = total_dipole;
 
-    // 5. 构造 1D 修正势
-    // 系数 factor = 4pi / Area * e^2
-    double factor = (ModuleBase::FOUR_PI / area) * ModuleBase::e2;
+    // 3. 计算离子修正能 (这是你想要加入的！)
+    double e_ion_corr = calc_energy_correction(cell, dir, net_charge, total_dipole, vacuum_center);
 
+    // 4. 构造 1D 修正势并叠加到 v_hartree
+    double factor = (ModuleBase::FOUR_PI / area) * ModuleBase::e2;
     int nrxx = rho_basis->nrxx; 
 
+    #ifdef _OPENMP
+    #pragma omp parallel for
+    #endif
     for (int ir = 0; ir < nrxx; ++ir)
     {
         int i = ir / (rho_basis->ny * rho_basis->nplane);
@@ -125,8 +121,14 @@ void ParabolicCorrection::apply_correction(const UnitCell& cell,
 
         v_hartree[ir] += v_corr;
     }
+
+    // 返回离子修正能，方便外部加到 Total Energy
+    return e_ion_corr;
 }
 
+// ---------------------------------------------------------
+// 辅助计算函数
+// ---------------------------------------------------------
 double ParabolicCorrection::calc_net_charge(const UnitCell& cell, double nelec)
 {
     double ion_charge = 0.0;
@@ -138,8 +140,8 @@ double ParabolicCorrection::calc_net_charge(const UnitCell& cell, double nelec)
 
 double ParabolicCorrection::calc_total_dipole(const UnitCell& cell, 
                                               const ModulePW::PW_Basis* rho_basis,
-                                              const double* const* rho_elec, // 【修正】类型匹配
-                                              int nspin,                     // 【修正】传递 nspin
+                                              const double* const* rho_elec,
+                                              int nspin,
                                               int dir,
                                               double center)
 {
@@ -170,20 +172,15 @@ double ParabolicCorrection::calc_ion_dipole(const UnitCell& cell, int dir, doubl
     return d;
 }
 
-// ---------------------------------------------------------
-// 3. 计算电子偶极矩 (包含多自旋求和)
-// ---------------------------------------------------------
 double ParabolicCorrection::calc_elec_dipole(const ModulePW::PW_Basis* rho_basis, 
-                                             const double* const* rho_elec, // 【修正】指针的指针
-                                             int nspin,                     // 【修正】接收 nspin
+                                             const double* const* rho_elec, 
+                                             int nspin,
                                              int dir, 
                                              double center,
                                              double omega)
 {
     double d = 0.0;
     int nrxx = rho_basis->nrxx;
-
-    // 【修正】Runtime 下，rho[0]=Up, rho[1]=Down，需要求和
     int n_components = (nspin == 2) ? 2 : 1;
 
     for (int ir = 0; ir < nrxx; ++ir)
@@ -203,15 +200,102 @@ double ParabolicCorrection::calc_elec_dipole(const ModulePW::PW_Basis* rho_basis
 
         double rho_val = 0.0;
         for(int is=0; is<n_components; ++is) {
-            rho_val += rho_elec[is][ir]; // 【修正】正确访问二维数组
+            rho_val += rho_elec[is][ir];
         }
 
         d += rho_val * dist;
-    } // 【修正】补回了丢失的括号
+    } // 【修复了这里丢失的括号】
 
     Parallel_Reduce::reduce_pool(d);
-
     d *= (omega / rho_basis->nxyz); 
-
     return d;
+}
+
+// ---------------------------------------------------------
+// 离子能量校正实现
+// ---------------------------------------------------------
+double ParabolicCorrection::calc_energy_correction(const UnitCell& cell, 
+                                                   int dir,
+                                                   double net_charge,
+                                                   double total_dipole,
+                                                   double vacuum_center)
+{
+    double lat_vec = 0.0;
+    if (dir == 0) lat_vec = cell.a1.norm() * cell.lat0;
+    else if (dir == 1) lat_vec = cell.a2.norm() * cell.lat0;
+    else lat_vec = cell.a3.norm() * cell.lat0;
+    
+    double area = cell.omega / lat_vec;
+    double factor = (ModuleBase::FOUR_PI / area) * ModuleBase::e2;
+    
+    double slab_center = vacuum_center + 0.5;
+    if(slab_center >= 1.0) slab_center -= 1.0;
+
+    double e_corr = 0.0;
+
+    for(int it=0; it<cell.ntype; ++it) {
+        double Z = cell.atoms[it].ncpp.zv;
+        for(int ia=0; ia<cell.atoms[it].na; ++ia) {
+            double pos = cell.atoms[it].taud[ia][dir];
+            
+            double dist_frac = pos - slab_center;
+            if (dist_frac > 0.5) dist_frac -= 1.0;
+            if (dist_frac < -0.5) dist_frac += 1.0;
+            double dist_bohr = dist_frac * lat_vec;
+
+            // V_corr at atom position
+            double v_at_atom = factor * ( -0.5 * net_charge * dist_bohr * dist_bohr 
+                                          + total_dipole * dist_bohr );
+            
+            e_corr += Z * v_at_atom;
+        }
+    }
+    return e_corr;
+}
+
+// ---------------------------------------------------------
+// 力校正实现 (保持独立，供 force 模块调用)
+// ---------------------------------------------------------
+void ParabolicCorrection::calc_force_correction(const UnitCell& cell, 
+                                                ModuleBase::matrix& force, 
+                                                int dir,
+                                                double net_charge,
+                                                double total_dipole,
+                                                double vacuum_center,
+                                                double area)
+{
+    double factor = (ModuleBase::FOUR_PI / area) * ModuleBase::e2;
+    double slab_center = vacuum_center + 0.5;
+    if(slab_center >= 1.0) slab_center -= 1.0;
+    
+    double lat_vec = cell.omega / area; // 反推 L
+
+    int iat = 0;
+    for(int it=0; it<cell.ntype; ++it) {
+        double Z = cell.atoms[it].ncpp.zv;
+        for(int ia=0; ia<cell.atoms[it].na; ++ia) {
+            
+            double pos = cell.atoms[it].taud[ia][dir];
+            double dist_frac = pos - slab_center;
+            if (dist_frac > 0.5) dist_frac -= 1.0;
+            if (dist_frac < -0.5) dist_frac += 1.0;
+            double dist = dist_frac * lat_vec;
+
+            // 1. Electric Field Force: -Z * dV/dz
+            double f_field = -Z * factor * (-net_charge * dist + total_dipole);
+
+            // 2. Dipole Response Force: Z * d(Total_Dipole)/dz * dE/dD
+            double f_dipole = factor * Z * (net_charge * dist - total_dipole);
+            
+            // 总修正力
+            double f_corr = f_field + f_dipole; 
+            // 注意：上面的推导中两项可能会相互抵消或合并，具体依赖于公式的变分导数
+            // 简单验证：Environ 中 f = (charge * pos - dipole) * fact * Z
+            // 这里的 f_corr 简化后确实类似。
+            
+            force(iat, dir) += f_corr;
+            
+            iat++;
+        }
+    }
 }