From 13496ca412920cc057db49301872a51fa300a206 Mon Sep 17 00:00:00 2001 From: esroamer Date: Thu, 4 Jul 2024 22:40:37 +0800 Subject: [PATCH 1/2] Length gague current --- source/Makefile.Objects | 1 + .../module_esolver/esolver_ks_lcao_tddft.cpp | 7 +- .../hamilt_lcaodft/hamilt_lcao.cpp | 3 +- .../operator_lcao/td_ekinetic_lcao.cpp | 33 +- .../operator_lcao/td_ekinetic_lcao.h | 6 +- .../operator_lcao/td_nonlocal_lcao.cpp | 4 +- .../module_tddft/CMakeLists.txt | 1 + .../module_tddft/td_current.cpp | 506 ++++++++++++++++++ .../module_tddft/td_current.h | 76 +++ .../module_tddft/td_velocity.cpp | 13 +- .../module_tddft/td_velocity.h | 5 +- source/module_io/td_current_io.cpp | 68 +-- source/module_io/td_current_io.h | 3 +- 13 files changed, 648 insertions(+), 78 deletions(-) create mode 100644 source/module_hamilt_lcao/module_tddft/td_current.cpp create mode 100644 source/module_hamilt_lcao/module_tddft/td_current.h diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 995919ce6b0..d6de50337ae 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -501,6 +501,7 @@ OBJS_LCAO=DM_gamma.o\ norm_psi.o\ propagator.o\ td_velocity.o\ + td_current.o\ snap_psibeta_half_tddft.o\ upsi.o\ FORCE_STRESS.o\ diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 94f5812578c..287750cf8b7 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -401,17 +401,14 @@ void ESolver_KS_LCAO_TDDFT::after_scf(const int istep) { } if (TD_Velocity::out_current == true) { - elecstate::DensityMatrix, double>* tmp_DM - = dynamic_cast>*>(this->pelec)->get_DM(); - hamilt::HamiltLCAO, double>* p_ham_lcao = dynamic_cast, double>*>(this->p_hamilt); + elecstate::DensityMatrix, double>* tmp_DM = dynamic_cast>*>(this->pelec)->get_DM(); ModuleIO::write_current(istep, this->psi, pelec, kv, two_center_bundle_.overlap_orb.get(), tmp_DM->get_paraV_pointer(), - this->RA, - p_ham_lcao->getSR()); // mohan add 2024-04-02 + this->RA); // mohan add 2024-04-02 } ESolver_KS_LCAO, double>::after_scf(istep); } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp index 0a71ebff5ea..5c5a8d21154 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.cpp @@ -339,11 +339,10 @@ HamiltLCAO::HamiltLCAO(Gint_Gamma* GG_in, // TDDFT_velocity_gague if (TD_Velocity::tddft_velocity) { - elecstate::H_TDDFT_pw::update_At(); + if(!TD_Velocity::init_vecpot_file) elecstate::H_TDDFT_pw::update_At(); Operator* td_ekinetic = new TDEkinetic>(LM_in, this->hR, &(this->getHk(LM_in)), - this->sR, kv, &GlobalC::ucell, &GlobalC::GridD, 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 fd740b3dcb3..9ede8ddb1d5 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 @@ -17,13 +17,12 @@ template TDEkinetic>::TDEkinetic(LCAO_Matrix* LM_in, hamilt::HContainer* hR_in, std::vector* hK_in, - hamilt::HContainer* SR_in, const K_Vectors* kv_in, const UnitCell* ucell_in, Grid_Driver* GridD_in, const Parallel_Orbitals* paraV, const TwoCenterIntegrator* intor) - : SR(SR_in), kv(kv_in), OperatorLCAO(LM_in, kv_in->kvec_d, hR_in, hK_in), intor_(intor) + : kv(kv_in), OperatorLCAO(LM_in, kv_in->kvec_d, hR_in, hK_in), intor_(intor) { this->LM = LM_in; this->ucell = ucell_in; @@ -44,18 +43,18 @@ TDEkinetic>::~TDEkinetic() } // term A^2*S template -void TDEkinetic>::td_ekinetic_scalar(std::complex* Hloc, TR* Sloc, int nnr) +void TDEkinetic>::td_ekinetic_scalar(std::complex* Hloc, TR overlap, int nnr) { return; } // term A^2*S template <> void TDEkinetic, double>>::td_ekinetic_scalar(std::complex* Hloc, - double* Sloc, + double overlap, int nnr) { // the correction term A^2/2 , 4.0 due to the unit transformation - std::complex tmp = {cart_At.norm2() * Sloc[nnr] / 4.0, 0}; + std::complex tmp = {cart_At.norm2() * overlap, 0}; Hloc[nnr] += tmp; return; } @@ -69,7 +68,7 @@ void TDEkinetic>::td_ekinetic_grad(std::complex* Hl //βˆ‡ refer to the integral βˆ«πœ™(π‘Ÿ)πœ•/πœ•π‘Ÿπœ™(π‘Ÿβˆ’π‘…)π‘‘π‘Ÿ,but abacus only provide the integral of βˆ«πœ™(π‘Ÿ)πœ•/πœ•Rπœ™(π‘Ÿβˆ’π‘…)π‘‘π‘Ÿ. An extra //minus must be counted in. The final term is iA dot βˆ‡R std::complex tmp = {0, grad_overlap * cart_At}; - Hloc[nnr] += tmp; + Hloc[nnr] += tmp * 2.0; return; } @@ -102,7 +101,6 @@ void TDEkinetic>::calculate_HR() ModuleBase::Vector3 dtau = this->ucell->cal_dtau(iat1, iat2, R_index2); hamilt::BaseMatrix>* tmp = this->hR_tmp->find_matrix(iat1, iat2, R_index2); - hamilt::BaseMatrix* tmps = this->SR->find_matrix(iat1, iat2, R_index2); if (tmp != nullptr) { if (TD_Velocity::out_current) @@ -110,14 +108,13 @@ void TDEkinetic>::calculate_HR() std::complex* tmp_c[3] = {nullptr, nullptr, nullptr}; for (int i = 0; i < 3; i++) { - tmp_c[i] - = td_velocity.get_current_term_pointer(i)->find_matrix(iat1, iat2, R_index2)->get_pointer(); + tmp_c[i] = td_velocity.get_current_term_pointer(i)->find_matrix(iat1, iat2, R_index2)->get_pointer(); } - this->cal_HR_IJR(iat1, iat2, paraV, dtau, tmp->get_pointer(), tmp_c, tmps->get_pointer()); + this->cal_HR_IJR(iat1, iat2, paraV, dtau, tmp->get_pointer(), tmp_c); } else { - this->cal_HR_IJR(iat1, iat2, paraV, dtau, tmp->get_pointer(), nullptr, tmps->get_pointer()); + this->cal_HR_IJR(iat1, iat2, paraV, dtau, tmp->get_pointer(), nullptr); } } else @@ -135,8 +132,7 @@ void TDEkinetic>::cal_HR_IJR(const int& iat1, const Parallel_Orbitals* paraV, const ModuleBase::Vector3& dtau, std::complex* data_pointer, - std::complex** data_pointer_c, - TR* s_pointer) + std::complex** data_pointer_c) { const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); // --------------------------------------------- @@ -171,6 +167,7 @@ void TDEkinetic>::cal_HR_IJR(const int& iat1, // calculate the Ekinetic matrix for each pair of orbitals // --------------------------------------------- double grad[3] = {0, 0, 0}; + double overlap = 0; auto row_indexes = paraV->get_indexes_row(iat1); auto col_indexes = paraV->get_indexes_col(iat2); const int step_trace = col_indexes.size() + 1; @@ -195,17 +192,16 @@ void TDEkinetic>::cal_HR_IJR(const int& iat1, int M2 = (m2 % 2 == 0) ? -m2 / 2 : (m2 + 1) / 2; // calculate , which equals to -. - intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, dtau * this->ucell->lat0, nullptr, grad); + intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, dtau * this->ucell->lat0, &overlap, grad); ModuleBase::Vector3 grad_overlap(grad[0], grad[1], grad[2]); for (int ipol = 0; ipol < npol; ipol++) { // key change - td_ekinetic_scalar(data_pointer, s_pointer, ipol * step_trace); + td_ekinetic_scalar(data_pointer, overlap, ipol * step_trace); td_ekinetic_grad(data_pointer, ipol * step_trace, grad_overlap); } data_pointer += npol; - s_pointer += npol; // current grad part if (data_pointer_c != nullptr) { @@ -216,13 +212,14 @@ void TDEkinetic>::cal_HR_IJR(const int& iat1, // part of Momentum operator, -iβˆ‡r,used to calculate the current // here is actually iβˆ‡R data_pointer_c[dir][ipol * step_trace] += std::complex(0, grad_overlap[dir]); + // part of Momentum operator, eA,used to calculate the current + data_pointer_c[dir][ipol * step_trace] += std::complex(overlap * cart_At[dir], 0); } data_pointer_c[dir] += npol; } } } data_pointer += (npol - 1) * col_indexes.size(); - s_pointer += (npol - 1) * col_indexes.size(); if (data_pointer_c != nullptr) { for (int dir = 0; dir < 3; dir++) @@ -236,7 +233,7 @@ void TDEkinetic>::init_td(void) { TD_Velocity::td_vel_op = &td_velocity; // calculate At in cartesian coorinates. - td_velocity.cal_cart_At(this->ucell->a1, this->ucell->a2, this->ucell->a3, elecstate::H_TDDFT_pw::At); + td_velocity.cal_cart_At(elecstate::H_TDDFT_pw::At); this->cart_At = td_velocity.cart_At; std::cout << "cart_At: " << cart_At[0] << " " << cart_At[1] << " " << cart_At[2] << std::endl; } diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h index 00365e2e562..78de85c263b 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_ekinetic_lcao.h @@ -36,7 +36,6 @@ class TDEkinetic> : public OperatorLCAO TDEkinetic>(LCAO_Matrix* LM_in, hamilt::HContainer* hR_in, std::vector* hK_in, - hamilt::HContainer* SR_in, const K_Vectors* kv_in, const UnitCell* ucell_in, Grid_Driver* GridD_in, @@ -70,8 +69,7 @@ class TDEkinetic> : public OperatorLCAO const Parallel_Orbitals* paraV, const ModuleBase::Vector3& dtau, std::complex* data_pointer, - std::complex** data_pointer_c, - TR* s_pointer); + std::complex** data_pointer_c); /** * @brief calculate the ekinetic matrix correction term in tddft with specific atom-pairs @@ -94,7 +92,7 @@ class TDEkinetic> : public OperatorLCAO const K_Vectors* kv; /// @brief correction term iAβ‹…βˆ‡ - void td_ekinetic_scalar(std::complex* Hloc, TR* Sloc, int nnr); + void td_ekinetic_scalar(std::complex* Hloc, TR Sloc, int nnr); /// @brief correction term A^2*S void td_ekinetic_grad(std::complex* Hloc, int nnr, ModuleBase::Vector3 grad_overlap); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp index e62b1e73d78..c2a8aa1ee64 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/operator_lcao/td_nonlocal_lcao.cpp @@ -198,7 +198,7 @@ void hamilt::TDNonlocal>::calculate_HR() atom1->iw2n[iw1], tau0 * this->ucell->lat0, T0, - cart_At / 2.0, + cart_At, TD_Velocity::out_current); for (int dir = 0; dir < nlm_dim; dir++) { @@ -412,7 +412,7 @@ void hamilt::TDNonlocal>::contributeHR() static_cast*>(this->next_sub_op)->set_HR_fixed(this->hR_tmp); } // calculate the values in hR_tmp - this->calculate_HR(); + //this->calculate_HR(); this->hR_tmp_done = true; } ModuleBase::timer::tick("TDNonlocal", "contributeHR"); diff --git a/source/module_hamilt_lcao/module_tddft/CMakeLists.txt b/source/module_hamilt_lcao/module_tddft/CMakeLists.txt index 12a76e0c740..2cde80aa228 100644 --- a/source/module_hamilt_lcao/module_tddft/CMakeLists.txt +++ b/source/module_hamilt_lcao/module_tddft/CMakeLists.txt @@ -8,6 +8,7 @@ if(ENABLE_LCAO) propagator.cpp upsi.cpp td_velocity.cpp + td_current.cpp snap_psibeta_half_tddft.cpp ) diff --git a/source/module_hamilt_lcao/module_tddft/td_current.cpp b/source/module_hamilt_lcao/module_tddft/td_current.cpp new file mode 100644 index 00000000000..cc823ee6685 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/td_current.cpp @@ -0,0 +1,506 @@ +#include "td_current.h" +#ifdef __LCAO +#include "module_base/timer.h" +#include "module_base/tool_title.h" +#include "module_hamilt_lcao/module_tddft/snap_psibeta_half_tddft.h" +#ifdef _OPENMP +#include +#endif + +TD_current::TD_current(const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const Parallel_Orbitals* paraV, + const TwoCenterIntegrator* intor) + : ucell(ucell_in), Grid(GridD_in), paraV(paraV) , intor_(intor) +{ + // for length gague, the A(t) = 0 for all the time. + this->cart_At = ModuleBase::Vector3(0,0,0); + this->initialize_vcomm_r(GridD_in, paraV); + this->initialize_grad_term(GridD_in, paraV); +} +TD_current::~TD_current() +{ + for (int dir=0;dir<3;dir++) + { + delete this->current_term[dir]; + } +} +//allocate space for current_term +void TD_current::initialize_vcomm_r(Grid_Driver* GridD, const Parallel_Orbitals* paraV) +{ + ModuleBase::TITLE("TD_current", "initialize_vcomm_r"); + ModuleBase::timer::tick("TD_current", "initialize_vcomm_r"); + for (int dir=0;dir<3;dir++) + { + if (this->current_term[dir] == nullptr) + this->current_term[dir] = new hamilt::HContainer>(paraV); + } + + this->adjs_vcommr.clear(); + this->adjs_vcommr.reserve(this->ucell->nat); + for (int iat0 = 0; iat0 < ucell->nat; iat0++) + { + auto tau0 = ucell->get_tau(iat0); + int T0, I0; + ucell->iat2iait(iat0, &I0, &T0); + AdjacentAtomInfo adjs; + GridD->Find_atom(*ucell, tau0, T0, I0, &adjs); + std::vector is_adj(adjs.adj_num + 1, false); + for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + { + const int T1 = adjs.ntype[ad1]; + const int I1 = adjs.natom[ad1]; + const int iat1 = ucell->itia2iat(T1, I1); + const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad1]; + const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; + // choose the real adjacent atoms + const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); + // Note: the distance of atoms should less than the cutoff radius, + // When equal, the theoretical value of matrix element is zero, + // but the calculated value is not zero due to the numerical error, which would lead to result changes. + if (this->ucell->cal_dtau(iat0, iat1, R_index1).norm() * this->ucell->lat0 + < orb.Phi[T1].getRcut() + this->ucell->infoNL.Beta[T0].get_rcut_max()) + { + is_adj[ad1] = true; + } + } + filter_adjs(is_adj, adjs); + this->adjs_vcommr.push_back(adjs); + for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + { + const int T1 = adjs.ntype[ad1]; + const int I1 = adjs.natom[ad1]; + const int iat1 = ucell->itia2iat(T1, I1); + const ModuleBase::Vector3& R_index1 = adjs.box[ad1]; + for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2) + { + const int T2 = adjs.ntype[ad2]; + const int I2 = adjs.natom[ad2]; + const int iat2 = ucell->itia2iat(T2, I2); + ModuleBase::Vector3& R_index2 = adjs.box[ad2]; + if (paraV->get_col_size(iat2) <= 0 || paraV->get_row_size(iat1) <= 0) + { + continue; + } + hamilt::AtomPair> tmp(iat1, + iat2, + R_index2.x - R_index1.x, + R_index2.y - R_index1.y, + R_index2.z - R_index1.z, + paraV); + for (int dir=0;dir<3;dir++) + { + this->current_term[dir]->insert_pair(tmp); + } + } + } + } + // allocate the memory of BaseMatrix in cal_vcomm_r_IJR, and set the new values to zero + for (int dir=0;dir<3;dir++) + { + this->current_term[dir]->allocate(nullptr, true); + } + ModuleBase::timer::tick("TD_current", "initialize_vcomm_r"); +} +void TD_current::initialize_grad_term(Grid_Driver* GridD, const Parallel_Orbitals* paraV) +{ + ModuleBase::TITLE("TD_current", "initialize_grad_term"); + ModuleBase::timer::tick("TD_current", "initialize_grad_term"); + + for (int dir=0;dir<3;dir++) + { + if (this->current_term[dir] == nullptr) + this->current_term[dir] = new hamilt::HContainer>(paraV); + } + for (int iat1 = 0; iat1 < ucell->nat; iat1++) + { + auto tau1 = ucell->get_tau(iat1); + int T1, I1; + ucell->iat2iait(iat1, &I1, &T1); + AdjacentAtomInfo adjs; + GridD->Find_atom(*ucell, tau1, T1, I1, &adjs); + std::vector is_adj(adjs.adj_num + 1, false); + for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + { + const int T2 = adjs.ntype[ad1]; + const int I2 = adjs.natom[ad1]; + const int iat2 = ucell->itia2iat(T2, I2); + if (paraV->get_row_size(iat1) <= 0 || paraV->get_col_size(iat2) <= 0) + { + continue; + } + const ModuleBase::Vector3& R_index2 = adjs.box[ad1]; + // choose the real adjacent atoms + const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); + // Note: the distance of atoms should less than the cutoff radius, + // When equal, the theoretical value of matrix element is zero, + // but the calculated value is not zero due to the numerical error, which would lead to result changes. + if (this->ucell->cal_dtau(iat1, iat2, R_index2).norm() * this->ucell->lat0 + < orb.Phi[T1].getRcut() + orb.Phi[T2].getRcut()) + { + is_adj[ad1] = true; + } + } + filter_adjs(is_adj, adjs); + this->adjs_grad.push_back(adjs); + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + int iat2 = ucell->itia2iat(T2, I2); + ModuleBase::Vector3& R_index = adjs.box[ad]; + hamilt::AtomPair> tmp(iat1, iat2, R_index.x, R_index.y, R_index.z, paraV); + for (int dir=0;dir<3;dir++) + { + this->current_term[dir]->insert_pair(tmp); + } + } + } + // allocate the memory of BaseMatrix in HR, and set the new values to zero + for (int dir=0;dir<3;dir++) + { + this->current_term[dir]->allocate(nullptr, true); + } + + ModuleBase::timer::tick("EkineticNew", "initialize_HR"); +} + +void TD_current::calculate_vcomm_r() +{ + ModuleBase::TITLE("TD_current", "calculate_vcomm_r"); + ModuleBase::timer::tick("TD_current", "calculate_vcomm_r"); + + const Parallel_Orbitals* paraV = this->current_term[0]->get_atom_pair(0).get_paraV(); + const int npol = this->ucell->get_npol(); + // 1. calculate for each pair of atoms +#ifdef _OPENMP +#pragma omp parallel +{ + std::unordered_set atom_row_list; + #pragma omp for + for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) + { + atom_row_list.insert(iat0); + } +#endif + for (int iat0 = 0; iat0 < this->ucell->nat; iat0++) + { + auto tau0 = ucell->get_tau(iat0); + int T0, I0; + ucell->iat2iait(iat0, &I0, &T0); + AdjacentAtomInfo& adjs = this->adjs_vcommr[iat0]; + + std::vector>>>> nlm_tot; + nlm_tot.resize(adjs.adj_num + 1); + for (int i = 0; i < adjs.adj_num + 1; i++) + { + nlm_tot[i].resize(4); + } + + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T1 = adjs.ntype[ad]; + const int I1 = adjs.natom[ad]; + const int iat1 = ucell->itia2iat(T1, I1); + const ModuleBase::Vector3& tau1 = adjs.adjacent_tau[ad]; + const Atom* atom1 = &ucell->atoms[T1]; + + const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); + auto all_indexes = paraV->get_indexes_row(iat1); +#ifdef _OPENMP + if(atom_row_list.find(iat1) == atom_row_list.end()) + { + all_indexes.clear(); + } +#endif + auto col_indexes = paraV->get_indexes_col(iat1); + // insert col_indexes into all_indexes to get universal set with no repeat elements + all_indexes.insert(all_indexes.end(), col_indexes.begin(), col_indexes.end()); + std::sort(all_indexes.begin(), all_indexes.end()); + all_indexes.erase(std::unique(all_indexes.begin(), all_indexes.end()), all_indexes.end()); + for (int iw1l = 0; iw1l < all_indexes.size(); iw1l += npol) + { + const int iw1 = all_indexes[iw1l] / npol; + std::vector>> nlm; + // nlm is a vector of vectors, but size of outer vector is only 1 when out_current is false + // and size of outer vector is 4 when out_current is true (3 for , 1 for + // ) inner loop : all projectors (L0,M0) + + // snap_psibeta_half_tddft() are used to calculate + // and as well if current are needed + + module_tddft::snap_psibeta_half_tddft(orb, + this->ucell->infoNL, + nlm, + tau1 * this->ucell->lat0, + T1, + atom1->iw2l[iw1], + atom1->iw2m[iw1], + atom1->iw2n[iw1], + tau0 * this->ucell->lat0, + T0, + this->cart_At, + true); + for (int dir = 0; dir < 4; dir++) + { + nlm_tot[ad][dir].insert({all_indexes[iw1l], nlm[dir]}); + } + } + } +// 2. calculate D for each pair of atoms + for (int ad1 = 0; ad1 < adjs.adj_num + 1; ++ad1) + { + const int T1 = adjs.ntype[ad1]; + const int I1 = adjs.natom[ad1]; + const int iat1 = ucell->itia2iat(T1, I1); +#ifdef _OPENMP + if(atom_row_list.find(iat1) == atom_row_list.end()) + { + continue; + } +#endif + ModuleBase::Vector3& R_index1 = adjs.box[ad1]; + for (int ad2 = 0; ad2 < adjs.adj_num + 1; ++ad2) + { + const int T2 = adjs.ntype[ad2]; + const int I2 = adjs.natom[ad2]; + const int iat2 = ucell->itia2iat(T2, I2); + ModuleBase::Vector3& R_index2 = adjs.box[ad2]; + ModuleBase::Vector3 R_vector(R_index2[0] - R_index1[0], + R_index2[1] - R_index1[1], + R_index2[2] - R_index1[2]); + std::complex* tmp_c[3] = {nullptr, nullptr, nullptr}; + for (int i = 0; i < 3; i++) + { + tmp_c[i] = this->current_term[i]->find_matrix(iat1, iat2, R_index2)->get_pointer(); + } + // if not found , skip this pair of atoms + if (tmp_c[0] != nullptr) + { + this->cal_vcomm_r_IJR(iat1, + iat2, + T0, + paraV, + nlm_tot[ad1], + nlm_tot[ad2], + tmp_c); + } + } + } + } +#ifdef _OPENMP +} +#endif + ModuleBase::timer::tick("TD_current", "calculate_vcomm_r"); +} + +// cal_HR_IJR() +void TD_current::cal_vcomm_r_IJR( + const int& iat1, + const int& iat2, + const int& T0, + const Parallel_Orbitals* paraV, + const std::vector>>>& nlm1_all, + const std::vector>>>& nlm2_all, + std::complex** data_pointer_c) +{ + // npol is the number of polarizations, + // 1 for non-magnetic (one Hamiltonian matrix only has spin-up or spin-down), + // 2 for magnetic (one Hamiltonian matrix has both spin-up and spin-down) + const int npol = this->ucell->get_npol(); + // --------------------------------------------- + // calculate the Nonlocal matrix for each pair of orbitals + // --------------------------------------------- + auto row_indexes = paraV->get_indexes_row(iat1); + auto col_indexes = paraV->get_indexes_col(iat2); + // step_trace = 0 for NSPIN=1,2; ={0, 1, local_col, local_col+1} for NSPIN=4 + std::vector step_trace(npol * npol, 0); + for (int is = 0; is < npol; is++) + { + for (int is2 = 0; is2 < npol; is2++) + { + step_trace[is * npol + is2] = col_indexes.size() * is + is2; + } + } + // calculate the local matrix + const std::complex* tmp_d = nullptr; + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + // const std::vector>* nlm1 = &(nlm1_all[0].find(row_indexes[iw1l])->second); + std::vector>*> nlm1; + for (int dir = 0; dir < 4; dir++) + { + nlm1.push_back(&(nlm1_all[dir].find(row_indexes[iw1l])->second)); + } + + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + std::vector>*> nlm2; + for (int dir = 0; dir < 4; dir++) + { + nlm2.push_back(&(nlm2_all[dir].find(col_indexes[iw2l])->second)); + } +#ifdef __DEBUG + assert(nlm1.size() == nlm2.size()); +#endif + for (int is = 0; is < npol * npol; ++is) + { + for (int dir = 0; dir < 3; dir++) + { + std::complex nlm_r_tmp = std::complex{0, 0}; + std::complex imag_unit = std::complex{0, 1}; + for (int no = 0; no < this->ucell->atoms[T0].ncpp.non_zero_count_soc[is]; no++) + { + const int p1 = this->ucell->atoms[T0].ncpp.index1_soc[is][no]; + const int p2 = this->ucell->atoms[T0].ncpp.index2_soc[is][no]; + this->ucell->atoms[T0].ncpp.get_d(is, p1, p2, tmp_d); + //- + // multiply d in the end + nlm_r_tmp += (nlm1[dir + 1]->at(p1) * std::conj(nlm2[0]->at(p2)) + - nlm1[0]->at(p1) * std::conj(nlm2[dir + 1]->at(p2))) + * (*tmp_d); + } + // -i[r,Vnl], 2.0 due to the unit transformation + data_pointer_c[dir][step_trace[is]] -= imag_unit * nlm_r_tmp / 2.0; + } + } + for (int dir = 0; dir < 3; dir++) + { + data_pointer_c[dir] += npol; + } + } + for (int dir = 0; dir < 3; dir++) + { + data_pointer_c[dir] += (npol - 1) * col_indexes.size(); + } + } +} + +void TD_current::calculate_grad_term() +{ + ModuleBase::TITLE("TD_current", "calculate_grad_term"); + if(this->current_term[0]==nullptr || this->current_term[0]->size_atom_pairs()<=0) + { + ModuleBase::WARNING_QUIT("TD_current::calculate_grad_term", "grad_term is nullptr or empty"); + } + ModuleBase::timer::tick("TD_current", "calculate_grad_term"); + + const Parallel_Orbitals* paraV = this->current_term[0]->get_atom_pair(0).get_paraV(); +#ifdef _OPENMP +#pragma omp parallel for +#endif + for (int iat1 = 0; iat1 < this->ucell->nat; iat1++) + { + auto tau1 = ucell->get_tau(iat1); + int T1, I1; + ucell->iat2iait(iat1, &I1, &T1); + AdjacentAtomInfo& adjs = this->adjs_grad[iat1]; + for (int ad = 0; ad < adjs.adj_num + 1; ++ad) + { + const int T2 = adjs.ntype[ad]; + const int I2 = adjs.natom[ad]; + const int iat2 = ucell->itia2iat(T2, I2); + const ModuleBase::Vector3& R_index2 = adjs.box[ad]; + ModuleBase::Vector3 dtau = this->ucell->cal_dtau(iat1, iat2, R_index2); + + std::complex* tmp_c[3] = {nullptr, nullptr, nullptr}; + for (int i = 0; i < 3; i++) + { + tmp_c[i] = this->current_term[i]->find_matrix(iat1, iat2, R_index2)->get_pointer(); + } + if (tmp_c[0] != nullptr) + { + this->cal_grad_IJR(iat1, iat2, paraV, dtau, tmp_c); + } + else + { + ModuleBase::WARNING_QUIT("TD_current::calculate_grad_term", "R_index not found in HR"); + } + } + } + ModuleBase::timer::tick("TD_current", "calculate_grad_term"); +} + +void TD_current::cal_grad_IJR(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const ModuleBase::Vector3& dtau, + std::complex** data_pointer_c) +{ + const LCAO_Orbitals& orb = LCAO_Orbitals::get_const_instance(); + // --------------------------------------------- + // get info of orbitals of atom1 and atom2 from ucell + // --------------------------------------------- + int T1, I1; + this->ucell->iat2iait(iat1, &I1, &T1); + int T2, I2; + this->ucell->iat2iait(iat2, &I2, &T2); + Atom& atom1 = this->ucell->atoms[T1]; + Atom& atom2 = this->ucell->atoms[T2]; + + // npol is the number of polarizations, + // 1 for non-magnetic (one Hamiltonian matrix only has spin-up or spin-down), + // 2 for magnetic (one Hamiltonian matrix has both spin-up and spin-down) + const int npol = this->ucell->get_npol(); + + const int* iw2l1 = atom1.iw2l; + const int* iw2n1 = atom1.iw2n; + const int* iw2m1 = atom1.iw2m; + const int* iw2l2 = atom2.iw2l; + const int* iw2n2 = atom2.iw2n; + const int* iw2m2 = atom2.iw2m; + // --------------------------------------------- + // get tau1 (in cell <0,0,0>) and tau2 (in cell R) + // in principle, only dtau is needed in this function + // snap_psipsi should be refactored to use dtau directly + // --------------------------------------------- + const ModuleBase::Vector3& tau1 = this->ucell->get_tau(iat1); + const ModuleBase::Vector3 tau2 = tau1 + dtau; + // --------------------------------------------- + // calculate the Ekinetic matrix for each pair of orbitals + // --------------------------------------------- + double grad[3] = {0, 0, 0}; + auto row_indexes = paraV->get_indexes_row(iat1); + auto col_indexes = paraV->get_indexes_col(iat2); + const int step_trace = col_indexes.size() + 1; + for(int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const int iw1 = row_indexes[iw1l] / npol; + const int L1 = iw2l1[iw1]; + const int N1 = iw2n1[iw1]; + const int m1 = iw2m1[iw1]; + // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) + int M1 = (m1 % 2 == 0) ? -m1 / 2 : (m1 + 1) / 2; + + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const int iw2 = col_indexes[iw2l] / npol; + const int L2 = iw2l2[iw2]; + const int N2 = iw2n2[iw2]; + const int m2 = iw2m2[iw2]; + // convert m (0,1,...2l) to M (-l, -l+1, ..., l-1, l) + int M2 = (m2 % 2 == 0) ? -m2 / 2 : (m2 + 1) / 2; + + // calculate , which equals to -. + intor_->calculate(T1, L1, N1, M1, T2, L2, N2, M2, dtau * this->ucell->lat0, nullptr, grad); + ModuleBase::Vector3 grad_overlap(grad[0], grad[1], grad[2]); + + for (int dir = 0; dir < 3; dir++) + { + for (int ipol = 0; ipol < npol; ipol++) + { + // part of Momentum operator, -iβˆ‡r,used to calculate the current + // here is actually iβˆ‡R + data_pointer_c[dir][ipol * step_trace] += std::complex(0, grad_overlap[dir]); + } + data_pointer_c[dir] += npol; + } + } + for (int dir = 0; dir < 3; dir++) + { + data_pointer_c[dir] += (npol - 1) * col_indexes.size(); + } + } +} + +#endif // __LCAO diff --git a/source/module_hamilt_lcao/module_tddft/td_current.h b/source/module_hamilt_lcao/module_tddft/td_current.h new file mode 100644 index 00000000000..a59ce93cc88 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/td_current.h @@ -0,0 +1,76 @@ +#ifndef TD_CURRENT_H +#define TD_CURRENT_H +#include +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_cell/module_neighbor/sltk_grid_driver.h" +#include "module_cell/unitcell.h" +#include "module_hamilt_lcao/module_hcontainer/hcontainer.h" +#include "module_basis/module_nao/two_center_integrator.h" +#include "td_velocity.h" +#include "module_base/vector3.h" + +#ifdef __LCAO +//design to calculate current for length gague +class TD_current +{ + public: + TD_current(const UnitCell* ucell_in, + Grid_Driver* GridD_in, + const Parallel_Orbitals* paraV, + const TwoCenterIntegrator* intor); + ~TD_current(); + + hamilt::HContainer>* get_current_term_pointer(const int& i)const + { + return this->current_term[i]; + } + + void calculate_vcomm_r(); + void calculate_grad_term(); + + private: + const UnitCell* ucell = nullptr; + + const Parallel_Orbitals* paraV = nullptr; + + Grid_Driver* Grid = nullptr; + /// @brief Store real space hamiltonian. TD term should include imaginary part, thus it has to be complex type. Only shared between TD operators. + std::vector>*> current_term = {nullptr, nullptr, nullptr}; + + const TwoCenterIntegrator* intor_ = nullptr; + + /** + * @brief initialize HR, search the nearest neighbor atoms + * HContainer is used to store the non-local pseudopotential matrix with specific atom-pairs + * the size of HR will be fixed after initialization + */ + void initialize_vcomm_r(Grid_Driver* GridD_in, const Parallel_Orbitals* paraV); + void initialize_grad_term(Grid_Driver* GridD_in, const Parallel_Orbitals* paraV); + + /** + * @brief calculate the HR local matrix of atom pair + */ + void cal_vcomm_r_IJR(const int& iat1, + const int& iat2, + const int& T0, + const Parallel_Orbitals* paraV, + const std::vector>>>& nlm1_all, + const std::vector>>>& nlm2_all, + std::complex** data_pointer_c); + void cal_grad_IJR(const int& iat1, + const int& iat2, + const Parallel_Orbitals* paraV, + const ModuleBase::Vector3& dtau, + std::complex** data_pointer_c); + + /// @brief exact the nearest neighbor atoms from all adjacent atoms + std::vector adjs_vcommr; + std::vector adjs_grad; + + /// @brief Store the vector potential for td_ekinetic term + ModuleBase::Vector3 cart_At; +}; + + +#endif // __LCAO +#endif // TD_CURRENT_H \ 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 8c4713a2f4f..e989aca91e1 100644 --- a/source/module_hamilt_lcao/module_tddft/td_velocity.cpp +++ b/source/module_hamilt_lcao/module_tddft/td_velocity.cpp @@ -60,7 +60,7 @@ void TD_Velocity::output_cart_At(const std::string& out_dir) // 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::scientific << std::setprecision(4) << std::setw(15) << cart_At[i] ; } ofs << std::endl; ofs.close(); @@ -68,10 +68,7 @@ void TD_Velocity::output_cart_At(const std::string& out_dir) return; } -void TD_Velocity::cal_cart_At(const ModuleBase::Vector3& a0, - const ModuleBase::Vector3& a1, - const ModuleBase::Vector3& a2, - const ModuleBase::Vector3& At) +void TD_Velocity::cal_cart_At(const ModuleBase::Vector3& At) { istep++; if (init_vecpot_file) @@ -80,8 +77,8 @@ void TD_Velocity::cal_cart_At(const ModuleBase::Vector3& a0, } 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]; + // transfrom into atomic unit + this->cart_At = At/2.0; } // output the vector potential if needed if (out_vecpot == true) @@ -132,8 +129,6 @@ void TD_Velocity::read_cart_At(void) } 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); } diff --git a/source/module_hamilt_lcao/module_tddft/td_velocity.h b/source/module_hamilt_lcao/module_tddft/td_velocity.h index f94c5afac5d..1a39811bb93 100644 --- a/source/module_hamilt_lcao/module_tddft/td_velocity.h +++ b/source/module_hamilt_lcao/module_tddft/td_velocity.h @@ -39,10 +39,7 @@ class TD_Velocity 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); + void cal_cart_At(const ModuleBase::Vector3& At); // allocate memory for current term. void initialize_current_term(const hamilt::HContainer>* HR, const Parallel_Orbitals* paraV); diff --git a/source/module_io/td_current_io.cpp b/source/module_io/td_current_io.cpp index d7544afb9c0..1d1e89b03fd 100644 --- a/source/module_io/td_current_io.cpp +++ b/source/module_io/td_current_io.cpp @@ -11,6 +11,7 @@ #include "module_elecstate/potentials/H_TDDFT_pw.h" #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_hamilt_lcao/module_tddft/td_velocity.h" +#include "module_hamilt_lcao/module_tddft/td_current.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #ifdef __LCAO @@ -121,23 +122,36 @@ void ModuleIO::write_current(const int istep, const K_Vectors& kv, const TwoCenterIntegrator* intor, const Parallel_Orbitals* pv, - Record_adj& ra, - hamilt::HContainer* sR) + Record_adj& ra) { ModuleBase::TITLE("ModuleIO", "write_current"); ModuleBase::timer::tick("ModuleIO", "write_current"); + TD_current* cal_current = nullptr; + std::vector>*> current_term = {nullptr, nullptr, nullptr}; if (!TD_Velocity::tddft_velocity) { - return; + cal_current = new TD_current(&GlobalC::ucell, &GlobalC::GridD, pv, intor); + cal_current->calculate_vcomm_r(); + cal_current->calculate_grad_term(); + for (int dir = 0; dir < 3; dir++) + { + current_term[dir]=cal_current->get_current_term_pointer(dir); + } } - if (TD_Velocity::td_vel_op == nullptr) + else { - ModuleBase::WARNING_QUIT("ModuleIO::write_current", "velocity gague infos is null!"); + if (TD_Velocity::td_vel_op == nullptr) + { + ModuleBase::WARNING_QUIT("ModuleIO::write_current", "velocity gague infos is null!"); + } + for (int dir = 0; dir < 3; dir++) + { + current_term[dir]=TD_Velocity::td_vel_op->get_current_term_pointer(dir); + } } - TD_Velocity* velocity_infos = TD_Velocity::td_vel_op; - ModuleBase::Vector3 cart_At = velocity_infos->cart_At; + // construct a DensityMatrix object // Since the function cal_dm_psi do not suport DMR in complex type, I replace it with two DMR in double type. Should // be refactored in the future. @@ -208,21 +222,14 @@ void ModuleIO::write_current(const int istep, // std::cout<< "iat1: " << iat1 << " iat2: " << iat2 << " Rx: " << Rx << " Ry: " << Ry << " Rz: // " << Rz << std::endl; // get BaseMatrix - hamilt::BaseMatrix* tmp_matrix_real - = DM_real.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix* tmp_matrix_imag - = DM_imag.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix* tmp_matrix_real = DM_real.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix* tmp_matrix_imag = DM_imag.get_DMR_pointer(is)->find_matrix(iat1, iat2, Rx, Ry, Rz); // refactor - hamilt::BaseMatrix>* tmp_m_rvx - = velocity_infos->get_current_term_pointer(0)->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix>* tmp_m_rvy - = velocity_infos->get_current_term_pointer(1)->find_matrix(iat1, iat2, Rx, Ry, Rz); - hamilt::BaseMatrix>* tmp_m_rvz - = velocity_infos->get_current_term_pointer(2)->find_matrix(iat1, iat2, Rx, Ry, Rz); - // get SR info - // SR are used to calculate the final part of momentum operator, eA. - hamilt::BaseMatrix* tmp_s = sR->find_matrix(iat1, iat2, Rx, Ry, Rz); - if (tmp_matrix_real == nullptr) { + hamilt::BaseMatrix>* tmp_m_rvx = current_term[0]->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix>* tmp_m_rvy = current_term[1]->find_matrix(iat1, iat2, Rx, Ry, Rz); + hamilt::BaseMatrix>* tmp_m_rvz = current_term[2]->find_matrix(iat1, iat2, Rx, Ry, Rz); + if (tmp_matrix_real == nullptr) + { continue; } int row_ap = pv->atom_begin_row[iat1]; @@ -235,7 +242,6 @@ void ModuleIO::write_current(const int istep, double dm2d1_real = 2.0 * tmp_matrix_real->get_value(mu, nu); double dm2d1_imag = 2.0 * tmp_matrix_imag->get_value(mu, nu); - double s_tmp = 0.0; std::complex rvx = {0, 0}; std::complex rvy = {0, 0}; std::complex rvz = {0, 0}; @@ -246,16 +252,10 @@ void ModuleIO::write_current(const int istep, rvy = tmp_m_rvy->get_value(mu, nu); rvz = tmp_m_rvz->get_value(mu, nu); } - if (tmp_s != nullptr) - { - s_tmp = tmp_s->get_value(mu, nu); - } - local_current_ik[0] -= dm2d1_real * rvx.real() - dm2d1_imag * rvx.imag() - + dm2d1_real * cart_At[0] * s_tmp / 2.0; - local_current_ik[1] -= dm2d1_real * rvy.real() - dm2d1_imag * rvy.imag() - + dm2d1_real * cart_At[1] * s_tmp / 2.0; - local_current_ik[2] -= dm2d1_real * rvz.real() - dm2d1_imag * rvz.imag() - + dm2d1_real * cart_At[2] * s_tmp / 2.0; + local_current_ik[0] -= dm2d1_real * rvx.real() - dm2d1_imag * rvx.imag(); + local_current_ik[1] -= dm2d1_real * rvy.real() - dm2d1_imag * rvy.imag(); + local_current_ik[2] -= dm2d1_real * rvz.real() - dm2d1_imag * rvz.imag(); + ++local_total_irr; ++irr; } // end kk @@ -304,6 +304,10 @@ void ModuleIO::write_current(const int istep, fout << istep << " " << current_total[0] << " " << current_total[1] << " " << current_total[2] << std::endl; fout.close(); } + if(!TD_Velocity::tddft_velocity) + { + delete cal_current; + } return; } #endif //__LCAO diff --git a/source/module_io/td_current_io.h b/source/module_io/td_current_io.h index f229535bdb7..17f34928d63 100644 --- a/source/module_io/td_current_io.h +++ b/source/module_io/td_current_io.h @@ -16,8 +16,7 @@ void write_current(const int istep, const K_Vectors& kv, const TwoCenterIntegrator* intor, const Parallel_Orbitals* pv, - Record_adj& ra, - hamilt::HContainer* sR); + Record_adj& ra); /// @brief calculate sum_n[𝜌_(π‘›π‘˜,πœ‡πœˆ)] for current calculation void cal_tmp_DM(elecstate::DensityMatrix, double>& DM_real, From 1c6253f5d9dad8cf159745be99af2bc3bf6099c2 Mon Sep 17 00:00:00 2001 From: esroamer Date: Thu, 4 Jul 2024 23:04:08 +0800 Subject: [PATCH 2/2] Add tests --- .../601_NO_TDDFT_H2_len_current/INPUT | 45 +++++++++++++++++++ .../integrate/601_NO_TDDFT_H2_len_current/KPT | 4 ++ .../601_NO_TDDFT_H2_len_current/STRU | 22 +++++++++ .../integrate/601_NO_TDDFT_H2_len_current/jd | 1 + .../601_NO_TDDFT_H2_len_current/result.ref | 6 +++ tests/integrate/601_NO_TDDFT_vel_H2/INPUT | 26 ++++++----- tests/integrate/601_NO_TDDFT_vel_H2/KPT | 2 +- tests/integrate/601_NO_TDDFT_vel_H2/jd | 2 +- .../integrate/601_NO_TDDFT_vel_H2/result.ref | 11 ++--- tests/integrate/601_NO_TDDFT_vel_Si/INPUT | 2 + tests/integrate/601_NO_TDDFT_vel_Si/jd | 2 +- .../integrate/601_NO_TDDFT_vel_Si/result.ref | 7 +-- tests/integrate/CASES_CPU.txt | 1 + tests/integrate/tools/catch_properties.sh | 8 ++++ 14 files changed, 116 insertions(+), 23 deletions(-) create mode 100644 tests/integrate/601_NO_TDDFT_H2_len_current/INPUT create mode 100644 tests/integrate/601_NO_TDDFT_H2_len_current/KPT create mode 100644 tests/integrate/601_NO_TDDFT_H2_len_current/STRU create mode 100644 tests/integrate/601_NO_TDDFT_H2_len_current/jd create mode 100644 tests/integrate/601_NO_TDDFT_H2_len_current/result.ref diff --git a/tests/integrate/601_NO_TDDFT_H2_len_current/INPUT b/tests/integrate/601_NO_TDDFT_H2_len_current/INPUT new file mode 100644 index 00000000000..556c7b6c025 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_len_current/INPUT @@ -0,0 +1,45 @@ +INPUT_PARAMETERS +#Parameters (General) +suffix autotest +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB +nbands 5 +calculation md +esolver_type tddft + +#Parameter (Accuracy) +ecutwfc 20 +scf_nmax 50 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +md_nstep 3 + +mixing_type broyden +mixing_beta 0.7 +mixing_gg0 0.0 +scf_thr 1.0e-6 + +cal_stress 1 +stress_thr 1e-6 +cal_force 1 +force_thr_ev 1.0e-3 + +md_type nve +md_dt 0.05 +init_vel 1 + +td_vext 1 +td_vext_dire 1 +td_stype 0 +td_ttype 0 +td_tstart 1 +td_tend 2 +td_gauss_freq 2.44 +td_gauss_phase 0.0 +td_gauss_sigma 0.5 +td_gauss_t0 0 +td_gauss_amp 0.1 + +out_current 1 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_current/KPT b/tests/integrate/601_NO_TDDFT_H2_len_current/KPT new file mode 100644 index 00000000000..c289c0158aa --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_len_current/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +1 1 1 0 0 0 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_current/STRU b/tests/integrate/601_NO_TDDFT_H2_len_current/STRU new file mode 100644 index 00000000000..af322391215 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_len_current/STRU @@ -0,0 +1,22 @@ +ATOMIC_SPECIES +H 1.00794 H.LDA.UPF + +NUMERICAL_ORBITAL +H_lda_8.0au_100Ry_2s1p.orb + +LATTICE_CONSTANT +15 + +LATTICE_VECTORS +1 0 0 #latvec1 +0 1 0 #latvec2 +0 0 1 #latvec3 + +ATOMIC_POSITIONS +Cartesian + +H #label +0 #magnetism +2 #number of atoms +0.502095122795 0.474071446432 0.509122850526 m 1 1 1 v -1.58317459417e-05 0.000196314331857 0.000681467807833 +0.497904877205 0.525928553568 0.592277149473 m 1 1 1 v 1.58317459417e-05 -0.000196314331857 -0.000681467807833 diff --git a/tests/integrate/601_NO_TDDFT_H2_len_current/jd b/tests/integrate/601_NO_TDDFT_H2_len_current/jd new file mode 100644 index 00000000000..c3541408ee7 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_len_current/jd @@ -0,0 +1 @@ +test TDDFT length-gague current output for Non-periodic system diff --git a/tests/integrate/601_NO_TDDFT_H2_len_current/result.ref b/tests/integrate/601_NO_TDDFT_H2_len_current/result.ref new file mode 100644 index 00000000000..b07509be358 --- /dev/null +++ b/tests/integrate/601_NO_TDDFT_H2_len_current/result.ref @@ -0,0 +1,6 @@ +etotref -30.91255703869586 +etotperatomref -15.4562785193 +totalforceref 0.459640 +totalstressref 0.945730 +CompareCurrent_pass 0 +totaltimeref 1.21 diff --git a/tests/integrate/601_NO_TDDFT_vel_H2/INPUT b/tests/integrate/601_NO_TDDFT_vel_H2/INPUT index 6785b062044..b83e3553ae5 100644 --- a/tests/integrate/601_NO_TDDFT_vel_H2/INPUT +++ b/tests/integrate/601_NO_TDDFT_vel_H2/INPUT @@ -14,30 +14,32 @@ scf_nmax 50 ks_solver scalapack_gvx basis_type lcao gamma_only 0 -md_nstep 2 +md_nstep 3 mixing_type broyden mixing_beta 0.7 mixing_gg0 0.0 scf_thr 1.0e-6 -cal_stress 1 -stress_thr 1e-6 -cal_force 1 -force_thr_ev 1.0e-3 +cal_stress 1 +stress_thr 1e-6 +cal_force 1 +force_thr_ev 1.0e-3 md_type nve md_dt 0.05 -init_vel 1 +init_vel 1 td_vext 1 -td_vext_dire 1 +td_vext_dire 1 td_stype 1 td_ttype 0 td_tstart 1 td_tend 2 -td_gauss_freq 2.44 -td_gauss_phase 0.0 -td_gauss_sigma 0.5 -td_gauss_t0 0.0 -td_gauss_amp 0.25 +td_gauss_freq 2.44 +td_gauss_phase 0.0 +td_gauss_sigma 0.5 +td_gauss_t0 0 +td_gauss_amp 0.1 + +out_current 1 diff --git a/tests/integrate/601_NO_TDDFT_vel_H2/KPT b/tests/integrate/601_NO_TDDFT_vel_H2/KPT index e769af76382..c289c0158aa 100644 --- a/tests/integrate/601_NO_TDDFT_vel_H2/KPT +++ b/tests/integrate/601_NO_TDDFT_vel_H2/KPT @@ -1,4 +1,4 @@ K_POINTS 0 Gamma -2 1 1 0 0 0 +1 1 1 0 0 0 diff --git a/tests/integrate/601_NO_TDDFT_vel_H2/jd b/tests/integrate/601_NO_TDDFT_vel_H2/jd index 771c8d3d17d..aac31f6ad5d 100644 --- a/tests/integrate/601_NO_TDDFT_vel_H2/jd +++ b/tests/integrate/601_NO_TDDFT_vel_H2/jd @@ -1 +1 @@ -test Non-periodic system with TDDFT velocity-gague +test TDDFT velocity-gague output for Non-periodic system diff --git a/tests/integrate/601_NO_TDDFT_vel_H2/result.ref b/tests/integrate/601_NO_TDDFT_vel_H2/result.ref index 69d8b035740..f71dd90e2a0 100644 --- a/tests/integrate/601_NO_TDDFT_vel_H2/result.ref +++ b/tests/integrate/601_NO_TDDFT_vel_H2/result.ref @@ -1,5 +1,6 @@ -etotref -30.91222553175344 -etotperatomref -15.4561127659 -totalforceref 0.574998 -totalstressref 1.137971 -totaltimeref 2.05 +etotref -30.91272078341994 +etotperatomref -15.4563603917 +totalforceref 0.459664 +totalstressref 0.945878 +CompareCurrent_pass 0 +totaltimeref 1.17 diff --git a/tests/integrate/601_NO_TDDFT_vel_Si/INPUT b/tests/integrate/601_NO_TDDFT_vel_Si/INPUT index 3914ee2b9e3..3baf2b0c624 100644 --- a/tests/integrate/601_NO_TDDFT_vel_Si/INPUT +++ b/tests/integrate/601_NO_TDDFT_vel_Si/INPUT @@ -36,3 +36,5 @@ td_gauss_phase 0.0 td_gauss_sigma 0.5 td_gauss_t0 1 td_gauss_amp 0.01 + +out_current 1 diff --git a/tests/integrate/601_NO_TDDFT_vel_Si/jd b/tests/integrate/601_NO_TDDFT_vel_Si/jd index 50c7315df5f..b2d6bb57ad9 100644 --- a/tests/integrate/601_NO_TDDFT_vel_Si/jd +++ b/tests/integrate/601_NO_TDDFT_vel_Si/jd @@ -1 +1 @@ -test periodic system with TDDFT velocity gague +test TDDFT velocity gague output for periodic system diff --git a/tests/integrate/601_NO_TDDFT_vel_Si/result.ref b/tests/integrate/601_NO_TDDFT_vel_Si/result.ref index 407d639b540..1aff915da03 100644 --- a/tests/integrate/601_NO_TDDFT_vel_Si/result.ref +++ b/tests/integrate/601_NO_TDDFT_vel_Si/result.ref @@ -1,3 +1,4 @@ -etotref -194.7716933059436 -etotperatomref -97.3858466530 -totaltimeref 9.23 +etotref -271.5879991254805 +etotperatomref -135.7939995627 +CompareCurrent_pass 0 +totaltimeref 0.99 diff --git a/tests/integrate/CASES_CPU.txt b/tests/integrate/CASES_CPU.txt index 664f7e28eea..416dbf9aec5 100644 --- a/tests/integrate/CASES_CPU.txt +++ b/tests/integrate/CASES_CPU.txt @@ -252,6 +252,7 @@ 601_NO_TDDFT_H2_len_trape 601_NO_TDDFT_H2_len_trigo 601_NO_TDDFT_O3 +601_NO_TDDFT_H2_len_current 601_NO_TDDFT_vel_H2 601_NO_TDDFT_vel_Si 701_LJ_MD_Anderson diff --git a/tests/integrate/tools/catch_properties.sh b/tests/integrate/tools/catch_properties.sh index d7f64619ccc..bd43ed9b902 100755 --- a/tests/integrate/tools/catch_properties.sh +++ b/tests/integrate/tools/catch_properties.sh @@ -68,6 +68,7 @@ out_chg=$(get_input_key_value "out_chg" "INPUT") base=$(get_input_key_value "basis_type" "INPUT") word="driver_line" symmetry=$(get_input_key_value "symmetry" "INPUT") +out_current=$(get_input_key_value "out_current" "INPUT") test -e $1 && rm $1 #-------------------------------------------- # if NOT non-self-consistent calculations @@ -457,6 +458,13 @@ if ! test -z "$symmetry" && [ $symmetry == 1 ]; then echo "nksibzref $nksibz" >>$1 fi +if ! test -z "$out_current" && [ $out_current ]; then + current1ref=refcurrent_total.dat + current1cal=OUT.autotest/current_total.dat + python3 ../tools/CompareFile.py $current1ref $current1cal 10 + echo "CompareCurrent_pass $?" >>$1 +fi + #echo $total_band ttot=`grep $word $running_path | awk '{print $3}'` echo "totaltimeref $ttot" >>$1