!
! 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:
!
!
! CDDL HEADER END
!

!
! Copyright (c) 2012, Mark R. Gilbert, CCFE Fusion Association.  All rights reserved.
!
! Contributors:
!    Mark R. Gilbert
!

!****************************************************************************
!**
!**  MODULE model_driver_PF_cubic_splines
!**
!**  EAM-like potential with cubic splines representing knot functions
!**  magnetic ability also available via B parameter
!**
!**  Language: Fortran 90
!**
!**
!**
!****************************************************************************


#include "KIM_API_status.h"
#define THIS_FILE_NAME __FILE__
#define TRUEFALSE(TRUTH) merge(1,0,(TRUTH))

module model_driver_pf_cubic_splines

use KIM_API
implicit none

save
private
public Compute_Energy_Forces, &
       model_cutoff,          &
       Destroy, &
       nknotv, nknotp,interpolate_num, &
       vknotcoeff,vknotpoint,linear_interpolate, r1,r2,A,B,Z,a_rho,a_inter, &
       pknotcoeff,pknotpoint

! Below are the definitions and values of all Model parameters
integer, parameter          :: DIM=3          ! dimensionality of space
integer, parameter          :: speccode = 1   ! internal species code
double precision:: model_cutoff  ! cutoff radius
                                                                    ! in angstroms
double precision :: model_cutsq

    double precision , parameter :: pc=1.0d0,pi=3.141592653589793d0
