Skip to content

Commit bfb00c0

Browse files
zgn-26714zgn-26714
committed
Reverted the process of removing PAW.
1 parent 480a7ef commit bfb00c0

File tree

1 file changed

+197
-9
lines changed

1 file changed

+197
-9
lines changed

‎source/module_hamilt_pw/hamilt_pwdft/forces.cpp‎

Lines changed: 197 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,9 @@
2121
#include <omp.h>
2222
#endif
2323

24+
#ifdef USE_PAW
25+
#include "module_cell/module_paw/paw_cell.h"
26+
#endif
2427

2528
template <typename FPTYPE, typename Device>
2629
void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
@@ -53,23 +56,108 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
5356
ModuleBase::matrix forceonsite(nat, 3);
5457

5558
// Force due to local ionic potential
56-
this->cal_force_loc(ucell,forcelc, rho_basis, locpp->vloc, chr);
59+
// For PAW, calculated together in paw_cell.calculate_force
60+
if (!PARAM.inp.use_paw)
61+
{
62+
this->cal_force_loc(ucell,forcelc, rho_basis, locpp->vloc, chr);
63+
}
64+
else
65+
{
66+
forcelc.zero_out();
67+
}
5768

5869
// Ewald
5970
this->cal_force_ew(ucell,forceion, rho_basis, p_sf);
6071

6172
// Force due to nonlocal part of pseudopotential
6273
if (wfc_basis != nullptr)
6374
{
64-
65-
this->npwx = wfc_basis->npwk_max;
66-
Forces::cal_force_nl(forcenl, wg, ekb, pkv, wfc_basis, p_sf, *p_nlpp, ucell, psi_in);
67-
68-
if (PARAM.globalv.use_uspp)
75+
if (!PARAM.inp.use_paw)
6976
{
70-
this->cal_force_us(forcenl, rho_basis, *p_nlpp, elec, ucell);
77+
this->npwx = wfc_basis->npwk_max;
78+
Forces::cal_force_nl(forcenl, wg, ekb, pkv, wfc_basis, p_sf, *p_nlpp, ucell, psi_in);
79+
if (PARAM.globalv.use_uspp)
80+
{
81+
this->cal_force_us(forcenl, rho_basis, *p_nlpp, elec, ucell);
82+
}
7183
}
84+
else
85+
{
86+
#ifdef USE_PAW
87+
for (int ik = 0; ik < wfc_basis->nks; ik++)
88+
{
89+
const int npw = wfc_basis->npwk[ik];
90+
ModuleBase::Vector3<double>* _gk = new ModuleBase::Vector3<double>[npw];
91+
for (int ig = 0; ig < npw; ig++)
92+
{
93+
_gk[ig] = wfc_basis->getgpluskcar(ik, ig);
94+
}
95+
double* kpt;
96+
kpt = new double[3];
97+
kpt[0] = wfc_basis->kvec_c[ik].x;
98+
kpt[1] = wfc_basis->kvec_c[ik].y;
99+
kpt[2] = wfc_basis->kvec_c[ik].z;
100+
101+
double** kpg;
102+
double** gcar;
103+
kpg = new double*[npw];
104+
gcar = new double*[npw];
105+
for (int ipw = 0; ipw < npw; ipw++)
106+
{
107+
kpg[ipw] = new double[3];
108+
kpg[ipw][0] = _gk[ipw].x;
109+
kpg[ipw][1] = _gk[ipw].y;
110+
kpg[ipw][2] = _gk[ipw].z;
111+
112+
gcar[ipw] = new double[3];
113+
gcar[ipw][0] = wfc_basis->getgcar(ik, ipw).x;
114+
gcar[ipw][1] = wfc_basis->getgcar(ik, ipw).y;
115+
gcar[ipw][2] = wfc_basis->getgcar(ik, ipw).z;
116+
}
117+
118+
GlobalC::paw_cell.set_paw_k(npw,
119+
wfc_basis->npwk_max,
120+
kpt,
121+
wfc_basis->get_ig2ix(ik).data(),
122+
wfc_basis->get_ig2iy(ik).data(),
123+
wfc_basis->get_ig2iz(ik).data(),
124+
(const double**)kpg,
125+
ucell.tpiba,
126+
(const double**)gcar);
127+
128+
delete[] kpt;
129+
for (int ipw = 0; ipw < npw; ipw++)
130+
{
131+
delete[] kpg[ipw];
132+
delete[] gcar[ipw];
133+
}
134+
delete[] kpg;
135+
delete[] gcar;
136+
137+
GlobalC::paw_cell.get_vkb();
138+
139+
GlobalC::paw_cell.set_currentk(ik);
72140

141+
psi_in[0].fix_k(ik);
142+
double *weight, *epsilon;
143+
weight = new double[PARAM.inp.nbands];
144+
epsilon = new double[PARAM.inp.nbands];
145+
for (int ib = 0; ib < PARAM.inp.nbands; ib++)
146+
{
147+
weight[ib] = wg(ik, ib);
148+
epsilon[ib] = ekb(ik, ib);
149+
}
150+
GlobalC::paw_cell.paw_nl_force(reinterpret_cast<std::complex<double>*>(psi_in[0].get_pointer()),
151+
epsilon,
152+
weight,
153+
PARAM.inp.nbands,
154+
forcenl.c);
155+
156+
delete[] weight;
157+
delete[] epsilon;
158+
}
159+
#endif
160+
}
73161
// DFT+U and DeltaSpin
74162
if(PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch)
75163
{
@@ -78,10 +166,26 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
78166
}
79167

80168
// non-linear core correction
81-
Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, ucell);
169+
// not relevant for PAW
170+
if (!PARAM.inp.use_paw)
171+
{
172+
Forces::cal_force_cc(forcecc, rho_basis, chr, locpp->numeric, ucell);
173+
}
174+
else
175+
{
176+
forcecc.zero_out();
177+
}
82178

