Index: src/Potentials/Specifics/SaturationPotential.cpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.cpp	(revision 4ffbb79d92e61ca142eecc9e02d8d3c22ff4e125)
+++ src/Potentials/Specifics/SaturationPotential.cpp	(revision 4ffbb79d92e61ca142eecc9e02d8d3c22ff4e125)
@@ -0,0 +1,245 @@
+/*
+ * Project: MoleCuilder
+ * Description: creates and alters molecular systems
+ * Copyright (C)  2012 University of Bonn. All rights reserved.
+ * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
+ * 
+ *
+ *   This file is part of MoleCuilder.
+ *
+ *    MoleCuilder is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 2 of the License, or
+ *    (at your option) any later version.
+ *
+ *    MoleCuilder is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with MoleCuilder.  If not, see <http://www.gnu.org/licenses/>. 
+ */
+
+/*
+ * SaturationPotential.cpp
+ *
+ *  Created on: Oct 11, 2012
+ *      Author: heber
+ */
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "CodePatterns/MemDebug.hpp"
+
+#include "SaturationPotential.hpp"
+
+#include "CodePatterns/Assert.hpp"
+
+#include "Potentials/helpers.hpp"
+
+SaturationPotential::SaturationPotential(
+    boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
+  energy_offset(0.),
+  triplefunction(_triplefunction),
+  saturation_cutoff(2.5)
+{}
+
+SaturationPotential::SaturationPotential(
+    const double _morse_spring_constant,
+    const double _morse_equilibrium_distance,
+    const double _morse_dissociation_energy,
+    const double _angle_spring_constant,
+    const double _angle_equilibrium_distance,
+    const double _all_energy_offset,
+    boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
+  energy_offset(_all_energy_offset),
+  triplefunction(_triplefunction),
+  saturation_cutoff(2.5)
+{
+  parameters_t morse_params(morse.getParameterDimension());
+  morse_params[PairPotential_Morse::spring_constant] = _morse_spring_constant;
+  morse_params[PairPotential_Morse::equilibrium_distance] = _morse_equilibrium_distance;
+  morse_params[PairPotential_Morse::dissociation_energy] = _morse_dissociation_energy;
+  morse_params[PairPotential_Morse::energy_offset] = 0.;
+  morse.setParameters(morse_params);
+  parameters_t angle_params(angle.getParameterDimension());
+  angle_params[PairPotential_Angle::spring_constant] = _angle_spring_constant;
+  angle_params[PairPotential_Angle::equilibrium_distance] = _angle_equilibrium_distance;
+  angle_params[PairPotential_Angle::energy_offset] = 0.;
+  angle.setParameters(angle_params);
+}
+
+void SaturationPotential::setParameters(const parameters_t &_params)
+{
+  const size_t paramsDim = _params.size();
+  ASSERT( paramsDim <= getParameterDimension(),
+      "SaturationPotential::setParameters() - we need not more than "
+      +toString(getParameterDimension())+" parameters.");
+//    LOG(1, "INFO: Setting new SaturationPotential params: " << _params);
+
+
+  // offsets
+  if (paramsDim > all_energy_offset)
+    energy_offset = _params[all_energy_offset];
+
+  // Morse
+  {
+    parameters_t morse_params(morse.getParameters());
+    if (paramsDim > morse_spring_constant)
+      morse_params[PairPotential_Morse::spring_constant] = _params[morse_spring_constant];
+    if (paramsDim > morse_equilibrium_distance)
+      morse_params[PairPotential_Morse::equilibrium_distance] = _params[morse_equilibrium_distance];
+    if (paramsDim > morse_dissociation_energy)
+      morse_params[PairPotential_Morse::dissociation_energy] = _params[morse_dissociation_energy];
+    morse_params[PairPotential_Morse::energy_offset] = 0.;
+    morse.setParameters(morse_params);
+  }
+
+  // Angle
+  {
+    parameters_t angle_params(angle.getParameters());
+    if (paramsDim > angle_spring_constant)
+      angle_params[PairPotential_Angle::spring_constant] = _params[angle_spring_constant];
+    if (paramsDim > angle_equilibrium_distance)
+      angle_params[PairPotential_Angle::equilibrium_distance] = _params[angle_equilibrium_distance];
+    angle_params[PairPotential_Angle::energy_offset] = 0.;
+    angle.setParameters(angle_params);
+  }
+#ifndef NDEBUG
+  parameters_t check_params(getParameters());
+  check_params.resize(paramsDim); // truncate to same size
+  ASSERT( check_params == _params,
+      "SaturationPotential::setParameters() - failed, mismatch in to be set "
+      +toString(_params)+" and set "+toString(check_params)+" params.");
+#endif
+}
+
+SaturationPotential::parameters_t SaturationPotential::getParameters() const
+{
+  parameters_t params(getParameterDimension());
+  const parameters_t morse_params = morse.getParameters();
+  const parameters_t angle_params = angle.getParameters();
+
+  params[all_energy_offset] = energy_offset;
+
+  params[morse_spring_constant] = morse_params[PairPotential_Morse::spring_constant];
+  params[morse_equilibrium_distance] = morse_params[PairPotential_Morse::equilibrium_distance];
+  params[morse_dissociation_energy] = morse_params[PairPotential_Morse::dissociation_energy];
+
+  params[angle_spring_constant] = angle_params[PairPotential_Angle::spring_constant];
+  params[angle_equilibrium_distance] = angle_params[PairPotential_Angle::equilibrium_distance];
+  return params;
+}
+
+SaturationPotential::results_t
+SaturationPotential::operator()(
+    const arguments_t &arguments
+    ) const
+{
+  double result = 0.;
+  for(arguments_t::const_iterator argiter = arguments.begin();
+      argiter != arguments.end();
+      ++argiter) {
+    const argument_t &r_ij = *argiter;
+    arguments_t args(1, r_ij);
+    if ((r_ij.indices.first == 0)) { // || (r_ij.indices.second == 0)) {
+      result += morse(args)[0];
+      std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
+      for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+          iter != triples.end(); ++iter) {
+        ASSERT( iter->size() == 2,
+            "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
+        const argument_t &r_ik = (*iter)[0];
+        const argument_t &r_jk = (*iter)[1];
+        arguments_t args;
+        args.push_back(r_ij);
+        args.push_back(r_ik);
+        args.push_back(r_jk);
+        result += .5 * angle(args)[0];
+      }
+    }
+  }
+  return std::vector<result_t>(1, energy_offset + result);
+}
+
+SaturationPotential::derivative_components_t
+SaturationPotential::derivative(
+    const arguments_t &arguments
+    ) const
+{
+  ASSERT( 0,
+      "SaturationPotential::operator() - not implemented.");
+  derivative_components_t result;
+  return result;
+}
+
+SaturationPotential::results_t
+SaturationPotential::parameter_derivative(
+    const arguments_t &arguments,
+    const size_t index
+    ) const
+{
+  double result = 0.;
+  for(arguments_t::const_iterator argiter = arguments.begin();
+      argiter != arguments.end();
+      ++argiter) {
+    const argument_t &r_ij = *argiter;
+    arguments_t args(1, r_ij);
+    if ((r_ij.indices.first == 0)) { // || (r_ij.indices.second == 0)) {
+      switch (index) {
+        case morse_spring_constant:
+        {
+          result += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
+          break;
+        }
+        case morse_equilibrium_distance:
+        {
+          result += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
+          break;
+        }
+        case morse_dissociation_energy:
+        {
+          result += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
+          break;
+        }
+        default:
+        {
+          std::vector<arguments_t> triples = triplefunction(r_ij, saturation_cutoff);
+          for (std::vector<arguments_t>::const_iterator iter = triples.begin();
+              iter != triples.end(); ++iter) {
+            ASSERT( iter->size() == 2,
+                "ManyBodyPotential_Tersoff::function_derivative_c() - the triples result must contain exactly two distances.");
+            const argument_t &r_ik = (*iter)[0];
+            const argument_t &r_jk = (*iter)[1];
+            arguments_t args;
+            args.push_back(r_ij);
+            args.push_back(r_ik);
+            args.push_back(r_jk);
+            switch (index) {
+              case angle_spring_constant:
+              {
+                result += .5*angle.parameter_derivative(args, PairPotential_Angle::spring_constant)[0];
+                break;
+              }
+              case angle_equilibrium_distance:
+              {
+                result += .5*angle.parameter_derivative(args, PairPotential_Angle::equilibrium_distance)[0];
+                break;
+              }
+              default:
+                break;
+            }
+          }
+        }
+      }
+    }
+  }
+  if (index == all_energy_offset)
+    result = 1.;
+  return SaturationPotential::results_t(1, result);
+}
Index: src/Potentials/Specifics/SaturationPotential.hpp
===================================================================
--- src/Potentials/Specifics/SaturationPotential.hpp	(revision 4ffbb79d92e61ca142eecc9e02d8d3c22ff4e125)
+++ src/Potentials/Specifics/SaturationPotential.hpp	(revision 4ffbb79d92e61ca142eecc9e02d8d3c22ff4e125)
@@ -0,0 +1,148 @@
+/*
+ * SaturationPotential.hpp
+ *
+ *  Created on: Oct 11, 2012
+ *      Author: heber
+ */
+
+#ifndef SATURATIONPOTENTIAL_HPP_
+#define SATURATIONPOTENTIAL_HPP_
+
+
+// include config.h
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <boost/function.hpp>
+#include <limits>
+
+#include "Potentials/EmpiricalPotential.hpp"
+#include "Potentials/Specifics/PairPotential_Angle.hpp"
+#include "Potentials/Specifics/PairPotential_Morse.hpp"
+#include "FunctionApproximation/FunctionModel.hpp"
+
+/** This is a combination of Morse and Angle potentials for saturated elements.
+ *
+ */
+class SaturationPotential : virtual public EmpiricalPotential, virtual public FunctionModel
+{
+  //!> grant unit test access to internal parts
+  friend class SaturationPotentialTest;
+  // some repeated typedefs to avoid ambiguities
+  typedef FunctionModel::arguments_t arguments_t;
+  typedef FunctionModel::result_t result_t;
+  typedef FunctionModel::results_t results_t;
+  typedef EmpiricalPotential::derivative_components_t derivative_components_t;
+  typedef FunctionModel::parameters_t parameters_t;
+public:
+  SaturationPotential(
+      boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction
+      );
+  SaturationPotential(
+      const double _morse_spring_constant,
+      const double _morse_equilibrium_distance,
+      const double _morse_dissociation_energy,
+      const double _angle_spring_constant,
+      const double _angle_equilibrium_distance,
+      const double _energy_offset,
+      boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction
+      );
+  virtual ~SaturationPotential() {}
+
+  /** Setter for parameters as required by FunctionModel interface.
+   *
+   * \param _params given set of parameters
+   */
+  void setParameters(const parameters_t &_params);
+
+  /** Getter for parameters as required by FunctionModel interface.
+   *
+   * \return set of parameters
+   */
+  parameters_t getParameters() const;
+
+  /** Getter for the number of parameters of this model function.
+   *
+   * \return number of parameters
+   */
+  size_t getParameterDimension() const
+  { return MAXPARAMS; }
+
+  /** Evaluates the harmonic potential function for the given arguments.
+   *
+   * @param arguments single distance
+   * @return value of the potential function
+   */
+  results_t operator()(const arguments_t &arguments) const;
+
+  /** Evaluates the derivative of the potential function.
+   *
+   * @param arguments single distance
+   * @return vector with derivative with respect to the input degrees of freedom
+   */
+  derivative_components_t derivative(const arguments_t &arguments) const;
+
+  /** Evaluates the derivative of the function with the given \a arguments
+   * with respect to a specific parameter indicated by \a index.
+   *
+   * \param arguments set of arguments as input variables to the function
+   * \param index derivative of which parameter
+   * \return result vector containing the derivative with respect to the given
+   *         input
+   */
+  results_t parameter_derivative(const arguments_t &arguments, const size_t index) const;
+
+  /** States whether lower and upper boundaries should be used to constraint
+   * the parameter search for this function model.
+   *
+   * \return true - constraints should be used, false - else
+   */
+  bool isBoxConstraint() const {
+    return true;
+  }
+
+  /** Returns a vector which are the lower boundaries for each parameter_t
+   * of this FunctionModel.
+   *
+   * \return vector of parameter_t resembling lowest allowed values
+   */
+  parameters_t getLowerBoxConstraints() const {
+    parameters_t lowerbounds(getParameterDimension(), -std::numeric_limits<double>::max());
+    lowerbounds[morse_equilibrium_distance] = 0.;
+    lowerbounds[angle_equilibrium_distance] = -1.;
+    return lowerbounds;
+  }
+
+  /** Returns a vector which are the upper boundaries for each parameter_t
+   * of this FunctionModel.
+   *
+   * \return vector of parameter_t resembling highest allowed values
+   */
+  parameters_t getUpperBoxConstraints() const {
+    parameters_t upperbounds(getParameterDimension(), std::numeric_limits<double>::max());
+    upperbounds[angle_equilibrium_distance] = 1.;
+    return upperbounds;
+  }
+
+  enum parameter_enum_t {
+    morse_spring_constant,
+    morse_equilibrium_distance,
+    morse_dissociation_energy,
+    angle_spring_constant,
+    angle_equilibrium_distance,
+    all_energy_offset,
+    MAXPARAMS
+  };
+
+private:
+  PairPotential_Morse morse;
+  PairPotential_Angle angle;
+  double energy_offset;
+
+  //!> bound function that obtains the triples for the internal coordinationb summation.
+  const boost::function< std::vector< arguments_t >(const argument_t &, const double)> &triplefunction;
+  const double saturation_cutoff;
+};
+
+#endif /* SATURATIONPOTENTIAL_HPP_ */
