Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
161 changes: 74 additions & 87 deletions source/source_esolver/esolver_ks_lcao_tddft.cpp
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
#include "esolver_ks_lcao_tddft.h"

#include "source_estate/elecstate_tools.h"
#include "source_io/cal_r_overlap_R.h"
#include "source_io/dipole_io.h"
#include "source_io/td_current_io.h"
#include "source_io/write_HS.h"
#include "source_io/write_HS_R.h"
#include "source_estate/elecstate_tools.h"

//--------------temporary----------------------------
#include "source_base/blas_connector.h"
Expand All @@ -17,17 +17,17 @@
#include "source_estate/module_dm/cal_edm_tddft.h"
#include "source_estate/module_dm/density_matrix.h"
#include "source_estate/occupy.h"
#include "source_io/print_info.h"
#include "source_lcao/module_rt/evolve_elec.h"
#include "source_lcao/module_rt/td_velocity.h"
#include "source_pw/module_pwdft/global.h"
#include "source_io/print_info.h"

//-----HSolver ElecState Hamilt--------
#include "module_parameter/parameter.h"
#include "source_estate/cal_ux.h"
#include "source_estate/elecstate_lcao.h"
#include "source_lcao/hamilt_lcaodft/hamilt_lcao.h"
#include "source_hsolver/hsolver_lcao.h"
#include "module_parameter/parameter.h"
#include "source_lcao/hamilt_lcaodft/hamilt_lcao.h"
#include "source_psi/psi.h"

//-----force& stress-------------------
Expand Down Expand Up @@ -87,52 +87,52 @@ void ESolver_KS_LCAO_TDDFT<Device>::before_all_runners(UnitCell& ucell, const In

template <typename Device>
void ESolver_KS_LCAO_TDDFT<Device>::hamilt2rho_single(UnitCell& ucell,
const int istep,
const int iter,
const double ethr)
const int istep,
const int iter,
const double ethr)
{
if (PARAM.inp.init_wfc == "file")
{
if (istep >= 1)
{
module_rt::Evolve_elec<Device>::solve_psi(istep,
PARAM.inp.nbands,
PARAM.globalv.nlocal,
kv.get_nks(),
this->p_hamilt,
this->pv,
this->psi,
this->psi_laststep,
this->Hk_laststep,
this->Sk_laststep,
this->pelec->ekb,
GlobalV::ofs_running,
td_htype,
PARAM.inp.propagator,
use_tensor,
use_lapack);
PARAM.inp.nbands,
PARAM.globalv.nlocal,
kv.get_nks(),
this->p_hamilt,
this->pv,
this->psi,
this->psi_laststep,
this->Hk_laststep,
this->Sk_laststep,
this->pelec->ekb,
GlobalV::ofs_running,
td_htype,
PARAM.inp.propagator,
use_tensor,
use_lapack);
this->weight_dm_rho();
}
this->weight_dm_rho();
}
else if (istep >= 2)
{
module_rt::Evolve_elec<Device>::solve_psi(istep,
PARAM.inp.nbands,
PARAM.globalv.nlocal,
kv.get_nks(),
this->p_hamilt,
this->pv,
this->psi,
this->psi_laststep,
this->Hk_laststep,
this->Sk_laststep,
this->pelec->ekb,
GlobalV::ofs_running,
td_htype,
PARAM.inp.propagator,
use_tensor,
use_lapack);
PARAM.inp.nbands,
PARAM.globalv.nlocal,
kv.get_nks(),
this->p_hamilt,
this->pv,
this->psi,
this->psi_laststep,
this->Hk_laststep,
this->Sk_laststep,
this->pelec->ekb,
GlobalV::ofs_running,
td_htype,
PARAM.inp.propagator,
use_tensor,
use_lapack);
this->weight_dm_rho();
}
else
Expand Down Expand Up @@ -163,18 +163,14 @@ void ESolver_KS_LCAO_TDDFT<Device>::hamilt2rho_single(UnitCell& ucell,
}

