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

#include "edip_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 "edip.hpp"
#include "edip_c.hpp"
#include "edip_multi.hpp"
#include "edip_single.hpp"
#include "helper.hpp"
#include "special.hpp"
#include "utility.hpp"

namespace {
static constexpr int kMaxNumberParameterFiles{2};

static constexpr int kForce{2};
static constexpr int kParticleEnergy{2};
static constexpr int kVirial{2};
static constexpr int kParticleVirial{2};

static constexpr int kMaxLineSize{1024};  // 2**10
static constexpr double kOne{1.0};        // 1
static constexpr double kThird{1.0 / 3.0};

static constexpr std::size_t kMaxNumberOfNeighbors{64};
}  // namespace

// public member functions

#ifdef KIM_LOGGER_OBJECT_NAME
#undef KIM_LOGGER_OBJECT_NAME
#endif
#define KIM_LOGGER_OBJECT_NAME model_driver_create
EDIPImplementation::EDIPImplementation(
    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) {
  // Everything is good
  *ier = false;

  if (!model_driver_create) {
    HELPER_LOG_ERROR("The model_driver_create pointer is not assigned.");
    *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("The EDIP model has 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 EDIPImplementation::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.";
      HELPER_LOG_ERROR(msg);
      for (int j = i; j--;) {
        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 EDIPImplementation::ProcessParameterFiles(
    KIM::ModelDriverCreate *const model_driver_create,
    int const number_parameter_files,
    std::FILE *const *parameter_file_pointers) {
  char next_line[kMaxLineSize];

  // create a utility object
  Utility ut;

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

  // Detect the EDIP type
  // Any of the `edip/c` or `edip/single` or `edip/multi`
  if (number_parameter_files == 1) {
    is_edip_c_ = 1;
  }

  if (is_edip_c_) {
    // Carbon EDIP
    element_name_.resize(1, "C");

    KIM::SpeciesName const species_name(element_name_[0]);

    int ier = model_driver_create->SetSpeciesCode(species_name, 0);
    if (ier) {
      LOG_ERROR("SetSpeciesCode failed to set the new species\n");
      return ier;
    }

    species_map[species_name] = 0;
  } else {
    // Read the EDIP element file

    // 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 EDIP element file.");
      return true;
    }

    // EDIP element names

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

    // # of EDIP element_name_
    if (element_name_.size() < 1) {
      HELPER_LOG_ERROR("Incorrect EDIP element file.\nNo element is provided.");
      return true;
    }

    for (std::size_t ielem = 0; ielem < element_name_.size(); ++ielem) {
      KIM::SpeciesName const species_name(element_name_[ielem]);

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

      // check for new species
      auto iter = species_map.find(species_name);
      if (iter == species_map.end()) {
        int elem = static_cast<int>(ielem);
        int ier = model_driver_create->SetSpeciesCode(species_name, elem);
        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 EDIP element file.\n";
        msg += "The Species '" + element_name_[ielem];
        msg += "' is already defined and exist.";
        HELPER_LOG_ERROR(msg);
        return true;
      }
    }  // Loop over ielem
  }

  // Detect the EDIP type
  // Any of the `edip/single` or `edip/multi`
  if (!is_edip_c_) {
    if (element_name_.size() == 1) {
      is_edip_single_ = 1;
    } else {
      is_edip_multi_ = 1;
    }
  }

  // Read EDIP parameter files.

  // edip/c
  if (is_edip_c_) {
    edip_c_.reset(new EDIPC);

    // Read EDIP parameter file. @parameter_file_pointers[0] if there is a
    // parameter file otherwise use the default parameters
    if (edip_c_->ProcessParameterFile(parameter_file_pointers[0],
                                      kMaxLineSize)) {
      HELPER_LOG_ERROR("Failed to read and parse the `edip/c` parameter file.");
      return true;
    }
  }
  // edip/single
  else if (is_edip_single_) {
    edip_single_.reset(new EDIPSingle);

    // Read EDIP parameter file. @parameter_file_pointers[1]
    if (edip_single_->ProcessParameterFile(parameter_file_pointers[1],
                                           kMaxLineSize, element_name_)) {
      HELPER_LOG_ERROR("Failed to read and parse the `edip` parameter file.");
      return true;
    }
  }
  // edip/multi
  else {
    edip_multi_.reset(new EDIPMulti);

    // Read EDIP parameter file. @parameter_file_pointers[1]
    if (edip_multi_->ProcessParameterFile(parameter_file_pointers[1],
                                          kMaxLineSize, element_name_)) {
      HELPER_LOG_ERROR("Failed to read and parse the `edip` parameter file.");
      return true;
    }
  }

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

void EDIPImplementation::CloseParameterFiles(
    int const number_parameter_files,
    std::FILE *const *parameter_file_pointers) {
  for (int i = number_parameter_files; i--;) {
    std::fclose(parameter_file_pointers[i]);
  }
}

#define KIM_LOGGER_OBJECT_NAME model_driver_create
int EDIPImplementation::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\n");
    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\n");
    return ier;
  }

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

    } else if (is_edip_single_) {
      edip_single_->ConvertUnit(convert_length, convert_energy);

    } else if (is_edip_multi_) {
      edip_multi_->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\n");
    return ier;
  }

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

template <class ModelObj>
int EDIPImplementation::SetRefreshMutableValues(ModelObj *const model_obj) {
  if (is_edip_c_) {
    edip_c_->CompleteSetup(&max_cutoff_);

    // In this model we need to request neighbors from non-contributing
    // particles
    model_will_not_request_neighbors_of_non_contributing_particles_ = 0;
    influence_distance_ = 2.0 * max_cutoff_;

  } else if (is_edip_single_) {
    edip_single_->CompleteSetup(&max_cutoff_);
    influence_distance_ = max_cutoff_;

  } else if (is_edip_multi_) {
    edip_multi_->CompleteSetup(&max_cutoff_);
    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_);

  // Set the Model's neighbor list data pointers
  model_obj->SetNeighborListPointers(
      1, &max_cutoff_,
      &model_will_not_request_neighbors_of_non_contributing_particles_);

  // Allocate local arrays which depenends on max number of possible neighbors
  if (is_edip_single_) {
    // Resize the arrays which depenends on max number of possible neighbors
    edip_single_->Resize(kMaxNumberOfNeighbors);

  } else if (is_edip_multi_) {
    edip_multi_->pre_force_coord_data_.resize(kMaxNumberOfNeighbors * 5);
  }

  // Everything is good
  return false;
}

int EDIPImplementation::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 EDIPImplementation::RegisterKIMParameters(
    KIM::ModelDriverCreate *const model_driver_create) {
  // Publish parameters
  int ier;

  if (is_edip_c_) {
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.eps, "eps", "Two-body parameter eps");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'eps'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.BB, "BB", "Two-body parameter B^4");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'BB'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.beta, "beta", "Two-body parameter beta");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'beta'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.sigma, "sigma", "Two-body parameter sigma");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'sigma'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(1, &edip_c_->params_.a, "a",
                                                   "Two-body parameter a");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'a'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.aprime, "aprime", "Two-body parameter aprime");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'aprime'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.Z0, "Z0", "Three-body parameter Z0");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Z0'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.lambda0, "lambda0",
        "Three-body parameter lambda0");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'lambda0'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.lambdaprime, "lambdaprime",
        "Three-body parameter lambdaprime");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'lambdaprime'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.gamma, "gamma", "Three-body parameter gamma");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'gamma'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(1, &edip_c_->params_.q, "q",
                                                   "Three-body parameter q");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'q'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.flow, "flow",
        "The three-parameter function lower limit");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'flow'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.fhigh, "fhigh",
        "The three-parameter function upper limit");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'fhigh'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.alpha, "alpha",
        "The three-parameter function  parameter alpha");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'alpha'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.plow, "plow",
        "The three-parameter function lower limit");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'plow'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.phigh, "phigh",
        "The three-parameter function upper limit");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'phigh'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.Zdih, "Zdih", "Coordination counting parameter");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Zdih'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.Zrep, "Zrep", "Coordination parameter");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Zrep'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_c_->params_.c0, "c0", "Repulsive cutoff limit");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'c0'\n");
      return ier;
    }
    if (edip_c_->lzbl_) {
      ier = model_driver_create->SetParameterPointer(
          1, &edip_c_->fermi_->cutoff_, "cutoff",
          "The Fermi switching cutoff distance");
      if (ier) {
        LOG_ERROR(
            "SetParameterPointer failed to set the pointer for 'cutoff'\n");
        return ier;
      }
      ier = model_driver_create->SetParameterPointer(
          1, &edip_c_->fermi_->skin_, "skin",
          "The Fermi switching skin distance (width of the switching region)");
      if (ier) {
        LOG_ERROR("SetParameterPointer failed to set the pointer for 'skin'\n");
        return ier;
      }
      ier = model_driver_create->SetParameterPointer(
          1, &edip_c_->fermi_->shift_, "shift",
          "The Fermi switching shift (delta)");
      if (ier) {
        LOG_ERROR(
            "SetParameterPointer failed to set the pointer for 'shift'\n");
        return ier;
      }
    }

  } else if (is_edip_single_) {
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.A, "A",
        "A coefficient for pair interaction I-J");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'A'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.B, "B",
        "B coefficient for pair interaction I-J");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'B'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.cutoffA, "cutoffA",
        "Cut-off distance for pair interaction I-J");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'cutoffA'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.cutoffC, "cutoffC",
        "Lower cut-off distance for calculating Z_I");
    if (ier) {
      LOG_ERROR(
          "SetParameterPointer failed to set the pointer for 'cutoffC'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.alpha, "alpha",
        "alpha coefficient for calculating Z_I");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'alpha'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.beta, "beta",
        "beta coefficient, for calculating in Gaussian function");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'beta'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.sigma, "sigma",
        "sigma coefficient for pair interaction I-J");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'sigma'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.rho, "rho",
        "rho coefficient for pair interaction I-J");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'rho'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.gamma, "gamma",
        "gamma coefficient for three-body interaction I-J-K");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'gamma'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.eta, "eta",
        "eta coefficient for function h(l,Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'eta'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.lambda, "lambda",
        "lambda coefficient for function h(l,Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'lambda'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.mu, "mu", "mu coefficient for function Q(Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'mu'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.Q0, "Q0", "Q0 coefficient for function Q(Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'Q0'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.u1, "u1",
        "u1 coefficient for function tau(Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'u1'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.u2, "u2",
        "u2 coefficient for function tau(Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'u2'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.u3, "u3",
        "u3 coefficient for function tau(Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'u3'\n");
      return ier;
    }
    ier = model_driver_create->SetParameterPointer(
        1, &edip_single_->params_.u4, "u4",
        "u4 coefficient for function tau(Z)");
    if (ier) {
      LOG_ERROR("SetParameterPointer failed to set the pointer for 'u4'\n");
      return ier;
    }

  } else if (is_edip_multi_) {
    std::string suffix;
    std::string param;
    std::size_t nelements = element_name_.size();
    for (std::size_t i = 0; i < nelements; ++i) {
      for (std::size_t j = 0; j < nelements; ++j) {
        for (std::size_t k = 0; k < nelements; ++k) {
          auto nparams = edip_multi_->elem3param_(i, j, k);
          suffix = element_name_[i] + element_name_[j] + element_name_[k];

          param = "A_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].A, param,
              "A coefficient for pair interaction I-J");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'A'\n");
            return ier;
          }
          param = "B_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].B, param,
              "B coefficient for pair interaction I-J");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'B'\n");
            return ier;
          }
          param = "cutoffA_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].cutoffA, param,
              "Cut-off distance for pair interaction I-J");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for "
                "'cutoffA'\n");
            return ier;
          }
          param = "cutoffC_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].cutoffC, param,
              "Lower cut-off distance for calculating Z_I");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for "
                "'cutoffC'\n");
            return ier;
          }
          param = "alpha_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].alpha, param,
              "alpha coefficient for calculating Z_I");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'alpha'\n");
            return ier;
          }
          param = "beta_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].beta, param,
              "beta coefficient, for calculating in Gaussian function");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'beta'\n");
            return ier;
          }
          param = "sigma_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].sigma, param,
              "sigma coefficient for pair interaction I-J");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'sigma'\n");
            return ier;
          }
          param = "rho_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].rho, param,
              "rho coefficient for pair interaction I-J");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'rho'\n");
            return ier;
          }
          param = "gamma_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].gamma, param,
              "gamma coefficient for three-body interaction I-J-K");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'gamma'\n");
            return ier;
          }
          param = "eta_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].eta, param,
              "eta coefficient for function h(l,Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'eta'\n");
            return ier;
          }
          param = "lambda_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].lambda, param,
              "lambda coefficient for function h(l,Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'lambda'\n");
            return ier;
          }
          param = "mu_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].mu, param,
              "mu coefficient for function Q(Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'mu'\n");
            return ier;
          }
          param = "Q0_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].Q0, param,
              "Q0 coefficient for function Q(Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'Q0'\n");
            return ier;
          }
          param = "u1_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].u1, param,
              "u1 coefficient for function tau(Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'u1'\n");
            return ier;
          }
          param = "u2_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].u2, param,
              "u2 coefficient for function tau(Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'u2'\n");
            return ier;
          }
          param = "u3_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].u3, param,
              "u3 coefficient for function tau(Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'u3'\n");
            return ier;
          }
          param = "u4_" + suffix;
          ier = model_driver_create->SetParameterPointer(
              1, &edip_multi_->params_[nparams].u4, param,
              "u4 coefficient for function tau(Z)");
          if (ier) {
            LOG_ERROR(
                "SetParameterPointer failed to set the pointer for 'u4'\n");
            return ier;
          }
        }  // Loop over k
      }    // Loop over j
    }      // Loop over i
  }

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

