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_; };
  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);

  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;

  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;
    // Relative to the fragment's own neutral reference state (the sum
    // of its atoms' nuclear charges) -- e.g. +1.0 means one electron
    // REMOVED (a cation). Converted to an absolute target electron
    // count once, inside BuildCDFTConstraint, matching CP2K's own
    // internal (absolute) TARGET convention exactly -- only the
    // user-facing options syntax is charge-relative, per the earlier
    // design discussion on this.
    double target_charge = 0.0;
    double initial_lambda = 0.0;
    // "warmstart" (default) or "fresh" -- see this field's own XML
    // help text (dftpackage.xml) for the full reasoning. Stored as
    // the raw string, not a bool, so an invalid value (a typo, say)
    // is caught by the XML schema's own choices="..." validation
    // rather than silently defaulting to one behavior or the other.
    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 read from options / actually used when initial_guess_ ==
  // "dimer_guess" -- see BuildDimerGuessFromMonomerFiles's own header
  // comment for what this guess does and why it exists.
  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;
 
  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;
  // Default false to preserve existing performance for callers not
  // using this feature -- computing forces adds real, non-trivial cost
  // (kinetic/nuclear-attraction/overlap derivatives, RI-J gradient, full
  // XC gradient, RI-K for hybrids) to every converged SCF, so this must
  // be explicit opt-in, not silently always-on. Settable via the
  // <xtpdft> options block, which flows through unmodified from
  // XTPDFT::RunDFT() (options_) straight into DFTEngine::Initialize --
  // confirmed directly by reading XTPDFT::ParseSpecificOptions, which
  // only extracts a single unrelated field (temporary_file) and does
  // not filter/transform anything else -- so this option is
  // automatically available through the full QMPackage/XTPDFT flow with
  // no changes needed there.
  bool compute_forces_ = false;
 
  // Empty by default -- the ONLY thing a standard, non-CDFT run needs
  // to know about this member is that it is empty, checked via a
  // single, cheap constraints_.empty() guard inside
  // EvaluateUKS's own Fock-matrix assembly (see that function's own
  // comment at the point the constraint potential term is added).
  // When empty, that guard means the added term is a complete no-op:
  // the Hamiltonian is built exactly as it always was, with no
  // measurable overhead and no change in behavior whatsoever for any
  // run that never touches this member. Populated only by the
  // (not yet implemented) outer Lagrange-multiplier optimization loop,
  // which is expected to modify each Constraint's own lambda field in
  // place between successive, warm-started calls into EvaluateUKS --
  // per the design discussion this grew out of (CP2K's own documented
  // approach: restart the inner SCF from the previous trial's
  // converged density at each new lambda, rather than a cold start).
  std::vector<HirshfeldPartition::Constraint> constraints_;
 
  // CDFT outer-loop (Lagrange-multiplier) control, used only by
  // RunCDFT below -- never read by the ordinary Evaluate/EvaluateUKS
  // path at all, so these have no bearing on any standard run either.
  // max_cdft_iterations_/cdft_population_tolerance_ are also settable
  // from options (see Initialize()'s own cdft.max_iterations/
  // cdft.population_tolerance parsing) when cdft.enabled=true; their
  // defaults here are what a directly-constructed RunCDFT call (e.g.
  // from a test, bypassing Initialize() entirely) gets instead.
  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