Fixed a bug in RT-TDDFT where the efield was not applied starting from the restart step (#5877)

Author: AsTonyshment

source/module_elecstate/potentials/H_TDDFT_pw.cpp

@@ -12,6 +12,7 @@ namespace elecstate
 {
 
 int H_TDDFT_pw::istep = -1;
+bool H_TDDFT_pw::is_initialized = false;
 
 double H_TDDFT_pw::amp;
 double H_TDDFT_pw::bmod;
@@ -76,6 +77,21 @@ int H_TDDFT_pw::heavi_count;
 std::vector<double> H_TDDFT_pw::heavi_t0;
 std::vector<double> H_TDDFT_pw::heavi_amp; // Ry/bohr
 
+void H_TDDFT_pw::current_step_info(const std::string& file_dir, int& istep)
+{
+    std::stringstream ssc;
+    ssc << file_dir << "Restart_md.dat";
+    std::ifstream file(ssc.str().c_str());
+
+    if (!file)
+    {
+        ModuleBase::WARNING_QUIT("H_TDDFT_pw::current_step_info", "No Restart_md.dat!");
+    }
+
+    file >> istep;
+    file.close();
+}
+
 void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo)
 {
     ModuleBase::TITLE("H_TDDFT_pw", "cal_fixed_v");

source/module_elecstate/potentials/H_TDDFT_pw.h

@@ -2,6 +2,7 @@
 #define H_TDDFT_PW_H
 
 #include "module_io/input_conv.h"
+#include "module_parameter/parameter.h" // PARAM.globalv.global_readin_dir, PARAM.inp.mdp.md_restart
 #include "pot_base.h"
 
 namespace elecstate
@@ -16,8 +17,25 @@ class H_TDDFT_pw : public PotBase
         this->fixed_mode = true;
 
         this->rho_basis_ = rho_basis_in;
+
+        // If it is the first time to create an H_TDDFT_pw instance and is restart calculation,
+        // initialize istep using current_step_info
+        if (!is_initialized && PARAM.inp.mdp.md_restart)
+        {
+            int restart_istep = -1;
+            std::string file_dir = PARAM.globalv.global_readin_dir;
+            current_step_info(file_dir, restart_istep);
+
+            if (restart_istep >= 0)
+            {
+                H_TDDFT_pw::istep = restart_istep - 1; // Update istep
+            }
+
+            is_initialized = true; // Mark as initialized, so that istep will not be initialized again
+        }
     }
-    ~H_TDDFT_pw(){};
+
+    ~H_TDDFT_pw() {};
 
     void cal_fixed_v(double* vl_pseudo) override;
 
@@ -42,17 +60,17 @@ class H_TDDFT_pw : public PotBase
     static int tstart;
     static int tend;
     static double dt;
-    //cut dt for integral
+    // cut dt for integral
     static double dt_int;
     static int istep_int;
 
     // space domain parameters
 
-    //length gauge
+    // length gauge
     static double lcut1;
     static double lcut2;
 
-    //velocity gauge, vector magnetic potential
+    // velocity gauge, vector magnetic potential
     static ModuleBase::Vector3<double> At;
 
     // time domain parameters
@@ -64,8 +82,8 @@ class H_TDDFT_pw : public PotBase
     static std::vector<double> gauss_sigma; // time(a.u.)
     static std::vector<double> gauss_t0;
     static std::vector<double> gauss_amp; // Ry/bohr
-    //add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough
-    //must be even, thus would get odd number of points for simpson integral
+    // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough
+    // must be even, thus would get odd number of points for simpson integral
     static std::vector<int> gauss_ncut;
 
     // trapezoid
@@ -76,7 +94,7 @@ class H_TDDFT_pw : public PotBase
     static std::vector<double> trape_t2;
     static std::vector<double> trape_t3;
     static std::vector<double> trape_amp; // Ry/bohr
-    //add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough
+    // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough
     static std::vector<int> trape_ncut;
 
     // Trigonometric
@@ -86,15 +104,15 @@ class H_TDDFT_pw : public PotBase
     static std::vector<double> trigo_phase1;
     static std::vector<double> trigo_phase2;
     static std::vector<double> trigo_amp; // Ry/bohr
-    //add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough
+    // add for velocity gauge, recut dt into n pieces to make sure the integral is accurate Enough
     static std::vector<int> trigo_ncut;
 
     // Heaviside
     static int heavi_count;
     static std::vector<double> heavi_t0;
     static std::vector<double> heavi_amp; // Ry/bohr
 
-    //update At for velocity gauge by intergral of E(t)dt
+    // update At for velocity gauge by intergral of E(t)dt
     static void update_At();
 
   private:
@@ -103,6 +121,7 @@ class H_TDDFT_pw : public PotBase
     // Vext will evolve by time, every time cal_fixed_v() is called, istep++
     //------------------------
     static int istep;
+    static bool is_initialized; // static flag variable, used to ensure initialization only once
 
     static double amp;
 
@@ -111,9 +130,12 @@ class H_TDDFT_pw : public PotBase
 
     const UnitCell* ucell_ = nullptr;
 
+    // Obtain the current MD step information, used for restart calculation
+    void current_step_info(const std::string& file_dir, int& istep);
+
     // potential of electric field in space domain : length gauge and velocity gauge
-    void cal_v_space(std::vector<double> &vext_space, int direc);
-    void cal_v_space_length(std::vector<double> &vext_space, int direc);
+    void cal_v_space(std::vector<double>& vext_space, int direc);
+    void cal_v_space_length(std::vector<double>& vext_space, int direc);
     double cal_v_space_length_potential(double i);
 
     // potential of electric field in time domain : Gauss , trapezoid, trigonometric, heaviside,  HHG
@@ -124,7 +146,7 @@ class H_TDDFT_pw : public PotBase
     static double cal_v_time_heaviside();
     // double cal_v_time_HHG();
 
-    //get ncut number for At integral
+    // get ncut number for At integral
     static int check_ncut(int t_type);
 
     void prepare(const ModuleBase::Matrix3& G, int& dir);

source/module_io/read_set_globalv.cpp

@@ -23,7 +23,7 @@ void ReadInput::set_globalv(Parameter& para)
         para.sys.global_matrix_dir = to_dir(para.sys.global_matrix_dir);
 
         /// get the global readin directory
-        para.sys.global_readin_dir = para.inp.read_file_dir + '/';
+        para.sys.global_readin_dir = para.inp.read_file_dir;
         para.sys.global_readin_dir = to_dir(para.sys.global_readin_dir);
 
         /// get the stru file for md restart case