//
// CDDL HEADER START
//
// The contents of this file are subject to the terms of the Common Development
// and Distribution License Version 1.0 (the "License").
//
// You can obtain a copy of the license at
// http://www.opensource.org/licenses/CDDL-1.0.  See the License for the
// specific language governing permissions and limitations under the License.
//
// When distributing Covered Code, include this CDDL HEADER in each file and
// include the License file in a prominent location with the name LICENSE.CDDL.
// If applicable, add the following below this CDDL HEADER, with the fields
// enclosed by brackets "[]" replaced with your own identifying information:
//
// Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
//
// CDDL HEADER END
//

//
// Copyright (c) 2013, Regents of the University of Minnesota.
// All rights reserved.
//
// Contributors:
//    Ryan S. Elliott
//    Stephen M. Whalen
//


#ifndef EAM_DYNAMO_IMPLEMENTATION_HPP_
#define EAM_DYNAMO_IMPLEMENTATION_HPP_

#include <map>
#include "EAM_Dynamo.hpp"
#include "KIM_API_status.h"

#define MAXLINE 1024
#define DIMENSION 3
#define ONE 1.0
#define TWO 2.0
#define HALF 0.5

#define MAX_PARAMETER_FILES 20
#define NUMBER_SETFL_COMMENT_LINES 3

#define NUMBER_SPLINE_COEFF 9
#define D2F_LINEAR 0
#define D2F_CONSTANT 1
#define DF_QUADRATIC 2
#define DF_LINEAR 3
#define DF_CONSTANT 4
#define F_CUBIC 5
#define F_QUADRATIC 6
#define F_LINEAR 7
#define F_CONSTANT 8

//==============================================================================
//
// Type definitions, enumerations, and helper function prototypes
//
//==============================================================================

// type declaration for get neighbor functions
typedef int (GetNeighborFunction)(void**, int*, int*, int*, int*, int**,
                                  double**);
// type declaration for vector of constant dimension
typedef double VectorOfSizeDIM[DIMENSION];
// type declaration for DYNAMO funcfl data
struct SetOfFuncflData
{
  int numberRhoPoints[MAX_PARAMETER_FILES];
  double deltaRho[MAX_PARAMETER_FILES];
  int numberRPoints[MAX_PARAMETER_FILES];
  double deltaR[MAX_PARAMETER_FILES];
  double cutoff[MAX_PARAMETER_FILES];
  double* embeddingData[MAX_PARAMETER_FILES];
  double* densityData[MAX_PARAMETER_FILES];
  double* ZData[MAX_PARAMETER_FILES];
};
// type declaration for deferred neighbors
typedef struct
{
  int index;
  VectorOfSizeDIM rij;
} neighbor;
// type declaration for iterating over deferred neighbor lists
typedef std::multimap<int, neighbor>::iterator deferredNeighborIterator;

// enumeration for the different types of NBC's
enum NBCTypeEnum {Neigh_Rvec, Neigh_Pure, Mi_Opbc, Cluster};
// enumeration for the different methods of computing rij
enum RijEnum {Coordinates, RVec, MI_OPBC};
// enumeration for EAMFileType
enum EAMFileType {Setfl, Funcfl, FinnisSinclair, Error};

//==============================================================================
//
// Helper class definitions
//
//==============================================================================

// Iterator object for CLUSTER NBC
class ClusterIterator
{
 private:
  int* list_;
  int baseconvert_;
  int const cachedNumberContributingParticles_;
  int request_;
 public:
  ClusterIterator(KIM_API_model* const pkim,
                  GetNeighborFunction* const get_neigh,
                  int const baseconvert,
                  int const cachedNumberContributingParticles,
                  int* const i,
                  int* const numnei,
                  int** const n1atom,
                  double** const pRij)
      : baseconvert_(baseconvert),
        cachedNumberContributingParticles_(cachedNumberContributingParticles),
        request_(0)
  {
    // allocate memory for neighbor list
    list_ = new int[cachedNumberContributingParticles_];
    for (int k = 0; k < cachedNumberContributingParticles_; ++k)
      list_[k] = k - baseconvert_;

    *i = request_;
    // CLUSTER always uses half-list behavior
    *numnei = cachedNumberContributingParticles_ - request_ - 1;
    *n1atom = &(list_[request_ + 1]);
    *pRij = NULL;
  }
  ~ClusterIterator()
  {
    delete [] list_;
  }
  bool done() const
  {
    return !(request_ < cachedNumberContributingParticles_);
  }
  int next(int* const i, int* const numnei, int** const n1atom,
           double** const pRij)
  {
    ++request_;

    *i = request_;
    *numnei = cachedNumberContributingParticles_ - request_ - 1;
    *n1atom = &(list_[request_ + 1]);
    *pRij = NULL;
    return KIM_STATUS_OK;
  }
};

// Iterator object for Locator mode access to neighbor list
class LocatorIterator
{
 private:
  KIM_API_model* const pkim_;
  GetNeighborFunction* const get_neigh_;
  int const baseconvert_;
  int const cachedNumberContributingParticles_;
  int request_;
  int const mode_;
 public:
  LocatorIterator(KIM_API_model* const pkim,
                  GetNeighborFunction* const get_neigh,
                  int const baseconvert,
                  int const cachedNumberContributingParticles,
                  int* const i,
                  int* const numnei,
                  int** const n1atom,
                  double** const pRij)
      : pkim_(pkim),
        get_neigh_(get_neigh),
        baseconvert_(baseconvert),
        cachedNumberContributingParticles_(cachedNumberContributingParticles),
        request_(-baseconvert_),  // set to first value (test-based indexing)
    mode_(1)  // locator mode
  {
    next(i, numnei, n1atom, pRij);
  }
  bool done() const
  {
    return !(request_ + baseconvert_ <= cachedNumberContributingParticles_);
  }
  int next(int* const i, int* const numnei, int** const n1atom,
           double** const pRij)
  {
    int ier;
    // Allow for request_ to be incremented to one more than contributing
    // without causing an error/warning from the openkim-api
    int req = std::min(request_,
                       cachedNumberContributingParticles_-baseconvert_-1);
    ier = (*get_neigh_)(
        reinterpret_cast<void**>(const_cast<KIM_API_model**>(&pkim_)),
        (int*) &mode_,
        &req,
        (int*) i,
        (int*) numnei,
        (int**) n1atom,
        (double**) pRij);
    *i += baseconvert_;  // adjust index of current particle

    ++request_;
    return ier;
  }
};

