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

#include "meam_implementation.hpp"

#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <map>
#include <numeric>
#include <string>

#include "KIM_LogMacros.hpp"
#include "KIM_ModelDriverHeaders.hpp"
#include "helper.hpp"
#include "meam.hpp"
#include "meam_c.hpp"
#include "meam_spline.hpp"
#include "meam_sw_spline.hpp"
#include "special.hpp"
#include "two_body.hpp"
#include "utility.hpp"

namespace {
/*!
 * \brief The various parameters in the MEAM formulas are listed in three files.
 *
 * 1 - MEAM element names between library & parameter file
 * 2 - MEAM library file in case of `meam/c` or a potential file in case of
 *     `meam/spline` or `meam/sw/spline` potentials.
 * 3 - MEAM parameter file in case of `meam/c` potential.
 *
 */
static constexpr int kMaxNumberParameterFiles = 3;

static constexpr int kForce = 2;
static constexpr int kParticleEnergy = 2;
static constexpr int kVirial = 2;
static constexpr int kParticleVirial = 2;
static constexpr int kKnotsOnRegularGrid = 2;

static constexpr int kMaxLineSize = 1024;                    // 2**10
static constexpr double kOne = 1.0;                          // 1
static constexpr double kThird = 0.33333333333333333333;     // 1/3
static constexpr double kTwoThird = 0.66666666666666666666;  // 2/3

/*! Voight notation index maps for 2D. x-y-z to Voigt-like index */
static constexpr int const kVoight2dIndex[9] = {0, 1, 2, 1, 3, 4, 2, 4, 5};

/*! Voight notation index maps for 3D. x-y-z to Voigt-like index */
static constexpr int const kVoight3dIndex[27] = {0, 1, 2, 1, 3, 4, 2, 4, 5,
                                                 1, 3, 4, 3, 6, 7, 4, 7, 8,
                                                 2, 4, 5, 4, 7, 8, 5, 8, 9};

/*!
 * \brief Array of factors to apply for Voight notation.
 *        multiplicity of Voigt index (i.e. [1] -> xy+yx = 2)
 */
static constexpr int const kVoight2dFactor[6] = {1, 2, 2, 1, 2, 1};

/*!
 * \brief Array of factors to apply for Voight notation.
 *        (multiplicity of Voigt index)
 */
static constexpr int const kVoight3dFactor[10] = {1, 3, 3, 3, 6, 3, 1, 3, 3, 1};
}  // namespace

// public member functions

#ifdef KIM_LOGGER_OBJECT_NAME
#undef KIM_LOGGER_OBJECT_NAME
#endif
#define KIM_LOGGER_OBJECT_NAME model_driver_create
MEAMImplementation::MEAMImplementation(
    KIM::ModelDriverCreate *const model_driver_create,
    KIM::LengthUnit const requested_length_unit,
    KIM::EnergyUnit const requested_energy_unit,
    KIM::ChargeUnit const requested_charge_unit,
    KIM::TemperatureUnit const requested_temperature_unit,
    KIM::TimeUnit const requested_time_unit, int *const ier)
    : meam_c_(nullptr), meam_spline_(nullptr), meam_sw_spline_(nullptr) {
  // Everything is good
  *ier = false;

  if (!model_driver_create) {
    HELPER_LOG_ERROR("The model_driver_create pointer is not assigned.\n");
    *ier = true;
    return;
  }

  {
    std::FILE *parameter_file_pointers[kMaxNumberParameterFiles];

    int number_parameter_files(0);

    // Getting the number of parameter files, open the files and process them
    model_driver_create->GetNumberOfParameterFiles(&number_parameter_files);

    if (number_parameter_files > kMaxNumberParameterFiles) {
      LOG_ERROR("Too many input parameter files!\n");
      *ier = true;
      return;
    }

    if (!number_parameter_files) {
      LOG_ERROR("There is no parameter file!\n");
      *ier = true;
      return;
    }

    *ier = OpenParameterFiles(model_driver_create, number_parameter_files,
                              parameter_file_pointers);
    if (*ier) {
      return;
    }

    *ier = ProcessParameterFiles(model_driver_create, number_parameter_files,
                                 parameter_file_pointers);

    CloseParameterFiles(number_parameter_files, parameter_file_pointers);

    if (*ier) {
      return;
    }
  }

  *ier = ConvertUnits(model_driver_create, requested_length_unit,
                      requested_energy_unit, requested_charge_unit,
                      requested_temperature_unit, requested_time_unit);
  if (*ier) {
    return;
  }

  *ier = SetRefreshMutableValues(model_driver_create);
  if (*ier) {
    return;
  }

  *ier = RegisterKIMModelSettings(model_driver_create);
  if (*ier) {
    return;
  }

  *ier = RegisterKIMParameters(model_driver_create);
  if (*ier) {
    return;
  }

  *ier = RegisterKIMFunctions(model_driver_create);
  if (*ier) {
    return;
  }
}
#undef KIM_LOGGER_OBJECT_NAME

// Private member functions

#define KIM_LOGGER_OBJECT_NAME model_driver_create
int MEAMImplementation::OpenParameterFiles(
    KIM::ModelDriverCreate *const model_driver_create,
    int const number_parameter_files, std::FILE **parameter_file_pointers) {
  std::string const *parameter_file_directory_name;
  model_driver_create->GetParameterFileDirectoryName(
      &parameter_file_directory_name);

  for (int i = 0; i < number_parameter_files; ++i) {
    std::string const *parameter_file_basename;

    if (model_driver_create->GetParameterFileBasename(
            i, &parameter_file_basename)) {
      LOG_ERROR("Unable to get the parameter file base name\n");
      return true;
    }

    std::string const parameter_file_name =
        *parameter_file_directory_name + "/" + *parameter_file_basename;

    parameter_file_pointers[i] = std::fopen(parameter_file_name.c_str(), "r");

    if (!parameter_file_pointers[i]) {
      std::string msg = "The parameter file (" + *parameter_file_basename;
      msg += ") can not be opened.\n";
      HELPER_LOG_ERROR(msg);
      for (int j = i - 1; i <= 0; --i) {
        std::fclose(parameter_file_pointers[j]);
      }
      return true;
    }
  }

  // everything is good
  return false;
}
#undef KIM_LOGGER_OBJECT_NAME

