diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index 104cd41955..9e87b52dcc 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -343,6 +343,8 @@ - [td\_heavi\_amp](#td_heavi_amp) - [out\_dipole](#out_dipole) - [out\_efield](#out_efield) + - [out\_vecpot](#out_vecpot) + - [init\_vecpot\_file](#init_vecpot_file) - [ocp](#ocp) - [ocp\_set](#ocp_set) - [Variables useful for debugging](#variables-useful-for-debugging) @@ -3262,6 +3264,22 @@ These variables are used to control berry phase and wannier90 interface paramete - False: do not output efield. - **Default**: False +### out_vecpot + +- **Type**: Boolean +- **Description**: output TDDFT Vector potential or not(a.u.) + - True: output Vector potential in file "OUT.suffix/At.dat" + - False: do not output Vector potential. +- **Default**: False + +### init_vecpot_file + +- **Type**: Boolean +- **Description**: Init vector potential through file or not + - True: init vector potential from file "At.dat".(a.u.) It consists of four columns, representing istep and vector potential on each direction. + - False: calculate vector potential by integral of Efield +- **Default**: False + ### ocp - **Type**: Boolean diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index 93601bac40..6f37cd87d3 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -38,7 +38,7 @@ double H_TDDFT_pw::lcut1; double H_TDDFT_pw::lcut2; //velocity gauge -double H_TDDFT_pw::At[3]={0.0,0.0,0.0}; +ModuleBase::Vector3 H_TDDFT_pw::At; // time domain parameters diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h index ea6b4f0ef9..8b3ffcaf42 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.h +++ b/source/module_elecstate/potentials/H_TDDFT_pw.h @@ -54,7 +54,7 @@ class H_TDDFT_pw : public PotBase static double lcut2; //velocity gauge, vector magnetic potential - static double At[3]; + static ModuleBase::Vector3 At; // time domain parameters diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp index 4102aae997..03bfdfdda5 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.cpp @@ -60,7 +60,7 @@ void TDEkinetic, double>>::td_ekinetic_scalar( template void TDEkinetic>::td_ekinetic_grad(std::complex* Hloc, int nnr, ModuleBase::Vector3 grad_overlap){ std::complex tmp= {0, grad_overlap*cart_At}; - Hloc[nnr] += tmp; + Hloc[nnr] -= tmp; return; } @@ -193,10 +193,9 @@ void TDEkinetic>::init_td(void) { TD_Velocity::td_vel_op = &td_velocity; //calculate At in cartesian coorinates. - double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()}; - double (&A)[3] = elecstate::H_TDDFT_pw::At; - cart_At = this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2]; - std::cout << "cart_At: " << cart_At[0] << " " <ucell->a1, this->ucell->a2, this->ucell->a3, elecstate::H_TDDFT_pw::At); + this->cart_At = td_velocity.cart_At; + std::cout<<"cart_At: "< void hamilt::TDNonlocal>::init_td(void) { //calculate At in cartesian coorinates. - double l_norm[3]={this->ucell->a1.norm() ,this->ucell->a2.norm() ,this->ucell->a3.norm()}; - double (&A)[3] = elecstate::H_TDDFT_pw::At; - cart_At = -(this->ucell->a1*A[0]/l_norm[0] + this->ucell->a2*A[1]/l_norm[1] + this->ucell->a3*A[2]/l_norm[2]); - std::cout << "cart_At: " << cart_At[0] << " " <cart_At=TD_Velocity::td_vel_op->cart_At; } // initialize_HR() template @@ -207,7 +204,7 @@ void hamilt::TDNonlocal>::calculate_HR() atom1->iw2n[iw1], tau0 * this->ucell->lat0, T0, - -cart_At, + -cart_At/2.0, 0); #else uot.snap_psibeta_half_tddft(orb, @@ -220,7 +217,7 @@ void hamilt::TDNonlocal>::calculate_HR() atom1->iw2n[iw1], tau0 * this->ucell->lat0, T0, - -cart_At, + -cart_At/2.0, 0); #endif nlm_tot[ad].insert({all_indexes[iw1l], nlm[0]}); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp index 8829bef774..162f7857e5 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.cpp @@ -376,4 +376,4 @@ void sparse_format::clear_zero_elements( } return; -} +} \ No newline at end of file diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h index d7b36f397f..d46a0304d3 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/spar_hsr.h @@ -46,4 +46,4 @@ namespace sparse_format } -#endif +#endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.cpp b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp index 5de3d35209..3056c3d414 100644 --- a/source/module_hamilt_lcao/module_tddft/td_velocity.cpp +++ b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp @@ -1,24 +1,145 @@ -#include "module_base/timer.h" #include "td_velocity.h" +#include "module_base/timer.h" +#include "module_elecstate/potentials/H_TDDFT_pw.h" bool TD_Velocity::tddft_velocity = false; bool TD_Velocity::out_mat_R = false; +bool TD_Velocity::out_vecpot = false; +bool TD_Velocity::init_vecpot_file = false; + TD_Velocity* TD_Velocity::td_vel_op = nullptr; +int TD_Velocity::istep = -1; +int TD_Velocity::max_istep = -1; +std::vector> TD_Velocity::At_from_file; + TD_Velocity::TD_Velocity() { - return; + if (init_vecpot_file && istep == -1) + { + this->read_cart_At(); + } + return; } TD_Velocity::~TD_Velocity() { - this->destroy_HS_R_td_sparse(); + this->destroy_HS_R_td_sparse(); +} + +void TD_Velocity::output_cart_At(const std::string& out_dir) +{ + if (GlobalV::MY_RANK == 0) + { + std::string out_file; + // generate the output file name + out_file = out_dir + "At.dat"; + std::ofstream ofs; + // output title + if (istep == 0) + { + ofs.open(out_file.c_str(), std::ofstream::out); + ofs << std::left << std::setw(8) << "#istep" << std::setw(15) << "A_x" << std::setw(15) << "A_y" + << std::setw(15) << "A_z" << std::endl; + } + else + { + ofs.open(out_file.c_str(), std::ofstream::app); + } + // output the vector potential + ofs << std::left << std::setw(8) << istep; + // divide by 2.0 to get the atomic unit + for (int i = 0; i < 3; i++) + { + ofs << std::scientific << std::setprecision(4) << std::setw(15) << cart_At[i] / 2.0; + } + ofs << std::endl; + ofs.close(); + } + return; } +void TD_Velocity::cal_cart_At(const ModuleBase::Vector3& a0, + const ModuleBase::Vector3& a1, + const ModuleBase::Vector3& a2, + const ModuleBase::Vector3& At) +{ + istep++; + if (init_vecpot_file) + { + this->cart_At = At_from_file[istep > max_istep ? max_istep : istep]; + } + else + { + const double l_norm[3] = {a0.norm(), a1.norm(), a2.norm()}; + this->cart_At = a0 * At[0] / l_norm[0] + a1 * At[1] / l_norm[1] + a2 * At[2] / l_norm[2]; + } + // output the vector potential if needed + if (out_vecpot == true) + { + this->output_cart_At(GlobalV::global_out_dir); + } +} + +void TD_Velocity::read_cart_At(void) +{ + std::string in_file; + // generate the input file name + in_file = "At.dat"; + std::ifstream ifs(in_file.c_str()); + // check if the file is exist + if (!ifs) + { + ModuleBase::WARNING_QUIT("TD_Velocity::read_cart_At", "Cannot open Vector potential file!"); + } + std::string line; + std::vector str_vec; + // use tmp to skip the istep number + int tmp = 0; + while (std::getline(ifs, line)) + { + // A tmporary vector3 to store the data of this line + ModuleBase::Vector3 At; + if (line[0] == '#') + { + continue; + } + std::istringstream iss(line); + // skip the istep number + if (!(iss >> tmp)) + { + ModuleBase::WARNING_QUIT("TD_Velocity::read_cart_At", "Error reading istep!"); + } + // read the vector potential + double component = 0; + // Read three components + for (int i = 0; i < 3; i++) + { + if (!(iss >> component)) + { + ModuleBase::WARNING_QUIT("TD_Velocity::read_cart_At", + "Error reading component " + std::to_string(i + 1) + " for istep " + + std::to_string(tmp) + "!"); + } + At[i] = component; + } + // unit transform ,change unit into Ry/bohr/e*t_a.u. + At *= 2.0; + // add the tmporary vector3 to the vector potential vector + At_from_file.push_back(At); + } + // set the max_istep + max_istep = At_from_file.size() - 1; + ifs.close(); + + return; +} void TD_Velocity::destroy_HS_R_td_sparse(void) { - std::map, std::map>>> empty_HR_sparse_td_vel_up; - std::map, std::map>>> empty_HR_sparse_td_vel_down; - HR_sparse_td_vel[0].swap(empty_HR_sparse_td_vel_up); - HR_sparse_td_vel[1].swap(empty_HR_sparse_td_vel_down); -} \ No newline at end of file + std::map, std::map>>> + empty_HR_sparse_td_vel_up; + std::map, std::map>>> + empty_HR_sparse_td_vel_down; + HR_sparse_td_vel[0].swap(empty_HR_sparse_td_vel_up); + HR_sparse_td_vel[1].swap(empty_HR_sparse_td_vel_down); +} diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.h b/source/module_hamilt_lcao/module_tddft/td_velocity.h index 0d274228dd..ddfaaf87cf 100644 --- a/source/module_hamilt_lcao/module_tddft/td_velocity.h +++ b/source/module_hamilt_lcao/module_tddft/td_velocity.h @@ -1,17 +1,19 @@ #ifndef TD_VELOCITY_H #define TD_VELOCITY_H -#include - -#include "module_base/timer.h" #include "module_base/abfs-vector3_order.h" -//Class to store TDDFT velocity gague infos. +#include "module_base/timer.h" + +#include +// Class to store TDDFT velocity gague infos. class TD_Velocity { public: TD_Velocity(); ~TD_Velocity(); - /// @brief Judge if in tddft calculation or not + void init(); + + /// @brief Judge if in tddft calculation or not static bool tddft_velocity; /// @brief switch to control the output of HR @@ -20,15 +22,42 @@ class TD_Velocity /// @brief pointer to the only TD_Velocity object itself static TD_Velocity* td_vel_op; - //For TDDFT velocity gague, to fix the output of HR - std::map, std::map>>> HR_sparse_td_vel[2]; + /// @brief switch to control the output of At + static bool out_vecpot; + + /// @brief switch to control the source of At + static bool init_vecpot_file; + /// @brief Store the vector potential for tddft calculation + ModuleBase::Vector3 cart_At; + + /// @brief calculate the At in cartesian coordinate + void cal_cart_At(const ModuleBase::Vector3& a0, + const ModuleBase::Vector3& a1, + const ModuleBase::Vector3& a2, + const ModuleBase::Vector3& At); + + // For TDDFT velocity gague, to fix the output of HR + std::map, std::map>>> HR_sparse_td_vel[2]; private: + /// @brief read At from output file + void read_cart_At(void); + + /// @brief output cart_At to output file + void output_cart_At(const std::string& out_dir); + + /// @brief store isteps now + static int istep; + + /// @brief total steps of read in At + static int max_istep; + + /// @brief store the read in At_data + static std::vector> At_from_file; /// @brief destory HSR data stored void destroy_HS_R_td_sparse(void); - }; -#endif \ No newline at end of file +#endif diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 0a70b615c5..64fae37c81 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -470,6 +470,8 @@ void Input::Default(void) out_dipole = false; out_efield = false; out_current = false; + out_vecpot = false; + init_vecpot_file = false; td_print_eij = -1.0; td_edm = 0; @@ -1813,6 +1815,14 @@ bool Input::Read(const std::string& fn) { read_value(ifs, out_efield); } + else if (strcmp("out_vecpot", word) == 0) + { + read_value(ifs, out_vecpot); + } + else if (strcmp("init_vecpot_file", word) == 0) + { + read_value(ifs, init_vecpot_file); + } else if (strcmp("td_print_eij", word) == 0) { read_value(ifs, td_print_eij); @@ -3752,6 +3762,8 @@ void Input::Bcast() Parallel_Common::bcast_bool(out_dipole); Parallel_Common::bcast_bool(out_efield); Parallel_Common::bcast_bool(out_current); + Parallel_Common::bcast_bool(out_vecpot); + Parallel_Common::bcast_bool(init_vecpot_file); Parallel_Common::bcast_double(td_print_eij); Parallel_Common::bcast_int(td_edm); Parallel_Common::bcast_bool(test_skip_ewald); diff --git a/source/module_io/input.h b/source/module_io/input.h index fedc85dcd0..3aad48b7ce 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -416,6 +416,8 @@ class Input bool out_dipole; // output the dipole or not bool out_efield; // output the efield or not bool out_current; //output the current or not + bool out_vecpot; // output the vector potential or not + bool init_vecpot_file; // initialize the vector potential, though file or integral double td_print_eij; // threshold to output Eij elements int td_edm; //0: new edm method 1: old edm method diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 72d49acc87..b2685ed702 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -550,6 +550,8 @@ void Input_Conv::Convert(void) module_tddft::Evolve_elec::out_current = INPUT.out_current; module_tddft::Evolve_elec::td_print_eij = INPUT.td_print_eij; module_tddft::Evolve_elec::td_edm = INPUT.td_edm; + TD_Velocity::out_vecpot = INPUT.out_vecpot; + TD_Velocity::init_vecpot_file = INPUT.init_vecpot_file; read_td_efield(); #endif diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h index d4449bdbf1..02d8c59484 100644 --- a/source/module_io/test/for_testing_input_conv.h +++ b/source/module_io/test/for_testing_input_conv.h @@ -16,6 +16,7 @@ #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_lcao/module_tddft/evolve_elec.h" +#include "module_hamilt_lcao/module_tddft/td_velocity.h" #include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h" #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" #include "module_hamilt_pw/hamilt_pwdft/wavefunc.h" @@ -41,6 +42,8 @@ std::vector module_tddft::Evolve_elec::td_vext_dire_case; bool module_tddft::Evolve_elec::out_dipole; bool module_tddft::Evolve_elec::out_efield; bool module_tddft::Evolve_elec::out_current; +bool TD_Velocity::out_vecpot; +bool TD_Velocity::init_vecpot_file; double module_tddft::Evolve_elec::td_print_eij; int module_tddft::Evolve_elec::td_edm; double elecstate::Gatefield::zgate = 0.5; diff --git a/source/module_io/write_HS_sparse.cpp b/source/module_io/write_HS_sparse.cpp index c7e7d5475a..2ca1d58049 100644 --- a/source/module_io/write_HS_sparse.cpp +++ b/source/module_io/write_HS_sparse.cpp @@ -849,4 +849,4 @@ template void ModuleIO::save_sparse>( const Parallel_Orbitals&, const std::string&, const int&, - const bool&); + const bool&); \ No newline at end of file