Merge develop branch to TDDFT branch

docs/advanced/input_files/input-main.md

@@ -312,9 +312,6 @@
 		- [out\_efield](#out_efield)
 		- [ocp](#ocp)
 		- [ocp\_set](#ocp_set)
-		- [td\_val\_elec\_01](#td_val_elec_01)
-		- [td\_val\_elec\_02](#td_val_elec_02)
-		- [td\_val\_elec\_03](#td_val_elec_03)
 	- [Variables useful for debugging](#variables-useful-for-debugging)
 		- [nurse](#nurse)
 		- [t\_in\_h](#t_in_h)
@@ -2588,24 +2585,6 @@ These variables are used to control berry phase and wannier90 interface paramete
 - **Description**: If ocp is true, the ocp_set is a string to set the number of occupancy, like 1 10 * 1 0 1 representing the 13 band occupancy, 12th band occupancy 0 and the rest 1, the code is parsing this string into an array through a regular expression.
 - **Default**: none
 
-### td_val_elec_01
-
-- **Type**: Integer
-- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the first element.
-- **Default**: 1.0
-
-### td_val_elec_02
-
-- **Type**: Integer
-- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the second element.
-- **Default**: 1.0
-
-### td_val_elec_03
-
-- **Type**: Integer
-- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the third element.
-- **Default**: 1.0
-
 [back to top](#full-list-of-input-keywords)
 
 ## Variables useful for debugging

source/module_elecstate/potentials/H_TDDFT_pw.cpp

@@ -1,9 +1,10 @@
 #include "H_TDDFT_pw.h"
 
-#include "module_io/input.h"
 #include "module_base/constants.h"
 #include "module_base/timer.h"
 #include "module_hamilt_lcao/module_tddft/ELEC_evolve.h"
+#include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_io/input.h"
 #include "module_io/input_conv.h"
 
 namespace elecstate
@@ -16,7 +17,7 @@ int H_TDDFT_pw::istep = -1;
 // fuxiang add in 2017-05
 //==========================================================
 
-void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo)
+void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
 {
     ModuleBase::TITLE("H_TDDFT_pw", "cal_fixed_v");
 
@@ -39,18 +40,19 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo)
     for (auto direc: ELEC_evolve::td_vext_dire_case)
     {
         std::vector<double> vext_space(this->rho_basis_->nrxx, 0.0);
-        double vext_time = cal_v_time(ttype[count]); 
+        double vext_time = cal_v_time(ttype[count]);
 
         if (ELEC_evolve::out_efield && GlobalV::MY_RANK == 0)
         {
             std::stringstream as;
-            as << GlobalV::global_out_dir << "efield_"<<count<<".dat";
+            as << GlobalV::global_out_dir << "efield_" << count << ".dat";
             std::ofstream ofs(as.str().c_str(), std::ofstream::app);
-            ofs << H_TDDFT_pw::istep*dt*ModuleBase::AU_to_FS << "\t" << vext_time <<endl;
+            ofs << H_TDDFT_pw::istep * dt * ModuleBase::AU_to_FS << "\t"
+                << vext_time * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << endl;
             ofs.close();
         }
 
-        cal_v_space(vext_space, direc); 
+        cal_v_space(vext_space, direc);
         for (size_t ir = 0; ir < this->rho_basis_->nrxx; ++ir)
             vl_pseudo[ir] += vext_space[ir] * vext_time;
         count++;
@@ -94,7 +96,7 @@ void H_TDDFT_pw::read_parameters(Input *in)
     gauss_count = 0;
     gauss_omega = set_parameters(in->td_gauss_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1
     gauss_phase = set_parameters(in->td_gauss_phase, 1.0);
-    gauss_sigma = set_parameters(in->td_gauss_sigma, 1/ModuleBase::AU_to_FS);
+    gauss_sigma = set_parameters(in->td_gauss_sigma, 1 / ModuleBase::AU_to_FS);
     gauss_t0 = set_parameters(in->td_gauss_t0, 1.0);
     gauss_amp = set_parameters(in->td_gauss_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr
 
@@ -163,6 +165,12 @@ void H_TDDFT_pw::cal_v_space_length(std::vector<double> &vext_space, int direc)
     ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space_length");
     ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space_length");
 
+    double bmod[3];
+    for (int i = 0; i < 3; i++)
+    {
+        bmod[i] = prepare(GlobalC::ucell, i);
+    }
+
     for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir)
     {
         int i = ir / (this->rho_basis_->ny * this->rho_basis_->nplane);
@@ -175,15 +183,15 @@ void H_TDDFT_pw::cal_v_space_length(std::vector<double> &vext_space, int direc)
         switch (direc)
         {
         case 1:
-            vext_space[ir] = cal_v_space_length_potential(x);
+            vext_space[ir] = cal_v_space_length_potential(x) / bmod[0];
             break;
 
         case 2:
-            vext_space[ir] = cal_v_space_length_potential(y);
+            vext_space[ir] = cal_v_space_length_potential(y) / bmod[1];
             break;
 
         case 3:
-            vext_space[ir] = cal_v_space_length_potential(z);
+            vext_space[ir] = cal_v_space_length_potential(z) / bmod[2];
             break;
 
         default:
@@ -198,18 +206,18 @@ void H_TDDFT_pw::cal_v_space_length(std::vector<double> &vext_space, int direc)
 
 double H_TDDFT_pw::cal_v_space_length_potential(double i)
 {
-    double vext_space=0.0;
-    if (i < this->rho_basis_->nx * lcut1)
+    double vext_space = 0.0;
+    if (i < lcut1)
     {
-        vext_space = ((i / this->rho_basis_->nx - lcut1)*(lcut2-lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0;
+        vext_space = ((i - lcut1) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0;
     }
-    else if (i >= this->rho_basis_->nx * lcut1 && i < this->rho_basis_->nx * lcut2)
+    else if (i >= lcut1 && i < lcut2)
     {
-        vext_space = -i / this->rho_basis_->nx * this->ucell_->lat0;
+        vext_space = -i * this->ucell_->lat0;
     }
-    else if (i >= this->rho_basis_->nx * lcut2)
+    else if (i >= lcut2)
     {
-        vext_space = ((i / this->rho_basis_->nx - lcut2)*(lcut2-lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0;
+        vext_space = ((i - lcut2) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0;
     }
     return vext_space;
 }
@@ -229,6 +237,10 @@ double H_TDDFT_pw::cal_v_time(int t_type)
         vext_time = cal_v_time_Gauss();
         break;
 
+    case 1:
+        vext_time = cal_v_time_trapezoid();
+        break;
+
     case 2:
         vext_time = cal_v_time_trigonometric();
         break;
@@ -237,9 +249,9 @@ double H_TDDFT_pw::cal_v_time(int t_type)
         vext_time = cal_v_time_heaviside();
         break;
 
-    // case 4:
-    //     vext_time = cal_v_time_HHG();
-    //     break;
+        // case 4:
+        //     vext_time = cal_v_time_HHG();
+        //     break;
 
     default:
         std::cout << "time_domain_type of electric field is wrong" << endl;
@@ -304,8 +316,7 @@ double H_TDDFT_pw::cal_v_time_trigonometric()
 
     const double timenow = istep * dt;
 
-    vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2)
-                * sin(omega2 * timenow + phase2);
+    vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) * sin(omega2 * timenow + phase2);
     trigo_count++;
 
     return vext_time;
@@ -345,4 +356,34 @@ double H_TDDFT_pw::cal_v_time_heaviside()
 //     return vext_time;
 // }
 
+double H_TDDFT_pw::prepare(const UnitCell &cell, int &dir)
+{
+    double bvec[3] = {0.0};
+    double bmod = 0.0;
+    if (dir == 0)
+    {
+        bvec[0] = cell.G.e11;
+        bvec[1] = cell.G.e12;
+        bvec[2] = cell.G.e13;
+    }
+    else if (dir == 1)
+    {
+        bvec[0] = cell.G.e21;
+        bvec[1] = cell.G.e22;
+        bvec[2] = cell.G.e23;
+    }
+    else if (dir == 2)
+    {
+        bvec[0] = cell.G.e31;
+        bvec[1] = cell.G.e32;
+        bvec[2] = cell.G.e33;
+    }
+    else
+    {
+        ModuleBase::WARNING_QUIT("H_TDDFT_pw::prepare", "direction is wrong!");
+    }
+    bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2));
+    return bmod;
+}
+
 } // namespace elecstate

source/module_elecstate/potentials/H_TDDFT_pw.h

@@ -109,6 +109,8 @@ class H_TDDFT_pw : public PotBase
     double cal_v_time_trigonometric();
     double cal_v_time_heaviside();
     // double cal_v_time_HHG();
+
+    double prepare(const UnitCell &cell, int &dir);
 };
 
 } // namespace elecstate

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -365,12 +365,15 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter)
                 = new psi::Psi<std::complex<double>>(GlobalC::kv.nks, GlobalV::NBANDS, GlobalV::NLOCAL, nullptr);
 #endif
 
-        std::complex<double> *p_psi = &psi[0](0,0,0);
-        std::complex<double> *p_psi_laststep = &psi_laststep[0](0,0,0);
-        for (int index = 0; index < psi[0].size(); ++index)
+        for (int ik = 0; ik < GlobalC::kv.nks; ++ik)
         {
-            p_psi_laststep[index] = p_psi[index];
+            this->psi->fix_k(ik);
+            this->psi_laststep->fix_k(ik);
+            int size0 = psi->get_nbands() * psi->get_nbasis();
+            for (int index = 0; index < size0; ++index)
+                psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index];
         }
+
         if (istep > 1 && ELEC_evolve::td_edm == 0)
             this->cal_edm_tddft();
     }

source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp

@@ -12,9 +12,6 @@ ELEC_evolve::ELEC_evolve(){};
 ELEC_evolve::~ELEC_evolve(){};
 
 double ELEC_evolve::td_force_dt;
-int ELEC_evolve::td_val_elec_01;
-int ELEC_evolve::td_val_elec_02;
-int ELEC_evolve::td_val_elec_03;
 bool ELEC_evolve::td_vext;
 std::vector<int> ELEC_evolve::td_vext_dire_case;
 bool ELEC_evolve::out_dipole;

source/module_hamilt_lcao/module_tddft/ELEC_evolve.h

@@ -30,9 +30,6 @@ class ELEC_evolve
     // fuxiang add 2021-05-25
 
     static double td_force_dt;
-    static int td_val_elec_01;
-    static int td_val_elec_02;
-    static int td_val_elec_03;
     static bool td_vext;
     static std::vector<int> td_vext_dire_case;
     static bool out_dipole;

source/module_io/dipole_io.h

@@ -6,11 +6,14 @@
 namespace ModuleIO
 {
 void write_dipole(const double *rho_save,
-                          const int &is,
-                          const int &istep,
-                          const std::string &fn,
-                          const int &precision = 11,
-                          const bool for_plot = false);
+                  const int &is,
+                  const int &istep,
+                  const std::string &fn,
+                  const int &precision = 11,
+                  const bool for_plot = false);
+
+double prepare(const UnitCell &cell, int &dir);
+
 } // namespace ModuleIO
 
 #endif

source/module_io/input.cpp

@@ -405,9 +405,6 @@ void Input::Default(void)
     // tddft
     //----------------------------------------------------------
     td_force_dt = 0.02;
-    td_val_elec_01 = 1;
-    td_val_elec_02 = 1;
-    td_val_elec_03 = 1;
     td_vext = false;
     td_vext_dire = "1";
 
@@ -1524,18 +1521,6 @@ bool Input::Read(const std::string &fn)
         {
             read_value(ifs, td_force_dt);
         }
-        else if (strcmp("td_val_elec_01", word) == 0)
-        {
-            read_value(ifs, td_val_elec_01);
-        }
-        else if (strcmp("td_val_elec_02", word) == 0)
-        {
-            read_value(ifs, td_val_elec_02);
-        }
-        else if (strcmp("td_val_elec_03", word) == 0)
-        {
-            read_value(ifs, td_val_elec_03);
-        }
         else if (strcmp("td_vext", word) == 0)
         {
             read_value(ifs, td_vext);
@@ -2969,9 +2954,6 @@ void Input::Bcast()
     Parallel_Common::bcast_int(vdw_cutoff_period.y);
     Parallel_Common::bcast_int(vdw_cutoff_period.z);
     // Fuxiang He add 2016-10-26
-    Parallel_Common::bcast_int(td_val_elec_01);
-    Parallel_Common::bcast_int(td_val_elec_02);
-    Parallel_Common::bcast_int(td_val_elec_03);
     Parallel_Common::bcast_double(td_force_dt);
     Parallel_Common::bcast_bool(td_vext);
     Parallel_Common::bcast_string(td_vext_dire);

source/module_io/input.h

@@ -372,9 +372,6 @@ class Input
     // Fuxiang He add 2016-10-26
     //==========================================================
     double td_force_dt; //"fs"
-    int td_val_elec_01; // valence electron 01
-    int td_val_elec_02; // valence electron 02
-    int td_val_elec_03; // valence electron 03
     bool td_vext; // add extern potential or not
     std::string td_vext_dire; // vext direction
     bool out_dipole; // output the dipole or not

source/module_io/input_conv.cpp

@@ -34,20 +34,32 @@ template <typename T> void Input_Conv::parse_expression(const std::string &fn, s
     int count = 0;
     std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)");
     std::vector<std::string> str;
-    std::string::size_type pos1, pos2;
-    std::string c = " ";
-    pos2 = fn.find(c);
-    pos1 = 0;
-    while (std::string::npos != pos2)
-    {
-        str.push_back(fn.substr(pos1, pos2 - pos1));
-        pos1 = pos2 + c.size();
-        pos2 = fn.find(c, pos1);
-    }
-    if (pos1 != fn.length())
-    {
-        str.push_back(fn.substr(pos1));
+    std::stringstream ss(fn);
+    std::string section;
+    while (ss >> section) {
+        int index = 0;
+        if (str.empty()) {
+            while (index < section.size() && std::isspace(section[index])) {
+                index++;
+            }
+        }
+        section.erase(0, index);
+        str.push_back(section);
     }
+    // std::string::size_type pos1, pos2;
+    // std::string c = " ";
+    // pos2 = fn.find(c);
+    // pos1 = 0;
+    // while (std::string::npos != pos2)
+    // {
+    //     str.push_back(fn.substr(pos1, pos2 - pos1));
+    //     pos1 = pos2 + c.size();
+    //     pos2 = fn.find(c, pos1);
+    // }
+    // if (pos1 != fn.length())
+    // {
+    //     str.push_back(fn.substr(pos1));
+    // }
     regex_t reg;
     regcomp(&reg, pattern.c_str(), REG_EXTENDED);
     regmatch_t pmatch[1];
@@ -329,9 +341,6 @@ void Input_Conv::Convert(void)
 //----------------------------------------------------------
 #ifdef __LCAO
     ELEC_evolve::td_force_dt = INPUT.td_force_dt;
-    ELEC_evolve::td_val_elec_01 = INPUT.td_val_elec_01;
-    ELEC_evolve::td_val_elec_02 = INPUT.td_val_elec_02;
-    ELEC_evolve::td_val_elec_03 = INPUT.td_val_elec_03;
     ELEC_evolve::td_vext = INPUT.td_vext;
     if (ELEC_evolve::td_vext)
     {

source/module_io/test/for_testing_input_conv.h

@@ -35,9 +35,6 @@ int hsolver::HSolverLCAO::out_mat_dh = 0;
 int Local_Orbital_Charge::out_dm = 0;
 int Local_Orbital_Charge::out_dm1 = 0;
 double ELEC_evolve::td_force_dt;
-int ELEC_evolve::td_val_elec_01;
-int ELEC_evolve::td_val_elec_02;
-int ELEC_evolve::td_val_elec_03;
 bool ELEC_evolve::td_vext;
 std::vector<int> ELEC_evolve::td_vext_dire_case;
 bool ELEC_evolve::out_dipole;