stochasticnetworkfit.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_STOCHASTICNETWORKFIT_H
#define VOTCA_XTP_STOCHASTICNETWORKFIT_H
 
// Standard includes
#include <cstdint>
#include <string>
#include <vector>
 
// VOTCA includes
#include <votca/tools/constants.h>
 
// Local VOTCA includes
#include "eigen.h"
#include "qmstate.h"
#include "stochasticnetwork.h"
#include "topology.h"
 
namespace votca {
namespace xtp {

class StochasticNetworkFit {
 public:
  struct TrapInput {
    Index carrier = 0;
    std::vector<Index> sites;   // deepest first
    double release_rate = 0.0;  // 1/s
  };
 
  struct Settings {
    std::vector<QMStateType> carriers;
    std::vector<TrapInput> traps;
    double bin = 0.05 * tools::conv::nm2bohr;
    double correlation_cutoff = 0.0;
    std::vector<Index> l_b;  // grids
    std::vector<Index> l_c;
    Index quantiles = 100;  // 0: Gaussian residuals
    bool fit_points = true;
    StochasticNetwork::PointProcess points;
    Index realizations = 4;
    Index max_evaluations = 150;
    std::uint64_t seed = 1;
    double candidates = 3.0;
    Index max_iterations = 200;
  };

  struct Data {
    Eigen::Vector3d box = Eigen::Vector3d::Zero();
    std::vector<Eigen::Vector3d> positions;  // wrapped into the box
    std::vector<QMStateType> carriers;
    std::vector<std::vector<double>> energies;  // per carrier, Hartree
    std::vector<double> U_xX_nN;                // per carrier, means
    std::vector<double> U_nX_nN;
    std::vector<double> U_xN_xX;
    struct Pair {
      Index i;
      Index j;
      double dist;
      std::vector<double> J2;  // per carrier, Hartree^2
    };
    std::vector<Pair> pairs;
  };
  static Data FromTopology(const Topology& top,
                           const std::vector<QMStateType>& carriers);
  static Data FromNetwork(const StochasticNetwork& net,
                          const std::vector<QMStateType>& carriers,
                          const Eigen::Vector3d& box);

  struct Curve {
    std::vector<double> r;  // bohr
    std::vector<double> value;
    std::vector<double> count;  // weights (pairs per bin), may be empty
  };
 
  struct EnergyFit {
    double mean = 0.0;
    double sigma = 0.0;        // sample standard deviation
    double sigma_model = 0.0;  // sigma / sqrt(1 - cbar) (see kappa)
    double cbar = 0.0;         // mean correlation over all pairs
    Index l_b = 1;
    Index l_c = 1;
    double w_a = 1.0;
    double w_b = 0.0;
    double w_c = 0.0;
    double rms = 0.0;  // rms deviation of the kappa fit
    Curve kappa;       // microscopic (sample mean and variance)
    Curve kappa_model;
    bool fitted = false;  // false: sigma zero, no correlation fit
    Curve l_c_scan;
    bool l_b_at_edge = false;  // the best l_b / l_c is the end of its grid
    bool l_c_at_edge = false;
  };
 
  struct CouplingFit {
    StochasticNetwork::CouplingModel model;
    Index used = 0;
    Index excluded = 0;  // J^2 <= 0
    double skewness = 0.0;
    double excess_kurtosis = 0.0;
  };
 
  struct PointFit {
    StochasticNetwork::PointProcess process;
    double objective = 0.0;
    Index evaluations = 0;
    Curve D, H, g;                    // microscopic
    Curve D_model, H_model, g_model;  // model, mean over the realizations
  };
 
  struct Result {
    StochasticNetwork::Options options;
    PointFit points;
    std::vector<EnergyFit> energies;
    std::vector<CouplingFit> couplings;
    Curve connection;
  };
 
  StochasticNetworkFit(const Data& data, const Settings& settings);
  Result Fit() const;
 
  // statistics, also used for the validation (all distances bohr)
  static std::vector<double> NearestDistances(
      const std::vector<Eigen::Vector3d>& points, const Eigen::Vector3d& box);
  static std::vector<double> ContactDistances(
      const std::vector<Eigen::Vector3d>& points, const Eigen::Vector3d& box,
      Index probes, std::uint64_t seed);
  static Curve Cdf(const std::vector<double>& x, const std::vector<double>& r);
  static Curve PairCorrelation(const std::vector<Eigen::Vector3d>& points,
                               const Eigen::Vector3d& box, double bin,
                               double rmax);
  static Curve Correlation(const std::vector<Eigen::Vector3d>& points,
                           const Eigen::Vector3d& box,
                           const std::vector<double>& z, double bin,
                           double rmax, const std::vector<char>& skip = {});
  static double NormalQuantile(double u);
  static std::vector<Curve> Overlaps(const std::vector<Eigen::Vector3d>& points,
                                     const Eigen::Vector3d& box,
                                     const std::vector<Index>& l, double bin,
                                     double rmax);
  static double FitWeights(const Curve& kappa, const Curve& Ob, const Curve& Oc,
                           double& w_b, double& w_c);
  static Curve PairsPerSite(const Data& data, double bin, double rmax);
  static Curve MeanLog10J2(const Data& data, Index carrier, double bin,
                           double rmax);
 
 private:
  const Data& data_;
  Settings settings_;
 
  PointFit FitPoints() const;
  Curve FitConnection() const;
  CouplingFit FitCoupling(Index carrier,
                          std::vector<double>& standardized) const;
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_STOCHASTICNETWORKFIT_H