Skip to content

Commit 3aaf22e

Browse files
sunml99sunml99
committed
add force correction of implicit solvation model
1 parent 815614a commit 3aaf22e

File tree

9 files changed

+228
-9
lines changed

9 files changed

+228
-9
lines changed

‎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+
force.cpp
1415
)

‎source/module_surchem/force.cpp‎

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

‎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

‎source/src_lcao/FORCE_STRESS.cpp‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66
#include "../src_pw/vdwd3.h"
77
#include "../module_base/timer.h"
88
#include "../module_surchem/efield.h" // liuyu add 2022-05-18
9+
#include "../module_surchem/surchem.h" //sunml add 2022-08-10
910
#ifdef __DEEPKS
1011
#include "../module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03
1112
#endif
@@ -183,6 +184,13 @@ void Force_Stress_LCAO::getForceStress(
183184
{
184185
fefield.create(nat, 3);
185186
Efield::compute_force(GlobalC::ucell, fefield);
187+
}
188+
//Force from implicit solvation model
189+
ModuleBase::matrix fsol;
190+
if(GlobalV::imp_sol&&isforce)
191+
{
192+
fsol.create(nat, 3);
193+
GlobalC::solvent_model.cal_force_sol(GlobalC::ucell,GlobalC::rhopw,fsol);
186194
}
187195
//Force contribution from DFT+U
188196
ModuleBase::matrix force_dftu;
@@ -255,6 +263,11 @@ void Force_Stress_LCAO::getForceStress(
255263
{
256264
fcs(iat, i) += fefield(iat, i);
257265
}
266+
//implicit solvation model
267+
if(GlobalV::imp_sol)
268+
{
269+
fcs(iat, i) += fsol(iat, i);
270+
}
258271
#ifdef __DEEPKS
259272
// mohan add 2021-08-04
260273
if (GlobalV::deepks_scf)
@@ -379,6 +392,11 @@ void Force_Stress_LCAO::getForceStress(
379392
f_pw.print("EFIELD FORCE", fefield,0);
380393
//this->print_force("EFIELD FORCE",fefield,1,ry);
381394
}
395+
if(GlobalV::imp_sol)
396+
{
397+
f_pw.print("IMP_SOL FORCE", fsol,0);
398+
//this->print_force("IMP_SOL FORCE",fsol,1,ry);
399+
}
382400
if(GlobalC::vdwd2_para.flag_vdwd2||GlobalC::vdwd3_para.flag_vdwd3)
383401
{
384402
f_pw.print("VDW FORCE", force_vdw,0);

‎source/src_pw/forces.cpp‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
#include "../src_parallel/parallel_reduce.h"
1010
#include "../module_base/timer.h"
1111
#include "../module_surchem/efield.h"
12+
#include "../module_surchem/surchem.h"
1213

1314
double Forces::output_acc = 1.0e-8; // (Ryd/angstrom).
1415

@@ -81,6 +82,17 @@ void Forces::init(ModuleBase::matrix& force, const psi::Psi<std::complex<double>
8182
}
8283
}
8384

85+
ModuleBase::matrix forcesol;
86+
if (GlobalV::imp_sol)
87+
{
88+
forcesol.create(GlobalC::ucell.nat, 3);
89+
GlobalC::solvent_model.cal_force_sol(GlobalC::ucell, GlobalC::rhopw, forcesol);
90+
if(GlobalV::TEST_FORCE)
91+
{
92+
Forces::print("IMP_SOL FORCE (Ry/Bohr)", forcesol);
93+
}
94+
}
95+
8496
//impose total force = 0
8597
int iat = 0;
8698
for (int ipol = 0; ipol < 3; ipol++)
@@ -109,6 +121,11 @@ void Forces::init(ModuleBase::matrix& force, const psi::Psi<std::complex<double>
109121
force(iat,ipol) = force(iat, ipol) + force_e(iat, ipol);
110122
}
111123

124+
if(GlobalV::imp_sol)
125+
{
126+
force(iat,ipol) = force(iat, ipol) + forcesol(iat, ipol);
127+
}
128+
112129
sum += force(iat, ipol);
113130

114131
iat++;
@@ -213,6 +230,7 @@ void Forces::init(ModuleBase::matrix& force, const psi::Psi<std::complex<double>
213230
Forces::print("ION FORCE (eV/Angstrom)", forceion,0);
214231
Forces::print("SCC FORCE (eV/Angstrom)", forcescc,0);
215232
if(GlobalV::EFIELD_FLAG) Forces::print("EFIELD FORCE (eV/Angstrom)", force_e,0);
233+
if(GlobalV::imp_sol) Forces::print("IMP_SOL FORCE (eV/Angstrom)", forcesol,0);
216234
}
217235
Forces::print(" TOTAL-FORCE (eV/Angstrom)", force,0);
218236

‎tests/integrate/115_PW_sol_H2O/INPUT‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ ntype 2
66
nbands 20
77
calculation scf
88
basis_type pw
9+
cal_force 1
910

1011
#Parameters (Accuracy)
1112
ecutwfc 20

‎tests/integrate/115_PW_sol_H2O/result.ref‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
1-
etotref -443.1237068197200
1+
etotref -443.1237068197164604
22
etotperatomref -147.7079022732
3-
esolelref -1.02249508527
4-
esolcavref +0.0325034001929
5-
totaltimeref 10.112
3+
totalforceref 16.865518
4+
esolelref -1.02249508526
5+
esolcavref +0.0325034001928
6+
totaltimeref 44.58741

‎tests/integrate/215_NO_sol_H2O/INPUT‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@ ntype 2
77
nbands 20
88
calculation scf
99
basis_type pw
10+
cal_force 1
1011

1112
#Parameters (Accuracy)
1213
ecutwfc 20

‎tests/integrate/215_NO_sol_H2O/result.ref‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
1-
etotref -443.1237068197200
1+
etotref -443.1237068197164604
22
etotperatomref -147.7079022732
3-
esolelref -1.02249508527
4-
esolcavref +0.0325034001929
5-
totaltimeref 10.042
3+
totalforceref 16.865518
4+
esolelref -1.02249508526
5+
esolcavref +0.0325034001928
6+
totaltimeref 43.79188

0 commit comments

Comments
 (0)