Skip to content

Commit babd750

Browse files
dyzhengdyzheng
authored
Merge pull request #1174 from sunml99/develop
Feature: force correction for implicit solvation model
2 parents 754a300 + 44d74fb commit babd750

File tree

15 files changed

+248
-29
lines changed

15 files changed

+248
-29
lines changed

‎source/Makefile.Objects‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -245,6 +245,7 @@ OBJS_SURCHEM=H_correction_pw.o\
245245
corrected_energy.o\
246246
minimize_cg.o\
247247
efield.o\
248+
sol_force.o\
248249

249250
OBJS_XC=xc_funct_corr_gga.o \
250251
xc_funct_corr_lda.o \

‎source/module_surchem/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,4 +11,5 @@ add_library(
1111
corrected_energy.cpp
1212
minimize_cg.cpp
1313
efield.cpp
14+
sol_force.cpp
1415
)

‎source/module_surchem/H_correction_pw.cpp‎

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ ModuleBase::matrix surchem::v_correction(const UnitCell &cell,
2626
complex<double> *Porter_g = new complex<double>[rho_basis->npw];
2727
ModuleBase::GlobalFunc::ZEROS(Porter_g, rho_basis->npw);
2828

29-
GlobalC::UFFT.ToReciSpace(Porter, Porter_g, rho_basis);
29+
rho_basis->real2recip(Porter, Porter_g);
3030

3131
complex<double> *N = new complex<double>[rho_basis->npw];
3232
complex<double> *TOTN = new complex<double>[rho_basis->npw];
@@ -149,7 +149,7 @@ ModuleBase::matrix surchem::v_compensating(const UnitCell &cell, ModulePW::PW_Ba
149149
// std::cout << " ecomp=" << ecomp << std::endl;
150150
comp_chg_energy = ecomp;
151151

152-
GlobalC::UFFT.ToRealSpace(phi_comp_G, phi_comp_R, rho_basis);
152+
rho_basis->recip2real(phi_comp_G, phi_comp_R);
153153

154154
ModuleBase::matrix v_comp(GlobalV::NSPIN, rho_basis->nrxx);
155155
if (GlobalV::NSPIN == 4)
@@ -205,7 +205,7 @@ void surchem::test_V_to_N(ModuleBase::matrix &v,
205205
phi_comp_R[ir] = v(0, ir);
206206
}
207207

208-
GlobalC::UFFT.ToReciSpace(phi_comp_R, phi_comp_G, rho_basis);
208+
rho_basis->real2recip(phi_comp_R, phi_comp_G);
209209
for (int ig = 0; ig < rho_basis->npw; ig++)
210210
{
211211
if (rho_basis->gg[ig] >= 1.0e-12) // LiuXh 20180410
@@ -214,7 +214,7 @@ void surchem::test_V_to_N(ModuleBase::matrix &v,
214214
comp_reci[ig] = phi_comp_G[ig] / fac;
215215
}
216216
}
217-
GlobalC::UFFT.ToRealSpace(comp_reci, N_real, rho_basis);
217+
rho_basis->recip2real(comp_reci, N_real);
218218

219219
complex<double> *vloc_g = new complex<double>[rho_basis->npw];
220220
complex<double> *ng = new complex<double>[rho_basis->npw];
@@ -225,8 +225,7 @@ void surchem::test_V_to_N(ModuleBase::matrix &v,
225225
for (int ir = 0; ir < rho_basis->nrxx; ir++)
226226
Porter[ir] = rho[0][ir];
227227

228-
GlobalC::UFFT.ToReciSpace(GlobalC::pot.vltot,
229-
vloc_g, rho_basis); // now n is vloc in Recispace
228+
rho_basis->real2recip(GlobalC::pot.vltot,vloc_g);// now n is vloc in Recispace
230229
for (int ig = 0; ig < rho_basis->npw; ig++) {
231230
if (rho_basis->gg[ig] >= 1.0e-12) // LiuXh 20180410
232231
{
@@ -237,7 +236,7 @@ void surchem::test_V_to_N(ModuleBase::matrix &v,
237236
}
238237
}
239238
double *nr = new double[rho_basis->nrxx];
240-
GlobalC::UFFT.ToRealSpace(ng, nr, rho_basis);
239+
rho_basis->recip2real(ng, nr);
241240

242241
double *diff = new double[rho_basis->nrxx];
243242
double *diff2 = new double[rho_basis->nrxx];

‎source/module_surchem/cal_totn.cpp‎

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,8 +7,7 @@ void surchem::cal_totn(const UnitCell &cell, ModulePW::PW_Basis* rho_basis,
77
complex<double> *vloc_g = new complex<double>[rho_basis->npw];
88
ModuleBase::GlobalFunc::ZEROS(vloc_g, rho_basis->npw);
99

10-
GlobalC::UFFT.ToReciSpace(GlobalC::pot.vltot,
11-
vloc_g, rho_basis); // now n is vloc in Recispace
10+
rho_basis->real2recip(GlobalC::pot.vltot, vloc_g); // now n is vloc in Recispace
1211
for (int ig = 0; ig < rho_basis->npw; ig++) {
1312
if(ig==rho_basis->ig_gge0)
1413
{

‎source/module_surchem/cal_vcav.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ void shape_gradn(const complex<double> *PS_TOTN, ModulePW::PW_Basis* rho_basis,
3636

3737
double *PS_TOTN_real = new double[rho_basis->nrxx];
3838
ModuleBase::GlobalFunc::ZEROS(PS_TOTN_real, rho_basis->nrxx);
39-
GlobalC::UFFT.ToRealSpace(PS_TOTN, PS_TOTN_real,rho_basis);
39+
rho_basis->recip2real(PS_TOTN, PS_TOTN_real);
4040

4141
double epr_c = 1.0 / sqrt(ModuleBase::TWO_PI) / GlobalV::sigma_k;
4242
double epr_z = 0;
@@ -119,8 +119,8 @@ void surchem::createcavity(const UnitCell &ucell, ModulePW::PW_Basis* rho_basis,
119119
// packs the real array into a complex one
120120
// to G space
121121
complex<double> *inv_gn = new complex<double>[rho_basis->npw];
122-
GlobalC::UFFT.ToReciSpace(sqrt_nablan_2, inv_gn,rho_basis);
123-
122+
rho_basis->real2recip(sqrt_nablan_2, inv_gn);
123+
124124
// \nabla(1 / |\nabla n|), ggn in real space
125125
ModuleBase::Vector3<double> *ggn = new ModuleBase::Vector3<double>[rho_basis->nrxx];
126126
XC_Functional::grad_rho(inv_gn, ggn, rho_basis);

‎source/module_surchem/cal_vel.cpp‎

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

6060
// double *TOTN_real = new double[pwb.nrxx];
61-
GlobalC::UFFT.ToRealSpace(TOTN, TOTN_real,rho_basis);
61+
rho_basis->recip2real(TOTN, TOTN_real);
6262

6363
// -4pi * TOTN(G)
6464
complex<double> *B = new complex<double>[rho_basis->npw];
@@ -69,7 +69,7 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
6969

7070
// Build a nrxx vector to DO FFT .
7171
double *PS_TOTN_real = new double[rho_basis->nrxx];
72-
GlobalC::UFFT.ToRealSpace(PS_TOTN, PS_TOTN_real,rho_basis);
72+
rho_basis->recip2real(PS_TOTN, PS_TOTN_real);
7373

7474
// build epsilon in real space (nrxx)
7575
double *epsilon = new double[rho_basis->nrxx];
@@ -95,8 +95,8 @@ ModuleBase::matrix surchem::cal_vel(const UnitCell &cell,
9595
double *phi_tilda_R0 = new double[rho_basis->nrxx];
9696
// double *delta_phi_R = new double[pwb.nrxx];
9797

98-
GlobalC::UFFT.ToRealSpace(Sol_phi, phi_tilda_R,rho_basis);
99-
GlobalC::UFFT.ToRealSpace(Sol_phi0, phi_tilda_R0,rho_basis);
98+
rho_basis->recip2real(Sol_phi, phi_tilda_R);
99+
rho_basis->recip2real(Sol_phi0, phi_tilda_R0);
100100

101101
// the 1st item of tmp_Vel
102102
for (int i = 0; i < rho_basis->nrxx; i++)

‎source/module_surchem/minimize_cg.cpp‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -202,7 +202,7 @@ void surchem::Leps2(const UnitCell &ucell,
202202
{
203203
grad_grad_phi[ir] = grad_phi[ir].x;
204204
}
205-
GlobalC::UFFT.ToReciSpace(grad_grad_phi, grad_grad_phi_G, rho_basis);
205+
rho_basis->real2recip(grad_grad_phi, grad_grad_phi_G);
206206
XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis);
207207
for (int ir = 0; ir < rho_basis->nrxx; ir++)
208208
{
@@ -217,7 +217,7 @@ void surchem::Leps2(const UnitCell &ucell,
217217
{
218218
grad_grad_phi[ir] = grad_phi[ir].y;
219219
}
220-
GlobalC::UFFT.ToReciSpace(grad_grad_phi, grad_grad_phi_G, rho_basis);
220+
rho_basis->real2recip(grad_grad_phi, grad_grad_phi_G);
221221
XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis);
222222
for (int ir = 0; ir < rho_basis->nrxx; ir++)
223223
{
@@ -232,7 +232,7 @@ void surchem::Leps2(const UnitCell &ucell,
232232
{
233233
grad_grad_phi[ir] = grad_phi[ir].z;
234234
}
235-
GlobalC::UFFT.ToReciSpace(grad_grad_phi, grad_grad_phi_G, rho_basis);
235+
rho_basis->real2recip(grad_grad_phi, grad_grad_phi_G);
236236
XC_Functional::grad_rho(grad_grad_phi_G, tmp_vector3, rho_basis);
237237
for (int ir = 0; ir < rho_basis->nrxx; ir++)
238238
{
@@ -243,7 +243,7 @@ void surchem::Leps2(const UnitCell &ucell,
243243
// cout << lp_real << [i] << endl;
244244
// }
245245

246-
GlobalC::UFFT.ToReciSpace(lp_real, lp, rho_basis);
246+
rho_basis->real2recip(lp_real, lp);
247247
// cout<<"lp: "<<endl;
248248
// test_print(lp, 10);
249249

‎source/module_surchem/sol_force.cpp‎

Lines changed: 177 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,177 @@
1+
#include "surchem.h"
2+
#include "../module_base/timer.h"
3+
4+
void force_cor_one(const UnitCell &cell, ModulePW::PW_Basis* rho_basis , ModuleBase::matrix& forcesol)
5+
{
6+
7+
8+
//delta phi multiply by the derivative of nuclear charge density with respect to the positions
9+
std::complex<double> *N = new std::complex<double>[rho_basis->npw];
10+
std::complex<double> *vloc_at = new std::complex<double>[rho_basis->npw];
11+
std::complex<double> *delta_phi_g = new complex<double>[rho_basis->npw];
12+
//ModuleBase::GlobalFunc::ZEROS(delta_phi_g, rho_basis->npw);
13+
14+
rho_basis->real2recip(GlobalC::solvent_model.delta_phi, delta_phi_g);
15+
//GlobalC::UFFT.ToReciSpace(GlobalC::solvent_model.delta_phi, delta_phi_g,rho_basis);
16+
double Ael=0;double Ael1 = 0;
17+
//ModuleBase::GlobalFunc::ZEROS(vg, ngmc);
18+
int iat = 0;
19+
20+
for (int it = 0;it < cell.ntype;it++)
21+
{
22+
for (int ia = 0;ia < cell.atoms[it].na ; ia++)
23+
{
24+
for (int ig = 0; ig < rho_basis->npw; ig++)
25+
{
26+
complex<double> phase = exp( ModuleBase::NEG_IMAG_UNIT *ModuleBase::TWO_PI * ( rho_basis->gcar[ig] * cell.atoms[it].tau[ia]));
27+
//vloc for each atom
28+
vloc_at[ig] = GlobalC::ppcell.vloc(it, rho_basis->ig2igg[ig]) * phase;
29+
if(rho_basis->ig_gge0 == ig)
30+
{
31+
N[ig] = GlobalC::ucell.atoms[it].zv / GlobalC::ucell.omega;
32+
}
33+
else
34+
{
35+
const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI /
36+
(cell.tpiba2 * rho_basis->gg[ig]);
37+
38+
N[ig] = -vloc_at[ig] / fac;
39+
}
40+
41+
//force for each atom
42+
forcesol(iat, 0) += rho_basis->gcar[ig][0] * imag(conj(delta_phi_g[ig]) * N[ig]);
43+
forcesol(iat, 1) += rho_basis->gcar[ig][1] * imag(conj(delta_phi_g[ig]) * N[ig]);
44+
forcesol(iat, 2) += rho_basis->gcar[ig][2] * imag(conj(delta_phi_g[ig]) * N[ig]);
45+
}
46+
47+
forcesol(iat, 0) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
48+
forcesol(iat, 1) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
49+
forcesol(iat, 2) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
50+
//unit Ry/Bohr
51+
forcesol(iat, 0) *= 2 ;
52+
forcesol(iat, 1) *= 2 ;
53+
forcesol(iat, 2) *= 2 ;
54+
55+
//cout<<"Force1"<<iat<<":"<<" "<<forcesol(iat, 0)<<" "<<forcesol(iat, 1)<<" "<<forcesol(iat, 2)<<endl;
56+
57+
++iat;
58+
}
59+
}
60+
61+
delete[] vloc_at;
62+
delete[] N;
63+
delete[] delta_phi_g;
64+
65+
}
66+
67+
void force_cor_two(const UnitCell &cell, ModulePW::PW_Basis* rho_basis , ModuleBase::matrix& forcesol)
68+
{
69+
70+
complex<double> *n_pseudo = new complex<double>[rho_basis->npw];
71+
ModuleBase::GlobalFunc::ZEROS(n_pseudo,rho_basis->npw);
72+
73+
//GlobalC::solvent_model.gauss_charge(cell, pwb, n_pseudo);
74+
75+
double *Vcav_sum = new double[rho_basis->nrxx];
76+
ModuleBase::GlobalFunc::ZEROS(Vcav_sum, rho_basis->nrxx);
77+
std::complex<double> *Vcav_g = new complex<double>[rho_basis->npw];
78+
std::complex<double> *Vel_g = new complex<double>[rho_basis->npw];
79+
ModuleBase::GlobalFunc::ZEROS(Vcav_g, rho_basis->npw);
80+
ModuleBase::GlobalFunc::ZEROS(Vel_g, rho_basis->npw);
81+
for(int is=0; is<GlobalV::NSPIN; is++)
82+
{
83+
for (int ir=0; ir<rho_basis->nrxx; ir++)
84+
{
85+
Vcav_sum[ir] += GlobalC::solvent_model.Vcav(is, ir);
86+
}
87+
}
88+
89+
rho_basis->real2recip(Vcav_sum, Vcav_g);
90+
rho_basis->real2recip(GlobalC::solvent_model.epspot, Vel_g);
91+
92+
int iat = 0;
93+
double Ael1 = 0;
94+
for (int it = 0;it < cell.ntype;it++)
95+
{
96+
double RCS = GlobalC::solvent_model.GetAtom.atom_RCS[cell.atoms[it].psd];
97+
double sigma_rc_k = RCS / 2.5;
98+
for (int ia = 0;ia < cell.atoms[it].na;ia++)
99+
{
100+
//cell.atoms[0].tau[0].z = 3.302;
101+
//cout<<cell.atoms[it].tau[ia]<<endl;
102+
ModuleBase::GlobalFunc::ZEROS(n_pseudo, rho_basis->npw);
103+
for (int ig = 0; ig < rho_basis->npw; ig++)
104+
{
105+
// G^2
106+
double gg = rho_basis->gg[ig];
107+
gg = gg * cell.tpiba2;
108+
complex<double> phase = exp( ModuleBase::NEG_IMAG_UNIT *ModuleBase::TWO_PI * ( rho_basis->gcar[ig] * cell.atoms[it].tau[ia]));
109+
110+
n_pseudo[ig].real((GlobalC::solvent_model.GetAtom.atom_Z[cell.atoms[it].psd] - cell.atoms[it].zv) * phase.real()
111+
* exp(-0.5 * gg * (sigma_rc_k * sigma_rc_k)));
112+
n_pseudo[ig].imag((GlobalC::solvent_model.GetAtom.atom_Z[cell.atoms[it].psd] - cell.atoms[it].zv) * phase.imag()
113+
* exp(-0.5 * gg * (sigma_rc_k * sigma_rc_k)));
114+
}
115+
116+
for (int ig = 0; ig < rho_basis->npw; ig++)
117+
{
118+
n_pseudo[ig] /= cell.omega;
119+
}
120+
for (int ig = 0; ig < rho_basis->npw; ig++)
121+
{
122+
forcesol(iat, 0) -= rho_basis->gcar[ig][0] * imag(conj(Vcav_g[ig]+Vel_g[ig]) * n_pseudo[ig]);
123+
forcesol(iat, 1) -= rho_basis->gcar[ig][1] * imag(conj(Vcav_g[ig]+Vel_g[ig]) * n_pseudo[ig]);
124+
forcesol(iat, 2) -= rho_basis->gcar[ig][2] * imag(conj(Vcav_g[ig]+Vel_g[ig]) * n_pseudo[ig]);
125+
}
126+
127+
forcesol(iat, 0) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
128+
forcesol(iat, 1) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
129+
forcesol(iat, 2) *= (GlobalC::ucell.tpiba * GlobalC::ucell.omega);
130+
//eV/Ang
131+
forcesol(iat, 0) *= 2 ;
132+
forcesol(iat, 1) *= 2 ;
133+
forcesol(iat, 2) *= 2 ;
134+
135+
//cout<<"Force2"<<iat<<":"<<" "<<forcesol(iat, 0)<<" "<<forcesol(iat, 1)<<" "<<forcesol(iat, 2)<<endl;
136+
137+
++iat;
138+
}
139+
}
140+
141+
delete[] n_pseudo;
142+
delete[] Vcav_sum;
143+
delete[] Vcav_g;
144+
delete[] Vel_g;
145+
146+
}
147+
148+
void surchem::cal_force_sol(const UnitCell &cell, ModulePW::PW_Basis* rho_basis , ModuleBase::matrix& forcesol)
149+
{
150+
ModuleBase::TITLE("surchem", "cal_force_sol");
151+
ModuleBase::timer::tick("surchem", "cal_force_sol");
152+
153+
int nat = GlobalC::ucell.nat;
154+
ModuleBase::matrix force1(nat, 3);
155+
ModuleBase::matrix force2(nat, 3);
156+
157+
force_cor_one(cell, rho_basis,force1);
158+
force_cor_two(cell, rho_basis,force2);
159+
160+
int iat = 0;
161+
for (int it = 0;it < GlobalC::ucell.ntype;it++)
162+
{
163+
for (int ia = 0;ia < GlobalC::ucell.atoms[it].na;ia++)
164+
{
165+
for(int ipol = 0; ipol < 3; ipol++)
166+
{
167+
forcesol(iat, ipol) = 0.5*force1(iat, ipol) + force2 (iat, ipol);
168+
}
169+
170+
++iat;
171+
}
172+
}
173+
174+
Parallel_Reduce::reduce_double_pool(forcesol.c, forcesol.nr * forcesol.nc);
175+
ModuleBase::timer::tick("surchem", "cal_force_sol");
176+
return;
177+
}

‎source/module_surchem/surchem.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -100,7 +100,9 @@ class surchem
100100
ModuleBase::matrix v_compensating(const UnitCell &cell, ModulePW::PW_Basis *pwb);
101101

102102
void test_V_to_N(ModuleBase::matrix &v, const UnitCell &cell, ModulePW::PW_Basis *rho_basis, const double *const *const rho);
103-
103+
104+
void cal_force_sol(const UnitCell &cell, ModulePW::PW_Basis* rho_basis , ModuleBase::matrix& forcesol);
105+
104106
private:
105107
};
106108

0 commit comments

Comments
 (0)