diff --git a/docs/advanced/input_files/input-main.md b/docs/advanced/input_files/input-main.md index c27eed5dcc..4f904f49d5 100644 --- a/docs/advanced/input_files/input-main.md +++ b/docs/advanced/input_files/input-main.md @@ -282,6 +282,7 @@ - [td\_edm](#td_edm) - [td\_print\_eij](#td_print_eij) - [td\_force\_dt](#td_force_dt) + - [propagator](#propagator) - [td\_vext](#td_vext) - [td\_vext\_dire](#td_vext_dire) - [td\_stype](#td_stype) @@ -312,9 +313,6 @@ - [out\_efield](#out_efield) - [ocp](#ocp) - [ocp\_set](#ocp_set) - - [td\_val\_elec\_01](#td_val_elec_01) - - [td\_val\_elec\_02](#td_val_elec_02) - - [td\_val\_elec\_03](#td_val_elec_03) - [Variables useful for debugging](#variables-useful-for-debugging) - [nurse](#nurse) - [t\_in\_h](#t_in_h) @@ -2337,6 +2335,16 @@ These variables are used to control berry phase and wannier90 interface paramete - **Description**: Time-dependent evolution force changes time step. (fs) - **Default**: 0.02 +### propagator + +- **Type**: Integer +- **Description**: + method of propagator. + - 0: Crank-Nicolson. + - 1: 4th Taylor expansions of exponential. + - 2: enforced time-reversal symmetry (ETRS). +- **Default**: 0 + ### td_vext - **Type**: Boolean @@ -2602,24 +2610,6 @@ These variables are used to control berry phase and wannier90 interface paramete - **Description**: If ocp is true, the ocp_set is a string to set the number of occupancy, like 1 10 * 1 0 1 representing the 13 band occupancy, 12th band occupancy 0 and the rest 1, the code is parsing this string into an array through a regular expression. - **Default**: none -### td_val_elec_01 - -- **Type**: Integer -- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the first element. -- **Default**: 1.0 - -### td_val_elec_02 - -- **Type**: Integer -- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the second element. -- **Default**: 1.0 - -### td_val_elec_03 - -- **Type**: Integer -- **Description**: Only useful when calculating the dipole. Specifies the number of valence electron associated with the third element. -- **Default**: 1.0 - [back to top](#full-list-of-input-keywords) ## Variables useful for debugging diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index a5ab2c513c..a915af5b7b 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -1,9 +1,10 @@ #include "H_TDDFT_pw.h" -#include "module_io/input.h" #include "module_base/constants.h" #include "module_base/timer.h" #include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/input.h" #include "module_io/input_conv.h" namespace elecstate @@ -16,7 +17,7 @@ int H_TDDFT_pw::istep = -1; // fuxiang add in 2017-05 //========================================================== -void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo) +void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) { ModuleBase::TITLE("H_TDDFT_pw", "cal_fixed_v"); @@ -39,18 +40,19 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo) for (auto direc: ELEC_evolve::td_vext_dire_case) { std::vector vext_space(this->rho_basis_->nrxx, 0.0); - double vext_time = cal_v_time(ttype[count]); + double vext_time = cal_v_time(ttype[count]); if (ELEC_evolve::out_efield && GlobalV::MY_RANK == 0) { std::stringstream as; - as << GlobalV::global_out_dir << "efield_"<rho_basis_->nrxx; ++ir) vl_pseudo[ir] += vext_space[ir] * vext_time; count++; @@ -94,7 +96,7 @@ void H_TDDFT_pw::read_parameters(Input *in) gauss_count = 0; gauss_omega = set_parameters(in->td_gauss_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 gauss_phase = set_parameters(in->td_gauss_phase, 1.0); - gauss_sigma = set_parameters(in->td_gauss_sigma, 1/ModuleBase::AU_to_FS); + gauss_sigma = set_parameters(in->td_gauss_sigma, 1 / ModuleBase::AU_to_FS); gauss_t0 = set_parameters(in->td_gauss_t0, 1.0); gauss_amp = set_parameters(in->td_gauss_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr @@ -163,6 +165,12 @@ void H_TDDFT_pw::cal_v_space_length(std::vector &vext_space, int direc) ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space_length"); ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space_length"); + double bmod[3]; + for (int i = 0; i < 3; i++) + { + bmod[i] = prepare(GlobalC::ucell, i); + } + for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir) { int i = ir / (this->rho_basis_->ny * this->rho_basis_->nplane); @@ -175,15 +183,15 @@ void H_TDDFT_pw::cal_v_space_length(std::vector &vext_space, int direc) switch (direc) { case 1: - vext_space[ir] = cal_v_space_length_potential(x); + vext_space[ir] = cal_v_space_length_potential(x) / bmod[0]; break; case 2: - vext_space[ir] = cal_v_space_length_potential(y); + vext_space[ir] = cal_v_space_length_potential(y) / bmod[1]; break; case 3: - vext_space[ir] = cal_v_space_length_potential(z); + vext_space[ir] = cal_v_space_length_potential(z) / bmod[2]; break; default: @@ -198,18 +206,18 @@ void H_TDDFT_pw::cal_v_space_length(std::vector &vext_space, int direc) double H_TDDFT_pw::cal_v_space_length_potential(double i) { - double vext_space=0.0; - if (i < this->rho_basis_->nx * lcut1) + double vext_space = 0.0; + if (i < lcut1) { - vext_space = ((i / this->rho_basis_->nx - lcut1)*(lcut2-lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0; + vext_space = ((i - lcut1) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0; } - else if (i >= this->rho_basis_->nx * lcut1 && i < this->rho_basis_->nx * lcut2) + else if (i >= lcut1 && i < lcut2) { - vext_space = -i / this->rho_basis_->nx * this->ucell_->lat0; + vext_space = -i * this->ucell_->lat0; } - else if (i >= this->rho_basis_->nx * lcut2) + else if (i >= lcut2) { - vext_space = ((i / this->rho_basis_->nx - lcut2)*(lcut2-lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0; + vext_space = ((i - lcut2) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0; } return vext_space; } @@ -229,6 +237,10 @@ double H_TDDFT_pw::cal_v_time(int t_type) vext_time = cal_v_time_Gauss(); break; + case 1: + vext_time = cal_v_time_trapezoid(); + break; + case 2: vext_time = cal_v_time_trigonometric(); break; @@ -237,9 +249,9 @@ double H_TDDFT_pw::cal_v_time(int t_type) vext_time = cal_v_time_heaviside(); break; - // case 4: - // vext_time = cal_v_time_HHG(); - // break; + // case 4: + // vext_time = cal_v_time_HHG(); + // break; default: std::cout << "time_domain_type of electric field is wrong" << endl; @@ -304,8 +316,7 @@ double H_TDDFT_pw::cal_v_time_trigonometric() const double timenow = istep * dt; - vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) - * sin(omega2 * timenow + phase2); + vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) * sin(omega2 * timenow + phase2); trigo_count++; return vext_time; @@ -345,4 +356,34 @@ double H_TDDFT_pw::cal_v_time_heaviside() // return vext_time; // } +double H_TDDFT_pw::prepare(const UnitCell &cell, int &dir) +{ + double bvec[3] = {0.0}; + double bmod = 0.0; + if (dir == 0) + { + bvec[0] = cell.G.e11; + bvec[1] = cell.G.e12; + bvec[2] = cell.G.e13; + } + else if (dir == 1) + { + bvec[0] = cell.G.e21; + bvec[1] = cell.G.e22; + bvec[2] = cell.G.e23; + } + else if (dir == 2) + { + bvec[0] = cell.G.e31; + bvec[1] = cell.G.e32; + bvec[2] = cell.G.e33; + } + else + { + ModuleBase::WARNING_QUIT("H_TDDFT_pw::prepare", "direction is wrong!"); + } + bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); + return bmod; +} + } // namespace elecstate diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h index df9ea2b3d3..e600219b24 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.h +++ b/source/module_elecstate/potentials/H_TDDFT_pw.h @@ -109,6 +109,8 @@ class H_TDDFT_pw : public PotBase double cal_v_time_trigonometric(); double cal_v_time_heaviside(); // double cal_v_time_HHG(); + + double prepare(const UnitCell &cell, int &dir); }; } // namespace elecstate diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index f766d4d9db..82a0486d88 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -1,11 +1,11 @@ #include "esolver_ks_lcao_tddft.h" #include "module_io/cal_r_overlap_R.h" +#include "module_io/dipole_io.h" #include "module_io/dm_io.h" #include "module_io/rho_io.h" -#include "module_io/dipole_io.h" -#include "module_io/write_HS_R.h" #include "module_io/write_HS.h" +#include "module_io/write_HS_R.h" //--------------temporary---------------------------- #include "module_base/blas_connector.h" @@ -41,6 +41,14 @@ ESolver_KS_LCAO_TDDFT::~ESolver_KS_LCAO_TDDFT() { // this->orb_con.clear_after_ions(GlobalC::UOT, GlobalC::ORB, GlobalV::deepks_setorb, GlobalC::ucell.infoNL.nproj); delete psi_laststep; + if (Hk_laststep != nullptr) + { + for (int ik = 0; ik < GlobalC::kv.nks; ++ik) + { + delete Hk_laststep[ik]; + } + delete Hk_laststep; + } } void ESolver_KS_LCAO_TDDFT::Init(Input& inp, UnitCell& ucell) @@ -188,7 +196,15 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr) if (GlobalV::ESOLVER_TYPE == "tddft" && istep >= 2 && !GlobalV::GAMMA_ONLY_LOCAL) { - ELEC_evolve::evolve_psi(istep, this->p_hamilt, this->LOWF, this->psi, this->psi_laststep, this->pelec_td->ekb); + ELEC_evolve::evolve_psi(istep, + this->p_hamilt, + this->LOWF, + this->psi, + this->psi_laststep, + this->Hk_laststep, + this->pelec_td->ekb, + td_htype, + INPUT.propagator); this->pelec_td->psiToRho_td(this->psi[0]); // this->pelec_td->psiToRho(this->psi[0]); } @@ -274,7 +290,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr) void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) { - //print Hamiltonian and Overlap matrix + // print Hamiltonian and Overlap matrix if (this->conv_elec) { if (!GlobalV::GAMMA_ONLY_LOCAL) @@ -283,7 +299,7 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) } for (int ik = 0; ik < GlobalC::kv.nks; ++ik) { - if(hsolver::HSolverLCAO::out_mat_hs) + if (hsolver::HSolverLCAO::out_mat_hs) { this->p_hamilt->updateHk(ik); } @@ -342,7 +358,8 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) if (!this->conv_elec) { - if(GlobalV::NSPIN==4) GlobalC::ucell.cal_ux(); + if (GlobalV::NSPIN == 4) + GlobalC::ucell.cal_ux(); this->pelec->pot->update_from_charge(this->pelec->charge, &GlobalC::ucell); GlobalC::en.delta_escf(this->pelec); } @@ -351,7 +368,7 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) GlobalC::en.cal_converged(this->pelec); } - // store wfc + // store wfc and Hk laststep if (istep >= 1 && this->conv_elec) { if (this->psi_laststep == nullptr) @@ -365,12 +382,37 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) = new psi::Psi>(GlobalC::kv.nks, GlobalV::NBANDS, GlobalV::NLOCAL, nullptr); #endif - std::complex *p_psi = &psi[0](0,0,0); - std::complex *p_psi_laststep = &psi_laststep[0](0,0,0); - for (int index = 0; index < psi[0].size(); ++index) + if (td_htype == 1) { - p_psi_laststep[index] = p_psi[index]; + if (this->Hk_laststep == nullptr) + { + this->Hk_laststep = new std::complex*[GlobalC::kv.nks]; + for (int ik = 0; ik < GlobalC::kv.nks; ++ik) + { + this->Hk_laststep[ik] = new std::complex[this->LOC.ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Hk_laststep[ik], this->LOC.ParaV->nloc); + } + } } + + for (int ik = 0; ik < GlobalC::kv.nks; ++ik) + { + this->psi->fix_k(ik); + this->psi_laststep->fix_k(ik); + int size0 = psi->get_nbands() * psi->get_nbasis(); + for (int index = 0; index < size0; ++index) + psi_laststep[0].get_pointer()[index] = psi[0].get_pointer()[index]; + + // store Hamiltonian + if (td_htype == 1) + { + this->p_hamilt->updateHk(ik); + hamilt::MatrixBlock> h_mat, s_mat; + this->p_hamilt->matrix(h_mat, s_mat); + BlasConnector::copy(this->LOC.ParaV->nloc, h_mat.p, 1, Hk_laststep[ik], 1); + } + } + if (istep > 1 && ELEC_evolve::td_edm == 0) this->cal_edm_tddft(); } @@ -413,7 +455,7 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep) const int precision = 3; std::stringstream ssc; ssc << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG.cube"; - double& ef_tmp = GlobalC::en.get_ef(is,GlobalV::TWO_EFERMI); + double& ef_tmp = GlobalC::en.get_ef(is, GlobalV::TWO_EFERMI); ModuleIO::write_rho( #ifdef __MPI GlobalC::bigpw->bz, @@ -527,7 +569,6 @@ void ESolver_KS_LCAO_TDDFT::afterscf(const int istep) else { r_matrix.out_rR(istep); - } } } @@ -547,14 +588,12 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() complex* tmp2 = new complex[this->LOC.ParaV->nloc]; complex* tmp3 = new complex[this->LOC.ParaV->nloc]; complex* tmp4 = new complex[this->LOC.ParaV->nloc]; - complex* tmp5 = new complex[this->LOC.ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Htmp, this->LOC.ParaV->nloc); ModuleBase::GlobalFunc::ZEROS(Sinv, this->LOC.ParaV->nloc); ModuleBase::GlobalFunc::ZEROS(tmp1, this->LOC.ParaV->nloc); ModuleBase::GlobalFunc::ZEROS(tmp2, this->LOC.ParaV->nloc); ModuleBase::GlobalFunc::ZEROS(tmp3, this->LOC.ParaV->nloc); ModuleBase::GlobalFunc::ZEROS(tmp4, this->LOC.ParaV->nloc); - ModuleBase::GlobalFunc::ZEROS(tmp5, this->LOC.ParaV->nloc); const int inc = 1; int nrow = this->LOC.ParaV->nrow; int ncol = this->LOC.ParaV->ncol; @@ -604,8 +643,8 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() const char N_char = 'N', T_char = 'T'; const complex one_float = {1.0, 0.0}, zero_float = {0.0, 0.0}; const complex half_float = {0.5, 0.0}; - pzgemm_(&T_char, - &T_char, + pzgemm_(&N_char, + &N_char, &GlobalV::NLOCAL, &GlobalV::NLOCAL, &GlobalV::NLOCAL, @@ -645,7 +684,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() this->LOC.ParaV->desc); pzgemm_(&N_char, - &T_char, + &N_char, &GlobalV::NLOCAL, &GlobalV::NLOCAL, &GlobalV::NLOCAL, @@ -665,7 +704,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() this->LOC.ParaV->desc); pzgemm_(&N_char, - &T_char, + &N_char, &GlobalV::NLOCAL, &GlobalV::NLOCAL, &GlobalV::NLOCAL, @@ -698,19 +737,7 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() &one_int, this->LOC.ParaV->desc); - pztranu_(&GlobalV::NLOCAL, - &GlobalV::NLOCAL, - &one_float, - tmp4, - &one_int, - &one_int, - this->LOC.ParaV->desc, - &zero_float, - tmp5, - &one_int, - &one_int, - this->LOC.ParaV->desc); - zcopy_(&this->LOC.ParaV->nloc, tmp5, &inc, this->LOC.edm_k_tddft[ik].c, &inc); + zcopy_(&this->LOC.ParaV->nloc, tmp4, &inc, this->LOC.edm_k_tddft[ik].c, &inc); delete[] Htmp; delete[] Sinv; @@ -718,7 +745,6 @@ void ESolver_KS_LCAO_TDDFT::cal_edm_tddft() delete[] tmp2; delete[] tmp3; delete[] tmp4; - delete[] tmp5; delete[] ipiv; #else this->LOC.edm_k_tddft[ik].create(this->LOC.ParaV->ncol, this->LOC.ParaV->nrow); diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index a943cebe5a..6df05bb38e 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -22,8 +22,10 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO void Init(Input& inp, UnitCell& cell) override; psi::Psi>* psi_laststep = nullptr; + std::complex** Hk_laststep = nullptr; //same as pelec elecstate::ElecStateLCAO_TDDFT* pelec_td = nullptr; + int td_htype = 1; protected: virtual void hamilt2density(const int istep, const int iter, const double ethr) override; diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp index ade328fce4..4b37b54197 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.cpp @@ -1,20 +1,17 @@ #include "ELEC_evolve.h" -#include "module_base/timer.h" +#include "LCAO_evolve.h" #include "module_base/parallel_reduce.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_base/timer.h" #include "module_elecstate/module_charge/symmetry_rho.h" -#include "LCAO_evolve.h" -#include "module_hamilt_lcao/module_dftu/dftu.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" +#include "module_hamilt_lcao/module_dftu/dftu.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" ELEC_evolve::ELEC_evolve(){}; ELEC_evolve::~ELEC_evolve(){}; double ELEC_evolve::td_force_dt; -int ELEC_evolve::td_val_elec_01; -int ELEC_evolve::td_val_elec_02; -int ELEC_evolve::td_val_elec_03; bool ELEC_evolve::td_vext; std::vector ELEC_evolve::td_vext_dire_case; bool ELEC_evolve::out_dipole; @@ -28,7 +25,10 @@ void ELEC_evolve::evolve_psi(const int& istep, Local_Orbital_wfc& lowf, psi::Psi>* psi, psi::Psi>* psi_laststep, - ModuleBase::matrix& ekb) + std::complex** Hk_laststep, + ModuleBase::matrix& ekb, + int htype, + int propagator) { ModuleBase::TITLE("ELEC_evolve", "eveolve_psi"); ModuleBase::timer::tick("ELEC_evolve", "evolve_psi"); @@ -42,7 +42,14 @@ void ELEC_evolve::evolve_psi(const int& istep, Evolve_LCAO_Matrix ELM(lowf.ParaV); psi->fix_k(ik); psi_laststep->fix_k(ik); - ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, &(ekb(ik, 0))); + if (Hk_laststep == nullptr) + { + ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, nullptr, &(ekb(ik, 0)), htype, propagator); + } + else + { + ELM.evolve_complex_matrix(ik, phm, psi, psi_laststep, Hk_laststep[ik], &(ekb(ik, 0)), htype, propagator); + } ModuleBase::timer::tick("Efficience", "evolve_k"); } // end k diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h index 86277853c6..dfa2b904eb 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h @@ -30,9 +30,6 @@ class ELEC_evolve // fuxiang add 2021-05-25 static double td_force_dt; - static int td_val_elec_01; - static int td_val_elec_02; - static int td_val_elec_03; static bool td_vext; static std::vector td_vext_dire_case; static bool out_dipole; @@ -47,7 +44,10 @@ class ELEC_evolve Local_Orbital_wfc& lowf, psi::Psi>* psi, psi::Psi>* psi_laststep, - ModuleBase::matrix& ekb); + std::complex** Hk_laststep, + ModuleBase::matrix& ekb, + int htype, + int propagator); }; #endif diff --git a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp index 66325ee7cf..cd7649f536 100644 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp +++ b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp @@ -1,12 +1,12 @@ #include "LCAO_evolve.h" -#include "module_io/input.h" -#include "module_base/scalapack_connector.h" -#include "module_base/lapack_connector.h" -#include "module_io/write_HS.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" #include "ELEC_evolve.h" +#include "module_base/lapack_connector.h" +#include "module_base/scalapack_connector.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" +#include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_io/input.h" +#include "module_io/write_HS.h" #include // fuxiang add 2016-10-28 @@ -27,7 +27,10 @@ void Evolve_LCAO_Matrix::evolve_complex_matrix(const int& ik, hamilt::Hamilt* p_hamilt, psi::Psi>* psi_k, psi::Psi>* psi_k_laststep, - double* ekb) const + std::complex* H_laststep, + double* ekb, + int htype, + int propagator) const { ModuleBase::TITLE("Evolve_LCAO_Matrix", "evolve_complex_matrix"); time_t time_start = time(NULL); @@ -36,7 +39,15 @@ void Evolve_LCAO_Matrix::evolve_complex_matrix(const int& ik, if (GlobalV::ESOLVER_TYPE == "tddft") { #ifdef __MPI - this->using_ScaLAPACK_complex(GlobalV::NBANDS, GlobalV::NLOCAL, psi_k_laststep[0].get_pointer(), p_hamilt, psi_k[0].get_pointer(), ekb); + this->using_ScaLAPACK_complex(GlobalV::NBANDS, + GlobalV::NLOCAL, + psi_k_laststep[0].get_pointer(), + H_laststep, + p_hamilt, + psi_k[0].get_pointer(), + ekb, + htype, + propagator); #else this->using_LAPACK_complex(ik, p_hamilt, psi_k[0].get_pointer(), psi_k_laststep[0].get_pointer(), ekb); #endif @@ -340,13 +351,15 @@ void Evolve_LCAO_Matrix::using_LAPACK_complex(const int& ik, #ifdef __MPI -void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( - const int nband, - const int nlocal, - const std::complex* psi_k_laststep, - hamilt::Hamilt* p_hamilt, - std::complex* psi_k, - double* ekb) const +void Evolve_LCAO_Matrix::using_ScaLAPACK_complex(const int nband, + const int nlocal, + const std::complex* psi_k_laststep, + const std::complex* H_laststep, + hamilt::Hamilt* p_hamilt, + std::complex* psi_k, + double* ekb, + int htype, + int propagator) const { ModuleBase::TITLE("Evolve_LCAO_Matrix", "using_ScaLAPACK_complex"); @@ -361,40 +374,69 @@ void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( complex* Htmp = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Htmp, this->ParaV->nloc); BlasConnector::copy(this->ParaV->nloc, h_mat.p, 1, Htmp, 1); - + + complex* Hold = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Hold, this->ParaV->nloc); + BlasConnector::copy(this->ParaV->nloc, h_mat.p, 1, Hold, 1); + complex* U_operator = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(U_operator, this->ParaV->nloc); -// (1)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // (1)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + /// @brief compute H(t+dt/2) + /// @input H_laststep, Htmp, print_matrix + /// @output Htmp + if (htype == 1 && propagator!=2) + { + half_Hmatrix(nband, nlocal, Htmp, H_laststep, print_matrix); + } + + // (2)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> /// @brief compute U_operator /// @input Stmp, Htmp, print_matrix /// @output U_operator - compute_U_operator(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + switch (propagator) + { + case 0: + compute_U_operator_CN2(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + break; + + case 1: + compute_U_operator_taylor(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + break; + + case 2: + compute_U_operator_etrs(nband, nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); -// (2)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + break; + + default: + std::cout << "method of propagator is wrong" << endl; + break; + } + + // (3)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> /// @brief apply U_operator to the wave function of the previous step for new wave function - /// @input U_operator, psi_k_laststep + /// @input U_operator, psi_k_laststep, print_matrix /// @output psi_k - U_to_wfc(nband, nlocal, U_operator, psi_k_laststep, psi_k); + U_to_wfc(nband, nlocal, U_operator, psi_k_laststep, psi_k, print_matrix); - - -// (3)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // (4)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> /// @brief normalize psi_k /// @input Stmp, psi_not_norm, psi_k, print_matrix /// @output psi_k norm_wfc(nband, nlocal, Stmp, psi_k, print_matrix); + // (5)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -// (4)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - /// @brief compute ekb /// @input Htmp, psi_k /// @output ekb - compute_ekb(nband, nlocal, Htmp, psi_k, ekb); + compute_ekb(nband, nlocal, Hold, psi_k, ekb); delete[] Stmp; delete[] Htmp; @@ -402,14 +444,80 @@ void Evolve_LCAO_Matrix::using_ScaLAPACK_complex( return; } +void Evolve_LCAO_Matrix::half_Hmatrix(const int nband, + const int nlocal, + std::complex* Htmp, + const std::complex* H_laststep, + const int print_matrix) const +{ + if (print_matrix) + { + GlobalV::ofs_running << std::setprecision(10); + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H(t+dt) :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" + << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H(t):" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << H_laststep[i * this->ParaV->ncol + j].real() << "+" + << H_laststep[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } + + complex alpha = {0.5, 0.0}; + complex beta = {0.5, 0.0}; + ScalapackConnector::geadd('N', + nlocal, + nlocal, + alpha, + H_laststep, + 1, + 1, + this->ParaV->desc, + beta, + Htmp, + 1, + 1, + this->ParaV->desc); + + if (print_matrix) + { + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H (t+dt/2) :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" + << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } +} -void Evolve_LCAO_Matrix::compute_U_operator( - const int nband, - const int nlocal, - const std::complex* Stmp, - const std::complex* Htmp, - std::complex* U_operator, - const int print_matrix) const +void Evolve_LCAO_Matrix::compute_U_operator_CN2(const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const { // (1) copy Htmp to Numerator & Denominator complex* Numerator = new complex[this->ParaV->nloc]; @@ -428,7 +536,8 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Stmp[i * this->ParaV->ncol + j].real() << "+" << Stmp[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Stmp[i * this->ParaV->ncol + j].real() << "+" + << Stmp[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -439,14 +548,15 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" << Numerator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" + << Numerator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } GlobalV::ofs_running << endl; } -// ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // (2) compute Numerator & Denominator by GEADD // Numerator = Stmp - i*para * Htmp; beta1 = - para = -0.25 * INPUT.mdp.md_dt // Denominator = Stmp + i*para * Htmp; beta2 = para = 0.25 * INPUT.mdp.md_dt @@ -454,128 +564,111 @@ void Evolve_LCAO_Matrix::compute_U_operator( complex beta1 = {0.0, -0.25 * INPUT.mdp.md_dt}; complex beta2 = {0.0, 0.25 * INPUT.mdp.md_dt}; - ScalapackConnector::geadd( - 'N', - nlocal, - nlocal, - alpha, - Stmp, - 1, - 1, - this->ParaV->desc, - beta1, - Numerator, - 1, - 1, - this->ParaV->desc - ); - ScalapackConnector::geadd( - 'N', - nlocal, - nlocal, - alpha, - Stmp, - 1, - 1, - this->ParaV->desc, - beta2, - Denominator, - 1, - 1, - this->ParaV->desc - ); + ScalapackConnector::geadd('N', + nlocal, + nlocal, + alpha, + Stmp, + 1, + 1, + this->ParaV->desc, + beta1, + Numerator, + 1, + 1, + this->ParaV->desc); + ScalapackConnector::geadd('N', + nlocal, + nlocal, + alpha, + Stmp, + 1, + 1, + this->ParaV->desc, + beta2, + Denominator, + 1, + 1, + this->ParaV->desc); if (print_matrix) { + GlobalV::ofs_running << " beta=" << beta1 << endl; GlobalV::ofs_running << " fenmu:" << endl; for (int i = 0; i < this->ParaV->nrow; i++) { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" << Denominator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" + << Denominator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } GlobalV::ofs_running << endl; } -//->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // (3) Next, invert Denominator int* ipiv = new int[this->ParaV->nloc]; int info = 0; // (3.1) compute ipiv - ScalapackConnector::getrf( - nlocal, - nlocal, - Denominator, - 1, - 1, - this->ParaV->desc, - ipiv, - &info - ); + ScalapackConnector::getrf(nlocal, nlocal, Denominator, 1, 1, this->ParaV->desc, ipiv, &info); int lwork = -1; int liwotk = -1; std::vector> work(1, 0); std::vector iwork(1, 0); // (3.2) compute work - ScalapackConnector::getri( - nlocal, - Denominator, - 1, - 1, - this->ParaV->desc, - ipiv, - work.data(), - &lwork, - iwork.data(), - &liwotk, - &info - ); + ScalapackConnector::getri(nlocal, + Denominator, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); lwork = work[0].real(); work.resize(lwork, 0); liwotk = iwork[0]; iwork.resize(liwotk, 0); - // (3.3) compute inverse matrix of matrix_A - ScalapackConnector::getri( - nlocal, - Denominator, - 1, - 1, - this->ParaV->desc, - ipiv, - work.data(), - &lwork, - iwork.data(), - &liwotk, - &info - ); + // (3.3) compute inverse matrix of Denominator + ScalapackConnector::getri(nlocal, + Denominator, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); assert(0 == info); -//->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> // (4) U_operator = Denominator * Numerator; - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nlocal, - nlocal, - 1.0, - Denominator, - 1, - 1, - this->ParaV->desc, - Numerator, - 1, - 1, - this->ParaV->desc, - 0.0, - U_operator, - 1, - 1, - this->ParaV->desc - ); + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + Denominator, + 1, + 1, + this->ParaV->desc, + Numerator, + 1, + 1, + this->ParaV->desc, + 0.0, + U_operator, + 1, + 1, + this->ParaV->desc); if (print_matrix) { @@ -584,7 +677,8 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" << Denominator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Denominator[i * this->ParaV->ncol + j].real() << "+" + << Denominator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -594,7 +688,8 @@ void Evolve_LCAO_Matrix::compute_U_operator( { for (int j = 0; j < this->ParaV->ncol; j++) { - GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" << Numerator[i * this->ParaV->ncol + j].imag() << "i "; + GlobalV::ofs_running << Numerator[i * this->ParaV->ncol + j].real() << "+" + << Numerator[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -617,103 +712,513 @@ void Evolve_LCAO_Matrix::compute_U_operator( } } -// cout << "U_operator Success!!!" <* U_operator, - const std::complex* psi_k_laststep, - std::complex* psi_k) const +void Evolve_LCAO_Matrix::compute_U_operator_taylor(const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const { + ModuleBase::GlobalFunc::ZEROS(U_operator, this->ParaV->nloc); + complex* A_matrix = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(A_matrix, this->ParaV->nloc); + complex* rank0 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank0, this->ParaV->nloc); + complex* rank2 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank2, this->ParaV->nloc); + complex* rank3 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank3, this->ParaV->nloc); + complex* rank4 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank4, this->ParaV->nloc); + complex* tmp1 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(tmp1, this->ParaV->nloc); + complex* tmp2 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(tmp2, this->ParaV->nloc); + complex* Sinv = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Sinv, this->ParaV->nloc); + BlasConnector::copy(this->ParaV->nloc, Stmp, 1, Sinv, 1); + + if (print_matrix) + { + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " S matrix :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Stmp[i * this->ParaV->ncol + j].real() << "+" + << Stmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " H matrix :" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * this->ParaV->ncol + j].real() << "+" + << Htmp[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } + + // set rank0 + int info; + int myid; + MPI_Comm_rank(this->ParaV->comm_2D, &myid); + int naroc[2]; // maximum number of row or column - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nlocal, - 1.0, - U_operator, - 1, - 1, - this->ParaV->desc, - psi_k_laststep, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc - ); + for (int iprow = 0; iprow < this->ParaV->dim0; ++iprow) + { + for (int ipcol = 0; ipcol < this->ParaV->dim1; ++ipcol) + { + const int coord[2] = {iprow, ipcol}; + int src_rank; + info = MPI_Cart_rank(this->ParaV->comm_2D, coord, &src_rank); + if (myid == src_rank) + { + naroc[0] = this->ParaV->nrow; + naroc[1] = this->ParaV->ncol; + for (int j = 0; j < naroc[1]; ++j) + { + int igcol = globalIndex(j, this->ParaV->nb, this->ParaV->dim1, ipcol); + if (igcol >= nlocal) + continue; + for (int i = 0; i < naroc[0]; ++i) + { + int igrow = globalIndex(i, this->ParaV->nb, this->ParaV->dim0, iprow); + if (igrow >= nlocal) + continue; + if (igcol == igrow) + { + rank0[j * naroc[0] + i] = {1.0, 0.0}; + } + else + { + rank0[j * naroc[0] + i] = {0.0, 0.0}; + } + } + } + } + } // loop ipcol + } // loop iprow + + // complex beta = {0.0, -0.5 * INPUT.mdp.md_dt}; + complex beta = {0.0, -0.25 * INPUT.mdp.md_dt}; // for ETRS + + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // invert Stmp + int* ipiv = new int[this->ParaV->nloc]; + // (3.1) compute ipiv + ScalapackConnector::getrf(nlocal, nlocal, Sinv, 1, 1, this->ParaV->desc, ipiv, &info); + int lwork = -1; + int liwotk = -1; + std::vector> work(1, 0); + std::vector iwork(1, 0); + // (3.2) compute work + ScalapackConnector::getri(nlocal, + Sinv, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); + lwork = work[0].real(); + work.resize(lwork, 0); + liwotk = iwork[0]; + iwork.resize(liwotk, 0); + ScalapackConnector::getri(nlocal, + Sinv, + 1, + 1, + this->ParaV->desc, + ipiv, + work.data(), + &lwork, + iwork.data(), + &liwotk, + &info); + assert(0 == info); + + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + // A_matrix = - idt S^-1 H ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + beta, + Sinv, + 1, + 1, + this->ParaV->desc, + Htmp, + 1, + 1, + this->ParaV->desc, + 0.0, + U_operator, + 1, + 1, + this->ParaV->desc); + + // rank2 = A^2 ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + U_operator, + 1, + 1, + this->ParaV->desc, + 0.0, + rank2, + 1, + 1, + this->ParaV->desc); + + // rank3 = A^3 ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + rank2, + 1, + 1, + this->ParaV->desc, + 0.0, + rank3, + 1, + 1, + this->ParaV->desc); + + // rank4 = A^4 ; + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + rank3, + 1, + 1, + this->ParaV->desc, + 0.0, + rank4, + 1, + 1, + this->ParaV->desc); + + complex p1 = {1.0, 0.0}; + complex p2 = {1.0 / 2.0, 0.0}; + complex p3 = {1.0 / 6.0, 0.0}; + complex p4 = {1.0 / 24.0, 0.0}; + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p1, + rank0, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p2, + rank2, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p3, + rank3, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + ScalapackConnector::geadd('N', + nlocal, + nlocal, + p4, + rank4, + 1, + 1, + this->ParaV->desc, + p1, + U_operator, + 1, + 1, + this->ParaV->desc); + + if (print_matrix) + { + GlobalV::ofs_running << " A_matrix:" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + GlobalV::ofs_running << A_matrix[i * this->ParaV->ncol + j].real() << "+" + << A_matrix[i * this->ParaV->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " U operator:" << endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + double aa, bb; + aa = U_operator[i * this->ParaV->ncol + j].real(); + bb = U_operator[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + } + + // cout << "U_operator Success!!!" <* Stmp, + const std::complex* Htmp, + const std::complex* H_laststep, + std::complex* U_operator, + const int print_matrix) const +{ + complex* U1 = new complex[this->ParaV->nloc]; + complex* U2 = new complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(U1, this->ParaV->nloc); + ModuleBase::GlobalFunc::ZEROS(U2, this->ParaV->nloc); + compute_U_operator_taylor(nband, nlocal, Stmp, Htmp, U1, print_matrix); + compute_U_operator_taylor(nband, nlocal, Stmp, H_laststep, U2, print_matrix); + ScalapackConnector::gemm('N', + 'N', + nlocal, + nlocal, + nlocal, + 1.0, + U1, + 1, + 1, + this->ParaV->desc, + U2, + 1, + 1, + this->ParaV->desc, + 0.0, + U_operator, + 1, + 1, + this->ParaV->desc); +} -void Evolve_LCAO_Matrix::norm_wfc( - const int nband, - const int nlocal, - const std::complex* Stmp, - std::complex* psi_k, - const int print_matrix) const +void Evolve_LCAO_Matrix::U_to_wfc(const int nband, + const int nlocal, + const std::complex* U_operator, + const std::complex* psi_k_laststep, + std::complex* psi_k, + const int print_matrix) const +{ + + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + U_operator, + 1, + 1, + this->ParaV->desc, + psi_k_laststep, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc); + + if (print_matrix) + { + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " psi_k:" << endl; + for (int i = 0; i < this->ParaV->ncol_bands; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + double aa, bb; + aa = psi_k[i * this->ParaV->ncol + j].real(); + bb = psi_k[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + GlobalV::ofs_running << " psi_k_laststep:" << endl; + for (int i = 0; i < this->ParaV->ncol_bands; i++) + { + for (int j = 0; j < this->ParaV->ncol; j++) + { + double aa, bb; + aa = psi_k_laststep[i * this->ParaV->ncol + j].real(); + bb = psi_k_laststep[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } +} + +void Evolve_LCAO_Matrix::norm_wfc(const int nband, + const int nlocal, + const std::complex* Stmp, + std::complex* psi_k, + const int print_matrix) const { complex* tmp1 = new complex[this->ParaV->nloc_wfc]; ModuleBase::GlobalFunc::ZEROS(tmp1, this->ParaV->nloc_wfc); - + complex* Cij = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Cij, this->ParaV->nloc); - - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nlocal, - 1.0, - Stmp, - 1, - 1, - this->ParaV->desc, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - tmp1, - 1, - 1, - this->ParaV->desc_wfc - ); - - ScalapackConnector::gemm( - 'C', - 'N', - nband, - nband, - nlocal, - 1.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - tmp1, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - Cij, - 1, - 1, - this->ParaV->desc - ); + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + Stmp, + 1, + 1, + this->ParaV->desc, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + tmp1, + 1, + 1, + this->ParaV->desc_wfc); + + ScalapackConnector::gemm('C', + 'N', + nband, + nband, + nlocal, + 1.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + tmp1, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + Cij, + 1, + 1, + this->ParaV->desc_Eij); + + if (print_matrix) + { + GlobalV::ofs_running << "original Cij :" << endl; + for (int i = 0; i < this->ParaV->ncol; i++) + { + for (int j = 0; j < this->ParaV->nrow; j++) + { + double aa, bb; + aa = Cij[i * this->ParaV->ncol + j].real(); + bb = Cij[i * this->ParaV->ncol + j].imag(); + if (abs(aa) < 1e-8) + aa = 0.0; + if (abs(bb) < 1e-8) + bb = 0.0; + GlobalV::ofs_running << aa << "+" << bb << "i "; + } + GlobalV::ofs_running << endl; + } + GlobalV::ofs_running << endl; + } int info; int myid; @@ -755,32 +1260,27 @@ void Evolve_LCAO_Matrix::norm_wfc( } // loop ipcol } // loop iprow - // std::cout << "nlocal" << nlocal << std::endl; - // std::cout << "GlobalV::NLOCAL" << GlobalV::NLOCAL << std::endl; BlasConnector::copy(this->ParaV->nloc_wfc, psi_k, 1, tmp1, 1); - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nband, - 1.0, - tmp1, - 1, - 1, - this->ParaV->desc_wfc, - Cij, - 1, - 1, - this->ParaV->desc, - 0.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc - ); - + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nband, + 1.0, + tmp1, + 1, + 1, + this->ParaV->desc_wfc, + Cij, + 1, + 1, + this->ParaV->desc_Eij, + 0.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc); if (print_matrix) { @@ -789,7 +1289,8 @@ void Evolve_LCAO_Matrix::norm_wfc( { for (int j = 0; j < this->ParaV->nrow; j++) { - GlobalV::ofs_running << Cij[i * this->ParaV->ncol + j].real() << "+" << Cij[i * this->ParaV->ncol_bands + j].imag() << "i "; + GlobalV::ofs_running << Cij[i * this->ParaV->ncol + j].real() << "+" + << Cij[i * this->ParaV->ncol + j].imag() << "i "; } GlobalV::ofs_running << endl; } @@ -812,7 +1313,7 @@ void Evolve_LCAO_Matrix::norm_wfc( GlobalV::ofs_running << endl; } GlobalV::ofs_running << endl; - GlobalV::ofs_running << " psi_k nlocal*nlocal:" << endl; + GlobalV::ofs_running << " psi_k before normalization:" << endl; for (int i = 0; i < this->ParaV->ncol; i++) { for (int j = 0; j < this->ParaV->ncol; j++) @@ -832,19 +1333,15 @@ void Evolve_LCAO_Matrix::norm_wfc( GlobalV::ofs_running << endl; } - delete[] tmp1; delete[] Cij; - } - -void Evolve_LCAO_Matrix::compute_ekb( - const int nband, - const int nlocal, - const std::complex* Htmp, - const std::complex* psi_k, - double* ekb) const +void Evolve_LCAO_Matrix::compute_ekb(const int nband, + const int nlocal, + const std::complex* Htmp, + const std::complex* psi_k, + double* ekb) const { complex* tmp1 = new complex[this->ParaV->nloc_wfc]; @@ -853,50 +1350,45 @@ void Evolve_LCAO_Matrix::compute_ekb( complex* Eij = new complex[this->ParaV->nloc]; ModuleBase::GlobalFunc::ZEROS(Eij, this->ParaV->nloc); - ScalapackConnector::gemm( - 'N', - 'N', - nlocal, - nband, - nlocal, - 1.0, - Htmp, - 1, - 1, - this->ParaV->desc, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - tmp1, - 1, - 1, - this->ParaV->desc_wfc - ); - - ScalapackConnector::gemm( - 'C', - 'N', - nband, - nband, - nlocal, - 1.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc, - tmp1, - 1, - 1, - this->ParaV->desc_wfc, - 0.0, - Eij, - 1, - 1, - this->ParaV->desc - ); - + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + Htmp, + 1, + 1, + this->ParaV->desc, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + tmp1, + 1, + 1, + this->ParaV->desc_wfc); + + ScalapackConnector::gemm('C', + 'N', + nband, + nband, + nlocal, + 1.0, + psi_k, + 1, + 1, + this->ParaV->desc_wfc, + tmp1, + 1, + 1, + this->ParaV->desc_wfc, + 0.0, + Eij, + 1, + 1, + this->ParaV->desc_Eij); if (ELEC_evolve::td_print_eij > 0.0) { @@ -904,9 +1396,9 @@ void Evolve_LCAO_Matrix::compute_ekb( << "------------------------------------------------------------------------------------------------" << endl; GlobalV::ofs_running << " Eij:" << endl; - for (int i = 0; i < this->ParaV->ncol; i++) + for (int i = 0; i < this->ParaV->nrow_bands; i++) { - for (int j = 0; j < this->ParaV->nrow; j++) + for (int j = 0; j < this->ParaV->ncol_bands; j++) { double aa, bb; aa = Eij[i * this->ParaV->ncol + j].real(); @@ -927,12 +1419,11 @@ void Evolve_LCAO_Matrix::compute_ekb( << endl; } - int info; int myid; int naroc[2]; MPI_Comm_rank(this->ParaV->comm_2D, &myid); - + double* Eii = new double[nband]; ModuleBase::GlobalFunc::ZEROS(Eii, nband); for (int iprow = 0; iprow < this->ParaV->dim0; ++iprow) @@ -966,11 +1457,10 @@ void Evolve_LCAO_Matrix::compute_ekb( } // loop ipcol } // loop iprow info = MPI_Allreduce(Eii, ekb, nband, MPI_DOUBLE, MPI_SUM, this->ParaV->comm_2D); - + delete[] tmp1; delete[] Eij; delete[] Eii; } - #endif diff --git a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h index dd13e45e2f..bcf470a771 100644 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h @@ -22,7 +22,10 @@ class Evolve_LCAO_Matrix hamilt::Hamilt* p_hamilt, psi::Psi>* psi_k, psi::Psi>* psi_k_laststep, - double* ekb) const; + std::complex* H_laststep, + double* ekb, + int htype, + int propagator) const; private: // LCAO_Matrix* LM; @@ -38,15 +41,42 @@ class Evolve_LCAO_Matrix const int nband, const int nlocal, const std::complex* psi_k_laststep, + const std::complex* H_laststep, hamilt::Hamilt* p_hamilt, std::complex* psi_k, - double* ekb) const; + double* ekb, + int htype, + int propagator) const; + + void half_Hmatrix( + const int nband, + const int nlocal, + std::complex* Htmp, + const std::complex* H_laststep, + const int print_matrix) const; + + void compute_U_operator_CN2( + const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const; - void compute_U_operator( + void compute_U_operator_taylor( + const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const; + + void compute_U_operator_etrs( const int nband, const int nlocal, const std::complex* Stmp, const std::complex* Htmp, + const std::complex* H_laststep, std::complex* U_operator, const int print_matrix) const; @@ -55,7 +85,8 @@ class Evolve_LCAO_Matrix const int nlocal, const std::complex* U_operator, const std::complex* psi_k_laststep, - std::complex* psi_k) const; + std::complex* psi_k, + const int print_matrix) const; void norm_wfc( const int nband, diff --git a/source/module_io/dipole_io.h b/source/module_io/dipole_io.h index fbd3d8f0a2..4817f0e3c1 100644 --- a/source/module_io/dipole_io.h +++ b/source/module_io/dipole_io.h @@ -6,11 +6,14 @@ namespace ModuleIO { void write_dipole(const double *rho_save, - const int &is, - const int &istep, - const std::string &fn, - const int &precision = 11, - const bool for_plot = false); + const int &is, + const int &istep, + const std::string &fn, + const int &precision = 11, + const bool for_plot = false); + +double prepare(const UnitCell &cell, int &dir); + } // namespace ModuleIO #endif diff --git a/source/module_io/input.cpp b/source/module_io/input.cpp index 585518fbb7..76f1967ae5 100644 --- a/source/module_io/input.cpp +++ b/source/module_io/input.cpp @@ -405,12 +405,11 @@ void Input::Default(void) // tddft //---------------------------------------------------------- td_force_dt = 0.02; - td_val_elec_01 = 1; - td_val_elec_02 = 1; - td_val_elec_03 = 1; td_vext = false; td_vext_dire = "1"; + propagator = 0; + out_dipole = false; out_efield = false; @@ -1524,18 +1523,6 @@ bool Input::Read(const std::string &fn) { read_value(ifs, td_force_dt); } - else if (strcmp("td_val_elec_01", word) == 0) - { - read_value(ifs, td_val_elec_01); - } - else if (strcmp("td_val_elec_02", word) == 0) - { - read_value(ifs, td_val_elec_02); - } - else if (strcmp("td_val_elec_03", word) == 0) - { - read_value(ifs, td_val_elec_03); - } else if (strcmp("td_vext", word) == 0) { read_value(ifs, td_vext); @@ -1560,6 +1547,10 @@ bool Input::Read(const std::string &fn) { read_value(ifs, td_edm); } + else if (strcmp("propagator", word) == 0) + { + read_value(ifs, propagator); + } else if (strcmp("td_stype", word) == 0) { read_value(ifs, td_stype); @@ -2982,12 +2973,10 @@ void Input::Bcast() Parallel_Common::bcast_int(vdw_cutoff_period.y); Parallel_Common::bcast_int(vdw_cutoff_period.z); // Fuxiang He add 2016-10-26 - Parallel_Common::bcast_int(td_val_elec_01); - Parallel_Common::bcast_int(td_val_elec_02); - Parallel_Common::bcast_int(td_val_elec_03); Parallel_Common::bcast_double(td_force_dt); Parallel_Common::bcast_bool(td_vext); Parallel_Common::bcast_string(td_vext_dire); + Parallel_Common::bcast_int(propagator); Parallel_Common::bcast_int(td_stype); Parallel_Common::bcast_string(td_ttype); Parallel_Common::bcast_int(td_tstart); diff --git a/source/module_io/input.h b/source/module_io/input.h index 99f75726da..df96fa3501 100644 --- a/source/module_io/input.h +++ b/source/module_io/input.h @@ -372,9 +372,6 @@ class Input // Fuxiang He add 2016-10-26 //========================================================== double td_force_dt; //"fs" - int td_val_elec_01; // valence electron 01 - int td_val_elec_02; // valence electron 02 - int td_val_elec_03; // valence electron 03 bool td_vext; // add extern potential or not std::string td_vext_dire; // vext direction bool out_dipole; // output the dipole or not @@ -383,6 +380,8 @@ class Input double td_print_eij; // threshold to output Eij elements int td_edm; //0: new edm method 1: old edm method + int propagator; // method of propagator + int td_stype ; //type of space domain 0 : length gauge 1: velocity gauge std::string td_ttype ; //type of time domain diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index e30e37281d..c11a0e5e3e 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -36,20 +36,32 @@ template void Input_Conv::parse_expression(const std::string &fn, s int count = 0; std::string pattern("([0-9]+\\*[0-9.]+|[0-9,.]+)"); std::vector str; - std::string::size_type pos1, pos2; - std::string c = " "; - pos2 = fn.find(c); - pos1 = 0; - while (std::string::npos != pos2) - { - str.push_back(fn.substr(pos1, pos2 - pos1)); - pos1 = pos2 + c.size(); - pos2 = fn.find(c, pos1); - } - if (pos1 != fn.length()) - { - str.push_back(fn.substr(pos1)); - } + std::stringstream ss(fn); + std::string section; + while (ss >> section) { + int index = 0; + if (str.empty()) { + while (index < section.size() && std::isspace(section[index])) { + index++; + } + } + section.erase(0, index); + str.push_back(section); + } + // std::string::size_type pos1, pos2; + // std::string c = " "; + // pos2 = fn.find(c); + // pos1 = 0; + // while (std::string::npos != pos2) + // { + // str.push_back(fn.substr(pos1, pos2 - pos1)); + // pos1 = pos2 + c.size(); + // pos2 = fn.find(c, pos1); + // } + // if (pos1 != fn.length()) + // { + // str.push_back(fn.substr(pos1)); + // } regex_t reg; regcomp(®, pattern.c_str(), REG_EXTENDED); regmatch_t pmatch[1]; @@ -357,9 +369,6 @@ void Input_Conv::Convert(void) //---------------------------------------------------------- #ifdef __LCAO ELEC_evolve::td_force_dt = INPUT.td_force_dt; - ELEC_evolve::td_val_elec_01 = INPUT.td_val_elec_01; - ELEC_evolve::td_val_elec_02 = INPUT.td_val_elec_02; - ELEC_evolve::td_val_elec_03 = INPUT.td_val_elec_03; ELEC_evolve::td_vext = INPUT.td_vext; if (ELEC_evolve::td_vext) { diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h index beae9bb6c8..c8031e85cd 100644 --- a/source/module_io/test/for_testing_input_conv.h +++ b/source/module_io/test/for_testing_input_conv.h @@ -35,9 +35,6 @@ int hsolver::HSolverLCAO::out_mat_dh = 0; int Local_Orbital_Charge::out_dm = 0; int Local_Orbital_Charge::out_dm1 = 0; double ELEC_evolve::td_force_dt; -int ELEC_evolve::td_val_elec_01; -int ELEC_evolve::td_val_elec_02; -int ELEC_evolve::td_val_elec_03; bool ELEC_evolve::td_vext; std::vector ELEC_evolve::td_vext_dire_case; bool ELEC_evolve::out_dipole; diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp index 181b3eed6d..0de993b1f5 100644 --- a/source/module_io/test/input_conv_test.cpp +++ b/source/module_io/test/input_conv_test.cpp @@ -127,9 +127,6 @@ TEST_F(InputConvTest, Conv) EXPECT_DOUBLE_EQ(elecstate::Gatefield::block_height,0.1); EXPECT_EQ(ELEC_evolve::td_force_dt,0.02); - EXPECT_EQ(ELEC_evolve::td_val_elec_01,1); - EXPECT_EQ(ELEC_evolve::td_val_elec_02,1); - EXPECT_EQ(ELEC_evolve::td_val_elec_03,1); EXPECT_EQ(ELEC_evolve::td_vext,false); EXPECT_EQ(ELEC_evolve::out_dipole,false); EXPECT_EQ(ELEC_evolve::out_efield,false); diff --git a/source/module_io/test/input_test.cpp b/source/module_io/test/input_test.cpp index f06ec30c54..561d28c0f5 100644 --- a/source/module_io/test/input_test.cpp +++ b/source/module_io/test/input_test.cpp @@ -252,11 +252,9 @@ TEST_F(InputTest, Default) EXPECT_DOUBLE_EQ(INPUT.soc_lambda,1.0); EXPECT_EQ(INPUT.input_error,0); EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); - EXPECT_EQ(INPUT.td_val_elec_01,1); - EXPECT_EQ(INPUT.td_val_elec_02,1); - EXPECT_EQ(INPUT.td_val_elec_03,1); EXPECT_FALSE(INPUT.td_vext); EXPECT_EQ(INPUT.td_vext_dire,"1"); + EXPECT_EQ(INPUT.propagator,0); EXPECT_EQ(INPUT.td_stype,0); EXPECT_EQ(INPUT.td_ttype,"0"); EXPECT_EQ(INPUT.td_tstart,1); @@ -589,11 +587,9 @@ TEST_F(InputTest, Read) EXPECT_DOUBLE_EQ(INPUT.soc_lambda,1.0); EXPECT_EQ(INPUT.input_error,0); EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); - EXPECT_EQ(INPUT.td_val_elec_01,1); - EXPECT_EQ(INPUT.td_val_elec_02,1); - EXPECT_EQ(INPUT.td_val_elec_03,1); EXPECT_EQ(INPUT.td_vext,0); // EXPECT_EQ(INPUT.td_vext_dire,"1"); + EXPECT_EQ(INPUT.propagator,0); EXPECT_EQ(INPUT.td_stype,0); // EXPECT_EQ(INPUT.td_ttype,"0"); EXPECT_EQ(INPUT.td_tstart,1); diff --git a/source/module_io/test/input_test_para.cpp b/source/module_io/test/input_test_para.cpp index f80762a374..315b560cec 100644 --- a/source/module_io/test/input_test_para.cpp +++ b/source/module_io/test/input_test_para.cpp @@ -258,11 +258,9 @@ TEST_F(InputParaTest,Bcast) EXPECT_DOUBLE_EQ(INPUT.soc_lambda,1.0); EXPECT_EQ(INPUT.input_error,0); EXPECT_DOUBLE_EQ(INPUT.td_force_dt,0.02); - EXPECT_EQ(INPUT.td_val_elec_01,1); - EXPECT_EQ(INPUT.td_val_elec_02,1); - EXPECT_EQ(INPUT.td_val_elec_03,1); EXPECT_FALSE(INPUT.td_vext); EXPECT_EQ(INPUT.td_vext_dire,"1"); + EXPECT_EQ(INPUT.propagator,0); EXPECT_EQ(INPUT.td_stype,0); EXPECT_EQ(INPUT.td_ttype,"0"); EXPECT_EQ(INPUT.td_tstart,1); diff --git a/source/module_io/test/support/INPUT b/source/module_io/test/support/INPUT index 9078dea961..3c57394565 100644 --- a/source/module_io/test/support/INPUT +++ b/source/module_io/test/support/INPUT @@ -282,11 +282,9 @@ exx_opt_orb_tolerence 0 # #Parameters (16.tddft) td_force_dt 0.02 #time of force change -td_val_elec_01 1 #td_val_elec_01 -td_val_elec_02 1 #td_val_elec_02 -td_val_elec_03 1 #td_val_elec_03 td_vext 0 #add extern potential or not td_vext_dire 1 #extern potential direction +propagator 0 # method of propagator td_stype 0 #space domain type td_ttype 0 #time domain type td_tstart 1 #the start step of electric field diff --git a/source/module_io/test/support/witestfile b/source/module_io/test/support/witestfile index c2e4563b46..7cb5e0f450 100644 --- a/source/module_io/test/support/witestfile +++ b/source/module_io/test/support/witestfile @@ -1,356 +1,354 @@ -INPUT_PARAMETERS -#Parameters (1.General) -suffix autotest #the name of main output directory -latname none #the name of lattice name -stru_file ./support/STRU #the filename of file containing atom positions -kpoint_file KPT #the name of file containing k points -pseudo_dir ../../PP_ORB/ #the directory containing pseudo files -orbital_dir ../../PP_ORB/ #the directory containing orbital files -pseudo_rcut 15 #cut-off radius for radial integration -pseudo_mesh 0 #0: use our own mesh to do radial renormalization; 1: use mesh as in QE -lmaxmax 2 #maximum of l channels used -dft_functional hse #exchange correlation functional -xc_temperature 0 #temperature for finite temperature functionals -calculation scf #test; scf; relax; nscf; ienvelope; istate -esolver_type ksdft #the energy solver: ksdft, sdft, ofdft, tddft, lj, dp -ntype 1 #atom species number -nspin 1 #1: single spin; 2: up and down spin; 4: noncollinear spin -kspacing 0 #unit in 1/bohr, should be > 0, default is 0 which means read KPT file -min_dist_coef 0.2 #factor related to the allowed minimum distance between two atoms -nbands 8 #number of bands -nbands_sto 256 #number of stochastic bands -nbands_istate 5 #number of bands around Fermi level for istate calulation -symmetry 1 #the control of symmetry -init_vel False #read velocity from STRU or not -symmetry_prec 1e-05 #accuracy for symmetry -nelec 0 #input number of electrons -out_mul 0 # mulliken charge or not -noncolin 0 #using non-collinear-spin -lspinorb 0 #consider the spin-orbit interaction -kpar 1 #devide all processors into kpar groups and k points will be distributed among each group -bndpar 1 #devide all processors into bndpar groups and bands will be distributed among each group -out_freq_elec 0 #the frequency ( >= 0) of electronic iter to output charge density and wavefunction. 0: output only when converged -dft_plus_dmft 0 #true:DFT+DMFT; false: standard DFT calcullation(default) -rpa 0 #true:generate output files used in rpa calculation; false:(default) -printe 100 #Print out energy for each band for every printe steps -mem_saver 0 #Only for nscf calculations. if set to 1, then a memory saving technique will be used for many k point calculations. -diago_proc 4 #the number of procs used to do diagonalization -nbspline -1 #the order of B-spline basis -wannier_card none #input card for wannier functions -soc_lambda 1 #The fraction of averaged SOC pseudopotential is given by (1-soc_lambda) -cal_force 0 #if calculate the force at the end of the electronic iteration -out_freq_ion 0 #the frequency ( >= 0 ) of ionic step to output charge density and wavefunction. 0: output only when ion steps are finished -device cpu #the computing device for ABACUS - -#Parameters (2.PW) -ecutwfc 20 ##energy cutoff for wave functions -pw_diag_thr 0.01 #threshold for eigenvalues is cg electron iterations -scf_thr 1e-08 #charge density error -init_wfc atomic #start wave functions are from 'atomic', 'atomic+random', 'random' or 'file' -init_chg atomic #start charge is from 'atomic' or file -chg_extrap atomic #atomic; first-order; second-order; dm:coefficients of SIA -out_chg FALSE #>0 output charge density for selected electron steps -out_pot 2 #output realspace potential -out_wfc_pw 0 #output wave functions -out_wfc_r 0 #output wave functions in realspace -out_dos 0 #output energy and dos -out_band false #output energy and band structure -out_proj_band FaLse #output projected band structure -restart_save f #print to disk every step for restart -restart_load F #restart from disk -read_file_dir auto #directory of files for reading -nx 0 #number of points along x axis for FFT grid -ny 0 #number of points along y axis for FFT grid -nz 0 #number of points along z axis for FFT grid -cell_factor 1.2 #used in the construction of the pseudopotential tables -pw_seed 1 #random seed for initializing wave functions - -#Parameters (3.Stochastic DFT) -method_sto 3 #1: slow and save memory, 2: fast and waste memory -npart_sto 1 #Reduce memory when calculating Stochastic DOS -nbands_sto 256 #number of stochstic orbitals -nche_sto 100 #Chebyshev expansion orders -emin_sto 0 #trial energy to guess the lower bound of eigen energies of the Hamitonian operator -emax_sto 0 #trial energy to guess the upper bound of eigen energies of the Hamitonian operator -seed_sto 0 #the random seed to generate stochastic orbitals -initsto_freq 0 #frequency to generate new stochastic orbitals when running md -cal_cond 0 #calculate electronic conductivities -cond_nche 20 #orders of Chebyshev expansions for conductivities -cond_dw 0.1 #frequency interval for conductivities -cond_wcut 10 #cutoff frequency (omega) for conductivities -cond_dt 0.07 #control the t interval -cond_dtbatch 2 #control dt batch -cond_fwhm 0.3 #FWHM for conductivities -cond_nonlocal 1 #Nonlocal effects for conductivities - -#Parameters (4.Relaxation) -ks_solver genelpa #cg; dav; lapack; genelpa; scalapack_gvx; cusolver -scf_nmax 50 ##number of electron iterations -relax_nmax 1 #number of ion iteration steps -out_stru 0 #output the structure files after each ion step -force_thr 0.001 #force threshold, unit: Ry/Bohr -force_thr_ev 0.0257112 #force threshold, unit: eV/Angstrom -force_thr_ev2 0 #force invalid threshold, unit: eV/Angstrom -relax_cg_thr 0.5 #threshold for switching from cg to bfgs, unit: eV/Angstrom -stress_thr 0.01 #stress threshold -press1 0 #target pressure, unit: KBar -press2 0 #target pressure, unit: KBar -press3 0 #target pressure, unit: KBar -relax_bfgs_w1 0.01 #wolfe condition 1 for bfgs -relax_bfgs_w2 0.5 #wolfe condition 2 for bfgs -relax_bfgs_rmax 0.8 #maximal trust radius, unit: Bohr -relax_bfgs_rmin 1e-05 #minimal trust radius, unit: Bohr -relax_bfgs_init 0.5 #initial trust radius, unit: Bohr -cal_stress 0 #calculate the stress or not -fixed_axes None #which axes are fixed -fixed_ibrav 0 #whether to preseve lattice type during relaxation -fixed_atoms 0 #whether to preseve direct coordinates of atoms during relaxation -relax_method cg #bfgs; sd; cg; cg_bfgs; -relax_new TRUE #whether to use the new relaxation method -relax_scale_force 0.5 #controls the size of the first CG step if relax_new is true -out_level ie #ie(for electrons); i(for ions); -out_dm 0 #>0 output density matrix -deepks_out_labels 0 #>0 compute descriptor for deepks -deepks_scf 0 #>0 add V_delta to Hamiltonian -deepks_bandgap 0 #>0 for bandgap label -deepks_out_unittest 0 #if set 1, prints intermediate quantities that shall be used for making unit test -deepks_model #file dir of traced pytorch model: 'model.ptg - -#Parameters (5.LCAO) -basis_type lcao #PW; LCAO in pw; LCAO -nb2d 0 #2d distribution of atoms -gamma_only F #Only for localized orbitals set and gamma point. If set to 1, a fast algorithm is used -search_radius -1 #input search radius (Bohr) -search_pbc 1 #input periodic boundary condition -lcao_ecut 20 #energy cutoff for LCAO -lcao_dk 0.01 #delta k for 1D integration in LCAO -lcao_dr 0.01 #delta r for 1D integration in LCAO -lcao_rmax 30 #max R for 1D two-center integration table -out_mat_hs 0 #output H and S matrix -out_mat_hs2 0 #output H(R) and S(R) matrix -out_interval 1 #interval for printing H(R) and S(R) matrix during MD -out_app_flag 0 #interval for printing H(R) and S(R) matrix during MD -out_element_info 0 #output (projected) wavefunction of each element -out_mat_r 0 #output r(R) matrix -out_wfc_lcao 0 #ouput LCAO wave functions -bx 2 #division of an element grid in FFT grid along x -by 2 #division of an element grid in FFT grid along y -bz 2 #division of an element grid in FFT grid along z - -bessel_nao_smooth 1 -bessel_nao_sigma 0.1 -bessel_nao_ecut default -bessel_nao_rcut 6.0 #-1.0 for forcing manual setting -bessel_nao_tolerence 1E-12 - -bessel_descriptor_lmax 2 # -1 for forcing manual setting -bessel_descriptor_smooth 1 -bessel_descriptor_sigma 0.1 -bessel_descriptor_ecut default -bessel_descriptor_rcut 6.0 #-1.0 for forcing manual setting -bessel_descriptor_tolerence 1E-12 - -#Parameters (6.Smearing) -smearing_method gauss #type of smearing_method: gauss; fd; fixed; mp; mp2; mv -smearing_sigma 0.002 #energy range for smearing - -#Parameters (7.Charge Mixing) -mixing_type pulay #plain; pulay; broyden -mixing_beta 0.7 #mixing parameter: 0 means no new charge -mixing_ndim 8 #mixing dimension in pulay -mixing_gg0 0 #mixing parameter in kerker - -#Parameters (8.DOS) -dos_emin_ev -15 #minimal range for dos -dos_emax_ev 15 #maximal range for dos -dos_edelta_ev 0.01 #delta energy for dos -dos_scale 0.01 #scale dos range by -dos_sigma 0.07 #gauss b coefficeinet(default=0.07) -dos_nche 100 #orders of Chebyshev expansions for dos - -#Parameters (9.Molecular dynamics) -md_type nvt #choose ensemble -md_thermostat nhc #choose thermostat -md_nstep 10 #md steps -md_dt 1 #time step -md_tchain 1 #number of Nose-Hoover chains -md_tfirst -1 #temperature first -md_tlast -1 #temperature last -md_dumpfreq 1 #The period to dump MD information -md_restartfreq 5 #The period to output MD restart information -md_seed -1 #random seed for MD -md_restart 0 #whether restart -lj_rcut 8.5 #cutoff radius of LJ potential -lj_epsilon 0.01032 #the value of epsilon for LJ potential -lj_sigma 3.405 #the value of sigma for LJ potential -pot_file graph.pb #the filename of potential files for CMD such as DP -msst_direction 2 #the direction of shock wave -msst_vel 0 #the velocity of shock wave -msst_vis 0 #artificial viscosity -msst_tscale 0.01 #reduction in initial temperature -msst_qmass -1 #mass of thermostat -md_tfreq 0 #oscillation frequency, used to determine qmass of NHC -md_damp 1 #damping parameter (time units) used to add force in Langevin method -md_nraise 1 #parameters used when md_type=nvt -md_tolerance 100 #tolerance for velocity rescaling (K) -md_pmode iso #NPT ensemble mode: iso, aniso, tri -md_pcouple none #whether couple different components: xyz, xy, yz, xz, none -md_pchain 1 #num of thermostats coupled with barostat -md_pfirst -1 #initial target pressure -md_plast -1 #final target pressure -md_pfreq 0 #oscillation frequency, used to determine qmass of thermostats coupled with barostat -dump_force 0 #output atomic forces into the file MD_dump or not. -dump_vel 0 #output atomic velocities into the file MD_dump or not -dump_virial 0 #output lattice virial into the file MD_dump or not - -#Parameters (10.Electric field and dipole correction) -efield_flag 0 #add electric field -dip_cor_flag 0 #dipole correction -efield_dir 2 #the direction of the electric field or dipole correction -efield_pos_max 0.5 #position of the maximum of the saw-like potential along crystal axis efield_dir -efield_pos_dec 0.1 #zone in the unit cell where the saw-like potential decreases -efield_amp 0 #amplitude of the electric field - -#Parameters (11.Gate field) -gate_flag 0 #compensating charge or not -zgate 0.5 #position of charged plate -relax 0 #allow relaxation along the specific direction -block 0 #add a block potential or not -block_down 0.45 #low bound of the block -block_up 0.55 #high bound of the block -block_height 0.1 #height of the block - -#Parameters (12.Test) -out_alllog F #output information for each processor, when parallel -nurse 0 #for coders -colour 0 #for coders, make their live colourful -t_in_h 1 #calculate the kinetic energy or not -vl_in_h 1 #calculate the local potential or not -vnl_in_h 1 #calculate the nonlocal potential or not -vh_in_h 1 #calculate the hartree potential or not -vion_in_h 1 #calculate the local ionic potential or not -test_force 0 #test the force -test_stress 0 #test the force -test_skip_ewald 0 #skip ewald energy - -#Parameters (13.vdW Correction) -vdw_method d2 #the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj -vdw_s6 default #scale parameter of d2/d3_0/d3_bj -vdw_s8 default #scale parameter of d3_0/d3_bj -vdw_a1 default #damping parameter of d3_0/d3_bj -vdw_a2 default #damping parameter of d3_bj -vdw_d 20 #damping parameter of d2 -vdw_abc 0 #third-order term? -vdw_C6_file default #filename of C6 -vdw_C6_unit Jnm6/mol #unit of C6, Jnm6/mol or eVA6 -vdw_R0_file default #filename of R0 -vdw_R0_unit A #unit of R0, A or Bohr -vdw_cutoff_type radius #expression model of periodic structure, radius or period -vdw_cutoff_radius default #radius cutoff for periodic structure -vdw_radius_unit Bohr #unit of radius cutoff for periodic structure -vdw_cn_thr 40 #radius cutoff for cn -vdw_cn_thr_unit Bohr #unit of cn_thr, Bohr or Angstrom -vdw_cutoff_period 3 3 3 #periods of periodic structure - -#Parameters (14.exx) -exx_hybrid_alpha default # -exx_hse_omega 0.11 # -exx_separate_loop 1 #0 or 1 -exx_hybrid_step 100 # -exx_lambda 0.3 # -exx_real_number default #0 or 1 -exx_pca_threshold 0 # -exx_c_threshold 0 # -exx_v_threshold 0 # -exx_dm_threshold 0 # -exx_schwarz_threshold 0 # -exx_cauchy_threshold 0 # -exx_c_grad_threshold 0 # -exx_v_grad_threshold 0 # -exx_ccp_threshold 1e-08 # -exx_ccp_rmesh_times default # -exx_distribute_type htime #htime or kmeans1 or kmeans2 -exx_opt_orb_lmax 0 # -exx_opt_orb_ecut 0 # -exx_opt_orb_tolerence 0 # - -#Parameters (16.tddft) -td_force_dt 0.02 #time of force change -td_val_elec_01 1 #td_val_elec_01 -td_val_elec_02 1 #td_val_elec_02 -td_val_elec_03 1 #td_val_elec_03 -td_vext 0 #add extern potential or not -td_vext_dire 1 #extern potential direction -td_stype 0 #space domain type -td_ttype 0 #time domain type -td_tstart 1 #the start step of electric field -td_tend 1000 #the start step of electric field -td_lcut1 0.05 # separator in length gauge -td_lcut2 0.95 # separator in length gauge -# parameters of Gauss electric field -td_gauss_freq 22.13 # fs^-1 -td_gauss_phase 0.0 -td_gauss_sigma 30.0 # fs -td_gauss_t0 100.0 -td_gauss_amp 0.25 # V/A -# parameters of Trapezoid electric field -td_trape_freq 1.60 # fs^-1 -td_trape_phase 0.0 -td_trape_t1 1875 -td_trape_t2 5625 -td_trape_t3 7500 -td_trape_amp 2.74 # V/A -#parameters of Trigonometric electric field -td_trigo_freq1 1.164656 # fs^-1 -td_trigo_freq2 0.029116 # fs^-1 -td_trigo_phase1 0.0 -td_trigo_phase2 0.0 -td_trigo_amp 2.74 # V/A -#parameters of Heaviside electric field -td_heavi_t0 100 -td_heavi_amp 1.0 # V/A -td_print_eij -1.0 # threshold to output Eij elements -td_edm 0 # 0: new edm method 1: old edm method -out_dipole 0 #output dipole or not -out_efield 0 #output efield or not -ocp 0 #change occupation or not -ocp_set none #set occupation - -#Parameters (17.berry_wannier) -berry_phase 0 #calculate berry phase or not -gdir 3 #calculate the polarization in the direction of the lattice vector -towannier90 0 #use wannier90 code interface or not -nnkpfile seedname.nnkp #the wannier90 code nnkp file name -wannier_spin up #calculate spin in wannier90 code interface - -#Parameters (18.implicit_solvation) -imp_sol 0 #calculate implicit solvation correction or not -eb_k 80 #the relative permittivity of the bulk solvent -tau 1.0798e-05 #the effective surface tension parameter -sigma_k 0.6 # the width of the diffuse cavity -nc_k 0.00037 # the cut-off charge density - -#Parameters (19.orbital free density functional theory) -of_kinetic vw #kinetic energy functional, such as tf, vw, wt -of_method tn #optimization method used in OFDFT, including cg1, cg2, tn (default) -of_conv energy #the convergence criterion, potential, energy (default), or both -of_tole 1e-06 #tolerance of the energy change (in Ry) for determining the convergence, default=2e-6 Ry -of_tolp 1e-05 #tolerance of potential for determining the convergence, default=1e-5 in a.u. -of_tf_weight 1 #weight of TF KEDF -of_vw_weight 1 #weight of vW KEDF -of_wt_alpha 0.833333 #parameter alpha of WT KEDF -of_wt_beta 0.833333 #parameter beta of WT KEDF -of_wt_rho0 1 #the average density of system, used in WT KEDF, in Bohr^-3 -of_hold_rho0 0 #If set to 1, the rho0 will be fixed even if the volume of system has changed, it will be set to 1 automaticly if of_wt_rho0 is not zero -of_full_pw 0 #If set to 1, ecut will be ignored when collect planewaves, so that all planewaves will be used -of_full_pw_dim 0 #If of_full_pw = true, dimention of FFT is testricted to be (0) either odd or even; (1) odd only; (2) even only -of_read_kernel 0 #If set to 1, the kernel of WT KEDF will be filled from file of_kernel_file, not from formula. Only usable for WT KEDF -of_kernel_file WTkernel.txt #The name of WT kernel file. - -#Parameters (19.dft+u) -dft_plus_u 0 #true:DFT+U correction; false: standard DFT calcullation(default) -yukawa_lambda -1 #default:0.0 -yukawa_potential 0 #default: false -omc 0 #the mode of occupation matrix control -hubbard_u 0 #Hubbard Coulomb interaction parameter U(ev) -orbital_corr -1 #which correlated orbitals need corrected ; d:2 ,f:3, do not need correction:-1 +INPUT_PARAMETERS +#Parameters (1.General) +suffix autotest #the name of main output directory +latname none #the name of lattice name +stru_file ./support/STRU #the filename of file containing atom positions +kpoint_file KPT #the name of file containing k points +pseudo_dir ../../PP_ORB/ #the directory containing pseudo files +orbital_dir ../../PP_ORB/ #the directory containing orbital files +pseudo_rcut 15 #cut-off radius for radial integration +pseudo_mesh 0 #0: use our own mesh to do radial renormalization; 1: use mesh as in QE +lmaxmax 2 #maximum of l channels used +dft_functional hse #exchange correlation functional +xc_temperature 0 #temperature for finite temperature functionals +calculation scf #test; scf; relax; nscf; ienvelope; istate +esolver_type ksdft #the energy solver: ksdft, sdft, ofdft, tddft, lj, dp +ntype 1 #atom species number +nspin 1 #1: single spin; 2: up and down spin; 4: noncollinear spin +kspacing 0 #unit in 1/bohr, should be > 0, default is 0 which means read KPT file +min_dist_coef 0.2 #factor related to the allowed minimum distance between two atoms +nbands 8 #number of bands +nbands_sto 256 #number of stochastic bands +nbands_istate 5 #number of bands around Fermi level for istate calulation +symmetry 1 #the control of symmetry +init_vel False #read velocity from STRU or not +symmetry_prec 1e-05 #accuracy for symmetry +nelec 0 #input number of electrons +out_mul 0 # mulliken charge or not +noncolin 0 #using non-collinear-spin +lspinorb 0 #consider the spin-orbit interaction +kpar 1 #devide all processors into kpar groups and k points will be distributed among each group +bndpar 1 #devide all processors into bndpar groups and bands will be distributed among each group +out_freq_elec 0 #the frequency ( >= 0) of electronic iter to output charge density and wavefunction. 0: output only when converged +dft_plus_dmft 0 #true:DFT+DMFT; false: standard DFT calcullation(default) +rpa 0 #true:generate output files used in rpa calculation; false:(default) +printe 100 #Print out energy for each band for every printe steps +mem_saver 0 #Only for nscf calculations. if set to 1, then a memory saving technique will be used for many k point calculations. +diago_proc 4 #the number of procs used to do diagonalization +nbspline -1 #the order of B-spline basis +wannier_card none #input card for wannier functions +soc_lambda 1 #The fraction of averaged SOC pseudopotential is given by (1-soc_lambda) +cal_force 0 #if calculate the force at the end of the electronic iteration +out_freq_ion 0 #the frequency ( >= 0 ) of ionic step to output charge density and wavefunction. 0: output only when ion steps are finished +device cpu #the computing device for ABACUS + +#Parameters (2.PW) +ecutwfc 20 ##energy cutoff for wave functions +pw_diag_thr 0.01 #threshold for eigenvalues is cg electron iterations +scf_thr 1e-08 #charge density error +init_wfc atomic #start wave functions are from 'atomic', 'atomic+random', 'random' or 'file' +init_chg atomic #start charge is from 'atomic' or file +chg_extrap atomic #atomic; first-order; second-order; dm:coefficients of SIA +out_chg FALSE #>0 output charge density for selected electron steps +out_pot 2 #output realspace potential +out_wfc_pw 0 #output wave functions +out_wfc_r 0 #output wave functions in realspace +out_dos 0 #output energy and dos +out_band false #output energy and band structure +out_proj_band FaLse #output projected band structure +restart_save f #print to disk every step for restart +restart_load F #restart from disk +read_file_dir auto #directory of files for reading +nx 0 #number of points along x axis for FFT grid +ny 0 #number of points along y axis for FFT grid +nz 0 #number of points along z axis for FFT grid +cell_factor 1.2 #used in the construction of the pseudopotential tables +pw_seed 1 #random seed for initializing wave functions + +#Parameters (3.Stochastic DFT) +method_sto 3 #1: slow and save memory, 2: fast and waste memory +npart_sto 1 #Reduce memory when calculating Stochastic DOS +nbands_sto 256 #number of stochstic orbitals +nche_sto 100 #Chebyshev expansion orders +emin_sto 0 #trial energy to guess the lower bound of eigen energies of the Hamitonian operator +emax_sto 0 #trial energy to guess the upper bound of eigen energies of the Hamitonian operator +seed_sto 0 #the random seed to generate stochastic orbitals +initsto_freq 0 #frequency to generate new stochastic orbitals when running md +cal_cond 0 #calculate electronic conductivities +cond_nche 20 #orders of Chebyshev expansions for conductivities +cond_dw 0.1 #frequency interval for conductivities +cond_wcut 10 #cutoff frequency (omega) for conductivities +cond_dt 0.07 #control the t interval +cond_dtbatch 2 #control dt batch +cond_fwhm 0.3 #FWHM for conductivities +cond_nonlocal 1 #Nonlocal effects for conductivities + +#Parameters (4.Relaxation) +ks_solver genelpa #cg; dav; lapack; genelpa; scalapack_gvx; cusolver +scf_nmax 50 ##number of electron iterations +relax_nmax 1 #number of ion iteration steps +out_stru 0 #output the structure files after each ion step +force_thr 0.001 #force threshold, unit: Ry/Bohr +force_thr_ev 0.0257112 #force threshold, unit: eV/Angstrom +force_thr_ev2 0 #force invalid threshold, unit: eV/Angstrom +relax_cg_thr 0.5 #threshold for switching from cg to bfgs, unit: eV/Angstrom +stress_thr 0.01 #stress threshold +press1 0 #target pressure, unit: KBar +press2 0 #target pressure, unit: KBar +press3 0 #target pressure, unit: KBar +relax_bfgs_w1 0.01 #wolfe condition 1 for bfgs +relax_bfgs_w2 0.5 #wolfe condition 2 for bfgs +relax_bfgs_rmax 0.8 #maximal trust radius, unit: Bohr +relax_bfgs_rmin 1e-05 #minimal trust radius, unit: Bohr +relax_bfgs_init 0.5 #initial trust radius, unit: Bohr +cal_stress 0 #calculate the stress or not +fixed_axes None #which axes are fixed +fixed_ibrav 0 #whether to preseve lattice type during relaxation +fixed_atoms 0 #whether to preseve direct coordinates of atoms during relaxation +relax_method cg #bfgs; sd; cg; cg_bfgs; +relax_new TRUE #whether to use the new relaxation method +relax_scale_force 0.5 #controls the size of the first CG step if relax_new is true +out_level ie #ie(for electrons); i(for ions); +out_dm 0 #>0 output density matrix +deepks_out_labels 0 #>0 compute descriptor for deepks +deepks_scf 0 #>0 add V_delta to Hamiltonian +deepks_bandgap 0 #>0 for bandgap label +deepks_out_unittest 0 #if set 1, prints intermediate quantities that shall be used for making unit test +deepks_model #file dir of traced pytorch model: 'model.ptg + +#Parameters (5.LCAO) +basis_type lcao #PW; LCAO in pw; LCAO +nb2d 0 #2d distribution of atoms +gamma_only F #Only for localized orbitals set and gamma point. If set to 1, a fast algorithm is used +search_radius -1 #input search radius (Bohr) +search_pbc 1 #input periodic boundary condition +lcao_ecut 20 #energy cutoff for LCAO +lcao_dk 0.01 #delta k for 1D integration in LCAO +lcao_dr 0.01 #delta r for 1D integration in LCAO +lcao_rmax 30 #max R for 1D two-center integration table +out_mat_hs 0 #output H and S matrix +out_mat_hs2 0 #output H(R) and S(R) matrix +out_interval 1 #interval for printing H(R) and S(R) matrix during MD +out_app_flag 0 #interval for printing H(R) and S(R) matrix during MD +out_element_info 0 #output (projected) wavefunction of each element +out_mat_r 0 #output r(R) matrix +out_wfc_lcao 0 #ouput LCAO wave functions +bx 2 #division of an element grid in FFT grid along x +by 2 #division of an element grid in FFT grid along y +bz 2 #division of an element grid in FFT grid along z + +bessel_nao_smooth 1 +bessel_nao_sigma 0.1 +bessel_nao_ecut default +bessel_nao_rcut 6.0 #-1.0 for forcing manual setting +bessel_nao_tolerence 1E-12 + +bessel_descriptor_lmax 2 # -1 for forcing manual setting +bessel_descriptor_smooth 1 +bessel_descriptor_sigma 0.1 +bessel_descriptor_ecut default +bessel_descriptor_rcut 6.0 #-1.0 for forcing manual setting +bessel_descriptor_tolerence 1E-12 + +#Parameters (6.Smearing) +smearing_method gauss #type of smearing_method: gauss; fd; fixed; mp; mp2; mv +smearing_sigma 0.002 #energy range for smearing + +#Parameters (7.Charge Mixing) +mixing_type pulay #plain; pulay; broyden +mixing_beta 0.7 #mixing parameter: 0 means no new charge +mixing_ndim 8 #mixing dimension in pulay +mixing_gg0 0 #mixing parameter in kerker + +#Parameters (8.DOS) +dos_emin_ev -15 #minimal range for dos +dos_emax_ev 15 #maximal range for dos +dos_edelta_ev 0.01 #delta energy for dos +dos_scale 0.01 #scale dos range by +dos_sigma 0.07 #gauss b coefficeinet(default=0.07) +dos_nche 100 #orders of Chebyshev expansions for dos + +#Parameters (9.Molecular dynamics) +md_type nvt #choose ensemble +md_thermostat nhc #choose thermostat +md_nstep 10 #md steps +md_dt 1 #time step +md_tchain 1 #number of Nose-Hoover chains +md_tfirst -1 #temperature first +md_tlast -1 #temperature last +md_dumpfreq 1 #The period to dump MD information +md_restartfreq 5 #The period to output MD restart information +md_seed -1 #random seed for MD +md_restart 0 #whether restart +lj_rcut 8.5 #cutoff radius of LJ potential +lj_epsilon 0.01032 #the value of epsilon for LJ potential +lj_sigma 3.405 #the value of sigma for LJ potential +pot_file graph.pb #the filename of potential files for CMD such as DP +msst_direction 2 #the direction of shock wave +msst_vel 0 #the velocity of shock wave +msst_vis 0 #artificial viscosity +msst_tscale 0.01 #reduction in initial temperature +msst_qmass -1 #mass of thermostat +md_tfreq 0 #oscillation frequency, used to determine qmass of NHC +md_damp 1 #damping parameter (time units) used to add force in Langevin method +md_nraise 1 #parameters used when md_type=nvt +md_tolerance 100 #tolerance for velocity rescaling (K) +md_pmode iso #NPT ensemble mode: iso, aniso, tri +md_pcouple none #whether couple different components: xyz, xy, yz, xz, none +md_pchain 1 #num of thermostats coupled with barostat +md_pfirst -1 #initial target pressure +md_plast -1 #final target pressure +md_pfreq 0 #oscillation frequency, used to determine qmass of thermostats coupled with barostat +dump_force 0 #output atomic forces into the file MD_dump or not. +dump_vel 0 #output atomic velocities into the file MD_dump or not +dump_virial 0 #output lattice virial into the file MD_dump or not + +#Parameters (10.Electric field and dipole correction) +efield_flag 0 #add electric field +dip_cor_flag 0 #dipole correction +efield_dir 2 #the direction of the electric field or dipole correction +efield_pos_max 0.5 #position of the maximum of the saw-like potential along crystal axis efield_dir +efield_pos_dec 0.1 #zone in the unit cell where the saw-like potential decreases +efield_amp 0 #amplitude of the electric field + +#Parameters (11.Gate field) +gate_flag 0 #compensating charge or not +zgate 0.5 #position of charged plate +relax 0 #allow relaxation along the specific direction +block 0 #add a block potential or not +block_down 0.45 #low bound of the block +block_up 0.55 #high bound of the block +block_height 0.1 #height of the block + +#Parameters (12.Test) +out_alllog F #output information for each processor, when parallel +nurse 0 #for coders +colour 0 #for coders, make their live colourful +t_in_h 1 #calculate the kinetic energy or not +vl_in_h 1 #calculate the local potential or not +vnl_in_h 1 #calculate the nonlocal potential or not +vh_in_h 1 #calculate the hartree potential or not +vion_in_h 1 #calculate the local ionic potential or not +test_force 0 #test the force +test_stress 0 #test the force +test_skip_ewald 0 #skip ewald energy + +#Parameters (13.vdW Correction) +vdw_method d2 #the method of calculating vdw (none ; d2 ; d3_0 ; d3_bj +vdw_s6 default #scale parameter of d2/d3_0/d3_bj +vdw_s8 default #scale parameter of d3_0/d3_bj +vdw_a1 default #damping parameter of d3_0/d3_bj +vdw_a2 default #damping parameter of d3_bj +vdw_d 20 #damping parameter of d2 +vdw_abc 0 #third-order term? +vdw_C6_file default #filename of C6 +vdw_C6_unit Jnm6/mol #unit of C6, Jnm6/mol or eVA6 +vdw_R0_file default #filename of R0 +vdw_R0_unit A #unit of R0, A or Bohr +vdw_cutoff_type radius #expression model of periodic structure, radius or period +vdw_cutoff_radius default #radius cutoff for periodic structure +vdw_radius_unit Bohr #unit of radius cutoff for periodic structure +vdw_cn_thr 40 #radius cutoff for cn +vdw_cn_thr_unit Bohr #unit of cn_thr, Bohr or Angstrom +vdw_cutoff_period 3 3 3 #periods of periodic structure + +#Parameters (14.exx) +exx_hybrid_alpha default # +exx_hse_omega 0.11 # +exx_separate_loop 1 #0 or 1 +exx_hybrid_step 100 # +exx_lambda 0.3 # +exx_real_number default #0 or 1 +exx_pca_threshold 0 # +exx_c_threshold 0 # +exx_v_threshold 0 # +exx_dm_threshold 0 # +exx_schwarz_threshold 0 # +exx_cauchy_threshold 0 # +exx_c_grad_threshold 0 # +exx_v_grad_threshold 0 # +exx_ccp_threshold 1e-08 # +exx_ccp_rmesh_times default # +exx_distribute_type htime #htime or kmeans1 or kmeans2 +exx_opt_orb_lmax 0 # +exx_opt_orb_ecut 0 # +exx_opt_orb_tolerence 0 # + +#Parameters (16.tddft) +td_force_dt 0.02 #time of force change +td_vext 0 #add extern potential or not +td_vext_dire 1 #extern potential direction +propagator 0 # method of propagator +td_stype 0 #space domain type +td_ttype 0 #time domain type +td_tstart 1 #the start step of electric field +td_tend 1000 #the start step of electric field +td_lcut1 0.05 # separator in length gauge +td_lcut2 0.95 # separator in length gauge +# parameters of Gauss electric field +td_gauss_freq 22.13 # fs^-1 +td_gauss_phase 0.0 +td_gauss_sigma 30.0 # fs +td_gauss_t0 100.0 +td_gauss_amp 0.25 # V/A +# parameters of Trapezoid electric field +td_trape_freq 1.60 # fs^-1 +td_trape_phase 0.0 +td_trape_t1 1875 +td_trape_t2 5625 +td_trape_t3 7500 +td_trape_amp 2.74 # V/A +#parameters of Trigonometric electric field +td_trigo_freq1 1.164656 # fs^-1 +td_trigo_freq2 0.029116 # fs^-1 +td_trigo_phase1 0.0 +td_trigo_phase2 0.0 +td_trigo_amp 2.74 # V/A +#parameters of Heaviside electric field +td_heavi_t0 100 +td_heavi_amp 1.0 # V/A +td_print_eij -1.0 # threshold to output Eij elements +td_edm 0 # 0: new edm method 1: old edm method +out_dipole 0 #output dipole or not +out_efield 0 #output efield or not +ocp 0 #change occupation or not +ocp_set none #set occupation + +#Parameters (17.berry_wannier) +berry_phase 0 #calculate berry phase or not +gdir 3 #calculate the polarization in the direction of the lattice vector +towannier90 0 #use wannier90 code interface or not +nnkpfile seedname.nnkp #the wannier90 code nnkp file name +wannier_spin up #calculate spin in wannier90 code interface + +#Parameters (18.implicit_solvation) +imp_sol 0 #calculate implicit solvation correction or not +eb_k 80 #the relative permittivity of the bulk solvent +tau 1.0798e-05 #the effective surface tension parameter +sigma_k 0.6 # the width of the diffuse cavity +nc_k 0.00037 # the cut-off charge density + +#Parameters (19.orbital free density functional theory) +of_kinetic vw #kinetic energy functional, such as tf, vw, wt +of_method tn #optimization method used in OFDFT, including cg1, cg2, tn (default) +of_conv energy #the convergence criterion, potential, energy (default), or both +of_tole 1e-06 #tolerance of the energy change (in Ry) for determining the convergence, default=2e-6 Ry +of_tolp 1e-05 #tolerance of potential for determining the convergence, default=1e-5 in a.u. +of_tf_weight 1 #weight of TF KEDF +of_vw_weight 1 #weight of vW KEDF +of_wt_alpha 0.833333 #parameter alpha of WT KEDF +of_wt_beta 0.833333 #parameter beta of WT KEDF +of_wt_rho0 1 #the average density of system, used in WT KEDF, in Bohr^-3 +of_hold_rho0 0 #If set to 1, the rho0 will be fixed even if the volume of system has changed, it will be set to 1 automaticly if of_wt_rho0 is not zero +of_full_pw 0 #If set to 1, ecut will be ignored when collect planewaves, so that all planewaves will be used +of_full_pw_dim 0 #If of_full_pw = true, dimention of FFT is testricted to be (0) either odd or even; (1) odd only; (2) even only +of_read_kernel 0 #If set to 1, the kernel of WT KEDF will be filled from file of_kernel_file, not from formula. Only usable for WT KEDF +of_kernel_file WTkernel.txt #The name of WT kernel file. + +#Parameters (19.dft+u) +dft_plus_u 0 #true:DFT+U correction; false: standard DFT calcullation(default) +yukawa_lambda -1 #default:0.0 +yukawa_potential 0 #default: false +omc 0 #the mode of occupation matrix control +hubbard_u 0 #Hubbard Coulomb interaction parameter U(ev) +orbital_corr -1 #which correlated orbitals need corrected ; d:2 ,f:3, do not need correction:-1 diff --git a/source/module_io/test/write_input_test.cpp b/source/module_io/test/write_input_test.cpp index e5c3a45106..e26ea020c5 100644 --- a/source/module_io/test/write_input_test.cpp +++ b/source/module_io/test/write_input_test.cpp @@ -294,9 +294,6 @@ TEST_F(write_input,print) EXPECT_THAT(output,testing::HasSubstr("")); EXPECT_THAT(output,testing::HasSubstr("#Parameters (16.tddft)")); EXPECT_THAT(output,testing::HasSubstr("td_force_dt 0.02 #time of force change")); - EXPECT_THAT(output,testing::HasSubstr("td_val_elec_01 1 #td_val_elec_01")); - EXPECT_THAT(output,testing::HasSubstr("td_val_elec_02 1 #td_val_elec_02")); - EXPECT_THAT(output,testing::HasSubstr("td_val_elec_03 1 #td_val_elec_03")); EXPECT_THAT(output,testing::HasSubstr("td_vext 0 #add extern potential or not")); EXPECT_THAT(output,testing::HasSubstr("td_vext_dire 1 #extern potential direction")); EXPECT_THAT(output,testing::HasSubstr("out_dipole 0 #output dipole or not")); diff --git a/source/module_io/write_dipole.cpp b/source/module_io/write_dipole.cpp index 2e34bc29fe..1b9dc884af 100644 --- a/source/module_io/write_dipole.cpp +++ b/source/module_io/write_dipole.cpp @@ -28,6 +28,12 @@ void ModuleIO::write_dipole(const double *rho_save, } } + double bmod[3]; + for (int i = 0; i < 3; i++) + { + bmod[i] = prepare(GlobalC::ucell, i); + } + #ifndef __MPI double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; for (int k = 0; k < GlobalC::rhopw->nz; k++) @@ -54,261 +60,73 @@ void ModuleIO::write_dipole(const double *rho_save, ofs << istep << " " << dipole_elec_x << " " << dipole_elec_y << dipole_elec_z; #else - double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; - // only do in the first pool. - if (GlobalV::MY_POOL == 0) - { - // num_z: how many planes on processor 'ip' - int *num_z = new int[GlobalV::NPROC_IN_POOL]; - ModuleBase::GlobalFunc::ZEROS(num_z, GlobalV::NPROC_IN_POOL); - for (int iz = 0; iz < GlobalC::bigpw->nbz; iz++) - { - int ip = iz % GlobalV::NPROC_IN_POOL; - num_z[ip] += GlobalC::bigpw->bz; - } + double dipole_elec[3] = {0.0, 0.0, 0.0}; - // start_z: start position of z in - // processor ip. - int *start_z = new int[GlobalV::NPROC_IN_POOL]; - ModuleBase::GlobalFunc::ZEROS(start_z, GlobalV::NPROC_IN_POOL); - for (int ip = 1; ip < GlobalV::NPROC_IN_POOL; ip++) - { - start_z[ip] = start_z[ip - 1] + num_z[ip - 1]; - } - - // which_ip: found iz belongs to which ip. - int *which_ip = new int[GlobalC::rhopw->nz]; - ModuleBase::GlobalFunc::ZEROS(which_ip, GlobalC::rhopw->nz); - for (int iz = 0; iz < GlobalC::rhopw->nz; iz++) - { - for (int ip = 0; ip < GlobalV::NPROC_IN_POOL; ip++) - { - if (iz >= start_z[GlobalV::NPROC_IN_POOL - 1]) - { - which_ip[iz] = GlobalV::NPROC_IN_POOL - 1; - break; - } - else if (iz >= start_z[ip] && iz < start_z[ip + 1]) - { - which_ip[iz] = ip; - break; - } - } - // GlobalV::ofs_running << "\n iz=" << iz << " ip=" << which_ip[iz]; - } - - // int count=0; - int nxy = GlobalC::rhopw->nx * GlobalC::rhopw->ny; - double *zpiece = new double[nxy]; - - // save the rho one z by one z. - for (int iz = 0; iz < GlobalC::rhopw->nz; iz++) - { - // std::cout << "\n iz=" << iz << std::endl; - // tag must be different for different iz. - ModuleBase::GlobalFunc::ZEROS(zpiece, nxy); - int tag = iz; - MPI_Status ierror; - - // case 1: the first part of rho in processor 0. - if (which_ip[iz] == 0 && GlobalV::RANK_IN_POOL == 0) - { - for (int ir = 0; ir < nxy; ir++) - { - // mohan change to rho_save on 2012-02-10 - // because this can make our next restart calculation lead - // to the same scf_thr as the one saved. - zpiece[ir] = rho_save[ir * GlobalC::rhopw->nplane + iz - GlobalC::rhopw->startz_current]; - // GlobalV::ofs_running << "\n get zpiece[" << ir << "]=" << zpiece[ir] << " - // ir*GlobalC::pw.nczp+iz=" << ir*GlobalC::pw.nczp+iz; - } - } - // case 2: > first part rho: send the rho to - // processor 0. - else if (which_ip[iz] == GlobalV::RANK_IN_POOL) - { - for (int ir = 0; ir < nxy; ir++) - { - // zpiece[ir] = rho[is][ir*num_z[GlobalV::RANK_IN_POOL]+iz]; - zpiece[ir] = rho_save[ir * GlobalC::rhopw->nplane + iz - GlobalC::rhopw->startz_current]; - // GlobalV::ofs_running << "\n get zpiece[" << ir << "]=" << zpiece[ir] << " - // ir*GlobalC::pw.nczp+iz=" << ir*GlobalC::pw.nczp+iz; - } - MPI_Send(zpiece, nxy, MPI_DOUBLE, 0, tag, POOL_WORLD); - } - - // case 2: > first part rho: processor 0 receive the rho - // from other processors - else if (GlobalV::RANK_IN_POOL == 0) - { - MPI_Recv(zpiece, nxy, MPI_DOUBLE, which_ip[iz], tag, POOL_WORLD, &ierror); - // GlobalV::ofs_running << "\n Receieve First number = " << zpiece[0]; - } - - // write data - if (GlobalV::MY_RANK == 0) - { - // ofs << "\niz=" << iz; - // mohan update 2011-03-30 - for (int iy = 0; iy < GlobalC::rhopw->ny; iy++) - { - for (int ix = 0; ix < GlobalC::rhopw->nx; ix++) - { - /* - if(ixny + iy] * ix * GlobalC::ucell.lat0 * 0.529177 - / GlobalC::rhopw->nx; - dipole_elec_y += zpiece[ix * GlobalC::rhopw->ny + iy] * iy * GlobalC::ucell.lat0 * 0.529177 - / GlobalC::rhopw->ny; - dipole_elec_z += zpiece[ix * GlobalC::rhopw->ny + iy] * iz * GlobalC::ucell.lat0 * 0.529177 - / GlobalC::rhopw->nz; - } - } - } - } // end iz + for (int ir = 0; ir < GlobalC::rhopw->nrxx; ++ir) + { + int i = ir / (GlobalC::rhopw->ny * GlobalC::rhopw->nplane); + int j = ir / GlobalC::rhopw->nplane - i * GlobalC::rhopw->ny; + int k = ir % GlobalC::rhopw->nplane + GlobalC::rhopw->startz_current; + double x = (double)i / GlobalC::rhopw->nx; + double y = (double)j / GlobalC::rhopw->ny; + double z = (double)k / GlobalC::rhopw->nz; + + dipole_elec[0] += rho_save[ir] * x; + dipole_elec[1] += rho_save[ir] * y; + dipole_elec[2] += rho_save[ir] * z; + } - delete[] zpiece; + Parallel_Reduce::reduce_double_pool(dipole_elec[0]); + Parallel_Reduce::reduce_double_pool(dipole_elec[1]); + Parallel_Reduce::reduce_double_pool(dipole_elec[2]); + for (int i = 0; i < 3; ++i) + { + dipole_elec[i] *= GlobalC::ucell.lat0 / bmod[i] * GlobalC::ucell.omega / GlobalC::rhopw->nxyz; + } - double nxyz = GlobalC::rhopw->nx * GlobalC::rhopw->ny * GlobalC::rhopw->nz; - dipole_elec_x *= GlobalC::ucell.omega / static_cast(nxyz); - dipole_elec_y *= GlobalC::ucell.omega / static_cast(nxyz); - dipole_elec_z *= GlobalC::ucell.omega / static_cast(nxyz); + std::cout << std::setprecision(8) << "dipole_elec_x: " << dipole_elec[0] << std::endl; + std::cout << std::setprecision(8) << "dipole_elec_y: " << dipole_elec[1] << std::endl; + std::cout << std::setprecision(8) << "dipole_elec_z: " << dipole_elec[2] << std::endl; - std::cout << std::setprecision(8) << "dipole_elec_x: " << dipole_elec_x << std::endl; - std::cout << std::setprecision(8) << "dipole_elec_y: " << dipole_elec_y << std::endl; - std::cout << std::setprecision(8) << "dipole_elec_z: " << dipole_elec_z << std::endl; + ofs << std::setprecision(8) << istep << " " << dipole_elec[0] << " " << dipole_elec[1] << " " << dipole_elec[2] + << std::endl; - ofs << istep << " " << dipole_elec_x << " " << dipole_elec_y << " " << dipole_elec_z << std::endl; + double dipole_ion[3] = {0.0}; + double dipole_sum = 0.0; - /* - double dipole_ion_x = 0.0, dipole_ion_y = 0.0, dipole_ion_z = 0.0, dipole_sum = 0.0; - if (GlobalC::ucell.ntype == 1) - { - for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - { - dipole_ion_x += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - dipole_ion_y += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - dipole_ion_z += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - } - } - else if (GlobalC::ucell.ntype == 2) - { - for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - { - dipole_ion_x += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - dipole_ion_y += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - dipole_ion_z += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - } - for (int ia = 0; ia < GlobalC::ucell.atoms[1].na; ia++) - { - dipole_ion_x += GlobalC::ucell.atoms[1].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_02; - dipole_ion_y += GlobalC::ucell.atoms[1].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_02; - dipole_ion_z += GlobalC::ucell.atoms[1].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_02; - } - } - else if (GlobalC::ucell.ntype == 3) + for (int i = 0; i < 3; ++i) + { + for (int it = 0; it < GlobalC::ucell.ntype; ++it) { - for (int ia = 0; ia < GlobalC::ucell.atoms[0].na; ia++) - { - dipole_ion_x += GlobalC::ucell.atoms[0].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - dipole_ion_y += GlobalC::ucell.atoms[0].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - dipole_ion_z += GlobalC::ucell.atoms[0].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_01; - } - for (int ia = 0; ia < GlobalC::ucell.atoms[1].na; ia++) - { - dipole_ion_x += GlobalC::ucell.atoms[1].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_02; - dipole_ion_y += GlobalC::ucell.atoms[1].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_02; - dipole_ion_z += GlobalC::ucell.atoms[1].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_02; - } - for (int ia = 0; ia < GlobalC::ucell.atoms[2].na; ia++) + double sum = 0; + for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ++ia) { - dipole_ion_x += GlobalC::ucell.atoms[2].taud[ia].x * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_03; - dipole_ion_y += GlobalC::ucell.atoms[2].taud[ia].y * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_03; - dipole_ion_z += GlobalC::ucell.atoms[2].taud[ia].z * GlobalC::ucell.lat0 * 0.529177 - * ELEC_evolve::td_val_elec_03; + + sum += GlobalC::ucell.atoms[it].taud[ia][i]; } + dipole_ion[i] += sum * GlobalC::ucell.atoms[it].ncpp.zv; } - else - { - std::cout << "atom ntype is too large!" << std::endl; - } - for(int it=1; it<(GlobalC::ucell.ntype); it++) - { - for(int ia=0; ia None: + self.dipolefile = dipolefile + self._indices, self.dipole_x, self.dipole_y, self.dipole_z = self.read( + self.dipolefile) + self.dt = dt # in fs + self.time = self._indices*dt + self.energies = Freq2eV*self._indices/dt/len(self._indices) + + @classmethod + def read(cls, filename: PathLike): + """Read dipole data file + + :params filename: string of dipole data file + """ + + data = np.loadtxt(filename, dtype=float) + indices = data[:, 0] + x = data[:, 1] + y = data[:, 2] + z = data[:, 3] + + return indices, x, y, z + + @property + def dipole_data(self): + return {0: self.dipole_x, 1: self.dipole_y, 2: self.dipole_z} + + def plot_dipole(self, fig: Figure, ax: axes.Axes, directions: list = [0, 1, 2], time_range: list = []): + """Plot dipole data in x,y,z directions + + :params fig: (matplotlib.figure.Figure) + :params ax: (matplotlib.axes.Axes) + :params directions: (list) 0->X, 1->Y, 2->Z + :params time_range: (list) range of time (in unit fs) to plot + """ + + labels = {0: 'X', 1: 'Y', 2: 'Z'} + + for direc in directions: + ax.plot(self.time, self.dipole_data[direc], label=labels[direc]) + + ax.set_xlabel('Times (fs)') + ax.set_ylabel('Dipole') + ax.legend() + if time_range: + ax.set_xlim(time_range) + + return fig, ax + + @property + def alpha_x(self): + return np.fft.fft(self.dipole_x) + + @property + def alpha_y(self): + return np.fft.fft(self.dipole_y) + + @property + def alpha_z(self): + return np.fft.fft(self.dipole_z) + + @property + def alpha_data(self): + return {0: self.alpha_x, 1: self.alpha_y, 2: self.alpha_z} + + def get_abs(self, direc: int): + S = np.abs(self.alpha_data[direc].imag) + return S + + def plot_abs(self, fig: Figure, ax: axes.Axes, directions: list = [0, 1, 2], x_range: list = [], unit: str = 'eV'): + """Plot Absportion Spectrum under Delta light field in x,y,z directions + + :params fig: (matplotlib.figure.Figure) + :params ax: (matplotlib.axes.Axes) + :params directions: (list) 0->X, 1->Y, 2->Z + :params x_range: (list) range of energies (in unit eV) to plot + :params unit: (str) + """ + + assert unit in ['eV', 'nm'] + labels = {0: 'X', 1: 'Y', 2: 'Z'} + x_data = self.energies if unit == 'eV' else sc.nu2lambda( + sc.eV/sc.h*self.energies)*1e9 + + #plot the adsorption spectra and output the data + adsorption_spectra_data = x_data[:, np.newaxis] + for direc in directions: + ax.plot(x_data, self.get_abs(direc), label=labels[direc]) + adsorption_spectra_data = np.concatenate((adsorption_spectra_data, self.get_abs(direc)[:, np.newaxis]),axis=1) + if direc != directions[-1]: + adsorption_spectra_data = np.concatenate((adsorption_spectra_data, x_data[:, np.newaxis]),axis=1) + np.savetxt('absorpation_spectra.dat', adsorption_spectra_data) + + xlabel = 'Energy (eV)' if unit == 'eV' else 'Wave Length (nm)' + ax.set_xlabel(xlabel) + ax.set_ylabel('Absportion') + ax.legend() + if x_range: + ax.set_xlim(x_range) + + return fig, ax + + +if __name__ == "__main__": + dipolefile = './SPIN1_DIPOLE' + dipole = Dipole(dipolefile, dt=0.0034) + + import matplotlib.pyplot as plt + fig1, ax1 = plt.subplots() + fig1, ax1 = dipole.plot_dipole(fig1, ax1) + fig1.savefig('dipole.png') + + fig2, ax2 = plt.subplots() + x_range = [0, 400] + fig2, ax2 = dipole.plot_abs( + fig2, ax2, x_range=x_range, unit='nm') + fig2.savefig('abs.png') diff --git a/tools/plot-tools/abacus_plot/utils.py b/tools/plot-tools/abacus_plot/utils.py index 7fd3679f1c..55d6ec6804 100644 --- a/tools/plot-tools/abacus_plot/utils.py +++ b/tools/plot-tools/abacus_plot/utils.py @@ -22,6 +22,7 @@ def remove_empty(a: list) -> list: while [] in a: a.remove([]) + def handle_data(data): data.remove('') @@ -31,6 +32,7 @@ def handle_elem(elem): return elist return list(map(handle_elem, data)) + def skip_notes(line: str) -> str: """Delete comments lines with '#' or '//' @@ -332,9 +334,10 @@ def parse_projected_data(orbitals, species: Union[Sequence[Any], Dict[Any, List[ return data, totnum + def key2int(species): """Convert keys of dict in str type to int""" - + new_species = {} if isinstance(species, dict): elements = list(map(int, species.keys())) @@ -352,5 +355,5 @@ def key2int(species): new_species[elem] = [] for l_index in l[i]: new_species[elem].append(int(l_index)) - - return new_species \ No newline at end of file + + return new_species