#define KIM_LOGGER_OBJECT_NAME model_driver_create
int MEAMImplementation::ProcessParameterFiles(
    KIM::ModelDriverCreate *const model_driver_create,
    int const number_parameter_files,
    std::FILE *const *parameter_file_pointers) {
  char next_line[kMaxLineSize];

  // Read the MEAM element file

  // create a utility object
  Utility ut;

  // read the first line, strip comment, skip line if blank
  if (ut.GetNextLine(parameter_file_pointers[0], next_line, kMaxLineSize)) {
    HELPER_LOG_ERROR("End of file while reading the MEAM element file.\n");
    return true;
  }

  // MEAM element names

  // element_name_ = list of unique element names
  element_name_ = ut.GetWordsInLine(next_line);

  // number_of_elements_ = # of MEAM element_name_
  number_of_elements_ = static_cast<int>(element_name_.size());

  if (number_of_elements_ < 1) {
    HELPER_LOG_ERROR("Incorrect MEAM element file.\nNo element is provided.\n");
    return true;
  }

  // keep track of known species
  std::map<KIM::SpeciesName const, int, KIM::SPECIES_NAME::Comparator>
      species_map;

  for (int ielem = 0; ielem < number_of_elements_; ++ielem) {
    KIM::SpeciesName const species_name(element_name_[ielem]);

    if (species_name.ToString() == "unknown") {
      std::string msg = "Incorrect format in the MEAM element file.\n";
      msg += "The species name '" + element_name_[ielem] + "' is unknown.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    // check for new species
    auto iter = species_map.find(species_name);
    if (iter == species_map.end()) {
      int ier = model_driver_create->SetSpeciesCode(species_name, ielem);
      if (ier) {
        LOG_ERROR("SetSpeciesCode failed to set the new species.\n");
        return ier;
      }
      species_map[species_name] = ielem;
    } else {
      std::string msg = "Incorrect format in the MEAM element file.\n";
      msg += "The Species '" + element_name_[ielem];
      msg += "' is already defined and exist.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }
  }

  // Detect the MEAM type
  // Any of the `meam/c` or `meam/spline` or `meam/sw/spline`

  do {
    if (number_parameter_files == 3) {
      // only meam/c might have three parameter files
      is_meam_c_ = 1;
      break;
    }

    // Skip the first line. The first line is always a comment
    if (ut.GetFirstLine(parameter_file_pointers[1], next_line, kMaxLineSize)) {
      std::string msg = "End of file while reading the first line ";
      msg += "from the input library/potential file.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    if (ut.GetNextLine(parameter_file_pointers[1], next_line, kMaxLineSize)) {
      std::string msg = "End of file while reading a line ";
      msg += "from the input library/potential file.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    char *word = std::strtok(next_line, " \t\n\r\f");

    if (std::strcmp(word, "meam/spline") == 0) {
      // the first word in the meam/spline new format
      // parameter file is the `meam/spline`
      is_meam_spline_ = 1;
      break;
    }

    // if it is not a meam/spline new format file, it can be:
    //
    //  - an old format meam/spline    : first word is the number of knots (int)
    //  - a             meam/sw/spline : first word is the number of knots (int)
    //  - a             meam/c         : first word is the species name (string)

    int number_knots = std::atoi(word);

    // if no conversion can be performed, ​0​ is returned.
    if (number_knots == 0) {
      is_meam_c_ = 1;
      break;
    }

    // it is not a meam/c, it can be:
    //    - an old format meam/spline    file
    //    - a             meam/sw/spline file

    // both files contain only a single species
    if (number_of_elements_ > 1) {
      std::string msg = "Unable to get the input file format.\n";
      msg += "At this point, input should be an old format ";
      msg += "`meam/spline` file, or a `meam/sw/spline` file. ";
      msg += "Either one supports only a single species.\nThere ";
      msg += "are more species defined in the element file.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    std::rewind(parameter_file_pointers[1]);

    // an old format `meam/spline` file, contains 5 sets of data,
    // while the `meam/sw/spline` file has 7.

    if (ut.GetFirstLine(parameter_file_pointers[1], next_line, kMaxLineSize)) {
      std::string msg = "End of file while reading the first line ";
      msg += "from the input library/potential file.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }

    int nset = 0;

    while (true) {
      if (ut.GetNextLine(parameter_file_pointers[1], next_line, kMaxLineSize)) {
        break;
      }

      word = std::strtok(next_line, " \t\n\r\f");

      number_knots = std::atoi(word);

      if (number_knots < 2) {
        std::string msg = "Invalid number of spline knots (n = ";
        msg += std::to_string(number_knots);
        msg += ") in the MEAM potential file.\n";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      if (ut.GetNextLine(parameter_file_pointers[1], next_line, kMaxLineSize)) {
        std::string msg = "End of file while reading a line from ";
        msg += "a potential file to detect MEAM type.\n";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      if (ut.GetNextLine(parameter_file_pointers[1], next_line, kMaxLineSize)) {
        std::string msg = "End of file while reading a line from ";
        msg += "a potential file to detect MEAM type.\n";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      for (int i = 0; i < number_knots; ++i) {
        if (ut.GetNextLine(parameter_file_pointers[1], next_line,
                           kMaxLineSize)) {
          std::string msg = "End of file while reading a line from ";
          msg += "a potential file to detect MEAM type.\n";
          HELPER_LOG_ERROR(msg);
          return true;
        }
      }

      ++nset;
    }

    // error if the potential file is empty
    if (nset == 0) {
      std::string msg = "The potential file is empty.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    } else if (nset == 5) {
      is_meam_spline_ = 1;
      break;
    } else if (nset == 7) {
      is_meam_sw_spline_ = 1;
      break;
    }

    HELPER_LOG_ERROR("Unable to detect the MEAM type from the input files.\n");
    return true;

  } while (false);

  // Read MEAM library/potential and/or parameter files.

  /*! \note In the current implementation all the input
   *        file comments should be started with \c # */

  // meam/c
  if (is_meam_c_) {
    meam_c_.reset(new MEAMC);

    // Read MEAM library file. @parameter_file_pointers[1]
    if (meam_c_->ProcessLibraryFile(parameter_file_pointers[1], kMaxLineSize,
                                    element_name_)) {
      HELPER_LOG_ERROR("Failed to read and parse the `meam/c` library file.\n");
      return true;
    }

    // done if there is no user paramater file
    if (number_parameter_files == 2) {
      // everything is good
      return false;
    }

    // Read MEAM parameter file. @parameter_file_pointers[2]
    if (meam_c_->ProcessParameterFile(parameter_file_pointers[2],
                                      kMaxLineSize)) {
      std::string msg = "Failed to read and parse the ";
      msg += "`meam/c` parameter file.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }
  }
  // meam/spline
  else if (is_meam_spline_) {
    meam_spline_.reset(new MEAMSpline);

    // Read meam/spline potential file. @parameter_file_pointers[1]
    if (meam_spline_->ProcessPotentialFile(parameter_file_pointers[1],
                                           kMaxLineSize, element_name_)) {
      std::string msg = "Failed to read and parse the ";
      msg += "`meam/spline` potential file.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }
  }
  // meam/sw/spline
  else if (is_meam_sw_spline_) {
    meam_sw_spline_.reset(new MEAMSWSpline);

    // Read meam/spline potential file. @parameter_file_pointers[1]
    if (meam_sw_spline_->ProcessPotentialFile(parameter_file_pointers[1],
                                              kMaxLineSize)) {
      std::string msg = "Failed to read and parse the ";
      msg += "`meam/sw/spline` potential file.\n";
      HELPER_LOG_ERROR(msg);
      return true;
    }
  }

  // everything is good
  return false;
}
#undef KIM_LOGGER_OBJECT_NAME

void MEAMImplementation::CloseParameterFiles(
    int const number_parameter_files,
    std::FILE *const *parameter_file_pointers) {
  for (int i = 0; i < number_parameter_files; ++i) {
    std::fclose(parameter_file_pointers[i]);
  }
}

#define KIM_LOGGER_OBJECT_NAME model_driver_create
int MEAMImplementation::ConvertUnits(
    KIM::ModelDriverCreate *const model_driver_create,
    KIM::LengthUnit const &requested_length_unit,
    KIM::EnergyUnit const &requested_energy_unit,
    KIM::ChargeUnit const &requested_charge_unit,
    KIM::TemperatureUnit const &requested_temperature_unit,
    KIM::TimeUnit const &requested_time_unit) {
  // Define metal unit as a default base units
  KIM::LengthUnit const fromLength = KIM::LENGTH_UNIT::A;
  KIM::EnergyUnit const fromEnergy = KIM::ENERGY_UNIT::eV;
  KIM::ChargeUnit const fromCharge = KIM::CHARGE_UNIT::e;
  KIM::TemperatureUnit const fromTemperature = KIM::TEMPERATURE_UNIT::K;
  KIM::TimeUnit const fromTime = KIM::TIME_UNIT::ps;

  // changing units of length
  double convert_length = kOne;
  int ier = model_driver_create->ConvertUnit(
      fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
      requested_length_unit, requested_energy_unit, requested_charge_unit,
      requested_temperature_unit, requested_time_unit, kOne, 0.0, 0.0, 0.0, 0.0,
      &convert_length);
  if (ier) {
    LOG_ERROR("Unable to convert length unit");
    return ier;
  }

  // Changing units of energy
  double convert_energy = kOne;
  ier = model_driver_create->ConvertUnit(
      fromLength, fromEnergy, fromCharge, fromTemperature, fromTime,
      requested_length_unit, requested_energy_unit, requested_charge_unit,
      requested_temperature_unit, requested_time_unit, 0.0, kOne, 0.0, 0.0, 0.0,
      &convert_energy);
  if (ier) {
    LOG_ERROR("Unable to convert energy unit");
    return ier;
  }

  // Convert to active units
  if (special::IsNotOne(convert_length) || special::IsNotOne(convert_energy)) {
    if (is_meam_c_) {
      meam_c_->ConvertUnit(convert_length, convert_energy);

    } else if (is_meam_spline_) {
      meam_spline_->ConvertUnit(convert_length, convert_energy);

    } else if (is_meam_sw_spline_) {
      meam_sw_spline_->ConvertUnit(convert_length, convert_energy);
    }
  }

  // Register units
  ier = model_driver_create->SetUnits(
      requested_length_unit, requested_energy_unit, KIM::CHARGE_UNIT::unused,
      KIM::TEMPERATURE_UNIT::unused, KIM::TIME_UNIT::unused);
  if (ier) {
    LOG_ERROR("Unable to set units to the requested values");
    return ier;
  }

  // Everything is good
  return false;
}
#undef KIM_LOGGER_OBJECT_NAME

template <class ModelObj>
int MEAMImplementation::SetRefreshMutableValues(ModelObj *const model_obj) {
  if (is_meam_c_) {
    meam_c_->CompleteSetup(&max_cutoff_);

  } else if (is_meam_spline_) {
    if (meam_spline_->CompleteSetup(&max_cutoff_)) {
      HELPER_LOG_ERROR("Failed to complete the setup.\n");
      return true;
    }

  } else if (is_meam_sw_spline_) {
    if (meam_sw_spline_->CompleteSetup(&max_cutoff_)) {
      HELPER_LOG_ERROR("Failed to complete the setup.\n");
      return true;
    }
  }

  influence_distance_ = max_cutoff_;
  max_cutoff_squared_ = max_cutoff_ * max_cutoff_;

  // Update the maximum cutoff or influence distance value in KIM API object
  model_obj->SetInfluenceDistancePointer(&influence_distance_);

  // In this model we do not need to request neighbors from non-contributing
  // particles model_will_not_request_neighbors_of_non_contributing_particles_
  // = 1 Update the cutoff value in KIM API object
  model_obj->SetNeighborListPointers(
      1, &influence_distance_,
      &model_will_not_request_neighbors_of_non_contributing_particles_);

  // Everything is good
  return false;
}

int MEAMImplementation::RegisterKIMModelSettings(
    KIM::ModelDriverCreate *const model_driver_create) const {
  // Register numbering
  return model_driver_create->SetModelNumbering(KIM::NUMBERING::zeroBased);
}

#define KIM_LOGGER_OBJECT_NAME model_driver_create
int MEAMImplementation::RegisterKIMParameters(
    KIM::ModelDriverCreate *const model_driver_create) {
  // Publish parameters
  int ier;

  if (is_meam_c_) {
    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->cutoff_radius_, "rc", "cutoff radius for cutoff function");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'rc'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->delr_, "delr",
        "length of smoothing distance for cutoff function");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'delr'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_, meam_c_->element_rho0_.data(), "rho0",
        "relative density for element I");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'rho0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_, meam_c_->element_Ec_.data(),
        "Ec",
        "cohesive energy of reference structure for I-J mixture.\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the cohesive energy of "
            "reference structure for I-J mixture use:\n"
            "index = ((I - 1) * number_of_elements + J); 'Ec(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Ec'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_,
        meam_c_->element_delta_.data(), "delta",
        "heat of formation for I-J alloy. (if Ec_IJ is input as zero, then "
        "sets Ec_IJ = (Ec_II + Ec_JJ)/2 - delta_IJ.\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the heat of formation for I-J alloy use:\n"
            "index = ((I - 1) * number_of_elements + J); 'delta(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'delta'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_,
        meam_c_->element_alpha_.data(), "alpha",
        "alpha parameter for pair potential between I and J. "
        "(can be computed from bulk modulus of reference structure.)\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the heat alpha parameter between I and J "
            "use:\nindex = ((I - 1) * number_of_elements + J); 'alpha(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'alpha'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_, meam_c_->element_re_.data(),
        "re",
        "equilibrium distance between I and J in the reference structure.\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the equilibrium distance between I and J "
            "use:\nindex = ((I - 1) * number_of_elements + J); 're(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 're'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_ * number_of_elements_,
        meam_c_->element_Cmax_.data(), "Cmax",
        "Cmax screening parameter when I-J pair is screened by K (I<=J).\n"
        "(Tensor of size N x N x N= " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the Cmax screening parameter when I-J pair "
            "screened by K use:\n"
            "index = (((I - 1) * number_of_elements + (J - 1)) * "
            "number_of_elements + k); 'Cmax(I,J,K)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Cmax'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_ * number_of_elements_,
        meam_c_->element_Cmin_.data(), "Cmin",
        "Cmin screening parameter when I-J pair is screened by K (I<=J).\n"
        "(Tensor of size N x N x N= " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the Cmin screening parameter when I-J pair "
            "screened by K use:\n"
            "index = (((I - 1) * number_of_elements + (J - 1)) * "
            "number_of_elements + k); 'Cmin(I,J,K)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Cmin'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_,
        reinterpret_cast<int *>(meam_c_->element_lattice_.data()), "lattce",
        "lattice structure of I-J reference structure.\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the lattice structure of I-J reference "
            "structure use:\nindex = ((I - 1) * number_of_elements + J); "
            "'lattce(I,J)'\n"
            " 0; fcc = face centered cubic\n"
            " 1; bcc = body centered cubic\n"
            " 2; hcp = hexagonal close-packed\n"
            " 3; dim = dimer\n"
            " 4; dia = diamond (interlaced fcc for alloy)\n"
            " 5; dia3= diamond structure with primary 1NN and secondary 3NN "
            "interaction\n"
            " 6; b1  = rock salt (NaCl structure)\n"
            " 7; c11 = MoSi2 structure\n"
            " 8; l12 = Cu3Au structure (lower case L, followed by 12)\n"
            " 9; b2  = CsCl structure (interpenetrating simple cubic)\n"
            "10; ch4 = methane-like structure, only for binary system\n"
            "11; lin = linear structure (180 degree angle)\n"
            "12; zig = zigzag structure with a uniform angle\n"
            "13; tri = H2O-like structure that has an angle \n");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'lattce'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_, meam_c_->element_nn2_.data(),
        "nn2",
        "turn on second-nearest neighbor MEAM formulation for I-J pair.\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find for I-J pair use:\n"
            "index = ((I - 1) * number_of_elements + J); 'nn2(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'nn2'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_,
        meam_c_->element_attrac_.data(), "attrac",
        "additional cubic attraction term in Rose energy I-J pair potential.\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the attraction term between I and J use:\n"
            "index = ((I - 1) * number_of_elements + J); 'attrac(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'attrac'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_,
        meam_c_->element_repuls_.data(), "repuls",
        "additional cubic repulsive term in Rose energy I-J pair potential.\n"
        "(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the repulsive term between I and J use:\n"
            "index = ((I - 1) * number_of_elements + J); 'repuls(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'repuls'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_, meam_c_->element_zbl_.data(),
        "zbl",
        "blend the MEAM I-J pair potential with the ZBL potential for small "
        "atom separations.\n(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the term between I and J use:\n"
            "index = ((I - 1) * number_of_elements + J); 'zbl(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'zbl'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        number_of_elements_ * number_of_elements_,
        meam_c_->element_theta_.data(), "theta",
        "angle between three atoms in line, zigzag, and trimer reference "
        "structures in degrees.\n(matrix of size N x N = " +
            std::to_string(number_of_elements_) + " x " +
            std::to_string(number_of_elements_) + ") " +
            "in row-major storage.\n"
            "For example, to find the term between I and J use:\n"
            "index = ((I - 1) * number_of_elements + J); 'theta(I,J)'");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'theta'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->gsmooth_factor_, "gsmooth_factor",
        "factor determining the length of the G-function smoothing region; "
        "only significant for ibar=0 or ibar=4.\n"
        "99.0 = short smoothing region, sharp step\n"
        "0.5  = long smoothing region, smooth step");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'gsmooth_factor'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->augment_t1_, "augt1",
        "integer flag for whether to augment t1 parameter by 3/5*t3 to account "
        "for old vs new meam formulations;\n0 = don't augment t1\n"
        "1 = augment t1");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'augt1'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->ialloy_, "ialloy",
        "integer flag to use alternative averaging rule for t parameters, "
        "for comparison with the DYNAMO MEAM code;\n"
        "0 = standard averaging (matches ialloy=0 in DYNAMO)\n"
        "1 = alternative averaging (matches ialloy=1 in DYNAMO)\n"
        "2 = no averaging of t (use single-element values)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'ialloy'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->mixing_rule_compute_t_, "mixture_ref_t",
        "integer flag to use mixture average of t to compute the background "
        "reference density for alloys, instead of the single-element values\n"
        "0 = do not use mixture averaging for t in the reference density\n"
        "1 = use mixture averaging for t in the reference density");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'mixture_ref_t'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->rose_function_form_, "erose_form",
        "integer value to select the form of the Rose energy function\n"
        "0: erose = -Ec*(1+astar+a3*(astar**3)/(r/re))*exp(-astar)\n"
        "1: erose = -Ec*(1+astar+(-attrac+repuls/r)*(astar**3))*exp(-astar)\n"
        "2: erose = -Ec*(1 +astar + a3*(astar**3))*exp(-astar)");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'erose_form'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->fcut_function_form_, "fcut_form",
        "integer value to select the form of the MEAM cut-off function\n"
        "0: has the (1-x) term to the fourth power\n"
        "1: has the (1-x) term to the sixth power");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'fcut_form'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->use_rhobar_linear_embedding_function_, "emb_lin_neg",
        "integer value to select embedding function for negative densities\n"
        "0 = F(rho)=0\n"
        "1 = F(rho) = -asub*esub*rho (linear in rho, matches DYNAMO)");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'emb_lin_neg'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, &meam_c_->dynamo_reference_density_, "bkgd_dyn",
        "integer value to select background density formula\n"
        "0 = rho_bkgd = rho_ref_meam(a) (as in the reference structure)\n"
        "1 = rho_bkgd = rho0_meam(a)*Z_meam(a) (matches DYNAMO)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'bkgd_dyn'");
      return ier;
    }
  }  // is_meam_c_
  else if (is_meam_spline_) {
    for (int i = 0; i < number_of_elements_; ++i) {
      int const _ind = i * number_of_elements_ - i * (i + 1) / 2;
      for (int j = i; j < number_of_elements_; ++j) {
        int const ind = _ind + j;

        std::string name =
            "Phi_x_" + ((i == j) ? element_name_[i]
                                 : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            meam_spline_->e_phi_[ind].GetNumberOfKnots(),
            meam_spline_->e_phi_[ind].x_data(), name,
            "Positions of spline knots");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }

        name = "Phi_" + ((i == j) ? element_name_[i]
                                  : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            meam_spline_->e_phi_[ind].GetNumberOfKnots(),
            meam_spline_->e_phi_[ind].value_data(), name,
            "The pair energy function tabulated over values of atomic "
            "separation distance");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }

        name =
            "Phi_d0_" + ((i == j) ? element_name_[i]
                                  : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            1, meam_spline_->e_phi_[ind].derivative_0_data(), name,
            "The first derivative of the pair energy function at the beginning "
            "of the spline");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }

        name =
            "Phi_dN_" + ((i == j) ? element_name_[i]
                                  : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            1, meam_spline_->e_phi_[ind].derivative_n_data(), name,
            "The first derivative of the pair energy function at the end of "
            "the spline");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }
      }  // Loop over j
    }    // Loop over i

    for (int i = 0; i < number_of_elements_; ++i) {
      std::string name = "U_x_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          meam_spline_->e_u_[i].GetNumberOfKnots(),
          meam_spline_->e_u_[i].x_data(), name, "Positions of spline knots");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "U_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          meam_spline_->e_u_[i].GetNumberOfKnots(),
          meam_spline_->e_u_[i].value_data(), name,
          "The embedding energy functional tabulated over values of the MEAM "
          "density");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "U_d0_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          1, meam_spline_->e_u_[i].derivative_0_data(), name,
          "The first derivative of the embedding energy functional at the "
          "beginning of the spline");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "U_dN_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          1, meam_spline_->e_u_[i].derivative_n_data(), name,
          "The first derivative of the embedding energy functional at the end "
          "of the spline");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }
    }  // Loop over i

    for (int i = 0; i < number_of_elements_; ++i) {
      std::string name = "rho_x_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          meam_spline_->rho_r_[i].GetNumberOfKnots(),
          meam_spline_->rho_r_[i].x_data(), name, "Positions of spline knots");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "rho_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          meam_spline_->rho_r_[i].GetNumberOfKnots(),
          meam_spline_->rho_r_[i].value_data(), name,
          "The rho function in computing the MEAM density tabulated over "
          "values of the atomic separation distance");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "rho_d0_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          1, meam_spline_->rho_r_[i].derivative_0_data(), name,
          "The first derivative of the rho function in computing the MEAM "
          "density at the beginning of the spline");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "rho_dN_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          1, meam_spline_->rho_r_[i].derivative_n_data(), name,
          "The first derivative of the rho function in computing the MEAM "
          "density at the end of the spline");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }
    }  // Loop over i

    for (int i = 0; i < number_of_elements_; ++i) {
      std::string name = "f_x_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          meam_spline_->rho_f_[i].GetNumberOfKnots(),
          meam_spline_->rho_f_[i].x_data(), name, "Positions of spline knots");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "f_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          meam_spline_->rho_f_[i].GetNumberOfKnots(),
          meam_spline_->rho_f_[i].value_data(), name,
          "The f function in computing the MEAM density tabulated over values "
          "of the atomic separation distance. It comprises the three-body "
          "contributions to the density");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "f_d0_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          1, meam_spline_->rho_f_[i].derivative_0_data(), name,
          "The first derivative of the f function in computing the MEAM "
          "density at the beginning of the spline");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }

      name = "f_dN_" + element_name_[i];
      ier = model_driver_create->SetParameterPointer(
          1, meam_spline_->rho_f_[i].derivative_n_data(), name,
          "The first derivative of the f function in computing the MEAM "
          "density at the end of the spline");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for '" + name +
                  "'");
        return ier;
      }
    }  // Loop over i

    for (int i = 0; i < number_of_elements_; ++i) {
      int const _ind = i * number_of_elements_ - i * (i + 1) / 2;
      for (int j = i; j < number_of_elements_; ++j) {
        int const ind = _ind + j;

        std::string name =
            "g_x_" + ((i == j) ? element_name_[i]
                               : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            meam_spline_->rho_g_[ind].GetNumberOfKnots(),
            meam_spline_->rho_g_[ind].x_data(), name,
            "Positions of spline knots");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }

        name = "g_" + ((i == j) ? element_name_[i]
                                : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            meam_spline_->rho_g_[ind].GetNumberOfKnots(),
            meam_spline_->rho_g_[ind].value_data(), name,
            "The g function in computing the MEAM density tabulated over "
            "values of the bond angle between atoms j, i, and k, centered "
            "on i. It comprises the three-body contributions to the density");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }

        name = "g_d0_" + ((i == j) ? element_name_[i]
                                   : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            1, meam_spline_->rho_g_[ind].derivative_0_data(), name,
            "The first derivative of the g function in computing the MEAM "
            "density at the beginning of the spline");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }

        name = "g_dN_" + ((i == j) ? element_name_[i]
                                   : element_name_[i] + "_" + element_name_[j]);
        ier = model_driver_create->SetParameterPointer(
            1, meam_spline_->rho_g_[ind].derivative_n_data(), name,
            "The first derivative of the g function in computing the MEAM "
            "density at the end of the spline");
        if (ier) {
          LOG_ERROR("SetParameterPointer failed to set the pointer for '" +
                    name + "'");
          return ier;
        }
      }  // Loop over j
    }    // Loop over i
  }      // is_meam_spline_
  else if (is_meam_sw_spline_) {
    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->e_phi_.GetNumberOfKnots(),
        meam_sw_spline_->e_phi_.x_data(), "Phi_x", "Positions of spline knots");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Phi'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->e_phi_.GetNumberOfKnots(),
        meam_sw_spline_->e_phi_.value_data(), "Phi",
        "The pair energy function tabulated over values of atomic "
        "separation distance");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Phi'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->e_phi_.derivative_0_data(), "Phi_d0",
        "The first derivative of the pair energy function at the beginning of "
        "the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Phi_d0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->e_phi_.derivative_n_data(), "Phi_dN",
        "The first derivative of the pair energy function at the end of the "
        "spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Phi_dN'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->e_u_.GetNumberOfKnots(),
        meam_sw_spline_->e_u_.x_data(), "U_x", "Positions of spline knots");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'U'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->e_u_.GetNumberOfKnots(),
        meam_sw_spline_->e_u_.value_data(), "U",
        "The embedding energy functional tabulated over values of the MEAM "
        "density");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'U'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->e_u_.derivative_0_data(), "U_d0",
        "The first derivative of the embedding energy functional at the "
        "beginning of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'U_d0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->e_u_.derivative_n_data(), "U_dN",
        "The first derivative of the embedding energy functional at the end of "
        "the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'U_dN'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->rho_r_.GetNumberOfKnots(),
        meam_sw_spline_->rho_r_.x_data(), "rho_x", "Positions of spline knots");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'rho'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->rho_r_.GetNumberOfKnots(),
        meam_sw_spline_->rho_r_.value_data(), "rho",
        "The rho function in computing the MEAM density tabulated over "
        "values of the atomic separation distance");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'rho'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->rho_r_.derivative_0_data(), "rho_d0",
        "The first derivative of the rho function in computing the MEAM "
        "density at the beginning of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'rho_d0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->rho_r_.derivative_n_data(), "rho_dN",
        "The first derivative of the rho function in computing the MEAM "
        "density at the end of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'rho_dN'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->rho_f_.GetNumberOfKnots(),
        meam_sw_spline_->rho_f_.x_data(), "f_x", "Positions of spline knots");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'f'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->rho_f_.GetNumberOfKnots(),
        meam_sw_spline_->rho_f_.value_data(), "f",
        "The f function in computing the MEAM density tabulated over values "
        "of the atomic separation distance. It comprises the three-body "
        "contributions to the density");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'f'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->rho_f_.derivative_0_data(), "f_d0",
        "The first derivative of the f function in computing the MEAM "
        "density at the beginning of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'f_d0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->rho_f_.derivative_n_data(), "f_dN",
        "The first derivative of the f function in computing the MEAM "
        "density at the end of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'f_dN'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->rho_g_.GetNumberOfKnots(),
        meam_sw_spline_->rho_g_.x_data(), "g_x", "Positions of spline knots");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'g'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->rho_g_.GetNumberOfKnots(),
        meam_sw_spline_->rho_g_.value_data(), "g",
        "The g function in computing the MEAM density tabulated over "
        "values of the bond angle between atoms j, i, and k, centered "
        "on i. It comprises the three-body contributions to the density");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'g'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->rho_g_.derivative_0_data(), "g_d0",
        "The first derivative of the g function in computing the MEAM "
        "density at the beginning of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'g_d0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->rho_g_.derivative_n_data(), "g_dN",
        "The first derivative of the g function in computing the MEAM "
        "density at the end of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'g_dN'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->esw_f_.GetNumberOfKnots(),
        meam_sw_spline_->esw_f_.x_data(), "F_SW_x",
        "Positions of spline knots");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'F_SW'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->esw_f_.GetNumberOfKnots(),
        meam_sw_spline_->esw_f_.value_data(), "F_SW",
        "The F_SW function in computing the MEAM energy tabulated over values "
        "of the atomic separation distance. It comprises the SW three-body "
        "contributions to the energy");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'F_SW'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->esw_f_.derivative_0_data(), "F_SW_d0",
        "The first derivative of the F_SW function in computing the MEAM "
        "energy at the beginning of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'F_SW_d0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->esw_f_.derivative_n_data(), "F_SW_dN",
        "The first derivative of the F_SW function in computing the MEAM "
        "energy at the end of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'F_SW_dN'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->esw_g_.GetNumberOfKnots(),
        meam_sw_spline_->esw_g_.x_data(), "G_SW_x",
        "Positions of spline knots");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'G_SW'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        meam_sw_spline_->esw_g_.GetNumberOfKnots(),
        meam_sw_spline_->esw_g_.value_data(), "G_SW",
        "The G_SW function in computing the MEAM energy tabulated over "
        "values of the bond angle between atoms j, i, and k, centered "
        "on i. It comprises the three-body contributions to the energy");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'G_SW'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->esw_g_.derivative_0_data(), "G_SW_d0",
        "The first derivative of the G_SW function in computing the MEAM "
        "energy at the beginning of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'G_SW_d0'");
      return ier;
    }

    ier = model_driver_create->SetParameterPointer(
        1, meam_sw_spline_->esw_g_.derivative_n_data(), "G_SW_dN",
        "The first derivative of the G_SW function in computing the MEAM "
        "energy at the end of the spline");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'G_SW_dN'");
      return ier;
    }
  }  // is_meam_sw_spline_

  // Everything is good
  return false;
}
#undef KIM_LOGGER_OBJECT_NAME

