| 1 | /*
 | 
|---|
| 2 |  * Project: MoleCuilder
 | 
|---|
| 3 |  * Description: creates and alters molecular systems
 | 
|---|
| 4 |  * Copyright (C)  2010 University of Bonn. All rights reserved.
 | 
|---|
| 5 |  * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 | 
|---|
| 6 |  */
 | 
|---|
| 7 | 
 | 
|---|
| 8 | /**
 | 
|---|
| 9 |  * \file potentials.dox
 | 
|---|
| 10 |  *
 | 
|---|
| 11 |  * Created on: Nov 28, 2012
 | 
|---|
| 12 |  *    Author: heber
 | 
|---|
| 13 |  */
 | 
|---|
| 14 | 
 | 
|---|
| 15 | /** \page potentials Empirical Potentials and FunctionModels
 | 
|---|
| 16 |  *
 | 
|---|
| 17 |  * On this page we explain what is meant with the Potentials sub folder.
 | 
|---|
| 18 |  *
 | 
|---|
| 19 |  * First, we are based on fragmenting a molecular system, i.e. dissecting its
 | 
|---|
| 20 |  * bond structure into connected subgraphs, calculating the energies of the
 | 
|---|
| 21 |  * fragments (ab-initio) and summing up to a good approximation of the total
 | 
|---|
| 22 |  * energy of the whole system.
 | 
|---|
| 23 |  * Second, having calculated these energies, there quickly comes up the thought
 | 
|---|
| 24 |  * that one actually calculates quite similar systems all time and if one could
 | 
|---|
| 25 |  * not cache results in an intelligent (i.e. interpolating) fashion ...
 | 
|---|
| 26 |  *
 | 
|---|
| 27 |  * That's where so-called empirical potentials come into play. They are
 | 
|---|
| 28 |  * functions depending on a number of "fitted" parameters and the variable
 | 
|---|
| 29 |  * distances within a molecular fragment (i.e. the bond lengths) in order to
 | 
|---|
| 30 |  * give a value for the total energy without the need to solve a complex
 | 
|---|
| 31 |  * ab-initio model.
 | 
|---|
| 32 |  *
 | 
|---|
| 33 |  * Empirical potentials have been thought of by fellows such as Lennard-Jones,
 | 
|---|
| 34 |  * Morse, Tersoff, Stillinger and Weber, etc. And in their honor, the
 | 
|---|
| 35 |  * potential form is named after its inventor. Hence, we speak e.g. of a
 | 
|---|
| 36 |  * Lennard-Jones potential.
 | 
|---|
| 37 |  *
 | 
|---|
| 38 |  * So, what we have to do in order to cache results is the following procedure:
 | 
|---|
| 39 |  * -# gather similar fragments
 | 
|---|
| 40 |  * -# perform a fit procedure to obtain the parameters for the empirical
 | 
|---|
| 41 |  *    potential
 | 
|---|
| 42 |  * -# evaluate the potential instead of an ab-initio calculation
 | 
|---|
| 43 |  *
 | 
|---|
| 44 |  * The terms we use, model the classes that are implemented:
 | 
|---|
| 45 |  * -# EmpiricalPotential: Contains the interface to a function that can be
 | 
|---|
| 46 |  *    evaluated given a number of arguments_t, i.e. distances. Also, one might
 | 
|---|
| 47 |  *    want to evaluate derivatives.
 | 
|---|
| 48 |  * -# FunctionModel: Is a function that can be fitted, i.e. that has internal
 | 
|---|
| 49 |  *    parameters to be set and got.
 | 
|---|
| 50 |  * -# argument_t: The Argument stores not only the distance but also the index
 | 
|---|
| 51 |  *    pair of the associated atoms and also their charges, to let the potential
 | 
|---|
| 52 |  *    check on validity.
 | 
|---|
| 53 |  * -# SerializablePotential: Eventually, one wants to store to or parse from
 | 
|---|
| 54 |  *    a file all potential parameters. This functionality is encoded in this
 | 
|---|
| 55 |  *    class.
 | 
|---|
| 56 |  * -# HomologyGraph: "Similar" fragments in our case have to have the same bond
 | 
|---|
| 57 |  *    graph. It is stored in the HomologyGraph that acts as representative
 | 
|---|
| 58 |  * -# HomologyContainer: This container combines, in multimap fashion, all
 | 
|---|
| 59 |  *    similar fragments with their energies together, with the HomologyGraph
 | 
|---|
| 60 |  *    as their "key".
 | 
|---|
| 61 |  * -# TrainingData: Here, we combine InputVector and OutputVector that contain
 | 
|---|
| 62 |  *    the set of distances required for the FunctionModel (e.g. only a single
 | 
|---|
| 63 |  *    distance/argument for a pair potential, three for an angle potential,
 | 
|---|
| 64 |  *    etc.) and also the expected OutputVector. This in combination with the
 | 
|---|
| 65 |  *    FunctionModel is the basis for the non-linear regression used for the
 | 
|---|
| 66 |  *    fitting procedure.
 | 
|---|
| 67 |  * -# Extractors: These set of functions yield the set of distances from a
 | 
|---|
| 68 |  *    given fragment that is stored in the HomologyContainer.
 | 
|---|
| 69 |  * -# FunctionApproximation: Contains the interface to the levmar package where
 | 
|---|
| 70 |  *    the Levenberg-Marquardt (Newton + Trust region) algorithm is used to
 | 
|---|
| 71 |  *    perform the fit.
 | 
|---|
| 72 |  *
 | 
|---|
| 73 |  * \section potentials-howto Howto use the potentials
 | 
|---|
| 74 |  *
 | 
|---|
| 75 |  *  We just give a brief run-down in terms of code on how to use the potentials.
 | 
|---|
| 76 |  *  Here, we just describe what to do in order to perform the fitting.
 | 
|---|
| 77 |  *
 | 
|---|
| 78 |  *  \code
 | 
|---|
| 79 |  *  // we need the homology container and the representative graph we want to
 | 
|---|
| 80 |  *  // fit to.
 | 
|---|
| 81 |  *  HomologyContainer homologies;
 | 
|---|
| 82 |  *  const HomologyGraph graph = getSomeGraph(homologies);
 | 
|---|
| 83 |  *  Fragment::charges_t h2o;
 | 
|---|
| 84 |  *  h2o += 8,1,1;
 | 
|---|
| 85 |  *  // TrainingData needs so called Extractors to get the required distances
 | 
|---|
| 86 |  *  // from the stored fragment. These are functions are bound.
 | 
|---|
| 87 |  *  TrainingData AngleData(
 | 
|---|
| 88 |  *      boost::bind(&Extractors::gatherDistancesFromFragment,
 | 
|---|
| 89 |  *          boost::bind(&Fragment::getPositions, _1),
 | 
|---|
| 90 |  *          boost::bind(&Fragment::getCharges, _1),
 | 
|---|
| 91 |  *          boost::cref(h2o),
 | 
|---|
| 92 |  *          _2)
 | 
|---|
| 93 |  *      );
 | 
|---|
| 94 |  *  // now we extract the distances and energies and store them
 | 
|---|
| 95 |  *  AngleData(homologies.getHomologousGraphs(graph));
 | 
|---|
| 96 |  *  // give ParticleTypes of this potential to make it unique
 | 
|---|
| 97 |  *  PairPotential_Angle::ParticleTypes_t types =
 | 
|---|
| 98 |  *        boost::assign::list_of<PairPotential_Angle::ParticleType_t>
 | 
|---|
| 99 |  *        (8)(1)(1)
 | 
|---|
| 100 |  *      ;
 | 
|---|
| 101 |  *  PairPotential_Angle angle(types);
 | 
|---|
| 102 |  *  // give initial parameter
 | 
|---|
| 103 |  *  FunctionModel::parameters_t params(PairPotential_Angle::MAXPARAMS, 0.);
 | 
|---|
| 104 |  *  ... set some initial parameters
 | 
|---|
| 105 |  *  angle.setParameters(params);
 | 
|---|
| 106 |  *
 | 
|---|
| 107 |  *  // use the potential as a FunctionModel along with prepared TrainingData
 | 
|---|
| 108 |  *  FunctionModel &model = angle;
 | 
|---|
| 109 |  *  FunctionApproximation approximator(AngleData, model);
 | 
|---|
| 110 |  *  approximator(FunctionApproximation::ParameterDerivative);
 | 
|---|
| 111 |  *
 | 
|---|
| 112 |  *  // obtain resulting parameters and check remaining L_2 and L_max error
 | 
|---|
| 113 |  *  angleparams = model.getParameters();
 | 
|---|
| 114 |  *  LOG(1, "INFO: L2sum = " << AngleData(model)
 | 
|---|
| 115 |  *      << ", LMax = " << AngleData(model) << ".");
 | 
|---|
| 116 |  *  \endcode
 | 
|---|
| 117 |  *
 | 
|---|
| 118 |  *  The evaluation of the fitted potential is then trivial, e.g.
 | 
|---|
| 119 |  *  \code
 | 
|---|
| 120 |  *  // constructed someplace
 | 
|---|
| 121 |  *  PairPotential_Angle angle(...);
 | 
|---|
| 122 |  *
 | 
|---|
| 123 |  *  // evaluate
 | 
|---|
| 124 |  *  FunctionModel::arguments_t args;
 | 
|---|
| 125 |  *  .. initialise args to the desired distances
 | 
|---|
| 126 |  *  const double value = angle(args)[0]; // output is a vector!
 | 
|---|
| 127 |  *  \endcode
 | 
|---|
| 128 |  *
 | 
|---|
| 129 |  * \date 2012-11-28
 | 
|---|
| 130 |  */
 | 
|---|