template <typename Device>
void ESolver_KS_LCAO_TDDFT<Device>::iter_finish(
UnitCell& ucell,
const int istep,
int& iter,
bool& conv_esolver)
void ESolver_KS_LCAO_TDDFT<Device>::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver)
{
// print occupation of each band
if (iter == 1 && istep <= 2)
{
GlobalV::ofs_running << " ---------------------------------------------------------"
<< std::endl;
GlobalV::ofs_running << " occupations of electrons" << std::endl;
GlobalV::ofs_running << " ----------------------------------------------------------" << std::endl;
GlobalV::ofs_running << " Occupations of electrons" << std::endl;
GlobalV::ofs_running << " ----------------------------------------------------------" << std::endl;
GlobalV::ofs_running << " k-point state occupation" << std::endl;
GlobalV::ofs_running << std::setiosflags(std::ios::showpoint);
GlobalV::ofs_running << std::left;
Expand All @@ -183,23 +179,21 @@ void ESolver_KS_LCAO_TDDFT<Device>::iter_finish(
{
for (int ib = 0; ib < PARAM.inp.nbands; ib++)
{
GlobalV::ofs_running << " " << std::setw(9)
<< ik+1 << std::setw(8) << ib + 1
<< std::setw(12) << this->pelec->wg(ik, ib) << std::endl;
GlobalV::ofs_running << " " << std::setw(9) << ik + 1 << std::setw(8) << ib + 1 << std::setw(12)
<< this->pelec->wg(ik, ib) << std::endl;
}
}
GlobalV::ofs_running << " ---------------------------------------------------------"
<< std::endl;
GlobalV::ofs_running << " ----------------------------------------------------------" << std::endl;
}

ESolver_KS_LCAO<std::complex<double>, double>::iter_finish(ucell, istep, iter, conv_esolver);
}

template <typename Device>
void ESolver_KS_LCAO_TDDFT<Device>::update_pot(UnitCell& ucell,
const int istep,
const int iter,
const bool conv_esolver)
void ESolver_KS_LCAO_TDDFT<Device>::update_pot(UnitCell& ucell,
const int istep,
const int iter,
const bool conv_esolver)
{
// Calculate new potential according to new Charge Density
if (!conv_esolver)
Expand Down Expand Up @@ -234,7 +228,6 @@ void ESolver_KS_LCAO_TDDFT<Device>::update_pot(UnitCell& ucell,
nrow_tmp = nlocal;
#endif
this->psi_laststep = new psi::Psi<std::complex<double>>(kv.get_nks(), ncol_tmp, nrow_tmp, kv.ngk, true);

}

// allocate memory for Hk_laststep and Sk_laststep
Expand Down Expand Up @@ -282,8 +275,8 @@ void ESolver_KS_LCAO_TDDFT<Device>::update_pot(UnitCell& ucell,
if (td_htype == 1)
{
this->p_hamilt->updateHk(ik);
hamilt::MatrixBlock <std::complex<double>> h_mat;
hamilt::MatrixBlock <std::complex<double>> s_mat;
hamilt::MatrixBlock<std::complex<double>> h_mat;
hamilt::MatrixBlock<std::complex<double>> s_mat;
this->p_hamilt->matrix(h_mat, s_mat);

if (use_tensor && use_lapack)
Expand Down Expand Up @@ -323,31 +316,31 @@ void ESolver_KS_LCAO_TDDFT<Device>::update_pot(UnitCell& ucell,
}

// print "eigen value" for tddft
// it seems uncessary to print out E_ii because the band energies are printed
/*
if (conv_esolver)
{
GlobalV::ofs_running << "----------------------------------------------------------"
<< std::endl;
GlobalV::ofs_running << " Print E=<psi_i|H|psi_i> " << std::endl;
GlobalV::ofs_running << " k-point state energy (eV)" << std::endl;
GlobalV::ofs_running << "----------------------------------------------------------"
<< std::endl;
GlobalV::ofs_running << std::setprecision(6);
GlobalV::ofs_running << std::setiosflags(std::ios::showpoint);

for (int ik = 0; ik < kv.get_nks(); ik++)
// it seems unnecessary to print out E_ii because the band energies are printed
/*
if (conv_esolver)
{
for (int ib = 0; ib < PARAM.inp.nbands; ib++)
GlobalV::ofs_running << "----------------------------------------------------------"
<< std::endl;
GlobalV::ofs_running << " Print E=<psi_i|H|psi_i> " << std::endl;
GlobalV::ofs_running << " k-point state energy (eV)" << std::endl;
GlobalV::ofs_running << "----------------------------------------------------------"
<< std::endl;
GlobalV::ofs_running << std::setprecision(6);
GlobalV::ofs_running << std::setiosflags(std::ios::showpoint);

for (int ik = 0; ik < kv.get_nks(); ik++)
{
GlobalV::ofs_running << " " << std::setw(7) << ik + 1
<< std::setw(7) << ib + 1
<< std::setw(10) << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV
<< std::endl;
for (int ib = 0; ib < PARAM.inp.nbands; ib++)
{
GlobalV::ofs_running << " " << std::setw(7) << ik + 1
<< std::setw(7) << ib + 1
<< std::setw(10) << this->pelec->ekb(ik, ib) * ModuleBase::Ry_to_eV
<< std::endl;
}
}
}
}
*/
*/
}

template <typename Device>
Expand All @@ -365,16 +358,11 @@ void ESolver_KS_LCAO_TDDFT<Device>::after_scf(UnitCell& ucell, const int istep,
{
std::stringstream ss_dipole;
ss_dipole << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_DIPOLE";
ModuleIO::write_dipole(ucell,
this->chr.rho_save[is],
this->chr.rhopw,
is,
istep,
ss_dipole.str());
ModuleIO::write_dipole(ucell, this->chr.rho_save[is], this->chr.rhopw, is, istep, ss_dipole.str());
}
}

// (2) write current information
// (2) write current information
if (TD_Velocity::out_current == true)
{
elecstate::DensityMatrix<std::complex<double>, double>* tmp_DM
Expand All @@ -392,7 +380,6 @@ void ESolver_KS_LCAO_TDDFT<Device>::after_scf(UnitCell& ucell, const int istep,
this->RA);
}


ModuleBase::timer::tick("ESolver_LCAO_TDDFT", "after_scf");
}

Expand All @@ -410,7 +397,7 @@ void ESolver_KS_LCAO_TDDFT<Device>::weight_dm_rho()
}

// calculate Eband energy
elecstate::calEBand(this->pelec->ekb,this->pelec->wg,this->pelec->f_en);
elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en);

