Skip to content

Commit b414508

Browse files
ouqi0711ouqi0711
authored
Merge pull request #961 from Asuna981002/develop
feature: add a function of compensating charge to the surchem module.
2 parents 4b9442b + ed1beb8 commit b414508

File tree

14 files changed

+275
-50
lines changed

14 files changed

+275
-50
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -240,6 +240,7 @@ OBJS_ESOLVER=esolver.o\
240240

241241
OBJS_SURCHEM=H_correction_pw.o\
242242
cal_epsilon.o\
243+
surchem.o\
243244
cal_pseudo.o\
244245
cal_totn.o\
245246
cal_vcav.o\
@@ -415,7 +416,6 @@ pzt2s.o\
415416
pdtrsm.o\
416417
pzhtrsm.o\
417418

418-
419419
PDIAG_MR_0=dcopy.o\
420420
dlae2.o\
421421
dlaebz.o\

‎source/input.cpp‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -420,6 +420,11 @@ void Input::Default(void)
420420
sigma_k = 0.6;
421421
nc_k = 0.00037;
422422

423+
comp_q = 0.0;
424+
comp_l = 1.0;
425+
comp_center = 0.0;
426+
comp_dim = 2;
427+
423428
return;
424429
}
425430

@@ -1495,6 +1500,22 @@ bool Input::Read(const std::string &fn)
14951500
{
14961501
read_value(ifs, nc_k);
14971502
}
1503+
else if (strcmp("comp_q", word) == 0)
1504+
{
1505+
read_value(ifs, comp_q);
1506+
}
1507+
else if (strcmp("comp_l", word) == 0)
1508+
{
1509+
read_value(ifs, comp_l);
1510+
}
1511+
else if (strcmp("comp_center", word) == 0)
1512+
{
1513+
read_value(ifs, comp_center);
1514+
}
1515+
else if (strcmp("comp_dim", word) == 0)
1516+
{
1517+
read_value(ifs, comp_dim);
1518+
}
14981519
//----------------------------------------------------------------------------------
14991520
else
15001521
{

‎source/input.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -399,6 +399,11 @@ class Input
399399
double tau;
400400
double sigma_k;
401401
double nc_k;
402+
// gauss charge
403+
double comp_q;
404+
double comp_l;
405+
double comp_center;
406+
int comp_dim;
402407

403408
//==========================================================
404409
// variables for test only

‎source/input_conv.cpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@
1414
#include "src_ions/ions_move_basic.h"
1515
#include "src_pw/global.h"
1616
#include "src_pw/occupy.h"
17+
#include "module_surchem/surchem.h"
1718
#ifdef __EXX
1819
#include "src_ri/exx_abfs-jle.h"
1920
#endif
@@ -501,6 +502,11 @@ void Input_Conv::Convert(void)
501502
GlobalV::tau = INPUT.tau;
502503
GlobalV::sigma_k = INPUT.sigma_k;
503504
GlobalV::nc_k = INPUT.nc_k;
505+
506+
GlobalC::solvent_model.comp_q = INPUT.comp_q;
507+
GlobalC::solvent_model.comp_l = INPUT.comp_l;
508+
GlobalC::solvent_model.comp_center = INPUT.comp_center;
509+
GlobalC::solvent_model.comp_dim = INPUT.comp_dim;
504510
ModuleBase::timer::tick("Input_Conv", "Convert");
505511
return;
506512
}

