//
// edip_single.cpp
//
// LGPL Version 2.1 HEADER START
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library 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
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
// MA 02110-1301  USA
//
// LGPL Version 2.1 HEADER END
//

//
// Copyright (c) 2021, Regents of the University of Minnesota.
// All rights reserved.
//
// Contributors:
//    Yaser Afshar
//

#include "edip_single.hpp"

#include <cmath>
#include <cstdio>
#include <cstring>
#include <limits>
#include <memory>
#include <string>
#include <vector>

#include "helper.hpp"
#include "special.hpp"
#include "utility.hpp"

namespace edip_single {
static constexpr double kGridDensityInverse{1.0 /
                                            static_cast<double>(kGridDensity)};
static constexpr double kLeftLimitToZero{-std::numeric_limits<double>::min() *
                                         1000.0};

/*!
 * \brief EDIP number of parameters per line.
 *
 * EDIP parameter file. One set of params can span multiple lines.
 * Only store the 1st entry.
 *
 * format of an entry (one or more lines) in a parameter file
 *
 *     element1   element2   element3
 *     A          B          cutoffA   cutoffC   alpha   beta   eta
 *     gamma      lambda     mu        rho       sigma   Q0
 *     u1         u2         u3        u4
 *
 * units for each parameters:
 *
 *     A, lambda -> are in eV
 *     B, cutoffA, cutoffC, gamma, sigma, -> are in Angstrom
 *     alpha, beta, eta, mu, rho, Q0, u1-u4, -> are pure numbers
 *
 */
static constexpr int kParamsPerLine{20};
}  // namespace edip_single

