23#include <boost/algorithm/string.hpp>
24#include <boost/format.hpp>
47 Index ignored_corelevels = 0;
50 basis.
Load(
"corelevels");
51 Index coreElectrons = 0;
52 for (
const auto& atom :
orbitals_.QMAtoms()) {
55 ignored_corelevels = coreElectrons / 2;
57 return ignored_corelevels;
74 std::string ranges = options.
get(
".ranges").
as<std::string>();
78 if (ranges ==
"factor") {
80 double rpamaxfactor = options.
get(
".rpamax").
as<
double>();
81 rpamax =
Index(rpamaxfactor *
double(num_of_levels)) -
84 double qpminfactor = options.
get(
".qpmin").
as<
double>();
86 num_of_occlevels -
Index(qpminfactor *
double(num_of_occlevels)) - 1;
88 double qpmaxfactor = options.
get(
".qpmax").
as<
double>();
90 num_of_occlevels +
Index(qpmaxfactor *
double(num_of_occlevels)) - 1;
92 double bseminfactor = options.
get(
".bsemin").
as<
double>();
94 num_of_occlevels -
Index(bseminfactor *
double(num_of_occlevels)) - 1;
96 double bsemaxfactor = options.
get(
".bsemax").
as<
double>();
98 num_of_occlevels +
Index(bsemaxfactor *
double(num_of_occlevels)) - 1;
100 }
else if (ranges ==
"explicit") {
105 bse_vmin = options.
get(
".bsemin").
as<
Index>();
106 bse_cmax = options.
get(
".bsemax").
as<
Index>();
107 }
else if (ranges ==
"default") {
108 rpamax = num_of_levels - 1;
110 qpmax = 3 * homo + 1;
112 bse_cmax = 3 * homo + 1;
113 }
else if (ranges ==
"full") {
114 rpamax = num_of_levels - 1;
116 qpmax = num_of_levels - 1;
118 bse_cmax = num_of_levels - 1;
120 std::string ignore_corelevels =
121 options.
get(
".ignore_corelevels").
as<std::string>();
123 if (ignore_corelevels ==
"RPA" || ignore_corelevels ==
"GW" ||
124 ignore_corelevels ==
"BSE") {
126 if (ignore_corelevels ==
"RPA") {
127 rpamin = ignored_corelevels;
129 if (ignore_corelevels ==
"GW" || ignore_corelevels ==
"RPA") {
130 if (qpmin < ignored_corelevels) {
131 qpmin = ignored_corelevels;
134 if (ignore_corelevels ==
"GW" || ignore_corelevels ==
"RPA" ||
135 ignore_corelevels ==
"BSE") {
136 if (bse_vmin < ignored_corelevels) {
137 bse_vmin = ignored_corelevels;
142 <<
TimeStamp() <<
" Ignoring " << ignored_corelevels
143 <<
" core levels for " << ignore_corelevels <<
" and beyond." << flush;
147 if (rpamax > num_of_levels) {
148 rpamax = num_of_levels - 1;
150 if (qpmax > num_of_levels) {
151 qpmax = num_of_levels - 1;
153 if (bse_cmax > num_of_levels) {
154 bse_cmax = num_of_levels - 1;
179 orbitals_.setBSEindices(bse_vmin, bse_cmax);
182 Index bse_vmax = homo;
183 Index bse_cmin = homo + 1;
184 Index bse_vtotal = bse_vmax - bse_vmin + 1;
185 Index bse_ctotal = bse_cmax - bse_cmin + 1;
186 Index bse_size = bse_vtotal * bse_ctotal;
189 <<
":" << rpamax <<
"]" << flush;
191 <<
":" << qpmax <<
"]" << flush;
193 <<
TimeStamp() <<
" BSE level range occ[" << bse_vmin <<
":" << bse_vmax
194 <<
"] virt[" << bse_cmin <<
":" << bse_cmax <<
"]" << flush;
204 options.
get(
"bse.davidson.correction").
as<std::string>();
207 options.
get(
"bse.davidson.tolerance").
as<std::string>();
210 options.
get(
"bse.davidson.update").
as<std::string>();
222 Index full_bse_size = (
bseopt_.useTDA) ? bse_size : 2 * bse_size;
224 << full_bse_size <<
"x" << full_bse_size << flush;
226 bseopt_.use_Hqp_offdiag = options.
get(
"bse.use_Hqp_offdiag").
as<
bool>();
228 if (!
bseopt_.use_Hqp_offdiag) {
230 <<
" BSE without Hqp offdiagonal elements" << flush;
233 <<
" BSE with Hqp offdiagonal elements" << flush;
237 bseopt_.dyn_tolerance = options.
get(
"bse.dyn_screen_tol").
as<
double>();
238 if (
bseopt_.max_dyn_iter > 0) {
248 }
else if (options.
exists(
".auxbasisset")) {
255 }
catch (std::runtime_error&) {
257 "There is no auxbasis from the dftcalculation nor did you specify an "
258 "auxbasisset for the gwbse calculation. Also no auxiliary basisset "
263 <<
" Could not find an auxbasisset using " <<
auxbasis_name_ << flush;
266 std::string mode = options.
get(
"gw.mode").
as<std::string>();
267 if (mode ==
"G0W0") {
268 gwopt_.gw_sc_max_iterations = 1;
269 }
else if (mode ==
"evGW") {
277 gwopt_.shift = options.
get(
"gw.scissor_shift").
as<
double>();
279 options.
get(
"gw.qp_sc_limit").
as<
double>();
282 gwopt_.g_sc_max_iterations =
288 if (mode ==
"evGW") {
293 options.
get(
"gw.sc_limit").
as<
double>();
296 <<
" qp_sc_limit [Hartree]: " <<
gwopt_.g_sc_limit << flush;
297 if (
gwopt_.gw_sc_max_iterations > 1) {
299 <<
" gw_sc_limit [Hartree]: " <<
gwopt_.gw_sc_limit << flush;
301 bseopt_.min_print_weight = options.
get(
"bse.print_weight").
as<
double>();
305 std::string tasks_string = options.
get(
".tasks").
as<std::string>();
306 boost::algorithm::to_lower(tasks_string);
307 if (tasks_string.find(
"all") != std::string::npos) {
312 if (tasks_string.find(
"gw") != std::string::npos) {
315 if (tasks_string.find(
"singlets") != std::string::npos) {
318 if (tasks_string.find(
"triplets") != std::string::npos) {
337 if (options.
exists(
"bse.fragments")) {
338 std::vector<tools::Property*> prop_region =
339 options.
Select(
"bse.fragments.fragment");
342 std::string indices = prop->get(
"indices").as<std::string>();
348 gwopt_.sigma_integration =
349 options.
get(
"gw.sigma_integrator").
as<std::string>();
351 <<
" Sigma integration: " <<
gwopt_.sigma_integration << flush;
354 if (
gwopt_.sigma_integration ==
"exact") {
356 <<
" RPA Hamiltonian size: " << (homo + 1 - rpamin) * (rpamax - homo)
359 if (
gwopt_.sigma_integration ==
"cda") {
362 <<
" Quadrature integration order : " <<
gwopt_.order << flush;
363 gwopt_.quadrature_scheme =
364 options.
get(
"gw.quadrature_scheme").
as<std::string>();
366 <<
" Quadrature integration scheme : " <<
gwopt_.quadrature_scheme
368 gwopt_.alpha = options.
get(
"gw.alpha").
as<
double>();
370 <<
" Alpha smoothing parameter : " <<
gwopt_.alpha << flush;
372 gwopt_.qp_solver = options.
get(
"gw.qp_solver").
as<std::string>();
375 if (
gwopt_.qp_solver ==
"grid") {
377 gwopt_.qp_grid_spacing = options.
get(
"gw.qp_grid_spacing").
as<
double>();
379 <<
" QP grid steps: " <<
gwopt_.qp_grid_steps << flush;
381 <<
" QP grid spacing: " <<
gwopt_.qp_grid_spacing << flush;
389 gwopt_.gw_mixing_alpha = options.
get(
"gw.mixing_alpha").
as<
double>();
391 if (mode ==
"evGW") {
392 if (
gwopt_.gw_mixing_order == 0) {
394 }
else if (
gwopt_.gw_mixing_order == 1) {
396 <<
gwopt_.gw_mixing_alpha << std::flush;
399 <<
gwopt_.gw_mixing_order <<
" using alpha "
400 <<
gwopt_.gw_mixing_alpha << std::flush;
404 if (options.
exists(
"gw.sigma_plot")) {
409 options.
get(
"gw.sigma_plot.filename").
as<std::string>();
430 (format(
"%1$+1.6f ") % (
orbitals_.getDFTTotalEnergy() * hrt2ev)).str());
441 (format(
"%1$+1.6f ") %
446 (format(
"%1$+1.6f ") % (
orbitals_.QPpertEnergies()(state) * hrt2ev))
449 level_summary.
add(
"qp_energy",
450 (format(
"%1$+1.6f ") %
451 (
orbitals_.QPdiag().eigenvalues()(state) * hrt2ev))
461 "omega", (format(
"%1$+1.6f ") %
462 (
orbitals_.BSESinglets().eigenvalues()(state) * hrt2ev))
466 const Eigen::Vector3d& dipoles = (
orbitals_.TransitionDipoles())[state];
467 double f = 2 * dipoles.squaredNorm() *
468 orbitals_.BSESinglets().eigenvalues()(state) / 3.0;
470 level_summary.
add(
"f", (format(
"%1$+1.6f ") % f).str());
472 "Trdipole", (format(
"%1$+1.4f %2$+1.4f %3$+1.4f") % dipoles.x() %
473 dipoles.y() % dipoles.z())
486 "omega", (format(
"%1$+1.6f ") %
487 (
orbitals_.BSETriplets().eigenvalues()(state) * hrt2ev))
505 if (
orbitals_.getXCFunctionalName().empty()) {
509 throw std::runtime_error(
"Functionals from DFT " +
510 orbitals_.getXCFunctionalName() +
" GWBSE " +
518 <<
TimeStamp() <<
" Setup grid for integration with gridsize: " <<
grid_
519 <<
" with " << grid.
getGridSize() <<
" points, divided into "
526 Eigen::MatrixXd DMAT =
orbitals_.DensityMatrixGroundState();
533 Eigen::MatrixXd mos =
536 Eigen::MatrixXd vxc = mos.transpose() * e_vxc_ao.
matrix() * mos;
538 <<
TimeStamp() <<
" Calculated exchange-correlation expectation values "
552 <<
TimeStamp() <<
" Using MKL overload for Eigen " << flush;
556 <<
" Using native Eigen implementation, no BLAS overload " << flush;
561 <<
TimeStamp() <<
" Using CUDA support for tensor multiplication with "
562 << nogpus <<
" GPUs." << flush;
566 <<
TimeStamp() <<
" Molecule Coordinates [A] " << flush;
568 std::string output = (boost::format(
"%5d"
570 " %1.4f %1.4f %1.4f") %
571 atom.getId() % atom.getElement() %
580 std::string dft_package =
orbitals_.getQMpackage();
582 <<
TimeStamp() <<
" DFT data was created by " << dft_package << flush;
604 <<
" size:" << frag.size() << flush;
609 throw std::runtime_error(
610 "You want no GW calculation but the orb file has no QPcoefficients for "
620 <<
" Calculating Mmn_beta (3-center-repulsion x orbitals) " << flush;
621 Mmn.
Fill(auxbasis, dftbasis,
orbitals_.MOs().eigenvectors());
624 <<
" functions from Aux Coulomb matrix to avoid near linear dependencies"
627 <<
TimeStamp() <<
" Calculated Mmn_beta (3-center-repulsion x orbitals) "
633 std::chrono::time_point<std::chrono::system_clock> start =
634 std::chrono::system_clock::now();
652 <<
TimeStamp() <<
" Calculating offdiagonal part of Sigma " << flush;
655 <<
TimeStamp() <<
" Calculated offdiagonal part of Sigma " << flush;
659 Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es =
661 if (es.info() == Eigen::ComputationInfo::Success) {
663 <<
TimeStamp() <<
" Diagonalized QP Hamiltonian " << flush;
666 orbitals_.QPdiag().eigenvectors() = es.eigenvectors();
667 orbitals_.QPdiag().eigenvalues() = es.eigenvalues();
668 std::chrono::duration<double> elapsed_time =
669 std::chrono::system_clock::now() - start;
671 << elapsed_time.count() <<
" seconds." << flush;
678 throw std::runtime_error(
679 "The ranges for GW and RPA do not agree with the ranges from the "
680 ".orb file, rerun your GW calculation");
682 const Eigen::MatrixXd& qpcoeff =
orbitals_.QPdiag().eigenvectors();
684 Hqp = qpcoeff *
orbitals_.QPdiag().eigenvalues().asDiagonal() *
691 std::chrono::time_point<std::chrono::system_clock> start =
692 std::chrono::system_clock::now();
698 Eigen::VectorXd Hd_static_contrib_triplet;
699 Eigen::VectorXd Hd_static_contrib_singlet;
704 <<
TimeStamp() <<
" Solved BSE for triplets " << flush;
711 <<
TimeStamp() <<
" Solved BSE for singlets " << flush;
729 std::chrono::duration<double> elapsed_time =
730 std::chrono::system_clock::now() - start;
732 << elapsed_time.count() <<
" seconds." << flush;
735 <<
TimeStamp() <<
" GWBSE calculation finished " << flush;
Container to hold Basisfunctions for all atoms.
Index AOBasisSize() const
void Solve_singlets(Orbitals &orb) const
void Analyze_singlets(std::vector< QMFragment< BSE_Population > > fragments, const Orbitals &orb) const
void Perturbative_DynamicalScreening(const QMStateType &type, Orbitals &orb)
void Analyze_triplets(std::vector< QMFragment< BSE_Population > > fragments, const Orbitals &orb) const
void configure(const options &opt, const Eigen::VectorXd &RPAEnergies, const Eigen::MatrixXd &Hqp_in)
void Solve_triplets(Orbitals &orb) const
void Load(const std::string &name)
void Load(const std::string &name)
const ECPElement & getElement(std::string element_type) const
std::vector< QMFragment< BSE_Population > > fragments_
std::string sigma_plot_states_
void addoutput(tools::Property &summary)
std::string auxbasis_name_
Eigen::MatrixXd CalculateVXC(const AOBasis &dftbasis)
std::string dftbasis_name_
void Initialize(tools::Property &options)
double sigma_plot_spacing_
std::string sigma_plot_filename_
bool do_dynamical_screening_bse_
void PlotSigma(std::string filename, Index steps, double spacing, std::string states) const
Eigen::SelfAdjointEigenSolver< Eigen::MatrixXd > DiagonalizeQPHamiltonian() const
void configure(const options &opt)
Eigen::MatrixXd getHQP() const
Eigen::VectorXd getGWAResults() const
Eigen::VectorXd RPAInputEnergies() const
void CalculateGWPerturbation()
Eigen::MatrixXd & matrix()
void Fill(const AOBasis &auxbasis, const AOBasis &dftbasis, const Eigen::MatrixXd &dft_orbitals)
void Initialize(Index basissize, Index mmin, Index mmax, Index nmin, Index nmax)
Index Removedfunctions() const
Timestamp returns the current time as a string Example: cout << TimeStamp()
Index getGridSize() const
void GridSetup(const std::string &type, const QMMolecule &atoms, const AOBasis &basis)
Index getBoxesSize() const
Mat_p_Energy IntegrateVXC(const Eigen::MatrixXd &density_matrix) const
void setXCfunctional(const std::string &functional)
#define XTP_LOG(level, log)
Charge transport classes.
bool XTP_HAS_MKL_OVERLOAD()
base class for all analysis tools