Skip to content

Commit

Permalink
key: ULR Hamilt & tear down HSolverLR
Browse files Browse the repository at this point in the history
  • Loading branch information
@maki49
maki49 committed Oct 17, 2024
1 parent afb1e4a commit 8bb65ea
Show file tree
Hide file tree
Showing 16 changed files with 451 additions and 261 deletions.
1 change: 0 additions & 1 deletion source/module_lr/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,6 @@ if(ENABLE_LCAO)
operator_casida/operator_lr_hxc.cpp
operator_casida/operator_lr_exx.cpp
potentials/pot_hxc_lrtd.cpp
hsolver_lrtd.cpp
lr_spectrum.cpp
hamilt_casida.cpp
esolver_lrtd_lcao.cpp)
Expand Down
81 changes: 50 additions & 31 deletions source/module_lr/esolver_lrtd_lcao.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,9 @@
#include "utils/gint_move.hpp"
#include "utils/lr_util.h"
#include "hamilt_casida.h"
#include "hamilt_ulr.hpp"
#include "module_lr/potentials/pot_hxc_lrtd.h"
#include "module_lr/hsolver_lrtd.h"
#include "module_lr/hsolver_lrtd.hpp"
#include "module_lr/lr_spectrum.h"
#include <memory>
#include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
Expand Down Expand Up @@ -141,7 +142,7 @@ void LR::ESolver_LR<T, TR>::reset_dim_spin2()
// npairs = { nocc[0] * nvirt[0], nocc[1] * nvirt[1] };
std::cout << "** Solve the spin-up and spin-down states separately for open-shell system. **" << std::endl;
}
if (nstates > std::min(npairs[0], npairs[1]) * nk) { throw std::invalid_argument("ESolver_LR: nstates > nocc*nvirt*nk"); }
if (nstates > (npairs[0] + npairs[1]) * nk) { throw std::invalid_argument("ESolver_LR: nstates > nocc*nvirt*nk"); }
}

template <typename T, typename TR>
Expand Down Expand Up @@ -431,33 +432,63 @@ void LR::ESolver_LR<T, TR>::runner(int istep, UnitCell& cell)
this->setup_eigenvectors_X();
this->pelec->ekb.create(nspin, this->nstates);

auto efile = [&](const std::string& label)->std::string {return PARAM.globalv.global_out_dir + "Excitation_Energy_" + label + ".dat";};
auto vfile = [&](const std::string& label)->std::string {return PARAM.globalv.global_out_dir + "Excitation_Amplitude_" + label + "_" + std::to_string(GlobalV::MY_RANK) + ".dat";};
if (this->input.lr_solver != "spectrum")
{
auto write_states = [&](const std::string& label, const Real<T>* e, const T* v, const int& dim, const int& nst, const int& prec = 8)->void
{
if (GlobalV::MY_RANK == 0) { assert(nst == LR_Util::write_value(efile(label), prec, e, nst)); }
assert(nst * dim == LR_Util::write_value(vfile(label), prec, v, nst, dim));
};
// allocate and initialize A matrix and density matrix
auto spin_types = (nspin == 2 && !openshell) ? std::vector<std::string>({ "singlet", "triplet" }) : std::vector<std::string>({ "updown" });
for (int is = 0;is < X.size();++is)
if (openshell)
{
if (nspin == 2) { std::cout << "Calculating " << spin_types[is] << " excitations" << std::endl; }
HamiltLR<T> hlr(xc_kernel, nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, orb_cutoff_, GlobalC::GridD, this->psi_ks, this->eig_ks,
std::cout << "Solving spin-conserving excitation for open-shell system." << std::endl;
HamiltULR<T> hulr(xc_kernel, nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, orb_cutoff_, GlobalC::GridD, *this->psi_ks, this->eig_ks,
#ifdef __EXX
this->exx_lri, this->exx_info.info_global.hybrid_alpha,
#endif
this->gint_, this->pot[is], this->kv, this->paraX_, this->paraC_, this->paraMat_,
spin_types[is], input.ri_hartree_benchmark, (input.ri_hartree_benchmark == "aims" ? input.aims_nbasis : std::vector<int>({})));
// solve the Casida equation
HSolverLR<T> hsol(nk, this->npairs[is], is, std::max(1e-13, this->input.lr_thr), this->input.out_wfc_lr);
hsol.solve(hlr, this->X[is].template data<T>(), nloc_per_band, nstates, this->pelec->ekb, this->input.lr_solver/*,
!std::set<std::string>({ "hf", "hse" }).count(this->xc_kernel)*/); //whether the kernel is Hermitian
this->gint_, this->pot, this->kv, this->paraX_, this->paraC_, this->paraMat_);
LR::HSolver::solve(hulr, this->X[0].template data<T>(), nloc_per_band, nstates, this->pelec->ekb.c, this->input.lr_solver, this->input.lr_thr);
if (input.out_wfc_lr) { write_states("openshell", this->pelec->ekb.c, this->X[0].template data<T>(), nloc_per_band, nstates); }
}
else
{
auto spin_types = std::vector<std::string>({ "singlet", "triplet" });
for (int is = 0;is < nspin;++is)
{
std::cout << "Calculating " << spin_types[is] << " excitations" << std::endl;
HamiltLR<T> hlr(xc_kernel, nspin, this->nbasis, this->nocc, this->nvirt, this->ucell, orb_cutoff_, GlobalC::GridD, *this->psi_ks, this->eig_ks,
#ifdef __EXX
this->exx_lri, this->exx_info.info_global.hybrid_alpha,
#endif
this->gint_, this->pot[is], this->kv, this->paraX_, this->paraC_, this->paraMat_,
spin_types[is], input.ri_hartree_benchmark, (input.ri_hartree_benchmark == "aims" ? input.aims_nbasis : std::vector<int>({})));
// solve the Casida equation
LR::HSolver::solve(hlr, this->X[is].template data<T>(), nloc_per_band, nstates,
this->pelec->ekb.c + is * nstates, this->input.lr_solver, this->input.lr_thr/*,
!std::set<std::string>({ "hf", "hse" }).count(this->xc_kernel)*/); //whether the kernel is Hermitian
if (input.out_wfc_lr) { write_states(spin_types[is], this->pelec->ekb.c + is * nstates, this->X[is].template data<T>(), nloc_per_band, nstates); }
}
}
}
else // read the eigenvalues
{
std::ifstream ifs(PARAM.globalv.global_readin_dir + "Excitation_Energy.dat");
std::cout << "reading the excitation energies from file: \n";
for (int is = 0;is < nspin;++is)
auto read_states = [&](const std::string& label, Real<T>* e, T* v, const int& dim, const int& nst)->void
{
if (GlobalV::MY_RANK == 0) { assert(nst == LR_Util::read_value(efile(label), e, nst)); }
assert(nst * dim == LR_Util::read_value(vfile(label), v, nst, dim));
};
std::cout << "reading the excitation amplitudes from file: \n";
if (openshell)
{
for (int i = 0;i < this->nstates;++i) { ifs >> this->pelec->ekb(is, i); }
for (int i = 0;i < this->nstates;++i) { std::cout << this->pelec->ekb(is, i) << " "; }
read_states("openshell", this->pelec->ekb.c, this->X[0].template data<T>(), nloc_per_band, nstates);
}
else
{
auto spin_types = std::vector<std::string>({ "singlet", "triplet" });
for (int is = 0;is < nspin;++is) { read_states(spin_types[is], this->pelec->ekb.c + is * nstates, this->X[is].template data<T>(), nloc_per_band, nstates); }
}
}
return;
Expand Down Expand Up @@ -513,19 +544,7 @@ void LR::ESolver_LR<T, TR>::setup_eigenvectors_X()

