source: src/Potentials/Specifics/ManyBodyPotential_Tersoff.hpp

/*
 * 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"

class PotentialFactory;
class TrainingData;

/** 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 :
  public EmpiricalPotential
{
  //!> grant unit test access to internal parts
  friend class ManyBodyPotential_TersoffTest;
  //!> grant PotentialFactory access to default cstor
  friend class PotentialFactory;
  // some repeated typedefs to avoid ambiguities
  typedef FunctionModel::list_of_arguments_t list_of_arguments_t;
  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;
private:
  /** Private default constructor.
   *
   * This prevents creation of potential without set ParticleTypes_t.
   *
   */
  ManyBodyPotential_Tersoff();

  /** Creates the binding model specific for this potential.
   *
   * Private because this is used internally some of the constructors.
   *
   * \param _ParticleTypes particle type for the potential
   * \return binding model
   */
  BindingModel generateBindingModel(const EmpiricalPotential::ParticleTypes_t &_ParticleTypes) const;

public:
  /** Constructor for class ManyBodyPotential_Tersoff.
   *
   * \param _ParticleTypes particle types for this potential
   */
  ManyBodyPotential_Tersoff(
      const ParticleTypes_t &_ParticleTypes
      );

  /** 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 _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 ParticleTypes_t &_ParticleTypes,
      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);

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

  /** Evaluates the Tersoff potential for the given arguments.
   *
   * @param listarguments list of distances
   * @return value of the potential function
   */
  results_t operator()(const list_of_arguments_t &listarguments) const;

  /** Evaluates the derivative of the Tersoff potential with respect to the
   * input variables.
   *
   * @param listarguments list of distances
   * @return vector with components of the derivative
   */
  derivative_components_t derivative(const list_of_arguments_t &listarguments) const;

  /** Evaluates the derivative of the function with the given \a arguments
   * with respect to a specific parameter indicated by \a index.
   *
   * \param listarguments list of distances
   * \param index derivative of which parameter
   * \return result vector containing the derivative with respect to the given
   *         input
   */
  results_t parameter_derivative(const list_of_arguments_t &listarguments, const size_t index) const;

  /** Returns the functor that converts argument_s into the
   * internal coordinate described by this potential function.
   *
   * \return coordinator functor
   */
  Coordinator::ptr getCoordinator() const
  { return coordinator; }

  /** Return the token name of this specific potential.
   *
   * \return token name of the potential
   */
  const std::string& getToken() const
  { return potential_token; }

  /** Returns a vector of parameter names.
   *
   * This is required from the specific implementation
   *
   * \return vector of strings containing parameter names
   */
  const ParameterNames_t& getParameterNames() const
  { return ParameterNames; }

  /** 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());
  }

  /** Returns a bound function to be used with TrainingData, extracting distances
   * from a Fragment.
   *
   * \return bound function extracting distances from a fragment
   */
  FunctionModel::filter_t getSpecificFilter() const;

  /** Returns the number of arguments the underlying function requires.
   *
   * \return number of arguments of the function
   */
  size_t getSpecificArgumentCount() const
  { return 1; }

  /** Sets the magic triple function that we use for getting angle distances.
   *
   * @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)
   */
  void setTriplefunction(triplefunction_t &_triplefunction)
  { triplefunction = _triplefunction;  }

  /** Getter for the graph specifying the binding model of the potential.
   *
   * \return BindingModel ref of the binding model
   */
  const BindingModel& getBindingModel() const
  { return bindingmodel; }

  /**
   * Returns the number of particle types associated with the potential.
   *
   * \return number of particle types
   */
  unsigned int getParticleTypeNumber() const
  { return 2; }

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,
//    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;
  }

  /** Sets the parameter randomly within the sensible range of each parameter.
   *
   * \param data container with training data for guesstimating range
   */
  void setParametersToRandomInitialValues(const TrainingData &data);

  /** 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.
  boost::function< std::vector< arguments_t >(const argument_t &, const double)> triplefunction;

  //!> static definitions of the parameter name for this potential
  static const ParameterNames_t ParameterNames;

  //!> static token of this potential type
  static const std::string potential_token;

  //!> internal coordinator object for converting arguments_t
  static Coordinator::ptr coordinator;

  //!> binding model for this potential
  const BindingModel bindingmodel;
};


#endif /* MANYBODYPOTENTIAL_TERSOFF_HPP_ */