diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 001bb4e5e3..52ba9f8ee7 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -256,6 +256,9 @@ - [of\_ml\_chi\_pnl](#of_ml_chi_pnl) - [of\_ml\_chi\_qnl](#of_ml_chi_qnl) - [of\_ml\_local\_test](#of_ml_local_test) + - [TD-OFDFT: time dependent orbital free density functional theory](#tdofdft-time-dependent-orbital-free-density-functional-theory) + - [of\_cd](#of_cd) + - [of\_mcd\_alpha](#of_mcd_alpha) - [Electric Field and Dipole Correction](#electric-field-and-dipole-correction) - [efield\_flag](#efield_flag) - [dip\_cor\_flag](#dip_cor_flag) @@ -528,6 +531,7 @@ These variables are used to control general system parameters. - **Description**: choose the energy solver. - ksdft: Kohn-Sham density functional theory - ofdft: orbital-free density functional theory + - tdofdft: time-dependent orbital-free density functional theory - sdft: [stochastic density functional theory](#electronic-structure-sdft) - tddft: real-time time-dependent density functional theory (TDDFT) - lj: Leonard Jones potential @@ -2747,6 +2751,25 @@ Warning: this function is not robust enough for the current version. Please try [back to top](#full-list-of-input-keywords) +## TDOFDFT: time dependent orbital free density functional theory + +### of_cd + +- **Type**: Boolean +- **Availability**: TDOFDFT +- **Type**: Boolean +- **Description**: Added the current dependent(CD) potential. (https://doi.org/10.1103/PhysRevB.98.144302) + - True: Added the CD potential. + - False: Not added the CD potential. +- **Default**: False + +### of_mcd_alpha + +- **Type**: Real +- **Availability**: TDOFDFT +- **Description**: The value of the parameter alpha in modified CD potential method. mCDPotenial=alpha*CDPotenial(proposed in paper PhysRevB.98.144302) +- **Default**: 1.0 + ## Electric field and dipole correction These variables are relevant to electric field and dipole correction diff --git a/source/source_io/module_parameter/input_parameter.h b/source/source_io/module_parameter/input_parameter.h index a8dd6d63bb..b1b2c08fce 100644 --- a/source/source_io/module_parameter/input_parameter.h +++ b/source/source_io/module_parameter/input_parameter.h @@ -687,5 +687,10 @@ struct Input_para // src/gga_c_pbe.c std::vector xc_corr_ext = { 130, 0.06672455060314922, 0.031090690869654895034, 1.00000}; + + // ============== #Parameters (24.td-ofdft) =========================== + bool of_cd = false; ///< add CD potential or not https://doi.org/10.1103/PhysRevB.98.144302 + double of_mCD_alpha = 1.0; /// parameter of modified CD Potential + }; #endif diff --git a/source/source_io/read_input.cpp b/source/source_io/read_input.cpp index c59c31d8c6..14f25d115d 100644 --- a/source/source_io/read_input.cpp +++ b/source/source_io/read_input.cpp @@ -101,6 +101,7 @@ ReadInput::ReadInput(const int& rank) this->item_sdft(); this->item_deepks(); this->item_rt_tddft(); + this->item_tdofdft(); this->item_lr_tddft(); this->item_output(); this->item_postprocess(); diff --git a/source/source_io/read_input.h b/source/source_io/read_input.h index 8f1f2e0ec0..92bda04bdb 100644 --- a/source/source_io/read_input.h +++ b/source/source_io/read_input.h @@ -108,6 +108,8 @@ class ReadInput // items for real time tddft void item_rt_tddft(); // items for linear response tddft + void item_tdofdft(); + // items for td-ofdft void item_lr_tddft(); // items for output void item_output(); diff --git a/source/source_io/read_input_item_tddft.cpp b/source/source_io/read_input_item_tddft.cpp index 89751fd224..66aa47164e 100644 --- a/source/source_io/read_input_item_tddft.cpp +++ b/source/source_io/read_input_item_tddft.cpp @@ -305,6 +305,22 @@ void ReadInput::item_rt_tddft() } +} +void ReadInput::item_tdofdft() +{ + // TD-OFDFT + { + Input_Item item("of_cd"); + item.annotation = "add CD Potential or not"; + read_sync_bool(input.of_cd); + this->add_item(item); + } + { + Input_Item item("of_mcd_alpha"); + item.annotation = "parameter of modified CD Potential"; + read_sync_double(input.of_mCD_alpha); + this->add_item(item); + } } void ReadInput::item_lr_tddft() { diff --git a/source/source_io/test/read_input_ptest.cpp b/source/source_io/test/read_input_ptest.cpp index a8cdfde867..a25bc2a926 100644 --- a/source/source_io/test/read_input_ptest.cpp +++ b/source/source_io/test/read_input_ptest.cpp @@ -364,6 +364,8 @@ TEST_F(InputParaTest, ParaRead) EXPECT_EQ(param.inp.of_full_pw_dim, 0); EXPECT_FALSE(param.inp.of_read_kernel); EXPECT_EQ(param.inp.of_kernel_file, "WTkernel.txt"); + EXPECT_FALSE(param.inp.of_cd); + EXPECT_DOUBLE_EQ(param.inp.of_mCD_alpha,1.0); EXPECT_DOUBLE_EQ(param.inp.of_xwm_kappa, 1.); EXPECT_DOUBLE_EQ(param.inp.of_xwm_rho_ref, 1.); EXPECT_EQ(param.inp.device, "cpu"); diff --git a/source/source_io/test/support/INPUT b/source/source_io/test/support/INPUT index d69b67a6b3..46840bd9da 100644 --- a/source/source_io/test/support/INPUT +++ b/source/source_io/test/support/INPUT @@ -387,3 +387,7 @@ nsc_min 4 #Minimum number of spin-constrained iteration sc_scf_nmin 4 #Minimum number of outer scf loop before initializing lambda loop alpha_trial 0.02 #Initial trial step size for lambda in eV/uB^2 sccut 4 #Maximal step size for lambda in eV/uB + +#Parameters (23. Time-dependent orbital-free DFT) +of_cd 0 #0: no CD potential; 1: add CD potential +of_mCD_alpha 1.0 # parameter of modified CD potential \ No newline at end of file diff --git a/source/source_pw/module_ofdft/evolve_ofdft.cpp b/source/source_pw/module_ofdft/evolve_ofdft.cpp index 4c158998cc..f4aaea2cec 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.cpp +++ b/source/source_pw/module_ofdft/evolve_ofdft.cpp @@ -25,7 +25,11 @@ void Evolve_OFDFT::cal_Hpsi(elecstate::ElecState* pelec, } pelec->pot->update_from_charge(&chr, &ucell); // Hartree + XC + external - this->cal_tf_potential(chr.rho,pw_rho ,pelec->pot->get_effective_v()); // TF potential + this->cal_tf_potential(chr.rho, pw_rho, pelec->pot->get_effective_v()); // TF potential + if (PARAM.inp.of_cd) + { + this->cal_CD_potential(psi_, pw_rho, pelec->pot->get_effective_v(), PARAM.inp.of_mCD_alpha); // CD potential + } #ifdef _OPENMP #pragma omp parallel for @@ -72,6 +76,9 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, ModulePW::PW_Basis* pw_rho, std::vector> Hpsi) { + if (PARAM.inp.nspin <= 0) { + ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive"); + } std::complex** rLapPhi = new std::complex*[PARAM.inp.nspin]; #ifdef _OPENMP #pragma omp parallel for @@ -118,13 +125,96 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector> pphi, void Evolve_OFDFT::cal_CD_potential(std::vector> psi_, ModulePW::PW_Basis* pw_rho, - ModuleBase::matrix& rpot) + ModuleBase::matrix& rpot, + double mCD_para) { + std::complex imag(0.0,1.0); + + if (PARAM.inp.nspin <= 0) { + ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive"); + } + std::complex** recipPhi = new std::complex*[PARAM.inp.nspin]; + std::complex** rPhi = new std::complex*[PARAM.inp.nspin]; +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int is = 0; is < PARAM.inp.nspin; ++is) { + rPhi[is] = new std::complex[pw_rho->nrxx]; + for (int ir = 0; ir < pw_rho->nrxx; ++ir) + { + rPhi[is][ir]=psi_[is * pw_rho->nrxx + ir]; + } + } + +#ifdef _OPENMP +#pragma omp parallel for +#endif for (int is = 0; is < PARAM.inp.nspin; ++is) { - //recipCurrent = new std::complex[pw_rho->npw]; - //delete[] recipCurrent; + 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]; + 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); + } + + pw_rho->real2recip(rPhi[is], recipPhi[is]); + for (int ik = 0; ik < pw_rho->npw; ++ik) + { + recipCurrent_x[ik]=imag*pw_rho->gcar[ik].x*recipPhi[is][ik]* pw_rho->tpiba; + 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); + 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); + 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]; + } + pw_rho->recip2real(recipCDPotential,rCDPotential); + + 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)); + } + delete[] recipCurrent_x; + delete[] recipCurrent_y; + delete[] recipCurrent_z; + delete[] rCurrent_x; + delete[] rCurrent_y; + delete[] rCurrent_z; } + +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int is = 0; is < PARAM.inp.nspin; ++is) + { + delete[] recipPhi[is]; + delete[] rPhi[is]; + } + delete[] recipPhi; + delete[] rPhi; } void Evolve_OFDFT::propagate_psi(elecstate::ElecState* pelec, diff --git a/source/source_pw/module_ofdft/evolve_ofdft.h b/source/source_pw/module_ofdft/evolve_ofdft.h index d3cd6d358b..45ffffe7ad 100644 --- a/source/source_pw/module_ofdft/evolve_ofdft.h +++ b/source/source_pw/module_ofdft/evolve_ofdft.h @@ -48,7 +48,8 @@ class Evolve_OFDFT std::vector> Hpsi); // -1/2 \nabla^2 \phi void cal_CD_potential(std::vector> psi_, ModulePW::PW_Basis* pw_rho, - ModuleBase::matrix& rpot); + ModuleBase::matrix& rpot, + double mCD_para); }; #endif \ No newline at end of file diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/INPUT b/tests/07_OFDFT/30_TDOFDFT_Al/INPUT new file mode 100644 index 0000000000..17fa540e80 --- /dev/null +++ b/tests/07_OFDFT/30_TDOFDFT_Al/INPUT @@ -0,0 +1,36 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation md +esolver_type tdofdft + +pseudo_dir ../../PP_ORB +pseudo_rcut 16 +cal_force 1 +cal_stress 1 + +#Parameters (2.Iteration) +ecutwfc 17 +scf_nmax 100 + +#OFDFT +of_kinetic tf+ +of_method tn +of_full_pw 1 +of_full_pw_dim 1 +of_cd 0 + +#Parameters (3.Basis) +basis_type pw + +init_vel 1 + +md_restart 0 +md_type nve +md_nstep 2 +md_dt 0.25 +md_tfirst 58022.52706 +md_dumpfreq 10 +md_tfreq 1.08 +md_tchain 1 +nbspline 10 # be sure fft dimension is odd diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/KPT b/tests/07_OFDFT/30_TDOFDFT_Al/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/07_OFDFT/30_TDOFDFT_Al/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/README b/tests/07_OFDFT/30_TDOFDFT_Al/README new file mode 100644 index 0000000000..be31e1c21f --- /dev/null +++ b/tests/07_OFDFT/30_TDOFDFT_Al/README @@ -0,0 +1 @@ +Test case without current dependent potential in TDOFDFT diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/STRU b/tests/07_OFDFT/30_TDOFDFT_Al/STRU new file mode 100644 index 0000000000..1498bc2cc1 --- /dev/null +++ b/tests/07_OFDFT/30_TDOFDFT_Al/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Al 26.98 al.lda.lps blps + +LATTICE_CONSTANT +7.50241114482312 // add lattice constant + +LATTICE_VECTORS +0.000000000000 0.500000000000 0.500000000000 +0.500000000000 0.000000000000 0.500000000000 +0.500000000000 0.500000000000 0.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +2 +0.000000000000 0.000000000000 0.000000000000 1 1 1 v 0 0 0.5 +0.500000000000 0.500000000000 0.500000000000 1 1 1 v 0 0 -0.5 diff --git a/tests/07_OFDFT/30_TDOFDFT_Al/result.ref b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref new file mode 100644 index 0000000000..404887187c --- /dev/null +++ b/tests/07_OFDFT/30_TDOFDFT_Al/result.ref @@ -0,0 +1,5 @@ +etotref -100.6027239728818614 +etotperatomref -50.3013619864 +totalforceref 5.786358 +totalstressref 11367.508779 +totaltimeref 0.10 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT b/tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT new file mode 100644 index 0000000000..b14e49962e --- /dev/null +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/INPUT @@ -0,0 +1,37 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation md +esolver_type tdofdft + +pseudo_dir ../../PP_ORB +pseudo_rcut 16 +cal_force 1 +cal_stress 1 + +#Parameters (2.Iteration) +ecutwfc 17 +scf_nmax 100 + +#OFDFT +of_kinetic tf+ +of_method tn +of_full_pw 1 +of_full_pw_dim 1 +of_cd 1 +of_mcd_alpha 1.0 + +#Parameters (3.Basis) +basis_type pw + +init_vel 1 + +md_restart 0 +md_type nve +md_nstep 2 +md_dt 0.25 +md_tfirst 58022.52706 +md_dumpfreq 10 +md_tfreq 1.08 +md_tchain 1 +nbspline 10 # be sure fft dimension is odd diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/KPT b/tests/07_OFDFT/31_TDOFDFT_Al_CD/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/README b/tests/07_OFDFT/31_TDOFDFT_Al_CD/README new file mode 100644 index 0000000000..d09f6a4ab9 --- /dev/null +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/README @@ -0,0 +1 @@ +Test case with current dependent potential in TDOFDFT diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/STRU b/tests/07_OFDFT/31_TDOFDFT_Al_CD/STRU new file mode 100644 index 0000000000..1498bc2cc1 --- /dev/null +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Al 26.98 al.lda.lps blps + +LATTICE_CONSTANT +7.50241114482312 // add lattice constant + +LATTICE_VECTORS +0.000000000000 0.500000000000 0.500000000000 +0.500000000000 0.000000000000 0.500000000000 +0.500000000000 0.500000000000 0.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +2 +0.000000000000 0.000000000000 0.000000000000 1 1 1 v 0 0 0.5 +0.500000000000 0.500000000000 0.500000000000 1 1 1 v 0 0 -0.5 diff --git a/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref new file mode 100644 index 0000000000..b94d694de2 --- /dev/null +++ b/tests/07_OFDFT/31_TDOFDFT_Al_CD/result.ref @@ -0,0 +1,5 @@ +etotref -100.6027239728818614 +etotperatomref -50.3013619864 +totalforceref 5.786358 +totalstressref 11367.508779 +totaltimeref 0.13 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT new file mode 100644 index 0000000000..b4c1593e9f --- /dev/null +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/INPUT @@ -0,0 +1,37 @@ +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest +calculation md +esolver_type tdofdft + +pseudo_dir ../../PP_ORB +pseudo_rcut 16 +cal_force 1 +cal_stress 1 + +#Parameters (2.Iteration) +ecutwfc 17 +scf_nmax 100 + +#OFDFT +of_kinetic tf+ +of_method tn +of_full_pw 1 +of_full_pw_dim 1 +of_cd 1 +of_mcd_alpha 2.0 + +#Parameters (3.Basis) +basis_type pw + +init_vel 1 + +md_restart 0 +md_type nve +md_nstep 2 +md_dt 0.25 +md_tfirst 58022.52706 +md_dumpfreq 10 +md_tfreq 1.08 +md_tchain 1 +nbspline 10 # be sure fft dimension is odd diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/KPT b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/KPT new file mode 100644 index 0000000000..c289c0158a --- /dev/null +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/README b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/README new file mode 100644 index 0000000000..cedbc26afd --- /dev/null +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/README @@ -0,0 +1 @@ +Test case with modified current dependent potential (alpha=2.0) in TDOFDFT diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/STRU b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/STRU new file mode 100644 index 0000000000..1498bc2cc1 --- /dev/null +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/STRU @@ -0,0 +1,19 @@ +ATOMIC_SPECIES +Al 26.98 al.lda.lps blps + +LATTICE_CONSTANT +7.50241114482312 // add lattice constant + +LATTICE_VECTORS +0.000000000000 0.500000000000 0.500000000000 +0.500000000000 0.000000000000 0.500000000000 +0.500000000000 0.500000000000 0.000000000000 + +ATOMIC_POSITIONS +Direct + +Al +0 +2 +0.000000000000 0.000000000000 0.000000000000 1 1 1 v 0 0 0.5 +0.500000000000 0.500000000000 0.500000000000 1 1 1 v 0 0 -0.5 diff --git a/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref new file mode 100644 index 0000000000..3b61275a93 --- /dev/null +++ b/tests/07_OFDFT/32_TDOFDFT_Al_mCD/result.ref @@ -0,0 +1,5 @@ +etotref -100.6027239728818614 +etotperatomref -50.3013619864 +totalforceref 5.786358 +totalstressref 11367.508779 +totaltimeref 0.11 diff --git a/tests/07_OFDFT/CASES_CPU.txt b/tests/07_OFDFT/CASES_CPU.txt index 67e095d710..e94290c382 100644 --- a/tests/07_OFDFT/CASES_CPU.txt +++ b/tests/07_OFDFT/CASES_CPU.txt @@ -26,4 +26,7 @@ 26_OF_MD_LibxcPBE 27_OF_CR_LDA 28_OF_KE_XWM -29_OF_XWM_para \ No newline at end of file +29_OF_XWM_para +30_TDOFDFT_Al +31_TDOFDFT_Al_CD +32_TDOFDFT_Al_mCD