Skip to content

Commit 22c0d77

Browse files
YuLiu98YuLiu98
committed
feature : dipole correction for pw
1 parent d2d579f commit 22c0d77

File tree

10 files changed

+260
-219
lines changed

10 files changed

+260
-219
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -244,6 +244,7 @@ OBJS_SURCHEM=H_correction_pw.o\
244244
cal_vel.o\
245245
corrected_energy.o\
246246
minimize_cg.o\
247+
dipole.o\
247248

248249
OBJS_XC=xc_funct_corr_gga.o \
249250
xc_funct_corr_lda.o \
@@ -295,7 +296,6 @@ write_HS_R.o\
295296
write_dm.o\
296297
write_wfc_realspace.o\
297298
efield.o \
298-
dipole.o\
299299
magnetism.o\
300300
optical.o\
301301
Cell_PW.o\

‎source/module_surchem/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,4 +9,5 @@ add_library(
99
cal_vel.cpp
1010
corrected_energy.cpp
1111
minimize_cg.cpp
12+
dipole.cpp
1213
)

‎source/module_surchem/dipole.cpp‎

Lines changed: 198 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,198 @@
1+
#include "dipole.h"
2+
#include "../module_base/constants.h"
3+
#include "../module_base/timer.h"
4+
#include "../module_base/global_variable.h"
5+
#include "../src_parallel/parallel_reduce.h"
6+
7+
int Dipole::dir = 2;
8+
double Dipole::dipole_energy = 0.0;
9+
double Dipole::max_pos = 0.5;
10+
double Dipole::de_reg = 0.1;
11+
12+
Dipole::Dipole(){}
13+
14+
Dipole::~Dipole(){}
15+
16+
//=======================================================
17+
// calculate dipole potential in surface calculations
18+
//=======================================================
19+
ModuleBase::matrix Dipole::v_dipole(const UnitCell &cell,
20+
PW_Basis &pwb,
21+
const int &nspin,
22+
const double *const *const rho)
23+
{
24+
ModuleBase::TITLE("Dipole", "v_dipole");
25+
ModuleBase::timer::tick("Dipole", "v_dipole");
26+
27+
double h_inv; // inverse of height
28+
if(dir == 0)
29+
{
30+
h_inv = sqrt(cell.G.e11 * cell.G.e11 + cell.G.e12 * cell.G.e12 + cell.G.e13 * cell.G.e13);
31+
}
32+
else if(dir == 1)
33+
{
34+
h_inv = sqrt(cell.G.e21 * cell.G.e21 + cell.G.e22 * cell.G.e22 + cell.G.e23 * cell.G.e23);
35+
}
36+
else if(dir = 2)
37+
{
38+
h_inv = sqrt(cell.G.e31 * cell.G.e31 + cell.G.e32 * cell.G.e32 + cell.G.e33 * cell.G.e33);
39+
}
40+
else
41+
{
42+
ModuleBase::WARNING_QUIT("Dipole::ion_dipole", "dipole direction is wrong!");
43+
}
44+
45+
double ion_dipole = cal_ion_dipole(cell, h_inv);
46+
double elec_dipole = cal_elec_dipole(cell, pwb, nspin, rho, h_inv);
47+
double tot_dipole = ion_dipole - elec_dipole;
48+
49+
// dipole energy correction
50+
dipole_energy = 0.5 * ModuleBase::e2 * tot_dipole * tot_dipole * cell.omega / ModuleBase::FOUR_PI;
51+
52+
// dipole potential
53+
ModuleBase::matrix v(nspin, pwb.nrxx);
54+
const int nspin0 = (nspin == 2) ? 2 : 1;
55+
56+
for (int ir = 0; ir < pwb.nrxx; ++ir)
57+
{
58+
int i = ir / (pwb.ncy * pwb.nczp);
59+
int j = ir / pwb.nczp - i * pwb.ncy;
60+
int k = ir % pwb.nczp + pwb.nczp_start;
61+
double x = (double)i / pwb.ncx;
62+
double y = (double)j / pwb.ncy;
63+
double z = (double)k / pwb.ncz;
64+
ModuleBase::Vector3<double> pos(x, y, z);
65+
66+
double saw = saw_function(max_pos, de_reg, pos[dir]);
67+
68+
for (int is = 0; is < nspin0; is++)
69+
{
70+
v(is, ir) = saw;
71+
}
72+
}
73+
74+
double fac = - ModuleBase::e2 * tot_dipole * cell.lat0 / h_inv;
75+
76+
ModuleBase::timer::tick("Dipole", "v_dipole");
77+
return v * fac;
78+
}
79+
80+
81+
//=======================================================
82+
// calculate dipole density in surface calculations
83+
//=======================================================
84+
// double Dipole::dipole(const UnitCell &cell,
85+
// PW_Basis &pwb,
86+
// const double *const *const rho,
87+
// complex<double> *TOTN)
88+
// {
89+
// double *Porter = new double[pwb.nrxx];
90+
91+
// for (int ir = 0; ir < pwb.nrxx; ++ir)
92+
// {
93+
// int i = ir / (pwb.ncy * pwb.nczp);
94+
// int j = ir / pwb.nczp - i * pwb.ncy;
95+
// int k = ir % pwb.nczp + pwb.nczp_start;
96+
// double x = (double)i / pwb.ncx;
97+
// double y = (double)j / pwb.ncy;
98+
// double z = (double)k / pwb.ncz;
99+
// ModuleBase::Vector3<double> pos(x, y, z);
100+
101+
// pos = cell.latvec * pos;
102+
103+
// Porter[ir] = cell.lat0 * pos[dir];
104+
// }
105+
106+
// complex<double> *Porter_g = new complex<double>[pwb.ngmc];
107+
// GlobalC::UFFT.ToReciSpace(Porter, Porter_g);
108+
109+
// double m = 0;
110+
// for (int ig = 0; ig < pwb.ngmc; ig++)
111+
// {
112+
// m += (conj(TOTN[ig]) * Porter_g[ig]).real();
113+
// }
114+
115+
// Parallel_Reduce::reduce_double_pool(m);
116+
117+
// // height
118+
// const double h = cell.latvec.to_matrix()(dir, dir) * cell.lat0;
119+
120+
// delete[] Porter, Porter_g;
121+
122+
// return -m * h * ModuleBase::e2;
123+
// }
124+
125+
double Dipole::cal_ion_dipole(const UnitCell &cell, const double &h_inv)
126+
{
127+
double ion_dipole = 0;
128+
for(int it=0; it<cell.ntype; ++it)
129+
{
130+
double sum = 0;
131+
for(int ia=0; ia<cell.atoms[it].na; ++ia)
132+
{
133+
sum += saw_function(max_pos, de_reg, cell.atoms[it].taud[ia][dir]);
134+
}
135+
ion_dipole += sum * cell.atoms[it].zv;
136+
}
137+
ion_dipole *= cell.lat0 / h_inv * ModuleBase::FOUR_PI / cell.omega;
138+
139+
// std::cout << "ion_dipole = " << ion_dipole << std::endl;
140+
141+
return ion_dipole;
142+
}
143+
144+
double Dipole::cal_elec_dipole(const UnitCell &cell,
145+
PW_Basis &pwb,
146+
const int &nspin,
147+
const double *const *const rho,
148+
const double &h_inv)
149+
{
150+
double elec_dipole = 0;
151+
const int nspin0 = (nspin == 2) ? 2 : 1;
152+
153+
for (int ir = 0; ir < pwb.nrxx; ++ir)
154+
{
155+
int i = ir / (pwb.ncy * pwb.nczp);
156+
int j = ir / pwb.nczp - i * pwb.ncy;
157+
int k = ir % pwb.nczp + pwb.nczp_start;
158+
double x = (double)i / pwb.ncx;
159+
double y = (double)j / pwb.ncy;
160+
double z = (double)k / pwb.ncz;
161+
ModuleBase::Vector3<double> pos(x, y, z);
162+
163+
double saw = saw_function(max_pos, de_reg, pos[dir]);
164+
165+
for (int is = 0; is < nspin0; is++)
166+
{
167+
elec_dipole += rho[is][ir] * saw;
168+
}
169+
}
170+
171+
Parallel_Reduce::reduce_double_pool(elec_dipole);
172+
elec_dipole *= cell.lat0 / h_inv * ModuleBase::FOUR_PI / pwb.ncxyz;
173+
174+
// std::cout << "elec_dipole = " << elec_dipole << std::endl;
175+
176+
return elec_dipole;
177+
}
178+
179+
double Dipole::saw_function(const double &a, const double &b, const double &x)
180+
{
181+
assert(x>=0);
182+
assert(x<=1);
183+
184+
const double fac = 1 - b;
185+
186+
if( x < a )
187+
{
188+
return x - a + 0.5 * fac;
189+
}
190+
else if( x > (a+b))
191+
{
192+
return x - a - 1 + 0.5 * fac;
193+
}
194+
else
195+
{
196+
return 0.5 * fac - fac * (x - a) / b;
197+
}
198+
}

