diff --git a/source/Makefile.Objects b/source/Makefile.Objects index d3761489ed..5c094792f2 100644 --- a/source/Makefile.Objects +++ b/source/Makefile.Objects @@ -382,14 +382,19 @@ OBJS_IO_LCAO=cal_r_overlap_R.o\ OBJS_LCAO=DM_gamma.o\ DM_k.o\ - ELEC_evolve.o\ + evolve_elec.o\ + evolve_psi.o\ + bandenergy.o\ + middle_hamilt.o\ + norm_psi.o\ + propagator.o\ + upsi.o\ FORCE_STRESS.o\ FORCE_gamma.o\ FORCE_gamma_edm.o\ FORCE_gamma_tvnl.o\ FORCE_gamma_vl.o\ FORCE_k.o\ - LCAO_evolve.o\ LCAO_gen_fixedH.o\ LCAO_hamilt.o\ LCAO_matrix.o\ diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp index a915af5b7b..c1b67cc64d 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp +++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp @@ -2,7 +2,7 @@ #include "module_base/constants.h" #include "module_base/timer.h" -#include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" +#include "module_hamilt_lcao/module_tddft/evolve_elec.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/input.h" #include "module_io/input_conv.h" @@ -11,11 +11,61 @@ namespace elecstate { int H_TDDFT_pw::istep = -1; -//========================================================== -// this function aims to add external time-dependent potential -// (eg: linear potential) used in tddft -// fuxiang add in 2017-05 -//========================================================== + +double H_TDDFT_pw::amp; +double H_TDDFT_pw::bmod; +double H_TDDFT_pw::bvec[3]; + +int H_TDDFT_pw::stype; // 0 : length gauge 1: velocity gauge + +std::vector H_TDDFT_pw::ttype; +// 0 Gauss type function. +// 1 trapezoid type function. +// 2 Trigonometric functions, sin^2. +// 3 heaviside function. +// 4 HHG function. + +int H_TDDFT_pw::tstart; +int H_TDDFT_pw::tend; +double H_TDDFT_pw::dt; + +// space domain parameters + +// length gauge +double H_TDDFT_pw::lcut1; +double H_TDDFT_pw::lcut2; + +// time domain parameters + +// Gauss +int H_TDDFT_pw::gauss_count; +std::vector H_TDDFT_pw::gauss_omega; // time(a.u.)^-1 +std::vector H_TDDFT_pw::gauss_phase; +std::vector H_TDDFT_pw::gauss_sigma; // time(a.u.) +std::vector H_TDDFT_pw::gauss_t0; +std::vector H_TDDFT_pw::gauss_amp; // Ry/bohr + +// trapezoid +int H_TDDFT_pw::trape_count; +std::vector H_TDDFT_pw::trape_omega; // time(a.u.)^-1 +std::vector H_TDDFT_pw::trape_phase; +std::vector H_TDDFT_pw::trape_t1; +std::vector H_TDDFT_pw::trape_t2; +std::vector H_TDDFT_pw::trape_t3; +std::vector H_TDDFT_pw::trape_amp; // Ry/bohr + +// Trigonometric +int H_TDDFT_pw::trigo_count; +std::vector H_TDDFT_pw::trigo_omega1; // time(a.u.)^-1 +std::vector H_TDDFT_pw::trigo_omega2; // time(a.u.)^-1 +std::vector H_TDDFT_pw::trigo_phase1; +std::vector H_TDDFT_pw::trigo_phase2; +std::vector H_TDDFT_pw::trigo_amp; // Ry/bohr + +// Heaviside +int H_TDDFT_pw::heavi_count; +std::vector H_TDDFT_pw::heavi_t0; +std::vector H_TDDFT_pw::heavi_amp; // Ry/bohr void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) { @@ -24,10 +74,8 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) // time evolve H_TDDFT_pw::istep++; - read_parameters(&INPUT); - // judgement to skip vext - if (!ELEC_evolve::td_vext || istep > tend || istep < tstart) + if (!module_tddft::Evolve_elec::td_vext || istep > tend || istep < tstart) { return; } @@ -36,13 +84,17 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) ModuleBase::timer::tick("H_TDDFT_pw", "cal_fixed_v"); int count = 0; + gauss_count = 0; + trape_count = 0; + trigo_count = 0; + heavi_count = 0; - for (auto direc: ELEC_evolve::td_vext_dire_case) + for (auto direc: module_tddft::Evolve_elec::td_vext_dire_case) { std::vector vext_space(this->rho_basis_->nrxx, 0.0); double vext_time = cal_v_time(ttype[count]); - if (ELEC_evolve::out_efield && GlobalV::MY_RANK == 0) + if (module_tddft::Evolve_elec::out_efield && GlobalV::MY_RANK == 0) { std::stringstream as; as << GlobalV::global_out_dir << "efield_" << count << ".dat"; @@ -62,80 +114,6 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo) return; } -std::vector H_TDDFT_pw::set_parameters(std::string params, double c) -{ - std::vector params_ori; - std::vector params_out; - Input_Conv::parse_expression(params, params_ori); - for (auto param: params_ori) - params_out.emplace_back(param * c); - - return params_out; -} - -void H_TDDFT_pw::read_parameters(Input *in) -{ - stype = in->td_stype; - - Input_Conv::parse_expression(in->td_ttype, ttype); - - tstart = in->td_tstart; - tend = in->td_tend; - - dt = in->mdp.md_dt; - - // space domain parameters - - // length gauge - lcut1 = in->td_lcut1; - lcut2 = in->td_lcut2; - - // time domain parameters - - // Gauss - 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_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 - - // trapezoid - trape_count = 0; - trape_omega = set_parameters(in->td_trape_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 - trape_phase = set_parameters(in->td_trape_phase, 1.0); - trape_t1 = set_parameters(in->td_trape_t1, 1.0); - trape_t2 = set_parameters(in->td_trape_t2, 1.0); - trape_t3 = set_parameters(in->td_trape_t3, 1.0); - trape_amp = set_parameters(in->td_trape_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr - - // Trigonometric - trigo_count = 0; - trigo_omega1 = set_parameters(in->td_trigo_freq1, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 - trigo_omega2 = set_parameters(in->td_trigo_freq2, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 - trigo_phase1 = set_parameters(in->td_trigo_phase1, 1.0); - trigo_phase2 = set_parameters(in->td_trigo_phase2, 1.0); - trigo_amp = set_parameters(in->td_trigo_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr - - // Heaviside - heavi_count = 0; - heavi_t0 = set_parameters(in->td_heavi_t0, 1.0); - heavi_amp = set_parameters(in->td_heavi_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr - - // HHG - // hhg_count = 0; - // hhg_amp1 = set_parameters(in->td_hhg_amp1, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr - // hhg_amp2 = set_parameters(in->td_hhg_amp2, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr - // hhg_omega1 = set_parameters(in->td_hhg_freq1, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 - // hhg_omega2 = set_parameters(in->td_hhg_freq2, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 - // hhg_phase1 = set_parameters(in->td_hhg_phase1, 1.0); - // hhg_phase2 = set_parameters(in->td_hhg_phase2, 1.0); - // hhg_t0 = set_parameters(in->td_hhg_t0, 1.0); - // hhg_sigma = set_parameters(in->td_hhg_sigma, 1/ModuleBase::AU_to_FS); - - return; -} - void H_TDDFT_pw::cal_v_space(std::vector &vext_space, int direc) { ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space"); @@ -165,11 +143,7 @@ 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); - } + prepare(GlobalC::ucell, direc); for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir) { @@ -183,15 +157,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) / bmod[0]; + vext_space[ir] = cal_v_space_length_potential(x) / bmod; break; case 2: - vext_space[ir] = cal_v_space_length_potential(y) / bmod[1]; + vext_space[ir] = cal_v_space_length_potential(y) / bmod; break; case 3: - vext_space[ir] = cal_v_space_length_potential(z) / bmod[2]; + vext_space[ir] = cal_v_space_length_potential(z) / bmod; break; default: @@ -336,30 +310,8 @@ double H_TDDFT_pw::cal_v_time_heaviside() return vext_time; } -// double H_TDDFT_pw::cal_v_time_HHG() -// { -// double vext_time = 0.0; -// double t0 = *(hhg_t0.begin() + hhg_count); -// double omega1 = *(hhg_omega1.begin() + hhg_count); -// double phase1 = *(hhg_phase1.begin() + hhg_count); -// double omega2 = *(hhg_omega2.begin() + hhg_count); -// double phase2 = *(hhg_phase2.begin() + hhg_count); -// double amp1 = *(hhg_amp1.begin() + hhg_count); -// double amp2 = *(hhg_amp2.begin() + hhg_count); -// double sigma = *(trigo_amp2.begin() + trigo_count); - -// double hhg_t = (istep - t0) * dt; -// vext_time = (cos(omega1 * hhg_t + phase1) * amp1 + cos(omega2 * hhg_t + phase2) * amp2) -// * exp(-hhg_t * hhg_t * 0.5 / (sigma * sigma)); -// hhg_count++; - -// return vext_time; -// } - -double H_TDDFT_pw::prepare(const UnitCell &cell, int &dir) +void 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; @@ -383,7 +335,22 @@ double H_TDDFT_pw::prepare(const UnitCell &cell, int &dir) 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; +} + +void H_TDDFT_pw ::compute_force(const UnitCell& cell, ModuleBase::matrix& fe) +{ + int iat = 0; + for (int it = 0; it < cell.ntype; ++it) + { + for (int ia = 0; ia < cell.atoms[it].na; ++ia) + { + for (int jj = 0; jj < 3; ++jj) + { + fe(iat, jj) = ModuleBase::e2 * amp * cell.atoms[it].ncpp.zv * bvec[jj] / bmod; + } + ++iat; + } + } } } // namespace elecstate diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h index e600219b24..15b97f9afd 100644 --- a/source/module_elecstate/potentials/H_TDDFT_pw.h +++ b/source/module_elecstate/potentials/H_TDDFT_pw.h @@ -1,8 +1,9 @@ #ifndef H_TDDFT_PW_H #define H_TDDFT_PW_H -#include "pot_base.h" #include "module_io/input.h" +#include "module_io/input_conv.h" +#include "pot_base.h" namespace elecstate { @@ -21,80 +22,79 @@ class H_TDDFT_pw : public PotBase void cal_fixed_v(double* vl_pseudo) override; - private: - // internal time-step, - //-------hypothesis------- - // Vext will evolve by time, every time cal_fixed_v() is called, istep++ - //------------------------ - static int istep; + /** + * @brief compute force of electric field + * + * @param[in] cell information of cell + * @param[out] fe force of electric field F=qE + */ + static void compute_force(const UnitCell& cell, ModuleBase::matrix& fe); // parameters - int stype ; // 0 : length gauge 1: velocity gauge + static int stype; // 0 : length gauge 1: velocity gauge - std::vector ttype ; + static std::vector ttype; // 0 Gauss type function. // 1 trapezoid type function. // 2 Trigonometric functions, sin^2. // 3 heaviside function. // 4 HHG function. - int tstart; - int tend; - double dt; + static int tstart; + static int tend; + static double dt; // space domain parameters //length gauge - double lcut1; - double lcut2; + static double lcut1; + static double lcut2; // time domain parameters // Gauss - int gauss_count; - std::vector gauss_omega; // time(a.u.)^-1 - std::vector gauss_phase ; - std::vector gauss_sigma ; // time(a.u.) - std::vector gauss_t0 ; - std::vector gauss_amp ; // Ry/bohr + static int gauss_count; + static std::vector gauss_omega; // time(a.u.)^-1 + static std::vector gauss_phase; + static std::vector gauss_sigma; // time(a.u.) + static std::vector gauss_t0; + static std::vector gauss_amp; // Ry/bohr // trapezoid - int trape_count; - std::vector trape_omega ; // time(a.u.)^-1 - std::vector trape_phase ; - std::vector trape_t1 ; - std::vector trape_t2 ; - std::vector trape_t3 ; - std::vector trape_amp ; // Ry/bohr + static int trape_count; + static std::vector trape_omega; // time(a.u.)^-1 + static std::vector trape_phase; + static std::vector trape_t1; + static std::vector trape_t2; + static std::vector trape_t3; + static std::vector trape_amp; // Ry/bohr // Trigonometric - int trigo_count; - std::vector trigo_omega1 ; // time(a.u.)^-1 - std::vector trigo_omega2 ; // time(a.u.)^-1 - std::vector trigo_phase1 ; - std::vector trigo_phase2 ; - std::vector trigo_amp ; // Ry/bohr + static int trigo_count; + static std::vector trigo_omega1; // time(a.u.)^-1 + static std::vector trigo_omega2; // time(a.u.)^-1 + static std::vector trigo_phase1; + static std::vector trigo_phase2; + static std::vector trigo_amp; // Ry/bohr // Heaviside - int heavi_count; - std::vector heavi_t0; - std::vector heavi_amp; // Ry/bohr - - // HHG - // int hhg_count = 0; - // std::vector hhg_amp1; // Ry/bohr - // std::vector hhg_amp2; // Ry/bohr - // std::vector hhg_omega1; // time(a.u.)^-1 - // std::vector hhg_omega2; // time(a.u.)^-1 - // std::vector hhg_phase1; - // std::vector hhg_phase2; - // std::vector hhg_t0; - // std::vector hhg_sigma; // time(a.u.) + static int heavi_count; + static std::vector heavi_t0; + static std::vector heavi_amp; // Ry/bohr - const UnitCell* ucell_ = nullptr; + private: + // internal time-step, + //-------hypothesis------- + // Vext will evolve by time, every time cal_fixed_v() is called, istep++ + //------------------------ + static int istep; + + static double amp; - std::vector set_parameters(std::string params, double c); - void read_parameters(Input* in); + static double bmod; + static double bvec[3]; + + const UnitCell* ucell_ = nullptr; // potential of electric field in space domain : length gauge and velocity gauge void cal_v_space(std::vector &vext_space, int direc); @@ -110,7 +110,7 @@ class H_TDDFT_pw : public PotBase double cal_v_time_heaviside(); // double cal_v_time_HHG(); - double prepare(const UnitCell &cell, int &dir); + void 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 6901233ed7..5e518c27dc 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp @@ -13,7 +13,7 @@ #include "module_base/scalapack_connector.h" #include "module_elecstate/module_charge/symmetry_rho.h" #include "module_elecstate/occupy.h" -#include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" +#include "module_hamilt_lcao/module_tddft/evolve_elec.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/print_info.h" @@ -114,75 +114,6 @@ void ESolver_KS_LCAO_TDDFT::Init(Input& inp, UnitCell& ucell) this->pelec_td = dynamic_cast(this->pelec); } -void ESolver_KS_LCAO_TDDFT::eachiterinit(const int istep, const int iter) -{ - // mohan add 2010-07-16 - // used for pulay mixing. - if (iter == 1) - GlobalC::CHR_MIX.reset(); - - // mohan update 2012-06-05 - this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(); - - if (GlobalC::wf.init_wfc == "file") - { - if (iter == 1) - { - std::cout << " WAVEFUN -> CHARGE " << std::endl; - - // The occupation should be read in together. - // Occupy::calculate_weights(); //mohan add 2012-02-15 - - // calculate the density matrix using read in wave functions - // and the ncalculate the charge density on grid. - if (this->psi != nullptr) - { - if (istep >= 2) - { - this->pelec_td->psiToRho_td(this->psi[0]); - } - else - { - this->pelec_td->psiToRho(this->psi[0]); - } - } - else - { - this->pelec_td->psiToRho(this->psid[0]); - } - - // calculate the local potential(rho) again. - // the grid integration will do in later grid integration. - - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // a puzzle remains here. - // if I don't renew potential, - // The scf_thr is very small. - // OneElectron, Hartree and - // Exc energy are all correct - // except the band energy. - // - // solved by mohan 2010-09-10 - // there are there rho here: - // rho1: formed by read in orbitals. - // rho2: atomic rho, used to construct H - // rho3: generated by after diagonalize - // here converged because rho3 and rho1 - // are very close. - // so be careful here, make sure - // rho1 and rho2 are the same rho. - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - this->pelec->pot->init_pot(istep, this->pelec->charge); - this->pelec->f_en.descf = this->pelec->cal_delta_escf(); - } - } - - if (!GlobalV::GAMMA_ONLY_LOCAL) - { - this->UHM.GK.renew(); - } -} - void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr) { @@ -190,18 +121,19 @@ 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->Hk_laststep, - this->pelec_td->ekb, - td_htype, - INPUT.propagator, - GlobalC::kv.nks); + module_tddft::Evolve_elec::solve_psi(istep, + GlobalV::NBANDS, + GlobalV::NLOCAL, + this->p_hamilt, + this->LOWF, + this->psi, + this->psi_laststep, + this->Hk_laststep, + this->pelec_td->ekb, + td_htype, + INPUT.propagator, + GlobalC::kv.nks); this->pelec_td->psiToRho_td(this->psi[0]); - // this->pelec_td->psiToRho(this->psi[0]); } // using HSolverLCAO::solve() else if (this->phsol != nullptr) @@ -223,6 +155,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr) ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "HSolver has not been initialed!"); } + // print occupation of each band if (iter == 1 && istep <= 2) { GlobalV::ofs_running @@ -246,19 +179,15 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr) << endl; } - // (3) sum bands to calculate charge density - // if (istep <= 1 ) Occupy::calculate_weights(); - for (int ik = 0; ik < GlobalC::kv.nks; ++ik) { this->pelec_td->print_band(ik, INPUT.printe, iter); } - // (4) mohan add 2010-06-24 // using new charge density. this->pelec->cal_energies(1); - // (5) symmetrize the charge density + // symmetrize the charge density only for ground state if (istep <= 1) { Symmetry_rho srho; @@ -269,7 +198,10 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(int istep, int iter, double ethr) } // (6) compute magnetization, only for spin==2 - GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, this->pelec->charge->nxyz, this->pelec->charge->rho, pelec->nelec_spin.data()); + GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx, + this->pelec->charge->nxyz, + this->pelec->charge->rho, + pelec->nelec_spin.data()); // (7) calculate delta energy this->pelec->f_en.deband = this->pelec->cal_delta_eband(); @@ -341,8 +273,7 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) elecstate::ElecStateLCAO::out_wfc_flag = 0; } - // (9) Calculate new potential according to new Charge Density. - + // Calculate new potential according to new Charge Density if (!this->conv_elec) { if (GlobalV::NSPIN == 4) @@ -400,10 +331,12 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) } } - if (istep > 1 && ELEC_evolve::td_edm == 0) + // calculate energy density matrix for tddft + if (istep > 1 && module_tddft::Evolve_elec::td_edm == 0) this->cal_edm_tddft(); } + // print "eigen value" for tddft if (this->conv_elec) { GlobalV::ofs_running @@ -430,156 +363,17 @@ void ESolver_KS_LCAO_TDDFT::updatepot(const int istep, const int iter) void ESolver_KS_LCAO_TDDFT::afterscf(const int istep) { - // if (this->conv_elec || iter == GlobalV::SCF_NMAX) - // { - //-------------------------------------- - // 1. output charge density for converged, - // 0 means don't need to consider iter, - //-------------------------------------- - for (int is = 0; is < GlobalV::NSPIN; is++) { - const int precision = 3; - std::stringstream ssc; - ssc << GlobalV::global_out_dir << "SPIN" << is + 1 << "_CHG.cube"; - const double ef_tmp = this->pelec->eferm.get_efval(is); - ModuleIO::write_rho( -#ifdef __MPI - GlobalC::bigpw->bz, - GlobalC::bigpw->nbz, - GlobalC::rhopw->nplane, - GlobalC::rhopw->startz_current, -#endif - pelec->charge->rho_save[is], - is, - GlobalV::NSPIN, - 0, - ssc.str(), - GlobalC::rhopw->nx, - GlobalC::rhopw->ny, - GlobalC::rhopw->nz, - ef_tmp, - &(GlobalC::ucell), - precision); - - if (ELEC_evolve::out_dipole == 1) + if (module_tddft::Evolve_elec::out_dipole == 1) { std::stringstream ss_dipole; ss_dipole << GlobalV::global_out_dir << "SPIN" << is + 1 << "_DIPOLE"; ModuleIO::write_dipole(pelec->charge->rho_save[is], pelec->charge->rhopw, is, istep, ss_dipole.str()); } - - std::stringstream ssd; - if (GlobalV::GAMMA_ONLY_LOCAL) - { - ssd << GlobalV::global_out_dir << "SPIN" << is + 1 << "_DM"; - } - else - { - ssd << GlobalV::global_out_dir << "SPIN" << is + 1 << "_DM_R"; - } - - ModuleIO::write_dm( -#ifdef __MPI - GlobalC::GridT.trace_lo, -#endif - is, - 0, - ssd.str(), - precision, - this->LOC.out_dm, - this->LOC.DM, - ef_tmp, - &(GlobalC::ucell)); - - if (GlobalV::out_pot == 1) // LiuXh add 20200701 - { - std::stringstream ssp; - ssp << GlobalV::global_out_dir << "SPIN" << is + 1 << "_POT.cube"; - this->pelec->pot->write_potential( -#ifdef __MPI - GlobalC::bigpw->bz, - GlobalC::bigpw->nbz, - this->pw_rho->nplane, - this->pw_rho->startz_current, -#endif - is, - 0, - ssp.str(), - this->pw_rho->nx, - this->pw_rho->ny, - this->pw_rho->nz, - this->pelec->pot->get_effective_v(), - precision); - } - } - - if (this->conv_elec) - { - GlobalV::ofs_running << "\n charge density convergence is achieved" << std::endl; - GlobalV::ofs_running << " final etot is " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" - << std::endl; } - if (GlobalV::OUT_LEVEL != "m") - { - // this->pelec->print_eigenvalue(GlobalV::ofs_running); - } - - if (this->conv_elec) - { - // xiaohui add "OUT_LEVEL", 2015-09-16 - if (GlobalV::OUT_LEVEL != "m") - GlobalV::ofs_running << std::setprecision(16); - if (GlobalV::OUT_LEVEL != "m") - GlobalV::ofs_running << " EFERMI = " << this->pelec->eferm.ef * ModuleBase::Ry_to_eV << " eV" << std::endl; - } - else - { - GlobalV::ofs_running << " !! convergence has not been achieved @_@" << std::endl; - if (GlobalV::OUT_LEVEL == "ie" || GlobalV::OUT_LEVEL == "m") // xiaohui add "m" option, 2015-09-16 - std::cout << " !! CONVERGENCE HAS NOT BEEN ACHIEVED !!" << std::endl; - } - - if( GlobalV::CALCULATION != "md" || (istep % GlobalV::out_interval == 0)) - { - if (hsolver::HSolverLCAO::out_mat_hsR) - { - ModuleIO::output_HS_R(istep, - this->pelec->pot->get_effective_v(), - this->UHM, - GlobalC::kv); // LiuXh add 2019-07-15 - } - - if (hsolver::HSolverLCAO::out_mat_t) - { - ModuleIO::output_T_R(istep, this->UHM); // LiuXh add 2019-07-15 - } - - if (hsolver::HSolverLCAO::out_mat_dh) - { - ModuleIO::output_dH_R(istep, - this->pelec->pot->get_effective_v(), - this->UHM, - GlobalC::kv); // LiuXh add 2019-07-15 - } - - // add by jingan for out r_R matrix 2019.8.14 - if (INPUT.out_mat_r) - { - cal_r_overlap_R r_matrix; - r_matrix.init(*this->LOWF.ParaV); - - if (hsolver::HSolverLCAO::out_mat_hsR) - { - r_matrix.out_rR_other(istep, this->LM.output_R_coor); - } - else - { - r_matrix.out_rR(istep); - } - } - } + ESolver_KS_LCAO::afterscf(istep); } // use the original formula (Hamiltonian matrix) to calculate energy density matrix diff --git a/source/module_esolver/esolver_ks_lcao_tddft.h b/source/module_esolver/esolver_ks_lcao_tddft.h index 6df05bb38e..5121a2b224 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.h +++ b/source/module_esolver/esolver_ks_lcao_tddft.h @@ -29,7 +29,6 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO protected: virtual void hamilt2density(const int istep, const int iter, const double ethr) override; - virtual void eachiterinit(const int istep, const int iter) override; virtual void updatepot(const int istep, const int iter) override; virtual void afterscf(const int istep) override; void cal_edm_tddft(); diff --git a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp index 6c3a83f3b6..bff2d0fba8 100644 --- a/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp +++ b/source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp @@ -1,15 +1,16 @@ #include "FORCE_STRESS.h" -#include "module_hamilt_pw/hamilt_pwdft/global.h" + #include "module_hamilt_lcao/hamilt_lcaodft/global_fp.h" -#include "module_hamilt_lcao/module_dftu/dftu.h" //Quxin add for DFT+U on 20201029 +#include "module_hamilt_lcao/module_dftu/dftu.h" //Quxin add for DFT+U on 20201029 +#include "module_hamilt_pw/hamilt_pwdft/global.h" // new #include "module_base/timer.h" -#include "module_elecstate/potentials/efield.h" // liuyu add 2022-05-18 -#include "module_hamilt_general/module_surchem/surchem.h" //sunml add 2022-08-10 +#include "module_elecstate/potentials/efield.h" // liuyu add 2022-05-18 #include "module_elecstate/potentials/gatefield.h" // liuyu add 2022-09-13 +#include "module_hamilt_general/module_surchem/surchem.h" //sunml add 2022-08-10 #include "module_hamilt_general/module_vdw/vdw.h" #ifdef __DEEPKS -#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03 +#include "module_hamilt_lcao/module_deepks/LCAO_deepks.h" //caoyu add for deepks 2021-06-03 #endif double Force_Stress_LCAO::force_invalid_threshold_ev = 0.00; @@ -33,20 +34,20 @@ void Force_Stress_LCAO::getForceStress( ModuleBase::matrix &scs, const K_Vectors& kv) { - ModuleBase::TITLE("Force_Stress_LCAO","getForceStress"); - ModuleBase::timer::tick("Force_Stress_LCAO","getForceStress"); - - if(!isforce&&!isstress) - { - ModuleBase::timer::tick("Force_Stress_LCAO","getForceStress"); - return; - } + ModuleBase::TITLE("Force_Stress_LCAO", "getForceStress"); + ModuleBase::timer::tick("Force_Stress_LCAO", "getForceStress"); + + if (!isforce && !isstress) + { + ModuleBase::timer::tick("Force_Stress_LCAO", "getForceStress"); + return; + } - const int nat = GlobalC::ucell.nat; + const int nat = GlobalC::ucell.nat; - //total force : ModuleBase::matrix fcs; + // total force : ModuleBase::matrix fcs; - // part of total force + // part of total force ModuleBase::matrix foverlap; ModuleBase::matrix ftvnl_dphi; ModuleBase::matrix fvnl_dbeta; @@ -54,21 +55,21 @@ void Force_Stress_LCAO::getForceStress( ModuleBase::matrix fvl_dvl; ModuleBase::matrix fewalds; ModuleBase::matrix fcc; - ModuleBase::matrix fscc; - - fvl_dphi.create (nat, 3);//must do it now, update it later, noted by zhengdy - - if(isforce) - { - fcs.create (nat, 3); - foverlap.create (nat, 3); - ftvnl_dphi.create (nat, 3); - fvnl_dbeta.create (nat, 3); - fvl_dvl.create (nat, 3); - fewalds.create (nat, 3); - fcc.create (nat, 3); - fscc.create (nat, 3); - //calculate basic terms in Force, same method with PW base + ModuleBase::matrix fscc; + + fvl_dphi.create(nat, 3); // must do it now, update it later, noted by zhengdy + + if (isforce) + { + fcs.create(nat, 3); + foverlap.create(nat, 3); + ftvnl_dphi.create(nat, 3); + fvnl_dbeta.create(nat, 3); + fvl_dvl.create(nat, 3); + fewalds.create(nat, 3); + fcc.create(nat, 3); + fscc.create(nat, 3); + // calculate basic terms in Force, same method with PW base this->calForcePwPart(fvl_dvl, fewalds, fcc, @@ -170,6 +171,15 @@ void Force_Stress_LCAO::getForceStress( fefield.create(nat, 3); elecstate::Efield::compute_force(GlobalC::ucell, fefield); } + + // implement force from E-field of tddft + ModuleBase::matrix fefield_tddft; + if (GlobalV::ESOLVER_TYPE == "TDDFT" && isforce) + { + fefield_tddft.create(nat, 3); + elecstate::Efield::compute_force(GlobalC::ucell, fefield_tddft); + } + // implement force from gate field ModuleBase::matrix fgate; if (GlobalV::GATE_FLAG && isforce) @@ -278,6 +288,11 @@ void Force_Stress_LCAO::getForceStress( { fcs(iat, i) += fefield(iat, i); } + // E-field force of tddft + if (GlobalV::ESOLVER_TYPE == "TDDFT") + { + fcs(iat, i) += fefield_tddft(iat, i); + } // Gate field force if (GlobalV::GATE_FLAG) { @@ -419,6 +434,11 @@ void Force_Stress_LCAO::getForceStress( f_pw.print("EFIELD FORCE", fefield, 0); // this->print_force("EFIELD FORCE",fefield,1,ry); } + if (GlobalV::ESOLVER_TYPE == "TDDFT") + { + f_pw.print("EFIELD_TDDFT FORCE", fefield_tddft, 0); + // this->print_force("EFIELD_TDDFT FORCE",fefield_tddft,1,ry); + } if (GlobalV::GATE_FLAG) { f_pw.print("GATEFIELD FORCE", fgate, 0); @@ -694,7 +714,7 @@ void Force_Stress_LCAO::calForceStressIntegralPart(const bool isGammaOnly, svl_dphi, svnl_dalpha, #else - svl_dphi, + svl_dphi, #endif uhm); } @@ -763,59 +783,83 @@ void Force_Stress_LCAO::calStressPwPart(ModuleBase::matrix& sigmadvl, { sigmaxc(i, i) = -etxc / GlobalC::ucell.omega; } - //Exchange-correlation for PBE + // Exchange-correlation for PBE sc_pw.stress_gga(sigmaxc, GlobalC::rhopw, chr); return; } #include "module_base/mathzone.h" -//do symmetry for total force +// do symmetry for total force void Force_Stress_LCAO::forceSymmetry(ModuleBase::matrix& fcs) { - double *pos; - double d1,d2,d3; - pos = new double[GlobalC::ucell.nat*3]; - ModuleBase::GlobalFunc::ZEROS(pos, GlobalC::ucell.nat*3); - int iat = 0; - for(int it = 0;it < GlobalC::ucell.ntype;it++) - { - for(int ia =0;ia< GlobalC::ucell.atoms[it].na;ia++) - { - pos[3*iat ] = GlobalC::ucell.atoms[it].taud[ia].x ; - pos[3*iat+1] = GlobalC::ucell.atoms[it].taud[ia].y ; - pos[3*iat+2] = GlobalC::ucell.atoms[it].taud[ia].z; - for(int k=0; k<3; ++k) - { - GlobalC::symm.check_translation( pos[iat*3+k], -floor(pos[iat*3+k])); - GlobalC::symm.check_boundary( pos[iat*3+k] ); - } - iat++; - } - } - - for(int iat=0; iat ELEC_evolve::td_vext_dire_case; -bool ELEC_evolve::out_dipole; -bool ELEC_evolve::out_efield; -double ELEC_evolve::td_print_eij; // the threshold to output Eij elements -int ELEC_evolve::td_edm; // 0: new edm method 1: old edm method - -// this routine only serves for TDDFT using LCAO basis set -void ELEC_evolve::evolve_psi(const int& istep, - hamilt::Hamilt* phm, - Local_Orbital_wfc& lowf, - psi::Psi>* psi, - psi::Psi>* psi_laststep, - std::complex** Hk_laststep, - ModuleBase::matrix& ekb, - int htype, - int propagator, - const int& nks) -{ - ModuleBase::TITLE("ELEC_evolve", "eveolve_psi"); - ModuleBase::timer::tick("ELEC_evolve", "evolve_psi"); - - // pool parallization in future -- mohan note 2021-02-09 - for (int ik = 0; ik < nks; ik++) - { - phm->updateHk(ik); - - ModuleBase::timer::tick("Efficience", "evolve_k"); - Evolve_LCAO_Matrix ELM(lowf.ParaV); - psi->fix_k(ik); - psi_laststep->fix_k(ik); - 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 - - ModuleBase::timer::tick("ELEC_evolve", "evolve_psi"); - return; -} diff --git a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp b/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp deleted file mode 100644 index cd7649f536..0000000000 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.cpp +++ /dev/null @@ -1,1466 +0,0 @@ -#include "LCAO_evolve.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 - -Evolve_LCAO_Matrix::~Evolve_LCAO_Matrix() -{ -} - -inline int globalIndex(int localindex, int nblk, int nprocs, int myproc) -{ - int iblock, gIndex; - iblock = localindex / nblk; - gIndex = (iblock * nprocs + myproc) * nblk + localindex % nblk; - return gIndex; -} - -void Evolve_LCAO_Matrix::evolve_complex_matrix(const int& ik, - hamilt::Hamilt* p_hamilt, - psi::Psi>* psi_k, - psi::Psi>* psi_k_laststep, - 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); - GlobalV::ofs_running << " Start Time : " << ctime(&time_start); - - if (GlobalV::ESOLVER_TYPE == "tddft") - { -#ifdef __MPI - 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 - } - else - { - ModuleBase::WARNING_QUIT("Evolve_LCAO_Matrix::evolve_complex_matrix", - "only esolver_type == tddft cando evolve"); - } - - time_t time_end = time(NULL); - ModuleBase::GlobalFunc::OUT_TIME("evolve(std::complex)", time_start, time_end); - - return; -} - -void Evolve_LCAO_Matrix::using_LAPACK_complex(const int& ik, - hamilt::Hamilt* p_hamilt, - std::complex* psi_k, - std::complex* psi_k_laststep, - double* ekb) const -{ - ModuleBase::TITLE("Evolve_LCAO_Matrix", "using_LAPACK_complex"); - - // Calculate the U operator - - ModuleBase::ComplexMatrix Htmp(GlobalV::NLOCAL, GlobalV::NLOCAL); - ModuleBase::ComplexMatrix Stmp(GlobalV::NLOCAL, GlobalV::NLOCAL); - hamilt::MatrixBlock> h_mat, s_mat; - p_hamilt->matrix(h_mat, s_mat); - for (int i = 0; i < GlobalV::NLOCAL; i++) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - Htmp(i, j) = h_mat.p[i * GlobalV::NLOCAL + j]; - // Htmp(i,j) = (this->LM->Hloc2[i*GlobalV::NLOCAL+j] +this->LM->Hloc2_laststep[i*GlobalV::NLOCAL+j])/2.0; - Stmp(i, j) = s_mat.p[i * GlobalV::NLOCAL + j]; - } - } - - ModuleBase::ComplexMatrix wfc_tmp(GlobalV::NBANDS, GlobalV::NLOCAL, true); - ModuleBase::ComplexMatrix wfc_laststep_tmp(GlobalV::NBANDS, GlobalV::NLOCAL, true); - // wfc_laststep_tmp.c = psi_k_laststep; - - for (int i = 0; i < GlobalV::NBANDS; i++) - { - for (int j = 0; j < GlobalV::NLOCAL; j++) - { - wfc_laststep_tmp.c[i * GlobalV::NLOCAL + j] = psi_k_laststep[i * GlobalV::NLOCAL + j]; - } - } - - /* - GlobalV::ofs_running << " Htmp: " < * work = new std::complex[lwork]; - ModuleBase::GlobalFunc::ZEROS(work, lwork); - int IPIV[GlobalV::NLOCAL]; - - LapackConnector::zgetrf( GlobalV::NLOCAL, GlobalV::NLOCAL, Stmp, GlobalV::NLOCAL, IPIV, &INFO); - LapackConnector::zgetri( GlobalV::NLOCAL, Stmp, GlobalV::NLOCAL, IPIV, work, lwork, &INFO); - */ - /* - std::cout << " S^-1: " <( -S_plus_H(i,j).imag(), S_plus_H(i,j).real() ); - } - } - - ModuleBase::ComplexMatrix Idmat(GlobalV::NLOCAL,GlobalV::NLOCAL); - for(int i=0; i(1.0, 0.0); - else Idmat(i,j) = std::complex(0.0, 0.0); - } - } - */ - ModuleBase::ComplexMatrix Numerator(GlobalV::NLOCAL, GlobalV::NLOCAL); - // Numerator = Idmat - 0.5*INPUT.mdp.md_dt*41.34*Denominator; - // Denominator = Idmat + 0.5*INPUT.mdp.md_dt*41.34*Denominator; - - complex para = {0.0, 1.0}; - para = para * 0.25 * INPUT.mdp.md_dt; - Numerator = Stmp - para * Htmp; - Denominator = Stmp + para * Htmp; - - int info; - int lwork = 3 * GlobalV::NLOCAL - 1; // tmp - std::complex* work = new std::complex[lwork]; - ModuleBase::GlobalFunc::ZEROS(work, lwork); - int* ipiv = new int[GlobalV::NLOCAL]; - - LapackConnector::zgetrf(GlobalV::NLOCAL, GlobalV::NLOCAL, Denominator, GlobalV::NLOCAL, ipiv, &info); - LapackConnector::zgetri(GlobalV::NLOCAL, Denominator, GlobalV::NLOCAL, ipiv, work, lwork, &info); - - ModuleBase::ComplexMatrix U_operator(GlobalV::NLOCAL, GlobalV::NLOCAL); - /* - GlobalV::ofs_running << " Numerator: " <* 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"); - - int print_matrix = 0; - hamilt::MatrixBlock> h_mat, s_mat; - p_hamilt->matrix(h_mat, s_mat); - - complex* Stmp = new complex[this->ParaV->nloc]; - ModuleBase::GlobalFunc::ZEROS(Stmp, this->ParaV->nloc); - BlasConnector::copy(this->ParaV->nloc, s_mat.p, 1, Stmp, 1); - - 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)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - - /// @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 - 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); - - 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, print_matrix - /// @output psi_k - U_to_wfc(nband, nlocal, U_operator, psi_k_laststep, psi_k, print_matrix); - - // (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)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> - - /// @brief compute ekb - /// @input Htmp, psi_k - /// @output ekb - compute_ekb(nband, nlocal, Hold, psi_k, ekb); - - delete[] Stmp; - delete[] Htmp; - delete[] U_operator; - 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_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]; - ModuleBase::GlobalFunc::ZEROS(Numerator, this->ParaV->nloc); - BlasConnector::copy(this->ParaV->nloc, Htmp, 1, Numerator, 1); - - complex* Denominator = new complex[this->ParaV->nloc]; - ModuleBase::GlobalFunc::ZEROS(Denominator, this->ParaV->nloc); - BlasConnector::copy(this->ParaV->nloc, Htmp, 1, Denominator, 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 << 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 - complex alpha = {1.0, 0.0}; - 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); - - 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 << 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); - 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); - lwork = work[0].real(); - work.resize(lwork, 0); - liwotk = iwork[0]; - iwork.resize(liwotk, 0); - // (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); - - if (print_matrix) - { - GlobalV::ofs_running << " fenmu^-1:" << 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 << endl; - } - GlobalV::ofs_running << endl; - GlobalV::ofs_running << " fenzi:" << endl; - for (int i = 0; i < this->ParaV->nrow; i++) - { - 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 << 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, - 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 - - 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::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_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; - MPI_Comm_rank(this->ParaV->comm_2D, &myid); - int naroc[2]; // maximum number of row or column - - 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 >= nband) - continue; - for (int i = 0; i < naroc[0]; ++i) - { - int igrow = globalIndex(i, this->ParaV->nb, this->ParaV->dim0, iprow); - if (igrow >= nband) - continue; - if (igcol == igrow) - { - Cij[j * naroc[0] + i] = {1.0 / sqrt(Cij[j * naroc[0] + i].real()), 0.0}; - } - else - { - Cij[j * naroc[0] + i] = {0.0, 0.0}; - } - } - } - } - } // loop ipcol - } // loop iprow - - 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_Eij, - 0.0, - psi_k, - 1, - 1, - this->ParaV->desc_wfc); - - if (print_matrix) - { - GlobalV::ofs_running << " Cij:" << endl; - for (int i = 0; i < this->ParaV->ncol; i++) - { - for (int j = 0; j < this->ParaV->nrow; j++) - { - GlobalV::ofs_running << Cij[i * this->ParaV->ncol + j].real() << "+" - << Cij[i * this->ParaV->ncol + j].imag() << "i "; - } - GlobalV::ofs_running << endl; - } - GlobalV::ofs_running << endl; - 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 before normalization:" << endl; - for (int i = 0; i < this->ParaV->ncol; i++) - { - for (int j = 0; j < this->ParaV->ncol; j++) - { - double aa, bb; - aa = tmp1[i * this->ParaV->ncol + j].real(); - bb = tmp1[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 << 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 -{ - - complex* tmp1 = new complex[this->ParaV->nloc_wfc]; - ModuleBase::GlobalFunc::ZEROS(tmp1, this->ParaV->nloc_wfc); - - 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_Eij); - - if (ELEC_evolve::td_print_eij > 0.0) - { - GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" - << endl; - GlobalV::ofs_running << " Eij:" << endl; - for (int i = 0; i < this->ParaV->nrow_bands; i++) - { - for (int j = 0; j < this->ParaV->ncol_bands; j++) - { - double aa, bb; - aa = Eij[i * this->ParaV->ncol + j].real(); - bb = Eij[i * this->ParaV->ncol + j].imag(); - if (abs(aa) < ELEC_evolve::td_print_eij) - aa = 0.0; - if (abs(bb) < ELEC_evolve::td_print_eij) - bb = 0.0; - if (aa > 0.0 || bb > 0.0) - { - GlobalV::ofs_running << i << " " << j << " " << aa << "+" << bb << "i " << endl; - } - } - } - GlobalV::ofs_running << endl; - GlobalV::ofs_running - << "------------------------------------------------------------------------------------------------" - << 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) - { - 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 >= nband) - continue; - for (int i = 0; i < naroc[0]; ++i) - { - int igrow = globalIndex(i, this->ParaV->nb, this->ParaV->dim0, iprow); - if (igrow >= nband) - continue; - if (igcol == igrow) - { - Eii[igcol] = Eij[j * naroc[0] + i].real(); - } - } - } - } - } // 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 deleted file mode 100644 index bcf470a771..0000000000 --- a/source/module_hamilt_lcao/module_tddft/LCAO_evolve.h +++ /dev/null @@ -1,108 +0,0 @@ -#ifndef EVOLVE_LCAO_MATRIX_H -#define EVOLVE_LCAO_MATRIX_H - -#include "module_base/complexmatrix.h" -#include "module_base/global_function.h" -#include "module_base/global_variable.h" -#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" -#include "module_psi/psi.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_matrix.h" -#include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_wfc.h" - -class Evolve_LCAO_Matrix -{ - public: - Evolve_LCAO_Matrix(const Parallel_Orbitals* pv) - { - this->ParaV = pv; - } - ~Evolve_LCAO_Matrix(); - - void evolve_complex_matrix(const int& ik, - hamilt::Hamilt* p_hamilt, - psi::Psi>* psi_k, - psi::Psi>* psi_k_laststep, - std::complex* H_laststep, - double* ekb, - int htype, - int propagator) const; - - private: - // LCAO_Matrix* LM; - const Parallel_Orbitals* ParaV; - - void using_LAPACK_complex(const int& ik, - hamilt::Hamilt* p_hamilt, - std::complex* psi_k, - std::complex* psi_k_laststep, - double* ekb) const; -#ifdef __MPI - void 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; - - 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_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; - - void 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; - - void norm_wfc( - const int nband, - const int nlocal, - const std::complex* Stmp, - std::complex* psi_k, - const int print_matrix) const; - - - void compute_ekb( - const int nband, - const int nlocal, - const std::complex* Htmp, - const std::complex* psi_k, - double* ekb) const; - -#endif -}; -#endif diff --git a/source/module_hamilt_lcao/module_tddft/bandenergy.cpp b/source/module_hamilt_lcao/module_tddft/bandenergy.cpp new file mode 100644 index 0000000000..d09f1129c8 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/bandenergy.cpp @@ -0,0 +1,151 @@ +#include "bandenergy.h" + +#include +#include + +#include "evolve_elec.h" +#include "module_base/lapack_connector.h" +#include "module_base/scalapack_connector.h" + +namespace module_tddft +{ +#ifdef __MPI + +inline int globalIndex(int localindex, int nblk, int nprocs, int myproc) +{ + int iblock, gIndex; + iblock = localindex / nblk; + gIndex = (iblock * nprocs + myproc) * nblk + localindex % nblk; + return gIndex; +} + +void compute_ekb(const Parallel_Orbitals* pv, + const int nband, + const int nlocal, + const std::complex* Htmp, + const std::complex* psi_k, + double* ekb) +{ + + std::complex* tmp1 = new std::complex[pv->nloc_wfc]; + ModuleBase::GlobalFunc::ZEROS(tmp1, pv->nloc_wfc); + + std::complex* Eij = new std::complex[pv->nloc]; + ModuleBase::GlobalFunc::ZEROS(Eij, pv->nloc); + + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + Htmp, + 1, + 1, + pv->desc, + psi_k, + 1, + 1, + pv->desc_wfc, + 0.0, + tmp1, + 1, + 1, + pv->desc_wfc); + + ScalapackConnector::gemm('C', + 'N', + nband, + nband, + nlocal, + 1.0, + psi_k, + 1, + 1, + pv->desc_wfc, + tmp1, + 1, + 1, + pv->desc_wfc, + 0.0, + Eij, + 1, + 1, + pv->desc_Eij); + + if (Evolve_elec::td_print_eij > 0.0) + { + GlobalV::ofs_running + << "------------------------------------------------------------------------------------------------" + << std::endl; + GlobalV::ofs_running << " Eij:" << std::endl; + for (int i = 0; i < pv->nrow_bands; i++) + { + for (int j = 0; j < pv->ncol_bands; j++) + { + double aa, bb; + aa = Eij[i * pv->ncol + j].real(); + bb = Eij[i * pv->ncol + j].imag(); + if (abs(aa) < Evolve_elec::td_print_eij) + aa = 0.0; + if (abs(bb) < Evolve_elec::td_print_eij) + bb = 0.0; + if (aa > 0.0 || bb > 0.0) + { + GlobalV::ofs_running << i << " " << j << " " << aa << "+" << bb << "i " << std::endl; + } + } + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running + << "------------------------------------------------------------------------------------------------" + << std::endl; + } + + int info; + int myid; + int naroc[2]; + MPI_Comm_rank(pv->comm_2D, &myid); + + double* Eii = new double[nband]; + ModuleBase::GlobalFunc::ZEROS(Eii, nband); + for (int iprow = 0; iprow < pv->dim0; ++iprow) + { + for (int ipcol = 0; ipcol < pv->dim1; ++ipcol) + { + const int coord[2] = {iprow, ipcol}; + int src_rank; + info = MPI_Cart_rank(pv->comm_2D, coord, &src_rank); + if (myid == src_rank) + { + naroc[0] = pv->nrow; + naroc[1] = pv->ncol; + for (int j = 0; j < naroc[1]; ++j) + { + int igcol = globalIndex(j, pv->nb, pv->dim1, ipcol); + if (igcol >= nband) + continue; + for (int i = 0; i < naroc[0]; ++i) + { + int igrow = globalIndex(i, pv->nb, pv->dim0, iprow); + if (igrow >= nband) + continue; + if (igcol == igrow) + { + Eii[igcol] = Eij[j * naroc[0] + i].real(); + } + } + } + } + } // loop ipcol + } // loop iprow + info = MPI_Allreduce(Eii, ekb, nband, MPI_DOUBLE, MPI_SUM, pv->comm_2D); + + delete[] tmp1; + delete[] Eij; + delete[] Eii; +} + +#endif + +} // namespace module_tddft \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/bandenergy.h b/source/module_hamilt_lcao/module_tddft/bandenergy.h new file mode 100644 index 0000000000..e7f361f658 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/bandenergy.h @@ -0,0 +1,32 @@ +/** + * @file bandenegy.h + * @brief compute band energy ekb + * This file originally belonged to file LCAO_evolve.cpp + */ +#ifndef BANDENERGY_H +#define BANDENERGY_H + +#include "module_basis/module_ao/parallel_orbitals.h" + +namespace module_tddft +{ +#ifdef __MPI +/** + * @brief compute band energy ekb + * + * @param[in] pv information of parallel + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] Htmp Hamiltonian + * @param[in] psi_k psi of this step + * @param[out] ekb band energy + */ +void compute_ekb(const Parallel_Orbitals* pv, + const int nband, + const int nlocal, + const std::complex* Htmp, + const std::complex* psi_k, + double* ekb); +#endif +} // namespace module_tddft +#endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp b/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp new file mode 100644 index 0000000000..d4139ba6a6 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/evolve_elec.cpp @@ -0,0 +1,85 @@ +#include "evolve_elec.h" + +#include "evolve_psi.h" +#include "module_base/parallel_reduce.h" +#include "module_base/timer.h" +#include "module_elecstate/module_charge/symmetry_rho.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" + +namespace module_tddft +{ +Evolve_elec::Evolve_elec(){}; +Evolve_elec::~Evolve_elec(){}; + +double Evolve_elec::td_force_dt; +bool Evolve_elec::td_vext; +std::vector Evolve_elec::td_vext_dire_case; +bool Evolve_elec::out_dipole; +bool Evolve_elec::out_efield; +double Evolve_elec::td_print_eij; // the threshold to output Eij elements +int Evolve_elec::td_edm; // 0: new edm method 1: old edm method + +// this routine only serves for TDDFT using LCAO basis set +void Evolve_elec::solve_psi(const int& istep, + const int nband, + const int nlocal, + hamilt::Hamilt* phm, + Local_Orbital_wfc& lowf, + psi::Psi>* psi, + psi::Psi>* psi_laststep, + std::complex** Hk_laststep, + ModuleBase::matrix& ekb, + int htype, + int propagator, + const int& nks) +{ + ModuleBase::TITLE("Evolve_elec", "eveolve_psi"); + ModuleBase::timer::tick("Evolve_elec", "evolve_psi"); + + for (int ik = 0; ik < nks; ik++) + { + phm->updateHk(ik); + + ModuleBase::timer::tick("Efficience", "evolve_k"); + psi->fix_k(ik); + psi_laststep->fix_k(ik); + if (htype == 0) + { + evolve_psi(nband, + nlocal, + lowf.ParaV, + phm, + psi[0].get_pointer(), + psi_laststep[0].get_pointer(), + nullptr, + &(ekb(ik, 0)), + htype, + propagator); + } + else if (htype == 1) + { + evolve_psi(nband, + nlocal, + lowf.ParaV, + phm, + psi[0].get_pointer(), + psi_laststep[0].get_pointer(), + Hk_laststep[ik], + &(ekb(ik, 0)), + htype, + propagator); + } + else + { + std::cout << "method of htype is wrong" << std::endl; + } + + ModuleBase::timer::tick("Efficience", "evolve_k"); + } // end k + + ModuleBase::timer::tick("Evolve_elec", "evolve_psi"); + return; +} +} // namespace module_tddft \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h b/source/module_hamilt_lcao/module_tddft/evolve_elec.h similarity index 55% rename from source/module_hamilt_lcao/module_tddft/ELEC_evolve.h rename to source/module_hamilt_lcao/module_tddft/evolve_elec.h index a40e4391fd..03bb49e078 100644 --- a/source/module_hamilt_lcao/module_tddft/ELEC_evolve.h +++ b/source/module_hamilt_lcao/module_tddft/evolve_elec.h @@ -1,13 +1,13 @@ -#ifndef ELEC_EVOLVE_H -#define ELEC_EVOLVE_H +#ifndef EVOLVE_ELEC_H +#define EVOLVE_ELEC_H #include "module_base/global_function.h" #include "module_base/global_variable.h" -#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h" #include "module_esolver/esolver_ks_lcao.h" #include "module_esolver/esolver_ks_lcao_tddft.h" -#include "module_psi/psi.h" +#include "module_hamilt_lcao/hamilt_lcaodft/LCAO_hamilt.h" #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" +#include "module_psi/psi.h" //----------------------------------------------------------- // mohan add 2021-02-09 @@ -16,7 +16,9 @@ // k is the index for the points in the first Brillouin zone //----------------------------------------------------------- -class ELEC_evolve +namespace module_tddft +{ +class Evolve_elec { friend class ELEC_scf; @@ -24,10 +26,8 @@ class ELEC_evolve friend class ModuleESolver::ESolver_KS_LCAO_TDDFT; public: - ELEC_evolve(); - ~ELEC_evolve(); - - // fuxiang add 2021-05-25 + Evolve_elec(); + ~Evolve_elec(); static double td_force_dt; static bool td_vext; @@ -36,19 +36,21 @@ class ELEC_evolve static bool out_efield; static double td_print_eij; // the threshold to output Eij elements - static int td_edm; // 0: new edm method 1: old edm method + static int td_edm; // 0: new edm method 1: old edm method private: - static void evolve_psi(const int& istep, - hamilt::Hamilt* phm, - Local_Orbital_wfc& lowf, - psi::Psi>* psi, - psi::Psi>* psi_laststep, - std::complex** Hk_laststep, - ModuleBase::matrix& ekb, - int htype, - int propagator, - const int& nks); + static void solve_psi(const int& istep, + const int nband, + const int nlocal, + hamilt::Hamilt* phm, + Local_Orbital_wfc& lowf, + psi::Psi>* psi, + psi::Psi>* psi_laststep, + std::complex** Hk_laststep, + ModuleBase::matrix& ekb, + int htype, + int propagator, + const int& nks); }; - +} // namespace module_tddft #endif diff --git a/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp b/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp new file mode 100644 index 0000000000..ebd578efdb --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/evolve_psi.cpp @@ -0,0 +1,105 @@ +#include "evolve_psi.h" + +#include + +#include "bandenergy.h" +#include "middle_hamilt.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 "norm_psi.h" +#include "propagator.h" +#include "upsi.h" + +namespace module_tddft +{ +void evolve_psi(const int nband, + const int nlocal, + const Parallel_Orbitals* pv, + hamilt::Hamilt* p_hamilt, + std::complex* psi_k, + std::complex* psi_k_laststep, + std::complex* H_laststep, + double* ekb, + int htype, + int propagator) +{ + ModuleBase::TITLE("Evolve_psi", "evolve_psi"); + time_t time_start = time(NULL); + GlobalV::ofs_running << " Start Time : " << ctime(&time_start); + +#ifdef __MPI + + int print_matrix = 0; + hamilt::MatrixBlock> h_mat, s_mat; + p_hamilt->matrix(h_mat, s_mat); + + std::complex* Stmp = new std::complex[pv->nloc]; + ModuleBase::GlobalFunc::ZEROS(Stmp, pv->nloc); + BlasConnector::copy(pv->nloc, s_mat.p, 1, Stmp, 1); + + std::complex* Htmp = new std::complex[pv->nloc]; + ModuleBase::GlobalFunc::ZEROS(Htmp, pv->nloc); + BlasConnector::copy(pv->nloc, h_mat.p, 1, Htmp, 1); + + std::complex* Hold = new std::complex[pv->nloc]; + ModuleBase::GlobalFunc::ZEROS(Hold, pv->nloc); + BlasConnector::copy(pv->nloc, h_mat.p, 1, Hold, 1); + + std::complex* U_operator = new std::complex[pv->nloc]; + ModuleBase::GlobalFunc::ZEROS(U_operator, pv->nloc); + + // (1)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + /// @brief compute H(t+dt/2) + /// @input H_laststep, Htmp, print_matrix + /// @output Htmp + if (htype == 1 && propagator != 2) + { + half_Hmatrix(pv, nband, nlocal, Htmp, H_laststep, print_matrix); + } + + // (2)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + /// @brief compute U_operator + /// @input Stmp, Htmp, print_matrix + /// @output U_operator + Propagator prop(propagator, pv); + prop.compute_propagator(nband, nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); + + // (3)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + /// @brief apply U_operator to the wave function of the previous step for new wave function + /// @input U_operator, psi_k_laststep, print_matrix + /// @output psi_k + upsi(pv, nband, nlocal, U_operator, psi_k_laststep, psi_k, print_matrix); + + // (4)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + /// @brief normalize psi_k + /// @input Stmp, psi_not_norm, psi_k, print_matrix + /// @output psi_k + norm_psi(pv, nband, nlocal, Stmp, psi_k, print_matrix); + + // (5)->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + + /// @brief compute ekb + /// @input Htmp, psi_k + /// @output ekb + compute_ekb(pv, nband, nlocal, Hold, psi_k, ekb); + + delete[] Stmp; + delete[] Htmp; + delete[] Hold; + delete[] U_operator; + +#endif + + time_t time_end = time(NULL); + ModuleBase::GlobalFunc::OUT_TIME("evolve(std::complex)", time_start, time_end); + + return; +} +} // namespace module_tddft \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/evolve_psi.h b/source/module_hamilt_lcao/module_tddft/evolve_psi.h new file mode 100644 index 0000000000..e905b837a5 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/evolve_psi.h @@ -0,0 +1,26 @@ +/** + * @file evolve_psi.h + * @brief evolve the wave function + * This file originally belonged to file LCAO_evolve.cpp + */ +#ifndef ELEC_PSI_H +#define ELEC_PSI_H + +#include "module_basis/module_ao/parallel_orbitals.h" +#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h" + +namespace module_tddft +{ +void evolve_psi(const int nband, + const int nlocal, + const Parallel_Orbitals* pv, + hamilt::Hamilt* p_hamilt, + std::complex* psi_k, + std::complex* psi_k_laststep, + std::complex* H_laststep, + double* ekb, + int htype, + int propagator); +} + +#endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/middle_hamilt.cpp b/source/module_hamilt_lcao/module_tddft/middle_hamilt.cpp new file mode 100644 index 0000000000..72ed8b68e0 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/middle_hamilt.cpp @@ -0,0 +1,68 @@ +#include "middle_hamilt.h" + +#include +#include + +#include "module_base/lapack_connector.h" +#include "module_base/scalapack_connector.h" + +namespace module_tddft +{ +#ifdef __MPI + +void half_Hmatrix(const Parallel_Orbitals* pv, + const int nband, + const int nlocal, + std::complex* Htmp, + const std::complex* H_laststep, + const int print_matrix) +{ + if (print_matrix) + { + GlobalV::ofs_running << std::setprecision(10); + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " H(t+dt) :" << std::endl; + for (int i = 0; i < pv->nrow; i++) + { + for (int j = 0; j < pv->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * pv->ncol + j].real() << "+" << Htmp[i * pv->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " H(t):" << std::endl; + for (int i = 0; i < pv->nrow; i++) + { + for (int j = 0; j < pv->ncol; j++) + { + GlobalV::ofs_running << H_laststep[i * pv->ncol + j].real() << "+" + << H_laststep[i * pv->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << std::endl; + } + GlobalV::ofs_running << std::endl; + } + + std::complex alpha = {0.5, 0.0}; + std::complex beta = {0.5, 0.0}; + ScalapackConnector::geadd('N', nlocal, nlocal, alpha, H_laststep, 1, 1, pv->desc, beta, Htmp, 1, 1, pv->desc); + + if (print_matrix) + { + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " H (t+dt/2) :" << std::endl; + for (int i = 0; i < pv->nrow; i++) + { + for (int j = 0; j < pv->ncol; j++) + { + GlobalV::ofs_running << Htmp[i * pv->ncol + j].real() << "+" << Htmp[i * pv->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << std::endl; + } + GlobalV::ofs_running << std::endl; + } +} +#endif +} // namespace module_tddft \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/middle_hamilt.h b/source/module_hamilt_lcao/module_tddft/middle_hamilt.h new file mode 100644 index 0000000000..6f6da7e8f8 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/middle_hamilt.h @@ -0,0 +1,34 @@ +/** + * @file middle_hamilt.h + * @brief compute H(t+dt/2) + * This file originally belonged to file LCAO_evolve.cpp + */ +#ifndef MIDDLE_HAMILT_H +#define MIDDLE_HAMILT_H + +#include "module_basis/module_ao/parallel_orbitals.h" + +namespace module_tddft +{ +#ifdef __MPI +/** + * @brief compute H(t+dt/2) + * + * @param[in] pv information of parallel + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] Htmp H(t+dt) + * @param[in] H_laststep H(t) + * @param[in] print_matirx print internal matrix or not + * @param[out] Htmp H(t+dt/2) + */ +void half_Hmatrix(const Parallel_Orbitals* pv, + const int nband, + const int nlocal, + std::complex* Htmp, + const std::complex* H_laststep, + const int print_matrix); +#endif +} // namespace module_tddft + +#endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/norm_psi.cpp b/source/module_hamilt_lcao/module_tddft/norm_psi.cpp new file mode 100644 index 0000000000..19cd0a5f84 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/norm_psi.cpp @@ -0,0 +1,211 @@ +#include "norm_psi.h" + +#include +#include + +#include "module_base/lapack_connector.h" +#include "module_base/scalapack_connector.h" + +namespace module_tddft +{ +#ifdef __MPI + +inline int globalIndex(int localindex, int nblk, int nprocs, int myproc) +{ + int iblock, gIndex; + iblock = localindex / nblk; + gIndex = (iblock * nprocs + myproc) * nblk + localindex % nblk; + return gIndex; +} + +void norm_psi(const Parallel_Orbitals* pv, + const int nband, + const int nlocal, + const std::complex* Stmp, + std::complex* psi_k, + const int print_matrix) +{ + std::complex* tmp1 = new std::complex[pv->nloc_wfc]; + ModuleBase::GlobalFunc::ZEROS(tmp1, pv->nloc_wfc); + + std::complex* Cij = new std::complex[pv->nloc]; + ModuleBase::GlobalFunc::ZEROS(Cij, pv->nloc); + + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + Stmp, + 1, + 1, + pv->desc, + psi_k, + 1, + 1, + pv->desc_wfc, + 0.0, + tmp1, + 1, + 1, + pv->desc_wfc); + + ScalapackConnector::gemm('C', + 'N', + nband, + nband, + nlocal, + 1.0, + psi_k, + 1, + 1, + pv->desc_wfc, + tmp1, + 1, + 1, + pv->desc_wfc, + 0.0, + Cij, + 1, + 1, + pv->desc_Eij); + + if (print_matrix) + { + GlobalV::ofs_running << "original Cij :" << std::endl; + for (int i = 0; i < pv->ncol; i++) + { + for (int j = 0; j < pv->nrow; j++) + { + double aa, bb; + aa = Cij[i * pv->ncol + j].real(); + bb = Cij[i * pv->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 << std::endl; + } + GlobalV::ofs_running << std::endl; + } + + int info; + int myid; + MPI_Comm_rank(pv->comm_2D, &myid); + int naroc[2]; // maximum number of row or column + + for (int iprow = 0; iprow < pv->dim0; ++iprow) + { + for (int ipcol = 0; ipcol < pv->dim1; ++ipcol) + { + const int coord[2] = {iprow, ipcol}; + int src_rank; + info = MPI_Cart_rank(pv->comm_2D, coord, &src_rank); + if (myid == src_rank) + { + naroc[0] = pv->nrow; + naroc[1] = pv->ncol; + for (int j = 0; j < naroc[1]; ++j) + { + int igcol = globalIndex(j, pv->nb, pv->dim1, ipcol); + if (igcol >= nband) + continue; + for (int i = 0; i < naroc[0]; ++i) + { + int igrow = globalIndex(i, pv->nb, pv->dim0, iprow); + if (igrow >= nband) + continue; + if (igcol == igrow) + { + Cij[j * naroc[0] + i] = {1.0 / sqrt(Cij[j * naroc[0] + i].real()), 0.0}; + } + else + { + Cij[j * naroc[0] + i] = {0.0, 0.0}; + } + } + } + } + } // loop ipcol + } // loop iprow + + BlasConnector::copy(pv->nloc_wfc, psi_k, 1, tmp1, 1); + + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nband, + 1.0, + tmp1, + 1, + 1, + pv->desc_wfc, + Cij, + 1, + 1, + pv->desc_Eij, + 0.0, + psi_k, + 1, + 1, + pv->desc_wfc); + + if (print_matrix) + { + GlobalV::ofs_running << " Cij:" << std::endl; + for (int i = 0; i < pv->ncol; i++) + { + for (int j = 0; j < pv->nrow; j++) + { + GlobalV::ofs_running << Cij[i * pv->ncol + j].real() << "+" << Cij[i * pv->ncol + j].imag() << "i "; + } + GlobalV::ofs_running << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " psi_k:" << std::endl; + for (int i = 0; i < pv->ncol_bands; i++) + { + for (int j = 0; j < pv->ncol; j++) + { + double aa, bb; + aa = psi_k[i * pv->ncol + j].real(); + bb = psi_k[i * pv->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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " psi_k before normalization:" << std::endl; + for (int i = 0; i < pv->ncol; i++) + { + for (int j = 0; j < pv->ncol; j++) + { + double aa, bb; + aa = tmp1[i * pv->ncol + j].real(); + bb = tmp1[i * pv->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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << std::endl; + } + + delete[] tmp1; + delete[] Cij; +} +#endif +} // namespace module_tddft diff --git a/source/module_hamilt_lcao/module_tddft/norm_psi.h b/source/module_hamilt_lcao/module_tddft/norm_psi.h new file mode 100644 index 0000000000..69287a1532 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/norm_psi.h @@ -0,0 +1,35 @@ +/** + * @file norm_psi.h + * @brief normalize the wave function + * This file originally belonged to file LCAO_evolve.cpp + */ +#ifndef NORM_PSI_H +#define NORM_PSI_H + +#include "module_basis/module_ao/parallel_orbitals.h" + +namespace module_tddft +{ +#ifdef __MPI +/** + * @brief normalize the wave function + * + * @param[in] pv information of parallel + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] Stmp overlap matrix + * @param[in] psi_k psi of this step + * @param[in] print_matirx print internal matrix or not + * @param[out] psi_k psi of this step after normalization + */ +void norm_psi(const Parallel_Orbitals* pv, + const int nband, + const int nlocal, + const std::complex* Stmp, + std::complex* psi_k, + const int print_matrix); + +#endif +} // namespace module_tddft + +#endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/propagator.cpp b/source/module_hamilt_lcao/module_tddft/propagator.cpp new file mode 100644 index 0000000000..f2711ed2f8 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/propagator.cpp @@ -0,0 +1,620 @@ +#include "propagator.h" + +#include +#include + +#include "module_base/lapack_connector.h" +#include "module_base/scalapack_connector.h" +#include "module_io/input.h" + +namespace module_tddft +{ +Propagator::~Propagator() +{ +} + +#ifdef __MPI + +inline int globalIndex(int localindex, int nblk, int nprocs, int myproc) +{ + int iblock, gIndex; + iblock = localindex / nblk; + gIndex = (iblock * nprocs + myproc) * nblk + localindex % nblk; + return gIndex; +} + +void Propagator::compute_propagator(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 +{ + int tag; + switch (ptype) + { + case 0: + compute_propagator_cn2(nband, nlocal, Stmp, Htmp, U_operator, print_matrix); + break; + + case 1: + tag = 1; + compute_propagator_taylor(nband, nlocal, Stmp, Htmp, U_operator, print_matrix, tag); + break; + + case 2: + compute_propagator_etrs(nband, nlocal, Stmp, Htmp, H_laststep, U_operator, print_matrix); + + break; + + default: + std::cout << "method of propagator is wrong" << std::endl; + break; + } +} + +void Propagator::compute_propagator_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 + std::complex* Numerator = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Numerator, this->ParaV->nloc); + BlasConnector::copy(this->ParaV->nloc, Htmp, 1, Numerator, 1); + + std::complex* Denominator = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(Denominator, this->ParaV->nloc); + BlasConnector::copy(this->ParaV->nloc, Htmp, 1, Denominator, 1); + + if (print_matrix) + { + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " S matrix :" << std::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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " H matrix :" << std::endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + 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 << std::endl; + } + GlobalV::ofs_running << std::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 + std::complex alpha = {1.0, 0.0}; + std::complex beta1 = {0.0, -0.25 * INPUT.mdp.md_dt}; + std::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); + + if (print_matrix) + { + GlobalV::ofs_running << " beta=" << beta1 << std::endl; + GlobalV::ofs_running << " fenmu:" << std::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 << std::endl; + } + GlobalV::ofs_running << std::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); + 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); + lwork = work[0].real(); + work.resize(lwork, 0); + liwotk = iwork[0]; + iwork.resize(liwotk, 0); + // (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); + + if (print_matrix) + { + GlobalV::ofs_running << " fenmu^-1:" << std::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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " fenzi:" << std::endl; + for (int i = 0; i < this->ParaV->nrow; i++) + { + 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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " U operator:" << std::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 << std::endl; + } + } + + delete[] Numerator; + delete[] Denominator; + delete[] ipiv; +} + +void Propagator::compute_propagator_taylor(const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix, + const int tag) const +{ + ModuleBase::GlobalFunc::ZEROS(U_operator, this->ParaV->nloc); + std::complex* A_matrix = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(A_matrix, this->ParaV->nloc); + std::complex* rank0 = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank0, this->ParaV->nloc); + std::complex* rank2 = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank2, this->ParaV->nloc); + std::complex* rank3 = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank3, this->ParaV->nloc); + std::complex* rank4 = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(rank4, this->ParaV->nloc); + std::complex* tmp1 = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(tmp1, this->ParaV->nloc); + std::complex* tmp2 = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(tmp2, this->ParaV->nloc); + std::complex* Sinv = new std::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 << std::endl; + GlobalV::ofs_running << " S matrix :" << std::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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " H matrix :" << std::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 << std::endl; + } + GlobalV::ofs_running << std::endl; + } + + // set rank0 + int info; + int myid; + MPI_Comm_rank(this->ParaV->comm_2D, &myid); + int naroc[2]; // maximum number of row or column + + 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 + + std::complex beta = {0.0, -0.5 * INPUT.mdp.md_dt / tag}; // for ETRS tag=2 , for taylor tag=1 + + //->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + // 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); + + std::complex p1 = {1.0, 0.0}; + std::complex p2 = {1.0 / 2.0, 0.0}; + std::complex p3 = {1.0 / 6.0, 0.0}; + std::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:" << std::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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " U operator:" << std::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 << std::endl; + } + } + + delete[] rank0; + delete[] rank2; + delete[] rank3; + delete[] rank4; + delete[] tmp1; + delete[] tmp2; + delete[] ipiv; +} + +void Propagator::compute_propagator_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 +{ + std::complex* U1 = new std::complex[this->ParaV->nloc]; + std::complex* U2 = new std::complex[this->ParaV->nloc]; + ModuleBase::GlobalFunc::ZEROS(U1, this->ParaV->nloc); + ModuleBase::GlobalFunc::ZEROS(U2, this->ParaV->nloc); + int tag = 2; + compute_propagator_taylor(nband, nlocal, Stmp, Htmp, U1, print_matrix, tag); + compute_propagator_taylor(nband, nlocal, Stmp, H_laststep, U2, print_matrix, tag); + 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); +} + +#endif +} // namespace module_tddft \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/propagator.h b/source/module_hamilt_lcao/module_tddft/propagator.h new file mode 100644 index 0000000000..7d9e65b2c4 --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/propagator.h @@ -0,0 +1,108 @@ +/** + * @file propagator.h + * @brief compute propagtor to evolve the wave function + * This file originally belonged to file LCAO_evolve.cpp + */ +#ifndef PROPAGATOR_H +#define PROPAGATOR_H + +#include "module_basis/module_ao/parallel_orbitals.h" + +namespace module_tddft +{ +class Propagator +{ + public: + Propagator(const int ptype, const Parallel_Orbitals* pv) + { + this->ptype = ptype; + this->ParaV = pv; + } + ~Propagator(); + +#ifdef __MPI + /** + * @brief compute propagator + * + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] Stmp overlap matrix + * @param[in] Htmp H(t+dt/2) or H(t+dt) + * @param[in] H_laststep H(t) + * @param[in] print_matirx print internal matrix or not + * @param[out] U_operator operator of propagator + */ + void compute_propagator(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; +#endif + + private: + int ptype; // type of propagator + const Parallel_Orbitals* ParaV; + +#ifdef __MPI + + /** + * @brief compute propagator of method Crank-Nicolson + * + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] Stmp overlap matrix + * @param[in] Htmp H(t+dt/2) or H(t+dt) + * @param[in] print_matirx print internal matrix or not + * @param[out] U_operator operator of propagator + */ + void compute_propagator_cn2(const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix) const; + + /** + * @brief compute propagator of method 4th Taylor + * + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] Stmp overlap matrix + * @param[in] Htmp H(t+dt/2) or H(t+dt) + * @param[in] print_matirx print internal matrix or not + * @param[in] tag a parametre different for 4th Taylor and ETRS + * @param[out] U_operator operator of propagator + */ + void compute_propagator_taylor(const int nband, + const int nlocal, + const std::complex* Stmp, + const std::complex* Htmp, + std::complex* U_operator, + const int print_matrix, + const int tag) const; + + /** + * @brief compute propagator of method ETRS + * + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] Stmp overlap matrix + * @param[in] Htmp H(t+dt/2) or H(t+dt) + * @param[in] H_laststep H(t) + * @param[in] print_matirx print internal matrix or not + * @param[out] U_operator operator of propagator + */ + void compute_propagator_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; +#endif +}; +} // namespace module_tddft + +#endif \ No newline at end of file diff --git a/source/module_hamilt_lcao/module_tddft/upsi.cpp b/source/module_hamilt_lcao/module_tddft/upsi.cpp new file mode 100644 index 0000000000..35c0b03b5e --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/upsi.cpp @@ -0,0 +1,82 @@ +#include "upsi.h" + +#include +#include + +#include "module_base/lapack_connector.h" +#include "module_base/scalapack_connector.h" + +namespace module_tddft +{ +#ifdef __MPI +void upsi(const Parallel_Orbitals* pv, + 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) +{ + + ScalapackConnector::gemm('N', + 'N', + nlocal, + nband, + nlocal, + 1.0, + U_operator, + 1, + 1, + pv->desc, + psi_k_laststep, + 1, + 1, + pv->desc_wfc, + 0.0, + psi_k, + 1, + 1, + pv->desc_wfc); + + if (print_matrix) + { + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " psi_k:" << std::endl; + for (int i = 0; i < pv->ncol_bands; i++) + { + for (int j = 0; j < pv->ncol; j++) + { + double aa, bb; + aa = psi_k[i * pv->ncol + j].real(); + bb = psi_k[i * pv->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 << std::endl; + } + GlobalV::ofs_running << std::endl; + GlobalV::ofs_running << " psi_k_laststep:" << std::endl; + for (int i = 0; i < pv->ncol_bands; i++) + { + for (int j = 0; j < pv->ncol; j++) + { + double aa, bb; + aa = psi_k_laststep[i * pv->ncol + j].real(); + bb = psi_k_laststep[i * pv->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 << std::endl; + } + GlobalV::ofs_running << std::endl; + } +} + +#endif +} // namespace module_tddft diff --git a/source/module_hamilt_lcao/module_tddft/upsi.h b/source/module_hamilt_lcao/module_tddft/upsi.h new file mode 100644 index 0000000000..7233eee77c --- /dev/null +++ b/source/module_hamilt_lcao/module_tddft/upsi.h @@ -0,0 +1,37 @@ +/** + * @file upsi.h + * @brief apply U_operator to the wave function of the previous step for new wave function + * This file originally belonged to file LCAO_evolve.cpp + */ + +#ifndef UPSI_H +#define UPSI_H + +#include "module_basis/module_ao/parallel_orbitals.h" + +namespace module_tddft +{ +#ifdef __MPI +/** + * @brief apply U_operator to the wave function of the previous step for new wave function + * + * @param[in] pv information of parallel + * @param[in] nband number of bands + * @param[in] nlocal number of orbitals + * @param[in] U_operator operator of propagator + * @param[in] psi_k_laststep psi of last step + * @param[in] print_matirx print internal matrix or not + * @param[out] psi_k psi of this step + */ +void upsi(const Parallel_Orbitals* pv, + 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); + +#endif +} // namespace module_tddft + +#endif \ No newline at end of file diff --git a/source/module_io/input_conv.cpp b/source/module_io/input_conv.cpp index 0b3a21c5df..2bc4345ead 100644 --- a/source/module_io/input_conv.cpp +++ b/source/module_io/input_conv.cpp @@ -15,22 +15,24 @@ #include "module_ri/exx_abfs-jle.h" #endif #ifdef __LCAO +#include "module_basis/module_ao/ORB_read.h" +#include "module_elecstate/potentials/H_TDDFT_pw.h" #include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h" #include "module_hamilt_lcao/hamilt_lcaodft/global_fp.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/module_dftu/dftu.h" -#include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" -#include "module_basis/module_ao/ORB_read.h" +#include "module_hamilt_lcao/module_tddft/evolve_elec.h" #endif #include "module_base/timer.h" #include "module_elecstate/elecstate_lcao.h" #include "module_elecstate/potentials/efield.h" #include "module_elecstate/potentials/gatefield.h" #include "module_hsolver/hsolver_lcao.h" -#include "module_psi/kernels/device.h" #include "module_md/md_func.h" +#include "module_psi/kernels/device.h" -template void Input_Conv::parse_expression(const std::string &fn, std::vector &vec) +template +void Input_Conv::parse_expression(const std::string &fn, std::vector &vec) { ModuleBase::TITLE("Input_Conv", "parse_expression"); int count = 0; @@ -110,6 +112,75 @@ template void Input_Conv::parse_expression(const std::string &fn, s regfree(®); } +#ifdef __LCAO +std::vector Input_Conv::convert_units(std::string params, double c) +{ + std::vector params_ori; + std::vector params_out; + parse_expression(params, params_ori); + for (auto param: params_ori) + params_out.emplace_back(param * c); + + return params_out; +} + +void Input_Conv::read_td_efield() +{ + elecstate::H_TDDFT_pw::stype = INPUT.td_stype; + + parse_expression(INPUT.td_ttype, elecstate::H_TDDFT_pw::ttype); + + elecstate::H_TDDFT_pw::tstart = INPUT.td_tstart; + elecstate::H_TDDFT_pw::tend = INPUT.td_tend; + + elecstate::H_TDDFT_pw::dt = INPUT.mdp.md_dt / ModuleBase::AU_to_FS; + + // space domain parameters + + // length gauge + elecstate::H_TDDFT_pw::lcut1 = INPUT.td_lcut1; + elecstate::H_TDDFT_pw::lcut2 = INPUT.td_lcut2; + + // time domain parameters + + // Gauss + elecstate::H_TDDFT_pw::gauss_omega + = convert_units(INPUT.td_gauss_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 + elecstate::H_TDDFT_pw::gauss_phase = convert_units(INPUT.td_gauss_phase, 1.0); + elecstate::H_TDDFT_pw::gauss_sigma = convert_units(INPUT.td_gauss_sigma, 1 / ModuleBase::AU_to_FS); + elecstate::H_TDDFT_pw::gauss_t0 = convert_units(INPUT.td_gauss_t0, 1.0); + elecstate::H_TDDFT_pw::gauss_amp + = convert_units(INPUT.td_gauss_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr + + // trapezoid + elecstate::H_TDDFT_pw::trape_omega + = convert_units(INPUT.td_trape_freq, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 + elecstate::H_TDDFT_pw::trape_phase = convert_units(INPUT.td_trape_phase, 1.0); + elecstate::H_TDDFT_pw::trape_t1 = convert_units(INPUT.td_trape_t1, 1.0); + elecstate::H_TDDFT_pw::trape_t2 = convert_units(INPUT.td_trape_t2, 1.0); + elecstate::H_TDDFT_pw::trape_t3 = convert_units(INPUT.td_trape_t3, 1.0); + elecstate::H_TDDFT_pw::trape_amp + = convert_units(INPUT.td_trape_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr + + // Trigonometric + elecstate::H_TDDFT_pw::trigo_omega1 + = convert_units(INPUT.td_trigo_freq1, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 + elecstate::H_TDDFT_pw::trigo_omega2 + = convert_units(INPUT.td_trigo_freq2, 2 * ModuleBase::PI * ModuleBase::AU_to_FS); // time(a.u.)^-1 + elecstate::H_TDDFT_pw::trigo_phase1 = convert_units(INPUT.td_trigo_phase1, 1.0); + elecstate::H_TDDFT_pw::trigo_phase2 = convert_units(INPUT.td_trigo_phase2, 1.0); + elecstate::H_TDDFT_pw::trigo_amp + = convert_units(INPUT.td_trigo_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr + + // Heaviside + elecstate::H_TDDFT_pw::heavi_t0 = convert_units(INPUT.td_heavi_t0, 1.0); + elecstate::H_TDDFT_pw::heavi_amp + = convert_units(INPUT.td_heavi_amp, ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV); // Ry/bohr + + return; +} +#endif + void Input_Conv::Convert(void) { ModuleBase::TITLE("Input_Conv", "Convert"); @@ -376,16 +447,17 @@ void Input_Conv::Convert(void) // Fuxiang He add 2016-10-26 //---------------------------------------------------------- #ifdef __LCAO - ELEC_evolve::td_force_dt = INPUT.td_force_dt; - ELEC_evolve::td_vext = INPUT.td_vext; - if (ELEC_evolve::td_vext) + module_tddft::Evolve_elec::td_force_dt = INPUT.td_force_dt; + module_tddft::Evolve_elec::td_vext = INPUT.td_vext; + if (module_tddft::Evolve_elec::td_vext) { - parse_expression(INPUT.td_vext_dire, ELEC_evolve::td_vext_dire_case); + parse_expression(INPUT.td_vext_dire, module_tddft::Evolve_elec::td_vext_dire_case); } - ELEC_evolve::out_dipole = INPUT.out_dipole; - ELEC_evolve::out_efield = INPUT.out_efield ; - ELEC_evolve::td_print_eij = INPUT.td_print_eij; - ELEC_evolve::td_edm = INPUT.td_edm; + module_tddft::Evolve_elec::out_dipole = INPUT.out_dipole; + module_tddft::Evolve_elec::out_efield = INPUT.out_efield; + module_tddft::Evolve_elec::td_print_eij = INPUT.td_print_eij; + module_tddft::Evolve_elec::td_edm = INPUT.td_edm; + read_td_efield(); #endif // setting for constrained DFT, jiyy add 2020.10.11 diff --git a/source/module_io/input_conv.h b/source/module_io/input_conv.h index f25ddb54f7..f553dea9e8 100644 --- a/source/module_io/input_conv.h +++ b/source/module_io/input_conv.h @@ -5,12 +5,13 @@ #ifndef INPUT_CONVERT_H #define INPUT_CONVERT_H -#include -#include -#include #include #include #include + +#include +#include +#include #include #include @@ -24,6 +25,23 @@ void Convert(void); // fn (string): expressions such as "3*1 0 2*0.5 3*0" // arr (vector): stores parsing results, for example, "3*1 0 2*0.5 1*1.5" can be parsed as [1, 1, 1, 0, 0.5, 0.5, 1.5] template void parse_expression(const std::string &fn, std::vector &arr); + +#ifdef __LCAO +/** + * @brief convert units of different parameters + * + * @param params input parameter + * @param c coefficients of unit conversion + * @return parame*c : parameter after unit vonversion + */ +std::vector convert_units(std::string params, double c); + +/** + * @brief read paramers of electric field for tddft and convert units + */ +void read_td_efield(); +#endif + } // namespace Input_Conv #endif // Input_Convert diff --git a/source/module_io/test/for_testing_input_conv.h b/source/module_io/test/for_testing_input_conv.h index eea79479d3..8048fdcdc0 100644 --- a/source/module_io/test/for_testing_input_conv.h +++ b/source/module_io/test/for_testing_input_conv.h @@ -5,12 +5,13 @@ #include "module_elecstate/elecstate_lcao.h" #include "module_elecstate/module_charge/charge_mixing.h" #include "module_elecstate/occupy.h" +#include "module_elecstate/potentials/H_TDDFT_pw.h" #include "module_elecstate/potentials/efield.h" #include "module_elecstate/potentials/gatefield.h" #include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h" #include "module_hamilt_lcao/hamilt_lcaodft/local_orbital_charge.h" #include "module_hamilt_lcao/module_dftu/dftu.h" -#include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" +#include "module_hamilt_lcao/module_tddft/evolve_elec.h" #include "module_hamilt_pw/hamilt_pwdft/VNL_in_pw.h" #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h" #include "module_hamilt_pw/hamilt_pwdft/wavefunc.h" @@ -23,7 +24,7 @@ #include "module_relax/relax_old/ions_move_cg.h" #include "module_relax/relax_old/lattice_change_basic.h" -bool berryphase::berry_phase_flag=false; +bool berryphase::berry_phase_flag = false; int elecstate::ElecStateLCAO::out_wfc_lcao = 0; bool elecstate::ElecStateLCAO::need_psi_grid = 1; int hsolver::HSolverLCAO::out_mat_hs = 0; @@ -32,13 +33,13 @@ int hsolver::HSolverLCAO::out_mat_t = 0; 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; -bool ELEC_evolve::td_vext; -std::vector ELEC_evolve::td_vext_dire_case; -bool ELEC_evolve::out_dipole; -bool ELEC_evolve::out_efield; -double ELEC_evolve::td_print_eij; -int ELEC_evolve::td_edm; +double module_tddft::Evolve_elec::td_force_dt; +bool module_tddft::Evolve_elec::td_vext; +std::vector module_tddft::Evolve_elec::td_vext_dire_case; +bool module_tddft::Evolve_elec::out_dipole; +bool module_tddft::Evolve_elec::out_efield; +double module_tddft::Evolve_elec::td_print_eij; +int module_tddft::Evolve_elec::td_edm; double elecstate::Gatefield::zgate = 0.5; bool elecstate::Gatefield::relax = false; bool elecstate::Gatefield::block = false; @@ -49,207 +50,326 @@ int elecstate::Efield::efield_dir; double elecstate::Efield::efield_pos_max; double elecstate::Efield::efield_pos_dec; double elecstate::Efield::efield_amp; + +// parameters of electric field for tddft + +int elecstate::H_TDDFT_pw::stype; // 0 : length gauge 1: velocity gauge + +std::vector elecstate::H_TDDFT_pw::ttype; +// 0 Gauss type function. +// 1 trapezoid type function. +// 2 Trigonometric functions, sin^2. +// 3 heaviside function. +// 4 HHG function. + +int elecstate::H_TDDFT_pw::tstart; +int elecstate::H_TDDFT_pw::tend; +double elecstate::H_TDDFT_pw::dt; + +// space domain parameters + +// length gauge +double elecstate::H_TDDFT_pw::lcut1; +double elecstate::H_TDDFT_pw::lcut2; + +// time domain parameters + +// Gauss +int elecstate::H_TDDFT_pw::gauss_count; +std::vector elecstate::H_TDDFT_pw::gauss_omega; // time(a.u.)^-1 +std::vector elecstate::H_TDDFT_pw::gauss_phase; +std::vector elecstate::H_TDDFT_pw::gauss_sigma; // time(a.u.) +std::vector elecstate::H_TDDFT_pw::gauss_t0; +std::vector elecstate::H_TDDFT_pw::gauss_amp; // Ry/bohr + +// trapezoid +int elecstate::H_TDDFT_pw::trape_count; +std::vector elecstate::H_TDDFT_pw::trape_omega; // time(a.u.)^-1 +std::vector elecstate::H_TDDFT_pw::trape_phase; +std::vector elecstate::H_TDDFT_pw::trape_t1; +std::vector elecstate::H_TDDFT_pw::trape_t2; +std::vector elecstate::H_TDDFT_pw::trape_t3; +std::vector elecstate::H_TDDFT_pw::trape_amp; // Ry/bohr + +// Trigonometric +int elecstate::H_TDDFT_pw::trigo_count; +std::vector elecstate::H_TDDFT_pw::trigo_omega1; // time(a.u.)^-1 +std::vector elecstate::H_TDDFT_pw::trigo_omega2; // time(a.u.)^-1 +std::vector elecstate::H_TDDFT_pw::trigo_phase1; +std::vector elecstate::H_TDDFT_pw::trigo_phase2; +std::vector elecstate::H_TDDFT_pw::trigo_amp; // Ry/bohr + +// Heaviside +int elecstate::H_TDDFT_pw::heavi_count; +std::vector elecstate::H_TDDFT_pw::heavi_t0; +std::vector elecstate::H_TDDFT_pw::heavi_amp; // Ry/bohr + double Force_Stress_LCAO::force_invalid_threshold_ev = 0.0; double BFGS_Basic::relax_bfgs_w1 = -1.0; double BFGS_Basic::relax_bfgs_w2 = -1.0; double Ions_Move_Basic::relax_bfgs_rmax = -1.0; double Ions_Move_Basic::relax_bfgs_rmin = -1.0; double Ions_Move_Basic::relax_bfgs_init = -1.0; -int Ions_Move_Basic::out_stru=0; -double Ions_Move_CG::RELAX_CG_THR =-1.0; +int Ions_Move_Basic::out_stru = 0; +double Ions_Move_CG::RELAX_CG_THR = -1.0; std::string Lattice_Change_Basic::fixed_axes = "None"; -int ModuleSymmetry::Symmetry::symm_flag=0; - -Charge_Mixing::Charge_Mixing(){} -Charge_Mixing::~Charge_Mixing(){} -pseudopot_cell_vnl::pseudopot_cell_vnl(){} -pseudopot_cell_vnl::~pseudopot_cell_vnl(){} -pseudopot_cell_vl::pseudopot_cell_vl(){} -pseudopot_cell_vl::~pseudopot_cell_vl(){} -ORB_gaunt_table::ORB_gaunt_table(){} -ORB_gaunt_table::~ORB_gaunt_table(){} -ModuleDFTU::DFTU::DFTU(){} -ModuleDFTU::DFTU::~DFTU(){} -Structure_Factor::Structure_Factor(){} -Structure_Factor::~Structure_Factor(){} -ModuleSymmetry::Symmetry::Symmetry(){} -ModuleSymmetry::Symmetry::~Symmetry(){} -ModuleSymmetry::Symmetry_Basic::Symmetry_Basic(){} +int ModuleSymmetry::Symmetry::symm_flag = 0; + +Charge_Mixing::Charge_Mixing() +{ +} +Charge_Mixing::~Charge_Mixing() +{ +} +pseudopot_cell_vnl::pseudopot_cell_vnl() +{ +} +pseudopot_cell_vnl::~pseudopot_cell_vnl() +{ +} +pseudopot_cell_vl::pseudopot_cell_vl() +{ +} +pseudopot_cell_vl::~pseudopot_cell_vl() +{ +} +ORB_gaunt_table::ORB_gaunt_table() +{ +} +ORB_gaunt_table::~ORB_gaunt_table() +{ +} +ModuleDFTU::DFTU::DFTU() +{ +} +ModuleDFTU::DFTU::~DFTU() +{ +} +Structure_Factor::Structure_Factor() +{ +} +Structure_Factor::~Structure_Factor() +{ +} +ModuleSymmetry::Symmetry::Symmetry() +{ +} +ModuleSymmetry::Symmetry::~Symmetry() +{ +} +ModuleSymmetry::Symmetry_Basic::Symmetry_Basic() +{ +} ModuleSymmetry::Symmetry_Basic::~Symmetry_Basic() { } -WF_atomic::WF_atomic(){} -WF_atomic::~WF_atomic(){} -wavefunc::wavefunc(){} -wavefunc::~wavefunc(){} -UnitCell::UnitCell(){ - if (GlobalV::test_unitcell) - ModuleBase::TITLE("unitcell", "Constructor"); - Coordinate = "Direct"; - latName = "none"; - lat0 = 0.0; - lat0_angstrom = 0.0; - - bool init_vel; - - ntype = 0; - nat = 0; - namax = 0; - nwmax = 0; - - iat2it = nullptr; - iat2ia = nullptr; - iwt2iat = nullptr; - iwt2iw = nullptr; - - itia2iat.create(1, 1); - lc = new int[3]; - - latvec = ModuleBase::Matrix3(); - latvec_supercell = ModuleBase::Matrix3(); - G = ModuleBase::Matrix3(); - GT = ModuleBase::Matrix3(); - GGT = ModuleBase::Matrix3(); - invGGT = ModuleBase::Matrix3(); - - tpiba = 0.0; - tpiba2 = 0.0; - omega = 0.0; - - atom_label = new string[1]; - atom_mass = nullptr; - pseudo_fn = new string[1]; - pseudo_type = new string[1]; - orbital_fn = new string[1]; - - set_atom_flag = false; -} -UnitCell::~UnitCell(){} +WF_atomic::WF_atomic() +{ +} +WF_atomic::~WF_atomic() +{ +} +wavefunc::wavefunc() +{ +} +wavefunc::~wavefunc() +{ +} +UnitCell::UnitCell() +{ + if (GlobalV::test_unitcell) + ModuleBase::TITLE("unitcell", "Constructor"); + Coordinate = "Direct"; + latName = "none"; + lat0 = 0.0; + lat0_angstrom = 0.0; + + bool init_vel; + + ntype = 0; + nat = 0; + namax = 0; + nwmax = 0; + + iat2it = nullptr; + iat2ia = nullptr; + iwt2iat = nullptr; + iwt2iw = nullptr; + + itia2iat.create(1, 1); + lc = new int[3]; + + latvec = ModuleBase::Matrix3(); + latvec_supercell = ModuleBase::Matrix3(); + G = ModuleBase::Matrix3(); + GT = ModuleBase::Matrix3(); + GGT = ModuleBase::Matrix3(); + invGGT = ModuleBase::Matrix3(); + + tpiba = 0.0; + tpiba2 = 0.0; + omega = 0.0; + + atom_label = new string[1]; + atom_mass = nullptr; + pseudo_fn = new string[1]; + pseudo_type = new string[1]; + orbital_fn = new string[1]; + + set_atom_flag = false; +} +UnitCell::~UnitCell() +{ +} #ifdef __LCAO -InfoNonlocal::InfoNonlocal(){} -InfoNonlocal::~InfoNonlocal(){} -LCAO_Orbitals::LCAO_Orbitals(){} -LCAO_Orbitals::~LCAO_Orbitals(){} +InfoNonlocal::InfoNonlocal() +{ +} +InfoNonlocal::~InfoNonlocal() +{ +} +LCAO_Orbitals::LCAO_Orbitals() +{ +} +LCAO_Orbitals::~LCAO_Orbitals() +{ +} #endif -Magnetism::Magnetism(){} -Magnetism::~Magnetism(){} - -void Charge_Mixing::set_mixing( - const std::string &mixing_mode_in, - const double &mixing_beta_in, - const int &mixing_ndim_in, - const double &mixing_gg0_in, - const bool &mixing_tau_in -){return;} -//void Charge_Mixing::need_auto_set(){} +Magnetism::Magnetism() +{ +} +Magnetism::~Magnetism() +{ +} + +void Charge_Mixing::set_mixing(const std::string& mixing_mode_in, + const double& mixing_beta_in, + const int& mixing_ndim_in, + const double& mixing_gg0_in, + const bool& mixing_tau_in) +{ + return; +} +// void Charge_Mixing::need_auto_set(){} void Charge_Mixing::need_auto_set() { - this->autoset = true; -} -void Occupy::decision(const std::string &name,const std::string &smearing_method,const double &smearing_sigma){return;} -//void UnitCell::setup(const std::string&,const int&,const int&,const bool&,const std::string&){return;} -void UnitCell::setup(const std::string &latname_in, - const int &ntype_in, - const int &lmaxmax_in, - const bool &init_vel_in, - const std::string &fixed_axes_in) -{ - this->latName = latname_in; - this->ntype = ntype_in; - this->lmaxmax = lmaxmax_in; - this->init_vel = init_vel_in; - // pengfei Li add 2018-11-11 - if (fixed_axes_in == "None") - { - this->lc[0] = 1; - this->lc[1] = 1; - this->lc[2] = 1; - } - else if (fixed_axes_in == "volume") - { - this->lc[0] = 1; - this->lc[1] = 1; - this->lc[2] = 1; - if(!GlobalV::relax_new) - { - ModuleBase::WARNING_QUIT("Input","there are bugs in the old implementation; set relax_new to be 1 for fixed_volume relaxation"); - } - } - else if (fixed_axes_in == "shape") - { - if(!GlobalV::relax_new) - { - ModuleBase::WARNING_QUIT("Input","set relax_new to be 1 for fixed_shape relaxation"); - } - this->lc[0] = 1; - this->lc[1] = 1; - this->lc[2] = 1; - } - else if (fixed_axes_in == "a") - { - this->lc[0] = 0; - this->lc[1] = 1; - this->lc[2] = 1; - } - else if (fixed_axes_in == "b") - { - this->lc[0] = 1; - this->lc[1] = 0; - this->lc[2] = 1; - } - else if (fixed_axes_in == "c") - { - this->lc[0] = 1; - this->lc[1] = 1; - this->lc[2] = 0; - } - else if (fixed_axes_in == "ab") - { - this->lc[0] = 0; - this->lc[1] = 0; - this->lc[2] = 1; - } - else if (fixed_axes_in == "ac") - { - this->lc[0] = 0; - this->lc[1] = 1; - this->lc[2] = 0; - } - else if (fixed_axes_in == "bc") - { - this->lc[0] = 1; - this->lc[1] = 0; - this->lc[2] = 0; - } - else if (fixed_axes_in == "abc") - { - this->lc[0] = 0; - this->lc[1] = 0; - this->lc[2] = 0; - } - else - { - ModuleBase::WARNING_QUIT("Input", "fixed_axes should be None,volume,shape,a,b,c,ab,ac,bc or abc!"); - } - return; -} -void Structure_Factor::set(const int&){return;} + this->autoset = true; +} +void Occupy::decision(const std::string& name, const std::string& smearing_method, const double& smearing_sigma) +{ + return; +} +// void UnitCell::setup(const std::string&,const int&,const int&,const bool&,const std::string&){return;} +void UnitCell::setup(const std::string& latname_in, + const int& ntype_in, + const int& lmaxmax_in, + const bool& init_vel_in, + const std::string& fixed_axes_in) +{ + this->latName = latname_in; + this->ntype = ntype_in; + this->lmaxmax = lmaxmax_in; + this->init_vel = init_vel_in; + // pengfei Li add 2018-11-11 + if (fixed_axes_in == "None") + { + this->lc[0] = 1; + this->lc[1] = 1; + this->lc[2] = 1; + } + else if (fixed_axes_in == "volume") + { + this->lc[0] = 1; + this->lc[1] = 1; + this->lc[2] = 1; + if (!GlobalV::relax_new) + { + ModuleBase::WARNING_QUIT( + "Input", + "there are bugs in the old implementation; set relax_new to be 1 for fixed_volume relaxation"); + } + } + else if (fixed_axes_in == "shape") + { + if (!GlobalV::relax_new) + { + ModuleBase::WARNING_QUIT("Input", "set relax_new to be 1 for fixed_shape relaxation"); + } + this->lc[0] = 1; + this->lc[1] = 1; + this->lc[2] = 1; + } + else if (fixed_axes_in == "a") + { + this->lc[0] = 0; + this->lc[1] = 1; + this->lc[2] = 1; + } + else if (fixed_axes_in == "b") + { + this->lc[0] = 1; + this->lc[1] = 0; + this->lc[2] = 1; + } + else if (fixed_axes_in == "c") + { + this->lc[0] = 1; + this->lc[1] = 1; + this->lc[2] = 0; + } + else if (fixed_axes_in == "ab") + { + this->lc[0] = 0; + this->lc[1] = 0; + this->lc[2] = 1; + } + else if (fixed_axes_in == "ac") + { + this->lc[0] = 0; + this->lc[1] = 1; + this->lc[2] = 0; + } + else if (fixed_axes_in == "bc") + { + this->lc[0] = 1; + this->lc[1] = 0; + this->lc[2] = 0; + } + else if (fixed_axes_in == "abc") + { + this->lc[0] = 0; + this->lc[1] = 0; + this->lc[2] = 0; + } + else + { + ModuleBase::WARNING_QUIT("Input", "fixed_axes should be None,volume,shape,a,b,c,ab,ac,bc or abc!"); + } + return; +} +void Structure_Factor::set(const int&) +{ + return; +} namespace MD_func { - double current_step(const int& my_rank, const std::string& file_dir){return 0;} +double current_step(const int& my_rank, const std::string& file_dir) +{ + return 0; } +} // namespace MD_func namespace GlobalC { - UnitCell ucell; - wavefunc wf; - ModuleSymmetry::Symmetry symm; - Structure_Factor sf; - ModuleDFTU::DFTU dftu; - Restart restart; - pseudopot_cell_vnl ppcell; - Charge_Mixing CHR_MIX; -} +UnitCell ucell; +wavefunc wf; +ModuleSymmetry::Symmetry symm; +Structure_Factor sf; +ModuleDFTU::DFTU dftu; +Restart restart; +pseudopot_cell_vnl ppcell; +Charge_Mixing CHR_MIX; +} // namespace GlobalC #undef private diff --git a/source/module_io/test/input_conv_test.cpp b/source/module_io/test/input_conv_test.cpp index d58f843fd5..9d6f911d0f 100644 --- a/source/module_io/test/input_conv_test.cpp +++ b/source/module_io/test/input_conv_test.cpp @@ -124,13 +124,13 @@ TEST_F(InputConvTest, Conv) EXPECT_DOUBLE_EQ(elecstate::Gatefield::block_down,0.45); EXPECT_DOUBLE_EQ(elecstate::Gatefield::block_up,0.55); EXPECT_DOUBLE_EQ(elecstate::Gatefield::block_height,0.1); - EXPECT_EQ(ELEC_evolve::td_force_dt,0.02); - EXPECT_EQ(ELEC_evolve::td_vext,false); - EXPECT_EQ(ELEC_evolve::out_dipole,false); - EXPECT_EQ(ELEC_evolve::out_efield,false); - EXPECT_EQ(ELEC_evolve::td_print_eij,-1.0); - EXPECT_EQ(ELEC_evolve::td_edm,0); - EXPECT_EQ(GlobalV::ocp,false); + EXPECT_EQ(module_tddft::Evolve_elec::td_force_dt, 0.02); + EXPECT_EQ(module_tddft::Evolve_elec::td_vext, false); + EXPECT_EQ(module_tddft::Evolve_elec::out_dipole, false); + EXPECT_EQ(module_tddft::Evolve_elec::out_efield, false); + EXPECT_EQ(module_tddft::Evolve_elec::td_print_eij, -1.0); + EXPECT_EQ(module_tddft::Evolve_elec::td_edm, 0); + EXPECT_EQ(GlobalV::ocp,false); EXPECT_EQ(GlobalV::ocp_set,INPUT.ocp_set); EXPECT_EQ(GlobalV::out_mul,0); EXPECT_EQ(GlobalC::ppcell.cell_factor,1.2); @@ -507,9 +507,9 @@ TEST_F(InputConvTest,parse ) INPUT.Default(); std::string input_file = "./support/INPUT"; INPUT.Read(input_file); - ELEC_evolve::td_vext=true; - Input_Conv::Convert(); - EXPECT_EQ(ELEC_evolve::td_vext_dire_case.size(),0); + module_tddft::Evolve_elec::td_vext = true; + Input_Conv::Convert(); + EXPECT_EQ(module_tddft::Evolve_elec::td_vext_dire_case.size(), 0); } TEST_F(InputConvTest,parse2 ) @@ -534,4 +534,63 @@ TEST_F(InputConvTest,ParseExpressionDouble) EXPECT_DOUBLE_EQ(vec[3],0.5); } +#ifdef __LCAO +TEST_F(InputConvTest, ConvertUnitsWithEmptyParams) +{ + std::string params = ""; + double c = 2.0; + std::vector expected = {}; + std::vector result = Input_Conv::convert_units(params, c); + EXPECT_EQ(result, expected); +} + +TEST_F(InputConvTest, ConvertUnitsWithSingleParam) +{ + std::string params = "1.23"; + double c = 2.0; + std::vector expected = {2.46}; + std::vector result = Input_Conv::convert_units(params, c); + EXPECT_EQ(result, expected); +} + +TEST_F(InputConvTest, ConvertUnitsWithMultipleParams) +{ + std::string params = "1.23 4.56 7.89"; + double c = 0.5; + std::vector expected = {0.615, 2.28, 3.945}; + std::vector result = Input_Conv::convert_units(params, c); + EXPECT_EQ(result, expected); +} + +TEST_F(InputConvTest, ReadTdEfieldTest) +{ + Input_Conv::read_td_efield(); + + EXPECT_EQ(elecstate::H_TDDFT_pw::stype, 0); + EXPECT_EQ(elecstate::H_TDDFT_pw::ttype[0], 0); + EXPECT_EQ(elecstate::H_TDDFT_pw::tstart, 1); + EXPECT_EQ(elecstate::H_TDDFT_pw::tend, 1000); + EXPECT_EQ(elecstate::H_TDDFT_pw::lcut1, 0.05); + EXPECT_EQ(elecstate::H_TDDFT_pw::lcut2, 0.95); + EXPECT_NEAR(elecstate::H_TDDFT_pw::gauss_omega[0], 22.13 * 2 * ModuleBase::PI * ModuleBase::AU_to_FS, 1e-8); + EXPECT_EQ(elecstate::H_TDDFT_pw::gauss_phase[0], 0.0); + EXPECT_NEAR(elecstate::H_TDDFT_pw::gauss_sigma[0], 30.0 / ModuleBase::AU_to_FS, 1e-8); + EXPECT_EQ(elecstate::H_TDDFT_pw::gauss_t0[0], 100.0); + EXPECT_NEAR(elecstate::H_TDDFT_pw::gauss_amp[0], 0.25 * ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV, 1e-8); + EXPECT_NEAR(elecstate::H_TDDFT_pw::trape_omega[0], 1.60 * 2 * ModuleBase::PI * ModuleBase::AU_to_FS, 1e-8); + EXPECT_EQ(elecstate::H_TDDFT_pw::trape_phase[0], 0.0); + EXPECT_EQ(elecstate::H_TDDFT_pw::trape_t1[0], 1875); + EXPECT_EQ(elecstate::H_TDDFT_pw::trape_t2[0], 5625); + EXPECT_EQ(elecstate::H_TDDFT_pw::trape_t3[0], 7500); + EXPECT_NEAR(elecstate::H_TDDFT_pw::trape_amp[0], 2.74 * ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV, 1e-8); + EXPECT_NEAR(elecstate::H_TDDFT_pw::trigo_omega1[0], 1.164656 * 2 * ModuleBase::PI * ModuleBase::AU_to_FS, 1e-8); + EXPECT_NEAR(elecstate::H_TDDFT_pw::trigo_omega2[0], 0.029116 * 2 * ModuleBase::PI * ModuleBase::AU_to_FS, 1e-8); + EXPECT_EQ(elecstate::H_TDDFT_pw::trigo_phase1[0], 0.0); + EXPECT_EQ(elecstate::H_TDDFT_pw::trigo_phase2[0], 0.0); + EXPECT_NEAR(elecstate::H_TDDFT_pw::trigo_amp[0], 2.74 * ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV, 1e-8); + EXPECT_EQ(elecstate::H_TDDFT_pw::heavi_t0[0], 100); + EXPECT_NEAR(elecstate::H_TDDFT_pw::heavi_amp[0], 1.00 * ModuleBase::BOHR_TO_A / ModuleBase::Ry_to_eV, 1e-8); +} +#endif + #undef private diff --git a/source/module_io/write_dipole.cpp b/source/module_io/write_dipole.cpp index a246050958..f2884f5638 100644 --- a/source/module_io/write_dipole.cpp +++ b/source/module_io/write_dipole.cpp @@ -1,6 +1,6 @@ #include "module_base/parallel_reduce.h" #include "module_elecstate/module_charge/charge.h" -#include "module_hamilt_lcao/module_tddft/ELEC_evolve.h" +#include "module_hamilt_lcao/module_tddft/evolve_elec.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" #include "module_io/dipole_io.h" diff --git a/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref b/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref index 689ce39369..a6d0fcd531 100644 --- a/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref +++ b/tests/integrate/601_NO_TDDFT_H2_taylor/result.ref @@ -1,5 +1,5 @@ -etotref -18.06759546636300 -etotperatomref -9.0337977332 -totalforceref 46.058216 -totalstressref 81.553158 -totaltimeref +3.9590 +etotref -17.65560405098635 +etotperatomref -8.8278020255 +totalforceref 26.351130 +totalstressref 46.577654 +totaltimeref +4.0140