//
// zbl.hpp
//
// 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) 2020-2021, Regents of the University of Minnesota.
// All rights reserved.
//
// Contributors:
//    Yaser Afshar
//

#ifndef ZBL_HPP
#define ZBL_HPP

#include <cmath>

#include "special.hpp"

/*!
 * \brief ZBL pair interaction style
 *
 * It computes the Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion
 * for describing high-energy collisions between atoms.
 *
 * The ZBL interaction is already smoothed to 0.0 at the cutoff.
 */
class ZBL {
 public:
  /*! Construct a new ZBL object */
  ZBL();

  /*! Set the coefficients */
  inline void SetCoeff();

  /*!
   * \brief Convert units of the parameters
   *
   * \param convert_length_factor length unit conversion factor
   * \param convert_energy_factor energy unit conversion factor
   */
  inline void ConvertUnit(double const convert_length_factor,
                          double const convert_energy_factor);

  /*!
   * \brief Compute ZBL pair energy
   *
   * \param r Pair atoms distance
   * \param r_inverse inverse of the pair atoms distance
   * \param dphidr derivative of the ZBL pair energy with respect to distance
   * \return double
   */
  inline double ComputePhi(double const r, double const r_inverse,
                           double &dphidr) const;

 private:
  /*! Carbon atomic number */
  double const carbon_atomic_number_{6.0};
  /*! Conversion of q^2/r to energy */
  double qqr2e_{14.399645354084361};
  /*! ZBL constant with length unit */
  double a0_{0.46850};
  /*! ZBL constants */
  double d1a_;
  double d2a_;
  double d3a_;
  double d4a_;
  double zze_;
};

inline ZBL::ZBL() { SetCoeff(); }

inline void ZBL::SetCoeff() {
  double const atomic_number_pow = std::pow(carbon_atomic_number_, 0.230);
  double const ainv = (atomic_number_pow + atomic_number_pow) / a0_;
  d1a_ = 0.201620 * ainv;
  d2a_ = 0.402900 * ainv;
  d3a_ = 0.942290 * ainv;
  d4a_ = 3.199800 * ainv;
  zze_ = carbon_atomic_number_ * carbon_atomic_number_ * qqr2e_;
}

inline void ZBL::ConvertUnit(double const convert_length_factor,
                             double const convert_energy_factor) {
  if (special::IsNotOne(convert_length_factor)) {
    qqr2e_ *= convert_length_factor;
    a0_ *= convert_length_factor;
  }
  if (special::IsNotOne(convert_energy_factor)) {
    qqr2e_ *= convert_energy_factor;
  }
  SetCoeff();
}

inline double ZBL::ComputePhi(double const r, double const r_inverse,
                              double &dphidr) const {
  double const sum1 = 0.028170 * std::exp(-d1a_ * r);
  double const sum2 = 0.280220 * std::exp(-d2a_ * r);
  double const sum3 = 0.509860 * std::exp(-d3a_ * r);
  double const sum4 = 0.181750 * std::exp(-d4a_ * r);
  double const sum = sum1 + sum2 + sum3 + sum4;
  double const dsum1 = -sum1 * d1a_;
  double const dsum2 = -sum2 * d2a_;
  double const dsum3 = -sum3 * d3a_;
  double const dsum4 = -sum4 * d4a_;
  double const dsum = dsum1 + dsum2 + dsum3 + dsum4;
  double const phi = zze_ * sum * r_inverse;
  dphidr = (zze_ * dsum - phi) * r_inverse;
  return phi;
}

#endif  // ZBL