‎source/module_surchem/dipole.h‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
#ifndef DIPOLE_H
2+
#define DIPOLE_H
3+
4+
#include "../src_pw/pw_basis.h"
5+
#include "../module_cell/unitcell.h"
6+
7+
class Dipole
8+
{
9+
public:
10+
Dipole();
11+
~Dipole();
12+
13+
static double dipole(const UnitCell &cell,
14+
PW_Basis &pwb,
15+
const double *const *const rho,
16+
complex<double> *TOTN);
17+
18+
static ModuleBase::matrix v_dipole(const UnitCell &cell,
19+
PW_Basis &pwb,
20+
const int &nspin,
21+
const double *const *const rho);
22+
23+
static double cal_elec_dipole(const UnitCell &cell,
24+
PW_Basis &pwb,
25+
const int &nspin,
26+
const double *const *const rho,
27+
const double &h_inv);
28+
29+
static double cal_ion_dipole(const UnitCell &cell, const double &h_inv);
30+
31+
static double saw_function(const double &a, const double &b, const double &x);
32+
33+
34+
35+
static double dipole_energy; // dipole energy
36+
static int dir; // 0, 1, 2 denotes x, y, z direction for dipole correction
37+
static double max_pos; // the maximum position of the saw function
38+
static double de_reg; // the decrease region length of the saw function
39+
};
40+
41+
42+
43+
44+
45+
#endif

