1919#include < complex>
2020#include < vector>
2121
22- namespace ModuleDFPT {
22+ namespace ModuleDFPT
23+ {
2324
2425DFPT_HamiltShift::DFPT_HamiltShift (const UnitCell& ucell,
2526 ModulePW::PW_Basis* pw_rho,
2627 ModulePW::PW_Basis_K* pw_wfc,
2728 const std::vector<double >& veff_r,
2829 const DFPT_Pert* pert)
29- : ucell_(&ucell),
30- pw_rho_ (pw_rho),
31- pw_wfc_(pw_wfc),
32- pert_(pert),
33- veff_r_(veff_r),
34- tpiba2_(ucell.tpiba2),
35- nrxx_(pw_rho != nullptr ? pw_rho->nrxx : 0 ) {
36- for (int it = 0 ; it < ucell_->ntype ; ++it) {
30+ : ucell_(&ucell), pw_rho_(pw_rho), pw_wfc_(pw_wfc), pert_(pert), veff_r_(veff_r), tpiba2_(ucell.tpiba2),
31+ nrxx_ (pw_rho != nullptr ? pw_rho->nrxx : 0 )
32+ {
33+ for (int it = 0 ; it < ucell_->ntype ; ++it)
34+ {
3735 const pseudo& ncpp = ucell_->atoms [it].ncpp ;
38- if (ncpp.tvanp || ncpp.has_so ) {
36+ if (ncpp.tvanp || ncpp.has_so )
37+ {
3938 ModuleBase::WARNING_QUIT (" DFPT_HamiltShift" ,
4039 " the shifted Sternheimer operator is implemented for "
4140 " normal-conserving separable pseudopotentials only." );
@@ -45,10 +44,13 @@ DFPT_HamiltShift::DFPT_HamiltShift(const UnitCell& ucell,
4544 std::vector<int > ib;
4645 std::vector<int > m;
4746 int mu = 0 ;
48- for (int ibeta = 0 ; ibeta < ncpp.nbeta ; ++ibeta) {
47+ for (int ibeta = 0 ; ibeta < ncpp.nbeta ; ++ibeta)
48+ {
4949 const int l = ncpp.lll [ibeta];
50- for (int im = 0 ; im < 2 * l + 1 ; ++im) {
51- if (mu < ncpp.nh ) {
50+ for (int im = 0 ; im < 2 * l + 1 ; ++im)
51+ {
52+ if (mu < ncpp.nh )
53+ {
5254 ib.push_back (ibeta);
5355 m.push_back (im);
5456 }
@@ -60,29 +62,36 @@ DFPT_HamiltShift::DFPT_HamiltShift(const UnitCell& ucell,
6062 }
6163}
6264
63- DFPT_HamiltShift::~DFPT_HamiltShift () {}
65+ DFPT_HamiltShift::~DFPT_HamiltShift ()
66+ {
67+ }
6468
65- void DFPT_HamiltShift::set_context (const ModuleBase::Vector3<double >& q_cart, int k_idx) {
69+ void DFPT_HamiltShift::set_context (const ModuleBase::Vector3<double >& q_cart, int k_idx)
70+ {
6671 kq_.init (pw_wfc_, pw_rho_, q_cart, k_idx);
6772 ik_cache_ = k_idx;
6873 const int npw = kq_.get_npwk ();
6974
7075 // k+q G index -> charge-grid G index (both bases share the FFT cell)
7176 kq2rho_.assign (npw, -1 );
72- for (int igl = 0 ; igl < npw; ++igl) {
77+ for (int igl = 0 ; igl < npw; ++igl)
78+ {
7379 kq2rho_[igl] = kq_.get_ig_rho (igl);
7480 }
7581
7682 // cache the beta projectors of every atom on the k+q list
7783 std::vector<ModuleBase::Vector3<double >> gk (npw);
78- for (int igl = 0 ; igl < npw; ++igl) {
84+ for (int igl = 0 ; igl < npw; ++igl)
85+ {
7986 gk[igl] = kq_.get_gpluskq (igl);
8087 }
8188 vkb_cache_.assign (ucell_->nat , std::vector<std::vector<std::complex <double >>>());
82- for (int iat = 0 ; iat < ucell_->nat ; ++iat) {
89+ for (int iat = 0 ; iat < ucell_->nat ; ++iat)
90+ {
8391 const int it = ucell_->iat2it [iat];
8492 const int ia = ucell_->iat2ia [iat];
85- if (ucell_->atoms [it].ncpp .nh == 0 ) {
93+ if (ucell_->atoms [it].ncpp .nh == 0 )
94+ {
8695 continue ;
8796 }
8897 pert_->build_vkb (it, ia, gk, vkb_cache_[iat]);
@@ -92,120 +101,152 @@ void DFPT_HamiltShift::set_context(const ModuleBase::Vector3<double>& q_cart, in
92101 x_r_.assign (nrxx_, std::complex <double >(0.0 , 0.0 ));
93102}
94103
95- void DFPT_HamiltShift::set_shift (double shift) {
104+ void DFPT_HamiltShift::set_shift (double shift)
105+ {
96106 shift_ = shift;
97107}
98108
99- int DFPT_HamiltShift::dimension () const {
109+ int DFPT_HamiltShift::dimension () const
110+ {
100111 return kq_.get_npwk ();
101112}
102113
103- void DFPT_HamiltShift::apply (const std::complex <double >* x, std::complex <double >* y) const {
114+ void DFPT_HamiltShift::apply (const std::complex <double >* x, std::complex <double >* y) const
115+ {
104116 const int npw = kq_.get_npwk ();
105- if (npw <= 0 || x == nullptr || y == nullptr ) {
117+ if (npw <= 0 || x == nullptr || y == nullptr )
118+ {
106119 return ;
107120 }
108121 // kinetic part minus the eigenvalue shift
109- for (int igl = 0 ; igl < npw; ++igl) {
122+ for (int igl = 0 ; igl < npw; ++igl)
123+ {
110124 y[igl] = (tpiba2_ * kq_.get_gk2 (igl) - shift_) * x[igl];
111125 }
112126 // local effective potential: phase-free FFT convolution on the shared
113127 // grid (the k+q Bloch phases cancel in the product, real_space_dv conv.)
114128 std::fill (x_recip_.begin (), x_recip_.end (), std::complex <double >(0.0 , 0.0 ));
115- for (int igl = 0 ; igl < npw; ++igl) {
116- if (kq2rho_[igl] >= 0 ) {
129+ for (int igl = 0 ; igl < npw; ++igl)
130+ {
131+ if (kq2rho_[igl] >= 0 )
132+ {
117133 x_recip_[kq2rho_[igl]] = x[igl];
118134 }
119135 }
120136 pw_rho_->recip2real (x_recip_.data (), x_r_.data ());
121- for (int ir = 0 ; ir < nrxx_; ++ir) {
137+ for (int ir = 0 ; ir < nrxx_; ++ir)
138+ {
122139 x_r_[ir] *= veff_r_[ir];
123140 }
124141 pw_rho_->real2recip (x_r_.data (), x_recip_.data ());
125- for (int igl = 0 ; igl < npw; ++igl) {
126- if (kq2rho_[igl] >= 0 ) {
142+ for (int igl = 0 ; igl < npw; ++igl)
143+ {
144+ if (kq2rho_[igl] >= 0 )
145+ {
127146 y[igl] += x_recip_[kq2rho_[igl]];
128147 }
129148 }
130149 // nonlocal part with the cached k+q projectors
131- for (int iat = 0 ; iat < ucell_->nat ; ++iat) {
150+ for (int iat = 0 ; iat < ucell_->nat ; ++iat)
151+ {
132152 const int it = ucell_->iat2it [iat];
133153 const int nh = ucell_->atoms [it].ncpp .nh ;
134- if (nh == 0 ) {
154+ if (nh == 0 )
155+ {
135156 continue ;
136157 }
137158 const std::vector<std::vector<std::complex <double >>>& vkb = vkb_cache_[iat];
138159 becp_.assign (nh, std::complex <double >(0.0 , 0.0 ));
139- for (int mu = 0 ; mu < nh; ++mu) {
140- for (int igl = 0 ; igl < npw; ++igl) {
160+ for (int mu = 0 ; mu < nh; ++mu)
161+ {
162+ for (int igl = 0 ; igl < npw; ++igl)
163+ {
141164 becp_[mu] += std::conj (vkb[mu][igl]) * x[igl];
142165 }
143166 }
144167 dbecp_.assign (nh, std::complex <double >(0.0 , 0.0 ));
145- for (int mu = 0 ; mu < nh; ++mu) {
146- for (int nu = 0 ; nu < nh; ++nu) {
147- if (mu_m_[it][mu] != mu_m_[it][nu]) {
168+ for (int mu = 0 ; mu < nh; ++mu)
169+ {
170+ for (int nu = 0 ; nu < nh; ++nu)
171+ {
172+ if (mu_m_[it][mu] != mu_m_[it][nu])
173+ {
148174 continue ;
149175 }
150176 dbecp_[mu] += ucell_->atoms [it].ncpp .dion (mu_ib_[it][mu], mu_ib_[it][nu]) * becp_[nu];
151177 }
152178 }
153- for (int mu = 0 ; mu < nh; ++mu) {
154- for (int igl = 0 ; igl < npw; ++igl) {
179+ for (int mu = 0 ; mu < nh; ++mu)
180+ {
181+ for (int igl = 0 ; igl < npw; ++igl)
182+ {
155183 y[igl] += vkb[mu][igl] * dbecp_[mu];
156184 }
157185 }
158186 }
159187}
160188
161- double DFPT_HamiltShift::debug_t_vnl (const std::vector<std::complex <double >>& x) const {
189+ double DFPT_HamiltShift::debug_t_vnl (const std::vector<std::complex <double >>& x) const
190+ {
162191 const int npw = kq_.get_npwk ();
163192 double ekin = 0.0 ;
164- for (int igl = 0 ; igl < npw; ++igl) {
193+ for (int igl = 0 ; igl < npw; ++igl)
194+ {
165195 ekin += tpiba2_ * kq_.get_gk2 (igl) * std::norm (x[igl]);
166196 }
167197 double vnl = 0.0 ;
168- for (int iat = 0 ; iat < ucell_->nat ; ++iat) {
198+ for (int iat = 0 ; iat < ucell_->nat ; ++iat)
199+ {
169200 const int it = ucell_->iat2it [iat];
170201 const int nh = ucell_->atoms [it].ncpp .nh ;
171- if (nh == 0 ) {
202+ if (nh == 0 )
203+ {
172204 continue ;
173205 }
174206 const std::vector<std::vector<std::complex <double >>>& vkb = vkb_cache_[iat];
175207 becp_.assign (nh, std::complex <double >(0.0 , 0.0 ));
176- for (int mu = 0 ; mu < nh; ++mu) {
177- for (int igl = 0 ; igl < npw; ++igl) {
208+ for (int mu = 0 ; mu < nh; ++mu)
209+ {
210+ for (int igl = 0 ; igl < npw; ++igl)
211+ {
178212 becp_[mu] += std::conj (vkb[mu][igl]) * x[igl];
179213 }
180214 }
181215 dbecp_.assign (nh, std::complex <double >(0.0 , 0.0 ));
182- for (int mu = 0 ; mu < nh; ++mu) {
183- for (int nu = 0 ; nu < nh; ++nu) {
184- if (mu_m_[it][mu] != mu_m_[it][nu]) {
216+ for (int mu = 0 ; mu < nh; ++mu)
217+ {
218+ for (int nu = 0 ; nu < nh; ++nu)
219+ {
220+ if (mu_m_[it][mu] != mu_m_[it][nu])
221+ {
185222 continue ;
186223 }
187224 dbecp_[mu] += ucell_->atoms [it].ncpp .dion (mu_ib_[it][mu], mu_ib_[it][nu]) * becp_[nu];
188225 }
189226 }
190- for (int mu = 0 ; mu < nh; ++mu) {
227+ for (int mu = 0 ; mu < nh; ++mu)
228+ {
191229 vnl += std::real (std::conj (becp_[mu]) * dbecp_[mu]);
192230 }
193231 }
194232 return ekin + vnl;
195233}
196234
197- double DFPT_HamiltShift::debug_v_wfc (const std::vector<std::complex <double >>& x) const {
235+ double DFPT_HamiltShift::debug_v_wfc (const std::vector<std::complex <double >>& x) const
236+ {
198237 const int npw = kq_.get_npwk ();
199238 std::vector<std::complex <double >> ur (nrxx_, std::complex <double >(0.0 , 0.0 ));
200239 pw_wfc_->recip2real (x.data (), ur.data (), ik_cache_);
201- for (int ir = 0 ; ir < nrxx_; ++ir) {
240+ for (int ir = 0 ; ir < nrxx_; ++ir)
241+ {
202242 ur[ir] *= veff_r_[ir];
203243 }
204244 std::vector<std::complex <double >> xg (pw_wfc_->npwk [ik_cache_], std::complex <double >(0.0 , 0.0 ));
205245 pw_wfc_->real2recip (ur.data (), xg.data (), ik_cache_);
206246 std::complex <double > dot (0.0 , 0.0 );
207247 const int n = std::min (static_cast <int >(xg.size ()), npw);
208- for (int igl = 0 ; igl < n; ++igl) {
248+ for (int igl = 0 ; igl < n; ++igl)
249+ {
209250 dot += std::conj (x[igl]) * xg[igl];
210251 }
211252 return dot.real ();
0 commit comments