Skip to content

Commit 67a0e88

Browse files
lyb9812claude
lyb9812
and
claude
authored
Refactor: clean up and fix TDOFDFT evolve_ofdft module (#7389)
* Remove debug output in renormalize_psi Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Fix: CD potential now applied to correct spin channel instead of always spin 0 Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Refactor: Replace raw new/delete with std::vector in cal_vw_potential_phi for automatic memory management Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Refactor: Replace raw new/delete with std::vector in cal_CD_potential for automatic memory management Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Refactor: Replace abs(x)*abs(x) with std::norm for clarity and consistency with RK2 Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Refactor: Remove unused <iostream> include Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> * Refactor: Remove dead nspin <= 0 checks that can never trigger Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 4.7 <noreply@anthropic.com>
1 parent 91c0c03 commit 67a0e88

1 file changed

Lines changed: 22 additions & 57 deletions

File tree

‎source/source_pw/module_ofdft/evolve_ofdft.cpp‎

Lines changed: 22 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
#include "evolve_ofdft.h"
22

33
#include "source_io/module_parameter/parameter.h"
4-
#include <iostream>
4+
#include <complex>
55

66
#include "source_base/parallel_reduce.h"
77

@@ -29,7 +29,7 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec,
2929
{
3030
for (int ir = 0; ir < nrxx; ++ir)
3131
{
32-
chr.rho[is][ir] = abs(psi_[is * nrxx + ir])*abs(psi_[is * nrxx + ir]);
32+
chr.rho[is][ir] = std::norm(psi_[is * nrxx + ir]);
3333
}
3434
}
3535
this->renormalize_psi(chr, pw_rho, psi_);
@@ -62,7 +62,7 @@ void Evolve_OFDFT::renormalize_psi(Charge& chr, ModulePW::PW_Basis* pw_rho, std:
6262
const int nspin = PARAM.inp.nspin;
6363
const int nrxx = pw_rho->nrxx;
6464

65-
std::cout<<"sr="<<sr<<" nelec="<<PARAM.inp.nelec<<" normalize_factor="<<normalize_factor<<std::endl;
65+
6666
#ifdef _OPENMP
6767
#pragma omp parallel for collapse(2)
6868
#endif
@@ -118,51 +118,33 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector<std::complex<double>>& pphi,
118118
const double* gg = pw_rho->gg;
119119
const double tpiba2 = pw_rho->tpiba2;
120120

121-
if (nspin <= 0) {
122-
ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive");
123-
}
124-
std::complex<double>** rLapPhi = new std::complex<double>*[nspin];
125-
#ifdef _OPENMP
126-
#pragma omp parallel for
127-
#endif
128-
for (int is = 0; is < nspin; ++is) {
129-
rLapPhi[is] = new std::complex<double>[nrxx];
121+
std::vector<std::vector<std::complex<double>>> rLapPhi(nspin, std::vector<std::complex<double>>(nrxx));
122+
std::vector<std::vector<std::complex<double>>> recipPhi(nspin, std::vector<std::complex<double>>(npw));
123+
124+
for (int is = 0; is < nspin; ++is)
125+
{
130126
for (int ir = 0; ir < nrxx; ++ir)
131127
{
132-
rLapPhi[is][ir]=pphi[is * nrxx + ir];
128+
rLapPhi[is][ir] = pphi[is * nrxx + ir];
133129
}
134130
}
135-
std::complex<double>** recipPhi = new std::complex<double>*[nspin];
136131

137132
#ifdef _OPENMP
138133
#pragma omp parallel for
139-
#endif
134+
#endif
140135
for (int is = 0; is < nspin; ++is)
141136
{
142-
recipPhi[is] = new std::complex<double>[npw];
143-
144-
pw_rho->real2recip(rLapPhi[is], recipPhi[is]);
137+
pw_rho->real2recip(rLapPhi[is].data(), recipPhi[is].data());
145138
for (int ik = 0; ik < npw; ++ik)
146139
{
147140
recipPhi[is][ik] *= gg[ik] * tpiba2;
148141
}
149-
pw_rho->recip2real(recipPhi[is], rLapPhi[is]);
142+
pw_rho->recip2real(recipPhi[is].data(), rLapPhi[is].data());
150143
for (int ir = 0; ir < nrxx; ++ir)
151144
{
152145
Hpsi[is * nrxx + ir] += rLapPhi[is][ir];
153146
}
154147
}
155-
156-
#ifdef _OPENMP
157-
#pragma omp parallel for
158-
#endif
159-
for (int is = 0; is < nspin; ++is)
160-
{
161-
delete[] recipPhi[is];
162-
delete[] rLapPhi[is];
163-
}
164-
delete[] recipPhi;
165-
delete[] rLapPhi;
166148
}
167149

168150
void Evolve_OFDFT::cal_CD_potential(std::vector<std::complex<double>>& psi_,
@@ -178,19 +160,14 @@ void Evolve_OFDFT::cal_CD_potential(std::vector<std::complex<double>>& psi_,
178160
const double* gg = pw_rho->gg;
179161
const ModuleBase::Vector3<double>* gcar = pw_rho->gcar;
180162

181-
if (nspin <= 0) {
182-
ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive");
183-
}
184-
std::complex<double>** recipPhi = new std::complex<double>*[nspin];
185-
std::complex<double>** rPhi = new std::complex<double>*[nspin];
186-
#ifdef _OPENMP
187-
#pragma omp parallel for
188-
#endif
189-
for (int is = 0; is < nspin; ++is) {
190-
rPhi[is] = new std::complex<double>[nrxx];
163+
std::vector<std::vector<std::complex<double>>> recipPhi(nspin, std::vector<std::complex<double>>(npw));
164+
std::vector<std::vector<std::complex<double>>> rPhi(nspin, std::vector<std::complex<double>>(nrxx));
165+
166+
for (int is = 0; is < nspin; ++is)
167+
{
191168
for (int ir = 0; ir < nrxx; ++ir)
192169
{
193-
rPhi[is][ir]=psi_[is * nrxx + ir];
170+
rPhi[is][ir] = psi_[is * nrxx + ir];
194171
}
195172
}
196173

@@ -208,14 +185,13 @@ void Evolve_OFDFT::cal_CD_potential(std::vector<std::complex<double>>& psi_,
208185
std::vector<std::complex<double>> rCurrent_z(nrxx);
209186
std::vector<std::complex<double>> kF_r(nrxx);
210187
std::vector<std::complex<double>> rCDPotential(nrxx);
211-
recipPhi[is] = new std::complex<double>[npw];
212188

213189
for (int ir = 0; ir < nrxx; ++ir)
214190
{
215191
kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1.0/3.0);
216192
}
217193

218-
pw_rho->real2recip(rPhi[is], recipPhi[is]);
194+
pw_rho->real2recip(rPhi[is].data(), recipPhi[is].data());
219195
for (int ik = 0; ik < npw; ++ik)
220196
{
221197
recipCurrent_x[ik]=imag*gcar[ik].x*recipPhi[is][ik]*tpiba;
@@ -250,25 +226,14 @@ void Evolve_OFDFT::cal_CD_potential(std::vector<std::complex<double>>& psi_,
250226

251227
for (int ir = 0; ir < nrxx; ++ir)
252228
{
253-
rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3)
229+
rpot(is, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3)
254230
/ (2.0*std::pow(std::real(kF_r[ir]),2));
255-
if (std::isnan(rpot(0, ir)))
231+
if (std::isnan(rpot(is, ir)))
256232
{
257-
rpot(0, ir)=0.0;
233+
rpot(is, ir)=0.0;
258234
}
259235
}
260236
}
261-
262-
#ifdef _OPENMP
263-
#pragma omp parallel for
264-
#endif
265-
for (int is = 0; is < nspin; ++is)
266-
{
267-
delete[] recipPhi[is];
268-
delete[] rPhi[is];
269-
}
270-
delete[] recipPhi;
271-
delete[] rPhi;
272237
}
273238

274239
void Evolve_OFDFT::propagate_psi_RK4(elecstate::ElecState* pelec,

0 commit comments

Comments
 (0)