int MEAMImplementation::RegisterKIMFunctions(
    KIM::ModelDriverCreate *const model_driver_create) const {
  // Register functions
  int ier =
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::Destroy, KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(MEAM::Destroy)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::Refresh, KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(MEAM::Refresh)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::WriteParameterizedModel,
          KIM::LANGUAGE_NAME::cpp, false,
          reinterpret_cast<KIM::Function *>(MEAM::WriteParameterizedModel)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::Compute, KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(MEAM::Compute)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::ComputeArgumentsCreate,
          KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(MEAM::ComputeArgumentsCreate)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::ComputeArgumentsDestroy,
          KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(MEAM::ComputeArgumentsDestroy));
  return ier;
}

#define KIM_LOGGER_OBJECT_NAME model_compute_arguments
int MEAMImplementation::SetComputeMutableValues(
    KIM::ModelComputeArguments const *const model_compute_arguments,
    bool &is_compute_energy, bool &is_compute_forces,
    bool &is_compute_particle_energy, bool &is_compute_virial,
    bool &is_compute_particle_virial, int const *&particle_species_codes,
    int const *&particle_contributing, VectorOfSizeDIM const *&coordinates,
    double *&energy, VectorOfSizeDIM *&forces, double *&particle_energy,
    VectorOfSizeSix *&virial, VectorOfSizeSix *&particle_virial) {
  // get compute flags
  int const *number_of_particles = nullptr;
  int ier = model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::numberOfParticles,
                &number_of_particles) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::particleSpeciesCodes,
                &particle_species_codes) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::particleContributing,
                &particle_contributing) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::coordinates,
                (double const **)&coordinates) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::partialEnergy, &energy) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::partialForces,
                (double const **)&forces) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::partialParticleEnergy,
                &particle_energy) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::partialVirial,
                (double const **)&virial) ||
            model_compute_arguments->GetArgumentPointer(
                KIM::COMPUTE_ARGUMENT_NAME::partialParticleVirial,
                (double const **)&particle_virial);
  if (ier) {
    LOG_ERROR("GetArgumentPointer return an error");
    return true;
  }

  is_compute_energy = (energy);
  is_compute_forces = (forces);
  is_compute_particle_energy = (particle_energy);
  is_compute_virial = (virial);
  is_compute_particle_virial = (particle_virial);

  // Update value
  cached_number_of_particles_ = *number_of_particles;

  // everything is good
  return false;
}
#undef KIM_LOGGER_OBJECT_NAME

int MEAMImplementation::GetComputeIndex(
    bool const is_compute_energy, bool const is_compute_forces,
    bool const is_compute_particle_energy, bool const is_compute_virial,
    bool const is_compute_particle_virial) const {
  int const is_compute_energy_int = static_cast<int>(is_compute_energy);
  int const is_compute_forces_int = static_cast<int>(is_compute_forces);
  int const is_compute_particle_energy_int =
      static_cast<int>(is_compute_particle_energy);
  int const is_compute_virial_int = static_cast<int>(is_compute_virial);
  int const is_compute_particle_virial_int =
      static_cast<int>(is_compute_particle_virial);

  // energy
  int index = is_compute_energy_int * kForce * kParticleEnergy * kVirial *
              kParticleVirial;
  // force
  index += is_compute_forces_int * kParticleEnergy * kVirial * kParticleVirial;
  // particle_energy
  index += is_compute_particle_energy_int * kVirial * kParticleVirial;
  // virial
  index += is_compute_virial_int * kParticleVirial;
  // particle_virial
  index += is_compute_particle_virial_int;
  return index;
}

int MEAMImplementation::GetComputeIndex(
    bool const is_compute_energy, bool const is_compute_forces,
    bool const is_compute_particle_energy, bool const is_compute_virial,
    bool const is_compute_particle_virial,
    bool const are_knots_on_regular_grid) const {
  int const is_compute_energy_int = static_cast<int>(is_compute_energy);
  int const is_compute_forces_int = static_cast<int>(is_compute_forces);
  int const is_compute_particle_energy_int =
      static_cast<int>(is_compute_particle_energy);
  int const is_compute_virial_int = static_cast<int>(is_compute_virial);
  int const is_compute_particle_virial_int =
      static_cast<int>(is_compute_particle_virial);
  int const are_knots_on_regular_grid_int =
      static_cast<int>(are_knots_on_regular_grid);

  // energy
  int index = is_compute_energy_int * kForce * kParticleEnergy * kVirial *
              kParticleVirial * kKnotsOnRegularGrid;
  // force
  index += is_compute_forces_int * kParticleEnergy * kVirial * kParticleVirial *
           kKnotsOnRegularGrid;
  // particle_energy
  index += is_compute_particle_energy_int * kVirial * kParticleVirial *
           kKnotsOnRegularGrid;
  // virial
  index += is_compute_virial_int * kParticleVirial * kKnotsOnRegularGrid;
  // particle_virial
  index += is_compute_particle_virial_int * kKnotsOnRegularGrid;
  // knots_are_on_regular_grid
  index += are_knots_on_regular_grid_int;
  return index;
}

#define KIM_LOGGER_OBJECT_NAME model_compute_arguments_create
int MEAMImplementation::RegisterKIMComputeArgumentsSettings(
    KIM::ModelComputeArgumentsCreate *const model_compute_arguments_create)
    const {
  // register arguments
  LOG_INFORMATION("Register argument support status");

  int ier = model_compute_arguments_create->SetArgumentSupportStatus(
                KIM::COMPUTE_ARGUMENT_NAME::partialEnergy,
                KIM::SUPPORT_STATUS::optional) ||
            model_compute_arguments_create->SetArgumentSupportStatus(
                KIM::COMPUTE_ARGUMENT_NAME::partialForces,
                KIM::SUPPORT_STATUS::optional) ||
            model_compute_arguments_create->SetArgumentSupportStatus(
                KIM::COMPUTE_ARGUMENT_NAME::partialParticleEnergy,
                KIM::SUPPORT_STATUS::optional) ||
            model_compute_arguments_create->SetArgumentSupportStatus(
                KIM::COMPUTE_ARGUMENT_NAME::partialVirial,
                KIM::SUPPORT_STATUS::optional) ||
            model_compute_arguments_create->SetArgumentSupportStatus(
                KIM::COMPUTE_ARGUMENT_NAME::partialParticleVirial,
                KIM::SUPPORT_STATUS::optional);
  return ier;
}
#undef KIM_LOGGER_OBJECT_NAME

// public method

int MEAMImplementation::Refresh(KIM::ModelRefresh *const model_refresh) {
  return SetRefreshMutableValues(model_refresh);
}