// Iterator object for Iterator mode access to neighbor list
class IteratorIterator
{
 private:
  KIM_API_model* const pkim_;
  GetNeighborFunction* const get_neigh_;
  int const baseconvert_;
  int request_;
  int const mode_;
  int ier_;
 public:
  IteratorIterator(KIM_API_model* const pkim,
                   GetNeighborFunction* const get_neigh,
                   int const baseconvert,
                   int const cachedNumberContributingParticles,
                   int* const i,
                   int* const numnei,
                   int** const n1atom,
                   double** const pRij)
      : pkim_(pkim),
        get_neigh_(get_neigh),
        baseconvert_(baseconvert),
        request_(0),  // reset iterator
        mode_(0),     // iterator mode
        ier_(KIM_STATUS_FAIL)
  {
    ier_ = (*get_neigh_)(
        reinterpret_cast<void**>(const_cast<KIM_API_model**>(&pkim_)),
        (int*) &mode_, &request_, i, numnei, n1atom, pRij);
    if (ier_ != KIM_STATUS_NEIGH_ITER_INIT_OK)
    {
      pkim->report_error(__LINE__, __FILE__, "iterator init failed", ier_);
      ier_ = KIM_STATUS_FAIL;
    }
    request_ = 1;  // increment iterator
    // return initial set of neighbors
    next(i, numnei, n1atom, pRij);
  }
  bool done() const
  {
    return ((ier_ == KIM_STATUS_FAIL) ||
            (ier_ == KIM_STATUS_NEIGH_ITER_PAST_END));
  }
  int next(int* const i, int* const numnei, int** const n1atom,
           double** const pRij)
  {
    ier_ = (*get_neigh_)(
        reinterpret_cast<void**>(const_cast<KIM_API_model**>(&pkim_)),
        (int*) &mode_,
        &request_,
        (int*) i,
        (int*) numnei,
        (int**) n1atom,
        (double**) pRij);
    *i += baseconvert_;  // adjust index of current particle

    return ier_;
  }
};

// helper routine declarations
void AllocateAndInitialize2DArray(double**& arrayPtr, int const extentZero,
                                  int const extentOne);
void Deallocate2DArray(double**& arrayPtr);
void AllocateAndInitialize3DArray(double***& arrayPtr, int const extentZero,
                                  int const extentOne, int const extentTwo);
void Deallocate3DArray(double***& arrayPtr);

//==============================================================================
//
// Declaration of EAM_DynamoImplementation class
//
//==============================================================================

//******************************************************************************
class EAM_DynamoImplementation
{
 public:
  EAM_DynamoImplementation(KIM_API_model* const pkim,
                           char const* const parameterFileNames,
                           int const parameterFileNameLength,
                           int const numberParameterFiles,
                           int* const ier);
  ~EAM_DynamoImplementation();  // no explicit Destroy() needed in this case

  int Reinit(KIM_API_model* pkim);
  int Compute(KIM_API_model* pkim);

 private:
  // Constant values that never change
  //
  //
  // KIM API Conventions
  int baseconvert_;
  NBCTypeEnum NBCType_;
  bool isHalf_;
  bool isLocatorMode_;
  //
  // KIM API Model Input indices
  int numberParticleTypesIndex_;
  int numberOfParticlesIndex_;
  int numberContributingParticlesIndex_;
  int particleTypesIndex_;
  int coordinatesIndex_;
  int boxSideLengthsIndex_;
  int get_neighIndex_;
  int process_dEdrIndex_;
  int process_d2Edr2Index_;
  //
  // KIM API Model Output indices
  int cutoffIndex_;
  int energyIndex_;
  int forcesIndex_;
  int particleEnergyIndex_;
  //
  // EAM_DynamoImplementation constants
  int numberModelTypes_;


  // Constant values that are read from the input files and never change
  //
  //
  // KIM API Model Fixed Parameters
  char* comments_ptr_[MAX_PARAMETER_FILES];
  char comments_[MAX_PARAMETER_FILES][MAXLINE];
  char particleNames_[MAXLINE];
  int* particleNumber_;
  double* particleMass_;
  double* latticeConstant_;
  char** latticeType_;
  int numberRhoPoints_;
  int numberRPoints_;
  //
  // KIM API Model Free Parameters whoses (pointer) values never change
  double** embeddingData_;
  double*** densityData_;
  double*** rPhiData_;


  // Mutable values that only change when reinit() executes
  //
  //
  // KIM API Model Free Parameters
  double cutoffParameter_;
  double deltaR_;
  double deltaRho_;
  //
  // EAM_DynamoImplementation values
  double cutoffSq_;
  double oneByDr_;
  double oneByDrho_;
  double** embeddingCoeff_;
  double*** densityCoeff_;
  double*** rPhiCoeff_;


  // Mutable values that can change with each call to Reinit() and Compute*()
  //
  //
  // EAM_DynamoImplementation values that change
  int cachedNumberOfParticles_;
  int cachedNumberContributingParticles_;
  double* densityValue_;
  double* embeddingDerivativeValue_;
  double* embeddingSecondDerivativeValue_;