// calculate the density matrix
ModuleBase::GlobalFunc::NOTE("Calculate the density matrix.");
Expand Down
5 changes: 2 additions & 3 deletions source/source_estate/module_dm/cal_edm_tddft.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
#include "source_base/scalapack_connector.h"
namespace elecstate
{

// use the original formula (Hamiltonian matrix) to calculate energy density matrix
void cal_edm_tddft(Parallel_Orbitals& pv,
elecstate::ElecState* pelec,
Expand Down Expand Up @@ -254,5 +253,5 @@ void cal_edm_tddft(Parallel_Orbitals& pv,
#endif
}
return;
}
} // namespace ModuleESolver
} // cal_edm_tddft
} // namespace elecstate
18 changes: 13 additions & 5 deletions source/source_lcao/module_rt/evolve_elec.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -83,8 +83,8 @@ void Evolve_elec<Device>::solve_psi(const int& istep,
propagator,
ofs_running,
print_matrix);
// std::cout << "Print ekb: " << std::endl;
// ekb.print(std::cout);
// GlobalV::ofs_running << "Print ekb: " << std::endl;
// ekb.print(GlobalV::ofs_running);
}
else
{
Expand Down Expand Up @@ -118,7 +118,7 @@ void Evolve_elec<Device>::solve_psi(const int& istep,
#ifdef __MPI
// Access the rank of the calling process in the communicator
int myid = 0;
int root_proc = 0;
const int root_proc = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &myid);

// Gather psi to the root process
Expand Down Expand Up @@ -203,8 +203,16 @@ void Evolve_elec<Device>::solve_psi(const int& istep,
len_HS_laststep);
syncmem_double_d2h_op()(&(ekb(ik, 0)), ekb_tensor.data<double>(), nband);

// std::cout << "Print ekb tensor: " << std::endl;
// ekb.print(std::cout);
#ifdef __MPI
const int root_proc = 0;
if (use_lapack)
{
// Synchronize ekb to all MPI processes
MPI_Bcast(&(ekb(ik, 0)), nband, MPI_DOUBLE, root_proc, MPI_COMM_WORLD);
}
#endif
// GlobalV::ofs_running << "Print ekb: " << std::endl;
// ekb.print(GlobalV::ofs_running);
}
}
else
Expand Down
2 changes: 1 addition & 1 deletion tests/15_rtTDDFT_GPU/CASES_GPU.txt
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
04_NO_CO_ocp_TDDFT_GPU
05_NO_cur_TDDFT_GPU
06_NO_dir_TDDFT_GPU
# 07_NO_EDM_TDDFT_GPU
07_NO_EDM_TDDFT_GPU
09_NO_HEAV_TDDFT_GPU
10_NO_HHG_TDDFT_GPU
11_NO_O3_TDDFT_GPU
Expand Down
Loading