source: src/documentation/constructs/potentials.dox@ c459c4

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/**
 * \file potentials.dox
 *
 * Created on: Nov 28, 2012
 *    Author: heber
 */

/** \page potentials Empirical Potentials and FunctionModels
 *
 * On this page we explain what is meant with the Potentials sub folder.
 *
 * First, we are based on fragmenting a molecular system, i.e. dissecting its
 * bond structure into connected subgraphs, calculating the energies of the
 * fragments (ab-initio) and summing up to a good approximation of the total
 * energy of the whole system, \sa fragmentation.
 * Second, having calculated these energies, there quickly comes up the idea
 * that one actually calculates quite similar systems all time and if one could
 * not cache results in an intelligent (i.e. interpolating) fashion ...
 *
 * That's where so-called empirical potentials come into play. They are
 * functions depending on a number of "fitted" parameters and the variable
 * distances within a molecular fragment (i.e. the bond lengths) in order to
 * give a value for the total energy without the need to solve a complex
 * ab-initio model (essentially, not solving the electronic Schrödinger equation
 * anymore). And they are accompanied by a specific binding model that
 * represents what kind of many-body force is represented by the potential.
 *
 * Empirical potentials have been thought of by fellows such as Lennard-Jones,
 * Morse, Tersoff, Stillinger and Weber, etc. And in their honor, most of the
 * potential forms are named after its inventor. Hence, we speak e.g. of a
 * Lennard-Jones potential.
 *
 * So, what we have to do in order to cache results is the following procedure:
 * -# gather similar fragments
 * -# perform a fit procedure to obtain the parameters for the empirical
 *    potential
 * -# evaluate the potential instead of an ab-initio calculation
 *
 * However, we need more: What are similar fragments? How do we perform the
 * fitting procedure? And as the potentials are mathematical functions, what
 * arguments do they depend on?
 *
 * Similar fragments are those that share the same bond graph, i.e. they have
 * the same number of nodes and the same number of edges. And each edge is
 * between the same two elements.
 *
 * The fitting procedure works by looking at a training set, i.e. a list of
 * elements where each contains an energy and a number of arguments, namely
 * pair-wise distances. The error is then the difference between the energies
 * in the set and all the energy values that we obtain when we feed the
 * arguments into the fitted potentials. This error is minimized in the
 * euclidian norm, i.e. least squares regression. But other norms might be
 * possible in the future, too.
 * 
 * And the pair-wise distances, we mentioned are the arguments.
 *
 * The terms, that we use, model the classes that are implemented:
 * -# EmpiricalPotential: Contains the interface to a function that can be
 *    evaluated given a number of arguments_t, i.e. distances. Also, one might
 *    want to evaluate derivatives.
 * -# FunctionModel: Is a function that can be fitted, i.e. it depends on a
 *    set of internal parameters that can be set and got.
 * -# argument_t: The Argument stores not only the distance but also the index
 *    pair of the associated atoms and also their charges, to let the potential
 *    check on validity.
 * -# SerializablePotential: Eventually, one wants to store to or parse from
 *    a file all potential parameters. This functionality is encoded in this
 *    class.
 * -# HomologyGraph: "Similar" fragments in our case have to have the same bond
 *    graph. It is stored in the HomologyGraph that acts as representative.
 * -# HomologyContainer: This container combines, in a ultimap fashion, all
 *    similar fragments with their energies together, with the HomologyGraph
 *    as their "key".
 * -# TrainingData: Here, we combine InputVector and OutputVector that contain
 *    the set of distances required for the FunctionModel (e.g. only a single
 *    distance/argument for a pair potential, three for an angle potential,
 *    etc.) and also the expected OutputVector, i.e. the energy of the specific
 *    configuration in our case. This in combination with the FunctionModel is
 *    the basis for the non-linear regression used for the fitting procedure.
 * -# Extractors: These set of functions yield the set of distances from a
 *    given fragment that is stored in the HomologyContainer.
 * -# FunctionApproximation: Contains the interface to the levmar package where
 *    the Levenberg-Marquardt (Newton + Trust region) algorithm is used to
 *    perform the fit.
 *
 * \section potentials-fit-potential-action What happens in FitPotentialAction.
 *
 *  First, charges and a potential type is used from the given options. This 
 *  is used to instantiate EmpiricalPotentials via the PotentialFactory, stored
 *  within the PotentialRegistry. This is the available set of potentials
 *  (without requiring any knowledge as to the nature of the fragment employed
 *  in fitting).
 *
 *  Second, the given fragment is used to get a suitable HomologyGraph from
 *  the World's HomologyContainer. This is given to a CompoundPotential, that in
 *  turn browses through the PotentialRegistry, picking out those
 *  EmpiricalPotential's that match with a subset of the FragmentNode's of the
 *  given graph. These are stored as a list of FunctionModel's within the
 *  CompoundPotential instance. Here comes the specific fragment into play,
 *  picking a subset from the available potentials.
 *
 *  Third, we need to setup the training data. For this we need vectors of input
 *  and output data that are obtained from the HomologyContainer with the
 *  HomologyGraph as key. The output vector in our case is simply a number
 *  (although the interface allows for more). The input vector is the set of
 *  distances. In order to pre-process the input data for the specific model
 *  a filter is required in the TrainingData's constructor. The purpose of the
 *  filter is to pick out the subset of distance arguments for each model one
 *  after the other and concatenate them to a list. On evaluation of the model
 *  this concatenated list of distances is given to the model and it may easily
 *  dissect the list and hand over each contained potential its subset of
 *  arguments. See Extractors for more information.
 *
 *  Afterwards, training may commence: The goal is to find a set of parameters
 *  for the model such that it as good as possible reproduces the output vector
 *  for a given input vector. This non-linear regression is contained in the
 *  levmar package and its functionality is wrapped in the FunctionApproximation
 *  class. An instance is initialized with both the gathered training data and
 *  the model containing a list of potentials. See
 *  [FunctionApproximation-details] for more details.
 *
 *  During the fitting procedure, first the derivatives of the model is checked
 *  for consistency, then the model is initialized with a sensible guess of
 *  starting parameters, and afterwards the Levenberg-Marquardt algorithm
 *  commences that makes numerous calls to evaluate the model and its derivative
 *  to find the minimum in the L2-norm.
 *
 *  This is done more than once as high-dimensional regression is sensitive to
 *  the starting values as there are possible numerous local minima. The lowest
 *  of the found minima is taken, either via a given threshold or the best of a
 *  given number of attempts.
 *
 *  Eventually, these parameters of the best model are streamed via
 *  PotentialSerializer back into a potential file. Each EmpiricalPotential in
 *  the CompoundPotential making up the whole model is also a
 *  SerializablePotential. Hence, each in turn writes a single line with its
 *  respective subset of parameters and particle types, describing this
 *  particular fit function.
 *
 * \section potentials-function-evaluation How does the model evaluation work
 *
 *  We now come to the question of how the model and its derivative are actually
 *  evaluated. We have an input vector from the training data and we have the
 *  model with a specific set of parameters.
 *
 *  FunctionModel is just an abstract interface that is implemented by the
 *  potential functions, such as CompoundPotential, that combines multiple
 *  potentials into a single function for fitting, or PairPotential_Harmonic,
 *  that is a specific fit function for harmonic bonds.
 *
 *  The main issue with the evaluation is picking the right set of distances from
 *  ones given in the input vector and feed it to each potential contained in
 *  CompoundPotential. 
 *
 *  Initially, the HomologyGraph only contains a list of configurations of a
 *  specific fragments (i.e. the position of each atom in the fragment) and an
 *  energy value. These first have to be converted into distances.
 *
 *  These distances are prepared by the TrainingData instantiation, i.e. a 
 *  fragment with all its atomic positions has already been converted to the
 *  set of all pair-wise interatomic distances.
 * 
 * \section potentials-distance-picking How does the distance picking work
 *
 *  Given a set of pair-wise distances, how do we pick the subset of distances
 *  needed by a particular potential.
 *  
 *  Let's make an example first: Imagine a water molecule, i.e. one oxygen and
 *  and two hydrogen atoms with two O-H bonds. Naturally, we obtain three pair-
 *  wise distances, OH1, OH2 and H1H2. Now, we would like to fit a Morse
 *  potential that just depends on a single interatomic distance. We would like
 *  it to represents the O-H bond energy. Hence, the distance picker, namely
 *  the Extractor function, needs to pick any subset of distance that contains
 *  a unique single O-H distance. In effect, it needs to return a list 
 *  containing OH1 and OH2 as the Morse potential needs to represent both bond
 *  energies together.
 *  
 *  Now, this is really still quite simple as long as the potential only 
 *  depends on a single distance. However, what if we continue and look at a
 *  angle potential, requiring three atoms, i.e. H-O-H?
 *  
 *  Or even more complicated: Imagine an ethane molecule (C2H6) and we would
 *  to represent the H-C-C angular interaction by a harmonic angle potential.
 *  Now, there are multiple of these at the same time, namely six angular
 *  interactions.
 *  
 *  What have to do is look for subgraphs inside a graph. Each potential comes
 *  with a small graph that represents the binding structure, in our terms 
 *  the bond model, that we expect. And we need to find the all matching
 *  subgraphs in the whole graph being the fragment itself. Then, for each
 *  subgraph the potential tells us in what order the pair-wise distances
 *  associated with the subgraph are required to be. All of these subset of
 *  distances are eventually concatenated and fed into the model on evaluation.
 *  
 * \section potentials-howto-use Howto use the potentials
 *
 *  We just give a brief run-down in terms of code on how to use the potentials.
 *  Here, we just describe what to do in order to perform the fitting. This is
 *  basically what is implemented in FragmentationFitPotentialAction.
 *
 *  \code
 *  // we need the homology container and the representative graph we want to
 *  // fit to.
 *  HomologyContainer homologies;
 *  const HomologyGraph graph = getSomeGraph(homologies);
 *  Fragment::atomicnumbers_t h2o;
 *  h2o += 8,1,1;
 *  // TrainingData needs so called Extractors to get the required distances
 *  // from the stored fragment. These functions are bound via boost::bind.
 *  TrainingData AngleData(
 *      boost::bind(&Extractors::gatherDistancesFromFragment,
 *          boost::bind(&Fragment::getPositions, _1),
 *          boost::bind(&Fragment::getAtomicNumbers, _1),
 *          boost::cref(h2o),
 *          _2)
 *      );
 *  // now we extract the distances and energies and store them
 *  AngleData(homologies.getHomologousGraphs(graph));
 *  // give ParticleTypes of this potential to make it unique
 *  PairPotential_Angle::ParticleTypes_t types =
 *        boost::assign::list_of<PairPotential_Angle::ParticleType_t>
 *        (8)(1)(1)
 *      ;
 *  PairPotential_Angle angle(types);
 *  // give initial parameter
 *  FunctionModel::parameters_t params(PairPotential_Angle::MAXPARAMS, 0.);
 *  // ... set some potential-specific initial parameters in params struct
 *  angle.setParameters(params);
 *
 *  // use the potential as a FunctionModel along with prepared TrainingData
 *  FunctionModel &model = angle;
 *  FunctionApproximation approximator(AngleData, model);
 *  approximator(FunctionApproximation::ParameterDerivative);
 *
 *  // obtain resulting parameters and check remaining L_2 and L_max error
 *  angleparams = model.getParameters();
 *  LOG(1, "INFO: L2sum = " << AngleData(model)
 *      << ", LMax = " << AngleData(model) << ".");
 *  \endcode
 *
 *  The evaluation of the fitted potential is then trivial, e.g.
 *  \code
 *  // constructed someplace
 *  PairPotential_Angle angle(...);
 *
 *  // evaluate
 *  FunctionModel::arguments_t args;
 *  // .. initialise args to the desired distances
 *  const double value = angle(args)[0]; // output is a vector!
 *  \endcode
 *
 * \section potentials-stability-of-fit note in stability of fit
 *
 *  As we always start from random initial parameters (within a certain sensible
 *  range at least), the non-linear fit does not always converge. Note that the
 *  random values are drawn from the defined distribution and the uniform distribution
 *  engine is obtained from the currently set, see \ref randomnumbers. Hence, you
 *  can manipulate both in order to get different results or to set the seed such that
 *  some "randomly" drawn value always work well (e.g. for testing). 
 *  
 *  In any case, the FragmentationFitPotentialAction has the option 
 *  "take-best-of" to allow for multiple fits where the best (in terms of l2 error) 
 *  is taken eventually. Furthermore, you can use the "set-threshold" option to
 *  repeat the fit procedure until the L2 error has dropped below the given
 *  threshold.
 *
 * \section potentials-howto-add Howto add new potentials
 *
 *  Adding a new potential requires the following:
 *  -# Add the new modules to Potentials/Specifics
 *  -# Add a unit test on the potential in Potentials/Specifics/unittests
 *  -# Give the potential a type name and add it to PotentialTypes.def. Note
 *     that the name must not contain white space.
 *  -# Add the potential name as case to PotentialFactory such that it knows
 *     how to instantiate your new potential when requested.
 *  -# Remember to use the the RandomNumberGenerator for getting random starting
 *     values!
 *
 * PotentialTypes.def contains a boost::preprocessor sequence of all
 * potential names. PotentialFactory uses this sequence to build its enum to
 * type map and inverse which the user sees when specifying the potential to
 * fit via PotentialTypeValidator.
 *
 *
 * \date 2017-05-14
 */