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

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


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

#include "CodePatterns/MemDebug.hpp"

#include "ManyBodyPotential_Tersoff.hpp"

#include <boost/bind.hpp>
#include <cmath>

#include "CodePatterns/Assert.hpp"

#include "Potentials/helpers.hpp"

ManyBodyPotential_Tersoff::results_t
ManyBodyPotential_Tersoff::operator()(
    const arguments_t &arguments
    ) const
{
  const double &distance = arguments[0].distance;
  const double cutoff = function_cutoff(distance);
  const double result = (cutoff == 0.) ? 0. : cutoff * (
      function_prefactor(
          manybodyparameter_alpha,
          boost::bind(&ManyBodyPotential_Tersoff::function_eta,
              boost::cref(*this),
              boost::cref(arguments[0])))
      * function_smoother(
          distance,
          manybodyparameter_A,
          manybodyparameter_lambda1)
      +
      function_prefactor(
          manybodyparameter_beta,
          boost::bind(&ManyBodyPotential_Tersoff::function_zeta,
              boost::cref(*this),
              boost::cref(arguments[0])))
      * function_smoother(
          distance,
          -manybodyparameter_B,
          manybodyparameter_lambda2)
  );
  return std::vector<result_t>(1, result);
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_cutoff(
    const double &distance
  ) const
{
  const double offset = (distance - cutoff_offset);
  if (offset < - cutoff_halfwidth)
    return 1.;
  else if (offset > cutoff_halfwidth)
    return 0.;
  else {
    return (0.5 - 0.5 * sin( .5 * M_PI * offset/cutoff_halfwidth));
  }
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_prefactor(
    const double &alpha,
    boost::function<result_t()> etafunction
  ) const
{
  return pow(
      (1. + Helpers::pow(alpha * etafunction(), manybodyparameter_n)),
      -1./(2.*manybodyparameter_n));
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_eta(
    const argument_t &r_ij
  ) const
{
  result_t result = 0.;

  // get all triples within the cutoff
  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
      iter != triples.end(); ++iter) {
    ASSERT( iter->size() == 2,
        "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain of exactly two distances.");
    const argument_t &r_ik = (*iter)[0];
    result += function_cutoff(r_ik.distance)
        * exp( Helpers::pow(manybodyparameter_lambda3 * (r_ij.distance - r_ik.distance) ,3));
  }

  return result;
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_zeta(
    const argument_t &r_ij
  ) const
{
  result_t result = 0.;

  // get all triples within the cutoff
  std::vector<arguments_t> triples = triplefunction(r_ij, cutoff_offset+cutoff_halfwidth);
  for (std::vector<arguments_t>::const_iterator iter = triples.begin();
      iter != triples.end(); ++iter) {
    ASSERT( iter->size() == 2,
        "ManyBodyPotential_Tersoff::function_zeta() - the triples result must contain exactly two distances.");
    const argument_t &r_ik = (*iter)[0];
    const argument_t &r_jk = (*iter)[1];
    result +=
        function_cutoff(r_ik.distance)
        * function_angle(r_ij.distance, r_ik.distance, r_jk.distance)
        * exp( Helpers::pow(manybodyparameter_lambda3 * (r_ij.distance - r_ik.distance) ,3));
  }

  return result;
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_angle(
    const double &r_ij,
    const double &r_ik,
    const double &r_jk
  ) const
{
  const double angle = Helpers::pow(r_ij,2) + Helpers::pow(r_ik,2) - Helpers::pow(r_jk,2);
  const double divisor = r_ij * r_ik;
  const double result =
      1.
      + (Helpers::pow(manybodyparameter_c, 2)/Helpers::pow(manybodyparameter_d, 2))
      - Helpers::pow(manybodyparameter_c, 2)/(Helpers::pow(manybodyparameter_d, 2) +
          Helpers::pow(manybodyparameter_h - cos(angle/divisor),2));
  return result;
}