From d46b2a1a3c31b39b851a2a1f6d32ec09801d356a Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Thu, 23 Jan 2025 14:43:08 +0800 Subject: [PATCH] Fixed a bug in RT-TDDFT where the efield was not applied starting from the restart step --- .../potentials/H_TDDFT_pw.cpp | 16 +++++++ .../module_elecstate/potentials/H_TDDFT_pw.h | 46 ++++++++++++++----- source/module_io/read_set_globalv.cpp | 2 +- 3 files changed, 51 insertions(+), 13 deletions(-) diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index cf0dbd2d40..4c8104e509 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/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 H_TDDFT_pw::heavi_t0; std::vector 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"); diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h index 8d7c42f9e5..0649487775 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.h +++ b/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 At; // time domain parameters @@ -64,8 +82,8 @@ class H_TDDFT_pw : public PotBase static std::vector gauss_sigma; // time(a.u.) static std::vector gauss_t0; static std::vector 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 gauss_ncut; // trapezoid @@ -76,7 +94,7 @@ class H_TDDFT_pw : public PotBase static std::vector trape_t2; static std::vector trape_t3; static std::vector 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 trape_ncut; // Trigonometric @@ -86,7 +104,7 @@ class H_TDDFT_pw : public PotBase static std::vector trigo_phase1; static std::vector trigo_phase2; static std::vector 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 trigo_ncut; // Heaviside @@ -94,7 +112,7 @@ class H_TDDFT_pw : public PotBase static std::vector heavi_t0; static std::vector 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 &vext_space, int direc); - void cal_v_space_length(std::vector &vext_space, int direc); + void cal_v_space(std::vector& vext_space, int direc); + void cal_v_space_length(std::vector& 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); diff --git a/source/module_io/read_set_globalv.cpp b/source/module_io/read_set_globalv.cpp index 1eb115fd9f..0b3434f170 100644 --- a/source/module_io/read_set_globalv.cpp +++ b/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