Skip to content

Fix: Fixed a bug in RT-TDDFT where the electric field was not applied starting from the restart step #5877

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jan 23, 2025
Merged

Fix: Fixed a bug in RT-TDDFT where the electric field was not applied starting from the restart step #5877

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 16 additions & 0 deletions source/module_elecstate/potentials/H_TDDFT_pw.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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");
Expand Down
46 changes: 34 additions & 12 deletions source/module_elecstate/potentials/H_TDDFT_pw.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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;

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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:
Expand All @@ -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;

Expand All @@ -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
Expand All @@ -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);
Expand Down
2 changes: 1 addition & 1 deletion source/module_io/read_set_globalv.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading