source: src/documentation/constructs/potentials.dox@ 0b6b77

/*
 * 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 thought
 * 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).
 *
 * Empirical potentials have been thought of by fellows such as Lennard-Jones,
 * Morse, Tersoff, Stillinger and Weber, etc. And in their honor, the
 * potential form is 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
 *
 * The terms 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. that has internal
 *    parameters to 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 multimap 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. 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-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::charges_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::getCharges, _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 distributionm
 *  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, For this 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
 *  stop restarting the fit procedure first when 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 2013-04-09
 */