int EDIPImplementation::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 *>(EDIP::Destroy)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::Refresh, KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(EDIP::Refresh)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::WriteParameterizedModel,
          KIM::LANGUAGE_NAME::cpp, false,
          reinterpret_cast<KIM::Function *>(EDIP::WriteParameterizedModel)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::Compute, KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(EDIP::Compute)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::ComputeArgumentsCreate,
          KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(EDIP::ComputeArgumentsCreate)) ||
      model_driver_create->SetRoutinePointer(
          KIM::MODEL_ROUTINE_NAME::ComputeArgumentsDestroy,
          KIM::LANGUAGE_NAME::cpp, true,
          reinterpret_cast<KIM::Function *>(EDIP::ComputeArgumentsDestroy));
  return ier;
}

#define KIM_LOGGER_OBJECT_NAME model_compute_arguments
int EDIPImplementation::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 returns an error\n");
    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;
  number_of_particles_ = static_cast<std::size_t>(*number_of_particles);

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

constexpr int EDIPImplementation::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) {
  return (((static_cast<int>(is_compute_energy) * kForce +
            static_cast<int>(is_compute_forces)) *
               kParticleEnergy +
           static_cast<int>(is_compute_particle_energy)) *
              kVirial +
          static_cast<int>(is_compute_virial)) *
             kParticleVirial +
         static_cast<int>(is_compute_particle_virial);
}

