Skip to content

Commit eb0cdb8

Browse files
YuLiu98YuLiu98
committed
feature : dipole correction for surface calculations
1 parent 3c80d55 commit eb0cdb8

File tree

7 files changed

+229
-2
lines changed

7 files changed

+229
-2
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -285,6 +285,7 @@ write_HS_R.o\
285285
write_dm.o\
286286
write_wfc_realspace.o\
287287
efield.o \
288+
dipole.o\
288289
magnetism.o\
289290
optical.o\
290291
Cell_PW.o\

‎source/src_pw/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@ list(APPEND objects
1010
charge_pulay.cpp
1111
diago_cg.cpp
1212
diago_david.cpp
13+
dipole.cpp
1314
efield.cpp
1415
energy.cpp
1516
forces.cpp

‎source/src_pw/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -72,6 +72,7 @@ charge_broyden.o \
7272
charge_extra.o \
7373
diago_david.o\
7474
diago_cg.o\
75+
dipole.o\
7576
electrons.o \
7677
efield.o \
7778
energy.o \

‎source/src_pw/dipole.cpp‎

Lines changed: 188 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
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+
10+
Dipole::Dipole(){}
11+
12+
Dipole::~Dipole(){}
13+
14+
//=======================================================
15+
// calculate dipole potential in surface calculations
16+
//=======================================================
17+
ModuleBase::matrix Dipole::v_dipole(const UnitCell &cell,
18+
PW_Basis &pwb,
19+
const int &nspin,
20+
const double *const *const rho)
21+
{
22+
ModuleBase::TITLE("Dipole", "v_dipole");
23+
ModuleBase::timer::tick("Dipole", "v_dipole");
24+
25+
const double m = dipole_density(cell, pwb, nspin, rho);
26+
const double fac = m * ModuleBase::FOUR_PI;
27+
ModuleBase::Vector3<int> gird;
28+
gird.set(pwb.ncx, pwb.ncy, pwb.ncz);
29+
30+
ModuleBase::matrix v(nspin, pwb.nrxx);
31+
double e_elec = 0;
32+
33+
for (int ir = 0; ir < pwb.nrxx; ++ir)
34+
{
35+
int i = ir / (pwb.ncy * pwb.nczp);
36+
int j = ir / pwb.nczp - i * pwb.ncy;
37+
int k = ir % pwb.nczp + pwb.nczp_start;
38+
ModuleBase::Vector3<int> index;
39+
index.set(i,j,k);
40+
41+
if (nspin == 4)
42+
{
43+
v(0, ir) = index[dir] / gird[dir] - 0.5;
44+
e_elec += v(0, ir) * rho[0][ir];
45+
}
46+
else
47+
{
48+
for (int is = 0; is < nspin; is++)
49+
{
50+
v(is, ir) = index[dir] / gird[dir] - 0.5;
51+
e_elec += v(is, ir) * rho[is][ir];
52+
}
53+
}
54+
}
55+
56+
Parallel_Reduce::reduce_double_pool(e_elec);
57+
e_elec *= cell.omega / pwb.ncxyz;
58+
59+
// which_ip: found iz belongs to which ip.
60+
int *which_ip = new int[pwb.ncz];
61+
ModuleBase::GlobalFunc::ZEROS(which_ip, pwb.ncz);
62+
63+
#ifdef __MPI
64+
// num_z: how many planes on processor 'ip'
65+
int *num_z = new int[GlobalV::NPROC_IN_POOL];
66+
ModuleBase::GlobalFunc::ZEROS(num_z, GlobalV::NPROC_IN_POOL);
67+
for (int iz=0;iz<pwb.nbz;iz++)
68+
{
69+
int ip = iz % GlobalV::NPROC_IN_POOL;
70+
num_z[ip] += pwb.bz;
71+
}
72+
73+
// start_z: start position of z in processor ip.
74+
int *start_z = new int[GlobalV::NPROC_IN_POOL];
75+
ModuleBase::GlobalFunc::ZEROS(start_z, GlobalV::NPROC_IN_POOL);
76+
for (int ip=1;ip<GlobalV::NPROC_IN_POOL;ip++)
77+
{
78+
start_z[ip] = start_z[ip-1]+num_z[ip-1];
79+
}
80+
81+
for(int iz=0; iz<pwb.ncz; iz++)
82+
{
83+
for(int ip=0; ip<GlobalV::NPROC_IN_POOL; ip++)
84+
{
85+
if(iz>=start_z[GlobalV::NPROC_IN_POOL-1])
86+
{
87+
which_ip[iz] = GlobalV::NPROC_IN_POOL-1;
88+
break;
89+
}
90+
else if(iz>=start_z[ip] && iz<start_z[ip+1])
91+
{
92+
which_ip[iz] = ip;
93+
break;
94+
}
95+
}
96+
}
97+
98+
delete[] num_z;
99+
delete[] start_z;
100+
#endif
101+
102+
double e_ion = 0;
103+
for(int it=0; it<cell.ntype; ++it)
104+
{
105+
for(int ia=0; ia<cell.atoms[it].na; ++ia)
106+
{
107+
int iz = cell.atoms[it].taud[ia][2] * pwb.ncz; // FFT divided by z axis
108+
if(GlobalV::RANK_IN_POOL == which_ip[iz])
109+
{
110+
int ix = cell.atoms[it].taud[ia][0] * pwb.ncx;
111+
int iy = cell.atoms[it].taud[ia][1] * pwb.ncy;
112+
iz -= pwb.nczp_start;
113+
int index = ix * pwb.ncy * pwb.nczp + iy * pwb.nczp + iz;
114+
115+
e_ion += cell.atoms[it].zv * v(0, index);
116+
}
117+
}
118+
}
119+
120+
Parallel_Reduce::reduce_double_pool(e_ion);
121+
122+
dipole_energy = 0.5 * fac * (e_ion - e_elec);
123+
v *= fac;
124+
125+
delete[] which_ip;
126+
127+
ModuleBase::timer::tick("Dipole", "v_dipole");
128+
return v;
129+
}
130+
131+
132+
//=======================================================
133+
// calculate dipole density in surface calculations
134+
//=======================================================
135+
double Dipole::dipole_density(const UnitCell &cell,
136+
PW_Basis &pwb,
137+
const int &nspin,
138+
const double *const *const rho)
139+
{
140+
// ion part
141+
double m_ion = 0;
142+
for(int it=0; it<cell.ntype; ++it)
143+
{
144+
double sum = 0;
145+
for(int ia=0; ia<cell.atoms[it].na; ++ia)
146+
{
147+
sum += cell.atoms[it].tau[ia][dir];
148+
}
149+
m_ion += sum * cell.atoms[it].zv;
150+
}
151+
m_ion *= cell.lat0 * ModuleBase::e2; // multiply ModuleBase::e2 to convert Ha to Ry ??
152+
153+
// height and surface area
154+
const double h = cell.latvec.to_matrix()(dir, dir) * cell.lat0;
155+
const double area = cell.omega / h;
156+
157+
m_ion /= area;
158+
159+
160+
// electron part
161+
double m_elec = 0;
162+
const int nspin0 = (nspin == 2) ? 2 : 1;
163+
164+
165+
for (int ir = 0; ir < pwb.nrxx; ++ir)
166+
{
167+
int i = ir / (pwb.ncy * pwb.nczp);
168+
int j = ir / pwb.nczp - i * pwb.ncy;
169+
int k = ir % pwb.nczp + pwb.nczp_start;
170+
ModuleBase::Vector3<int> index;
171+
index.set(i,j,k);
172+
173+
for (int is = 0; is < nspin0; is++)
174+
{
175+
m_elec += rho[is][ir] * index[dir];
176+
}
177+
}
178+
179+
Parallel_Reduce::reduce_double_pool(m_elec);
180+
181+
ModuleBase::Vector3<int> gird;
182+
gird.set(pwb.ncx, pwb.ncy, pwb.ncz);
183+
184+
m_elec *= h * h / pwb.ncxyz / gird[dir] * ModuleBase::e2;
185+
186+
return m_ion - m_elec;
187+
}
188+

‎source/src_pw/dipole.h‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
#ifndef DIPOLE_H
2+
#define DIPOLE_H
3+
4+
#include "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_density(const UnitCell &cell, PW_Basis &pwb, const int &nspin, const double *const *const rho);
14+
static ModuleBase::matrix v_dipole(const UnitCell &cell, PW_Basis &pwb, const int &nspin, const double *const *const rho);
15+
16+
17+
18+
static double dipole_energy; // dipole energy
19+
static int dir; // 0, 1, 2 denotes x, y, z direction for dipole correction
20+
};
21+
22+
23+
24+
25+
26+
#endif

