/*
 * ManyBodyPotential_Tersoff.hpp
 *
 *  Created on: Sep 26, 2012
 *      Author: heber
 */

#ifndef MANYBODYPOTENTIAL_TERSOFF_HPP_
#define MANYBODYPOTENTIAL_TERSOFF_HPP_

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

#include <boost/function.hpp>
#include <cmath>
#include <limits>

#include "Potentials/EmpiricalPotential.hpp"
#include "FunctionApproximation/FunctionModel.hpp"

/** This class is the implementation of the Tersoff potential function.
 *
 * \note The arguments_t argument list is here in the following order:
 * -# first \f$ r_{ij} \f$,
 * -# then all \f$ r_{ik} \f$ that are within the cutoff, i.e. \f$ r_{ik} < R + D\f$
 *
 */
class ManyBodyPotential_Tersoff : virtual public EmpiricalPotential, virtual public FunctionModel
{
  //!> grant unit test access to internal parts
  friend class ManyBodyPotential_TersoffTest;
  // 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:
  /** Constructor for class ManyBodyPotential_Tersoff.
   *
   * @param _triplefunction function that returns a list of triples (i.e. the
   *        two remaining distances) to a given pair of points (contained as
   *        indices within the argument)
   */
  ManyBodyPotential_Tersoff(
      boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction
      );

  /** Constructor for class ManyBodyPotential_Tersoff.
   *
   * @param _R offset for cutoff
   * @param _S halfwidth for cutoff relative to \a _R
   * @param A
   * @param B
   * @param lambda
   * @param mu
   * @param lambda3
   * @param alpha
   * @param beta
   * @param chi
   * @param omega
   * @param n
   * @param c
   * @param d
   * @param h
   * @param offset
   * @param _triplefunction function that returns a list of triples (i.e. the
   *        two remaining distances) to a given pair of points (contained as
   *        indices within the argument)
   */
  ManyBodyPotential_Tersoff(
      const double &_R,
      const double &_S,
      const double &_A,
      const double &_B,
      const double &_lambda,
      const double &_mu,
      const double &_lambda3,
      const double &_alpha,
      const double &_beta,
      const double &_chi,
      const double &_omega,
      const double &_n,
      const double &_c,
      const double &_d,
      const double &_h,
      const double &_offset,
      boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction);

  /** Destructor of class ManyBodyPotential_Tersoff.
   *
   */
  virtual ~ManyBodyPotential_Tersoff() {}

  /** Evaluates the Tersoff potential 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 Tersoff potential with respect to the
   * input variables.
   *
   * @param arguments single distance
   * @return vector with components of the derivative
   */
  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 lowerbound(getParameterDimension(), -std::numeric_limits<double>::max());
//    lowerbound[R] = 0.;
//    lowerbound[S] = 0.;
//    lowerbound[lambda3] = 0.;
//    lowerbound[alpha] = 0.;
    lowerbound[beta] = std::numeric_limits<double>::min();
    lowerbound[n] = std::numeric_limits<double>::min();
    lowerbound[c] = std::numeric_limits<double>::min();
    lowerbound[d] = std::numeric_limits<double>::min();
    return lowerbound;
  }

  /** 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 {
    return parameters_t(getParameterDimension(), std::numeric_limits<double>::max());
  }

private:
  /** Prohibit private default constructor.
   *
   * We essentially need the triplefunction, hence without this function cannot
   * be.
   */
  ManyBodyPotential_Tersoff();

private:
  /** This function represents the cutoff \f$ f_C \f$.
   *
   * @param distance variable of the function
   * @return a value in [0,1].
   */
  result_t function_cutoff(
      const double &distance
      ) const;
  /** This  function has the exponential feature from the Morse potential.
   *
   * @param prefactor prefactor parameter to exp function
   * @param lambda scale parameter of exp function's argument
   * @param distance variable of the function
   * @return
   */
  result_t function_smoother(
      const double &prefactor,
      const double &lambda,
      const double &distance
      ) const;

  /** This function represents \f$ (1 + \alpha^n \eta^n)^{-1/2n} \f$.
   *
   * @param alpha prefactor to eta function
   * @param r_ij distance argument
   * @param eta result value of eta or zeta
   * @return \f$ (1 + \alpha^n \eta^n)^{-1/2n} \f$
   */
  result_t function_prefactor(
      const double &alpha,
      const double &eta
        ) const;

  result_t
  function_eta(
      const argument_t &r_ij
    ) const;

  result_t
  function_zeta(
      const argument_t &r_ij
    ) const;

  result_t
  function_theta(
      const double &r_ij,
      const double &r_ik,
      const double &r_jk
    ) const;

  result_t
  function_angle(
      const double &r_ij,
      const double &r_ik,
      const double &r_jk
    ) const;

private:
  result_t
  function_derivative_c(
      const argument_t &r_ij
    ) const;

  result_t
  function_derivative_d(
      const argument_t &r_ij
    ) const;

  result_t
  function_derivative_h(
      const argument_t &r_ij
    ) const;

public:
  enum parameter_enum_t {
    A,
    B,
    lambda,
    mu,
    beta,
    n,
    c,
    d,
    h,
    offset,
//    R,
//    S,
//    lambda3,
//    alpha,
//    chi,
//    omega,
    MAXPARAMS
  };

private:
  //!> parameter vector with parameters as in enum parameter_enum_t
  parameters_t params;

public:
  // some internal parameters which are fixed
  const double R;
  const double S;
  const double lambda3;
  const double alpha;
  const double chi;
  const double omega;

public:
  /** 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
  {
    return params;
  }

  /** Getter for the number of parameters of this model function.
   *
   * \return number of parameters
   */
  size_t getParameterDimension() const
  {
    return MAXPARAMS;
  }

private:
  //!> bound function that obtains the triples for the internal coordinationb summation.
  const boost::function< std::vector< arguments_t >(const argument_t &, const double)> &triplefunction;
};


#endif /* MANYBODYPOTENTIAL_TERSOFF_HPP_ */