! electric constant and electron charge ( units of
    double precision , parameter :: E_0=8.854187817d-12,Ec=1.60217653d-19
!bohr constant in angstroms
   double precision , parameter :: a_B=0.5291772108d0

! constants of biersack-ziegler coulomb potential
    double precision , parameter :: p1=0.1818d0,p2=0.5099d0,p3=0.2802d0,p4=0.02817d0
    double precision , parameter :: Bt1=-3.2d0,Bt2=-0.9423d0,Bt3=-0.4029d0,Bt4=-0.2016d0
!
! new variables for general potential - to be read in - 12/10/06
    double precision       :: A, B, Z,r1,r2,a_inter(7),a_rho
    integer               :: nknotv,nknotp,interpolate_num
    LOGICAL              :: linear_interpolate
    double precision , allocatable :: vknotcoeff(:), vknotpoint(:), &
                                           pknotcoeff(:), pknotpoint(:)

contains

!-------------------------------------------------------------------------------
!
!  Calculate pair potential phi(r)
!
!-------------------------------------------------------------------------------
subroutine calc_phi(r,phi)
implicit none

!-- Transferred variables
double precision, intent(in)  :: r
double precision, intent(out) :: phi
double precision :: a_s,phi_r,B_Zpre
integer :: i

!-- Local variables
! <FILL place any local variable definitions here>

if (r .gt. model_cutoff) then
   ! Argument exceeds cutoff radius
   phi = 0.d0
else
    IF (r>=r2) THEN
      phi=0.d0
      DO i=1,nknotv
        IF (r<vknotpoint(i)) THEN
         phi=phi+vknotcoeff(i)*(vknotpoint(i)-r)**3
        END IF
      END DO
    ELSE IF (r<=r1) THEN
     a_s=0.88534d0*a_B/(sqrt(2.d0)*(Z**(1.d0/3.d0)))
     B_Zpre=Z**2*Ec*1.0d10/(4.d0*pi*E_0)
     phi_r=p1*exp(Bt1*r/a_s)+p2*exp(Bt2*r/a_s)+p3*exp(Bt3*r/a_s)+ &
               p4*exp(Bt4*r/a_s)
     phi=B_Zpre*phi_r/r
    ELSE
     IF (linear_interpolate) THEN
      phi=a_inter(1)+a_inter(2)*r+a_inter(3)*(r**2)+ &
       a_inter(4)*(r**3)+a_inter(5)*(r**4)+a_inter(6)*(r**5)+a_inter(7)*(r**6)
     ELSE
      phi=exp(a_inter(1)+a_inter(2)*r+a_inter(3)*(r**2)+ &
       a_inter(4)*(r**3)+a_inter(5)*(r**4)+a_inter(6)*(r**5)+a_inter(7)*(r**6))
     END IF
    END IF
endif

end subroutine calc_phi

!-------------------------------------------------------------------------------
!
!  Calculate pair potential phi(r) and its derivative dphi(r)
!
!-------------------------------------------------------------------------------
subroutine calc_phi_dphi(r,phi,dphi)
implicit none

!-- Transferred variables
double precision, intent(in)  :: r
double precision, intent(out) :: phi,dphi
double precision :: a_s,phi_r,B_Zpre,dphi_r
integer :: i


!-- Local variables
! <FILL place any local variable definitions here>

if (r .gt. model_cutoff) then
   ! Argument exceeds cutoff radius
   phi    = 0.d0
   dphi   = 0.d0
else
    IF (r>=r2) THEN
      dphi=0.d0
      phi=0.d0
      DO i=1,nknotv
        IF (r<vknotpoint(i)) THEN
         dphi=dphi-3*vknotcoeff(i)*(vknotpoint(i)-r)**2
         phi=phi+vknotcoeff(i)*(vknotpoint(i)-r)**3
        END IF
      END DO
    ELSE IF (r<=r1) THEN
     a_s=0.88534d0*a_B/(sqrt(2.d0)*(Z**(1.d0/3.d0)))
     B_Zpre=Z**2*Ec*1.0d10/(4.d0*pi*E_0)
     phi_r=p1*exp(Bt1*r/a_s)+p2*exp(Bt2*r/a_s)+p3*exp(Bt3*r/a_s)+ &
               p4*exp(Bt4*r/a_s)
     dphi_r=(p1*Bt1*exp(Bt1*r/a_s)+p2*Bt2*exp(Bt2*r/a_s)+p3*Bt3*exp(Bt3*r/a_s)+ &
               p4*Bt4*exp(Bt4*r/a_s))/a_s
     dphi=-B_Zpre*phi_r/(r**2)+B_Zpre*dphi_r/r
     phi=B_Zpre*phi_r/r
    ELSE
     IF (linear_interpolate) THEN
      dphi=a_inter(2)+2.d0*a_inter(3)*r+3.d0*a_inter(4)*(r**2)+ &
                      4.d0*a_inter(5)*(r**3)+5.d0*a_inter(6)*(r**4)+6.d0*a_inter(7)*(r**5)
      phi=a_inter(1)+a_inter(2)*r+a_inter(3)*(r**2)+ &
          a_inter(4)*(r**3)+a_inter(5)*(r**4)+a_inter(6)*(r**5)+a_inter(7)*(r**6)
     ELSE
       phi_r=exp(a_inter(1)+a_inter(2)*r+a_inter(3)*(r**2)+ &
        a_inter(4)*(r**3)+a_inter(5)*(r**4)+a_inter(6)*(r**5)+a_inter(7)*(r**6))
       dphi_r=(a_inter(2)+2.d0*a_inter(3)*r+3.d0*a_inter(4)*(r**2)+ &
                      4.d0*a_inter(5)*(r**3)+5.d0*a_inter(6)*(r**4)+6.d0*a_inter(7)*(r**5))
       dphi=phi_r*dphi_r
       phi=exp(a_inter(1)+a_inter(2)*r+a_inter(3)*(r**2)+ &
          a_inter(4)*(r**3)+a_inter(5)*(r**4)+a_inter(6)*(r**5)+a_inter(7)*(r**6))
     END IF

    END IF
endif

end subroutine calc_phi_dphi

!-------------------------------------------------------------------------------
!
!  Calculate electron density g(r)
!
!-------------------------------------------------------------------------------
subroutine calc_g(r,g)
implicit none

!-- Transferred variables
double precision, intent(in)  :: r
double precision, intent(out) :: g
integer :: i

!-- Local variables
! <FILL place any local variable definitions here>

if (r .gt. model_cutoff) then
   ! Argument exceeds cutoff radius
   g = 0.d0
else
   g = 0.d0
    DO i=1,nknotp
      IF (r<pknotpoint(i)) THEN
       g=g+pknotcoeff(i)*(pknotpoint(i)-r)**3
      END IF
    END DO
endif

end subroutine calc_g

!-------------------------------------------------------------------------------
!
!  Calculate electron density derivative dg(r)
!
!-------------------------------------------------------------------------------
subroutine calc_dg(r,dg)
implicit none

!-- Transferred variables
double precision, intent(in)  :: r
double precision, intent(out) :: dg
integer :: i

!-- Local variables
! <FILL place any local variable definitions here>

if (r .gt. model_cutoff) then
   ! Argument exceeds cutoff radius
   dg = 0.d0
else
   dg = 0.d0
    DO i=1,nknotp
      IF (r<pknotpoint(i)) THEN
       dg=dg-3.d0*pknotcoeff(i)*(pknotpoint(i)-r)**2
      END IF
    END DO
endif

end subroutine calc_dg

!-------------------------------------------------------------------------------
!
!  Calculate embedding function U(rho)
!
!-------------------------------------------------------------------------------
subroutine calc_U(rho,U)
implicit none

!-- Transferred variables
double precision, intent(in)  :: rho
double precision, intent(out) :: U

!-- Local variables
! <FILL place any local variable definitions here>

    IF (rho<pc) THEN
     U=-A*sqrt(rho)-(B/log(2.0d0))*(1.d0-sqrt(rho/pc))*log(2.0d0-rho/pc)+ &
                     a_rho*(rho**2)
    ELSE
      U=-A*sqrt(rho)+a_rho*(rho**2)
    END IF

end subroutine calc_U

!-------------------------------------------------------------------------------
!
!  Calculate embedding function U(rho) and first derivative dU(rho)
!
!-------------------------------------------------------------------------------
subroutine calc_U_dU(rho,U,dU)
implicit none

!-- Transferred variables
double precision, intent(in)  :: rho
double precision, intent(out) :: U,dU

!-- Local variables
! <FILL place any local variable definitions here>
 IF(rho.le.1e-10) THEN
    dU=0d0
    U=0.d0
 ELSE
    IF (rho<pc) THEN
     U=-A*sqrt(rho)-(B/log(2.0d0))*(1.d0-sqrt(rho/pc))*log(2.0d0-rho/pc)+ &
                     a_rho*(rho**2)
     dU=-A/(2.d0*sqrt(rho))-(B/log(2.d0))*( &
               (-1.d0/(2.d0*sqrt(rho*pc)))*log(2.d0-rho/pc)+ &
               (1.d0-sqrt(rho/pc))*(-1.d0/(2.d0*pc-rho)))+2.d0*a_rho*rho
    ELSE
      U=-A*sqrt(rho)+a_rho*(rho**2)
      dU=-A/(2.d0*sqrt(rho))+2.d0*a_rho*rho
    END IF
 END IF




end subroutine calc_U_dU


!-------------------------------------------------------------------------------
!
!  Calculate embedding function U(rho) and first +2nd derivative dU(rho)
!
!-------------------------------------------------------------------------------
subroutine calc_U_dU_ddU(rho,U,dU,ddU)
implicit none

!-- Transferred variables
double precision, intent(in)  :: rho
double precision, intent(out) :: U,dU,ddU

!-- Local variables
! <FILL place any local variable definitions here>
 IF(rho.le.1e-10) THEN
    dU=0d0
    U=0.d0
    ddU=0.d0
 ELSE
    IF (rho<pc) THEN
     U=-A*sqrt(rho)-(B/log(2.0d0))*(1.d0-sqrt(rho/pc))*log(2.0d0-rho/pc)+ &
                     a_rho*(rho**2)
     dU=-A/(2.d0*sqrt(rho))-(B/log(2.d0))*( &
               (-1.d0/(2.d0*sqrt(rho*pc)))*log(2.d0-rho/pc)+ &
               (1.d0-sqrt(rho/pc))*(-1.d0/(2.d0*pc-rho)))+2.d0*a_rho*rho
     ddU=A/(4.d0*(rho**(3.d0/2.d0)))-(B/log(2.d0))*( &
                (rho**(-3.d0/2.d0)/(4.d0*sqrt(pc)))*log(2.d0-rho/pc)+ &
                 (1.d0/sqrt(pc*rho))*(1.d0/(2.d0*pc-rho))- &
                  (1.d0-sqrt(rho/pc))*(1.d0/((2.d0*pc-rho)**2)))+2.d0*a_rho
    ELSE
      U=-A*sqrt(rho)+a_rho*(rho**2)
      dU=-A/(2.d0*sqrt(rho))+2.d0*a_rho*rho
      ddU=A/(4.d0*(rho**(3.d0/2.d0)))+2.d0*a_rho
    END IF
 END IF




end subroutine calc_U_dU_ddU




!-------------------------------------------------------------------------------
!
! Compute energy and forces on atoms from the positions.
!
!-------------------------------------------------------------------------------
integer function Compute_Energy_Forces(pkim)
implicit none

!-- Transferred variables
integer(kind=kim_intptr), intent(in)  :: pkim

!-- Local variables
double precision :: Rij(DIM)
double precision :: r,Rsqij,phi,dphi,g,dg,dU,U,dphieff
double precision :: dphii,dUi,Ei,dphij,dUj,Ej
integer :: ier
integer :: i,j,jj,numnei,comp_force,comp_enepot,comp_virial,comp_energy
double precision, allocatable :: rho(:),derU(:)
integer, allocatable, target :: nei1atom_substitute(:)
character*80 :: error_message

!-- KIM variables
integer N;                  pointer(pN,N)
double precision energy;              pointer(penergy,energy)
double precision coordum(DIM,1);      pointer(pcoor,coordum)
double precision forcedum(DIM,1);     pointer(pforce,forcedum)
double precision enepotdum(1);        pointer(penepot,enepotdum)
double precision boxSideLengths(DIM); pointer(pboxSideLengths,boxSideLengths)
double precision Rij_list(DIM,1);     pointer(pRij_list,Rij_list)
integer numContrib;         pointer(pnumContrib,numContrib)
integer nei1atom(1);        pointer(pnei1atom,nei1atom)
integer particleTypes(1);   pointer(pparticleTypes,particleTypes)
double precision virialdum(1);        pointer(pvirial,virialdum)
character*64 NBC_Method;    pointer(pNBC_Method,NBC_Method)
double precision, pointer :: coor(:,:),force(:,:),ene_pot(:),virial_global(:)
integer IterOrLoca
integer HalfOrFull
integer NBC
integer numberContrib
integer idum
integer atom_ret

! Determine neighbor list boundary condition (NBC)
! and half versus full mode:
! *****************************
! * HalfOrFull = 1 -- Half
! *            = 2 -- Full
! *****************************
!
model_cutsq   = model_cutoff**2
!
pNBC_Method = kim_api_get_nbc_method_f(pkim, ier)
if (ier.lt.KIM_STATUS_OK) then
   idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_get_nbc_method_f", ier)
   goto 42
endif
if (index(NBC_Method,"CLUSTER").eq.1) then
   NBC = 0
   HalfOrFull = 1
elseif (index(NBC_Method,"MI_OPBC_H").eq.1) then
   NBC = 1
   HalfOrFull = 1
elseif (index(NBC_Method,"MI_OPBC_F").eq.1) then
   NBC = 1
   HalfOrFull = 2
elseif (index(NBC_Method,"NEIGH_PURE_H").eq.1) then
   NBC = 2
   HalfOrFull = 1
elseif (index(NBC_Method,"NEIGH_PURE_F").eq.1) then
   NBC = 2
   HalfOrFull = 2
elseif (index(NBC_Method,"NEIGH_RVEC_F").eq.1) then
   NBC = 3
   HalfOrFull = 2
else
   ier = KIM_STATUS_FAIL
   idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                 "Unknown NBC method", ier)
   goto 42