‎source/src_pw/energy.cpp‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@
1616
//new
1717
#include "H_Ewald_pw.h"
1818
#include "H_Hartree_pw.h"
19+
#include "dipole.h" // liuyu add 2022-05-06
1920
#ifdef __DEEPKS
2021
#include "../module_deepks/LCAO_deepks.h"
2122
#endif
@@ -60,7 +61,9 @@ void energy::calculate_harris(const int &flag)
6061
+ demet
6162
+ exx
6263
+ Efield::etotefield
63-
+ evdw; // Peize Lin add evdw 2021.03.09
64+
+ evdw // Peize Lin add evdw 2021.03.09
65+
+ Dipole::dipole_energy; // liuyu add 2022-05-06
66+
6467
#ifdef __LCAO
6568
if(INPUT.dft_plus_u)
6669
{
@@ -90,7 +93,8 @@ void energy::calculate_etot(void)
9093
+ descf
9194
+ exx
9295
+ Efield::etotefield
93-
+ evdw; // Peize Lin add evdw 2021.03.09
96+
+ evdw // Peize Lin add evdw 2021.03.09
97+
+ Dipole::dipole_energy; // liuyu add 2022-05-06
9498

9599
//Quxin adds for DFT+U energy correction on 20201029
96100
/*

‎source/src_pw/potential.cpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,7 @@
1010
// new
1111
#include "../module_surchem/surchem.h"
1212
#include "H_Hartree_pw.h"
13+
#include "dipole.h"
1314
#ifdef __LCAO
1415
#include "../src_lcao/ELEC_evolve.h"
1516
#endif
@@ -362,6 +363,11 @@ ModuleBase::matrix Potential::v_of_rho(const double *const *const rho_in, const
362363
}
363364
}
364365

366+
if (!GlobalV::EFIELD && GlobalV::DIPOLE)
367+
{
368+
v += Dipole::v_dipole(GlobalC::ucell, GlobalC::pw, GlobalV::NSPIN, rho_in);
369+
}
370+
365371
// mohan add 2011-06-20
366372
if (GlobalV::EFIELD && GlobalV::DIPOLE)
367373
{

0 commit comments

Comments
 (0)