[Feature] Add current dependent potential for TD-OFDFT

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

source/source_io/module_parameter/input_parameter.h

@@ -687,5 +687,10 @@ struct Input_para
     // src/gga_c_pbe.c
     std::vector<double> 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

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();

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();

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()
 {

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");

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

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<std::complex<double>> pphi,
                                         ModulePW::PW_Basis* pw_rho, 
                                         std::vector<std::complex<double>> Hpsi)
 {
+    if (PARAM.inp.nspin <= 0) {
+        ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive");
+    }
     std::complex<double>** rLapPhi = new std::complex<double>*[PARAM.inp.nspin];
 #ifdef _OPENMP
 #pragma omp parallel for
@@ -118,13 +125,96 @@ void Evolve_OFDFT::cal_vw_potential_phi(std::vector<std::complex<double>> pphi,
 
 void Evolve_OFDFT::cal_CD_potential(std::vector<std::complex<double>> psi_, 
                                     ModulePW::PW_Basis* pw_rho, 
-                                    ModuleBase::matrix& rpot)
+                                    ModuleBase::matrix& rpot,
+                                    double mCD_para)
 {
+    std::complex<double> imag(0.0,1.0);
+
+    if (PARAM.inp.nspin <= 0) {
+        ModuleBase::WARNING_QUIT("Evolve_OFDFT","nspin must be positive");
+    }
+    std::complex<double>** recipPhi = new std::complex<double>*[PARAM.inp.nspin];
+    std::complex<double>** rPhi = new std::complex<double>*[PARAM.inp.nspin];
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+    for (int is = 0; is < PARAM.inp.nspin; ++is) {
+        rPhi[is] = new std::complex<double>[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<double>[pw_rho->npw];
-        //delete[] recipCurrent;
+        std::complex<double> *recipCurrent_x=new std::complex<double>[pw_rho->npw];
+        std::complex<double> *recipCurrent_y=new std::complex<double>[pw_rho->npw];
+        std::complex<double> *recipCurrent_z=new std::complex<double>[pw_rho->npw];
+        std::complex<double> *recipCDPotential=new std::complex<double>[pw_rho->npw];
+        std::complex<double> *rCurrent_x=new std::complex<double>[pw_rho->nrxx];
+        std::complex<double> *rCurrent_y=new std::complex<double>[pw_rho->nrxx];
+        std::complex<double> *rCurrent_z=new std::complex<double>[pw_rho->nrxx];
+        std::complex<double> *kF_r=new std::complex<double>[pw_rho->nrxx];
+        std::complex<double> *rCDPotential=new std::complex<double>[pw_rho->nrxx];
+        recipPhi[is] = new std::complex<double>[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, 

source/source_pw/module_ofdft/evolve_ofdft.h

@@ -48,7 +48,8 @@ class Evolve_OFDFT
                               std::vector<std::complex<double>> Hpsi); // -1/2 \nabla^2 \phi
     void cal_CD_potential(std::vector<std::complex<double>> psi_, 
                           ModulePW::PW_Basis* pw_rho, 
-                          ModuleBase::matrix& rpot);
+                          ModuleBase::matrix& rpot,
+                          double mCD_para);
 
 };
 #endif

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

tests/07_OFDFT/30_TDOFDFT_Al/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0

tests/07_OFDFT/30_TDOFDFT_Al/README

@@ -0,0 +1 @@
+Test case without current dependent potential in TDOFDFT

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

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

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

tests/07_OFDFT/31_TDOFDFT_Al_CD/KPT

@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0

tests/07_OFDFT/31_TDOFDFT_Al_CD/README

@@ -0,0 +1 @@
+Test case with current dependent potential in TDOFDFT

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

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