convergenceacc.cc

Go to the documentation of this file.

/*
 *            Copyright 2009-2026 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.
 *
 */
 
// Local VOTCA includes
#include "votca/xtp/convergenceacc.h"
 
namespace votca {
namespace xtp {

 
// Build the symmetric orthogonalization matrix X = S^{-1/2}. All Fock-like
// matrices are diagonalized in the orthogonal AO basis X^T F X.
void ConvergenceAcc::setOverlap(AOOverlap& S, double etol) {
  S_ = &S;
  Sminusahalf = S.Pseudo_InvSqrt(etol);
  XTP_LOG(Log::error, *log_)
      << TimeStamp() << " Smallest value of AOOverlap matrix is "
      << S_->SmallestEigenValue() << std::flush;
  XTP_LOG(Log::error, *log_)
      << TimeStamp() << " Removed " << S_->Removedfunctions()
      << " basisfunction from inverse overlap matrix" << std::flush;
  return;
}
 
// Perform one SCF acceleration step.
//
// The commutator residual
//
//   R = X^T (F P S - S P F) X
//
// vanishes at self-consistency in a non-orthogonal AO basis. Its maximum
// element is used as the DIIS error metric, while the history of Fock and
// density matrices is passed to either DIIS or ADIIS to construct the next
// extrapolated Fock matrix. When the extrapolation is deemed unsafe, linear
// density mixing is used instead.
Eigen::MatrixXd ConvergenceAcc::Iterate(const Eigen::MatrixXd& dmat,
                                        Eigen::MatrixXd& H,
                                        tools::EigenSystem& MOs, double totE) {
  Eigen::MatrixXd H_guess = Eigen::MatrixXd::Zero(H.rows(), H.cols());
 
  if (int(mathist_.size()) == opt_.histlength) {
    totE_.erase(totE_.begin() + maxerrorindex_);
    mathist_.erase(mathist_.begin() + maxerrorindex_);
    dmatHist_.erase(dmatHist_.begin() + maxerrorindex_);
  }
 
  totE_.push_back(totE);
  if (opt_.mode != KSmode::fractional && nocclevels_ > 0 &&
      nocclevels_ < MOs.eigenvalues().size()) {
    double gap =
        MOs.eigenvalues()(nocclevels_) - MOs.eigenvalues()(nocclevels_ - 1);
    if ((diiserror_ > opt_.levelshiftend && opt_.levelshift > 0.0) ||
        gap < 1e-6) {
      Levelshift(H, MOs.eigenvectors());
    }
  }
  const Eigen::MatrixXd& S = S_->Matrix();
  Eigen::MatrixXd errormatrix =
      Sminusahalf.transpose() * (H * dmat * S - S * dmat * H) * Sminusahalf;
  diiserror_ = errormatrix.cwiseAbs().maxCoeff();
 
  mathist_.push_back(H);
  dmatHist_.push_back(dmat);
 
  if (opt_.maxout) {
    if (diiserror_ > maxerror_) {
      maxerror_ = diiserror_;
      maxerrorindex_ = mathist_.size() - 1;
    }
  }
 
  diis_.Update(maxerrorindex_, errormatrix);
  bool diis_error = false;
  XTP_LOG(Log::error, *log_)
      << TimeStamp() << " DIIs error " << getDIIsError() << std::flush;
 
  XTP_LOG(Log::error, *log_)
      << TimeStamp() << " Delta Etot " << getDeltaE() << std::flush;
 
  if ((diiserror_ < opt_.adiis_start || diiserror_ < opt_.diis_start) &&
      opt_.usediis && mathist_.size() > 2) {
    Eigen::VectorXd coeffs;
    // use ADIIs if energy has risen a lot in current iteration
 
    if (diiserror_ > opt_.diis_start ||
        totE_.back() > 0.9 * totE_[totE_.size() - 2]) {
      coeffs = adiis_.CalcCoeff(dmatHist_, mathist_);
      diis_error = !adiis_.Info();
      XTP_LOG(Log::warning, *log_)
          << TimeStamp() << " Using ADIIS for next guess" << std::flush;
 
    } else {
      coeffs = diis_.CalcCoeff();
      diis_error = !diis_.Info();
      XTP_LOG(Log::warning, *log_)
          << TimeStamp() << " Using DIIS for next guess" << std::flush;
    }
    if (diis_error) {
      XTP_LOG(Log::warning, *log_)
          << TimeStamp() << " (A)DIIS failed using mixing instead"
          << std::flush;
      H_guess = H;
    } else {
      for (Index i = 0; i < coeffs.size(); i++) {
        if (std::abs(coeffs(i)) < 1e-8) {
          continue;
        }
        H_guess += coeffs(i) * mathist_[i];
      }
    }
 
  } else {
    H_guess = H;
  }
 
  MOs = SolveFockmatrix(H_guess);
 
  Eigen::MatrixXd dmatout = DensityMatrix(MOs);
  if (diiserror_ > opt_.adiis_start || !opt_.usediis || diis_error ||
      mathist_.size() <= 2) {
    usedmixing_ = true;
    dmatout =
        opt_.mixingparameter * dmat + (1.0 - opt_.mixingparameter) * dmatout;
    XTP_LOG(Log::warning, *log_)
        << TimeStamp() << " Using Mixing with alpha=" << opt_.mixingparameter
        << std::flush;
  } else {
    usedmixing_ = false;
  }
  return dmatout;
}
 
void ConvergenceAcc::PrintConfigOptions() const {
  XTP_LOG(Log::error, *log_)
      << TimeStamp() << " Convergence Options:" << std::flush;
  XTP_LOG(Log::error, *log_)
      << "\t\t Delta E [Ha]: " << opt_.Econverged << std::flush;
  XTP_LOG(Log::error, *log_)
      << "\t\t DIIS max error: " << opt_.error_converged << std::flush;
  if (opt_.usediis) {
    XTP_LOG(Log::error, *log_)
        << "\t\t DIIS histlength: " << opt_.histlength << std::flush;
    XTP_LOG(Log::error, *log_)
        << "\t\t ADIIS start: " << opt_.adiis_start << std::flush;
    XTP_LOG(Log::error, *log_)
        << "\t\t DIIS start: " << opt_.diis_start << std::flush;
    std::string del = "oldest";
    if (opt_.maxout) {
      del = "largest";
    }
    XTP_LOG(Log::error, *log_)
        << "\t\t Deleting " << del << " element from DIIS hist" << std::flush;
  }
  XTP_LOG(Log::error, *log_)
      << "\t\t Levelshift[Ha]: " << opt_.levelshift << std::flush;
  XTP_LOG(Log::error, *log_)
      << "\t\t Levelshift end: " << opt_.levelshiftend << std::flush;
  XTP_LOG(Log::error, *log_)
      << "\t\t Mixing Parameter alpha: " << opt_.mixingparameter << std::flush;
}
 
// Solve the generalized Roothaan-Hall problem
//
//   F C = S C eps
//
// by symmetric orthogonalization: diagonalize X^T F X with X = S^{-1/2} and
// back-transform the eigenvectors as C = X C' .
tools::EigenSystem ConvergenceAcc::SolveFockmatrix(
    const Eigen::MatrixXd& H) const {
  // transform to orthogonal for
  Eigen::MatrixXd H_ortho = Sminusahalf.transpose() * H * Sminusahalf;
  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(H_ortho);
 
  if (es.info() != Eigen::ComputationInfo::Success) {
    throw std::runtime_error("Matrix Diagonalisation failed. DiagInfo" +
                             std::to_string(es.info()));
  }
 
  tools::EigenSystem result;
  result.eigenvalues() = es.eigenvalues();
 
  result.eigenvectors() = Sminusahalf * es.eigenvectors();
  return result;
}
 
// Add a virtual-space level shift
//
//   F <- F + S C_virt Delta C_virt^T S,
//
// implemented by placing the scalar shift on the virtual diagonal in the MO
// basis built from the previous iteration. This leaves occupied orbitals
// unchanged while opening the HOMO-LUMO gap during difficult SCF phases.
void ConvergenceAcc::Levelshift(Eigen::MatrixXd& H,
                                const Eigen::MatrixXd& MOs_old) const {
  if (opt_.levelshift < 1e-9) {
    return;
  }
  Eigen::VectorXd virt = Eigen::VectorXd::Zero(H.rows());
  for (Index i = nocclevels_; i < H.rows(); i++) {
    virt(i) = opt_.levelshift;
  }
 
  XTP_LOG(Log::error, *log_)
      << TimeStamp() << " Using levelshift:" << opt_.levelshift << " Hartree"
      << std::flush;
  Eigen::MatrixXd vir = S_->Matrix() * MOs_old * virt.asDiagonal() *
                        MOs_old.transpose() * S_->Matrix();
  H += vir;
  return;
}
 
/*
Eigen::MatrixXd ConvergenceAcc::DensityMatrix(
    const tools::EigenSystem& MOs) const {
  Eigen::MatrixXd result;
  if (opt_.mode == KSmode::closed) {
    result = DensityMatrixGroundState(MOs.eigenvectors());
  } else if (opt_.mode == KSmode::open) {
    result = DensityMatrixGroundState_unres(MOs.eigenvectors());
  } else if (opt_.mode == KSmode::fractional) {
    result = DensityMatrixGroundState_frac(MOs);
  }
  return result;
} */
 
// Closed-shell AO density matrix
//
//   P = 2 C_occ C_occ^T,
//
// where each occupied spatial orbital contributes one alpha and one beta
// electron.
Eigen::MatrixXd ConvergenceAcc::DensityMatrixGroundState(
    const Eigen::MatrixXd& MOs) const {
  const Eigen::MatrixXd occstates = MOs.leftCols(nocclevels_);
  Eigen::MatrixXd dmatGS = 2.0 * occstates * occstates.transpose();
  return dmatGS;
}
 
// Spin-resolved unrestricted AO density matrix for one spin channel,
//
//   P^sigma = C_occ^sigma (C_occ^sigma)^T,
//
// with no factor of two because a single spin channel is represented.
Eigen::MatrixXd ConvergenceAcc::DensityMatrixGroundState_unres(
    const Eigen::MatrixXd& MOs) const {
  if (nocclevels_ == 0) {
    return Eigen::MatrixXd::Zero(MOs.rows(), MOs.rows());
  }
  Eigen::MatrixXd occstates = MOs.leftCols(nocclevels_);
  Eigen::MatrixXd dmatGS = occstates * occstates.transpose();
  return dmatGS;
}
 
// Fractionally occupied AO density matrix
//
//   P = C f C^T,
//
// where f is a diagonal matrix of orbital occupations assembled from the
// configured electron count.
// Fractional-occupation AO density matrix
//
//   P = C n C^T,
//
// where n is the diagonal matrix of orbital occupations provided in the
// EigenSystem container.
Eigen::MatrixXd ConvergenceAcc::DensityMatrixGroundState_frac(
    const tools::EigenSystem& MOs) const {
  if (opt_.numberofelectrons == 0) {
    return Eigen::MatrixXd::Zero(MOs.eigenvectors().rows(),
                                 MOs.eigenvectors().rows());
  }
 
  Eigen::VectorXd occupation = Eigen::VectorXd::Zero(MOs.eigenvalues().size());
  std::vector<std::vector<Index> > degeneracies;
  double buffer = 1e-4;
  degeneracies.push_back(std::vector<Index>{0});
  for (Index i = 1; i < occupation.size(); i++) {
    if (MOs.eigenvalues()(i) <
        MOs.eigenvalues()(degeneracies[degeneracies.size() - 1][0]) + buffer) {
      degeneracies[degeneracies.size() - 1].push_back(i);
    } else {
      degeneracies.push_back(std::vector<Index>{i});
    }
  }
  Index numofelec = opt_.numberofelectrons;
  for (const std::vector<Index>& deglevel : degeneracies) {
    Index numofpossibleelectrons = 2 * Index(deglevel.size());
    if (numofpossibleelectrons <= numofelec) {
      for (Index i : deglevel) {
        occupation(i) = 2;
      }
      numofelec -= numofpossibleelectrons;
    } else {
      double occ = double(numofelec) / double(deglevel.size());
      for (Index i : deglevel) {
        occupation(i) = occ;
      }
      break;
    }
  }
  Eigen::MatrixXd dmatGS = MOs.eigenvectors() * occupation.asDiagonal() *
                           MOs.eigenvectors().transpose();
  return dmatGS;
}
 
/*******************************************************
 * EXTENSION FOR SPIN-KS-DFT
 *******************************************************/
 
ConvergenceAcc::SpinDensity
    ConvergenceAcc::DensityMatrixGroundState_restricted_open(
        const Eigen::MatrixXd& MOs) const {
 
  const Index n_docc =
      std::min(opt_.number_alpha_electrons, opt_.number_beta_electrons);
  const Index n_socc_alpha = opt_.number_alpha_electrons - n_docc;
 
  SpinDensity result;
  result.alpha = Eigen::MatrixXd::Zero(MOs.rows(), MOs.rows());
  result.beta = Eigen::MatrixXd::Zero(MOs.rows(), MOs.rows());
 
  if (n_docc > 0) {
    const Eigen::MatrixXd docc = MOs.leftCols(n_docc);
    const Eigen::MatrixXd d_docc = docc * docc.transpose();
    result.alpha += d_docc;
    result.beta += d_docc;
  }
 
  if (n_socc_alpha > 0) {
    const Eigen::MatrixXd socc = MOs.middleCols(n_docc, n_socc_alpha);
    result.alpha += socc * socc.transpose();
  }
 
  return result;
}
 
// Construct spin-resolved densities according to the configured occupation
// model. For restricted open-shell cases the same spatial orbitals are split
// into doubly and singly occupied subsets using the stored alpha/beta counts.
ConvergenceAcc::SpinDensity ConvergenceAcc::DensityMatrixSpinResolved(
    const tools::EigenSystem& MOs) const {
 
  if (opt_.mode == KSmode::restricted_open) {
    return DensityMatrixGroundState_restricted_open(MOs.eigenvectors());
  } else if (opt_.mode == KSmode::closed) {
    Eigen::MatrixXd d = DensityMatrixGroundState(MOs.eigenvectors());
    return {0.5 * d, 0.5 * d};
  } else if (opt_.mode == KSmode::open) {
    Eigen::MatrixXd d = DensityMatrixGroundState_unres(MOs.eigenvectors());
    return {d, Eigen::MatrixXd::Zero(d.rows(), d.cols())};
  } else {
    Eigen::MatrixXd d = DensityMatrixGroundState_frac(MOs);
    return {0.5 * d, 0.5 * d};
  }
}
 
Eigen::MatrixXd ConvergenceAcc::DensityMatrix(
    const tools::EigenSystem& MOs) const {
  SpinDensity spin_dmat = DensityMatrixSpinResolved(MOs);
  return spin_dmat.total();
}
 
}  // namespace xtp
}  // namespace votca