83179
// force due to core charge
84-
this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, ucell);
180+
// For PAW, calculated together in paw_cell.calculate_force
181+
if (!PARAM.inp.use_paw)
182+
{
183+
this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, ucell);
184+
}
185+
else
186+
{
187+
forcescc.zero_out();
188+
}
85189

86190
ModuleBase::matrix stress_vdw_pw; //.create(3,3);
87191
ModuleBase::matrix force_vdw;
@@ -134,6 +238,52 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
134238
ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL FORCE (Ry/Bohr)", forcesol);
135239
}
136240
}
241+
#ifdef USE_PAW
242+
if (PARAM.inp.use_paw)
243+
{
244+
double* force_paw;
245+
double* rhor;
246+
rhor = new double[rho_basis->nrxx];
247+
for (int ir = 0; ir < rho_basis->nrxx; ir++)
248+
{
249+
rhor[ir] = 0.0;
250+
}
251+
for (int is = 0; is < PARAM.inp.nspin; is++)
252+
{
253+
for (int ir = 0; ir < rho_basis->nrxx; ir++)
254+
{
255+
rhor[ir] += chr->rho[is][ir] + chr->nhat[is][ir];
256+
}
257+
}
258+
259+
force_paw = new double[3 * this->nat];
260+
ModuleBase::matrix v_xc, v_effective;
261+
v_effective.create(PARAM.inp.nspin, rho_basis->nrxx);
262+
v_effective.zero_out();
263+
elec.pot->update_from_charge(elec.charge, &ucell);
264+
v_effective = elec.pot->get_effective_v();
265+
266+
v_xc.create(PARAM.inp.nspin, rho_basis->nrxx);
267+
v_xc.zero_out();
268+
const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v
269+
= XC_Functional::v_xc(rho_basis->nrxx, elec.charge, &ucell);
270+
v_xc = std::get<2>(etxc_vtxc_v);
271+
272+
GlobalC::paw_cell.calculate_force(v_effective.c, v_xc.c, rhor, force_paw);
273+
274+
for (int iat = 0; iat < this->nat; iat++)
275+
{
276+
// Ha to Ry
277+
forcepaw(iat, 0) = force_paw[3 * iat] * 2.0;
278+
forcepaw(iat, 1) = force_paw[3 * iat + 1] * 2.0;
279+
forcepaw(iat, 2) = force_paw[3 * iat + 2] * 2.0;
280+
}
281+
282+
delete[] force_paw;
283+
delete[] rhor;
284+
}
285+
#endif
286+
137287

