chargeresponsefit.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_CHARGERESPONSEFIT_H
#define VOTCA_XTP_CHARGERESPONSEFIT_H
 
// Standard includes
#include <array>
#include <string>
#include <vector>
 
// Local VOTCA includes
#include "aobasis.h"
#include "eigen.h"
#include "hirshfeldpartition.h"
#include "vxc_grid.h"
 
namespace votca {
namespace tools {
class Property;
}
namespace xtp {
class Logger;
class Orbitals;

class ChargeResponseFit {
 public:
  struct Result {
    Eigen::MatrixXd pi;         // symmetrized, charge conserving
    double pi_asymmetry = 0.0;  // ||Pi - Pi^T|| / ||Pi|| before
    Eigen::VectorXd pi_eigenvalues;
    Eigen::Matrix3d alpha;                         // molecular, a.u.
    Eigen::Matrix3d alpha_ct;                      // X^T Pi X
    Eigen::Matrix3d alpha_res;                     // alpha - alpha_ct
    Eigen::VectorXd res_weights;                   // sum 1
    std::vector<Eigen::Matrix3d> alpha_res_atoms;  // weights * alpha_res
    std::vector<Eigen::Matrix3d> alpha_atoms;      // Hirshfeld, a.u.
  };
 
  struct FitOptions {
    bool analytic = true;
    double tolerance = 1e-7;
    double step = 1e-3;
    bool use_inversion = true;
  };

  static Result Fit(const tools::Property& dftoptions, Orbitals& orb,
                    const FitOptions& fit, Logger& log);

  struct Inversion {
    std::vector<Index> atom_map;
    Eigen::MatrixXd T;
  };
  static Inversion FindInversion(const QMMolecule& mol, const AOBasis& basis,
                                 double tol = 1e-4);

  static std::vector<Eigen::MatrixXd> WeightMatrices(
      const std::vector<HirshfeldPartition::AtomicReference>& atoms,
      const AOBasis& basis, const Vxc_Grid& grid);

  static std::vector<Eigen::Matrix3d> AtomicPolarizabilities(
      const std::vector<HirshfeldPartition::AtomicReference>& atoms,
      const std::vector<Eigen::Vector3d>& positions, const AOBasis& basis,
      const Vxc_Grid& grid, const std::array<Eigen::MatrixXd, 3>& dP_field);

  static Result Analyse(const std::vector<Eigen::Vector3d>& positions,
                        const std::vector<Eigen::MatrixXd>& W,
                        const std::vector<Eigen::MatrixXd>& dP_potential,
                        const std::array<Eigen::MatrixXd, 3>& dipole,
                        const std::array<Eigen::MatrixXd, 3>& dP_field,
                        const std::vector<Eigen::Matrix3d>& alpha_atoms);
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_CHARGERESPONSEFIT_H