auto spin_types = (nspin == 2 && !openshell) ? std::vector<std::string>({ "singlet", "triplet" }) : std::vector<std::string>({ "updown" });
// if spectrum-only, read the LR-eigenstates from file and return
if (this->input.lr_solver == "spectrum")
{
std::cout << "reading the excitation amplitudes from file: \n";
for (int is = 0; is < this->X.size(); ++is)
{
std::ifstream ifs(PARAM.globalv.global_readin_dir + "Excitation_Amplitude_" + spin_types[is] + "_" + std::to_string(GlobalV::MY_RANK) + ".dat");
assert(nstates * nloc_per_band == LR_Util::read_value(ifs, X[is].data<T>(), nstates, nloc_per_band));
}
}
else
{
set_X_initial_guess();
}
if (this->input.lr_solver != "spectrum") { set_X_initial_guess(); }
}

template<typename T, typename TR>
Expand Down Expand Up @@ -564,7 +583,7 @@ void LR::ESolver_LR<T, TR>::set_X_initial_guess()
const int ik = ib / np;
const int xstart_b = ib * nloc_per_band; //start index of band ib
const int xstart_bs = (openshell && is == 1) ? xstart_b + nk * paraX_[0].get_local_size() : xstart_b; // start index of band ib, spin is
const int is_in_x = openshell ? 1 : is; // if openshell, spin-up and spin-down are put together
const int is_in_x = openshell ? 0 : is; // if openshell, spin-up and spin-down are put together
if (px.in_this_processor(virt_global, occ_global))
{
const int ipair_loc = px.global2local_col(occ_global) * px.get_row_size() + px.global2local_row(virt_global);
Expand Down
18 changes: 6 additions & 12 deletions source/module_lr/hamilt_casida.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,20 +8,21 @@ namespace LR
const int no = this->nocc[0];
const int nv = this->nvirt[0];
const auto& px = this->pX[0];
const int nloc_per_band = nk * px.get_local_size();
int npairs = no * nv;
std::vector<T> Amat_full(this->nk * npairs * this->nk * npairs, 0.0);
for (int ik = 0;ik < this->nk;++ik)
for (int j = 0;j < no;++j)
for (int b = 0;b < nv;++b)
{//calculate A^{ai} for each bj
int bj = j * nv + b;
int kbj = ik * npairs + bj;
int bj = j * nv + b; //global
int kbj = ik * npairs + bj; //global
psi::Psi<T> X_bj(1, 1, this->nk * px.get_local_size()); // k1-first, like in iterative solver
X_bj.zero_out();
// X_bj(0, 0, lj * px.get_row_size() + lb) = this->one();
int lj = px.global2local_col(j);
int lb = px.global2local_row(b);
if (px.in_this_processor(b, j)) X_bj(0, 0, ik * px.get_local_size() + lj * px.get_row_size() + lb) = this->one();
if (px.in_this_processor(b, j)) { X_bj(0, 0, ik * px.get_local_size() + lj * px.get_row_size() + lb) = this->one(); }
psi::Psi<T> A_aibj(1, 1, this->nk * px.get_local_size()); // k1-first
A_aibj.zero_out();

Expand All @@ -41,15 +42,8 @@ namespace LR
#endif
}
// output Amat
std::cout << "Amat_full:" << std::endl;
for (int i = 0;i < this->nk * npairs;++i)
{
for (int j = 0;j < this->nk * npairs;++j)
{
std::cout << Amat_full[i * this->nk * npairs + j] << " ";
}
std::cout << std::endl;
}
std::cout << "Full A matrix:" << std::endl;
LR_Util::print_value(Amat_full.data(), nk * npairs, nk * npairs);
return Amat_full;
}

