source: src/Potentials/PartialNucleiChargeFitter.hpp@ eaf08f

/*
 * PartialNucleiChargeFitter.hpp
 *
 *  Created on: 12.05.2013
 *      Author: heber
 */

#ifndef PARTIALNUCLEICHARGEFITTER_HPP_
#define PARTIALNUCLEICHARGEFITTER_HPP_

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <vector>

#include "LinearAlgebra/Vector.hpp"
#include "LinearAlgebra/VectorContent.hpp"

#include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"

class PartialNucleiChargeFitterTest;
class SamplingGrid;

/** This class is a functor that fits the magnitudes of a give number of
 * point charges to a sampled electrostatic potential distribution by
 * solving an over-determined system of equations.
 *
 * We require thus a sampled potential and the spatial positions of the
 * charges. As a result we offer the partial point charge at each position.
 *
 * Note that we fit to the potential of both nuclei and electrons.
 * And also note that open boundary conditions should be used as we compare
 * to the potential of spherical point charges and we do not do any Wolf
 * or Ewald summation.
 */
class PartialNucleiChargeFitter
{
  //!> grant unit test access to private parts.
  friend class PartialNucleiChargeFitterTest;
public:
  //!> typedef for the charge type to fit
  typedef double charge_t;
  //!> typedef for the number of charges to fit at the same time
  typedef std::vector<charge_t> charges_t;
  //!> typedef for specifying a position in 3D
  typedef Vector position_t;
  //!> typedef for specifying positions of all nuclei whose charges to fit
  typedef std::vector< position_t > positions_t;
  typedef size_t dimension_t;
  typedef std::vector<dimension_t> dimensions_t;

  /** Constructor for class PartialNucleiChargeFitter.
   *
   * \note Copies given parameters as it does not impact on overall
   * performance significantly.
   *
   * \warning the \a threshold parameter is important! Core electrons are very
   * strongly localized and can in the general case not be properly sampled on
   * a grid. As for the partial charges we are not interested in modelling the
   * potential around the nuclei but around the molecule as a whole, with this
   * parameter the inner potential is masked and excluded from the fit. This does
   * not mask any charge, it just avoids the strongly peaked electrons along with
   * "in short-range falsely" smeared-out nuclei charges. We recommend a value of
   * e.\,g. 1 or more. This should stabilize of finer grids and larger thresholds.
   *
   * \param grid grid with sampled (long-range) potential to match
   * \param _positions vector of positions
   * \param _threshold radius of spherical mask around nuclei
   */
  PartialNucleiChargeFitter(
      const SamplingGrid &grid,
      const positions_t &_positions,
      const double _threshold = 0.);

  /** Destructor for class PartialNucleiChargeFitter.
   *
   */
   ~PartialNucleiChargeFitter();

   /** Function to evaluate the over-determined system given by the
    * desired sampled potential and a number of charges at given positions.
    *
    * \return L_2-Fehler des Residuums (Ax-b)
    */
   double operator()();

   /** Getter for the problem matrix A.
    *
    * The matrix represents the sampled potential per charge/column vector
    * with unit charge.
    *
    * \return matrix A in system Ax=b or NLL if not constructed yet.
    */
   const MatrixContent &getMatrix() const
   {
     return *PotentialFromCharges;
   }

   /** Getter for solution as charges_t
    *
    * \return solution vector x as type charges_t
    */
   charges_t getSolutionAsCharges_t() const;

   /** Helper function to construct the discretized Coulomb operator for the
    * nuclei charges.
    *
    */
   void writeMatrix();

   /** Setter for solution as charges_t
    *
    * \param _charges solution vector x as type charges_t
    */
   void setCharges(const charges_t &_charges) const;

   /** Calculates the residuum for the currently set charges.
    *
    * \return residual for every charge
    */
   VectorContent calculateResiduum();

private:
   /** Helper function to construct the problem Matrix A.
    *
    */
   void constructMatrix();

   /** Internal function to calculate the discrete grid dimension per axis.
    *
    * \param grid grid whose dimensions are to be determind
    * \return vector of 3 grid dimensions
    */
   dimensions_t getGridDimensions(const SamplingGrid &grid) const;

   /** Helper function to mask out regions around nuclei.
    *
    * The electronic potential is too peaked around the nuclei, i.e.
    * cannot be sanely modelled by the partial point charges. Hence,
    * we mask out these regions
    *
    * \return true - point is at least threshold away from any nuclei, false - else
    */
   bool isGridPointSettable(
       const positions_t &_positions,
       const Vector &grid_position) const;

private:
   //!> grid dimensions per axis for internal use
   const dimensions_t total;
   //!> sampled potential as right hand side for Ax=b
   const VectorContent SampledPotential;
   //!> properties of grid for calculating potential on same grid
   const SamplingGridProperties grid_properties;
   //!> positions of nuclei whose charges are to be fitted
   const positions_t positions;
   //!> internal matrix representing potential by spherical charges (we use MatrixContent as there is only a 3x3 RealSpaceMatrix, no generic one)
   MatrixContent *PotentialFromCharges;
   //!> internal representation of solution vector x
   VectorContent *PartialCharges;
   //!> threshold for minimum distance to any nuclei in fitting
   const double threshold;
};

#endif /* PARTIALNUCLEICHARGEFITTER_HPP_ */