‎source/src_io/write_pot.cpp‎

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#include "../src_pw/potential.h"
22
#include "../src_pw/global.h"
3+
#include "../module_surchem/dipole.h"
34
#include "../module_base/timer.h"
45

56
// translate from write_rho in charge.cpp.
@@ -251,12 +252,22 @@ void Potential::write_elecstat_pot(const std::string &fn, const std::string &fn_
251252
//transform hartree potential to real space
252253
//==========================================
253254
GlobalC::pw.FFT_chg.FFT3D(Porter, 1);
255+
256+
//==========================================
257+
// Dipole correction
258+
//==========================================
259+
ModuleBase::matrix v_dip(GlobalV::NSPIN, GlobalC::pw.nrxx);
260+
if (!GlobalV::EFIELD && GlobalV::DIPOLE)
261+
{
262+
v_dip = Dipole::v_dipole(GlobalC::ucell, GlobalC::pw, GlobalV::NSPIN, GlobalC::CHR.rho);
263+
}
264+
254265
//==========================================
255266
//Add hartree potential and local pseudopot
256267
//==========================================
257268
for (int ir = 0;ir < GlobalC::pw.nrxx;ir++)
258269
{
259-
v_elecstat[ir] = Porter[ir].real() + this->vltot[ir];
270+
v_elecstat[ir] = Porter[ir].real() + this->vltot[ir] + v_dip(0, ir);
260271
}
261272

262273
//-------------------------------------------

‎source/src_pw/CMakeLists.txt‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,6 @@ list(APPEND objects
1010
charge_pulay.cpp
1111
diago_cg.cpp
1212
diago_david.cpp
13-
dipole.cpp
1413
efield.cpp
1514
energy.cpp
1615
forces.cpp

0 commit comments

Comments
 (0)