Expand Down
34 changes: 7 additions & 27 deletions source/module_lr/hamilt_casida.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ namespace LR
const UnitCell& ucell_in,
const std::vector<double>& orb_cutoff,
Grid_Driver& gd_in,
const psi::Psi<T>* psi_ks_in,
const psi::Psi<T>& psi_ks_in,
const ModuleBase::matrix& eig_ks,
#ifdef __EXX
std::weak_ptr<Exx_LRI<T>> exx_lri_in,
Expand All @@ -38,9 +38,7 @@ namespace LR
const Parallel_Orbitals& pmat_in,
const std::string& spin_type,
const std::string& ri_hartree_benchmark = "none",
const std::vector<int>& aims_nbasis = {}) : nspin(nspin), nocc(nocc), nvirt(nvirt), pX(pX_in),
nk(kv_in.get_nks() / nspin), openshell(spin_type == "up" || spin_type == "down"),
nloc_per_band(nk* (openshell ? pX_in[0].get_local_size() + pX_in[1].get_local_size() : pX_in[0].get_local_size()))
const std::vector<int>& aims_nbasis = {}) : nspin(nspin), nocc(nocc), nvirt(nvirt), pX(pX_in), nk(kv_in.get_nks() / nspin)
{
ModuleBase::TITLE("HamiltLR", "HamiltLR");
if (ri_hartree_benchmark != "aims") { assert(aims_nbasis.empty()); }
Expand Down Expand Up @@ -78,7 +76,7 @@ namespace LR
}
if (!std::set<std::string>({ "rpa", "hf" }).count(xc_kernel)) { throw std::runtime_error("ri_hartree_benchmark is only supported for xc_kernel rpa and hf"); }
RI_Benchmark::OperatorRIHartree<T>* ri_hartree_op
= new RI_Benchmark::OperatorRIHartree<T>(ucell_in, naos, nocc[0], nvirt[0], *psi_ks_in,
= new RI_Benchmark::OperatorRIHartree<T>(ucell_in, naos, nocc[0], nvirt[0], psi_ks_in,
Cs_read, Vs_read, ri_hartree_benchmark == "aims", aims_nbasis);
this->ops->add(ri_hartree_op);
}
Expand Down Expand Up @@ -124,28 +122,12 @@ namespace LR

virtual void hPsi(const T* psi_in, T* hpsi, const int ld_psi, const int& nband) const
{
assert(ld_psi == this->nloc_per_band);
assert(ld_psi == nk * pX[0].get_local_size());
hamilt::Operator<T>* node(this->ops);
if (openshell)
while (node != nullptr)
{
/// band-wise act (also works for close-shell, but not efficient)
for (int ib = 0;ib < nband;++ib)
{
const int offset = ib * ld_psi;
while (node != nullptr)
{
node->act(1, ld_psi, /*npol=*/1, psi_in + offset, hpsi + offset);
node = (hamilt::Operator<T>*)(node->next_op);
}
}
}
else
{
while (node != nullptr)
{
node->act(nband, ld_psi, /*npol=*/1, psi_in, hpsi);
node = (hamilt::Operator<T>*)(node->next_op);
}
node->act(nband, ld_psi, /*npol=*/1, psi_in, hpsi);
node = (hamilt::Operator<T>*)(node->next_op);
}
}

Expand All @@ -154,8 +136,6 @@ namespace LR
const std::vector<int>& nvirt;
const int nspin = 1;
const int nk = 1;
const bool openshell = false;
const int nloc_per_band = 1;
const std::vector<Parallel_2D>& pX;
T one()const;
/// transition density matrix in AO representation
Expand Down
Loading

0 comments on commit 8bb65ea

Please sign in to comment.