Merge branch 'eband' of https://github.com/sunliang98/abacus-develop into eband

Author: sunliang98

docs/advanced/input_files/input-main.md

@@ -255,6 +255,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)
@@ -527,6 +530,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
@@ -1072,7 +1076,20 @@ calculations.
 
 - **Type**: String
 - **Description**: In our package, the XC functional can either be set explicitly using the `dft_functional` keyword in `INPUT` file. If `dft_functional` is not specified, ABACUS will use the xc functional indicated in the pseudopotential file.
-  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names such as 'LDA', 'PBE', 'SCAN'. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/source_hamilt/module_xc/xc_functional.cpp). The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://libxc.gitlab.io/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**.
+  On the other hand, if dft_functional is specified, it will overwrite the functional from pseudopotentials and performs calculation with whichever functional the user prefers. We further offer two ways of supplying exchange-correlation functional. The first is using 'short-hand' names. A complete list of 'short-hand' expressions can be found in [the source code](../../../source/source_hamilt/module_xc/xc_functional.cpp). Supported density functionals are:
+  - LDA functionals
+    - LDA (equivalent with PZ and SLAPZNOGXNOGC), PWLDA
+  - GGA functionals
+    - PBE (equivalent with SLAPWPBXPBC), PBESOL, REVPBE, WC, BLYP, BP(referred to BP86), PW91, HCTH, OLYP, BLYP_LR
+  - meta-GGA functionals
+    - SCAN (require LIBXC)
+  - Hybrid functionals
+    - PBE0, HF
+    - If LIBXC is avaliale, additional short-hand names of hybrid functionals are supported: HSE(referred to HSE06), B3LYP, LC_PBE, LC_WPBE, LRC_WPBE, LRC_WPBEH, CAM_PBEH, WP22, CWP22, MULLER (equivalent with POWER)
+  - Hybrid meta-GGA functionals
+    - SCAN0 (require LIBXC)
+
+   The other way is only available when ***compiling with LIBXC***, and it allows for supplying exchange-correlation functionals as combinations of LIBXC keywords for functional components, joined by a plus sign, for example, dft_functional='LDA_X_1D_EXPONENTIAL+LDA_C_1D_CSC'. The list of LIBXC keywords can be found on its [website](https://libxc.gitlab.io/functionals/). In this way, **we support all the LDA,GGA and mGGA functionals provided by LIBXC**. Some popular functionals and their usage are: RPBE of [Hammer et al.](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.59.7413), set `dft_functional` to 'GGA_X_RPBE+GGA_C_PBE', and [r$^{2}$SCAN](https://pubs.acs.org/doi/10.1021/acs.jpclett.0c02405), set `dft_functional` to 'MGGA_X_R2SCAN+MGGA_C_R2SCAN'.
 
   Furthermore, the old INPUT parameter exx_hybrid_type for hybrid functionals has been absorbed into dft_functional. Options are `hf` (pure Hartree-Fock), `pbe0`(PBE0), `hse` (Note: in order to use HSE functional, LIBXC is required). Note also that HSE has been tested while PBE0 has NOT been fully tested yet, and the maximum CPU cores for running exx in parallel is $N(N+1)/2$, with N being the number of atoms.
 
@@ -2725,6 +2742,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_cell/update_cell.cpp

@@ -496,12 +496,13 @@ void periodic_boundary_adjustment(Atom* atoms,
                     atom->taud[ia][ik] -= 1.0;
                 }
             }
-            if (atom->taud[ia].x < 0 
-                || atom->taud[ia].y < 0
-                || atom->taud[ia].z < 0 
-                || atom->taud[ia].x >= 1.0
-                || atom->taud[ia].y >= 1.0 
-                || atom->taud[ia].z >= 1.0) 
+            const double eps = 1e-12;
+            if (atom->taud[ia].x < -eps
+                || atom->taud[ia].y < -eps
+                || atom->taud[ia].z < -eps
+                || atom->taud[ia].x >= 1.0+eps
+                || atom->taud[ia].y >= 1.0+eps
+                || atom->taud[ia].z >= 1.0+eps) 
             {
                 GlobalV::ofs_warning << " atom type=" << it + 1 << " atom index=" << ia + 1 << std::endl;
                 GlobalV::ofs_warning << " direct coordinate=" << atom->taud[ia].x << " "

source/source_io/module_parameter/input_parameter.h

@@ -686,5 +686,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