endif
call free(pNBC_Method) ! don't forget to release the memory...

! Determine neighbor list handling mode
!
if (NBC.ne.0) then
   !*****************************
   !* IterOrLoca = 1 -- Iterator
   !*            = 2 -- Locator
   !*****************************
   IterOrLoca = kim_api_get_neigh_mode_f(pkim, ier)
   if (ier.lt.KIM_STATUS_OK) then
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_get_neigh_mode_f", ier)
      goto 42
   endif
   if (IterOrLoca.ne.1 .and. IterOrLoca.ne.2) then
      ier = KIM_STATUS_FAIL
      write(error_message,'(a,i1)') &
         'Unsupported IterOrLoca mode = ',IterOrLoca
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    error_message, ier)
      goto 42
   endif
else
   IterOrLoca = 2   ! for CLUSTER NBC
endif

! Check to see if we have been asked to compute the forces, energyperatom,
! energy and virial
!
call kim_api_getm_compute_f(pkim, ier, &
     "energy",         comp_energy, 1, &
     "forces",         comp_force,  1, &
     "particleEnergy", comp_enepot, 1, &
     "virial",         comp_virial, 1)
if (ier.lt.KIM_STATUS_OK) then
   idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_getm_compute_f", ier)
   goto 42
endif

! Unpack data from KIM object
!
call kim_api_getm_data_f(pkim, ier, &
     "numberOfParticles",           pN,              1,                           &
     "particleTypes",               pparticleTypes,  1,                           &
     "coordinates",                 pcoor,           1,                           &
     "numberContributingParticles", pnumContrib,     TRUEFALSE(HalfOrFull.eq.1),  &
     "boxSideLengths",              pboxSideLengths, TRUEFALSE(NBC.eq.1),         &
     "energy",                      penergy,         TRUEFALSE(comp_energy.eq.1), &
     "forces",                      pforce,          TRUEFALSE(comp_force.eq.1),  &
     "particleEnergy",              penepot,         TRUEFALSE(comp_enepot.eq.1), &
     "virial",                      pvirial,         TRUEFALSE(comp_virial.eq.1)  &
     )