  // Helper methods
  //
  //
  // Related to constructor
  int SetConstantValues(KIM_API_model* const pkim);
  int DetermineNBCTypeAndHalf(KIM_API_model* const pkim);
  void AllocateFixedParameterMemory();
  static int OpenParameterFiles(
      KIM_API_model* const pkim,
      char const* const parameterFileNames,
      int const parameterFileNameLength,
      int const numberParameterFiles,
      FILE* parameterFilePointers[MAX_PARAMETER_FILES]);
  static EAMFileType DetermineParameterFileTypes(
      KIM_API_model* const pkim,
      int const numberModelTypes,
      FILE* const parameterFilePointers[MAX_PARAMETER_FILES],
      int const numberParameterFiles);
  static EAMFileType IsFuncflOrSetfl(FILE* const fptr);
  static EAMFileType IsSetflOrFinnisSinclair(KIM_API_model* const pkim,
                                             FILE* const fptr);
  int ProcessParameterFileHeaders(
      KIM_API_model* const pkim,
      EAMFileType const eamFileType,
      FILE* const parameterFilePointers[MAX_PARAMETER_FILES],
      int const numberParameterFiles, SetOfFuncflData& funcflData);
  int ReadDynamoSetflHeader(KIM_API_model* const pkim, FILE* const fptr);
  int ReadDynamoFuncflHeader(KIM_API_model* const pkim, FILE* const fptr,
                             int const fileIndex, int& numberRhoPoints,
                             double& deltaRho, int& numberRPoints,
                             double& deltaR, double& cutoffParameter);
  int SetParticleNamesForDynamoFuncflModels(KIM_API_model* const pkim);
  void AllocateFreeParameterMemory();
  int ProcessParameterFileData(
      KIM_API_model* const pkim,
      EAMFileType const eamFileType,
      FILE* const parameterFilePointers[MAX_PARAMETER_FILES],
      int const numberParameterFiles, SetOfFuncflData& funcflData);
  int ReadDynamoSetflData(KIM_API_model* const pkim, FILE* const fptr);
  int ReadDynamoFinnisSinclairData(KIM_API_model* const pkim, FILE* const fptr);
  static int ReadDynamoFuncflData(KIM_API_model* const pkim, FILE* const fptr,
                                  int const fileIndex,
                                  SetOfFuncflData& funcflData);
  static int GrabData(KIM_API_model* const pkim, FILE* const fptr, int const n,
                      double* const list);
  void ReinterpolateAndMix(SetOfFuncflData const& funcflData);
  static void CloseParameterFiles(
      FILE* const parameterFilePointers[MAX_PARAMETER_FILES],
      int const numberParameterFiles);
  int ConvertUnits(KIM_API_model* const pkim);
  int RegisterKIMParameters(KIM_API_model* const pkim,
                            EAMFileType const eamFileType) const;
  int RegisterKIMFunctions(KIM_API_model* const pkim) const;
  //
  // Related to Reinit()
  int SetReinitMutableValues(KIM_API_model* const pkim);
  void CubicSplineInterpolateAllData();
  static void CubicSplineInterpolate(double const* const dat,
                                     double const delta, int const n,
                                     double* const coe);
  //
  // Related to Compute()
  int SetComputeMutableValues(KIM_API_model* const pkim,
                              bool& isComputeProcess_dEdr,
                              bool& isComputeProcess_d2Edr2,
                              bool& isComputeEnergy,
                              bool& isComputeForces,
                              bool& isComputeParticleEnergy,
                              int const*& particleTypes,
                              GetNeighborFunction *& get_neigh,
                              double const*& boxSideLengths,
                              VectorOfSizeDIM const*& coordinates,
                              double*& energy,
                              double*& particleEnergy,
                              VectorOfSizeDIM*& forces);
  int CheckParticleTypes(KIM_API_model* const pkim,
                         int const* const particleTypes) const;
  int GetComputeIndex(const bool& isComputeProcess_dEdr,
                      const bool& isComputeProcess_d2Edr2,
                      const bool& isComputeEnergy,
                      const bool& isComputeForces,
                      const bool& isComputeParticleEnergy) const;
  static void ApplyMIOPBC(double const* const boxSideLengths, double* const dx);

  // compute functions
  template< class Iter, bool isHalf, RijEnum rijMethod,
            bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
            bool isComputeEnergy, bool isComputeForces,
            bool isComputeParticleEnergy >
  int Compute(KIM_API_model* const pkim,
              const int* const particleTypes,
              GetNeighborFunction* const get_neigh,
              const double* const boxSideLengths,
              const VectorOfSizeDIM* const coordinates,
              double* const energy,
              VectorOfSizeDIM* const forces,
              double* const particleEnergy);
};

//==============================================================================
//
// Definition of MACROs for improved efficiency
//
//==============================================================================

//******************************************************************************
// MACRO to compute parameters for cubic spline interpolation
// (used for efficiency)
//
// X - function argument
// H - 1/dX where dX is the spline knot spacing
// N - number of knots in the spline
// INDX = int(X/dX) * NUMBER_SPLINE_COEFF
// DELTAX = X/dX - int(X/dX)
#define GET_DELTAX_AND_INDEX(X, H, N, DELTAX, INDX)             \
  DELTAX  = X * H;                                              \
  INDX    = static_cast<int>(DELTAX);                           \
  INDX    = std::min(INDX, N - 2);                              \
  DELTAX -= static_cast<double>(INDX);                          \
  DELTAX  = std::min(DELTAX, 1.0);                              \
  INDX   *= NUMBER_SPLINE_COEFF;