‎source/module_surchem/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
add_library(
22
surchem
33
OBJECT
4+
surchem.cpp
45
H_correction_pw.cpp
56
cal_epsilon.cpp
67
cal_pseudo.cpp

‎source/module_surchem/cal_pseudo.cpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
#include "surchem.h"
22

3-
atom_in surchem::GetAtom;
3+
// atom_in surchem::GetAtom;
44

55
void surchem::gauss_charge(const UnitCell &cell, PW_Basis &pwb, complex<double> *N)
66
{

‎source/module_surchem/cal_totn.cpp‎

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,94 @@
11
#include "surchem.h"
22

3+
void add_comp_chg(const UnitCell &cell, PW_Basis &pwb, double q, double l, double center, complex<double> *NG, int dim)
4+
{
5+
// x dim
6+
if (dim == 0)
7+
{
8+
double L = cell.a1[0];
9+
q /= (cell.a2[1] * cell.a3[2] * l);
10+
ModuleBase::GlobalFunc::ZEROS(NG, pwb.ngmc);
11+
for (int ig = pwb.gstart; ig < pwb.ngmc; ig++)
12+
{
13+
double GX = pwb.get_G_cartesian_projection(ig, 0);
14+
double GY = pwb.get_G_cartesian_projection(ig, 1);
15+
double GZ = pwb.get_G_cartesian_projection(ig, 2);
16+
GX = GX * 2 * ModuleBase::PI;
17+
if (GY == 0 && GZ == 0 && GX != 0)
18+
{
19+
NG[ig] = exp(ModuleBase::NEG_IMAG_UNIT * GX * center) * complex<double>(2.0 * q * sin(GX * l / 2.0) / (L * GX), 0.0);
20+
}
21+
}
22+
NG[0] = complex<double>(q * l / L, 0.0);
23+
}
24+
// y dim
25+
else if (dim == 1)
26+
{
27+
double L = cell.a2[1];
28+
q /= (cell.a1[0] * cell.a3[2] * l);
29+
ModuleBase::GlobalFunc::ZEROS(NG, pwb.ngmc);
30+
for (int ig = pwb.gstart; ig < pwb.ngmc; ig++)
31+
{
32+
double GX = pwb.get_G_cartesian_projection(ig, 0);
33+
double GY = pwb.get_G_cartesian_projection(ig, 1);
34+
double GZ = pwb.get_G_cartesian_projection(ig, 2);
35+
GY = GY * 2 * ModuleBase::PI;
36+
if (GX == 0 && GZ == 0 && GY != 0)
37+
{
38+
NG[ig] = exp(ModuleBase::NEG_IMAG_UNIT * GY * center) * complex<double>(2.0 * q * sin(GY * l / 2.0) / (L * GY), 0.0);
39+
}
40+
}
41+
NG[0] = complex<double>(q * l / L, 0.0);
42+
}
43+
// z dim
44+
else if (dim == 2)
45+
{
46+
double L = cell.a3[2];
47+
q /= (cell.a1[0] * cell.a2[1] * l);
48+
ModuleBase::GlobalFunc::ZEROS(NG, pwb.ngmc);
49+
for (int ig = pwb.gstart; ig < pwb.ngmc; ig++)
50+
{
51+
double GX = pwb.get_G_cartesian_projection(ig, 0);
52+
double GY = pwb.get_G_cartesian_projection(ig, 1);
53+
double GZ = pwb.get_G_cartesian_projection(ig, 2);
54+
GZ = GZ * 2 * ModuleBase::PI;
55+
if (GX == 0 && GY == 0 && GZ != 0)
56+
{
57+
NG[ig] = exp(ModuleBase::NEG_IMAG_UNIT * GZ * center) * complex<double>(2.0 * q * sin(GZ * l / 2.0) / (L * GZ), 0.0);
58+
}
59+
}
60+
NG[0] = complex<double>(q * l / L, 0.0);
61+
}
62+
}
63+
64+
void test_sum(const UnitCell &cell, PW_Basis &pwb, double* NR)
65+
{
66+
double sum = 0.0;
67+
for (int i=0; i < pwb.nrxx; i++)
68+
{
69+
sum += NR[i];
70+
}
71+
cout << sum * cell.omega / pwb.nrxx << endl;
72+
}
73+
374
void surchem::cal_totn(const UnitCell &cell, PW_Basis &pwb,
475
const complex<double> *Porter_g, complex<double> *N,
576
complex<double> *TOTN) {
77+
/*
78+
// test compensating charge
79+
complex<double> *comp_reci = new complex<double>[pwb.ngmc];
80+
double *comp_real = new double[pwb.nrxx];
81+
add_comp_chg(cell, pwb, comp_q, comp_l, comp_center, comp_reci, comp_dim);
82+
GlobalC::UFFT.ToRealSpace(comp_reci, comp_real);
83+
for (int ir = 0; ir < pwb.nz; ir++)
84+
{
85+
cout << comp_real[ir] << endl;
86+
}
87+
cout << "sum" << endl;
88+
test_sum(cell, pwb, comp_real);
89+
delete[] comp_real;
90+
delete[] comp_reci;
91+
*/
692

793
// vloc to N
894
complex<double> *vloc_g = new complex<double>[pwb.ngmc];
@@ -23,11 +109,13 @@ void surchem::cal_totn(const UnitCell &cell, PW_Basis &pwb,
23109

24110
if (GlobalV::MY_RANK == 0) {
25111
N[0] = Porter_g[0];
112+
// cout << "Ng[0]" << N[0] << endl;
26113
}
27114

28115
for (int ig = 0; ig < pwb.ngmc; ig++) {
29116
TOTN[ig] = N[ig] - Porter_g[ig];
30117
}
118+
31119
delete[] vloc_g;
32120
return;
33121
}

‎source/module_surchem/cal_vcav.cpp‎

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -73,7 +73,7 @@ void shape_gradn(const complex<double> *PS_TOTN, PW_Basis &pw, double *eprime)
7373
delete[] PS_TOTN_real;
7474
}
7575

76-
void createcavity(const UnitCell &ucell, PW_Basis &pwb, const complex<double> *PS_TOTN, double *vwork)
76+
void surchem::createcavity(const UnitCell &ucell, PW_Basis &pwb, const complex<double> *PS_TOTN, double *vwork)
7777
{
7878
ModuleBase::Vector3<double> *nablan = new ModuleBase::Vector3<double>[pwb.nrxx];
7979
ModuleBase::GlobalFunc::ZEROS(nablan, pwb.nrxx);
@@ -122,7 +122,7 @@ void createcavity(const UnitCell &ucell, PW_Basis &pwb, const complex<double> *P
122122
// quantum surface area, integral of (gamma*A / n) * |\nabla n|
123123
//=term1 * sqrt_nablan_2
124124
//-------------------------------------------------------------
125-
double qs = 0;
125+
qs = 0;
126126

127127
for (int ir = 0; ir < pwb.nrxx; ir++)
128128
{
@@ -136,7 +136,7 @@ void createcavity(const UnitCell &ucell, PW_Basis &pwb, const complex<double> *P
136136
// cavitation energy
137137
//-------------------------------------------------------------
138138

139-
double Ael = surchem::cal_Acav(ucell, pwb, qs);
139+
// double Ael = cal_Acav(ucell, pwb);
140140

141141
// packs the real array into a complex one
142142
// to G space
@@ -176,17 +176,17 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell &ucell, PW_Basis &pwb, const
176176
ModuleBase::TITLE("surchem", "cal_vcav");
177177
ModuleBase::timer::tick("surchem", "cal_vcav");
178178

179-
double *Vcav = new double[pwb.nrxx];
180-
ModuleBase::GlobalFunc::ZEROS(Vcav, pwb.nrxx);
179+
double *tmp_Vcav = new double[pwb.nrxx];
180+
ModuleBase::GlobalFunc::ZEROS(tmp_Vcav, pwb.nrxx);
181181

182-
createcavity(ucell, pwb, PS_TOTN, Vcav);
182+
createcavity(ucell, pwb, PS_TOTN, tmp_Vcav);
183183

184-
ModuleBase::matrix v(nspin, pwb.nrxx);
184+
ModuleBase::GlobalFunc::ZEROS(Vcav.c, nspin * pwb.nrxx);
185185
if (nspin == 4)
186186
{
187187
for (int ir = 0; ir < pwb.nrxx; ir++)
188188
{
189-
v(0, ir) += Vcav[ir];
189+
Vcav(0, ir) += tmp_Vcav[ir];
190190
}
191191
}
192192
else
@@ -195,11 +195,12 @@ ModuleBase::matrix surchem::cal_vcav(const UnitCell &ucell, PW_Basis &pwb, const
195195
{
196196
for (int ir = 0; ir < pwb.nrxx; ir++)
197197
{
198-
v(is, ir) += Vcav[ir];
198+
Vcav(is, ir) += tmp_Vcav[ir];
199199
}
200200
}
201201
}
202-
delete[] Vcav;
202+
203+
delete[] tmp_Vcav;
203204
ModuleBase::timer::tick("surchem", "cal_vcav");
204-
return v;
205+
return Vcav;
205206
}

‎source/module_surchem/cal_vel.cpp‎

Lines changed: 20 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
5757
ModuleBase::TITLE("surchem", "cal_vel");
5858
ModuleBase::timer::tick("surchem", "cal_vel");
5959

60-
double *TOTN_real = new double[pwb.nrxx];
60+
// double *TOTN_real = new double[pwb.nrxx];
6161
GlobalC::UFFT.ToRealSpace(TOTN, TOTN_real);
6262

6363
// -4pi * TOTN(G)
@@ -80,8 +80,8 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
8080
complex<double> *Sol_phi0 = new complex<double>[pwb.ngmc];
8181
int ncgsol = 0;
8282

83-
double *Vel = new double[pwb.nrxx];
84-
ModuleBase::GlobalFunc::ZEROS(Vel, pwb.nrxx);
83+
double *tmp_Vel = new double[pwb.nrxx];
84+
ModuleBase::GlobalFunc::ZEROS(tmp_Vel, pwb.nrxx);
8585

8686
// Calculate Sol_phi with epsilon.
8787
ncgsol = 0;
@@ -93,39 +93,37 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
9393

9494
double *phi_tilda_R = new double[pwb.nrxx];
9595
double *phi_tilda_R0 = new double[pwb.nrxx];
96-
double *delta_phi_R = new double[pwb.nrxx];
96+
// double *delta_phi_R = new double[pwb.nrxx];
9797

9898
GlobalC::UFFT.ToRealSpace(Sol_phi, phi_tilda_R);
9999
GlobalC::UFFT.ToRealSpace(Sol_phi0, phi_tilda_R0);
100100

101-
// the 1st item of Vel
101+
// the 1st item of tmp_Vel
102102
for (int i = 0; i < pwb.nrxx; i++)
103103
{
104-
delta_phi_R[i] = phi_tilda_R[i] - phi_tilda_R0[i];
105-
Vel[i] += delta_phi_R[i];
104+
delta_phi[i] = phi_tilda_R[i] - phi_tilda_R0[i];
105+
tmp_Vel[i] += delta_phi[i];
106106
}
107107

108108
// calculate Ael
109-
double Ael = cal_Ael(cell, pwb, TOTN_real, delta_phi_R);
109+
// double Ael = cal_Ael(cell, pwb);
110110

111-
// the 2nd item of Vel
112-
double *Vel2 = new double[pwb.nrxx];
113-
ModuleBase::GlobalFunc::ZEROS(Vel2, pwb.nrxx);
114-
115-
eps_pot(PS_TOTN_real, Sol_phi, pwb, epsilon, Vel2);
111+
// the 2nd item of tmp_Vel
112+
eps_pot(PS_TOTN_real, Sol_phi, pwb, epsilon, epspot);
116113

117114
for (int i = 0; i < pwb.nrxx; i++)
118115
{
119-
Vel[i] += Vel2[i];
116+
tmp_Vel[i] += epspot[i];
120117
}
121118

122-
ModuleBase::matrix v(nspin, pwb.nrxx);
119+
// ModuleBase::matrix v(nspin, pwb.nrxx);
120+
ModuleBase::GlobalFunc::ZEROS(Vel.c, nspin * pwb.nrxx);
123121

124122
if (nspin == 4)
125123
{
126124
for (int ir = 0; ir < pwb.nrxx; ir++)
127125
{
128-
v(0, ir) += Vel[ir];
126+
Vel(0, ir) += tmp_Vel[ir];
129127
}
130128
}
131129
else
@@ -134,7 +132,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
134132
{
135133
for (int ir = 0; ir < pwb.nrxx; ir++)
136134
{
137-
v(is, ir) += Vel[ir];
135+
Vel(is, ir) += tmp_Vel[ir];
138136
}
139137
}
140138
}
@@ -145,13 +143,13 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
145143
delete[] B;
146144
delete[] epsilon;
147145
delete[] epsilon0;
148-
delete[] Vel;
149-
delete[] Vel2;
150-
delete[] TOTN_real;
146+
delete[] tmp_Vel;
147+
// delete[] Vel2;
148+
// delete[] TOTN_real;
151149
delete[] phi_tilda_R;
152150
delete[] phi_tilda_R0;
153-
delete[] delta_phi_R;
151+
// delete[] delta_phi_R;
154152

155153
ModuleBase::timer::tick("surchem", "cal_vel");
156-
return v;
154+
return Vel;
157155
}

0 commit comments

Comments
 (0)