#define KIM_LOGGER_OBJECT_NAME model_compute_arguments_create
int EDIPImplementation::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 EDIPImplementation::Refresh(KIM::ModelRefresh *const model_refresh) {
  return SetRefreshMutableValues(model_refresh);
}

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

  std::string const *path{nullptr};
  model_write_parameterized_model->GetPath(&path);

  std::string const *model_name{nullptr};
  model_write_parameterized_model->GetModelName(&model_name);

  if (is_edip_c_) {
    // EDIP/C parameter file
    std::string parameter_file;
    for (auto element_name : element_name_) {
      parameter_file += element_name;
    }
    parameter_file += ".parameter.edip";

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

    // EDIP/C parameter 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.");
      return true;
    } else {
      if (fs.fail()) {
        HELPER_LOG_ERROR(
            "An error has occurred from opening the file on the associated "
            "stream!");
        return true;
      }

      fs << "# EDIP/C parameter file\n\n";

      auto sigma = edip_c_->params_.sigma;
      auto a = edip_c_->params_.a;
      auto aprime = edip_c_->params_.aprime;
      auto gamma = edip_c_->params_.gamma;
      auto flow = edip_c_->params_.flow;
      auto fhigh = edip_c_->params_.fhigh;
      auto plow = edip_c_->params_.plow;
      auto phigh = edip_c_->params_.phigh;
      auto c0 = edip_c_->params_.c0;
      auto BB = edip_c_->params_.BB;
      auto eps = edip_c_->params_.eps;
      auto lambda0 = edip_c_->params_.lambda0;
      auto Zdih = edip_c_->params_.Zdih;
      auto Zrep = edip_c_->params_.Zrep;
      double cutoff{0.0};
      double skin{0.0};
      double shift{0.0};
      if (edip_c_->lzbl_) {
        cutoff = edip_c_->fermi_->cutoff_;
        skin = edip_c_->fermi_->skin_;
        shift = edip_c_->fermi_->shift_;
      }

      // (Angstroms in metal units).
      if (special::IsNotOne(edip_c_->convert_length_factor_)) {
        double const convert_length_factor_inverse =
            1.0 / edip_c_->convert_length_factor_;
        sigma *= convert_length_factor_inverse;
        a *= convert_length_factor_inverse;
        aprime *= convert_length_factor_inverse;
        gamma *= convert_length_factor_inverse;
        flow *= convert_length_factor_inverse;
        fhigh *= convert_length_factor_inverse;
        plow *= convert_length_factor_inverse;
        phigh *= convert_length_factor_inverse;
        c0 *= convert_length_factor_inverse;

        auto cst2 =
            convert_length_factor_inverse * convert_length_factor_inverse;
        auto cst4 = cst2 * cst2;
        BB *= cst4;

        cutoff *= convert_length_factor_inverse;
        skin *= convert_length_factor_inverse;
        shift *= convert_length_factor_inverse;
      }

      //  (eV in the case of metal units).
      if (special::IsNotOne(edip_c_->convert_energy_factor_)) {
        double const convert_energy_factor_inverse =
            1.0 / edip_c_->convert_energy_factor_;

        eps *= convert_energy_factor_inverse;
        lambda0 *= convert_energy_factor_inverse;
        Zdih *= convert_energy_factor_inverse;
        Zrep *= convert_energy_factor_inverse;
      }

      fs << eps << "  " << BB << "  " << edip_c_->params_.beta << "  " << sigma
         << "  " << a << "  " << aprime << "\n"
         << edip_c_->params_.Z0 << "  " << lambda0 << "  "
         << edip_c_->params_.lambdaprime << "  " << gamma << "  "
         << edip_c_->params_.q << "\n"
         << flow << "  " << fhigh << "  " << edip_c_->params_.alpha << "  "
         << plow << "  " << phigh << Zdih << "  " << Zrep << "  " << c0 << "  ";
      if (edip_c_->lzbl_) {
        fs << edip_c_->lzbl_ << "  " << cutoff << " " << skin << "  " << shift;
      } else {
        fs << edip_c_->lzbl_;
      }
      fs << "\n\n\n\n#\n# EDIP/C parameter file\n#\n";
      fs.close();
    }

  } else {
    // EDIP elements file
    {
      std::string elements_file(*model_name);

      elements_file += ".elements.edip";

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

      // EDIP 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.");
        return true;
      } else {
        if (fs.fail()) {
          std::string msg = "An error has occurred from opening the ";
          msg += "file on the associated stream!";
          HELPER_LOG_ERROR(msg);
          return true;
        }
        fs << "# EDIP elements \n";
        for (auto element_name : element_name_) {
          fs << element_name << "  ";
        }
        fs << "\n\n\n\n#\n# EDIP elements file\n#\n";
        fs.close();
      }
    }

    // EDIP parameter file
    {
      std::string parameter_file;
      for (auto element_name : element_name_) {
        parameter_file += element_name;
      }
      parameter_file += ".parameter.edip";

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

      // EDIP parameter 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.");
        return true;
      } else {
        if (fs.fail()) {
          std::string msg = "An error has occurred from opening the ";
          msg += "file on the associated stream!";
          HELPER_LOG_ERROR(msg);
          return true;
        }

        fs << "# EDIP parameter file";

        if (is_edip_single_) {
          auto B = edip_single_->params_.B;
          auto cutoffA = edip_single_->params_.cutoffA;
          auto cutoffC = edip_single_->params_.cutoffC;
          auto gamma = edip_single_->params_.gamma;
          auto sigma = edip_single_->params_.sigma;

          auto A = edip_single_->params_.A;
          auto lambda = edip_single_->params_.lambda;

          if (special::IsNotOne(edip_single_->convert_length_factor_)) {
            double const convert_length_factor_inverse =
                kOne / edip_single_->convert_length_factor_;
            B *= convert_length_factor_inverse;
            cutoffA *= convert_length_factor_inverse;
            cutoffC *= convert_length_factor_inverse;
            gamma *= convert_length_factor_inverse;
            sigma *= convert_length_factor_inverse;
          }

          if (special::IsNotOne(edip_single_->convert_energy_factor_)) {
            double const convert_energy_factor_inverse =
                kOne / edip_single_->convert_energy_factor_;
            A *= convert_energy_factor_inverse;
            lambda *= convert_energy_factor_inverse;
          }

          fs << "\n\n"
             << element_name_[0] << "  " << element_name_[0] << "  "
             << element_name_[0] << "\n"
             << A << "  " << B << "  " << cutoffA << "  " << cutoffC << "  "
             << edip_single_->params_.alpha << "  "
             << edip_single_->params_.beta << "  " << edip_single_->params_.eta
             << "\n"
             << gamma << "  " << lambda << "  " << edip_single_->params_.mu
             << "  " << edip_single_->params_.rho << "  " << sigma << "  "
             << edip_single_->params_.Q0 << "\n"
             << edip_single_->params_.u1 << "  " << edip_single_->params_.u2
             << "  " << edip_single_->params_.u3 << "  "
             << edip_single_->params_.u4;

        } else if (is_edip_multi_) {
          double const convert_length_factor{
              edip_multi_->convert_length_factor_};
          double const convert_length_factor_inverse =
              kOne / convert_length_factor;
          double const convert_energy_factor{
              edip_multi_->convert_energy_factor_};
          double const convert_energy_factor_inverse =
              kOne / convert_energy_factor;

          auto nelements = element_name_.size();
          for (std::size_t i = 0; i < nelements; ++i) {
            for (std::size_t j = 0; j < nelements; ++j) {
              for (std::size_t k = 0; k < nelements; ++k) {
                auto const nparams = edip_multi_->elem3param_(i, j, k);

                auto B = edip_multi_->params_[nparams].B;
                auto cutoffA = edip_multi_->params_[nparams].cutoffA;
                auto cutoffC = edip_multi_->params_[nparams].cutoffC;
                auto gamma = edip_multi_->params_[nparams].gamma;
                auto sigma = edip_multi_->params_[nparams].sigma;

                auto A = edip_multi_->params_[nparams].A;
                auto lambda = edip_multi_->params_[nparams].lambda;

                if (special::IsNotOne(convert_length_factor)) {
                  B *= convert_length_factor_inverse;
                  cutoffA *= convert_length_factor_inverse;
                  cutoffC *= convert_length_factor_inverse;
                  gamma *= convert_length_factor_inverse;
                  sigma *= convert_length_factor_inverse;
                }

                if (special::IsNotOne(convert_energy_factor)) {
                  A *= convert_energy_factor_inverse;
                  lambda *= convert_energy_factor_inverse;
                }

                fs << "\n\n"
                   << element_name_[i] << "  " << element_name_[j] << "  "
                   << element_name_[k] << "\n"
                   << A << "  " << B << "  " << cutoffA << "  " << cutoffC
                   << "  " << edip_multi_->params_[nparams].alpha << "  "
                   << edip_multi_->params_[nparams].beta << "  "
                   << edip_multi_->params_[nparams].eta << "\n"
                   << gamma << "  " << lambda << "  "
                   << edip_multi_->params_[nparams].mu << "  "
                   << edip_multi_->params_[nparams].rho << "  " << sigma << "  "
                   << edip_multi_->params_[nparams].Q0 << "\n"
                   << edip_multi_->params_[nparams].u1 << "  "
                   << edip_multi_->params_[nparams].u2 << "  "
                   << edip_multi_->params_[nparams].u3 << "  "
                   << edip_multi_->params_[nparams].u4;
              }  // Loop over k
            }    // Loop over j
          }      // Loop over i
        }        // is_edip_single_ || is_edip_multi_
      }          // fs.is_open()

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

  // Everything is good
  return false;
}

