trah.h

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.
 *
 */
 
#pragma once
#ifndef VOTCA_XTP_TRAH_H
#define VOTCA_XTP_TRAH_H
 
// Standard includes
#include <functional>
 
// Local VOTCA includes
#include "eigen.h"
 
namespace votca {
namespace xtp {

namespace trah {
 
struct StepOptions {
  // The micro-iterations stop when |(H - mu) x + g| <= eta |g| with the
  // forcing term eta = min(max_forcing, |g|), so the macro-iterations
  // converge quadratically close to the minimum.
  double max_forcing = 0.1;
  double residual_floor = 1e-10;
  Index max_iterations = 50;
  Index max_subspace = 60;
  double alpha_max = 1e4;
  double lower_fraction = 0.8;
};
 
struct StepResult {
  Eigen::VectorXd x;
  // g.x + 1/2 x.H.x of the step taken, and g.x alone
  double predicted = 0.0;
  double linear = 0.0;
  double alpha = 1.0;
  double mu = 0.0;
  Index iterations = 0;
  Index products = 0;
  // the step started from the subspace passed in
  bool reused = false;
  double residual = 0.0;
  bool converged = false;
};

using HessianProduct = std::function<Eigen::MatrixXd(const Eigen::MatrixXd&)>;

struct Subspace {
  Eigen::MatrixXd B;
  Eigen::MatrixXd HB;
  Eigen::VectorXd g;
  void Clear() {
    B.resize(0, 0);
    HB.resize(0, 0);
    g.resize(0);
  }
};
 
StepResult Step(const Eigen::VectorXd& g, const Eigen::VectorXd& diag,
                const HessianProduct& hessian, double trust_radius,
                const StepOptions& opt = StepOptions(),
                Subspace* subspace = nullptr);

class TrustRadius {
 public:
  explicit TrustRadius(double radius = 0.5, double max_radius = 1.0)
      : radius_(radius), max_radius_(max_radius) {}
  double Radius() const { return radius_; }
  double LastRatio() const { return ratio_; }
  bool Update(double predicted, double actual);
  bool Update(double predicted, double actual, double linear, double step_norm);
 
  static constexpr double noise_floor = 1e-10;
 
 private:
  double radius_;
  double max_radius_;
  double ratio_ = 0.0;
};

Eigen::MatrixXd RotateOrbitals(const Eigen::MatrixXd& C,
                               const Eigen::VectorXd& x, Index nocc);

Eigen::VectorXd OccVirt(const Eigen::MatrixXd& M, Index nocc);
 
}  // namespace trah
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_TRAH_H