//******************************************************************************
// MACRO to interpolate F(X) (used for efficiency)
//
// DX - DELTAX as computed in GET_DELTAX_AND_INDEX
// I  - INDX   as computed in GET_DELTAX_AND_INDEX
// F  - F(X)
#define INTERPOLATE_F(COEFF, DX, I, F)                    \
  F = COEFF[I + F_CUBIC] * DX + COEFF[I + F_QUADRATIC];   \
  F = F * DX + COEFF[I + F_LINEAR];                       \
  F = F * DX + COEFF[I + F_CONSTANT];

//******************************************************************************
// MACRO to interpolate dF(X)/dX (used for efficiency)
//
// DX - DELTAX as computed in GET_DELTAX_AND_INDEX
// I  - INDX   as computed in GET_DELTAX_AND_INDEX
// DF - dF(X)/dX
#define INTERPOLATE_DF(COEFF, DX, I, DF)                        \
  DF = COEFF[I + DF_QUADRATIC] * DX + COEFF[I + DF_LINEAR];     \
  DF = DF * DX + COEFF[I + DF_CONSTANT];

//******************************************************************************
// MACRO to interpolate d^2F(X)/dX^2 (used for efficiency)
//
// DX - DELTAX as computed in GET_DELTAX_AND_INDEX
// I  - INDX   as computed in GET_DELTAX_AND_INDEX
// D2F- d^2F(X)/dX^2
#define INTERPOLATE_D2F(COEFF, DX, I, D2F)                      \
  D2F = COEFF[I + D2F_LINEAR] * DX + COEFF[I + D2F_CONSTANT];

//==============================================================================
//
// Definition of EAM_DynamoImplementation::Compute functions
//
// NOTE: Here we rely on the compiler optimizations to prune dead code
//       after the template expansions.  This provides high efficiency
//       and easy maintenance.
//
//==============================================================================

//******************************************************************************
template< class Iter, bool isHalf, RijEnum rijMethod,
          bool isComputeProcess_dEdr, bool isComputeProcess_d2Edr2,
          bool isComputeEnergy, bool isComputeForces,
          bool isComputeParticleEnergy >