int EDIPImplementation::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!");
    return ier;
  }

#include "edip_implementation_compute_dispatch.cpp"

  return ier;
}

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

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

std::size_t EDIPImplementation::TotalNumberOfNeighbors(
    KIM::ModelComputeArguments const *const model_compute_arguments) 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) {
    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors,
                                             &neighbors_of_particle);
    total_number_of_neighbors += static_cast<std::size_t>(number_of_neighbors);
  }  // Loop over i
  return total_number_of_neighbors;
}

void EDIPImplementation::SetDistanceDirections(
    KIM::ModelComputeArguments const *const model_compute_arguments,
    const VectorOfSizeDIM *const coordinates) const {
  auto number_of_neighbors = &edip_c_->number_of_neighbors_[0];
  auto neighbors_offset = &edip_c_->neighbors_offset_[0];
  auto neigh_to_part = &edip_c_->neighbors_to_particles_[0];
  auto r = &edip_c_->r_[0];
  auto N = &edip_c_->N_[0];

  int const *neighbors_of_particle;
  int number_of_neighbors_i;

  std::size_t nn{0};
  std::size_t offset{0};

  for (int i = 0; i < cached_number_of_particles_; ++i) {
    // Get the neighbors
    model_compute_arguments->GetNeighborList(0, i, &number_of_neighbors_i,
                                             &neighbors_of_particle);

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

    // update offset
    *neighbors_offset++ = offset;

    nn = 0;
    for (int jn = 0; jn < number_of_neighbors_i; ++jn) {
      int const j = neighbors_of_particle[jn];

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

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

      if (rij2 < max_cutoff_squared_) {
        *neigh_to_part++ = static_cast<std::size_t>(j);
        double const rij = std::sqrt(rij2);
        *r++ = rij;
        *N++ = dx / rij;
        *N++ = dy / rij;
        *N++ = dz / rij;
        ++nn;
      }
    }  // Loop over jn neighbors

    offset += nn;
    *number_of_neighbors++ = nn;
  }  // Loop over i
}

// Main compute function

template <bool IsComputeEnergy, bool IsComputeForces,
          bool IsComputeParticleEnergy, bool IsComputeVirial,
          bool IsComputeParticleVirial>
int EDIPImplementation::EDIPCCompute(
    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;
    }
  }

  std::size_t const total_number_of_neighbors =
      TotalNumberOfNeighbors(model_compute_arguments);

  // Allocate memory
  edip_c_->Resize(number_of_particles_, total_number_of_neighbors);

  // Compute correct number of neighbors & distances and their directions
  SetDistanceDirections(model_compute_arguments, coordinates);

  // Complete the memory allocation
  edip_c_->Resize(number_of_particles_);

  // Calculate spherical coordination z_i for all particles.
  edip_c_->CutoffFunction(number_of_particles_);

  // Calculate switching functions for all particles.
  edip_c_->SwitchFunction(number_of_particles_);

  if (edip_c_->lzbl_) {
    for (std::size_t i = 0; i < number_of_particles_; ++i) {
      if (!particle_contributing[i]) {
        continue;
      }

      edip_c_->Coordination(i);

      edip_c_
          ->PairZbl<IsComputeEnergy, IsComputeForces, IsComputeParticleEnergy,
                    IsComputeVirial, IsComputeParticleVirial>(
              i, particle_contributing, energy, forces, particle_energy, virial,
              particle_virial);

      edip_c_->Triple<IsComputeEnergy, IsComputeForces, IsComputeParticleEnergy,
                      IsComputeVirial, IsComputeParticleVirial>(
          i, particle_contributing, energy, forces, particle_energy, virial,
          particle_virial);
    }  // Loop over i
  } else {
    for (std::size_t i = 0; i < number_of_particles_; ++i) {
      if (!particle_contributing[i]) {
        continue;
      }

      edip_c_->Coordination(i);

      edip_c_->Pair<IsComputeEnergy, IsComputeForces, IsComputeParticleEnergy,
                    IsComputeVirial, IsComputeParticleVirial>(
          i, particle_contributing, energy, forces, particle_energy, virial,
          particle_virial);

      edip_c_->Triple<IsComputeEnergy, IsComputeForces, IsComputeParticleEnergy,
                      IsComputeVirial, IsComputeParticleVirial>(
          i, particle_contributing, energy, forces, particle_energy, virial,
          particle_virial);
    }  // Loop over i
  }    // lzbl_

  // everything is good
  return false;
}

