33#include " source_io/module_parameter/parameter.h"
44#include " surchem.h"
55
6- void shape_gradn (const double * PS_TOTN_real, const ModulePW::PW_Basis* rho_basis, double * eprime)
6+
7+ static void local_grad_rho_vel (const UnitCell& ucell,
8+ const ModulePW::PW_Basis* rho_basis,
9+ const std::complex <double >* rho_G, // G
10+ ModuleBase::Vector3<double >* grad_R, // R
11+ std::complex <double >* aux_G, //
12+ double * aux_R) //
713{
14+ for (int i = 0 ; i < 3 ; ++i)
15+ {
16+ // 1. iG * rho(G)
17+ for (int ig = 0 ; ig < rho_basis->npw ; ig++) {
18+ aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar [ig][i];
19+ }
820
21+ // 2. FFT: G -> R
22+ rho_basis->recip2real (aux_G, aux_R);
23+
24+ // 3. 2pi/a
25+ for (int ir = 0 ; ir < rho_basis->nrxx ; ir++) {
26+ grad_R[ir][i] = aux_R[ir] * ucell.tpiba ;
27+ }
28+ }
29+ }
30+
31+ void shape_gradn (const double * PS_TOTN_real, const ModulePW::PW_Basis* rho_basis, double * eprime)
32+ {
933 double epr_c = 1.0 / sqrt (ModuleBase::TWO_PI ) / PARAM .inp .sigma_k ;
1034 double epr_z = 0 ;
1135 double min = 1e-10 ;
@@ -16,6 +40,25 @@ void shape_gradn(const double* PS_TOTN_real, const ModulePW::PW_Basis* rho_basis
1640 }
1741}
1842
43+ static void local_grad_rho_vel_tpiba (const double tpiba,
44+ const ModulePW::PW_Basis* rho_basis,
45+ const std::complex <double >* rho_G,
46+ ModuleBase::Vector3<double >* grad_R,
47+ std::complex <double >* aux_G,
48+ double * aux_R)
49+ {
50+ for (int i = 0 ; i < 3 ; ++i)
51+ {
52+ for (int ig = 0 ; ig < rho_basis->npw ; ig++) {
53+ aux_G[ig] = ModuleBase::IMAG_UNIT * rho_G[ig] * rho_basis->gcar [ig][i];
54+ }
55+ rho_basis->recip2real (aux_G, aux_R);
56+ for (int ir = 0 ; ir < rho_basis->nrxx ; ir++) {
57+ grad_R[ir][i] = aux_R[ir] * tpiba;
58+ }
59+ }
60+ }
61+
1962void eps_pot (const double * PS_TOTN_real,
2063 const double & tpiba,
2164 const std::complex <double >* phi,
@@ -36,8 +79,21 @@ void eps_pot(const double* PS_TOTN_real,
3679 ModuleBase::Vector3<double > *nabla_phi = new ModuleBase::Vector3<double >[rho_basis->nrxx ];
3780 double *phisq = new double [rho_basis->nrxx ];
3881
39- // nabla phi
40- XC_Functional::grad_rho (phi, nabla_phi, rho_basis, tpiba);
82+
83+
84+ std::complex <double > *aux_G = new std::complex <double >[rho_basis->npw ];
85+ double *aux_R = new double [rho_basis->nrxx ];
86+
87+
88+ local_grad_rho_vel_tpiba (tpiba, rho_basis, phi, nabla_phi, aux_G, aux_R);
89+
90+
91+ delete[] aux_G;
92+ delete[] aux_R;
93+
94+ // ========================================================
95+ // 修复结束
96+ // ========================================================
4197
4298 for (int ir = 0 ; ir < rho_basis->nrxx ; ir++)
4399 {
@@ -120,6 +176,7 @@ void surchem::cal_vel(const UnitCell& cell,
120176 this ->Ael *= cell.omega / rho_basis->nxyz ;
121177
122178 // the 2nd item of tmp_Vel
179+ // eps_pot internally calls shape_gradn and now the SAFE local_grad_rho_vel_tpiba
123180 eps_pot (PS_TOTN_real, cell.tpiba , Sol_phi, rho_basis, epsilon, epspot);
124181
125182 for (int i = 0 ; i < rho_basis->nrxx ; i++)
@@ -162,4 +219,4 @@ void surchem::cal_vel(const UnitCell& cell,
162219
163220 ModuleBase::timer::tick (" surchem" , " cal_vel" );
164221 return ;
165- }
222+ }
0 commit comments