138288
// impose total force = 0
139289
int iat = 0;
@@ -149,6 +299,11 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
149299
force(iat, ipol) = forcelc(iat, ipol) + forceion(iat, ipol) + forcenl(iat, ipol) + forcecc(iat, ipol)
150300
+ forcescc(iat, ipol);
151301

302+
if (PARAM.inp.use_paw)
303+
{
304+
force(iat, ipol) += forcepaw(iat, ipol);
305+
}
306+
152307
if (vdw_solver != nullptr) // linpz and jiyy added vdw force, modified by zhengdy
153308
{
154309
force(iat, ipol) += force_vdw(iat, ipol);
@@ -287,6 +442,15 @@ void Forces<FPTYPE, Device>::cal_force(UnitCell& ucell,
287442
ModuleIO::print_force(GlobalV::ofs_running, ucell, "ION FORCE (eV/Angstrom)", forceion, false);
288443
ModuleIO::print_force(GlobalV::ofs_running, ucell, "SCC FORCE (eV/Angstrom)", forcescc, false);
289444

445+
if (PARAM.inp.use_paw)
446+
{
447+
ModuleIO::print_force(GlobalV::ofs_running,
448+
ucell,
449+
"PAW FORCE (eV/Angstrom)",
450+
forcepaw,
451+
false);
452+
}
453+
290454
if (PARAM.inp.efield_flag)
291455
{
292456
ModuleIO::print_force(GlobalV::ofs_running, ucell, "EFIELD FORCE (eV/Angstrom)", force_e, false);
@@ -504,6 +668,13 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
504668
if (ucell.atoms[it].na != 0)
505669
{
506670
double dzv;
671+
if (PARAM.inp.use_paw)
672+
{
673+
#ifdef USE_PAW
674+
dzv = GlobalC::paw_cell.get_val(it);
675+
#endif
676+
}
677+
else
507678
{
508679
dzv = ucell.atoms[it].ncpp.zv;
509680
}
@@ -520,6 +691,13 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
520691
double charge = 0.0;
521692
for (int it = 0; it < ucell.ntype; it++)
522693
{
694+
if (PARAM.inp.use_paw)
695+
{
696+
#ifdef USE_PAW
697+
charge += ucell.atoms[it].na * GlobalC::paw_cell.get_val(it);
698+
#endif
699+
}
700+
else
523701
{
524702
charge += ucell.atoms[it].na * ucell.atoms[it].ncpp.zv; // mohan modify 2007-11-7
525703
}
@@ -687,6 +865,16 @@ void Forces<FPTYPE, Device>::cal_force_ew(const UnitCell& ucell,
687865
const double rr = sqrt(r2[n]) * ucell.lat0;
688866

689867
double factor;
868+
if (PARAM.inp.use_paw)
869+
{
870+
#ifdef USE_PAW
871+
factor = GlobalC::paw_cell.get_val(T1) * GlobalC::paw_cell.get_val(T2) * ModuleBase::e2
872+
/ (rr * rr)
873+
* (erfc(sqa * rr) / rr + sq8a_2pi * ModuleBase::libm::exp(-alpha * rr * rr))
874+
* ucell.lat0;
875+
#endif
876+
}
877+
else
690878
{
691879
factor = ucell.atoms[T1].ncpp.zv * ucell.atoms[T2].ncpp.zv
692880
* ModuleBase::e2 / (rr * rr)

0 commit comments

Comments
 (0)