if (ier.lt.KIM_STATUS_OK) then
   idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_getm_data_f", ier)
   goto 42
endif

call KIM_to_F90_real_array_2d(coordum,coor,DIM,N)
if (comp_force.eq.1)  call KIM_to_F90_real_array_2d(forcedum,force,DIM,N)
if (comp_enepot.eq.1) call KIM_to_F90_real_array_1d(enepotdum,ene_pot,N)
if (comp_virial.eq.1) call KIM_to_F90_real_array_1d(virialdum,virial_global,6)
if (HalfOrFull.eq.1) then
   if (NBC.ne.0) then ! non-CLUSTER cases
      numberContrib = numContrib
   else               ! CLUSTER cases
      numberContrib = N
   endif
endif

! Check to be sure that the atom types are correct
!
ier = KIM_STATUS_FAIL ! assume an error
do i = 1,N
   if (particleTypes(i).ne.speccode) then
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    "Unexpected species type detected", ier)
      goto 42
   endif
enddo
ier = KIM_STATUS_OK ! everything is ok

! Initialize potential energies, forces, virial term, electron density
!
! Note: that the variable `ene_pot' does not need to be initialized
!       because it's initial value is set during the embedding energy
!       calculation.
!
if (comp_energy.eq.1) energy         = 0.d0
if (comp_force.eq.1)  force(1:DIM,1:N) = 0.d0
if (comp_virial.eq.1) virial_global  = 0.d0
allocate( rho(N) )  ! pair functional electron density
rho(1:N) = 0.d0
if (comp_force.eq.1.or.comp_virial.eq.1) allocate( derU(N) )  ! EAM embedded energy deriv

! Initialize neighbor handling for CLUSTER NBC
!
if (NBC.eq.0) then
   allocate( nei1atom_substitute(N) )
   pnei1atom = loc(nei1atom_substitute)
endif

!
!  Compute energy and forces
!

! Reset iterator if one is being used
!
if (IterOrLoca.eq.1) then
   ier = kim_api_get_neigh_f(pkim,0,0,atom_ret,numnei,pnei1atom,pRij_list)
   if (ier.lt.KIM_STATUS_OK) then
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_get_neigh_f", ier)
      goto 42
   endif
endif

!  Loop over particles in the neighbor list a first time,
!  to compute electron density (=coordination)
!
i = 0
do

   ! Set up neighbor list for next atom for all NBC methods
   !
   call get_current_atom_neighbors(IterOrLoca,HalfOrFull,NBC,N,pkim,      &
                                   i,numnei,pnei1atom,pRij_list,ier)
   if (ier.eq.KIM_STATUS_NEIGH_ITER_PAST_END) exit  ! atom counter incremented past end of list
   if (ier.lt.KIM_STATUS_OK) then
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    "get_current_atom_neighbors", ier)
      goto 42
   endif

   ! Loop over the neighbors of atom i
   !
   do jj = 1, numnei

      j = nei1atom(jj)                            ! get neighbor ID

      ! compute relative position vector
      !
      if (NBC.ne.3) then                          ! all methods except NEIGH_RVEC_F
         Rij(:) = coor(:,j) - coor(:,i)           ! distance vector between i j
      else
         Rij(:) = Rij_list(:,jj)
      endif

      ! apply periodic boundary conditions if required
      !
      if (NBC.eq.1) then
         where ( abs(Rij) .gt. 0.5d0*boxSideLengths )  ! periodic boundary conditions
            Rij = Rij - sign(boxSideLengths,Rij)       ! applied where needed.
         end where                                !
      endif

      ! compute contribution to electron density
      !
      Rsqij = dot_product(Rij,Rij)                ! compute square distance
      if ( Rsqij .lt. model_cutsq ) then          ! particles are interacting?
         r = sqrt(Rsqij)                          ! compute distance
         call calc_g(r,g)                         ! compute electron density
         rho(i) = rho(i) + g                      ! accumulate electron density
         if ((HalfOrFull.eq.1) .and. &
             (j .le. numberContrib)) &            ! HALF mode
            rho(j) = rho(j) + g                   !      (add contrib to j)
      endif

   enddo  ! loop on jj

enddo  ! infinite do loop (terminated by exit statements above)

!  Now that we know the electron densities, calculate embedding part of energy
!  U and its derivative U' (derU)
!
do i = 1,N
   if (comp_force.eq.1.or.comp_virial.eq.1) then
      call calc_U_dU(rho(i),U,dU)                 ! compute embedding energy
                                                  !   and its derivative
      derU(i) = dU                                ! store dU for later use
   else
      call calc_U(rho(i),U)                       ! compute just embedding energy
   endif

   ! accumulate the embedding energy contribution
   !
   ! Assuming U(rho=0) = 0.0d0
   !
   if (comp_enepot.eq.1) then                     ! accumulate embedding energy contribution
      ene_pot(i) = U
   endif
   if (comp_energy.eq.1) then
      energy = energy + U
   endif

   if ((HalfOrFull.eq.1) .and. (i .gt. numberContrib)) exit
enddo

!  Loop over particles in the neighbor list a second time, to compute
!  the forces and complete energy calculation
!

! Reset iterator if one is being used
!
if (IterOrLoca.eq.1) then
   ier = kim_api_get_neigh_f(pkim,0,0,atom_ret,numnei,pnei1atom,pRij_list)
   if (ier.lt.KIM_STATUS_OK) then
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_get_neigh_f", ier)
      goto 42
   endif
endif

i = 0
do

   ! Set up neighbor list for next atom for all NBC methods
   !
   call get_current_atom_neighbors(IterOrLoca,HalfOrFull,NBC,N,pkim,      &
                                   i,numnei,pnei1atom,pRij_list,ier)

   if (ier.eq.KIM_STATUS_NEIGH_ITER_PAST_END) exit  ! atom counter incremented past end of list
   if (ier.lt.KIM_STATUS_OK) then
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    "get_current_atom_neighbors", ier)
      goto 42
   endif


   ! Loop over the neighbors of atom i
   !
   do jj = 1, numnei

      j = nei1atom(jj)                            ! get neighbor ID

      ! compute relative position vector
      !
      if (NBC.ne.3) then                          ! all methods except NEIGH_RVEC_F
         Rij(:) = coor(:,j) - coor(:,i)           ! distance vector between i j
      else
         Rij(:) = Rij_list(:,jj)
      endif

      ! apply periodic boundary conditions if required
      !
      if (NBC.eq.1) then
         where ( abs(Rij) .gt. 0.5d0*boxSideLengths )  ! periodic boundary conditions
            Rij = Rij - sign(boxSideLengths,Rij)       ! applied where needed.
         end where                                !
      endif

      ! compute energy and forces
      !
      Rsqij = dot_product(Rij,Rij)                ! compute square distance
      if ( Rsqij .lt. model_cutsq ) then          ! particles are interacting?

         r = sqrt(Rsqij)                          ! compute distance
         if (comp_force.eq.1.or.comp_virial.eq.1) then
            call calc_phi_dphi(r,phi,dphi)        ! compute pair potential
                                                  !   and it derivative
            call calc_dg(r,dg)                    ! compute elect dens first deriv
            if ((HalfOrFull.eq.1) .and. &
                (j .le. numberContrib)) then      ! HALF mode
               dphii  = 0.5d0*dphi
               dphij  = 0.5d0*dphi
               dUi    = derU(i)*dg
               dUj    = derU(j)*dg
            else                                  ! FULL mode
               dphii  = 0.5d0*dphi
               dphij  = 0.0d0
               dUi    = derU(i)*dg
               dUj    = 0.0d0
            endif
            dphieff = dphii + dphij + dUi + dUj
         else
            call calc_phi(r,phi)                  ! compute just pair potential
         endif
         if ((HalfOrFull.eq.1) .and. &
             (j .le. numberContrib)) then         ! HALF mode
            Ei     = 0.5d0*phi
            Ej     = 0.5d0*phi
         else                                  ! FULL mode
            Ei     = 0.5d0*phi
            Ej     = 0.0d0
         endif

         ! contribution to energy
         !
         if (comp_enepot.eq.1) then
            ene_pot(i) = ene_pot(i) + Ei          ! accumulate energy Ei
            ene_pot(j) = ene_pot(j) + Ej          ! accumulate energy Ej
         endif
         if (comp_energy.eq.1) then
            energy = energy + Ei                  ! accumulate energy
            energy = energy + Ej                  ! accumulate energy
         endif

         ! contribution to virial tensor
         !
         if (comp_virial.eq.1) then
            virial_global(1) = virial_global(1) + Rij(1)*Rij(1)*dphieff/r
            virial_global(2) = virial_global(2) + Rij(2)*Rij(2)*dphieff/r
            virial_global(3) = virial_global(3) + Rij(3)*Rij(3)*dphieff/r
            virial_global(4) = virial_global(4) + Rij(2)*Rij(3)*dphieff/r
            virial_global(5) = virial_global(5) + Rij(1)*Rij(3)*dphieff/r
            virial_global(6) = virial_global(6) + Rij(1)*Rij(2)*dphieff/r
         endif

         ! contribution to forces
         !
         if (comp_force.eq.1) then                        ! Ei contribution
            force(:,i) = force(:,i) + dphieff*Rij/r ! accumulate force on atom i
            force(:,j) = force(:,j) - dphieff*Rij/r ! accumulate force on atom j
         endif

      endif

   enddo  ! loop on jj

enddo  ! infinite do loop (terminated by exit statements above)

! Free temporary storage
!
if (NBC.eq.0) deallocate( nei1atom_substitute )
deallocate( rho )
if (comp_force.eq.1.or.comp_virial.eq.1) deallocate( derU )

! Everything is great
!
ier = KIM_STATUS_OK
42 continue
Compute_Energy_Forces = ier
return

end function Compute_Energy_Forces

!-------------------------------------------------------------------------------
!
! Get list of neighbors for current atom using all NBC methods
!
!-------------------------------------------------------------------------------
subroutine get_current_atom_neighbors(IterOrLoca,HalfOrFull,NBC,N,pkim,      &
                                      atom,numnei,pnei1atom,pRij_list,ier)
implicit none

!-- Transferred variables
integer,                  intent(in)    :: IterOrLoca
integer,                  intent(in)    :: HalfOrFull
integer,                  intent(in)    :: NBC
integer,                  intent(in)    :: N
integer(kind=kim_intptr), intent(in)    :: pkim
integer,                  intent(inout) :: atom
integer,                  intent(out)   :: numnei
integer,                  intent(out)   :: ier
integer nei1atom(1);    pointer(pnei1atom,nei1atom)
double precision Rij_list(DIM,1); pointer(pRij_list,Rij_list)

!-- Local variables
integer atom_ret, jj
integer idum

! Set up neighbor list for next atom for all NBC methods
!
if (IterOrLoca.eq.1) then    ! ITERATOR mode
   ier = kim_api_get_neigh_f(pkim,0,1,atom_ret,numnei, &
                             pnei1atom,pRij_list)
   if (ier.eq.KIM_STATUS_NEIGH_ITER_PAST_END) then
                          ! past end of the list, terminate loop in
      return              ! calling routine
   endif
   if (ier.lt.KIM_STATUS_OK) then     ! some sort of problem, exit
      idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_get_neigh_f", ier)
      return
   endif
   atom = atom_ret

else                         ! LOCATOR mode

   atom = atom + 1
   if (atom.gt.N) then                     ! incremented past end of list,
      ier = KIM_STATUS_NEIGH_ITER_PAST_END ! terminate loop in calling routine
      return
   endif

   if (NBC.eq.0) then ! CLUSTER NBC method
      numnei = N - atom   ! number of neighbors in list atom+1, ..., N
      nei1atom(1:numnei) = (/ (atom+jj, jj = 1,numnei) /)
      ier = KIM_STATUS_OK
   else
      ier = kim_api_get_neigh_f(pkim,1,atom,atom_ret,numnei,pnei1atom,pRij_list)
      if (ier.ne.KIM_STATUS_OK) then ! some sort of problem, exit
         idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                       "kim_api_get_neigh_f", ier)
         ier = KIM_STATUS_FAIL
         return
      endif
   endif
endif

return
end subroutine get_current_atom_neighbors
!-------------------------------------------------------------------------------
!
! Model driver destroy routine
!
!-------------------------------------------------------------------------------
integer function destroy(pkim)
use KIM_API
implicit none

!-- Transferred variables
integer(kind=kim_intptr), intent(in) :: pkim

!-- Local variables

DEALLOCATE(pknotcoeff,pknotpoint,vknotpoint,vknotcoeff)
!-- KIM variables

destroy = KIM_STATUS_OK
return

end function destroy

end module model_driver_pf_cubic_splines

!-------------------------------------------------------------------------------
!
! Model driver initialization routine (REQUIRED)
!
!-------------------------------------------------------------------------------
integer function model_driver_init(pkim, byte_paramfile, nmstrlen, numparamfiles)
use model_driver_pf_cubic_splines
use KIM_API
implicit none

!-- Transferred variables
integer(kind=kim_intptr), intent(in) :: pkim
integer,                  intent(in) :: nmstrlen
integer,                  intent(in) :: numparamfiles
byte,                     intent(in) :: byte_paramfile(nmstrlen*numparamfiles)

!-- Local variables
integer(kind=kim_intptr), parameter :: one=1
character(len=nmstrlen) paramfile_names(numparamfiles)
integer i, j, ier, idum, return_error
character*80 :: error_message
character (LEN=55) :: text
logical ldum

! define variables for all model parameters to be read in

!-- KIM variables
double precision  cutoff;          pointer(pcutoff,cutoff)
character*64 NBC_Method; pointer(pNBC_Method,NBC_Method)

! assume all is well
!
return_error = KIM_STATUS_OK

! generic code to process model parameter file names from byte string
do i=0,numparamfiles-1
   write(paramfile_names(i+1),'(1000a)')  &
        char(byte_paramfile(i*nmstrlen+1: &
                            i*nmstrlen+minloc(abs(byte_paramfile),dim=1)-1))
enddo

! store function pointers in KIM object
call kim_api_setm_method_f(pkim, ier, &
     "compute", one, loc(Compute_Energy_Forces), 1, &
     "destroy", one, loc(destroy),               1)
if (ier.lt.KIM_STATUS_OK) then
   idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_setm_method_f", ier)
   return_error = ier
   goto 1000
endif

! Read in model parameters from parameter file
!
open(10,file=paramfile_names(1),status="old")
    read(10,*,iostat=ier,err=100) model_cutoff
    read(10,*,iostat=ier,err=100) Z
    read(10,*,iostat=ier,err=100) r1, r2, linear_interpolate, interpolate_num, a_rho
    read(10,*,iostat=ier,err=100) nknotp, nknotv
    !
    ! for A,B, and knot functions there are numbers and logicals to skip
    read(10,*,iostat=ier,err=100) idum, A, ldum
    !
    !
    read(10,*,iostat=ier,err=100) idum, B, ldum

    allocate (pknotcoeff(nknotp),pknotpoint(nknotp))
    allocate (vknotcoeff(nknotv),vknotpoint(nknotv))

    i=1
    DO WHILE ((ier==0).AND.(i.LE.nknotp))
     read(10,*,iostat=ier,err=100) idum, pknotcoeff(i), ldum, pknotpoint(i)
     i=i+1
    END DO
    i=1
    DO WHILE ((ier==0).AND.(i.LE.nknotv))
     read(10,*,iostat=ier,err=100) idum, vknotcoeff(i), ldum, vknotpoint(i)
     i=i+1
    END DO
    i=1
    a_inter=0.0d0
    DO WHILE ((ier==0).AND.(i.LE.interpolate_num))
     read(10,*,iostat=ier,err=100) a_inter(i)
     i=i+1
    END DO
    !PRINT *,Z,r1,r2,A,B,a_rho,nknotp,nknotv,a_inter,model_cutoff,interpolate_num
    !PRINT *,vknotpoint,vknotcoeff,pknotpoint,pknotcoeff
    read(10,'(A55)',iostat=ier,err=100) text
close(10)


    write (*,*) 'Using ',TRIM(ADJUSTL(text))
    write (*,*) 'atomic number is: ',Z

goto 200
100 continue
! reading parameters failed
idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                              "Unable to read potential parameters, kimerror = ",KIM_STATUS_FAIL)
return_error = KIM_STATUS_FAIL
goto 1000
!PRINT *,Z,r1,r2,A,B,a_rho,nknotp,nknotv,a_inter,model_cutoff,interpolate_num
!PRINT *,vknotpoint,vknotcoeff,pknotpoint,pknotcoeff

200 continue

! store model cutoff in KIM object
pcutoff =  kim_api_get_data_f(pkim,"cutoff",ier)
if (ier.lt.KIM_STATUS_OK) then
   idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_get_data", ier)
   goto 1000
endif
cutoff = model_cutoff

1000 continue
model_driver_init = return_error
return

end function model_driver_init