stochasticnetwork.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_STOCHASTICNETWORK_H
#define VOTCA_XTP_STOCHASTICNETWORK_H
 
// Standard includes
#include <array>
#include <cstdint>
#include <functional>
#include <random>
#include <string>
#include <vector>
 
// Local VOTCA includes
#include "eigen.h"
#include "qmstate.h"
#include "topology.h"
 
namespace votca {
namespace tools {
class Property;
}
}  // namespace votca
 
namespace votca {
namespace xtp {

class ModelRandom {
 public:
  ModelRandom(std::uint64_t seed, std::uint64_t stream);
  double Uniform() {
    return (double(engine_() >> 11) + 1.0) * (1.0 / 9007199254740992.0);
  }
  double Normal();
  Index Poisson(double mean);
  Index Integer(Index n);
 
 private:
  std::mt19937_64 engine_;
  bool has_spare_ = false;
  double spare_ = 0.0;
};

class PeriodicCells {
 public:
  PeriodicCells(const Eigen::Vector3d& box, double radius,
                const std::vector<Eigen::Vector3d>& points);
  Eigen::Vector3d Connection(const Eigen::Vector3d& r1,
                             const Eigen::Vector3d& r2) const;
  void Within(const Eigen::Vector3d& p, double radius,
              const std::function<void(Index, const Eigen::Vector3d&, double)>&
                  f) const;
  std::vector<Index> Nearest(Index i, Index l) const;
  std::vector<Index> NearestTo(const Eigen::Vector3d& p, Index l) const;
  Index NumberOfPoints() const { return Index(points_.size()); }
 
 private:
  Eigen::Vector3d box_;
  std::array<Index, 3> n_;
  Eigen::Vector3d cell_;
  const std::vector<Eigen::Vector3d>& points_;
  std::vector<Index> start_;
  std::vector<Index> members_;
  Index CellIndex(const Eigen::Vector3d& p, std::array<Index, 3>& c) const;
  template <class F>
  void Shell(const std::array<Index, 3>& c, Index s, F&& f) const;
};

class StochasticNetwork {
 public:
  struct PointProcess {
    double hard_core = 0.0;        // r_h
    double radius_mean = 0.0;      // mean of X
    double radius_variance = 0.0;  // variance of X (0: X = mean)
    double candidates = 3.0;
    Index max_iterations = 200;
  };
 
  struct EnergyModel {
    double mean = 0.0;
    double sigma = 0.0;  // 0: all sites at mean
    double w_a = 1.0;
    double w_b = 0.0;
    Index l_b = 1;
    Index l_c = 1;
  };
 
  struct CouplingModel {
    // log10(J^2/eV^2) mean = mean_offset + mean_slope r[nm]
    double mean_offset = -6.0;
    double mean_slope = 0.0;
    // variance = var_offset + var_slope r[nm] (clipped at 0)
    double var_offset = 1.0;
    double var_slope = 0.0;
    std::vector<double> residual_quantiles;
  };
 
  struct Carrier {
    QMStateType state;
    EnergyModel energy;
    CouplingModel coupling;
    double U_xX_nN = 0.0;  // internal energy
    double U_nX_nN = 0.0;  // reorganization parts
    double U_xN_xX = 0.0;
  };

  struct Trap {
    Index carrier = 0;     // the carrier it traps (index into carriers)
    double density = 0.0;  // traps of this kind per site
    std::vector<std::vector<double>> offsets;
    double release_rate = 0.0;  // microscopic release rate (1/s, information)
  };
 
  struct Options {
    Eigen::Vector3d box = Eigen::Vector3d::Zero();  // orthorhombic, bohr
    Index sites = 0;
    PointProcess points;
    std::vector<double> connection_r;
    std::vector<double> connection_p;
    std::vector<Carrier> carriers;
    double energy_correlation = 0.0;
    double coupling_correlation = 0.0;
    std::uint64_t seed = 1;
    std::string segment_name = "site";
    std::vector<Trap> traps;
    double trap_factor = 1.0;
  };
 
  explicit StochasticNetwork(const Options& opt);
 
  Index NumberOfSites() const { return Index(positions_.size()); }
  const std::vector<Eigen::Vector3d>& Positions() const { return positions_; }
  const std::vector<double>& Radii() const { return radii_; }
  struct InsertedTrap {
    Index type;
    std::vector<Index> sites;
  };
  const std::vector<InsertedTrap>& Traps() const { return inserted_; }
  const std::vector<double>& SiteEnergies(Index c) const {
    return energies_[c];
  }
  struct Pair {
    Index i;
    Index j;
    Eigen::Vector3d R;       // shortest connection from i to j
    std::vector<double> J2;  // per carrier, Hartree^2
  };
  const std::vector<Pair>& Pairs() const { return pairs_; }
  Index Iterations() const { return iterations_; }

  void BuildTopology(Topology& top) const;

  static Options FromProperty(const tools::Property& p);
  static std::string ToXML(const Options& opt);

  static std::vector<double> MovingAverage(const PeriodicCells& cells,
                                           const EnergyModel& m, Index sites,
                                           const std::vector<double>& A,
                                           const std::vector<double>& Bf,
                                           const std::vector<double>& C);
  static Index SamplePoints(const PointProcess& pp, const Eigen::Vector3d& box,
                            Index sites, ModelRandom& random,
                            std::vector<Eigen::Vector3d>& positions,
                            std::vector<double>& radii);
  double Connection(double r) const;
  static double Residual(const CouplingModel& m, double g);
  static double Log10J2(const CouplingModel& m, double r, double g);
 
 private:
  Options opt_;
  std::vector<Eigen::Vector3d> positions_;
  std::vector<double> radii_;
  std::vector<std::vector<double>> energies_;
  std::vector<Pair> pairs_;
  Index iterations_ = 0;
  std::vector<InsertedTrap> inserted_;
 
  void Check() const;
  void GenerateEnergies(ModelRandom& random);
  void GeneratePairs(ModelRandom& random);
  void InsertTraps(ModelRandom& random);
};
 
}  // namespace xtp
}  // namespace votca
 
#endif  // VOTCA_XTP_STOCHASTICNETWORK_H