int EAM_DynamoImplementation::Compute(
    KIM_API_model* const pkim,
    const int* const particleTypes,
    GetNeighborFunction* const get_neigh,
    const double* const boxSideLengths,
    const VectorOfSizeDIM* const coordinates,
    double* const energy,
    VectorOfSizeDIM* const forces,
    double* const particleEnergy)
{
  int ier = KIM_STATUS_OK;

  // initialize electron density for each contributing particle
  for (int i = 0; i < cachedNumberContributingParticles_; ++i)
  {
    densityValue_[i] = 0.0;
    // no need to initialize embeddingDerivativeValue_
  }
  if (isComputeEnergy == true)
  {
    *energy = 0.0;
  }
  if (isComputeForces == true)
  {
    for (int i = 0; i < cachedNumberOfParticles_; ++i)
    {
      for (int j = 0; j < DIMENSION; ++j)
        forces[i][j] = 0.0;
    }
  }

  // compute electron density
  // Setup loop over contributing particles
  int i = 0;
  int numnei = 0;
  int* n1atom = 0;
  double* pRij = 0;
  for (Iter iterator(pkim, get_neigh, baseconvert_,
                     cachedNumberContributingParticles_, &i, &numnei, &n1atom,
                     &pRij);
       iterator.done() == false;
       ier = iterator.next(&i, &numnei, &n1atom, &pRij))
  {
    if (ier < KIM_STATUS_OK) // check that iterator.next was successful
    {
      pkim->report_error(__LINE__, __FILE__, "get_neigh", ier);
      return ier;
    }

    // Setup loop over neighbors of current particle
    for (int jj = 0; jj < numnei; ++jj)
    {  // adjust index of particle neighbor
      int const j = n1atom[jj] + baseconvert_;
      double* r_ij;
      double r_ijValue[DIMENSION];
      // Compute r_ij appropriately
      switch (rijMethod)
      {
        case Coordinates:
        {
          r_ij = r_ijValue;
          for (int k = 0; k < DIMENSION; ++k)
            r_ij[k] = coordinates[j][k] - coordinates[i][k];
          break;
        }
        case MI_OPBC:
        {
          r_ij = r_ijValue;
          for (int k = 0; k < DIMENSION; ++k)
            r_ij[k] = coordinates[j][k] - coordinates[i][k];
          // apply minimum image convention
          ApplyMIOPBC(boxSideLengths, r_ij);
          break;
        }
        case RVec:
        {
          r_ij = &pRij[jj * DIMENSION];
          break;
        }
      }

      // compute distance squared
      double rij2 = 0.0;
      for (int k = 0; k < DIMENSION; ++k)
        rij2 += r_ij[k] * r_ij[k];

      if (rij2 <= cutoffSq_)
      { // compute contribution to electron density
        double rijOffset;
        int rijIndex;
        double const rij = sqrt(rij2);

        // compute rijOffset and rijIndex
        GET_DELTAX_AND_INDEX(rij, oneByDr_, numberRPoints_, rijOffset,
                             rijIndex);

        // interpolate value of rho_beta(r_ij)
        double densityBetaValue;
        double const* const densityBetaCoeff
            = densityCoeff_[particleTypes[i]][particleTypes[j]];
        INTERPOLATE_F(densityBetaCoeff, rijOffset, rijIndex, densityBetaValue);
        densityValue_[i] += densityBetaValue;

        // compute ED contribution from neighbor if half list
        if (isHalf == true)
        {
          if (j < cachedNumberContributingParticles_)
          { // if using half list and j contributes, add its contribution.
            // interpolate value of rho_alpha(r_ij)
            double densityAlphaValue;
            double const* const densityAlphaCoeff
                = densityCoeff_[particleTypes[j]][particleTypes[i]];
            INTERPOLATE_F(densityAlphaCoeff, rijOffset, rijIndex,
                          densityAlphaValue);
            densityValue_[j] += densityAlphaValue;
          }
        }
      }
    }  // end of loop over neighbors
  }  // end of loop over contributing particles

  // calculate embedding function and its derivative
  for (int i = 0; i < cachedNumberContributingParticles_; ++i)
  {
    double densityOffset;
    int densityIndex;
    // compute densityOffset and densityIndex
    GET_DELTAX_AND_INDEX(densityValue_[i], oneByDrho_, numberRhoPoints_,
                         densityOffset, densityIndex);
    // interpolate F_i(rho_i)
    double embeddingValue;
    double const* const embeddingAlphaCoeff
        = embeddingCoeff_[particleTypes[i]];
    INTERPOLATE_F(embeddingAlphaCoeff, densityOffset, densityIndex,
                  embeddingValue);
    // Contribute embedding term to Energy
    if (isComputeEnergy == true)
    {
      *energy += embeddingValue;
    }
    // Contribute embedding term to ParticleEnergy
    if (isComputeParticleEnergy == true)
    {
      particleEnergy[i] = embeddingValue;
    }
    // Compute embedding derivative
    if ((isComputeForces == true) || (isComputeProcess_dEdr == true) ||
        (isComputeProcess_d2Edr2 == true))
    {
      // interpolate dF_i(rho_i)/d(rho_i)
      INTERPOLATE_DF(embeddingAlphaCoeff, densityOffset, densityIndex,
                     embeddingDerivativeValue_[i]);
    }
    // Compute embedding second derivative
    if (isComputeProcess_d2Edr2 == true)
    {
      // interpolate d^2F_i(rho_i)/d(rho_i)^2
      INTERPOLATE_D2F(embeddingAlphaCoeff, densityOffset, densityIndex,
                      embeddingSecondDerivativeValue_[i]);
    }
  }

  // calculate contribution from electron density to the force part
  // and from pair function
  //
  // Setup loop over contributing particles
  for (Iter iterator(pkim, get_neigh, baseconvert_,
                     cachedNumberContributingParticles_, &i, &numnei, &n1atom,
                     &pRij);
       iterator.done() == false;
       iterator.next(&i, &numnei, &n1atom, &pRij))
  {
    // Setup loop over neighbors of current particle
    for (int jj = 0; jj < numnei; ++jj)
    { // adjust index of particle neighbor
      int const j = n1atom[jj] + baseconvert_;
      double* r_ij;
      double r_ijValue[DIMENSION];
      // Compute r_ij appropriately
      switch (rijMethod)
      {
        case Coordinates:
        {
          r_ij = r_ijValue;
          for (int k = 0; k < DIMENSION; ++k)
            r_ij[k] = coordinates[j][k] - coordinates[i][k];
          break;
        }
        case MI_OPBC:
        {
          r_ij = r_ijValue;
          for (int k = 0; k < DIMENSION; ++k)
            r_ij[k] = coordinates[j][k] - coordinates[i][k];
          // apply minimum image convention
          ApplyMIOPBC(boxSideLengths, r_ij);
          break;
        }
        case RVec:
        {
          r_ij = &pRij[jj * DIMENSION];
          break;
        }
      }

      // compute distance squared
      double rij2 = 0.0;
      for (int k = 0; k < DIMENSION; ++k)
        rij2 += r_ij[k] * r_ij[k];

      if (rij2 <= cutoffSq_)
      { // compute contribution to energy and force
        double rijOffset;
        int rijIndex;
        double const rij = sqrt(rij2);

        // compute rijOffset and rijIndex
        GET_DELTAX_AND_INDEX(rij, oneByDr_, numberRPoints_, rijOffset,
                             rijIndex);

        // interpolate r_ij*phi(r_ij)
        double rijPhiValue;
        double const* const rijPhiAlphaBetaCoeff
            = rPhiCoeff_[particleTypes[i]][particleTypes[j]];
        INTERPOLATE_F(rijPhiAlphaBetaCoeff, rijOffset, rijIndex, rijPhiValue);

        // find phi(r_ij)
        double const oneByRij = ONE / rij;
        double const pairPotentialValue = rijPhiValue * oneByRij;

        // Contribute pair term to Energy as half or full
        if (isComputeEnergy == true)
        {
          if (isHalf == true)
          {
            if (j < cachedNumberContributingParticles_)
            {
              *energy += pairPotentialValue;
            }
            else
            {
              *energy += HALF * pairPotentialValue;
            }
          }
          else
          {
            *energy += HALF * pairPotentialValue;
          }
        }

        // Contribute pair term to Particle Energy
        if (isComputeParticleEnergy == true)
        {
          if (isHalf == true)
          {
            if (i < cachedNumberContributingParticles_)
            {
              particleEnergy[i] += HALF * pairPotentialValue;
            }
            if (j < cachedNumberContributingParticles_)
            {
              particleEnergy[j] += HALF * pairPotentialValue;
            }
          }
          else
          {
            particleEnergy[i] += HALF * pairPotentialValue;
          }
        }

        // Compute dEdrByR terms as half or full
        double dEdrByRij = 0.0;
        if ((isComputeForces == true) || (isComputeProcess_dEdr == true))
        {
          if (isHalf == true)
          { // interpolate derivative of r_ij*phi(r_ij) function
            double rijPhiDerivativeValue;
            INTERPOLATE_DF(rijPhiAlphaBetaCoeff, rijOffset, rijIndex,
                           rijPhiDerivativeValue);

            // interpolate derivative of rho_beta(r_ij)
            double densityBetaDerivativeValue;
            double const* const densityBetaCoeff
                = densityCoeff_[particleTypes[i]][particleTypes[j]];
            INTERPOLATE_DF(densityBetaCoeff, rijOffset, rijIndex,
                           densityBetaDerivativeValue);

            // compute dEdr contribution
            if (j < cachedNumberContributingParticles_)
            { // interpolate derivative of rho_alpha(r_ij)
              double densityAlphaDerivativeValue;
              double const* const densityAlphaCoeff
                  = densityCoeff_[particleTypes[j]][particleTypes[i]];
              INTERPOLATE_DF(densityAlphaCoeff, rijOffset, rijIndex,
                             densityAlphaDerivativeValue);

              // embedding contribution to dEdr
              double const embeddingContribution
                  = ((embeddingDerivativeValue_[i] * densityBetaDerivativeValue)
                     +
                     (embeddingDerivativeValue_[j] *
                      densityAlphaDerivativeValue));

              // pair potential contribution to dEdr
              double const pairPotentialContribution
                  = (rijPhiDerivativeValue - pairPotentialValue) * oneByRij;

              // divide by r so we can multiply by r_ij below
              dEdrByRij = (embeddingContribution + pairPotentialContribution)
                  * oneByRij;
            }
            else
            { // embedding contribution to dEdr
              double const embeddingContribution
                  = (embeddingDerivativeValue_[i] * densityBetaDerivativeValue);

              // pair potential contribution
              double const pairPotentialContribution
                  = HALF * (rijPhiDerivativeValue - pairPotentialValue)
                  * oneByRij;

              // divide by r so we can multiply by r_ij below
              dEdrByRij = (embeddingContribution + pairPotentialContribution)
                  * oneByRij;
            }
          }
          else
          { // interpolate derivative of r_ij*phi(r_ij) function
            double rijPhiDerivativeValue;
            INTERPOLATE_DF(rijPhiAlphaBetaCoeff, rijOffset, rijIndex,
                           rijPhiDerivativeValue);

            // interpolate derivative of rho_beta(r_ij)
            double densityBetaDerivativeValue;
            double const* const densityBetaCoeff
                = densityCoeff_[particleTypes[i]][particleTypes[j]];
            INTERPOLATE_DF(densityBetaCoeff, rijOffset, rijIndex,
                           densityBetaDerivativeValue);

            // compute dEdr contribution
            // embedding contribution to dEdr
            double const embeddingContribution
                = (embeddingDerivativeValue_[i] * densityBetaDerivativeValue);

            // pair potential contribution
            double const pairPotentialContribution
                = HALF * (rijPhiDerivativeValue - pairPotentialValue)
                * oneByRij;

            // divide by r so we can multiply by r_ij below
            dEdrByRij = (embeddingContribution + pairPotentialContribution)
                * oneByRij;
          }
        }

        // Contribute dEdrByR to forces
        if (isComputeForces == true)
        {
          for (int k = 0; k < DIMENSION; ++k)
          {
            forces[i][k] += dEdrByRij * r_ij[k];
            forces[j][k] -= dEdrByRij * r_ij[k];
          }
        }

        // Call process_dEdr
        if (isComputeProcess_dEdr == true)
        {
          double const dEdr = dEdrByRij * rij;
          ier = pkim->process_dEdr(const_cast<KIM_API_model**>(&pkim),
                                   const_cast<double*>(&dEdr),
                                   const_cast<double*>(&rij),
                                   &r_ij,
                                   &i,
                                   const_cast<int*>(&j));
          if (ier < KIM_STATUS_OK)
          {
            pkim->report_error(__LINE__, __FILE__, "process_dEdr", ier);
            return ier;
          }
        }
      }  // if particles i and j interact
    }  // end of first neighbor loop
  }  // end of loop over contributing particles

  // Separate loop nest for process_d2Edr2
  if (isComputeProcess_d2Edr2 == true)
  { // Separate implementations for half and full cases
    //
    // For this potential, the second derivative
    // d^2E/dr_{ij}dr_{kl} is nonzero only if i,j,k,l are
    // not distinct. As such, we need only address all
    // derivatives d^2E/dr_{ij}dr_{ik}.
    //
    // If we have full lists, we can simply iterate over
    // particle i's neighbor list in a doubly-nested
    // triangular loop.
    //
    // If we have half lists, then the presence of j in
    // i's neighbor list implies that i is also a neighbor
    // of j, so that we should be processing all derivs
    // d^2E/dr_{ij}dr_{ik} and d^2E/dr_{ji}dr_{jl}. The
    // hurdle to overcome here is the need to access the
    // neighbor list of particle j (and also particle k).
    //
    // We avoid this problem by taking a different
    // approach. For each neighbor (j or k) of i, we save
    // i in a "deferred neighbor list" for each of j and
    // k. For each particle, we process neighbors both in
    // its neighbor list and its deferred neighbor list.
    // (In effect, we are reconstructing full lists for
    // each particle.)

    // Containers for deferred neighbor lists for process_d2Edr2
    std::multimap<int, neighbor> deferredNeighborMap;

    // Iterators over deferred neighbor lists
    std::pair<deferredNeighborIterator, deferredNeighborIterator>
        deferredNeighborRange;
    deferredNeighborIterator deferredNeighborPtrJ;

    // Setup loop over contributing particles
    for (Iter iterator(pkim, get_neigh, baseconvert_,
                       cachedNumberContributingParticles_, &i, &numnei, &n1atom,
                       &pRij);
         iterator.done() == false;
         iterator.next(&i, &numnei, &n1atom, &pRij))
    { // Number of neighbors to visit. May be greater than
      // numnei if previous particles inserted deferred
      // neighbors to this particle, or if this particle
      // is its own neighbor.
      int extendedNumnei;
      extendedNumnei = numnei;

      // Begin by deferring i to each of its neighbors.
      // It is tempting to try to include this in the
      // subsequent loop that also does the calculations,
      // but this cannot be done. The complete deferred
      // list must be available the first time we go
      // through the inner, second-neighbor, loop.
      if (isHalf == true)
      {
        for (int jj = 0; jj < numnei; ++jj)
        {
          neighbor reflectedNeighJ;

          // Index of particle neighbor
          int const atomj = n1atom[jj];
          // Adjust index of particle neighbor
          int const j = atomj + baseconvert_;

          // Save "unadjusted" particle index
          reflectedNeighJ.index = i - baseconvert_;

          // Save displacement if necessary
          if (RVec == rijMethod)
          {
            for (int k = 0; k < DIMENSION; ++k)
              reflectedNeighJ.rij[k] = -pRij[jj * DIMENSION + k];
          }

          deferredNeighborMap.insert(std::make_pair(j, reflectedNeighJ));
        }

        // No further particles will contribute to i's neighbor list.
        // Safe to finalize the neighbor count.
        extendedNumnei += deferredNeighborMap.count(i);
        deferredNeighborRange = deferredNeighborMap.equal_range(i);
        deferredNeighborPtrJ  = deferredNeighborRange.first;
      }

      // Setup loop over neighbors of current particle
      for (int jj = 0; jj < extendedNumnei; ++jj)
      {
        int atomj;

        if (isHalf == true && jj >= numnei)
        { // Look in deferred neighbor list
          atomj = deferredNeighborPtrJ->second.index;
        }
        else
        { // Not a deferred neighbor
          atomj = n1atom[jj];
        }

        // adjust index of particle neighbor
        int const j = atomj + baseconvert_;

        // Declare enough space to hold two displacements,
        // one for each neighbor. Ensuring that the two
        // displacements are contiguous in memory allows
        // the use of r_ijValue in the call to process_d2Edr2.
        double r_ijValue[2*DIMENSION];

        // Pointer to the first displacement.
        // Perhaps confusingly, this might be set to point
        // into the RVEC displacements list (pRij), rather
        // than always pointing into r_ijValue.
        double* r_ij;

        // Compute r_ij appropriately
        switch (rijMethod)
        {
          case Coordinates:
          {
            r_ij = r_ijValue;
            for (int k = 0; k < DIMENSION; ++k)
              r_ij[k] = coordinates[j][k] - coordinates[i][k];
            break;
          }
          case MI_OPBC:
          {
            r_ij = r_ijValue;
            for (int k = 0; k < DIMENSION; ++k)
              r_ij[k] = coordinates[j][k] - coordinates[i][k];
            // apply minimum image convention
            ApplyMIOPBC(boxSideLengths, r_ij);
            break;
          }
          case RVec:
          {
            if (isHalf == true && jj >= numnei)
            { // Look in deferred neighbor list
              r_ij = r_ijValue;
              for (int k = 0; k < DIMENSION; ++k)
                r_ij[k] = deferredNeighborPtrJ->second.rij[k];
            }
            else
            {
              r_ij = &pRij[jj * DIMENSION];
            }
            break;
          }
        }

        // compute distance squared
        double rij2 = 0.0;
        for (int k = 0; k < DIMENSION; ++k)
          rij2 += r_ij[k] * r_ij[k];

        if (rij2 <= cutoffSq_)
        {
          double rijOffset;
          int rijIndex;
          double rij[2] = { sqrt(rij2), 0.0 };  // rij[1] used for ik pair

          // compute rijOffset and rijIndex
          GET_DELTAX_AND_INDEX(rij[0], oneByDr_, numberRPoints_, rijOffset,
                               rijIndex);

          // interpolate r_ij*phi(r_ij)
          double rijPhiValue;
          double const* const rijPhiAlphaBetaCoeff
              = rPhiCoeff_[particleTypes[i]][particleTypes[j]];
          INTERPOLATE_F(rijPhiAlphaBetaCoeff, rijOffset, rijIndex, rijPhiValue);

          // find phi(r_ij)
          double const oneByRij = ONE / rij[0];
          double const pairPotentialValue = rijPhiValue * oneByRij;

          // interpolate derivative of r_ij*phi(r_ij) function
          double rijPhiDerivativeValue;
          INTERPOLATE_DF(rijPhiAlphaBetaCoeff, rijOffset, rijIndex,
                         rijPhiDerivativeValue);

          // find derivative of phi(r_ij)
          double const pairPotentialDerivativeValue
              = (rijPhiDerivativeValue - pairPotentialValue) * oneByRij;

          // interpolate derivative of rho_beta(r_ij)
          double densityBetaDerivativeValue;
          double const* const densityBetaCoeff
              = densityCoeff_[particleTypes[i]][particleTypes[j]];
          INTERPOLATE_DF(densityBetaCoeff, rijOffset, rijIndex,
                         densityBetaDerivativeValue);

          // Setup second loop over neighbors of current particle
          deferredNeighborIterator
              deferredNeighborPtrK = deferredNeighborPtrJ;
          for (int kk = jj; kk < extendedNumnei; ++kk)
          {
            int atomk;
            if (isHalf == true && kk >= numnei)
            { // Look in deferred neighbor list
              atomk = deferredNeighborPtrK->second.index;
            }
            else
            {
              atomk = n1atom[kk];
            }

            // adjust index of particle neighbor
            int const k = atomk + baseconvert_;
            double* r_ik;

            // Pointer to the second displacement vector
            r_ik = &(r_ijValue[DIMENSION]);
            switch (rijMethod)
            {
              case Coordinates:
              {
                for (int d = 0; d < DIMENSION; ++d)
                  r_ik[d] = coordinates[k][d] - coordinates[i][d];
                break;
              }
              case MI_OPBC:
              {
                for (int d = 0; d < DIMENSION; ++d)
                  r_ik[d] = coordinates[k][d] - coordinates[i][d];
                // apply minimum image convention
                ApplyMIOPBC(boxSideLengths, r_ik);
                break;
              }
              case RVec:
              { // In this case, we could have earlier set r_ij
                // to point into pRij, and not into r_ijValue.
                // In other words, r_ijValue might not yet
                // contain the components of the first neighbor's
                // displacement. We need them to be there.
                for (int d = 0; d < DIMENSION; ++d)
                  r_ijValue[d] = r_ij[d];

                r_ij = r_ijValue;

                if (isHalf == true && kk >= numnei)
                { // Look in deferred neighbor list
                  for (int d = 0; d < DIMENSION; ++d)
                    r_ik[d] = deferredNeighborPtrK->second.rij[d];
                }
                else
                {
                  for (int d = 0; d < DIMENSION; ++d)
                    r_ik[d] = pRij[kk * DIMENSION + d];
                }

                break;
              }
            }

            // compute distance squared
            double rik2 = 0.0;
            for (int d = 0; d < DIMENSION; ++d)
              rik2 += r_ik[d] * r_ik[d];

            if (rik2 <= cutoffSq_)
            {
              double rikOffset;
              int rikIndex;
              rij[1] = sqrt(rik2);

              // compute rikOffset and rikIndex
              GET_DELTAX_AND_INDEX(rij[1], oneByDr_, numberRPoints_, rikOffset,
                                   rikIndex);

              // interpolate derivative of rho_gamma(r_ik)
              double densityGammaDerivativeValue;
              double const* const densityGammaCoeff
                  = densityCoeff_[particleTypes[i]][particleTypes[k]];
              INTERPOLATE_DF(densityGammaCoeff, rikOffset, rikIndex,
                             densityGammaDerivativeValue);

              // mixed-index embedding contribution to d2Edr2
              double const mixedEmbeddingContribution
                  = embeddingSecondDerivativeValue_[i]
                  * densityBetaDerivativeValue
                  * densityGammaDerivativeValue;

              double d2Edr2 = mixedEmbeddingContribution;

              if (kk == jj)
              { // interpolate second derivative of r_ij*phi(r_ij) function
                double rijPhiSecondDerivativeValue;
                INTERPOLATE_D2F(rijPhiAlphaBetaCoeff, rijOffset, rijIndex,
                                rijPhiSecondDerivativeValue);

                // interpolate second derivative of rho_beta(r_ij)
                double densityBetaSecondDerivativeValue;
                double const* const densityBetaCoeff
                    = densityCoeff_[particleTypes[i]][particleTypes[j]];
                INTERPOLATE_D2F(densityBetaCoeff, rijOffset, rijIndex,
                                densityBetaSecondDerivativeValue);

                // pair potential contribution
                double const pairPotentialSecondDerivativeValue
                    = (rijPhiSecondDerivativeValue
                       - TWO * pairPotentialDerivativeValue) * oneByRij;
                double const pairPotentialContribution
                    = HALF * pairPotentialSecondDerivativeValue;

                // second derivative of embedding contribution
                double const embeddingContribution
                    = embeddingDerivativeValue_[i]
                    * densityBetaSecondDerivativeValue;

                d2Edr2 += pairPotentialContribution + embeddingContribution;
              }
              else
              { // Process transpose pair
                const int iis[2] = { i, i };
                const int jjs[2] = { k, j };
                int const* const piis = &iis[0];
                int const* const pjjs = &jjs[0];
                const double rikrij[2] = { rij[1], rij[0] };
                double const* const prikrij = &rikrij[0];
                const double r_ikr_ij[6] = { r_ij[3], r_ij[4], r_ij[5],
                                             r_ij[0], r_ij[1], r_ij[2] };
                double const* const pr_ikr_ij = &r_ikr_ij[0];
                ier = pkim->process_d2Edr2(
                    const_cast<KIM_API_model**>(&pkim),
                    &d2Edr2,
                    const_cast<double**>(&prikrij),
                    const_cast<double**>(&pr_ikr_ij),
                    const_cast<int**>(&piis),
                    const_cast<int**>(&pjjs));
                if (ier < KIM_STATUS_OK)
                {
                  pkim->report_error(__LINE__, __FILE__, "process_d2Edr2", ier);
                  return ier;
                }
              }

              int const iis[2] = { i, i };
              int const jjs[2] = { j, k };
              int const* const piis = &iis[0];
              int const* const pjjs = &jjs[0];
              double const* const prij = &rij[0];
              ier = pkim->process_d2Edr2(
                  const_cast<KIM_API_model**>(&pkim),
                  &d2Edr2,
                  const_cast<double**>(&prij),
                  &r_ij,
                  const_cast<int**>(&piis),
                  const_cast<int**>(&pjjs));
              if (ier < KIM_STATUS_OK)
              {
                pkim->report_error(__LINE__, __FILE__, "process_d2Edr2", ier);
                return ier;
              }
            }  // if particles i and k interact

            if (isHalf == true && kk >= numnei)
            {
              ++deferredNeighborPtrK;
            }
          }  // end of second neighbor loop
        }  // if particles i and j interact

        if (isHalf == true && jj >= numnei)
        {
          ++deferredNeighborPtrJ;
        }
      }  // end of first neighbor loop

      deferredNeighborMap.erase(i);
    }  // end of loop over contributing particles
  }  // if (isComputeProcess_d2Edr2 == true)

  // everything is good
  ier = KIM_STATUS_OK;
  return ier;
}

#endif  // EAM_DYNAMO_IMPLEMENTATION_HPP_