template <bool IsComputeEnergy, bool IsComputeForces,
          bool IsComputeParticleEnergy, bool IsComputeVirial,
          bool IsComputeParticleVirial>
int EDIPImplementation::EDIPSingleCompute(
    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;
    }
  }

  double const cutoffA{edip_single_->params_.cutoffA};
  double const cutoffC{edip_single_->params_.cutoffC};
  double const sigma{edip_single_->params_.sigma};
  double const gamma{edip_single_->params_.gamma};
  double const cutsq{edip_single_->params_.cutsq};
  double const beta{edip_single_->params_.beta};
  double const rho{edip_single_->params_.rho};
  double const lambda{edip_single_->params_.lambda};
  double const eta{edip_single_->params_.eta};
  double const mu{edip_single_->params_.mu};

  double const GridStart{edip_single_->GridStart_};

  auto exp3B = &edip_single_->exp3B_[0];
  auto exp2B = &edip_single_->exp2B_[0];
  auto pow2B = &edip_single_->pow2B_[0];
  auto cutoffFunction = &edip_single_->cutoffFunction_[0];
  auto cutoffFunctionDerived = &edip_single_->cutoffFunctionDerived_[0];
  auto expMinusBetaZeta_iZeta_iGrid =
      &edip_single_->expMinusBetaZeta_iZeta_iGrid_[0];
  auto qFunctionGrid = &edip_single_->qFunctionGrid_[0];
  auto tauFunctionGrid = &edip_single_->tauFunctionGrid_[0];
  auto tauFunctionDerivedGrid = &edip_single_->tauFunctionDerivedGrid_[0];

  auto pre_force_rij_inverse = &edip_single_->pre_force_rij_inverse_[0];
  auto pre_force_Exp3B_ij = &edip_single_->pre_force_Exp3B_ij_[0];
  auto pre_force_Exp3BDerived_ij = &edip_single_->pre_force_Exp3BDerived_ij_[0];
  auto pre_force_Exp2B_ij = &edip_single_->pre_force_Exp2B_ij_[0];
  auto pre_force_Exp2BDerived_ij = &edip_single_->pre_force_Exp2BDerived_ij_[0];
  auto pre_force_Pow2B_ij = &edip_single_->pre_force_Pow2B_ij_[0];
  auto pre_force_coord_data = &edip_single_->pre_force_coord_data_[0];

  int number_of_neighbors;
  int const *neighbors_of_particle;
  double y1;
  double ydiff;

  // Loop over full neighbor list of my atoms
  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];

    double zeta_i{0.0};
    int number_pre_force_coord_pairs{0};

    // Setup loop over neighbors of current particle
    double *pre_force_coord = &pre_force_coord_data[0];

    // pre-loop to compute environment coordination f(Z)
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];

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

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

      if (rij2 > cutsq) {
        continue;
      }

      double const rij = std::sqrt(rij2);
      double const rij_inverse = 1.0 / rij;
      pre_force_rij_inverse[jn] = rij_inverse;

      double const rij_ca = rij - cutoffA;
      double const rij_ca_inverse = 1.0 / rij_ca;
      double const sigma_rij_ca_inverse = sigma * rij_ca_inverse;
      double const gamma_rij_ca_inverse = gamma * rij_ca_inverse;

      double const deltax_start = rij > GridStart ? rij - GridStart : 0.0;
      double const deltax_grid_tmp = deltax_start * edip_single::kGridDensity;
      std::size_t const index = static_cast<std::size_t>(deltax_grid_tmp);
      double const deltax = deltax_grid_tmp - static_cast<double>(index);

      y1 = exp3B[index];
      ydiff = exp3B[index + 1] - y1;

      double const exp3B_ij = y1 + ydiff * deltax;
      double const exp3BDerived_ij =
          -exp3B_ij * gamma_rij_ca_inverse * rij_ca_inverse;

      pre_force_Exp3B_ij[jn] = exp3B_ij;
      pre_force_Exp3BDerived_ij[jn] = exp3BDerived_ij;

      y1 = exp2B[index];
      ydiff = exp2B[index + 1] - y1;

      double const exp2B_ij = y1 + ydiff * deltax;
      double const exp2BDerived_ij =
          -exp2B_ij * sigma_rij_ca_inverse * rij_ca_inverse;

      pre_force_Exp2B_ij[jn] = exp2B_ij;
      pre_force_Exp2BDerived_ij[jn] = exp2BDerived_ij;

      y1 = pow2B[index];
      ydiff = pow2B[index + 1] - y1;

      double const pow2B_ij = y1 + ydiff * deltax;

      pre_force_Pow2B_ij[jn] = pow2B_ij;

      // zeta and its derivative

      double zeta_iDerivedInvR_ij{};

      if (rij < cutoffC) {
        zeta_i += 1.0;
      } else {
        y1 = cutoffFunction[index];
        ydiff = cutoffFunction[index + 1] - y1;

        double const cutoffFunction_ij = y1 + ydiff * deltax;

        zeta_i += cutoffFunction_ij;

        y1 = cutoffFunctionDerived[index];
        ydiff = cutoffFunctionDerived[index + 1] - y1;

        double const zeta_iDerived = y1 + ydiff * deltax;
        zeta_iDerivedInvR_ij = zeta_iDerived * rij_inverse;

        pre_force_coord[0] = zeta_iDerivedInvR_ij;
        pre_force_coord[1] = dx;
        pre_force_coord[2] = dy;
        pre_force_coord[3] = dz;
        pre_force_coord[4] = static_cast<double>(j);

        pre_force_coord += 5;
        ++number_pre_force_coord_pairs;
      }  // rij cutoffC
    }    // Loop over jn

    // quantities depending on zeta_i

    double const interpolTMP = (zeta_i * edip_single::kGridDensity);
    std::size_t const interpolIDX = static_cast<std::size_t>(interpolTMP);
    double const interpolTMPIDX =
        interpolTMP - static_cast<double>(interpolIDX);

    y1 = expMinusBetaZeta_iZeta_iGrid[interpolIDX];
    ydiff = expMinusBetaZeta_iZeta_iGrid[interpolIDX + 1] - y1;
    double const expMinusBetaZeta_iZeta_i = y1 + ydiff * interpolTMPIDX;

    y1 = qFunctionGrid[interpolIDX];
    ydiff = qFunctionGrid[interpolIDX + 1] - y1;
    double const qFunction = y1 + ydiff * interpolTMPIDX;

    y1 = tauFunctionGrid[interpolIDX];
    ydiff = tauFunctionGrid[interpolIDX + 1] - y1;
    double const tauFunction = y1 + ydiff * interpolTMPIDX;

    y1 = tauFunctionDerivedGrid[interpolIDX];
    ydiff = tauFunctionDerivedGrid[interpolIDX + 1] - y1;
    double const tauFunctionDerived = y1 + ydiff * interpolTMPIDX;

    double const forceModCoord_factor =
        2.0 * beta * zeta_i * expMinusBetaZeta_iZeta_i;

    double forceModCoord{0.0};

    // two-body interactions, skip half of them
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];

      VectorOfSizeDIM dr_ij;
      dr_ij[0] = coordinates[j][0] - xi;
      dr_ij[1] = coordinates[j][1] - yi;
      dr_ij[2] = coordinates[j][2] - zi;

      double const rij2 =
          dr_ij[0] * dr_ij[0] + dr_ij[1] * dr_ij[1] + dr_ij[2] * dr_ij[2];

      if (rij2 > cutsq) {
        continue;
      }

      int const contributing_j =
          IsComputeParticleEnergy ? particle_contributing[j] : false;

      double const rij_inverse = pre_force_rij_inverse[jn];

      double const pow2B_ij = pre_force_Pow2B_ij[jn];
      double const exp2BDerived_ij = pre_force_Exp2BDerived_ij[jn];

      double const exp2B_ij = pre_force_Exp2B_ij[jn];
      double const pow2BDerived_ij = -rho * rij_inverse * pow2B_ij;

      double const potential2B_factor = pow2B_ij - expMinusBetaZeta_iZeta_i;
      forceModCoord += forceModCoord_factor * exp2B_ij;

      double const forceMod2B =
          exp2BDerived_ij * potential2B_factor + exp2B_ij * pow2BDerived_ij;

      double const directorCos_ij_x = rij_inverse * dr_ij[0];
      double const directorCos_ij_y = rij_inverse * dr_ij[1];
      double const directorCos_ij_z = rij_inverse * dr_ij[2];

      double const exp3B_ij = pre_force_Exp3B_ij[jn];
      double const exp3BDerived_ij = pre_force_Exp3BDerived_ij[jn];

      if (IsComputeEnergy || IsComputeParticleEnergy) {
        double const phi = exp2B_ij * potential2B_factor;

        // Contribution to energy
        if (IsComputeEnergy) {
          *energy += phi;
        }  // IsComputeEnergy

        // Contribution to particle_energy
        if (IsComputeParticleEnergy) {
          if (contributing_j) {
            double const phihalf = 0.5 * phi;
            particle_energy[i] += phihalf;
            particle_energy[j] += phihalf;
          } else {
            particle_energy[i] += phi;
          }
        }  // IsComputeParticleEnergy
      }    // IsComputeEnergy || IsComputeParticleEnergy

      if (IsComputeForces || IsComputeVirial || IsComputeParticleVirial) {
        VectorOfSizeDIM f_ij;
        f_ij[0] = forceMod2B * directorCos_ij_x;
        f_ij[1] = forceMod2B * directorCos_ij_y;
        f_ij[2] = forceMod2B * directorCos_ij_z;

        if (IsComputeForces) {
          forces[i][0] += f_ij[0];
          forces[i][1] += f_ij[1];
          forces[i][2] += f_ij[2];

          forces[j][0] -= f_ij[0];
          forces[j][1] -= f_ij[1];
          forces[j][2] -= f_ij[2];
        }  // IsComputeForces

        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] = -dr_ij[0] * f_ij[0];
          v[1] = -dr_ij[1] * f_ij[1];
          v[2] = -dr_ij[2] * f_ij[2];
          v[3] = -dr_ij[1] * f_ij[2];
          v[4] = -dr_ij[0] * f_ij[2];
          v[5] = -dr_ij[0] * f_ij[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] *= 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
      }      // IsComputeForces || IsComputeVirial || IsComputeParticleVirial

      // three-body Forces
      for (int kn = jn + 1; kn < number_of_neighbors; ++kn) {
        int const k = neighbors_of_particle[kn];

        VectorOfSizeDIM dr_ik;
        dr_ik[0] = coordinates[k][0] - xi;
        dr_ik[1] = coordinates[k][1] - yi;
        dr_ik[2] = coordinates[k][2] - zi;

        double const rik2 =
            dr_ik[0] * dr_ik[0] + dr_ik[1] * dr_ik[1] + dr_ik[2] * dr_ik[2];

        if (rik2 > cutsq) {
          continue;
        }

        double const rik_inverse = pre_force_rij_inverse[kn];

        double const directorCos_ik_x = rik_inverse * dr_ik[0];
        double const directorCos_ik_y = rik_inverse * dr_ik[1];
        double const directorCos_ik_z = rik_inverse * dr_ik[2];

        double const cosTeta = directorCos_ij_x * directorCos_ik_x +
                               directorCos_ij_y * directorCos_ik_y +
                               directorCos_ij_z * directorCos_ik_z;

        double const cosTetaDiff = cosTeta + tauFunction;
        double const cosTetaDiffCosTetaDiff = cosTetaDiff * cosTetaDiff;

        double const qFunctionCosTetaDiffCosTetaDiff =
            cosTetaDiffCosTetaDiff * qFunction;
        double const expMinusQFunctionCosTetaDiffCosTetaDiff =
            std::exp(-qFunctionCosTetaDiffCosTetaDiff);

        double const potentia3B_factor =
            lambda * ((1.0 - expMinusQFunctionCosTetaDiffCosTetaDiff) +
                      eta * qFunctionCosTetaDiffCosTetaDiff);

        double const exp3B_ik = pre_force_Exp3B_ij[kn];
        double const exp3BDerived_ik = pre_force_Exp3BDerived_ij[kn];

        double const forceMod3B_factor1_ij =
            -exp3BDerived_ij * exp3B_ik * potentia3B_factor;
        double const forceMod3B_factor2 =
            2.0 * lambda * exp3B_ij * exp3B_ik * qFunction * cosTetaDiff *
            (eta + expMinusQFunctionCosTetaDiffCosTetaDiff);
        double const forceMod3B_factor2_ij = forceMod3B_factor2 * rij_inverse;

        if (IsComputeEnergy || IsComputeParticleEnergy) {
          // potential energy
          double const phi = exp3B_ij * exp3B_ik * potentia3B_factor;

          // Contribution to energy
          if (IsComputeEnergy) {
            *energy += phi;
          }  // IsComputeEnergy

          // Contribution to particle_energy
          if (IsComputeParticleEnergy) {
            double const phithird = kThird * phi;
            particle_energy[i] += phithird;
            if (contributing_j) {
              particle_energy[j] += phithird;
            } else {
              particle_energy[i] += phithird;
            }
            if (particle_contributing[k]) {
              particle_energy[k] += phithird;
            } else {
              particle_energy[i] += phithird;
            }
          }  // IsComputeParticleEnergy
        }    // IsComputeEnergy || IsComputeParticleEnergy

        if (IsComputeForces || IsComputeVirial || IsComputeParticleVirial) {
          VectorOfSizeDIM f_ij;
          f_ij[0] = forceMod3B_factor1_ij * directorCos_ij_x +
                    forceMod3B_factor2_ij *
                        (cosTeta * directorCos_ij_x - directorCos_ik_x);
          f_ij[1] = forceMod3B_factor1_ij * directorCos_ij_y +
                    forceMod3B_factor2_ij *
                        (cosTeta * directorCos_ij_y - directorCos_ik_y);
          f_ij[2] = forceMod3B_factor1_ij * directorCos_ij_z +
                    forceMod3B_factor2_ij *
                        (cosTeta * directorCos_ij_z - directorCos_ik_z);

          double const forceMod3B_factor1_ik =
              -exp3BDerived_ik * exp3B_ij * potentia3B_factor;
          double const forceMod3B_factor2_ik = forceMod3B_factor2 * rik_inverse;

          VectorOfSizeDIM f_ik;
          f_ik[0] = forceMod3B_factor1_ik * directorCos_ik_x +
                    forceMod3B_factor2_ik *
                        (cosTeta * directorCos_ik_x - directorCos_ij_x);
          f_ik[1] = forceMod3B_factor1_ik * directorCos_ik_y +
                    forceMod3B_factor2_ik *
                        (cosTeta * directorCos_ik_y - directorCos_ij_y);
          f_ik[2] = forceMod3B_factor1_ik * directorCos_ik_z +
                    forceMod3B_factor2_ik *
                        (cosTeta * directorCos_ik_z - directorCos_ij_z);

          forceModCoord += (forceMod3B_factor2 *
                            (tauFunctionDerived - 0.5 * mu * cosTetaDiff));

          if (IsComputeForces) {
            forces[j][0] += f_ij[0];
            forces[j][1] += f_ij[1];
            forces[j][2] += f_ij[2];

            forces[k][0] += f_ik[0];
            forces[k][1] += f_ik[1];
            forces[k][2] += f_ik[2];

            forces[i][0] -= f_ij[0] + f_ik[0];
            forces[i][1] -= f_ij[1] + f_ik[1];
            forces[i][2] -= f_ij[2] + f_ik[2];
          }  // IsComputeForces

          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] = dr_ij[0] * f_ij[0] + dr_ik[0] * f_ik[0];
            v[1] = dr_ij[1] * f_ij[1] + dr_ik[1] * f_ik[1];
            v[2] = dr_ij[2] * f_ij[2] + dr_ik[2] * f_ik[2];
            v[3] = dr_ij[1] * f_ij[2] + dr_ik[1] * f_ik[2];
            v[4] = dr_ij[0] * f_ij[2] + dr_ik[0] * f_ik[2];
            v[5] = dr_ij[0] * f_ij[1] + dr_ik[0] * f_ik[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
        }      // IsComputeForces || IsComputeVirial || IsComputeParticleVirial
      }        // Loop over kn
    }          // Loop over jn

    // forces due to environment coordination f(Z)
    if (IsComputeForces || IsComputeVirial || IsComputeParticleVirial) {
      pre_force_coord = &pre_force_coord_data[0];

      for (int idx = 0; idx < number_pre_force_coord_pairs; ++idx) {
        double const dzetair = pre_force_coord[0] * forceModCoord;
        double const dx = pre_force_coord[1];
        double const dy = pre_force_coord[2];
        double const dz = pre_force_coord[3];
        int const j = static_cast<int>(pre_force_coord[4]);

        pre_force_coord += 5;

        double const dx_dzetair = dx * dzetair;
        double const dy_dzetair = dy * dzetair;
        double const dz_dzetair = dz * dzetair;

        // Contribution to forces
        if (IsComputeForces) {
          forces[i][0] -= dx_dzetair;
          forces[i][1] -= dy_dzetair;
          forces[i][2] -= dz_dzetair;

          forces[j][0] += dx_dzetair;
          forces[j][1] += dy_dzetair;
          forces[j][2] += dz_dzetair;
        }  //

        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 * dx_dzetair;
          v[1] = -dy * dy_dzetair;
          v[2] = -dz * dz_dzetair;
          v[3] = -dz * dy_dzetair;
          v[4] = -dz * dx_dzetair;
          v[5] = -dx * dy_dzetair;

          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
      }      // Loop over idx
    }        // IsComputeForces || IsComputeVirial || IsComputeParticleVirial
  }          // Loop over i

  // everything is good
  return false;
}