int MEAMImplementation::WriteParameterizedModel(
    KIM::ModelWriteParameterizedModel const
        *const model_write_parameterized_model) const {
  if (!model_write_parameterized_model) {
    std::string msg = "The model_write_parameterized_model ";
    msg += "pointer is not assigned.\n";
    HELPER_LOG_ERROR(msg);
    return true;
  }

  std::string const *path = nullptr;
  std::string const *model_name = nullptr;

  model_write_parameterized_model->GetPath(&path);
  model_write_parameterized_model->GetModelName(&model_name);

  // MEAM elements file
  {
    std::string elements_file(*model_name);

    elements_file += ".elements.meam";

    // Set the file name for the next parameter file.
    model_write_parameterized_model->SetParameterFileName(elements_file);

    // MEAM elements file
    std::string buffer = *path + "/" + elements_file;

    std::fstream fs;

    fs.open(buffer.c_str(), std::fstream::out);

    if (!fs.is_open()) {
      HELPER_LOG_ERROR("Unable to open the elements file for writing.\n");
      return true;
    } else {
      if (fs.fail()) {
        std::string msg = "An error has occurred from opening the ";
        msg += "file on the associated stream!\n";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      fs << "# MEAM elements ";
      for (int ielem = 0; ielem < number_of_elements_; ++ielem) {
        fs << element_name_[ielem] << " ";
      }

      fs << "\n\n\n\n#\n# MEAM elements file\n#\n";
      fs.close();
    }
  }

  if (is_meam_c_) {
    std::string parameter_file;
    for (auto i : element_name_) {
      parameter_file += i;
    }
    parameter_file += ".parameter.meam";

    // Set the file name for the next parameter file.
    model_write_parameterized_model->SetParameterFileName(parameter_file);

    // MEAM elements file
    std::string buffer = *path + "/" + parameter_file;

    std::fstream fs;

    fs.open(buffer.c_str(), std::fstream::out);

    if (!fs.is_open()) {
      HELPER_LOG_ERROR("Unable to open the parameter file for writing.\n");
      return true;
    } else {
      if (fs.fail()) {
        std::string msg = "An error has occurred from opening the ";
        msg += "file on the associated stream!\n";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      fs << "# MEAM parameter file ";
      if (!meam_c_->augment_t1_) {
        fs << "augt1 = 0";
      }
      if (meam_c_->ialloy_) {
        fs << "ialloy = " << meam_c_->ialloy_;
      }
      if (meam_c_->mixing_rule_compute_t_) {
        fs << "mixture_ref_t = 1";
      }
      if (meam_c_->rose_function_form_) {
        fs << "erose_form = " << meam_c_->rose_function_form_;
      }
      if (meam_c_->fcut_function_form_) {
        fs << "fcut_form = " << meam_c_->fcut_function_form_;
      }
      if (meam_c_->use_rhobar_linear_embedding_function_) {
        fs << "emb_lin_neg = "
           << meam_c_->use_rhobar_linear_embedding_function_;
      }
      if (meam_c_->dynamo_reference_density_) {
        fs << "bkgd_dyn = " << meam_c_->dynamo_reference_density_;
      }
      if (special::IsNotZero(meam_c_->cutoff_radius_ - 4.0)) {
        fs << "rc = " << meam_c_->cutoff_radius_;
      }
      if (special::IsNotZero(meam_c_->delr_ - 0.1)) {
        fs << "delr = " << meam_c_->delr_;
      }
      if (special::IsNotZero(meam_c_->gsmooth_factor_ - 99.0)) {
        fs << "gsmooth_factor = " << meam_c_->gsmooth_factor_;
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        fs << "rho0(" << i + 1 << ") = " << meam_c_->element_rho0_[i];
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (meam_c_->element_lattice_(i, j) != Lattice::FCC) {
            fs << "lattce(" << i + 1 << "," << j + 1 << ") = '"
               << meam_c_->LatticeToString(meam_c_->element_lattice_(i, j))
               << "'";
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (meam_c_->element_nn2_(i, j)) {
            fs << "nn2(" << i + 1 << "," << j + 1 << ") = 1";
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (!meam_c_->element_zbl_(i, j)) {
            fs << "zbl(" << i + 1 << "," << j + 1 << ") = 0";
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (special::IsNotZero(meam_c_->element_alpha_(i, j))) {
            fs << "alpha(" << i + 1 << "," << j + 1
               << ") = " << meam_c_->element_alpha_(i, j);
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (special::IsNotZero(meam_c_->element_re_(i, j))) {
            fs << "re(" << i + 1 << "," << j + 1
               << ") = " << meam_c_->element_re_(i, j);
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (special::IsNotZero(meam_c_->element_Ec_(i, j))) {
            fs << "Ec(" << i + 1 << "," << j + 1
               << ") = " << meam_c_->element_Ec_(i, j);
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (special::IsNotZero(meam_c_->element_delta_(i, j))) {
            fs << "delta(" << i + 1 << "," << j + 1
               << ") = " << meam_c_->element_delta_(i, j);
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (special::IsNotZero(meam_c_->element_attrac_(i, j))) {
            fs << "attrac(" << i + 1 << "," << j + 1
               << ") = " << meam_c_->element_attrac_(i, j);
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (special::IsNotZero(meam_c_->element_repuls_(i, j))) {
            fs << "repuls(" << i + 1 << "," << j + 1
               << ") = " << meam_c_->element_repuls_(i, j);
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          if (special::IsNotZero(meam_c_->element_theta_(i, j) - 180.0)) {
            fs << "theta(" << i + 1 << "," << j + 1
               << ") = " << meam_c_->element_theta_(i, j);
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          for (int k = 0; k < number_of_elements_; ++k) {
            if (special::IsNotZero(meam_c_->element_Cmin_(i, j, k) - 2.0)) {
              fs << "Cmin(" << i + 1 << "," << j + 1 << "," << k + 1
                 << ") = " << meam_c_->element_Cmin_(i, j, k);
            }
          }
        }
      }
      for (int i = 0; i < number_of_elements_; ++i) {
        for (int j = i; j < number_of_elements_; ++j) {
          for (int k = 0; k < number_of_elements_; ++k) {
            if (special::IsNotZero(meam_c_->element_Cmax_(i, j, k) - 2.8)) {
              fs << "Cmax(" << i + 1 << "," << j + 1 << "," << k + 1
                 << ") = " << meam_c_->element_Cmax_(i, j, k);
            }
          }
        }
      }

      fs << "\n\n\n\n#\n# MEAM parameter file\n#\n";
      fs.close();
    }

  } else if (is_meam_spline_) {
    std::string potential_file;
    for (auto i : element_name_) {
      potential_file += i;
    }
    potential_file += ".potential.meam.spline";

    // Set the file name for the next parameter file.
    model_write_parameterized_model->SetParameterFileName(potential_file);

    // MEAM elements file
    std::string buffer = *path + "/" + potential_file;

    std::fstream fs;

    fs.open(buffer.c_str(), std::fstream::out);

    if (!fs.is_open()) {
      HELPER_LOG_ERROR("Unable to open the potential file for writing.\n");
      return true;
    } else {
      if (fs.fail()) {
        std::string msg = "An error has occurred from opening the ";
        msg += "file on the associated stream!\n";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      fs << "# MEAM SPLINE potential ";
      {
        std::string unique_names;
        for (auto i : element_name_) {
          unique_names += " " + i;
        }
        fs << "meam/spline " << number_of_elements_ << unique_names;
      }
      for (std::size_t i = 0; i < meam_spline_->e_phi_.size(); ++i) {
        fs << "spline3eq";
        int const number_of_knots = meam_spline_->e_phi_[i].GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_spline_->e_phi_[i].derivative_0_data() << " "
           << *meam_spline_->e_phi_[i].derivative_n_data();
        auto x = meam_spline_->e_phi_[i].x_data();
        auto v = meam_spline_->e_phi_[i].value_data();
        auto d = meam_spline_->e_phi_[i].second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      for (std::size_t i = 0; i < meam_spline_->rho_r_.size(); ++i) {
        fs << "spline3eq";
        int const number_of_knots = meam_spline_->rho_r_[i].GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_spline_->rho_r_[i].derivative_0_data() << " "
           << *meam_spline_->rho_r_[i].derivative_n_data();
        auto x = meam_spline_->rho_r_[i].x_data();
        auto v = meam_spline_->rho_r_[i].value_data();
        auto d = meam_spline_->rho_r_[i].second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      for (std::size_t i = 0; i < meam_spline_->e_u_.size(); ++i) {
        fs << "spline3eq";
        int const number_of_knots = meam_spline_->e_u_[i].GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_spline_->e_u_[i].derivative_0_data() << " "
           << *meam_spline_->e_u_[i].derivative_n_data();
        auto x = meam_spline_->e_u_[i].x_data();
        auto v = meam_spline_->e_u_[i].value_data();
        auto d = meam_spline_->e_u_[i].second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      for (std::size_t i = 0; i < meam_spline_->rho_f_.size(); ++i) {
        fs << "spline3eq";
        int const number_of_knots = meam_spline_->rho_f_[i].GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_spline_->rho_f_[i].derivative_0_data() << " "
           << *meam_spline_->rho_f_[i].derivative_n_data();
        auto x = meam_spline_->rho_f_[i].x_data();
        auto v = meam_spline_->rho_f_[i].value_data();
        auto d = meam_spline_->rho_f_[i].second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      for (std::size_t i = 0; i < meam_spline_->rho_g_.size(); ++i) {
        fs << "spline3eq";
        int const number_of_knots = meam_spline_->rho_g_[i].GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_spline_->rho_g_[i].derivative_0_data() << " "
           << *meam_spline_->rho_g_[i].derivative_n_data();
        auto x = meam_spline_->rho_g_[i].x_data();
        auto v = meam_spline_->rho_g_[i].value_data();
        auto d = meam_spline_->rho_g_[i].second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }

      fs << "\n\n\n\n#\n# MEAM SPLINE potential file\n#\n";
      fs.close();
    }

  } else if (is_meam_sw_spline_) {
    std::string potential_file;
    for (auto i : element_name_) {
      potential_file += i;
    }
    potential_file += ".potential.meam.sw.spline";

    // Set the file name for the next parameter file.
    model_write_parameterized_model->SetParameterFileName(potential_file);

    // MEAM elements file
    std::string buffer = *path + "/" + potential_file;

    std::fstream fs;

    fs.open(buffer.c_str(), std::fstream::out);

    if (!fs.is_open()) {
      HELPER_LOG_ERROR("Unable to open the potential file for writing.\n");
      return true;
    } else {
      if (fs.fail()) {
        std::string msg = "An error has occurred from opening the ";
        msg += "file on the associated stream!\n";
        HELPER_LOG_ERROR(msg);
        return true;
      }

      fs << "# MEAM SW SPLINE potential ";

      // e_phi_
      {
        int const number_of_knots = meam_sw_spline_->e_phi_.GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_sw_spline_->e_phi_.derivative_0_data() << " "
           << *meam_sw_spline_->e_phi_.derivative_n_data();
        fs << meam_sw_spline_->e_phi_.OldFormatString();
        auto x = meam_sw_spline_->e_phi_.x_data();
        auto v = meam_sw_spline_->e_phi_.value_data();
        auto d = meam_sw_spline_->e_phi_.second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      // esw_f_
      {
        int const number_of_knots = meam_sw_spline_->esw_f_.GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_sw_spline_->esw_f_.derivative_0_data() << " "
           << *meam_sw_spline_->esw_f_.derivative_n_data();
        fs << meam_sw_spline_->esw_f_.OldFormatString();
        auto x = meam_sw_spline_->esw_f_.x_data();
        auto v = meam_sw_spline_->esw_f_.value_data();
        auto d = meam_sw_spline_->esw_f_.second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      // esw_g_
      {
        int const number_of_knots = meam_sw_spline_->esw_g_.GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_sw_spline_->esw_g_.derivative_0_data() << " "
           << *meam_sw_spline_->esw_g_.derivative_n_data();
        fs << meam_sw_spline_->esw_g_.OldFormatString();
        auto x = meam_sw_spline_->esw_g_.x_data();
        auto v = meam_sw_spline_->esw_g_.value_data();
        auto d = meam_sw_spline_->esw_g_.second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      // rho_r_
      {
        int const number_of_knots = meam_sw_spline_->rho_r_.GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_sw_spline_->rho_r_.derivative_0_data() << " "
           << *meam_sw_spline_->rho_r_.derivative_n_data();
        fs << meam_sw_spline_->rho_r_.OldFormatString();
        auto x = meam_sw_spline_->rho_r_.x_data();
        auto v = meam_sw_spline_->rho_r_.value_data();
        auto d = meam_sw_spline_->rho_r_.second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      // e_u_
      {
        int const number_of_knots = meam_sw_spline_->e_u_.GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_sw_spline_->e_u_.derivative_0_data() << " "
           << *meam_sw_spline_->e_u_.derivative_n_data();
        fs << meam_sw_spline_->e_u_.OldFormatString();
        auto x = meam_sw_spline_->e_u_.x_data();
        auto v = meam_sw_spline_->e_u_.value_data();
        auto d = meam_sw_spline_->e_u_.second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      // rho_f_
      {
        int const number_of_knots = meam_sw_spline_->rho_f_.GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_sw_spline_->rho_f_.derivative_0_data() << " "
           << *meam_sw_spline_->rho_f_.derivative_n_data();
        fs << meam_sw_spline_->rho_f_.OldFormatString();
        auto x = meam_sw_spline_->rho_f_.x_data();
        auto v = meam_sw_spline_->rho_f_.value_data();
        auto d = meam_sw_spline_->rho_f_.second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }
      // rho_g_
      {
        int const number_of_knots = meam_sw_spline_->rho_g_.GetNumberOfKnots();
        fs << number_of_knots;
        fs << *meam_sw_spline_->rho_g_.derivative_0_data() << " "
           << *meam_sw_spline_->rho_g_.derivative_n_data();
        fs << meam_sw_spline_->rho_g_.OldFormatString();
        auto x = meam_sw_spline_->rho_g_.x_data();
        auto v = meam_sw_spline_->rho_g_.value_data();
        auto d = meam_sw_spline_->rho_g_.second_derivative_value_data();
        for (int n = 0; n < number_of_knots; ++n) {
          fs << x[n] << " " << v[n] << " " << d[n];
        }
      }

      fs << "\n\n\n\n#\n# MEAM SW SPLINE potential file\n#\n";
      fs.close();
    }
  }

  return false;
}

int MEAMImplementation::Compute(
    KIM::ModelCompute const *const model_compute,
    KIM::ModelComputeArguments const *const model_compute_arguments) {
  // KIM API Model Output compute flags
  bool is_compute_energy = false;
  bool is_compute_forces = false;
  bool is_compute_particle_energy = false;
  bool is_compute_virial = false;
  bool is_compute_particle_virial = false;

  // KIM API Model Input
  int const *particle_species_codes = nullptr;
  int const *particle_contributing = nullptr;

  VectorOfSizeDIM const *coordinates = nullptr;

  // KIM API Model Output
  double *energy = nullptr;
  VectorOfSizeDIM *forces = nullptr;
  double *particle_energy = nullptr;
  VectorOfSizeSix *virial = nullptr;
  VectorOfSizeSix *particle_virial = nullptr;

  int ier = SetComputeMutableValues(
      model_compute_arguments, is_compute_energy, is_compute_forces,
      is_compute_particle_energy, is_compute_virial, is_compute_particle_virial,
      particle_species_codes, particle_contributing, coordinates, energy,
      forces, particle_energy, virial, particle_virial);
  if (ier) {
    HELPER_LOG_ERROR("SetComputeMutableValues fails.\n");
    return true;
  }

#include "meam_implementation_compute_dispatch.cpp"

  return ier;
}

int MEAMImplementation::ComputeArgumentsCreate(
    KIM::ModelComputeArgumentsCreate *const model_compute_arguments_create)
    const {
  return RegisterKIMComputeArgumentsSettings(model_compute_arguments_create);
}

int MEAMImplementation::ComputeArgumentsDestroy(
    KIM::ModelComputeArgumentsDestroy *const) const {
  // Nothing to do for this case
  // Everything is good
  return false;
}

std::size_t MEAMImplementation::TotalNumberOfNeighbors(
    KIM::ModelComputeArguments const *const model_compute_arguments,
    int const *const particle_contributing) const {
  int const *neighbors_of_particle;
  std::size_t total_number_of_neighbors = 0;
  for (int i = 0, number_of_neighbors; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    // Setup loop over neighbors of current particle
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];
      if (!(particle_contributing[j] && j < i)) {
        ++total_number_of_neighbors;
      }  // Effective half list
    }    // Loop over j
  }      // Loop over i
  return total_number_of_neighbors;
}

int MEAMImplementation::MaxNumberOfNeighbors(
    KIM::ModelComputeArguments const *const model_compute_arguments,
    int const *const particle_contributing) const {
  int const *neighbors_of_particle;
  int max_number_of_neighbors = 0;
  for (int i = 0, number_of_neighbors; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    if (number_of_neighbors > max_number_of_neighbors) {
      max_number_of_neighbors = number_of_neighbors;
    }
  }  // Loop over i
  return max_number_of_neighbors;
}

// Main compute function

#define KIM_LOGGER_OBJECT_NAME model_compute_arguments
template <bool IsComputeEnergy, bool IsComputeForces,
          bool IsComputeParticleEnergy, bool IsComputeVirial,
          bool IsComputeParticleVirial>
int MEAMImplementation::MeamCCompute(
    KIM::ModelCompute const *const,
    KIM::ModelComputeArguments const *const model_compute_arguments,
    int const *const particle_species_codes,
    int const *const particle_contributing,
    const VectorOfSizeDIM *const coordinates, double *const energy,
    VectorOfSizeDIM *const forces, double *const particle_energy,
    VectorOfSizeSix virial, VectorOfSizeSix *const particle_virial) const {
  // Everything is good
  int ier = false;

  // Nothing to compute
  if (!IsComputeEnergy && !IsComputeForces && !IsComputeParticleEnergy &&
      !IsComputeVirial && !IsComputeParticleVirial) {
    // Everything is good
    return ier;
  }

  // Initialize energy
  if (IsComputeEnergy) {
    *energy = 0.0;
  }

  // Initialize forces
  if (IsComputeForces) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      forces[i][0] = 0.0;
      forces[i][1] = 0.0;
      forces[i][2] = 0.0;
    }
  }

  // Initialize particle energy
  if (IsComputeParticleEnergy) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      particle_energy[i] = 0.0;
    }
  }

  // Initialize virial
  if (IsComputeVirial) {
    virial[0] = 0.0;
    virial[1] = 0.0;
    virial[2] = 0.0;
    virial[3] = 0.0;
    virial[4] = 0.0;
    virial[5] = 0.0;
  }

  // Initialize particle virial
  if (IsComputeParticleVirial) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      particle_virial[i][0] = 0.0;
      particle_virial[i][1] = 0.0;
      particle_virial[i][2] = 0.0;
      particle_virial[i][3] = 0.0;
      particle_virial[i][4] = 0.0;
      particle_virial[i][5] = 0.0;
    }
  }

  // Set the density arrays, and initialize them if necessary
  meam_c_->ResizeDenistyArrays(cached_number_of_particles_);

  // Grow local arrays scrfcn_, and dscrfcn_ if necessary

  // Compute the total number of neighbors of all contributing atoms
  std::size_t const total_number_of_neighbors =
      TotalNumberOfNeighbors(model_compute_arguments, particle_contributing);

  meam_c_->ResizeScreeningArrays(total_number_of_neighbors);

  int number_of_neighbors;
  int const *neighbors_of_particle = nullptr;

  // Initialize the density calculation
  for (int i = 0, offset = 0; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    meam_c_->InitializeDensityCalculation(
        i, number_of_neighbors, neighbors_of_particle, offset, coordinates,
        particle_species_codes, particle_contributing);
  }  // Loop over i

  // Finalize the density calculation
  for (int i = 0; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    int const species_i = particle_species_codes[i];

    double embedding_energy;
    meam_c_->FinalizeDensityCalculation(i, species_i, embedding_energy);

    if (IsComputeEnergy) {
      *energy += embedding_energy;
    }

    if (IsComputeParticleEnergy) {
      particle_energy[i] += embedding_energy;
    }
  }  // Loop over i

  // Temporary sum
  double ts;

  for (int i = 0, offset = 0; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    double const *const scrfcn = &meam_c_->scrfcn_[offset];
    double const *const dscrfcn = &meam_c_->dscrfcn_[offset];

    int const species_i = particle_species_codes[i];

    double const xi = coordinates[i][0];
    double const yi = coordinates[i][1];
    double const zi = coordinates[i][2];

    int nth_neighbor_of_i = -1;

    // Setup loop over neighbors of current particle
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];
      int const contributing_j = particle_contributing[j];

      // Effective half list
      if (!(contributing_j && j < i)) {
        ++nth_neighbor_of_i;

        // Screening function
        double const sij = scrfcn[nth_neighbor_of_i];

        if (special::IsNotZero(sij)) {
          int const species_j = particle_species_codes[j];

          double const xj = coordinates[j][0];
          double const yj = coordinates[j][1];
          double const zj = coordinates[j][2];

          double const dx = xj - xi;
          double const dy = yj - yi;
          double const dz = zj - zi;

          double const rij2 = dx * dx + dy * dy + dz * dz;

          if (rij2 < max_cutoff_squared_) {
            double const rij = std::sqrt(rij2);
            double const rij_inv = 1.0 / rij;

            // Compute pair densities and derivatives

            double rhoa0_i, drhoa0_i;
            double rhoa1_i, drhoa1_i;
            double rhoa2_i, drhoa2_i;
            double rhoa3_i, drhoa3_i;

            double rhoa0_j, drhoa0_j;
            double rhoa1_j, drhoa1_j;
            double rhoa2_j, drhoa2_j;
            double rhoa3_j, drhoa3_j;

            // Eq.4.8
            meam_c_->ComputeAtomicElectronDensities(
                species_i, species_j, rij, rhoa0_i, drhoa0_i, rhoa1_i, drhoa1_i,
                rhoa2_i, drhoa2_i, rhoa3_i, drhoa3_i, rhoa0_j, drhoa0_j,
                rhoa1_j, drhoa1_j, rhoa2_j, drhoa2_j, rhoa3_j, drhoa3_j);

            double const dxdxv = dx * dx * kVoight2dFactor[0];
            double const dxdyv = dx * dy * kVoight2dFactor[1];
            double const dxdzv = dx * dz * kVoight2dFactor[2];
            double const dydyv = dy * dy * kVoight2dFactor[3];
            double const dydzz = dy * dz * kVoight2dFactor[4];
            double const dzdzv = dz * dz * kVoight2dFactor[5];

            double const dxdxdxv = dx * dx * dx * kVoight3dFactor[0];
            double const dxdxdyv = dx * dx * dy * kVoight3dFactor[1];
            double const dxdxdzv = dx * dx * dz * kVoight3dFactor[2];
            double const dxdydyv = dx * dy * dy * kVoight3dFactor[3];
            double const dxdydzv = dx * dy * dz * kVoight3dFactor[4];
            double const dxdzdzv = dx * dz * dz * kVoight3dFactor[5];
            double const dydydyv = dy * dy * dy * kVoight3dFactor[6];
            double const dydydzv = dy * dy * dz * kVoight3dFactor[7];
            double const dydzdzv = dy * dz * dz * kVoight3dFactor[8];
            double const dzdzdzv = dz * dz * dz * kVoight3dFactor[9];

            double const *const arho1_1d_i = &meam_c_->arho1_(i, 0);

            // In Eq.4.30a & 4.30b
            double const y1r_i =
                arho1_1d_i[0] * dx + arho1_1d_i[1] * dy + arho1_1d_i[2] * dz;

            double const *const arho2_1d_i = &meam_c_->arho2_(i, 0);

            // In Eq.4.30d & 4.30e
            ts = arho2_1d_i[0] * dxdxv;
            ts += arho2_1d_i[1] * dxdyv;
            ts += arho2_1d_i[2] * dxdzv;
            ts += arho2_1d_i[3] * dydyv;
            ts += arho2_1d_i[4] * dydzz;
            ts += arho2_1d_i[5] * dzdzv;
            double const y2rr_i = ts;

            double const *const arho3_1d_i = &meam_c_->arho3_(i, 0);

            // In Eq.4.30g & 4.30h
            ts = arho3_1d_i[0] * dxdxdxv;
            ts += arho3_1d_i[1] * dxdxdyv;
            ts += arho3_1d_i[2] * dxdxdzv;
            ts += arho3_1d_i[3] * dxdydyv;
            ts += arho3_1d_i[4] * dxdydzv;
            ts += arho3_1d_i[5] * dxdzdzv;
            ts += arho3_1d_i[6] * dydydyv;
            ts += arho3_1d_i[7] * dydydzv;
            ts += arho3_1d_i[8] * dydzdzv;
            ts += arho3_1d_i[9] * dzdzdzv;
            double const y3rrr_i = ts;

            double const *const arho3b_1d_i = &meam_c_->arho3b_(i, 0);

            // In Eq.4.30g & 4.30h
            double const w3r_i =
                arho3b_1d_i[0] * dx + arho3b_1d_i[1] * dy + arho3b_1d_i[2] * dz;

            // rho0_ terms Eq.4.26a
            double const drho0dr_i = drhoa0_j * sij;

            // rho0_ terms Eq.4.26a
            double const drho0dr_j = drhoa0_i * sij;

            // In Eq.4.30a
            double const c2sr = 2.0 * sij * rij_inv;

            // Eq.4.30a
            double const drho1dr_i =
                c2sr * (drhoa1_j - rhoa1_j * rij_inv) * y1r_i;

            // Eq.4.30c
            double const drho1drm_i0 = c2sr * rhoa1_j * arho1_1d_i[0];
            double const drho1drm_i1 = c2sr * rhoa1_j * arho1_1d_i[1];
            double const drho1drm_i2 = c2sr * rhoa1_j * arho1_1d_i[2];

            // In Eq.4.30d
            double const c2sr2 = 2.0 * sij / rij2;

            // Eq.4.30d
            double const drho2dr_i =
                c2sr2 * (drhoa2_j - 2 * rhoa2_j * rij_inv) * y2rr_i -
                kTwoThird * meam_c_->arho2b_[i] * drhoa2_j * sij;

            // In Eq.4.30f
            double const c4sr2 = 4 * sij / rij2;

            // Eq.4.30f
            ts = arho2_1d_i[kVoight2dIndex[0]] * dx;   // [0][0]
            ts += arho2_1d_i[kVoight2dIndex[1]] * dy;  // [0][1]
            ts += arho2_1d_i[kVoight2dIndex[2]] * dz;  // [0][2]
            double const drho2drm_i0 = ts * c4sr2 * rhoa2_j;

            ts = arho2_1d_i[kVoight2dIndex[3]] * dx;   // [1][0]
            ts += arho2_1d_i[kVoight2dIndex[4]] * dy;  // [1][1]
            ts += arho2_1d_i[kVoight2dIndex[5]] * dz;  // [1][2]
            double const drho2drm_i1 = ts * c4sr2 * rhoa2_j;

            ts = arho2_1d_i[kVoight2dIndex[6]] * dx;   // [2][0]
            ts += arho2_1d_i[kVoight2dIndex[7]] * dy;  // [2][1]
            ts += arho2_1d_i[kVoight2dIndex[8]] * dz;  // [2][2]
            double const drho2drm_i2 = ts * c4sr2 * rhoa2_j;

            // In Eq.4.30g
            double const rij3 = rij * rij2;
            double const c2sr3 = 2 * sij / rij3;
            double const c65sr = 6.0 / 5.0 * sij * rij_inv;

            // Eq.4.30g
            double const drho3dr_i =
                c2sr3 * (drhoa3_j - 3 * rhoa3_j * rij_inv) * y3rrr_i -
                c65sr * (drhoa3_j - rhoa3_j * rij_inv) * w3r_i;

            // In Eq.4.30i
            double const c6sr3 = 6 * sij / rij3;

            // Eq.4.30i
            ts = arho3_1d_i[kVoight3dIndex[0]] * dxdxv;   // [0][0][0]
            ts += arho3_1d_i[kVoight3dIndex[1]] * dxdyv;  // [0][0][1]
            ts += arho3_1d_i[kVoight3dIndex[2]] * dxdzv;  // [0][0][2]
            ts += arho3_1d_i[kVoight3dIndex[4]] * dydyv;  // [0][1][1]
            ts += arho3_1d_i[kVoight3dIndex[5]] * dydzz;  // [0][1][2]
            ts += arho3_1d_i[kVoight3dIndex[8]] * dzdzv;  // [0][2][2]
            double const drho3drm_i0 =
                (c6sr3 * ts - c65sr * arho3b_1d_i[0]) * rhoa3_j;

            ts = arho3_1d_i[kVoight3dIndex[9]] * dxdxv;    // [1][0][0]
            ts += arho3_1d_i[kVoight3dIndex[10]] * dxdyv;  // [1][0][1]
            ts += arho3_1d_i[kVoight3dIndex[11]] * dxdzv;  // [1][0][2]
            ts += arho3_1d_i[kVoight3dIndex[13]] * dydyv;  // [1][1][1]
            ts += arho3_1d_i[kVoight3dIndex[14]] * dydzz;  // [1][1][2]
            ts += arho3_1d_i[kVoight3dIndex[17]] * dzdzv;  // [1][2][2]
            double const drho3drm_i1 =
                (c6sr3 * ts - c65sr * arho3b_1d_i[1]) * rhoa3_j;

            ts = arho3_1d_i[kVoight3dIndex[18]] * dxdxv;   // [2][0][0]
            ts += arho3_1d_i[kVoight3dIndex[19]] * dxdyv;  // [2][0][1]
            ts += arho3_1d_i[kVoight3dIndex[20]] * dxdzv;  // [2][0][2]
            ts += arho3_1d_i[kVoight3dIndex[22]] * dydyv;  // [2][1][1]
            ts += arho3_1d_i[kVoight3dIndex[23]] * dydzz;  // [2][1][2]
            ts += arho3_1d_i[kVoight3dIndex[26]] * dzdzv;  // [2][2][2]
            double const drho3drm_i2 =
                (c6sr3 * ts - c65sr * arho3b_1d_i[2]) * rhoa3_j;

            double y1r_j;
            double y2rr_j;
            double y3rrr_j;
            double w3r_j;
            double drho1dr_j;
            double drho1drm_j0;
            double drho1drm_j1;
            double drho1drm_j2;
            double drho2dr_j;
            double drho2drm_j0;
            double drho2drm_j1;
            double drho2drm_j2;
            double drho3dr_j;
            double drho3drm_j0;
            double drho3drm_j1;
            double drho3drm_j2;

            if (contributing_j) {
              double const *const arho1_1d_j = &meam_c_->arho1_(j, 0);

              // In Eq.4.30a & 4.30b
              y1r_j = -arho1_1d_j[0] * dx;
              y1r_j -= arho1_1d_j[1] * dy;
              y1r_j -= arho1_1d_j[2] * dz;

              double const *const arho2_1d_j = &meam_c_->arho2_(j, 0);

              // In Eq.4.30d & 4.30e
              y2rr_j = arho2_1d_j[0] * dxdxv;
              y2rr_j += arho2_1d_j[1] * dxdyv;
              y2rr_j += arho2_1d_j[2] * dxdzv;
              y2rr_j += arho2_1d_j[3] * dydyv;
              y2rr_j += arho2_1d_j[4] * dydzz;
              y2rr_j += arho2_1d_j[5] * dzdzv;

              double const *const arho3_1d_j = &meam_c_->arho3_(j, 0);

              // In Eq.4.30g & 4.30h
              y3rrr_j = -arho3_1d_j[0] * dxdxdxv;
              y3rrr_j -= arho3_1d_j[1] * dxdxdyv;
              y3rrr_j -= arho3_1d_j[2] * dxdxdzv;
              y3rrr_j -= arho3_1d_j[3] * dxdydyv;
              y3rrr_j -= arho3_1d_j[4] * dxdydzv;
              y3rrr_j -= arho3_1d_j[5] * dxdzdzv;
              y3rrr_j -= arho3_1d_j[6] * dydydyv;
              y3rrr_j -= arho3_1d_j[7] * dydydzv;
              y3rrr_j -= arho3_1d_j[8] * dydzdzv;
              y3rrr_j -= arho3_1d_j[9] * dzdzdzv;

              double const *const arho3b_1d_j = &meam_c_->arho3b_(j, 0);

              // In Eq.4.30g & 4.30h
              w3r_j = -arho3b_1d_j[0] * dx;
              w3r_j -= arho3b_1d_j[1] * dy;
              w3r_j -= arho3b_1d_j[2] * dz;

              // Eq.4.30a
              drho1dr_j = c2sr * (drhoa1_i - rhoa1_i * rij_inv) * y1r_j;

              // Eq.4.30c
              drho1drm_j0 = -c2sr * rhoa1_i * arho1_1d_j[0];
              drho1drm_j1 = -c2sr * rhoa1_i * arho1_1d_j[1];
              drho1drm_j2 = -c2sr * rhoa1_i * arho1_1d_j[2];

              // Eq.4.30d
              drho2dr_j = c2sr2 * (drhoa2_i - 2 * rhoa2_i * rij_inv) * y2rr_j -
                          kTwoThird * meam_c_->arho2b_[j] * drhoa2_i * sij;

              // Eq.4.30f
              ts = -arho2_1d_j[kVoight2dIndex[0]] * dx;  // [0][0]
              ts -= arho2_1d_j[kVoight2dIndex[1]] * dy;  // [0][1]
              ts -= arho2_1d_j[kVoight2dIndex[2]] * dz;  // [0][2]
              drho2drm_j0 = ts * (-c4sr2 * rhoa2_i);

              ts = -arho2_1d_j[kVoight2dIndex[3]] * dx;  // [1][0]
              ts -= arho2_1d_j[kVoight2dIndex[4]] * dy;  // [1][1]
              ts -= arho2_1d_j[kVoight2dIndex[5]] * dz;  // [1][2]
              drho2drm_j1 = ts * (-c4sr2 * rhoa2_i);

              ts = -arho2_1d_j[kVoight2dIndex[6]] * dx;  // [2][0]
              ts -= arho2_1d_j[kVoight2dIndex[7]] * dy;  // [2][1]
              ts -= arho2_1d_j[kVoight2dIndex[8]] * dz;  // [2][2]
              drho2drm_j2 = ts * (-c4sr2 * rhoa2_i);

              // Eq.4.30g
              drho3dr_j = c2sr3 * (drhoa3_i - 3 * rhoa3_i * rij_inv) * y3rrr_j -
                          c65sr * (drhoa3_i - rhoa3_i * rij_inv) * w3r_j;

              // Eq.4.30i
              ts = arho3_1d_j[kVoight3dIndex[0]] * dxdxv;   // [0][0][0]
              ts += arho3_1d_j[kVoight3dIndex[1]] * dxdyv;  // [0][0][1]
              ts += arho3_1d_j[kVoight3dIndex[2]] * dxdzv;  // [0][0][2]
              ts += arho3_1d_j[kVoight3dIndex[4]] * dydyv;  // [0][1][1]
              ts += arho3_1d_j[kVoight3dIndex[5]] * dydzz;  // [0][1][2]
              ts += arho3_1d_j[kVoight3dIndex[8]] * dzdzv;  // [0][2][2]
              drho3drm_j0 = (-c6sr3 * ts + c65sr * arho3b_1d_j[0]) * rhoa3_i;

              ts = arho3_1d_j[kVoight3dIndex[9]] * dxdxv;    // [1][0][0]
              ts += arho3_1d_j[kVoight3dIndex[10]] * dxdyv;  // [1][0][1]
              ts += arho3_1d_j[kVoight3dIndex[11]] * dxdzv;  // [1][0][2]
              ts += arho3_1d_j[kVoight3dIndex[13]] * dydyv;  // [1][1][1]
              ts += arho3_1d_j[kVoight3dIndex[14]] * dydzz;  // [1][1][2]
              ts += arho3_1d_j[kVoight3dIndex[17]] * dzdzv;  // [1][2][2]
              drho3drm_j1 = (-c6sr3 * ts + c65sr * arho3b_1d_j[1]) * rhoa3_i;

              ts = arho3_1d_j[kVoight3dIndex[18]] * dxdxv;   // [2][0][0]
              ts += arho3_1d_j[kVoight3dIndex[19]] * dxdyv;  // [2][0][1]
              ts += arho3_1d_j[kVoight3dIndex[20]] * dxdzv;  // [2][0][2]
              ts += arho3_1d_j[kVoight3dIndex[22]] * dydyv;  // [2][1][1]
              ts += arho3_1d_j[kVoight3dIndex[23]] * dydzz;  // [2][1][2]
              ts += arho3_1d_j[kVoight3dIndex[26]] * dzdzv;  // [2][2][2]
              drho3drm_j2 = (-c6sr3 * ts + c65sr * arho3b_1d_j[2]) * rhoa3_i;
            }  // particle_contributing

            // Compute derivatives of weighting functions t wrt rij

            double const t1m_i = meam_c_->element_t1_[species_i];
            double const t1m_j = meam_c_->element_t1_[species_j];

            double const t2m_i = meam_c_->element_t2_[species_i];
            double const t2m_j = meam_c_->element_t2_[species_j];

            double const t3m_i = meam_c_->element_t3_[species_i];
            double const t3m_j = meam_c_->element_t3_[species_j];

            double const t1_i = meam_c_->t_ave_(i, 0);
            double const t1_j = meam_c_->t_ave_(j, 0);

            double const t2_i = meam_c_->t_ave_(i, 1);
            double const t2_j = meam_c_->t_ave_(j, 1);

            double const t3_i = meam_c_->t_ave_(i, 2);
            double const t3_j = meam_c_->t_ave_(j, 2);

            double dt1dr_i;
            double dt2dr_i;
            double dt3dr_i;

            if (meam_c_->ialloy_ == 1) {
              double const *const tsq_ave_1d_i = &meam_c_->tsq_ave_(i, 0);

              double const a1i =
                  special::FloatDivZero(drhoa0_j * sij, tsq_ave_1d_i[0]);
              double const a2i =
                  special::FloatDivZero(drhoa0_j * sij, tsq_ave_1d_i[1]);
              double const a3i =
                  special::FloatDivZero(drhoa0_j * sij, tsq_ave_1d_i[2]);

              dt1dr_i = a1i * (t1m_j - t1_i * special::Square(t1m_j));
              dt2dr_i = a2i * (t2m_j - t2_i * special::Square(t2m_j));
              dt3dr_i = a3i * (t3m_j - t3_i * special::Square(t3m_j));
            } else if (meam_c_->ialloy_ == 2) {
              dt1dr_i = 0.0;
              dt2dr_i = 0.0;
              dt3dr_i = 0.0;
            } else {
              double const ai =
                  special::FloatDivZero(drhoa0_j * sij, meam_c_->rho0_[i]);
              dt1dr_i = ai * (t1m_j - t1_i);
              dt2dr_i = ai * (t2m_j - t2_i);
              dt3dr_i = ai * (t3m_j - t3_i);
            }

            // Compute derivatives of total density wrt rij, sij and rij(3)

            VectorOfSizeDIM shape_factors_i;

            meam_c_->GetShapeFactors(
                meam_c_->element_lattice_(species_i, species_i),
                meam_c_->element_stheta_(species_i, species_i),
                meam_c_->element_ctheta_(species_i, species_i),
                shape_factors_i);

            // Eq.4.36a
            double const drhodr_i =
                meam_c_->dgamma1_[i] * drho0dr_i +
                meam_c_->dgamma2_[i] *
                    (dt1dr_i * meam_c_->rho1_[i] + t1_i * drho1dr_i +
                     dt2dr_i * meam_c_->rho2_[i] + t2_i * drho2dr_i +
                     dt3dr_i * meam_c_->rho3_[i] + t3_i * drho3dr_i) -
                meam_c_->dgamma3_[i] * (shape_factors_i[0] * dt1dr_i +
                                        shape_factors_i[1] * dt2dr_i +
                                        shape_factors_i[2] * dt3dr_i);

            // Eq.4.36c
            double const drhodrm_i0 =
                meam_c_->dgamma2_[i] *
                (t1_i * drho1drm_i0 + t2_i * drho2drm_i0 + t3_i * drho3drm_i0);
            double const drhodrm_i1 =
                meam_c_->dgamma2_[i] *
                (t1_i * drho1drm_i1 + t2_i * drho2drm_i1 + t3_i * drho3drm_i1);
            double const drhodrm_i2 =
                meam_c_->dgamma2_[i] *
                (t1_i * drho1drm_i2 + t2_i * drho2drm_i2 + t3_i * drho3drm_i2);

            VectorOfSizeDIM shape_factors_j;
            double drhodr_j{};
            double drhodrm_j0{};
            double drhodrm_j1{};
            double drhodrm_j2{};

            if (contributing_j) {
              double dt1dr_j;
              double dt2dr_j;
              double dt3dr_j;

              if (meam_c_->ialloy_ == 1) {
                double const *const tsq_ave_1d_j = &meam_c_->tsq_ave_(j, 0);

                double const a1j =
                    special::FloatDivZero(drhoa0_i * sij, tsq_ave_1d_j[0]);
                double const a2j =
                    special::FloatDivZero(drhoa0_i * sij, tsq_ave_1d_j[1]);
                double const a3j =
                    special::FloatDivZero(drhoa0_i * sij, tsq_ave_1d_j[2]);

                dt1dr_j = a1j * t1m_i * (1.0 - t1_j * t1m_i);
                dt2dr_j = a2j * t2m_i * (1.0 - t2_j * t2m_i);
                dt3dr_j = a3j * t3m_i * (1.0 - t3_j * t3m_i);

              } else if (meam_c_->ialloy_ == 2) {
                dt1dr_j = 0.0;
                dt2dr_j = 0.0;
                dt3dr_j = 0.0;

              } else {
                double const aj =
                    special::FloatDivZero(drhoa0_i * sij, meam_c_->rho0_[j]);
                dt1dr_j = aj * (t1m_i - t1_j);
                dt2dr_j = aj * (t2m_i - t2_j);
                dt3dr_j = aj * (t3m_i - t3_j);
              }

              meam_c_->GetShapeFactors(
                  meam_c_->element_lattice_(species_j, species_j),
                  meam_c_->element_stheta_(species_i, species_i),
                  meam_c_->element_ctheta_(species_i, species_i),
                  shape_factors_j);

              // Eq.4.36a
              drhodr_j = meam_c_->dgamma1_[j] * drho0dr_j +
                         meam_c_->dgamma2_[j] *
                             (dt1dr_j * meam_c_->rho1_[j] + t1_j * drho1dr_j +
                              dt2dr_j * meam_c_->rho2_[j] + t2_j * drho2dr_j +
                              dt3dr_j * meam_c_->rho3_[j] + t3_j * drho3dr_j) -
                         meam_c_->dgamma3_[j] * (shape_factors_j[0] * dt1dr_j +
                                                 shape_factors_j[1] * dt2dr_j +
                                                 shape_factors_j[2] * dt3dr_j);

              // Eq.4.36c
              drhodrm_j0 = meam_c_->dgamma2_[j] *
                           (t1_j * drho1drm_j0 + t2_j * drho2drm_j0 +
                            t3_j * drho3drm_j0);
              drhodrm_j1 = meam_c_->dgamma2_[j] *
                           (t1_j * drho1drm_j1 + t2_j * drho2drm_j1 +
                            t3_j * drho3drm_j1);
              drhodrm_j2 = meam_c_->dgamma2_[j] *
                           (t1_j * drho1drm_j2 + t2_j * drho2drm_j2 +
                            t3_j * drho3drm_j2);
            }  // particle_contributing

            double dphi;
            double const phi =
                meam_c_->GetPhiAndDerivative(species_i, species_j, rij, dphi);

            if (contributing_j) {
              // Contribute pair term to Energy
              if (IsComputeEnergy) {
                *energy += phi * sij;
              }

              // Contribute pair term to Particle Energy
              if (IsComputeParticleEnergy) {
                particle_energy[i] += 0.5 * phi * sij;
                particle_energy[j] += 0.5 * phi * sij;
              }
            }
            // particle j is not contributing
            else {
              // Contribute pair term to Energy
              if (IsComputeEnergy) {
                *energy += 0.5 * phi * sij;
              }

              // Contribute pair term to Particle Energy
              if (IsComputeParticleEnergy) {
                particle_energy[i] += 0.5 * phi * sij;
              }
            }

            // Derivative of the screening function
            double const dsij = dscrfcn[nth_neighbor_of_i];

            double dEdsij = 0.0;

            // Compute derivatives wrt sij, but only if necessary
            if (special::IsNotZero(dsij)) {
              // Eq.4.26b
              double const drho0ds_i = rhoa0_j;

              // Eq.4.26b
              double const drho0ds_j = rhoa0_i;

              // in Eq.4.30b
              double const c2r = 2.0 * rij_inv;

              // Eq.4.30b
              double const drho1ds_i = c2r * rhoa1_j * y1r_i;

              // in Eq.4.30e
              double const c2r2 = 2.0 / rij2;

              // Eq.4.30e
              double const drho2ds_i =
                  c2r2 * rhoa2_j * y2rr_i -
                  kTwoThird * meam_c_->arho2b_[i] * rhoa2_j;

              // in Eq.4.30h
              double const c2r3 = 2.0 / rij3;
              double const c65r = 6.0 / (5.0 * rij);

              // Eq.4.30h
              double const drho3ds_i =
                  c2r3 * rhoa3_j * y3rrr_i - c65r * rhoa3_j * w3r_i;

              double dt1ds_i;
              double dt2ds_i;
              double dt3ds_i;

              if (meam_c_->ialloy_ == 1) {
                double const *const tsq_ave_1d_i = &meam_c_->tsq_ave_(i, 0);

                double const a1i =
                    special::FloatDivZero(rhoa0_j, tsq_ave_1d_i[0]);
                double const a2i =
                    special::FloatDivZero(rhoa0_j, tsq_ave_1d_i[1]);
                double const a3i =
                    special::FloatDivZero(rhoa0_j, tsq_ave_1d_i[2]);

                dt1ds_i = a1i * t1m_j * (1.0 - t1_i * t1m_j);
                dt2ds_i = a2i * t2m_j * (1.0 - t2_i * t2m_j);
                dt3ds_i = a3i * t3m_j * (1.0 - t3_i * t3m_j);

              } else if (meam_c_->ialloy_ == 2) {
                dt1ds_i = 0.0;
                dt2ds_i = 0.0;
                dt3ds_i = 0.0;

              } else {
                double const ai =
                    special::FloatDivZero(rhoa0_j, meam_c_->rho0_[i]);

                dt1ds_i = ai * (t1m_j - t1_i);
                dt2ds_i = ai * (t2m_j - t2_i);
                dt3ds_i = ai * (t3m_j - t3_i);
              }

              // Eq.4.36b
              double const drhods_i =
                  meam_c_->dgamma1_[i] * drho0ds_i +
                  meam_c_->dgamma2_[i] *
                      (dt1ds_i * meam_c_->rho1_[i] + t1_i * drho1ds_i +
                       dt2ds_i * meam_c_->rho2_[i] + t2_i * drho2ds_i +
                       dt3ds_i * meam_c_->rho3_[i] + t3_i * drho3ds_i) -
                  meam_c_->dgamma3_[i] * (shape_factors_i[0] * dt1ds_i +
                                          shape_factors_i[1] * dt2ds_i +
                                          shape_factors_i[2] * dt3ds_i);

              if (contributing_j) {
                // Eq.4.30b
                double const drho1ds_j = c2r * rhoa1_i * y1r_j;

                // Eq.4.30e
                double const drho2ds_j =
                    c2r2 * rhoa2_i * y2rr_j -
                    kTwoThird * meam_c_->arho2b_[j] * rhoa2_i;

                // Eq.4.30h
                double const drho3ds_j =
                    c2r3 * rhoa3_i * y3rrr_j - c65r * rhoa3_i * w3r_j;

                double dt1ds_j;
                double dt2ds_j;
                double dt3ds_j;

                if (meam_c_->ialloy_ == 1) {
                  double const *const tsq_ave_1d_j = &meam_c_->tsq_ave_(j, 0);

                  double const a1j =
                      special::FloatDivZero(rhoa0_i, tsq_ave_1d_j[0]);
                  double const a2j =
                      special::FloatDivZero(rhoa0_i, tsq_ave_1d_j[1]);
                  double const a3j =
                      special::FloatDivZero(rhoa0_i, tsq_ave_1d_j[2]);

                  dt1ds_j = a1j * t1m_i * (1.0 - t1_j * t1m_i);
                  dt2ds_j = a2j * t2m_i * (1.0 - t2_j * t2m_i);
                  dt3ds_j = a3j * t3m_i * (1.0 - t3_j * t3m_i);

                } else if (meam_c_->ialloy_ == 2) {
                  dt1ds_j = 0.0;
                  dt2ds_j = 0.0;
                  dt3ds_j = 0.0;

                } else {
                  double const aj =
                      special::FloatDivZero(rhoa0_i, meam_c_->rho0_[j]);

                  dt1ds_j = aj * (t1m_i - t1_j);
                  dt2ds_j = aj * (t2m_i - t2_j);
                  dt3ds_j = aj * (t3m_i - t3_j);
                }

                // Eq.4.36b
                double const drhods_j =
                    meam_c_->dgamma1_[j] * drho0ds_j +
                    meam_c_->dgamma2_[j] *
                        (dt1ds_j * meam_c_->rho1_[j] + t1_j * drho1ds_j +
                         dt2ds_j * meam_c_->rho2_[j] + t2_j * drho2ds_j +
                         dt3ds_j * meam_c_->rho3_[j] + t3_j * drho3ds_j) -
                    meam_c_->dgamma3_[j] * (shape_factors_j[0] * dt1ds_j +
                                            shape_factors_j[1] * dt2ds_j +
                                            shape_factors_j[2] * dt3ds_j);

                // Eq.4.41b
                dEdsij = phi + meam_c_->frhop_[i] * drhods_i +
                         meam_c_->frhop_[j] * drhods_j;
              }
              // particle j is not contributing
              else {
                // Eq.4.41b
                dEdsij = 0.5 * phi + meam_c_->frhop_[i] * drhods_i;
              }  // particle_contributing
            }    // dscrfcn is not zero

            if (IsComputeForces || IsComputeVirial || IsComputeParticleVirial) {
              // Compute derivatives of energy wrt rij, sij and rij[3]

              // Eq.4.41a
              double dEdrij;

              // Eq.4.41d
              double dEdrijm0;
              double dEdrijm1;
              double dEdrijm2;

              if (contributing_j) {
                // Eq.4.41a
                dEdrij = dphi * sij + meam_c_->frhop_[i] * drhodr_i +
                         meam_c_->frhop_[j] * drhodr_j;

                // Eq.4.41d
                dEdrijm0 = meam_c_->frhop_[i] * drhodrm_i0 +
                           meam_c_->frhop_[j] * drhodrm_j0;
                dEdrijm1 = meam_c_->frhop_[i] * drhodrm_i1 +
                           meam_c_->frhop_[j] * drhodrm_j1;
                dEdrijm2 = meam_c_->frhop_[i] * drhodrm_i2 +
                           meam_c_->frhop_[j] * drhodrm_j2;

              }
              // particle j is not contributing
              else {
                // Eq.4.41a
                dEdrij = 0.5 * dphi * sij + meam_c_->frhop_[i] * drhodr_i;

                // Eq.4.41d
                dEdrijm0 = meam_c_->frhop_[i] * drhodrm_i0;
                dEdrijm1 = meam_c_->frhop_[i] * drhodrm_i1;
                dEdrijm2 = meam_c_->frhop_[i] * drhodrm_i2;
              }

              // Add the part of the force due to dEdrij and dEdsij
              double const force = dEdrij * rij_inv + dEdsij * dsij;

              double const fi0 = dx * force + dEdrijm0;
              double const fi1 = dy * force + dEdrijm1;
              double const fi2 = dz * force + dEdrijm2;

              if (IsComputeForces) {
                forces[i][0] += fi0;
                forces[i][1] += fi1;
                forces[i][2] += fi2;

                forces[j][0] -= fi0;
                forces[j][1] -= fi1;
                forces[j][2] -= fi2;
              }

              if (IsComputeVirial || IsComputeParticleVirial) {
                // Virial has 6 components and is stored as a 6-element
                // vector in the following order: xx, yy, zz, yz, xz, xy.
                VectorOfSizeSix v;
                v[0] = dx * fi0;
                v[1] = dy * fi1;
                v[2] = dz * fi2;
                v[3] = 0.5 * (dy * fi2 + dz * fi1);
                v[4] = 0.5 * (dx * fi2 + dz * fi0);
                v[5] = 0.5 * (dx * fi1 + dy * fi0);

                if (IsComputeVirial) {
                  virial[0] += v[0];
                  virial[1] += v[1];
                  virial[2] += v[2];
                  virial[3] += v[3];
                  virial[4] += v[4];
                  virial[5] += v[5];
                }

                if (IsComputeParticleVirial) {
                  v[0] *= 0.5;
                  v[1] *= 0.5;
                  v[2] *= 0.5;
                  v[3] *= 0.5;
                  v[4] *= 0.5;
                  v[5] *= 0.5;

                  particle_virial[i][0] += v[0];
                  particle_virial[i][1] += v[1];
                  particle_virial[i][2] += v[2];
                  particle_virial[i][3] += v[3];
                  particle_virial[i][4] += v[4];
                  particle_virial[i][5] += v[5];

                  particle_virial[j][0] += v[0];
                  particle_virial[j][1] += v[1];
                  particle_virial[j][2] += v[2];
                  particle_virial[j][3] += v[3];
                  particle_virial[j][4] += v[4];
                  particle_virial[j][5] += v[5];
                }  // IsComputeParticleVirial
              }    // IsComputeVirial || IsComputeParticleVirial

              // Now compute forces on other atoms k due to the change in sij

              if (special::IsZero(sij) || special::IsOne(sij)) {
                // Loop over j
                continue;
              }

              double dx_jk;
              double dy_jk;
              double dz_jk;

              double dx_ik;
              double dy_ik;
              double dz_ik;

              // Setup loop over neighbors of current particle
              for (int kn = 0; kn < number_of_neighbors; ++kn) {
                int const k = neighbors_of_particle[kn];
                if (k == j) {
                  continue;
                }

                double dsij1 = 0.0;
                double dsij2 = 0.0;

                if (special::IsNotZero(sij) && special::IsNotOne(sij)) {
                  double const rbound =
                      meam_c_->element_ebound_(species_i, species_j) * rij2;

                  double const xk = coordinates[k][0];
                  double const yk = coordinates[k][1];
                  double const zk = coordinates[k][2];

                  dx_jk = xk - xj;
                  dy_jk = yk - yj;
                  dz_jk = zk - zj;

                  double const rjk2 =
                      dx_jk * dx_jk + dy_jk * dy_jk + dz_jk * dz_jk;

                  // lies within the screening function ellipse,
                  // so we need to calculate its effects
                  if (rjk2 <= rbound) {
                    dx_ik = xk - xi;
                    dy_ik = yk - yi;
                    dz_ik = zk - zi;

                    double const rik2 =
                        dx_ik * dx_ik + dy_ik * dy_ik + dz_ik * dz_ik;

                    // lies within the screening function ellipse,
                    // so we need to calculate its effects
                    if (rik2 <= rbound) {
                      double const xik = rik2 / rij2;
                      double const xjk = rjk2 / rij2;

                      double a = 1 - (xik - xjk) * (xik - xjk);

                      if (special::IsNotZero(a)) {
                        double cijk = (2.0 * (xik + xjk) + a - 2.0) / a;

                        int const species_k = particle_species_codes[k];

                        double const cmax = meam_c_->element_Cmax_(
                            species_i, species_j, species_k);
                        double const cmin = meam_c_->element_Cmin_(
                            species_i, species_j, species_k);

                        if (cijk >= cmin && cijk <= cmax) {
                          double const delc = cmax - cmin;

                          cijk = (cijk - cmin) / delc;

                          double dfc_drij;
                          double const sijk = meam_c_->RadialCutoff(
                              cijk, dfc_drij, meam_c_->fcut_function_form_);

                          // Eq.4.17b
                          double dcijk_drik;
                          // Eq.4.17c
                          double dcijk_drjk;

                          meam_c_->DCijkDRikDRjk(rij2, rik2, rjk2, dcijk_drik,
                                                 dcijk_drjk);

                          a = sij / delc * dfc_drij / sijk;

                          // Eq.4.22b
                          dsij1 = a * dcijk_drik;
                          // Eq.4.22c
                          dsij2 = a * dcijk_drjk;
                        }  // cijk >= cmin && cijk <= cmax
                      }    // IsNotZero(a)
                    }      // rik2 <= rbound
                  }        // rjk2 <= rbound
                }          //  IsNotZero(sij) && IsNotOne(sij)

                if (special::IsNotZero(dsij1) || special::IsNotZero(dsij2)) {
                  // In Eq.4.40
                  double const force1 = dEdsij * dsij1;
                  double const fi0 = force1 * dx_ik;
                  double const fi1 = force1 * dy_ik;
                  double const fi2 = force1 * dz_ik;

                  // In Eq.4.40
                  double const force2 = dEdsij * dsij2;
                  double const fj0 = force2 * dx_jk;
                  double const fj1 = force2 * dy_jk;
                  double const fj2 = force2 * dz_jk;

                  if (IsComputeForces) {
                    forces[i][0] += fi0;
                    forces[i][1] += fi1;
                    forces[i][2] += fi2;

                    forces[j][0] += fj0;
                    forces[j][1] += fj1;
                    forces[j][2] += fj2;

                    forces[k][0] -= (fi0 + fj0);
                    forces[k][1] -= (fi1 + fj1);
                    forces[k][2] -= (fi2 + fj2);
                  }

                  if (IsComputeVirial || IsComputeParticleVirial) {
                    // Tabulate per-atom virial as symmetrized stress
                    // tensor

                    // Virial has 6 components and is stored as a
                    // 6 - element vector in the following order:
                    // xx, yy, zz, yz, xz, xy.

                    VectorOfSizeSix v;
                    v[0] = kTwoThird * (dx_ik * fi0 + dx_jk * fj0);
                    v[1] = kTwoThird * (dy_ik * fi1 + dy_jk * fj1);
                    v[2] = kTwoThird * (dz_ik * fi2 + dz_jk * fj2);
                    v[3] = kThird * (dy_ik * fi2 + dy_jk * fj2 + dz_ik * fi1 +
                                     dz_jk * fj1);
                    v[4] = kThird * (dx_ik * fi2 + dx_jk * fj2 + dz_ik * fi0 +
                                     dz_jk * fj0);
                    v[5] = kThird * (dx_ik * fi1 + dx_jk * fj1 + dy_ik * fi0 +
                                     dy_jk * fj0);

                    if (IsComputeVirial) {
                      virial[0] += v[0];
                      virial[1] += v[1];
                      virial[2] += v[2];
                      virial[3] += v[3];
                      virial[4] += v[4];
                      virial[5] += v[5];
                    }  // IsComputeVirial

                    if (IsComputeParticleVirial) {
                      v[0] *= 0.5;
                      v[1] *= 0.5;
                      v[2] *= 0.5;
                      v[3] *= 0.5;
                      v[4] *= 0.5;
                      v[5] *= 0.5;

                      particle_virial[i][0] += v[0];
                      particle_virial[i][1] += v[1];
                      particle_virial[i][2] += v[2];
                      particle_virial[i][3] += v[3];
                      particle_virial[i][4] += v[4];
                      particle_virial[i][5] += v[5];

                      particle_virial[j][0] += v[0];
                      particle_virial[j][1] += v[1];
                      particle_virial[j][2] += v[2];
                      particle_virial[j][3] += v[3];
                      particle_virial[j][4] += v[4];
                      particle_virial[j][5] += v[5];

                      particle_virial[k][0] += v[0];
                      particle_virial[k][1] += v[1];
                      particle_virial[k][2] += v[2];
                      particle_virial[k][3] += v[3];
                      particle_virial[k][4] += v[4];
                      particle_virial[k][5] += v[5];
                    }  // IsComputeParticleVirial
                  }    // IsComputeVirial || IsComputeParticleVirial
                }      // IsNotZero(dsij1) || IsNotZero(dsij2)
              }        // Loop over k
            }          // IsComputeForces || IsComputeVirial ||
                       // IsComputeParticleVirial
          }            // rij2 < max_cutoff_squared_
        }              // scrfcn is not zero
      }                // Effective half list
    }                  // Loop over j

    if (nth_neighbor_of_i > -1) {
      offset += nth_neighbor_of_i + 1;
    }
  }  // Loop over i

  // everything is good
  return false;
}

#define KIM_LOGGER_OBJECT_NAME model_compute_arguments
template <bool IsComputeEnergy, bool IsComputeForces,
          bool IsComputeParticleEnergy, bool IsComputeVirial,
          bool IsComputeParticleVirial, bool AreKnotsOnRegularGid>
int MEAMImplementation::MeamSplineCompute(
    KIM::ModelCompute const *const,
    KIM::ModelComputeArguments const *const model_compute_arguments,
    int const *const particle_species_codes,
    int const *const particle_contributing,
    const VectorOfSizeDIM *const coordinates, double *const energy,
    VectorOfSizeDIM *const forces, double *const particle_energy,
    VectorOfSizeSix virial, VectorOfSizeSix *const particle_virial) const {
  // Everything is good
  int ier = false;

  // Nothing to compute
  if (!IsComputeEnergy && !IsComputeForces && !IsComputeParticleEnergy &&
      !IsComputeVirial && !IsComputeParticleVirial) {
    // Everything is good
    return ier;
  }

  // Initialize energy
  if (IsComputeEnergy) {
    *energy = 0.0;
  }

  // Initialize forces
  if (IsComputeForces) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      forces[i][0] = 0.0;
      forces[i][1] = 0.0;
      forces[i][2] = 0.0;
    }
  }

  // Initialize particle energy
  if (IsComputeParticleEnergy) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      particle_energy[i] = 0.0;
    }
  }

  // Initialize virial
  if (IsComputeVirial) {
    virial[0] = 0.0;
    virial[1] = 0.0;
    virial[2] = 0.0;
    virial[3] = 0.0;
    virial[4] = 0.0;
    virial[5] = 0.0;
  }

  // Initialize particle virial
  if (IsComputeParticleVirial) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      particle_virial[i][0] = 0.0;
      particle_virial[i][1] = 0.0;
      particle_virial[i][2] = 0.0;
      particle_virial[i][3] = 0.0;
      particle_virial[i][4] = 0.0;
      particle_virial[i][5] = 0.0;
    }
  }

  // Grow per-atom array if necessary
  meam_spline_->ResizeLocalArray(cached_number_of_particles_);

  // Determine the maximum number of neighbors a single atom has
  int const max_number_of_neighbors =
      MaxNumberOfNeighbors(model_compute_arguments, particle_contributing);

  // Temporary array
  std::vector<TwoBody> two_body(max_number_of_neighbors);

  int number_of_neighbors;
  int const *neighbors_of_particle;

  // Temporary variable
  double prime;

  // Sum three-body contributions to charge density and
  // the embedding energy
  for (int i = 0; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    int const species_i = particle_species_codes[i];

    double const xi = coordinates[i][0];
    double const yi = coordinates[i][1];
    double const zi = coordinates[i][2];

    // Compute number_bonds & charge density
    int number_bonds = 0;

    double rho_value = 0.0;

    // Setup loop over neighbors of the current particle
    auto bond_j = two_body.begin();
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];

      double const dx = coordinates[j][0] - xi;
      double const dy = coordinates[j][1] - yi;
      double const dz = coordinates[j][2] - zi;

      double const rij2 = dx * dx + dy * dy + dz * dz;

      if (rij2 < max_cutoff_squared_) {
        double const rij = std::sqrt(rij2);

        int const species_j = particle_species_codes[j];

        bond_j->index = j;
        bond_j->species = species_j;
        bond_j->r = rij;
        bond_j->f = meam_spline_->rho_f_[species_j].Eval<AreKnotsOnRegularGid>(
            rij, prime);
        bond_j->f_prime = prime;
        bond_j->dx = dx / rij;
        bond_j->dy = dy / rij;
        bond_j->dz = dz / rij;

        int const species_ind =
            species_j * number_of_elements_ - species_j * (species_j + 1) / 2;

        double partial_sum = 0.0;

        auto bond_k = two_body.begin();
        for (int kb = 0; kb < number_bonds; ++kb, ++bond_k) {
          int const ind = bond_k->species + species_ind;

          double const cos_theta = bond_j->dx * bond_k->dx +
                                   bond_j->dy * bond_k->dy +
                                   bond_j->dz * bond_k->dz;

          double const rho_g =
              meam_spline_->rho_g_[ind].Eval<AreKnotsOnRegularGid>(cos_theta);

          partial_sum += bond_k->f * rho_g;
        }

        rho_value += bond_j->f * partial_sum;

        double const rho_r =
            meam_spline_->rho_r_[species_j].Eval<AreKnotsOnRegularGid>(rij);
        rho_value += rho_r;

        ++number_bonds;
        ++bond_j;
      }  // rij2 < max_cutoff_squared_
    }    // Loop over j

    // Compute embedding energy and its derivative
    double e_u_prime;

    double const e_u = meam_spline_->e_u_[species_i].Eval<AreKnotsOnRegularGid>(
                           rho_value, e_u_prime) -
                       meam_spline_->zero_atom_energies_[species_i];

    meam_spline_->e_u_prime_[i] = e_u_prime;

    if (IsComputeEnergy) {
      *energy += e_u;
    }

    if (IsComputeParticleEnergy) {
      particle_energy[i] += e_u;
    }

    // Compute three-body contributions to force
    VectorOfSizeDIM f_i = {};

    bond_j = two_body.begin();
    for (int jb = 0; jb < number_bonds; ++jb, ++bond_j) {
      int const j = bond_j->index;
      int const species_j = bond_j->species;
      double const r_ij = bond_j->r;
      double const f_ij = bond_j->f;
      double const f_prime_ij = bond_j->f_prime;

      int const species_ind =
          species_j * number_of_elements_ - species_j * (species_j + 1) / 2;

      VectorOfSizeDIM f_j = {};

      auto bond_k = two_body.begin();
      for (int kb = 0; kb < jb; ++kb, ++bond_k) {
        int const k = bond_k->index;
        double const r_ik = bond_k->r;
        double const f_ik = bond_k->f;
        double const f_prime_ik = bond_k->f_prime;

        int const ind = bond_k->species + species_ind;

        double const cos_theta = bond_j->dx * bond_k->dx +
                                 bond_j->dy * bond_k->dy +
                                 bond_j->dz * bond_k->dz;

        double const rho_g =
            meam_spline_->rho_g_[ind].Eval<AreKnotsOnRegularGid>(cos_theta,
                                                                 prime);

        double const prefactor = e_u_prime * f_ij * prime * f_ik;

        double const prefactor_ij = prefactor / r_ij;
        double const prefactor_ik = prefactor / r_ik;

        double const fij =
            -e_u_prime * rho_g * f_prime_ij * f_ik + prefactor_ij * cos_theta;

        double const fj0 = bond_j->dx * fij - bond_k->dx * prefactor_ij;
        double const fj1 = bond_j->dy * fij - bond_k->dy * prefactor_ij;
        double const fj2 = bond_j->dz * fij - bond_k->dz * prefactor_ij;

        f_j[0] += fj0;
        f_j[1] += fj1;
        f_j[2] += fj2;

        double const fik =
            -e_u_prime * rho_g * f_prime_ik * f_ij + prefactor_ik * cos_theta;

        double const fk0 = bond_k->dx * fik - bond_j->dx * prefactor_ik;
        double const fk1 = bond_k->dy * fik - bond_j->dy * prefactor_ik;
        double const fk2 = bond_k->dz * fik - bond_j->dz * prefactor_ik;

        f_i[0] -= fk0;
        f_i[1] -= fk1;
        f_i[2] -= fk2;

        if (IsComputeForces) {
          forces[k][0] += fk0;
          forces[k][1] += fk1;
          forces[k][2] += fk2;
        }

        if (IsComputeVirial || IsComputeParticleVirial) {
          // Tabulate per-atom virial as symmetrized stress
          // tensor

          // Virial has 6 components and is stored as a
          // 6 - element vector in the following order:
          // xx, yy, zz, yz, xz, xy.

          double const dx_ij = bond_j->dx * r_ij;
          double const dy_ij = bond_j->dy * r_ij;
          double const dz_ij = bond_j->dz * r_ij;

          double const dx_ik = bond_k->dx * r_ik;
          double const dy_ik = bond_k->dy * r_ik;
          double const dz_ik = bond_k->dz * r_ik;

          VectorOfSizeSix v;
          v[0] = dx_ij * fj0 + dx_ik * fk0;
          v[1] = dy_ij * fj1 + dy_ik * fk1;
          v[2] = dz_ij * fj2 + dz_ik * fk2;
          v[3] = dy_ij * fj2 + dy_ik * fk2;
          v[4] = dx_ij * fj2 + dx_ik * fk2;
          v[5] = dx_ij * fj1 + dx_ik * fk1;

          if (IsComputeVirial) {
            virial[0] += v[0];
            virial[1] += v[1];
            virial[2] += v[2];
            virial[3] += v[3];
            virial[4] += v[4];
            virial[5] += v[5];
          }  // IsComputeVirial

          if (IsComputeParticleVirial) {
            v[0] *= kThird;
            v[1] *= kThird;
            v[2] *= kThird;
            v[3] *= kThird;
            v[4] *= kThird;
            v[5] *= kThird;

            particle_virial[i][0] += v[0];
            particle_virial[i][1] += v[1];
            particle_virial[i][2] += v[2];
            particle_virial[i][3] += v[3];
            particle_virial[i][4] += v[4];
            particle_virial[i][5] += v[5];

            particle_virial[j][0] += v[0];
            particle_virial[j][1] += v[1];
            particle_virial[j][2] += v[2];
            particle_virial[j][3] += v[3];
            particle_virial[j][4] += v[4];
            particle_virial[j][5] += v[5];

            particle_virial[k][0] += v[0];
            particle_virial[k][1] += v[1];
            particle_virial[k][2] += v[2];
            particle_virial[k][3] += v[3];
            particle_virial[k][4] += v[4];
            particle_virial[k][5] += v[5];
          }  // IsComputeParticleVirial
        }    // IsComputeVirial || IsComputeParticleVirial
      }      // Loop over k bonds

      if (IsComputeForces) {
        forces[i][0] -= f_j[0];
        forces[i][1] -= f_j[1];
        forces[i][2] -= f_j[2];

        forces[j][0] += f_j[0];
        forces[j][1] += f_j[1];
        forces[j][2] += f_j[2];
      }
    }  // Loop over j bonds

    if (IsComputeForces) {
      forces[i][0] += f_i[0];
      forces[i][1] += f_i[1];
      forces[i][2] += f_i[2];
    }
  }  // Loop over i

  // Compute two-body pair interactions
  for (int i = 0; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    int const species_i = particle_species_codes[i];
    int const species_ind =
        species_i * number_of_elements_ - species_i * (species_i + 1) / 2;

    double const xi = coordinates[i][0];
    double const yi = coordinates[i][1];
    double const zi = coordinates[i][2];

    // Setup loop over neighbors of current particle
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];
      int const contributing_j = particle_contributing[j];

      // Effective half list
      if (!(contributing_j && j < i)) {
        double const dx = coordinates[j][0] - xi;
        double const dy = coordinates[j][1] - yi;
        double const dz = coordinates[j][2] - zi;

        double const rij2 = dx * dx + dy * dy + dz * dz;

        if (rij2 < max_cutoff_squared_) {
          double const rij = std::sqrt(rij2);

          int const species_j = particle_species_codes[j];

          int const ind = species_j + species_ind;

          double fij;
          double const phi =
              meam_spline_->e_phi_[ind].Eval<AreKnotsOnRegularGid>(rij, fij);

          meam_spline_->rho_r_[species_j].Eval<AreKnotsOnRegularGid>(rij,
                                                                     prime);

          if (contributing_j) {
            fij += prime * meam_spline_->e_u_prime_[i];

            if (species_i != species_j) {
              meam_spline_->rho_r_[species_i].Eval<AreKnotsOnRegularGid>(rij,
                                                                         prime);
            }

            fij += prime * meam_spline_->e_u_prime_[j];

            // Contribute pair term to Energy
            if (IsComputeEnergy) {
              *energy += phi;
            }

            // Contribute pair term to Particle Energy
            if (IsComputeParticleEnergy) {
              particle_energy[i] += 0.5 * phi;
              particle_energy[j] += 0.5 * phi;
            }

          }
          // particle j is not contributing
          else {
            fij = 0.5 * fij + prime * meam_spline_->e_u_prime_[i];

            // Contribute pair term to Energy
            if (IsComputeEnergy) {
              *energy += 0.5 * phi;
            }

            // Contribute pair term to Particle Energy
            if (IsComputeParticleEnergy) {
              particle_energy[i] += 0.5 * phi;
            }
          }

          // Divide by r_ij to get forces from gradient
          fij /= rij;

          if (IsComputeForces) {
            forces[i][0] += dx * fij;
            forces[i][1] += dy * fij;
            forces[i][2] += dz * fij;

            forces[j][0] -= dx * fij;
            forces[j][1] -= dy * fij;
            forces[j][2] -= dz * fij;
          }  // IsComputeForces

          if (IsComputeVirial || IsComputeParticleVirial) {
            // Tabulate per-atom virial as symmetrized stress
            // tensor

            // Virial has 6 components and is stored as a
            // 6 - element vector in the following order:
            // xx, yy, zz, yz, xz, xy.

            VectorOfSizeSix v;
            v[0] = dx * dx * fij;
            v[1] = dy * dy * fij;
            v[2] = dz * dz * fij;
            v[3] = dy * dz * fij;
            v[4] = dx * dz * fij;
            v[5] = dx * dy * fij;

            if (IsComputeVirial) {
              virial[0] += v[0];
              virial[1] += v[1];
              virial[2] += v[2];
              virial[3] += v[3];
              virial[4] += v[4];
              virial[5] += v[5];
            }  // IsComputeVirial

            if (IsComputeParticleVirial) {
              v[0] *= 0.5;
              v[1] *= 0.5;
              v[2] *= 0.5;
              v[3] *= 0.5;
              v[4] *= 0.5;
              v[5] *= 0.5;

              particle_virial[i][0] += v[0];
              particle_virial[i][1] += v[1];
              particle_virial[i][2] += v[2];
              particle_virial[i][3] += v[3];
              particle_virial[i][4] += v[4];
              particle_virial[i][5] += v[5];

              particle_virial[j][0] += v[0];
              particle_virial[j][1] += v[1];
              particle_virial[j][2] += v[2];
              particle_virial[j][3] += v[3];
              particle_virial[j][4] += v[4];
              particle_virial[j][5] += v[5];
            }  // IsComputeParticleVirial
          }    // IsComputeVirial || IsComputeParticleVirial
        }      // rij2 < max_cutoff_squared_
      }        // Effective half list
    }          // Loop over j
  }            // Loop over i

  // everything is good
  return false;
}

#define KIM_LOGGER_OBJECT_NAME model_compute_arguments
template <bool IsComputeEnergy, bool IsComputeForces,
          bool IsComputeParticleEnergy, bool IsComputeVirial,
          bool IsComputeParticleVirial, bool AreKnotsOnRegularGid>
int MEAMImplementation::MeamSWSplineCompute(
    KIM::ModelCompute const *const,
    KIM::ModelComputeArguments const *const model_compute_arguments,
    int const *const, int const *const particle_contributing,
    const VectorOfSizeDIM *const coordinates, double *const energy,
    VectorOfSizeDIM *const forces, double *const particle_energy,
    VectorOfSizeSix virial, VectorOfSizeSix *const particle_virial) const {
  // Everything is good
  int ier = false;

  // Nothing to compute
  if (!IsComputeEnergy && !IsComputeForces && !IsComputeParticleEnergy &&
      !IsComputeVirial && !IsComputeParticleVirial) {
    // Everything is good
    return ier;
  }

  // Initialize energy
  if (IsComputeEnergy) {
    *energy = 0.0;
  }

  // Initialize forces
  if (IsComputeForces) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      forces[i][0] = 0.0;
      forces[i][1] = 0.0;
      forces[i][2] = 0.0;
    }
  }

  // Initialize particle energy
  if (IsComputeParticleEnergy) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      particle_energy[i] = 0.0;
    }
  }

  // Initialize virial
  if (IsComputeVirial) {
    virial[0] = 0.0;
    virial[1] = 0.0;
    virial[2] = 0.0;
    virial[3] = 0.0;
    virial[4] = 0.0;
    virial[5] = 0.0;
  }

  // Initialize particle virial
  if (IsComputeParticleVirial) {
    for (int i = 0; i < cached_number_of_particles_; ++i) {
      particle_virial[i][0] = 0.0;
      particle_virial[i][1] = 0.0;
      particle_virial[i][2] = 0.0;
      particle_virial[i][3] = 0.0;
      particle_virial[i][4] = 0.0;
      particle_virial[i][5] = 0.0;
    }
  }

  // Grow per-atom array if necessary
  meam_sw_spline_->ResizeLocalArray(cached_number_of_particles_);

  // Determine the maximum number of neighbors a single atom has
  int const max_number_of_neighbors =
      MaxNumberOfNeighbors(model_compute_arguments, particle_contributing);

  // Temporary array
  std::vector<SWTwoBody> two_body(max_number_of_neighbors);

  int number_of_neighbors;
  int const *neighbors_of_particle;

  // Temporary variable
  double prime;

  // Sum three-body contributions to charge density and
  // the embedding energy
  for (int i = 0; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    double const xi = coordinates[i][0];
    double const yi = coordinates[i][1];
    double const zi = coordinates[i][2];

    // Compute number_bonds & charge density
    int number_bonds = 0;

    double rho_value = 0.0;
    double rho_sw_value = 0.0;

    // Setup loop over neighbors of the current particle
    auto bond_j = two_body.begin();
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];

      double const dx = coordinates[j][0] - xi;
      double const dy = coordinates[j][1] - yi;
      double const dz = coordinates[j][2] - zi;

      double const rij2 = dx * dx + dy * dy + dz * dz;

      if (rij2 < max_cutoff_squared_) {
        double const rij = std::sqrt(rij2);

        bond_j->index = j;
        bond_j->r = rij;
        bond_j->f =
            meam_sw_spline_->rho_f_.Eval<AreKnotsOnRegularGid>(rij, prime);
        bond_j->f_prime = prime;
        bond_j->sw_f =
            meam_sw_spline_->esw_f_.Eval<AreKnotsOnRegularGid>(rij, prime);
        bond_j->sw_f_prime = prime;
        bond_j->dx = dx / rij;
        bond_j->dy = dy / rij;
        bond_j->dz = dz / rij;

        double partial_sum = 0.0;
        double partial_sum2 = 0.0;

        auto bond_k = two_body.begin();
        for (int kb = 0; kb < number_bonds; ++kb, ++bond_k) {
          double const cos_theta =
              (bond_j->dx * bond_k->dx + bond_j->dy * bond_k->dy +
               bond_j->dz * bond_k->dz);

          double const rho_g =
              meam_sw_spline_->rho_g_.Eval<AreKnotsOnRegularGid>(cos_theta);
          partial_sum += bond_k->f * rho_g;

          double const sw_g =
              meam_sw_spline_->esw_g_.Eval<AreKnotsOnRegularGid>(cos_theta);
          partial_sum2 += bond_k->sw_f * sw_g;
        }

        rho_value += bond_j->f * partial_sum;
        rho_sw_value += bond_j->sw_f * partial_sum2;

        double const rho_r =
            meam_sw_spline_->rho_r_.Eval<AreKnotsOnRegularGid>(rij);
        rho_value += rho_r;

        ++number_bonds;
        ++bond_j;
      }  // rij2 < max_cutoff_squared_
    }    // Loop over j

    // Compute embedding energy and its derivative
    double e_u_prime;
    double const e_u =
        meam_sw_spline_->e_u_.Eval<AreKnotsOnRegularGid>(rho_value, e_u_prime) -
        meam_sw_spline_->zero_atom_energy_;

    meam_sw_spline_->e_u_prime_[i] = e_u_prime;

    double const e_sw_prime = 1.0;
    double const e_sw = rho_sw_value;

    if (IsComputeEnergy) {
      *energy += e_u + e_sw;
    }

    if (IsComputeParticleEnergy) {
      particle_energy[i] += e_u + e_sw;
    }

    // Compute three-body contributions to force
    VectorOfSizeDIM f_i = {};

    bond_j = two_body.begin();
    for (int jb = 0; jb < number_bonds; ++jb, ++bond_j) {
      int const j = bond_j->index;
      double const r_ij = bond_j->r;
      double const f_ij = bond_j->f;
      double const f_prime_ij = bond_j->f_prime;
      double const sw_f_ij = bond_j->sw_f;
      double const sw_f_prime_ij = bond_j->sw_f_prime;

      VectorOfSizeDIM f_j = {};

      auto bond_k = two_body.begin();
      for (int kb = 0; kb < jb; ++kb, ++bond_k) {
        int const k = bond_k->index;
        double const r_ik = bond_k->r;
        double const f_ik = bond_k->f;
        double const f_prime_ik = bond_k->f_prime;
        double const sw_f_ik = bond_k->sw_f;
        double const sw_f_prime_ik = bond_k->sw_f_prime;

        double const cos_theta =
            (bond_j->dx * bond_k->dx + bond_j->dy * bond_k->dy +
             bond_j->dz * bond_k->dz);

        double const rho_g = meam_sw_spline_->rho_g_.Eval<AreKnotsOnRegularGid>(
            cos_theta, prime);
        double const prefactor = e_u_prime * f_ij * prime * f_ik;
        double const prefactor_ij = prefactor / r_ij;
        double const prefactor_ik = prefactor / r_ik;
        double const fij =
            -e_u_prime * rho_g * f_prime_ij * f_ik + prefactor_ij * cos_theta;
        double const fik =
            -e_u_prime * rho_g * f_prime_ik * f_ij + prefactor_ik * cos_theta;

        double const sw_g = meam_sw_spline_->esw_g_.Eval<AreKnotsOnRegularGid>(
            cos_theta, prime);
        double const prefactor2 = e_sw_prime * sw_f_ij * prime * sw_f_ik;
        double const prefactor2_ij = prefactor2 / r_ij;
        double const prefactor2_ik = prefactor2 / r_ik;
        double const sw_fij = -e_sw_prime * sw_g * sw_f_prime_ij * sw_f_ik +
                              prefactor2_ij * cos_theta;
        double const sw_fik = -e_sw_prime * sw_g * sw_f_prime_ik * sw_f_ij +
                              prefactor2_ik * cos_theta;

        VectorOfSizeDIM fj;

        fj[0] = bond_j->dx * fij - bond_k->dx * prefactor_ij;
        fj[1] = bond_j->dy * fij - bond_k->dy * prefactor_ij;
        fj[2] = bond_j->dz * fij - bond_k->dz * prefactor_ij;

        fj[0] += bond_j->dx * sw_fij - bond_k->dx * prefactor2_ij;
        fj[1] += bond_j->dy * sw_fij - bond_k->dy * prefactor2_ij;
        fj[2] += bond_j->dz * sw_fij - bond_k->dz * prefactor2_ij;

        f_j[0] += fj[0];
        f_j[1] += fj[1];
        f_j[2] += fj[2];

        VectorOfSizeDIM fk;

        fk[0] = bond_k->dx * fik - bond_j->dx * prefactor_ik;
        fk[1] = bond_k->dy * fik - bond_j->dy * prefactor_ik;
        fk[2] = bond_k->dz * fik - bond_j->dz * prefactor_ik;

        fk[0] += bond_k->dx * sw_fik - bond_j->dx * prefactor2_ik;
        fk[1] += bond_k->dy * sw_fik - bond_j->dy * prefactor2_ik;
        fk[2] += bond_k->dz * sw_fik - bond_j->dz * prefactor2_ik;

        f_i[0] -= fk[0];
        f_i[1] -= fk[1];
        f_i[2] -= fk[2];

        if (IsComputeForces) {
          forces[k][0] += fk[0];
          forces[k][1] += fk[1];
          forces[k][2] += fk[2];
        }

        if (IsComputeVirial || IsComputeParticleVirial) {
          // Tabulate per-atom virial as symmetrized stress
          // tensor

          // Virial has 6 components and is stored as a
          // 6 - element vector in the following order:
          // xx, yy, zz, yz, xz, xy.

          double const dx_ij = bond_j->dx * r_ij;
          double const dy_ij = bond_j->dy * r_ij;
          double const dz_ij = bond_j->dz * r_ij;

          double const dx_ik = bond_k->dx * r_ik;
          double const dy_ik = bond_k->dy * r_ik;
          double const dz_ik = bond_k->dz * r_ik;

          VectorOfSizeSix v;
          v[0] = dx_ij * fj[0] + dx_ik * fk[0];
          v[1] = dy_ij * fj[1] + dy_ik * fk[1];
          v[2] = dz_ij * fj[2] + dz_ik * fk[2];
          v[3] = dy_ij * fj[2] + dy_ik * fk[2];
          v[4] = dx_ij * fj[2] + dx_ik * fk[2];
          v[5] = dx_ij * fj[1] + dx_ik * fk[1];

          if (IsComputeVirial) {
            virial[0] += v[0];
            virial[1] += v[1];
            virial[2] += v[2];
            virial[3] += v[3];
            virial[4] += v[4];
            virial[5] += v[5];
          }  // IsComputeVirial

          if (IsComputeParticleVirial) {
            v[0] *= kThird;
            v[1] *= kThird;
            v[2] *= kThird;
            v[3] *= kThird;
            v[4] *= kThird;
            v[5] *= kThird;

            particle_virial[i][0] += v[0];
            particle_virial[i][1] += v[1];
            particle_virial[i][2] += v[2];
            particle_virial[i][3] += v[3];
            particle_virial[i][4] += v[4];
            particle_virial[i][5] += v[5];

            particle_virial[j][0] += v[0];
            particle_virial[j][1] += v[1];
            particle_virial[j][2] += v[2];
            particle_virial[j][3] += v[3];
            particle_virial[j][4] += v[4];
            particle_virial[j][5] += v[5];

            particle_virial[k][0] += v[0];
            particle_virial[k][1] += v[1];
            particle_virial[k][2] += v[2];
            particle_virial[k][3] += v[3];
            particle_virial[k][4] += v[4];
            particle_virial[k][5] += v[5];
          }  // IsComputeParticleVirial
        }    // IsComputeVirial || IsComputeParticleVirial
      }      // Loop over k bonds

      if (IsComputeForces) {
        forces[i][0] -= f_j[0];
        forces[i][1] -= f_j[1];
        forces[i][2] -= f_j[2];

        forces[j][0] += f_j[0];
        forces[j][1] += f_j[1];
        forces[j][2] += f_j[2];
      }
    }  // Loop over j bonds

    if (IsComputeForces) {
      forces[i][0] += f_i[0];
      forces[i][1] += f_i[1];
      forces[i][2] += f_i[2];
    }
  }  // Loop over i

  // Compute two-body pair interactions
  for (int i = 0; i < cached_number_of_particles_; ++i) {
    if (!particle_contributing[i]) {
      continue;
    }

    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);

    double const xi = coordinates[i][0];
    double const yi = coordinates[i][1];
    double const zi = coordinates[i][2];

    // Setup loop over neighbors of current particle
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];
      int const contributing_j = particle_contributing[j];

      // Effective half list
      if (!(contributing_j && j < i)) {
        double const dx = coordinates[j][0] - xi;
        double const dy = coordinates[j][1] - yi;
        double const dz = coordinates[j][2] - zi;

        double const rij2 = dx * dx + dy * dy + dz * dz;

        if (rij2 < max_cutoff_squared_) {
          double const rij = std::sqrt(rij2);

          double fij;
          double const phi =
              meam_sw_spline_->e_phi_.Eval<AreKnotsOnRegularGid>(rij, fij);

          meam_sw_spline_->rho_r_.Eval<AreKnotsOnRegularGid>(rij, prime);

          if (contributing_j) {
            prime *= (meam_sw_spline_->e_u_prime_[i] +
                      meam_sw_spline_->e_u_prime_[j]);
            fij += prime;

            // Contribute pair term to Energy
            if (IsComputeEnergy) {
              *energy += phi;
            }

            // Contribute pair term to Particle Energy
            if (IsComputeParticleEnergy) {
              particle_energy[i] += 0.5 * phi;
              particle_energy[j] += 0.5 * phi;
            }
          }
          // particle j is not contributing
          else {
            prime *= meam_sw_spline_->e_u_prime_[i];
            fij = 0.5 * fij + prime;

            // Contribute pair term to Energy
            if (IsComputeEnergy) {
              *energy += 0.5 * phi;
            }

            // Contribute pair term to Particle Energy
            if (IsComputeParticleEnergy) {
              particle_energy[i] += 0.5 * phi;
            }
          }

          // Divide by r_ij to get forces from gradient
          fij /= rij;

          if (IsComputeForces) {
            forces[i][0] += dx * fij;
            forces[i][1] += dy * fij;
            forces[i][2] += dz * fij;

            forces[j][0] -= dx * fij;
            forces[j][1] -= dy * fij;
            forces[j][2] -= dz * fij;
          }  // IsComputeForces

          if (IsComputeVirial || IsComputeParticleVirial) {
            // Tabulate per-atom virial as symmetrized stress
            // tensor

            // Virial has 6 components and is stored as a
            // 6 - element vector in the following order:
            // xx, yy, zz, yz, xz, xy.

            VectorOfSizeSix v;
            v[0] = dx * dx * fij;
            v[1] = dy * dy * fij;
            v[2] = dz * dz * fij;
            v[3] = dy * dz * fij;
            v[4] = dx * dz * fij;
            v[5] = dx * dy * fij;

            if (IsComputeVirial) {
              virial[0] += v[0];
              virial[1] += v[1];
              virial[2] += v[2];
              virial[3] += v[3];
              virial[4] += v[4];
              virial[5] += v[5];
            }  // IsComputeVirial

            if (IsComputeParticleVirial) {
              v[0] *= 0.5;
              v[1] *= 0.5;
              v[2] *= 0.5;
              v[3] *= 0.5;
              v[4] *= 0.5;
              v[5] *= 0.5;

              particle_virial[i][0] += v[0];
              particle_virial[i][1] += v[1];
              particle_virial[i][2] += v[2];
              particle_virial[i][3] += v[3];
              particle_virial[i][4] += v[4];
              particle_virial[i][5] += v[5];

              particle_virial[j][0] += v[0];
              particle_virial[j][1] += v[1];
              particle_virial[j][2] += v[2];
              particle_virial[j][3] += v[3];
              particle_virial[j][4] += v[4];
              particle_virial[j][5] += v[5];
            }  // IsComputeParticleVirial
          }    // IsComputeVirial || IsComputeParticleVirial
        }      // rij2 < max_cutoff_squared_
      }        // Effective half list
    }          // Loop over j
  }            // Loop over i

  // everything is good
  return false;
}

#undef KIM_LOGGER_OBJECT_NAME