dftengine.h

Go to the documentation of this file.

/*
 *            Copyright 2009-2023 The VOTCA Development Team
 *                       (http://www.votca.org)
 *
 *      Licensed under the Apache License, Version 2.0 (the "License")
 *
 * You may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *              http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
 
#pragma once
#ifndef VOTCA_XTP_DFTENGINE_H
#define VOTCA_XTP_DFTENGINE_H
 
// Standard includes
#include <map>
 
// VOTCA includes
#include <votca/tools/property.h>
 
// Local VOTCA includes
#include "ERIs.h"
#include "convergenceacc.h"
#include "dfttimings.h"
#include "hirshfeldpartition.h"
#include "uks_convergenceacc.h"
 
#include "ecpaobasis.h"
#include "extended_hueckel.h"
#include "logger.h"
#include "qmmolecule.h"
#include "staticsite.h"
#include "vxc_grid.h"
#include "vxc_potential.h"
 
namespace votca {
namespace xtp {
class Orbitals;
struct DFTSetupCache {
  std::string eris_key;
  bool has_eris = false;
  ERIs eris;
  std::string ewald_key;
  Eigen::MatrixXd ewald_matrix;
  // AO potential of the permanent moments of the external sites; the induced
  // dipoles change between QM/MM iterations and are integrated every run.
  std::string multipole_key;
  Eigen::MatrixXd multipole_sites;  // one row per site: pos, rank, Q
  Eigen::MatrixXd multipole_matrix;
};
 
class DFTEngineTestAccess;

bool HasLibint2DerivativeSupport();

 
class DFTEngine {
 public:
  void Initialize(tools::Property& options);

  void setLogger(Logger* pLog) { pLog_ = pLog; }

  void setExternalcharges(
      std::vector<std::unique_ptr<StaticSite> >* externalsites) {
    externalsites_ = externalsites;
  }
 
  void setEwaldgrid(const Vxc_Grid& ewaldgrid) {
    external_ewaldgrid_ = ewaldgrid;
    has_ewaldgrid_ = true;
  }
 
  // The nuclei's share of the external Ewald potential, sum_A Z_A
  // phi(R_A). Added to E0 alongside the grid's contribution to H0 -- the
  // grid integrates against the DENSITY and so reaches the electrons
  // only, exactly as IntegrateExternalMultipoles pairs its AO matrix with
  // ExternalRepulsion for the multipole route.
  void setEwaldNuclearEnergy(double energy) { ewald_nuclear_energy_ = energy; }

  void setWarmStart(bool warm_start) { warm_start_ = warm_start; }

  void setSetupCache(DFTSetupCache* cache) { setup_cache_ = cache; }

  void setSCFToleranceFloor(double energy, double error);
  void ReturnSetupCache();

  bool Evaluate(Orbitals& orb);

  bool RunCDFT(Orbitals& orb, HirshfeldPartition::Constraint& constraint);

  bool EvaluateActiveRegion(Orbitals& orb);
  bool EvaluateTruncatedActiveRegion(Orbitals& trunc_orb);

  std::string getDFTBasisName() const { return dftbasis_name_; };

  std::vector<Eigen::MatrixXd> StaticResponse(
      const Orbitals& orb, const std::vector<Eigen::MatrixXd>& perturbations,
      double tolerance = 1e-7, Index max_iterations = 100);

  const std::string& getGridName() const { return grid_name_; }

  void setExternalPotentialMatrix(const Eigen::MatrixXd& v) {
    ext_potential_matrix_ = v;
  }

  std::vector<HirshfeldPartition::AtomicReference> HirshfeldReferences(
      const QMMolecule& mol) const;
  struct SpinDensity {
    Eigen::MatrixXd alpha;
    Eigen::MatrixXd beta;

    Eigen::MatrixXd total() const { return alpha + beta; }

    Eigen::MatrixXd spin() const { return alpha - beta; }
  };

  bool IsRestrictedOpenShell() const {
    return num_alpha_electrons_ != num_beta_electrons_;
  }

  Index NumberOfRestrictedOccupiedOrbitals() const {
    return std::max(num_alpha_electrons_, num_beta_electrons_);
  }

  SpinDensity BuildSpinDensity(const tools::EigenSystem& MOs) const;

  ConvergenceAcc::options BuildConvergenceOptions() const;
 
  // development path RKS vs UKS
  void setForceUKSPath(bool force) { force_uks_path_ = force; }
 
 private:
  friend class DFTEngineTestAccess;

  void Prepare(Orbitals& orb, Index numofelectrons = -1);

  Vxc_Potential<Vxc_Grid> SetupVxc(const QMMolecule& mol);

  UKSConvergenceAcc::SpinFock BuildUKSFock(
      const Eigen::MatrixXd& H0, const Vxc_Potential<Vxc_Grid>& vxcpotential,
      const Eigen::MatrixXd& alpha_new, const Eigen::MatrixXd& beta_new);
  Eigen::MatrixXd RKSOrbitalHessian(const Vxc_Potential<Vxc_Grid>& vxc,
                                    const tools::EigenSystem& mos,
                                    const Eigen::MatrixXd& F,
                                    const Eigen::MatrixXd& X);
  Eigen::MatrixXd UKSOrbitalHessian(const Vxc_Potential<Vxc_Grid>& vxc,
                                    const tools::EigenSystem& mos_alpha,
                                    const tools::EigenSystem& mos_beta,
                                    const Eigen::MatrixXd& F_alpha,
                                    const Eigen::MatrixXd& F_beta,
                                    const Eigen::MatrixXd& X);
 
  std::vector<Eigen::MatrixXd> CoupledResponse(
      const Vxc_Potential<Vxc_Grid>& vxc, const Eigen::MatrixXd& P0,
      const Eigen::MatrixXd& C_occ, const Eigen::MatrixXd& C_virt,
      const std::vector<Eigen::MatrixXd>& U);

  Eigen::MatrixXd OrthogonalizeGuess(const Eigen::MatrixXd& GuessMOs) const;
  void PrintMOs(const Eigen::VectorXd& MOEnergies, Log::Level level);
  void PrintMOsUKS(const Eigen::VectorXd& alpha_energies,
                   const Eigen::VectorXd& beta_energies,
                   Log::Level level) const;
  void CalcElDipole(const Orbitals& orb) const;

  std::array<Eigen::MatrixXd, 2> CalcERIs_EXX(const Eigen::MatrixXd& MOCoeff,
                                              const Eigen::MatrixXd& Dmat,
                                              double error) const;

  Eigen::MatrixXd CalcERIs(const Eigen::MatrixXd& Dmat, double error) const;

  void ComputePolarizability(Orbitals& orb);

  bool EvaluateAndTime(Orbitals& orb);
  void ReportDimensionsAndMemory() const;

  void ConfigOrbfile(Orbitals& orb);
  void SetupInvariantMatrices();
  Eigen::MatrixXd McWeenyPurification(Eigen::MatrixXd& Dmat_in,
                                      AOOverlap& overlap);

  Mat_p_Energy SetupH0(const QMMolecule& mol) const;
  Mat_p_Energy IntegrateExternalMultipoles(
      const QMMolecule& mol,
      const std::vector<std::unique_ptr<StaticSite> >& multipoles) const;
  Mat_p_Energy IntegrateExternalDensity(const QMMolecule& mol,
                                        const Orbitals& extdensity) const;
  Eigen::MatrixXd IntegrateExternalField(const QMMolecule& mol) const;

  tools::EigenSystem IndependentElectronGuess(const Mat_p_Energy& H0) const;
  tools::EigenSystem ModelPotentialGuess(
      const Mat_p_Energy& H0, const QMMolecule& mol,
      const Vxc_Potential<Vxc_Grid>& vxcpotential) const;

  Eigen::MatrixXd AtomicGuess(const QMMolecule& mol) const;

  Eigen::VectorXd BuildEHTOrbitalEnergies(const QMMolecule& mol) const;
  Eigen::MatrixXd BuildEHTHamiltonian(const QMMolecule& mol) const;
  tools::EigenSystem ExtendedHuckelGuess(const QMMolecule& mol) const;
  tools::EigenSystem ExtendedHuckelDFTGuess(
      const Mat_p_Energy& H0, const QMMolecule& mol,
      const Vxc_Potential<Vxc_Grid>& vxcpotential) const;

  Orbitals BuildDimerGuessFromMonomerFiles(const QMMolecule& dimer_mol) const;

  Eigen::MatrixXd RunAtomicDFT_unrestricted(
      const QMAtom& uniqueAtom, bool use_hunds_rule_occupation = false) const;

  std::map<std::string, Eigen::MatrixXd> ComputeHirshfeldReferenceDensities(
      const QMMolecule& mol) const;

  struct CDFTConstraintSpec {
    std::vector<Index> atom_indices;
    // Charge relative to the neutral fragment (sum of its nuclear charges);
    // +1.0 means one electron removed. BuildCDFTConstraint converts it to an
    // absolute electron count (CP2K's TARGET convention).
    double target_charge = 0.0;
    double initial_lambda = 0.0;
    // "warmstart" (default) or "fresh", see dftpackage.xml. Kept as a string
    // so that invalid values are rejected by the XML choices validation.
    std::string guess_strategy = "warmstart";
  };

  HirshfeldPartition::Constraint BuildCDFTConstraint(
      const QMMolecule& mol, const CDFTConstraintSpec& spec) const;

  double NuclearRepulsion(const QMMolecule& mol) const;
  double ExternalRepulsion(
      const QMMolecule& mol,
      const std::vector<std::unique_ptr<StaticSite> >& multipoles) const;
  Eigen::MatrixXd SphericalAverageShells(const Eigen::MatrixXd& dmat,
                                         const AOBasis& dftbasis) const;

  void TruncateBasis(Orbitals& orb, std::vector<Index>& activeatoms,
                     Mat_p_Energy& H0,
                     Eigen::MatrixXd InitialActiveDensityMatrix,
                     Eigen::MatrixXd v_embedding,
                     Eigen::MatrixXd InitialInactiveMOs);

  void TruncMOsFullBasis(Orbitals& orb, std::vector<Index> activeatoms,
                         std::vector<Index> numfuncpatom);
  Eigen::MatrixXd InsertZeroCols(Eigen::MatrixXd MOsMatrix, Index startidx,
                                 Index numofzerocols);
  Eigen::MatrixXd InsertZeroRows(Eigen::MatrixXd MOsMatrix, Index startidx,
                                 Index numofzerorows);

  bool EvaluateClosedShell(Orbitals& orb, const Mat_p_Energy& H0,
                           const Vxc_Potential<Vxc_Grid>& vxcpotential);

  bool EvaluateUKS(Orbitals& orb, const Mat_p_Energy& H0,
                   const Vxc_Potential<Vxc_Grid>& vxcpotential);

  void ComputeAndStoreForces(Orbitals& orb, const Eigen::MatrixXd& Dmat,
                             const Vxc_Potential<Vxc_Grid>& vxcpotential) const;

  Eigen::MatrixXd ComputeOverlapPulayGradientUKS(
      const QMMolecule& mol, const tools::EigenSystem& MOs_alpha,
      const tools::EigenSystem& MOs_beta) const;

  Eigen::MatrixXd ComputeNonXCGradientUKS(
      const QMMolecule& mol, const UKSConvergenceAcc::SpinDensity& Dspin,
      const tools::EigenSystem& MOs_alpha,
      const tools::EigenSystem& MOs_beta) const;

  void ComputeAndStoreForcesUKS(
      Orbitals& orb, const UKSConvergenceAcc::SpinDensity& Dspin,
      const tools::EigenSystem& MOs_alpha, const tools::EigenSystem& MOs_beta,
      const Vxc_Potential<Vxc_Grid>& vxcpotential) const;
 
  Logger* pLog_;
 
  // basis sets
  std::string auxbasis_name_;
  std::string dftbasis_name_;
  std::string ecp_name_;
  AOBasis dftbasis_;
  AOBasis auxbasis_;
  ECPAOBasis ecp_;
 
  Index fock_matrix_reset_;
  // Pre-screening
  double screening_eps_;
  // Eigenvalues of the AO overlap below this are removed from S^-1/2
  // (xtpdft.overlap_tolerance).
  double overlap_tolerance_ = 1e-8;
  // Shell pairs whose RI three-centre integrals are bounded by less than this
  // are not stored (xtpdft.ri_pair_threshold).
  double ri_pair_threshold_ = 1e-10;
 
  // numerical integration Vxc
  std::string grid_name_;
 
  // AO Matrices
  AOOverlap dftAOoverlap_;
 
  std::string initial_guess_;
 
  // Only used if initial_guess_ == "dimer_guess", see
  // BuildDimerGuessFromMonomerFiles.
  std::string dimer_guess_orbA_name_;
  std::string dimer_guess_orbB_name_;
 
  // Convergence
  Index numofelectrons_ = 0;
  Index max_iter_;
  ConvergenceAcc::options conv_opt_;
  // DIIS variables
  ConvergenceAcc conv_accelerator_;
  // Electron repulsion integrals
  ERIs ERIs_;
  // Wall-clock accounting, reported at the end of Evaluate()
  mutable DFTTimings timings_;
  bool warm_start_ = false;
  bool warm_started_ = false;
  DFTSetupCache* setup_cache_ = nullptr;
  std::string eris_key_;  // key of the RI integrals in ERIs_
  std::string RISetupKey() const;
  bool UsableAsWarmStart(const Orbitals& orb, const std::string& previous_basis,
                         std::string& reason) const;
 
  // external charges
  std::vector<std::unique_ptr<StaticSite> >* externalsites_ = nullptr;
 
  // exchange and correlation
  double ScaHFX_ = 0.0;
  std::string xc_functional_name_;
 
  bool integrate_ext_density_ = false;
  // integrate external density
  std::string orbfilename_;
  std::string gridquality_;
  std::string state_;
 
  QMMolecule activemol_ =
      QMMolecule("molecule made of atoms participating in Active region", 1);
 
  Eigen::Vector3d extfield_ = Eigen::Vector3d::Zero();
  bool integrate_ext_field_ = false;
  // xtpdft.polarizability: static polarizability after the SCF
  bool compute_polarizability_ = false;
  Eigen::MatrixXd ext_potential_matrix_;
 
  std::string active_atoms_as_string_;
  double active_threshold_;
  double levelshift_;
 
  // truncation
  Eigen::MatrixXd H0_trunc_;
  Eigen::MatrixXd InitialActiveDmat_trunc_;
  Eigen::MatrixXd v_embedding_trunc_;
  bool truncate_;
  Index active_electrons_;
  double Total_E_full_;
  double E_nuc_;
  double truncation_threshold_;
  std::vector<Index> active_and_border_atoms_;
  std::vector<Index> numfuncpatom_;
 
  // QMEwald
  Vxc_Grid external_ewaldgrid_;
  bool has_ewaldgrid_ = false;
  double ewald_nuclear_energy_ = 0.0;
  // Spin-DFT Extension
  Index num_alpha_electrons_ = 0;
  Index num_beta_electrons_ = 0;
  Index num_docc_ = 0;
  Index num_socc_alpha_ = 0;
  Index spin_ = 1;
  Index charge_ = 0;
  bool force_uks_path_ = false;
  // xtpdft.compute_forces: compute the nuclear gradient after the SCF.
  // Off by default because of its cost (one-electron and overlap
  // derivatives, RI-J, XC and, for hybrids, RI-K gradients).
  bool compute_forces_ = false;
 
  // CDFT constraints added to the UKS Fock matrix; empty (no-op) for
  // ordinary runs. Set by RunCDFT, which updates each lambda between
  // warm-started EvaluateUKS calls.
  std::vector<HirshfeldPartition::Constraint> constraints_;
 
  // CDFT outer-loop control, used only by RunCDFT. Read from
  // cdft.max_iterations/cdft.population_tolerance in Initialize() if
  // cdft.enabled; the defaults apply when RunCDFT is called directly.
  Index max_cdft_iterations_ = 50;
  double cdft_population_tolerance_ = 1.e-4;
 
  bool cdft_enabled_ = false;
  CDFTConstraintSpec cdft_constraint_spec_;
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_DFTENGINE_H