int EDIPSingle::ProcessParameterFile(
    std::FILE *const parameter_file_pointer, int const max_line_size,
    std::vector<std::string> const &element_name) {
  std::unique_ptr<char> next_line_(new char[max_line_size]);
  char *next_line = next_line_.get();

  // Create a utility object
  Utility ut;

  // Read each set of parameters from the EDIP parameter file.
  // One set of params can span multiple lines. Only store the 1st entry.
  char *words[edip_single::kParamsPerLine + 1];

  while (true) {
    if (ut.GetNextLine(parameter_file_pointer, next_line, max_line_size)) {
      std::string msg{"End of file while reading a line from the `edip` "};
      msg += "parameter file.\n'" + element_name[0];
      msg += "' parameters are not provided.";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    int nwords = ut.GetNumberOfWordsInLine(next_line);

    // concatenate additional lines until have kParamsPerLine words
    while (nwords < edip_single::kParamsPerLine) {
      std::size_t const n = std::strlen(next_line);

      if (ut.GetNextLine(parameter_file_pointer, &next_line[n],
                         max_line_size - static_cast<int>(n))) {
        std::string msg{"End of file while reading a line from the `edip` "};
        msg += "parameter file.\nIncorrect format in the `edip` parameter ";
        msg += "file.\n" + std::to_string(nwords) + " parameters are ";
        msg += "provided while the model expects 20 parameters.";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      nwords = ut.GetNumberOfWordsInLine(next_line);
    }  // nwords < kParamsPerLine

    if (nwords != edip_single::kParamsPerLine) {
      std::string msg{"Incorrect format in the `edip` parameter file.\n"};
      msg += std::to_string(nwords) + " parameters are provided while the ";
      msg += "model expects 20 parameters.";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    nwords = 0;
    words[nwords++] = std::strtok(next_line, "' \t\n\r\f");
    while ((words[nwords] = std::strtok(nullptr, "' \t\n\r\f"))) {
      ++nwords;
    }

    // ielement,jelement,kelement = 1st args
    // if all 3 args are in element list, then parse this line
    // else skip to the next entry in file
    std::string ielement(words[0]);
    std::string jelement(words[1]);
    std::string kelement(words[2]);
    if (ielement != element_name[0] || jelement != element_name[0] ||
        kelement != element_name[0]) {
      continue;
    }

    for (int i = 3; i < edip_single::kParamsPerLine; ++i) {
      if (!ut.IsDouble(words[i])) {
        std::string msg(words[i]);
        msg += " is not a valid floating-point number.";
        HELPER_LOG_ERROR(msg);
        return true;
      }
    }  // Loop over i

    params_.A = std::atof(words[3]);
    params_.B = std::atof(words[4]);
    params_.cutoffA = std::atof(words[5]);
    params_.cutoffC = std::atof(words[6]);
    params_.alpha = std::atof(words[7]);
    params_.beta = std::atof(words[8]);
    params_.eta = std::atof(words[9]);
    params_.gamma = std::atof(words[10]);
    params_.lambda = std::atof(words[11]);
    params_.mu = std::atof(words[12]);
    params_.rho = std::atof(words[13]);
    params_.sigma = std::atof(words[14]);
    params_.Q0 = std::atof(words[15]);
    params_.u1 = std::atof(words[16]);
    params_.u2 = std::atof(words[17]);
    params_.u3 = std::atof(words[18]);
    params_.u4 = std::atof(words[19]);

    if (params_.A < 0.0 || params_.B < 0.0 || params_.cutoffA < 0.0 ||
        params_.cutoffC < 0.0 || params_.alpha < 0.0 || params_.beta < 0.0 ||
        params_.eta < 0.0 || params_.gamma < 0.0 || params_.lambda < 0.0 ||
        params_.mu < 0.0 || params_.rho < 0.0 || params_.sigma < 0.0) {
      HELPER_LOG_ERROR("Illegal EDIP parameter.");
      return true;
    }

    break;
  }  // loop

  // everything is good
  return false;
}

void EDIPSingle::ConvertUnit(double const convert_length_factor,
                             double const convert_energy_factor) {
  // (Angstroms in metal units).
  if (special::IsNotOne(convert_length_factor)) {
    convert_length_factor_ = convert_length_factor;

    // B, cutoffA, cutoffC, gamma, sigma, -> are in Angstrom
    params_.B *= convert_length_factor;
    params_.cutoffA *= convert_length_factor;
    params_.cutoffC *= convert_length_factor;
    params_.gamma *= convert_length_factor;
    params_.sigma *= convert_length_factor;

    if (convert_length_factor < 1.0) {
      // This conversion is based on some rule of thumbs
      GridStart_ = params_.cutoffC / 15.0;
    }
  }

  // (eV in the case of metal units).
  if (special::IsNotOne(convert_energy_factor)) {
    convert_energy_factor_ = convert_energy_factor;

    // A, lambda -> are in eV
    params_.A *= convert_energy_factor;
    params_.lambda *= convert_energy_factor;
  }
}

void EDIPSingle::CompleteSetup(double *max_cutoff) {
  // set cutoff square
  params_.cutsq = params_.cutoffA * params_.cutoffA;
  cutmax_ = params_.cutoffA;
  *max_cutoff = cutmax_;

  // Grow local arrays
  Resize();
}

void EDIPSingle::Resize() {
  double const A{params_.A};
  double const B{params_.B};
  double const cutoffA{params_.cutoffA};
  double const cutoffC{params_.cutoffC};
  double const alpha{params_.alpha};
  double const sigma{params_.sigma};
  double const rho{params_.rho};
  double const gamma{params_.gamma};

  // cutoffFunction
  double const cg = cutoffC - GridStart_;
  double const ac = cutoffA - cutoffC;

  std::size_t const number_grid_points_one_cuttoff =
      static_cast<std::size_t>(cg * edip_single::kGridDensity);

  std::size_t const number_grid_points_not_one_cuttoff =
      static_cast<std::size_t>(ac * edip_single::kGridDensity);

  std::size_t const number_grid_points =
      number_grid_points_one_cuttoff + number_grid_points_not_one_cuttoff + 2;

  cutoffFunction_.resize(number_grid_points);
  cutoffFunctionDerived_.resize(number_grid_points);

  // Init cutoffFunction
  {
    for (std::size_t l = 0; l < number_grid_points_one_cuttoff; ++l) {
      cutoffFunction_[l] = 1.0;
    }

    for (std::size_t l = 0; l < number_grid_points_one_cuttoff; ++l) {
      cutoffFunctionDerived_[l] = 0.0;
    }

    double r = GridStart_ + number_grid_points_one_cuttoff *
                                edip_single::kGridDensityInverse;

    for (auto l = number_grid_points_one_cuttoff; l < number_grid_points; ++l) {
      double const ac_rc = ac / (r - cutoffC);
      double const ac_rc3 = ac_rc * ac_rc * ac_rc;
      double const ac_rc4 = ac_rc3 * ac_rc;
      double const onemac_rc3 = 1.0 - ac_rc3;
      double const alpha_1cst3 = alpha / onemac_rc3;
      double const ealpha_1cst3 = std::exp(alpha_1cst3);
      double const c1 = -3.0 * alpha / ac;
      double const c2 = ac_rc4 / (onemac_rc3 * onemac_rc3);
      cutoffFunction_[l] = ealpha_1cst3;
      cutoffFunctionDerived_[l] = c1 * c2 * ealpha_1cst3;
      r += edip_single::kGridDensityInverse;
    }
  }

  // pow2B
  double const cst = cutoffA + edip_single::kLeftLimitToZero - GridStart_;
  std::size_t const number_grid_points_R =
      static_cast<std::size_t>(cst * edip_single::kGridDensity) + 2;

  pow2B_.resize(number_grid_points_R);
  exp2B_.resize(number_grid_points_R);
  exp3B_.resize(number_grid_points_R);

  // Init pow2B
  {
    std::size_t l;
    double r = GridStart_;
    for (l = 0; l < number_grid_points_R - 2; ++l) {
      double const br = B / r;
      pow2B_[l] = std::pow(br, rho);
      double const ra = r - cutoffA;
      double const sr = sigma / ra;
      exp2B_[l] = A * std::exp(sr);
      double const gr = gamma / ra;
      exp3B_[l] = std::exp(gr);
      r += edip_single::kGridDensityInverse;
    }

    double const br1 = B / r;
    pow2B_[l] = std::pow(br1, rho);
    exp2B_[l] = 0.0;
    exp3B_[l] = 0.0;
    ++l;

    r += edip_single::kGridDensityInverse;

    double const br2 = B / r;
    pow2B_[l] = std::pow(br2, rho);
    exp2B_[l] = 0.0;
    exp3B_[l] = 0.0;
  }
}

void EDIPSingle::Resize(std::size_t const max_number_of_neighbors) {
  std::size_t const number_grid_points =
      max_number_of_neighbors * edip_single::kGridDensity + 2;

  if (number_grid_points <= tauFunctionGrid_.size()) {
    return;
  }

  // tauFunctionGrid
  tauFunctionGrid_.resize(number_grid_points);
  tauFunctionDerivedGrid_.resize(number_grid_points);

  // Init tauFunctionGrid
  {
    double const u1{params_.u1};
    double const u2{params_.u2};
    double const u3{params_.u3};
    double const u4{params_.u4};
    double const u2u3 = u2 * u3;
    double const u2u4 = u2 * u4;
    double const twou2u4 = 2.0 * u2u4;
    double const u2u3u4 = -u2u3 * u4;
    double r{0.0};
    for (std::size_t l = 0; l < number_grid_points; ++l) {
      double const u4r = -u4 * r;
      double const eu4r = std::exp(u4r);
      double const twou4r = 2.0 * u4r;
      double const etwou4r = std::exp(twou4r);
      double const t = u1 + u2u3 * eu4r - u2 * etwou4r;
      double const dt = u2u3u4 * eu4r + twou2u4 * etwou4r;
      tauFunctionGrid_[l] = t;
      tauFunctionDerivedGrid_[l] = dt;
      r += edip_single::kGridDensityInverse;
    }
  }

  // expMinusBetaZeta_iZeta_iGrid
  expMinusBetaZeta_iZeta_iGrid_.resize(number_grid_points);

  // Init expMinusBetaZeta_iZeta_iGrid
  {
    double const beta{-params_.beta};
    double r{0.0};
    for (std::size_t l = 0; l < number_grid_points; ++l) {
      double const cst = beta * r * r;
      expMinusBetaZeta_iZeta_iGrid_[l] = std::exp(cst);
      r += edip_single::kGridDensityInverse;
    }
  }

  // qFunctionGrid
  qFunctionGrid_.resize(number_grid_points);

  // Init qFunctionGrid
  {
    double const Q0{params_.Q0};
    double const mu{-params_.mu};
    double r{0.0};
    for (std::size_t l = 0; l < qFunctionGrid_.size(); ++l) {
      qFunctionGrid_[l] = Q0 * std::exp(mu * r);
      r += edip_single::kGridDensityInverse;
    }
  }

  // Temporary arrays
  pre_force_rij_inverse_.resize(max_number_of_neighbors);
  pre_force_Exp3B_ij_.resize(max_number_of_neighbors);
  pre_force_Exp3BDerived_ij_.resize(max_number_of_neighbors);
  pre_force_Exp2B_ij_.resize(max_number_of_neighbors);
  pre_force_Exp2BDerived_ij_.resize(max_number_of_neighbors);
  pre_force_Pow2B_ij_.resize(max_number_of_neighbors);
  pre_force_coord_data_.resize(max_number_of_neighbors * 5);
}