template <bool IsComputeEnergy, bool IsComputeForces,
          bool IsComputeParticleEnergy, bool IsComputeVirial,
          bool IsComputeParticleVirial>
int EDIPImplementation::EDIPMultiCompute(
    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;
    }
  }

  int number_of_neighbors;
  int const *neighbors_of_particle;

  // Loop over full neighbor list of my atoms
  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);

    std::size_t const species_i =
        static_cast<std::size_t>(particle_species_codes[i]);

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

    double zeta_i{0.0};
    int number_pre_force_coord_pairs{0};

    // Setup loop over neighbors of current particle
    double *pre_force_coord = &edip_multi_->pre_force_coord_data_[0];

    // pre-loop to compute environment coordination f(Z)
    for (int jn = 0; jn < number_of_neighbors; ++jn) {
      int const j = neighbors_of_particle[jn];
      std::size_t const species_j =
          static_cast<std::size_t>(particle_species_codes[j]);

      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;

      std::size_t const ijparam =
          edip_multi_->elem3param_(species_i, species_j, species_j);

      if (rij2 > edip_multi_->params_[ijparam].cutsq) {
        continue;
      }

      double const rij = std::sqrt(rij2);

      // zeta and its derivative dZ/dr
      if (rij < edip_multi_->params_[ijparam].cutoffC) {
        zeta_i += 1.0;
      } else {
        double fdr;
        double const f = edip_multi_->EdipFc(rij, ijparam, fdr);

        zeta_i += f;

        pre_force_coord[0] = -fdr / rij;
        pre_force_coord[1] = dx;
        pre_force_coord[2] = dy;
        pre_force_coord[3] = dz;
        pre_force_coord[4] = static_cast<double>(j);

        pre_force_coord += 5;
        ++number_pre_force_coord_pairs;
      }  // rij cutoffC
    }    // Loop over jn

    // two-body interactions
    double dpairZ{0.0};
    double dtripleZ{0.0};

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

      std::size_t const species_j =
          static_cast<std::size_t>(particle_species_codes[j]);

      VectorOfSizeDIM dr_ij;
      dr_ij[0] = coordinates[j][0] - xi;
      dr_ij[1] = coordinates[j][1] - yi;
      dr_ij[2] = coordinates[j][2] - zi;

      double const rij2 =
          dr_ij[0] * dr_ij[0] + dr_ij[1] * dr_ij[1] + dr_ij[2] * dr_ij[2];

      std::size_t const ijparam =
          edip_multi_->elem3param_(species_i, species_j, species_j);

      if (rij2 > edip_multi_->params_[ijparam].cutsq) {
        continue;
      }

      double const rij = std::sqrt(rij2);

      int const contributing_j =
          IsComputeParticleEnergy ? particle_contributing[j] : false;

      // potential energy and force
      // since pair i-j is different from pair j-i, double counting is
      // already considered in constructing the potential
      double fdr;
      double fdZ;
      double const phi = edip_multi_->EdipPair(rij, zeta_i, ijparam, fdr, fdZ);

      double fpair = -fdr / rij;
      dpairZ += fdZ;

      // Contribution to energy
      if (IsComputeEnergy) {
        *energy += phi;
      }

      // Contribution to particle_energy
      if (IsComputeParticleEnergy) {
        if (contributing_j) {
          double const phihalf = 0.5 * phi;
          particle_energy[i] += phihalf;
          particle_energy[j] += phihalf;
        } else {
          particle_energy[i] += phi;
        }
      }  // IsComputeParticleEnergy

      if (IsComputeForces || IsComputeVirial || IsComputeParticleVirial) {
        VectorOfSizeDIM f_ij;
        f_ij[0] = fpair * dr_ij[0];
        f_ij[1] = fpair * dr_ij[1];
        f_ij[2] = fpair * dr_ij[2];

        // Contribution to forces
        if (IsComputeForces) {
          forces[i][0] -= f_ij[0];
          forces[i][1] -= f_ij[1];
          forces[i][2] -= f_ij[2];

          forces[j][0] += f_ij[0];
          forces[j][1] += f_ij[1];
          forces[j][2] += f_ij[2];
        }

        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] = dr_ij[0] * f_ij[0];
          v[1] = dr_ij[1] * f_ij[1];
          v[2] = dr_ij[2] * f_ij[2];
          v[3] = dr_ij[1] * f_ij[2];
          v[4] = dr_ij[0] * f_ij[2];
          v[5] = dr_ij[0] * f_ij[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] *= 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
      }      // IsComputeForces || IsComputeVirial || IsComputeParticleVirial

      // three-body Forces
      for (int kn = jn + 1; kn < number_of_neighbors; ++kn) {
        int const k = neighbors_of_particle[kn];
        std::size_t const species_k =
            static_cast<std::size_t>(particle_species_codes[k]);

        VectorOfSizeDIM dr_ik;
        dr_ik[0] = coordinates[k][0] - xi;
        dr_ik[1] = coordinates[k][1] - yi;
        dr_ik[2] = coordinates[k][2] - zi;

        double const rik2 =
            dr_ik[0] * dr_ik[0] + dr_ik[1] * dr_ik[1] + dr_ik[2] * dr_ik[2];

        std::size_t const ikparam =
            edip_multi_->elem3param_(species_i, species_k, species_k);

        if (rik2 > edip_multi_->params_[ikparam].cutsq) {
          continue;
        }

        double const rik = std::sqrt(rik2);

        std::size_t const ijkparam =
            edip_multi_->elem3param_(species_i, species_j, species_k);

        double const costheta = special::Dot3(dr_ij, dr_ik) / rij / rik;

        double v2;
        double const v1 = edip_multi_->EdipFcut3(rij, ijparam, v2);
        double v4;
        double const v3 = edip_multi_->EdipFcut3(rik, ikparam, v4);
        double v6;
        double v7;
        double const v5 =
            edip_multi_->EdipH(costheta, zeta_i, ijkparam, v6, v7);

        // forces
        dtripleZ += v1 * v3 * v7;

        if (IsComputeEnergy || IsComputeParticleEnergy) {
          // potential energy
          double const phi3 = v1 * v3 * v5;

          // Contribution to energy
          if (IsComputeEnergy) {
            *energy += phi3;
          }  // IsComputeEnergy

          // Contribution to particle_energy
          if (IsComputeParticleEnergy) {
            double const phithird = kThird * phi3;
            particle_energy[i] += phithird;
            if (contributing_j) {
              particle_energy[j] += phithird;
            } else {
              particle_energy[i] += phithird;
            }
            if (particle_contributing[k]) {
              particle_energy[k] += phithird;
            } else {
              particle_energy[i] += phithird;
            }
          }  // IsComputeParticleEnergy
        }    // IsComputeEnergy || IsComputeParticleEnergy

        double const dhl = v1 * v3 * v6;

        VectorOfSizeDIM dri;
        VectorOfSizeDIM drj;
        VectorOfSizeDIM drk;
        special::CosTheta(dr_ij, rij, dr_ik, rik, dri, drj, drk);

        if (IsComputeForces || IsComputeVirial || IsComputeParticleVirial) {
          VectorOfSizeDIM f_ij;
          f_ij[0] = -dhl * drj[0];
          f_ij[1] = -dhl * drj[1];
          f_ij[2] = -dhl * drj[2];

          VectorOfSizeDIM f_ik;
          f_ik[0] = -dhl * drk[0];
          f_ik[1] = -dhl * drk[1];
          f_ik[2] = -dhl * drk[2];

          double dfr = v2 * v3 * v5;
          fpair = -dfr / rij;

          f_ij[0] += fpair * dr_ij[0];
          f_ij[1] += fpair * dr_ij[1];
          f_ij[2] += fpair * dr_ij[2];

          dfr = v1 * v4 * v5;
          fpair = -dfr / rik;

          f_ik[0] += fpair * dr_ik[0];
          f_ik[1] += fpair * dr_ik[1];
          f_ik[2] += fpair * dr_ik[2];

          // Contribution to forces
          if (IsComputeForces) {
            forces[j][0] += f_ij[0];
            forces[j][1] += f_ij[1];
            forces[j][2] += f_ij[2];

            forces[k][0] += f_ik[0];
            forces[k][1] += f_ik[1];
            forces[k][2] += f_ik[2];

            forces[i][0] -= f_ij[0] + f_ik[0];
            forces[i][1] -= f_ij[1] + f_ik[1];
            forces[i][2] -= f_ij[2] + f_ik[2];
          }  // IsComputeForces

          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] = dr_ij[0] * f_ij[0] + dr_ik[0] * f_ik[0];
            v[1] = dr_ij[1] * f_ij[1] + dr_ik[1] * f_ik[1];
            v[2] = dr_ij[2] * f_ij[2] + dr_ik[2] * f_ik[2];
            v[3] = dr_ij[1] * f_ij[2] + dr_ik[1] * f_ik[2];
            v[4] = dr_ij[0] * f_ij[2] + dr_ik[0] * f_ik[2];
            v[5] = dr_ij[0] * f_ij[1] + dr_ik[0] * f_ik[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
        }      // IsComputeForces || IsComputeVirial || IsComputeParticleVirial
      }        // Loop over kn
    }          // Loop over jn

    // forces due to environment coordination f(Z)
    if (IsComputeForces || IsComputeVirial || IsComputeParticleVirial) {
      double const dpairZ_plus_dtripleZ = dpairZ + dtripleZ;
      pre_force_coord = &edip_multi_->pre_force_coord_data_[0];

      for (int idx = 0; idx < number_pre_force_coord_pairs; ++idx) {
        double const dzetair = pre_force_coord[0] * dpairZ_plus_dtripleZ;
        double const dx = pre_force_coord[1];
        double const dy = pre_force_coord[2];
        double const dz = pre_force_coord[3];
        int const j = static_cast<int>(pre_force_coord[4]);

        pre_force_coord += 5;

        double const dx_dzetair = dx * dzetair;
        double const dy_dzetair = dy * dzetair;
        double const dz_dzetair = dz * dzetair;

        // Contribution to forces
        if (IsComputeForces) {
          forces[i][0] -= dx_dzetair;
          forces[i][1] -= dy_dzetair;
          forces[i][2] -= dz_dzetair;

          forces[j][0] += dx_dzetair;
          forces[j][1] += dy_dzetair;
          forces[j][2] += dz_dzetair;
        }  //

        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 * dx_dzetair;
          v[1] = dy * dy_dzetair;
          v[2] = dz * dz_dzetair;
          v[3] = dz * dy_dzetair;
          v[4] = dz * dx_dzetair;
          v[5] = dx * dy_dzetair;

          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
      }      // Loop over idx
    }        // IsComputeForces || IsComputeVirial || IsComputeParticleVirial
  }          // Loop over i

  // everything is good
  return false;
}

#undef KIM_LOGGER_OBJECT_NAME