diff --git a/source/source_esolver/esolver_of_tddft.cpp b/source/source_esolver/esolver_of_tddft.cpp index a978a1267c..daeda628cb 100644 --- a/source/source_esolver/esolver_of_tddft.cpp +++ b/source/source_esolver/esolver_of_tddft.cpp @@ -16,6 +16,8 @@ //-----stress------------------ #include "source_pw/module_ofdft/of_stress_pw.h" +#include + namespace ModuleESolver { @@ -45,12 +47,13 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) this->iter_time = std::chrono::system_clock::now(); #endif - if (istep==0) + if (this->phi_td.empty()) { - this->phi_td.resize(PARAM.inp.nspin*this->pw_rho->nrxx); + const int size = PARAM.inp.nspin * this->pw_rho->nrxx; + this->phi_td.resize(size, std::complex(0.0, 0.0)); } - if ((istep<1) && PARAM.inp.init_chg != "file") + if ((istep==0) && PARAM.inp.init_chg != "file") { while (true) { @@ -91,7 +94,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } } } - else if ((istep<1) && PARAM.inp.init_chg == "file") + else if ((istep==0) && PARAM.inp.init_chg == "file") { #ifdef _OPENMP #pragma omp parallel for collapse(2) @@ -107,7 +110,7 @@ void ESolver_OF_TDDFT::runner(UnitCell& ucell, const int istep) } else { - this->evolve_ofdft->propagate_psi(this->pelec, this->chr, ucell, this->phi_td, this->pw_rho); + this->evolve_ofdft->propagate_psi_RK4(this->pelec, this->chr, ucell, this->phi_td, this->pw_rho); #ifdef _OPENMP #pragma omp parallel for collapse(2) #endif diff --git a/source/source_esolver/esolver_of_tddft.h b/source/source_esolver/esolver_of_tddft.h index 99a4238c0f..85293b1761 100644 --- a/source/source_esolver/esolver_of_tddft.h +++ b/source/source_esolver/esolver_of_tddft.h @@ -15,7 +15,7 @@ class ESolver_OF_TDDFT : public ESolver_OF virtual void runner(UnitCell& ucell, const int istep) override; protected: - std::vector> phi_td; // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi(). + std::vector> phi_td; // time dependent wavefunction Evolve_OFDFT* evolve_ofdft=nullptr; }; } // namespace ModuleESolver diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 05a058671d..44c0f8b5c6 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -6,11 +6,11 @@ #include "source_base/parallel_reduce.h" void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, - const Charge& chr, + Charge& chr, UnitCell& ucell, - std::vector> psi_, + std::vector>& psi_, ModulePW::PW_Basis* pw_rho, - std::vector> Hpsi) + std::vector>& Hpsi) { // update rho #ifdef _OPENMP @@ -23,6 +23,7 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, chr.rho[is][ir] = abs(psi_[is * pw_rho->nrxx + ir])*abs(psi_[is * pw_rho->nrxx + ir]); } } + this->renormalize_psi(chr, pw_rho, psi_); pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_eff_v()); // TF potential @@ -45,6 +46,26 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, this->cal_vw_potential_phi(psi_, pw_rho, Hpsi); } +void Evolve_OFDFT::renormalize_psi(Charge& chr, ModulePW::PW_Basis* pw_rho, std::vector>& pphi_) +{ + const double sr = chr.sum_rho(); + const double normalize_factor = PARAM.inp.nelec / sr; + + std::cout<<"sr="<> pphi, pw_rho->recip2real(recipPhi[is], rLapPhi[is]); for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - Hpsi[is * pw_rho->nrxx + ir]+=rLapPhi[is][ir]; + Hpsi[is * pw_rho->nrxx + ir] += rLapPhi[is][ir]; } } @@ -123,7 +144,7 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, delete[] rLapPhi; } -void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, +void Evolve_OFDFT::cal_CD_potential(std::vector>& psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot, double mCD_para) @@ -151,20 +172,20 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, #endif for (int is = 0; is < PARAM.inp.nspin; ++is) { - std::complex *recipCurrent_x=new std::complex[pw_rho->npw]; - std::complex *recipCurrent_y=new std::complex[pw_rho->npw]; - std::complex *recipCurrent_z=new std::complex[pw_rho->npw]; - std::complex *recipCDPotential=new std::complex[pw_rho->npw]; - std::complex *rCurrent_x=new std::complex[pw_rho->nrxx]; - std::complex *rCurrent_y=new std::complex[pw_rho->nrxx]; - std::complex *rCurrent_z=new std::complex[pw_rho->nrxx]; - std::complex *kF_r=new std::complex[pw_rho->nrxx]; - std::complex *rCDPotential=new std::complex[pw_rho->nrxx]; + std::vector> recipCurrent_x(pw_rho->npw); + std::vector> recipCurrent_y(pw_rho->npw); + std::vector> recipCurrent_z(pw_rho->npw); + std::vector> recipCDPotential(pw_rho->npw); + std::vector> rCurrent_x(pw_rho->nrxx); + std::vector> rCurrent_y(pw_rho->nrxx); + std::vector> rCurrent_z(pw_rho->nrxx); + std::vector> kF_r(pw_rho->nrxx); + std::vector> rCDPotential(pw_rho->nrxx); recipPhi[is] = new std::complex[pw_rho->npw]; for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1/3); + kF_r[ir]=std::pow(3*std::pow(ModuleBase::PI*std::abs(rPhi[is][ir]),2),1.0/3.0); } pw_rho->real2recip(rPhi[is], recipPhi[is]); @@ -174,38 +195,40 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, recipCurrent_y[ik]=imag*pw_rho->gcar[ik].y*recipPhi[is][ik]* pw_rho->tpiba; recipCurrent_z[ik]=imag*pw_rho->gcar[ik].z*recipPhi[is][ik]* pw_rho->tpiba; } - pw_rho->recip2real(recipCurrent_x,rCurrent_x); - pw_rho->recip2real(recipCurrent_y,rCurrent_y); - pw_rho->recip2real(recipCurrent_z,rCurrent_z); + pw_rho->recip2real(recipCurrent_x.data(),rCurrent_x.data()); + pw_rho->recip2real(recipCurrent_y.data(),rCurrent_y.data()); + pw_rho->recip2real(recipCurrent_z.data(),rCurrent_z.data()); for (int ir = 0; ir < pw_rho->nrxx; ++ir) { rCurrent_x[ir]=std::imag(rCurrent_x[ir]*std::conj(rPhi[is][ir])); rCurrent_y[ir]=std::imag(rCurrent_y[ir]*std::conj(rPhi[is][ir])); rCurrent_z[ir]=std::imag(rCurrent_z[ir]*std::conj(rPhi[is][ir])); } - pw_rho->real2recip(rCurrent_x,recipCurrent_x); - pw_rho->real2recip(rCurrent_y,recipCurrent_y); - pw_rho->real2recip(rCurrent_z,recipCurrent_z); + pw_rho->real2recip(rCurrent_x.data(),recipCurrent_x.data()); + pw_rho->real2recip(rCurrent_y.data(),recipCurrent_y.data()); + pw_rho->real2recip(rCurrent_z.data(),recipCurrent_z.data()); for (int ik = 0; ik < pw_rho->npw; ++ik) { - recipCDPotential[ik]=recipCurrent_x[ik]*pw_rho->gcar[ik].x - +recipCurrent_y[ik]*pw_rho->gcar[ik].y - +recipCurrent_z[ik]*pw_rho->gcar[ik].z; - recipCDPotential[ik]*=imag/pw_rho->gg[ik]; + recipCDPotential[ik]=recipCurrent_x[ik]*pw_rho->gcar[ik].x+recipCurrent_y[ik]*pw_rho->gcar[ik].y+recipCurrent_z[ik]*pw_rho->gcar[ik].z; + if (pw_rho->gg[ik]==0) + { + recipCDPotential[ik]=0.0; + } + else + { + recipCDPotential[ik]*=imag/sqrt(pw_rho->gg[ik]); + } } - pw_rho->recip2real(recipCDPotential,rCDPotential); + pw_rho->recip2real(recipCDPotential.data(),rCDPotential.data()); for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) - / (2.0*std::pow(std::real(kF_r[ir]),2)); + rpot(0, ir) -= mCD_para*2.0*std::real(rCDPotential[ir])*std::pow(ModuleBase::PI,3) / (2.0*std::pow(std::real(kF_r[ir]),2)); + if (isnan(rpot(0, ir))) + { + rpot(0, ir)=0.0; + } } - delete[] recipCurrent_x; - delete[] recipCurrent_y; - delete[] recipCurrent_z; - delete[] rCurrent_x; - delete[] rCurrent_y; - delete[] rCurrent_z; } #ifdef _OPENMP @@ -220,15 +243,16 @@ void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, delete[] rPhi; } -void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, - const Charge& chr, UnitCell& ucell, - std::vector> pphi_, +void Evolve_OFDFT::propagate_psi_RK4(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, ModulePW::PW_Basis* pw_rho) { - ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK4"); std::complex imag(0.0,1.0); - double dt=PARAM.inp.mdp.md_dt; + double dt=PARAM.inp.mdp.md_dt / ModuleBase::AU_to_FS; const int nspin = PARAM.inp.nspin; const int nrxx = pw_rho->nrxx; const int total_size = nspin * nrxx; @@ -247,7 +271,7 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K1[is * nrxx + ir]=-1.0*K1[is * nrxx + ir]*dt*imag; + K1[is * nrxx + ir]=-0.5*K1[is * nrxx + ir]*dt*imag; // 0.5 convert Ry to Hartree psi1[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K1[is * nrxx + ir]; } } @@ -258,7 +282,7 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K2[is * nrxx + ir]=-1.0*K2[is * nrxx + ir]*dt*imag; + K2[is * nrxx + ir]=-0.5*K2[is * nrxx + ir]*dt*imag; psi2[is * nrxx + ir]=pphi_[is * nrxx + ir]+0.5*K2[is * nrxx + ir]; } } @@ -269,7 +293,7 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K3[is * nrxx + ir]=-1.0*K3[is * nrxx + ir]*dt*imag; + K3[is * nrxx + ir]=-0.5*K3[is * nrxx + ir]*dt*imag; psi3[is * nrxx + ir]=pphi_[is * nrxx + ir]+K3[is * nrxx + ir]; } } @@ -280,12 +304,80 @@ void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, for (int is = 0; is < PARAM.inp.nspin; ++is){ for (int ir = 0; ir < pw_rho->nrxx; ++ir) { - K4[is * nrxx + ir]=-1.0*K4[is * nrxx + ir]*dt*imag; - pphi_[is * nrxx + ir]+=1.0/6.0*(K1[is * nrxx + ir] - +2.0*K2[is * nrxx + ir]+2.0*K3[is * nrxx + ir] - +K4[is * nrxx + ir]); + K4[is * nrxx + ir]=-0.5*K4[is * nrxx + ir]*dt*imag; + pphi_[is * nrxx + ir]+=1.0/6.0*(K1[is * nrxx + ir]+2.0*K2[is * nrxx + ir]+2.0*K3[is * nrxx + ir]+K4[is * nrxx + ir]); + } + } + +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + chr.rho[is][ir] = abs(pphi_[is * pw_rho->nrxx + ir])*abs(pphi_[is * pw_rho->nrxx + ir]); + } + } + this->renormalize_psi(chr, pw_rho, pphi_); + + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK4"); +} + +void Evolve_OFDFT::propagate_psi_RK2(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, + ModulePW::PW_Basis* pw_rho) +{ + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK2"); + + const std::complex imag(0.0, 1.0); + double dt=PARAM.inp.mdp.md_dt / ModuleBase::AU_to_FS; + const int nspin = PARAM.inp.nspin; + const int nrxx = pw_rho->nrxx; + const int total_size = nspin * nrxx; + + std::vector> K1(total_size); + std::vector> K2(total_size); + std::vector> psi_mid(total_size); + + cal_Hpsi(pelec, chr, ucell, pphi_, pw_rho, K1); + +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < nspin; ++is) { + for (int ir = 0; ir < nrxx; ++ir) { + const int idx = is * nrxx + ir; + K1[idx] = -0.5 * K1[idx] * dt * imag; + psi_mid[idx] = pphi_[idx] + 0.5 * K1[idx]; + } + } + + cal_Hpsi(pelec, chr, ucell, psi_mid, pw_rho, K2); + +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < nspin; ++is) { + for (int ir = 0; ir < nrxx; ++ir) { + const int idx = is * nrxx + ir; + K2[idx] = -0.5 * K2[idx] * dt * imag; + pphi_[idx] += K2[idx]; } } - ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagte_psi"); +#ifdef _OPENMP +#pragma omp parallel for collapse(2) +#endif + for (int is = 0; is < nspin; ++is) { + for (int ir = 0; ir < nrxx; ++ir) { + chr.rho[is][ir] = std::norm(pphi_[is * nrxx + ir]); + } + } + + this->renormalize_psi(chr, pw_rho, pphi_); + + ModuleBase::timer::tick("ESolver_OF_TDDFT", "propagate_psi_RK2"); } diff --git a/source/source_pw/module_ofdft/evolve_ofdft.h b/source/source_pw/module_ofdft/evolve_ofdft.h index 45ffffe7ad..8d592a56a4 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.h +++ b/source/source_pw/module_ofdft/evolve_ofdft.h @@ -24,10 +24,19 @@ class Evolve_OFDFT ~Evolve_OFDFT() { } - void propagate_psi(elecstate::ElecState* pelec, - const Charge& chr, UnitCell& ucell, - std::vector> pphi_, + void propagate_psi_RK4(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, ModulePW::PW_Basis* pw_rho); + + void propagate_psi_RK2(elecstate::ElecState* pelec, + Charge& chr, + UnitCell& ucell, + std::vector>& pphi_, + ModulePW::PW_Basis* pw_rho); + + void renormalize_psi(Charge& chr, ModulePW::PW_Basis* pw_rho, std::vector>& pphi_); private: const double c_tf_ @@ -35,18 +44,18 @@ class Evolve_OFDFT * 2; // 10/3*(3*pi^2)^{2/3}, multiply by 2 to convert unit from Hartree to Ry, finally in Ry*Bohr^(-2) void cal_Hpsi(elecstate::ElecState* pelec, - const Charge& chr, + Charge& chr, UnitCell& ucell, - std::vector> psi_, + std::vector>& psi_, ModulePW::PW_Basis* pw_rho, - std::vector> Hpsi); + std::vector>& Hpsi); void cal_tf_potential(const double* const* prho, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot); - void cal_vw_potential_phi(std::vector> pphi, + void cal_vw_potential_phi(std::vector>& pphi, ModulePW::PW_Basis* pw_rho, - std::vector> Hpsi); // -1/2 \nabla^2 \phi - void cal_CD_potential(std::vector> psi_, + std::vector>& Hpsi); // -1/2 \nabla^2 \phi + void cal_CD_potential(std::vector>& psi_, ModulePW::PW_Basis* pw_rho, ModuleBase::matrix& rpot, double mCD_para); diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/result.ref b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref index 404887187c..b53afa1a10 100644 --- a/tests/07_OFDFT/30_TDOFDFT_Al/result.ref +++ b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref @@ -1,5 +1,5 @@ -etotref -100.6027239728818614 -etotperatomref -50.3013619864 -totalforceref 5.786358 -totalstressref 11367.508779 -totaltimeref 0.10 +etotref 226.9979568201541724 +etotperatomref 113.4989784101 +totalforceref 562.524006 +totalstressref 194931.772028 +totaltimeref 0.45 diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/threshold b/tests/07_OFDFT/30_TDOFDFT_Al/threshold new file mode 100644 index 0000000000..a137b5cae7 --- /dev/null +++ b/tests/07_OFDFT/30_TDOFDFT_Al/threshold @@ -0,0 +1,4 @@ +threshold 0.0001 +force_threshold 0.001 +stress_threshold 0.01 +fatal_threshold 1 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref index b94d694de2..286a081395 100644 --- a/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref @@ -1,5 +1,5 @@ -etotref -100.6027239728818614 -etotperatomref -50.3013619864 -totalforceref 5.786358 -totalstressref 11367.508779 -totaltimeref 0.13 +etotref 225.6175831066692297 +etotperatomref 112.8087915533 +totalforceref 488.530072 +totalstressref 185999.479726 +totaltimeref 0.09 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold b/tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold new file mode 100644 index 0000000000..a137b5cae7 --- /dev/null +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/threshold @@ -0,0 +1,4 @@ +threshold 0.0001 +force_threshold 0.001 +stress_threshold 0.01 +fatal_threshold 1 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref index 3b61275a93..a6e3cab7f3 100644 --- a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref @@ -1,5 +1,5 @@ -etotref -100.6027239728818614 -etotperatomref -50.3013619864 -totalforceref 5.786358 -totalstressref 11367.508779 -totaltimeref 0.11 +etotref 223.1193797466232809 +etotperatomref 111.5596898733 +totalforceref 386.437466 +totalstressref 181323.311310 +totaltimeref 0.08 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold new file mode 100644 index 0000000000..a137b5cae7 --- /dev/null +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/threshold @@ -0,0 +1,4 @@ +threshold 0.0001 +force_threshold 0.001 +stress_threshold 0.01 +fatal_threshold 1