/*
 * 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 "CodePatterns/Info.hpp"

#include "Potentials/helpers.hpp"

ManyBodyPotential_Tersoff::ManyBodyPotential_Tersoff(
    boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction
    ) :
    params(parameters_t(MAXPARAMS, 0.)),
    lambda3(0.),
    alpha(0.),
    chi(1.),
    omega(1.),
    triplefunction(_triplefunction)
{}

ManyBodyPotential_Tersoff::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,
    boost::function< std::vector<arguments_t>(const argument_t &, const double)> &_triplefunction) :
  params(parameters_t(MAXPARAMS, 0.)),
  lambda3(_lambda3),
  alpha(_alpha),
  chi(_chi),
  omega(_mu),
  triplefunction(_triplefunction)
{
//  Info info(__func__);
  params[R] = _R;
  params[S] = _S;
  params[A] = _A;
  params[B] = _B;
  params[lambda] = _lambda;
  params[mu] = _mu;
//  lambda3 = _lambda3;
//  alpha = _alpha;
  params[beta] = _beta;
//  chi = _chi;
//  omega = _omega;
  params[n] = _n;
  params[c] = _c;
  params[d] = _d;
  params[h] = _h;
}

ManyBodyPotential_Tersoff::results_t
ManyBodyPotential_Tersoff::operator()(
    const arguments_t &arguments
    ) const
{
//  Info info(__func__);
  const argument_t &r_ij = arguments[0];
  const double cutoff = function_cutoff(r_ij.distance);
  const double result = (cutoff == 0.) ?
      0. :
      cutoff * (
          function_prefactor(
              alpha,
              function_eta(r_ij))
          * function_smoother(
              params[A],
              params[lambda],
              r_ij.distance)
          +
          function_prefactor(
              params[beta],
              function_zeta(r_ij))
          * function_smoother(
              -params[B],
              params[mu],
              r_ij.distance)
      );
//  LOG(2, "DEBUG: operator()(" << r_ij.distance << ") = " << result);
  return std::vector<result_t>(1, result);
}

ManyBodyPotential_Tersoff::derivative_components_t
ManyBodyPotential_Tersoff::derivative(
    const arguments_t &arguments
    ) const
{
//  Info info(__func__);
  return ManyBodyPotential_Tersoff::derivative_components_t();
}

ManyBodyPotential_Tersoff::results_t
ManyBodyPotential_Tersoff::parameter_derivative(
    const arguments_t &arguments,
    const size_t index
    ) const
{
//  Info info(__func__);
//  ASSERT( arguments.size() == 1,
//      "PairPotential_Harmonic::parameter_derivative() - requires exactly one argument.");
  const argument_t &r_ij = arguments[0];
  switch (index) {
    case R:
    {
      const double result = 0.;
      return results_t(1, result);
      break;
    }
    case S:
    {
      const double result = 0.;
      return results_t(1, result);
      break;
    }
    case A:
    {
      const double result = 0.;
      return results_t(1, result);
      break;
    }
    case B:
    {
      const double result = 0.;
      return results_t(1, result);
      break;
    }
    case lambda:
    {
      const double cutoff = function_cutoff(r_ij.distance);
      const double result = (cutoff == 0.) ?
          0. :
          -r_ij.distance * cutoff * params[lambda] * (
          function_prefactor(
              alpha,
              function_eta(r_ij))
          * function_smoother(
              params[A],
              params[lambda],
              r_ij.distance)
          );
      return results_t(1, result);
      break;
    }
    case mu:
    {
      const double cutoff = function_cutoff(r_ij.distance);
      const double result = (cutoff == 0.) ?
          0. :
          -r_ij.distance * cutoff * params[mu] *(
          function_prefactor(
              params[beta],
              function_zeta(r_ij))
          * function_smoother(
              -params[B],
              params[mu],
              r_ij.distance)
      );
      return results_t(1, result);
      break;
    }
//    case lambda3:
//    {
//      const double result = 0.;
//      return results_t(1, result);
//      break;
//    }
//    case alpha:
//    {
//      const double temp =
//          pow(alpha*function_eta(r_ij), params[n]);
//      const double cutoff = function_cutoff(r_ij.distance);
//      const double result = (cutoff == 0.) || (alpha == 0. )?
//          0. :
//          function_smoother(
//              -params[A],
//              params[lambda],
//              r_ij.distance)
//          * (-.5) * alpha * (temp/alpha)
//          / (1. + temp)
//          ;
//      return results_t(1, result);
//      break;
//    }
//    case chi:
//    {
//      const double result = 0.;
//      return results_t(1, result);
//      break;
//    }
//    case omega:
//    {
//      const double result = 0.;
//      return results_t(1, result);
//      break;
//    }
    case beta:
    {
      const double temp =
          pow(params[beta]*function_zeta(r_ij), params[n]);
      const double cutoff = function_cutoff(r_ij.distance);
      const double result = (cutoff == 0.) || (params[beta] == 0. )?
          0. : cutoff *
          function_smoother(
              -params[B],
              params[mu],
              r_ij.distance)
          * (-.5) * params[beta] * (temp/params[beta])
          / (1. + temp)
          ;
      return results_t(1, result);
      break;
    }
    case n:
    {
      const double temp =
          pow(params[beta]*function_zeta(r_ij), params[n]);
      const double cutoff = function_cutoff(r_ij.distance);
      const double result = (cutoff == 0.) ?
          0. : cutoff *
          function_smoother(
              -params[B],
              params[mu],
              r_ij.distance)
          * params[B]
          * ( log(1.+temp)/(params[n]*params[n]) - temp
              * (log(function_zeta(r_ij)) + log(params[beta]))
              /(params[n]*(1.+temp)))
          ;
      return results_t(1, result);
      break;
    }
    case c:
    {
      const double result = 0.;
      return results_t(1, result);
      break;
    }
    case d:
    {
      const double result = 0.;
      return results_t(1, result);
      break;
    }
    case h:
    {
      const double result = 0.;
      return results_t(1, result);
      break;
    }
    default:
      break;
  }
  return results_t(1, 0.);
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_cutoff(
    const double &distance
  ) const
{
//  Info info(__func__);
  double result = 0.;
  if (distance < params[R])
    result = 1.;
  else if (distance > params[S])
    result = 0.;
  else {
    result = (0.5 + 0.5 * cos( M_PI * (distance - params[R])/(params[S]-params[R])));
  }
//  LOG(2, "DEBUG: function_cutoff(" << distance << ") = " << result);
  return result;
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_prefactor(
    const double &alpha,
    const double &eta
  ) const
{
//  Info info(__func__);
  const double result = chi * pow(
      (1. + pow(alpha * eta, params[n])),
      -1./(2.*params[n]));
//  LOG(2, "DEBUG: function_prefactor(" << alpha << "," << eta << ") = " << result);
  return result;
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_smoother(
    const double &prefactor,
    const double &lambda,
    const double &distance
    ) const
{
//  Info info(__func__);
  const double result = prefactor * exp(-lambda * distance);
//  LOG(2, "DEBUG: function_smoother(" << prefactor << "," << lambda << "," << distance << ") = " << result);
  return result;
}

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

  // get all triples within the cutoff
  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
  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(lambda3 * (r_ij.distance - r_ik.distance) ,3));
  }

//  LOG(2, "DEBUG: function_eta(" << r_ij.distance << ") = " << result);
  return result;
}

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

  // get all triples within the cutoff
  std::vector<arguments_t> triples = triplefunction(r_ij, params[S]);
  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)
        * omega
        * function_angle(r_ij.distance, r_ik.distance, r_jk.distance)
        * exp( Helpers::pow(lambda3 * (r_ij.distance - r_ik.distance) ,3));
  }

//  LOG(2, "DEBUG: function_zeta(" << r_ij.distance << ") = " << result);
  return result;
}

ManyBodyPotential_Tersoff::result_t
ManyBodyPotential_Tersoff::function_angle(
    const double &r_ij,
    const double &r_ik,
    const double &r_jk
  ) const
{
//  Info info(__func__);
  const double angle = Helpers::pow(r_ij,2) + Helpers::pow(r_ik,2) - Helpers::pow(r_jk,2);
  const double divisor = 2.* r_ij * r_ik;
//  LOG(2, "DEBUG: cos(theta)= " << angle/divisor);
  if (divisor != 0.) {
    const double result =
        1.
        + (Helpers::pow(params[c]/params[d], 2))
        - Helpers::pow(params[c], 2)/(Helpers::pow(params[d], 2) +
            Helpers::pow(params[h] - angle/divisor,2));

//    LOG(2, "DEBUG: function_angle(" << r_ij << "," << r_ik << "," << r_jk << ") = " << result);
    return result;
  } else
    return 0.;
}