utility_forces_numer_deriv.f90

Go to the documentation of this file.

!
! KIM-API: An API for interatomic models
! Copyright (c) 2013--2022, Regents of the University of Minnesota.
! All rights reserved.
!
! Contributors:
!    Ellad B. Tadmor
!    Ryan S. Elliott
!    Stephen M. Whalen
!
! SPDX-License-Identifier: LGPL-2.1-or-later
!
! 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
!

module error
  use, intrinsic :: iso_c_binding
  implicit none
  public

contains
  recursive subroutine my_error(message)
    implicit none
    character(len=*, kind=c_char), intent(in) :: message

    print *, "* Error : ", trim(message)
    stop
  end subroutine my_error

  recursive subroutine my_warning(message)
    implicit none
    character(len=*, kind=c_char), intent(in) :: message

    print *, "* Warning : ", trim(message)
  end subroutine my_warning
end module error

!-------------------------------------------------------------------------------
!
! module mod_neighborlist :
!
!    Module contains type and routines related to neighbor list calculation
!
!-------------------------------------------------------------------------------

module mod_neighborlist

  use, intrinsic :: iso_c_binding

  public get_neigh

  type, bind(c) :: neighobject_type
    real(c_double) :: cutoff
    integer(c_int) :: number_of_particles
    type(c_ptr) :: neighbor_list_pointer
  end type neighobject_type
contains

!-------------------------------------------------------------------------------
!
! get_neigh neighbor list access function
!
! This function implements Locator and Iterator mode
!
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
!
! get_neigh neighbor list access function
!
! This function implements Locator and Iterator mode
!
!-------------------------------------------------------------------------------
  recursive subroutine get_neigh(data_object, number_of_neighbor_lists, &
                                 cutoffs, neighbor_list_index, request, &
                                 numnei, pnei1part, ierr) bind(c)
    use error
    implicit none

    !-- Transferred variables
    type(c_ptr), value, intent(in) :: data_object
    integer(c_int), value, intent(in) :: number_of_neighbor_lists
    real(c_double), intent(in) :: cutoffs(*)
    integer(c_int), value, intent(in) :: neighbor_list_index
    integer(c_int), value, intent(in)  :: request
    integer(c_int), intent(out) :: numnei
    type(c_ptr), intent(out) :: pnei1part
    integer(c_int), intent(out) :: ierr

    !-- Local variables
    integer(c_int) numberofparticles
    type(neighobject_type), pointer :: neighobject
    integer(c_int), pointer :: neighborlist(:, :)

    if (number_of_neighbor_lists > 1) then
      call my_warning("Model requires too many neighbor lists")
      ierr = 1
      return
    end if

    call c_f_pointer(data_object, neighobject)
    numberofparticles = neighobject%number_of_particles
    call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
                     [numberofparticles + 1, numberofparticles])

    if (cutoffs(neighbor_list_index) > neighobject%cutoff) then
      call my_warning("neighbor list cutoff too small for model cutoff")
      ierr = 1
      return
    end if

    if ((request > numberofparticles) .or. (request < 1)) then
      print *, request
      call my_warning("Invalid part ID in get_neigh")
      ierr = 1
      return
    end if

    ! set the returned number of neighbors for the returned part
    numnei = neighborlist(1, request)

    ! set the location for the returned neighbor list
    pnei1part = c_loc(neighborlist(2, request))

    ierr = 0
    return
  end subroutine get_neigh

end module mod_neighborlist

module mod_utilities
  use kim_simulator_headers_module
  implicit none
  public

contains

!-------------------------------------------------------------------------------
!
!  Check if we are compatible with the model
!
!-------------------------------------------------------------------------------
  recursive subroutine check_model_compatibility( &
    compute_arguments_handle, forces_optional, particle_energy_supported, &
    particle_energy_optional, model_is_compatible, ierr)
    use, intrinsic :: iso_c_binding
    use error
    implicit none

    !-- Transferred variables
    type(kim_compute_arguments_handle_type), intent(in) :: &
      compute_arguments_handle
    logical, intent(out) :: forces_optional
    logical, intent(out) :: particle_energy_supported
    logical, intent(out) :: particle_energy_optional
    logical, intent(out) :: model_is_compatible
    integer(c_int), intent(out) :: ierr

    !-- Local variables
    integer(c_int) i
    integer(c_int) number_of_argument_names
    integer(c_int) number_of_callback_names
    type(kim_compute_argument_name_type) argument_name
    type(kim_support_status_type) support_status
    type(kim_compute_callback_name_type) callback_name

    ! assume fail
    model_is_compatible = .false.
    particle_energy_supported = .false.
    particle_energy_optional = .false.
    forces_optional = .false.
    ierr = 0

    ! check arguments
    call kim_get_number_of_compute_argument_names( &
      number_of_argument_names)
    do i = 1, number_of_argument_names
      call kim_get_compute_argument_name(i, argument_name, &
                                         ierr)
      if (ierr /= 0) then
        call my_warning("can't get argument name")
        return
      end if
      call kim_get_argument_support_status( &
        compute_arguments_handle, argument_name, support_status, ierr)
      if (ierr /= 0) then
        call my_warning("can't get argument support_status")
        return
      end if

      ! can only handle energy, particle_energy and forces as required args
      if (support_status == kim_support_status_required) then
        if (.not. ( &
            (argument_name == kim_compute_argument_name_partial_energy) .or. &
            (argument_name == &
             kim_compute_argument_name_partial_particle_energy) .or. &
            (argument_name == kim_compute_argument_name_partial_forces))) then
          call my_warning("unsupported required argument")
          ierr = 0
          return
        end if
      end if

      ! need both energy and forces not "notSupported"
      if ((argument_name == kim_compute_argument_name_partial_energy) .and. &
          (support_status == kim_support_status_not_supported)) then
        call my_warning("model does not support energy")
        ierr = 0
        return
      end if
      if (argument_name == kim_compute_argument_name_partial_forces) then
        if (support_status == kim_support_status_not_supported) then
          call my_warning("model does not support forces")
          ierr = 0
          return
        else if (support_status == kim_support_status_required) then
          forces_optional = .false.
        else if (support_status == kim_support_status_optional) then
          forces_optional = .true.
        else
          call my_warning("unknown support_status for forces")
          ierr = 0
          return
        end if
      end if

      ! check support for particle_energy
      if (argument_name == &
          kim_compute_argument_name_partial_particle_energy) then
        if (support_status == kim_support_status_not_supported) then
          call my_warning("model does not support partial_particle_energy.  &
            &The associated checks will be disabled.")
          particle_energy_supported = .false.
          particle_energy_optional = .false.
        else if (support_status == kim_support_status_required) then
          particle_energy_supported = .true.
          particle_energy_optional = .false.
        else if (support_status == kim_support_status_optional) then
          particle_energy_supported = .true.
          particle_energy_optional = .true.
        else
          call my_warning("unknown support_status for particle energy")
          ierr = 0
          return
        end if
      end if
    end do

    ! check call backs
    call kim_get_number_of_compute_callback_names( &
      number_of_callback_names)
    do i = 1, number_of_callback_names
      call kim_get_compute_callback_name(i, callback_name, &
                                         ierr)
      if (ierr /= 0) then
        call my_warning("can't get call back name")
        return
      end if
      call kim_get_callback_support_status( &
        compute_arguments_handle, callback_name, support_status, ierr)
      if (ierr /= 0) then
        call my_warning("can't get call back support_status")
        return
      end if

      ! cannot handle any "required" call backs
      if (support_status == kim_support_status_required) then
        call my_warning("unsupported required call back")
        ierr = 0
        return
      end if
    end do

    ! got to here, then everything must be OK
    model_is_compatible = .true.
    ierr = 0
    return
  end subroutine check_model_compatibility

  !-----------------------------------------------------------------------------
  !
  !  Get number and identities of particle species supported by
  !  KIM Model `modelname'
  !
  !-----------------------------------------------------------------------------
  recursive subroutine get_model_supported_species( &
    model_handle, max_species, model_species, num_species, ier)
    use, intrinsic :: iso_c_binding
    implicit none

    !-- Transferred variables
    type(kim_model_handle_type), intent(in)   :: model_handle
    integer(c_int), intent(in)   :: max_species
    type(kim_species_name_type), intent(out) :: model_species(max_species)
    integer(c_int), intent(out)  :: num_species
    integer(c_int), intent(out)  :: ier

    !-- Local variables
    integer(c_int) i
    integer(c_int) total_num_species
    type(kim_species_name_type) :: species_name
    integer(c_int) species_is_supported
    integer(c_int) code

    ! Initialize error flag
    ier = 1

    num_species = 0 ! initialize

    call kim_get_number_of_species_names(total_num_species)

    if (total_num_species > max_species) return

    num_species = 0
    do i = 1, total_num_species
      call kim_get_species_name(i, species_name, ier)
      call kim_get_species_support_and_code(model_handle, species_name, &
                                            species_is_supported, code, ier)
      if ((ier == 0) .and. (species_is_supported /= 0)) then
        num_species = num_species + 1
        model_species(num_species) = species_name
      end if
    end do

    ier = 0
    return

  end subroutine get_model_supported_species

  recursive subroutine update_neighborlist(DIM, N, coords, cutoff, cutpad, &
                                           do_update_list, coordsave, &
                                           neighObject, ierr)
    use, intrinsic :: iso_c_binding
    use mod_neighborlist
    implicit none
    integer(c_int), parameter :: cd = c_double ! used for literal constants

    !-- Transferred variables
    integer(c_int), intent(in)    :: DIM
    integer(c_int), intent(in)    :: N
    real(c_double), intent(in)    :: coords(dim, n)
    real(c_double), intent(in)    :: cutoff
    real(c_double), intent(in)    :: cutpad
    logical, intent(inout) :: do_update_list
    real(c_double), intent(inout) :: coordsave(dim, n)
    type(neighobject_type), intent(inout) :: neighObject
    integer(c_int), intent(out)   :: ierr

    !-- Local variables
    real(c_double) disp, disp1, disp2, cutrange, dispvec(dim)
    integer(c_int) i

    ! Initialize error code
    !
    ierr = 0

    ! Update neighbor lists if necessary
    !
    if (.not. do_update_list) then   ! if update not requested

      ! check whether a neighbor list update is necessary even if it hasn't been
      ! requested using the "two max sum" criterion
      disp1 = 0.0_cd
      disp2 = 0.0_cd
      do i = 1, n
        dispvec(1:dim) = coords(1:dim, i) - coordsave(1:dim, i)
        disp = sqrt(dot_product(dispvec, dispvec))
        if (disp >= disp1) then        !  1st position taken
          disp2 = disp1               !  push current 1st into 2nd place
          disp1 = disp                !  and put this one into current 1st
        else if (disp >= disp2) then   !  2nd position taken
          disp2 = disp
        end if
      end do
      do_update_list = (disp1 + disp2 > cutpad)

    end if

    if (do_update_list) then

      ! save current coordinates
      coordsave(1:dim, 1:n) = coords(1:dim, 1:n)

      ! compute neighbor lists
      cutrange = cutoff + cutpad
      call neigh_pure_cluster_neighborlist(.false., n, coords, cutrange, &
                                           neighobject)

      ! neighbor list uptodate, no need to compute again for now
      do_update_list = .false.
    end if

    return

  end subroutine update_neighborlist

  !-----------------------------------------------------------------------------
  !
  ! NEIGH_PURE_cluster_neighborlist
  !
  !-----------------------------------------------------------------------------
  recursive subroutine neigh_pure_cluster_neighborlist( &
    half, numberOfParticles, coords, cutoff, neighObject)
    use, intrinsic :: iso_c_binding
    use mod_neighborlist
    implicit none

    !-- Transferred variables
    logical, intent(in)            :: half
    integer(c_int), intent(in)            :: numberOfParticles
    real(c_double), dimension(3, numberOfParticles), &
      intent(in)            :: coords
    real(c_double), intent(in)            :: cutoff
    type(neighobject_type), intent(inout) :: neighObject

    !-- Local variables
    integer(c_int), pointer :: neighborList(:, :)
    integer(c_int) i, j, a
    real(c_double) dx(3)
    real(c_double) r2
    real(c_double) cutoff2

    call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
                     [numberofparticles + 1, numberofparticles])

    neighobject%cutoff = cutoff

    cutoff2 = cutoff**2

    do i = 1, numberofparticles
      a = 1
      do j = 1, numberofparticles
        dx(:) = coords(:, j) - coords(:, i)
        r2 = dot_product(dx, dx)
        if (r2 <= cutoff2) then
          ! part j is a neighbor of part i
          if ((j > i) .OR. ((.not. half) .AND. (i /= j))) then
            a = a + 1
            neighborlist(a, i) = j
          end if
        end if
      end do
      ! part i has a-1 neighbors
      neighborlist(1, i) = a - 1
    end do

    return

  end subroutine neigh_pure_cluster_neighborlist

  !-----------------------------------------------------------------------------
  !
  ! create_FCC_configuration subroutine
  !
  !  creates a cubic configuration of FCC particles with lattice spacing
  !  `FCCspacing' and `nCellsPerSide' cells along each direction.
  !
  !  With periodic==.true. this will result in a total number of particles equal
  !  to 4*(nCellsPerSide)**3 + 6*(nCellsPerSide)**2 + 3*(nCellsPerSide) + 1
  !
  !  With periodic==.false. this will result in a total number of particles equal
  !  to 4*(nCellsPerSide)**3
  !
  !  Returns the Id of the particle situated in the middle of the configuration
  !  (this particle will have the most neighbors.)
  !
  !-----------------------------------------------------------------------------
  recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, &
                                                periodic, coords, MiddlePartId)
    use, intrinsic :: iso_c_binding
    implicit none
    integer(c_int), parameter :: cd = c_double ! used for literal constants

    !-- Transferred variables
    real(c_double), intent(in)  :: FCCspacing
    integer(c_int), intent(in)  :: nCellsPerSide
    logical, intent(in)  :: periodic
    real(c_double), intent(out) :: coords(3, *)
    integer(c_int), intent(out) :: MiddlePartId
    !
    ! cluster setup variables
    !
    real(c_double) FCCshifts(3, 4)
    real(c_double) latVec(3)
    integer(c_int) a, i, j, k, m

    ! Create a cubic FCC cluster
    !
    fccshifts(1, 1) = 0.0_cd
    fccshifts(2, 1) = 0.0_cd
    fccshifts(3, 1) = 0.0_cd
    fccshifts(1, 2) = 0.5_cd * fccspacing
    fccshifts(2, 2) = 0.5_cd * fccspacing
    fccshifts(3, 2) = 0.0_cd
    fccshifts(1, 3) = 0.5_cd * fccspacing
    fccshifts(2, 3) = 0.0_cd
    fccshifts(3, 3) = 0.5_cd * fccspacing
    fccshifts(1, 4) = 0.0_cd
    fccshifts(2, 4) = 0.5_cd * fccspacing
    fccshifts(3, 4) = 0.5_cd * fccspacing

    middlepartid = 1 ! Always put middle particle as #1
    a = 1            ! leave space for middle particle as particle #1
    do i = 1, ncellsperside
      latvec(1) = (i - 1) * fccspacing
      do j = 1, ncellsperside
        latvec(2) = (j - 1) * fccspacing
        do k = 1, ncellsperside
          latvec(3) = (k - 1) * fccspacing
          do m = 1, 4
            a = a + 1
            coords(:, a) = latvec + fccshifts(:, m)
            if ((i == ncellsperside / 2 + 1) &
                .and. (j == ncellsperside / 2 + 1) &
                .and. (k == ncellsperside / 2 + 1) &
                .and. (m == 1)) &
              then
              ! put middle particle as #1
              coords(:, 1) = latvec + fccshifts(:, m)
              a = a - 1
            end if
          end do
        end do
        if (.not. periodic) then
          ! Add in the remaining three faces
          ! pos-x face
          latvec(1) = ncellsperside * fccspacing
          latvec(2) = (i - 1) * fccspacing
          latvec(3) = (j - 1) * fccspacing
          a = a + 1; coords(:, a) = latvec
          a = a + 1; coords(:, a) = latvec + fccshifts(:, 4)
          ! pos-y face
          latvec(1) = (i - 1) * fccspacing
          latvec(2) = ncellsperside * fccspacing
          latvec(3) = (j - 1) * fccspacing
          a = a + 1; coords(:, a) = latvec
          a = a + 1; coords(:, a) = latvec + fccshifts(:, 3)
          ! pos-z face
          latvec(1) = (i - 1) * fccspacing
          latvec(2) = (j - 1) * fccspacing
          latvec(3) = ncellsperside * fccspacing
          a = a + 1; coords(:, a) = latvec
          a = a + 1; coords(:, a) = latvec + fccshifts(:, 2)
        end if
      end do
      if (.not. periodic) then
        ! Add in the remaining three edges
        latvec(1) = (i - 1) * fccspacing
        latvec(2) = ncellsperside * fccspacing
        latvec(3) = ncellsperside * fccspacing
        a = a + 1; coords(:, a) = latvec
        latvec(1) = ncellsperside * fccspacing
        latvec(2) = (i - 1) * fccspacing
        latvec(3) = ncellsperside * fccspacing
        a = a + 1; coords(:, a) = latvec
        latvec(1) = ncellsperside * fccspacing
        latvec(2) = ncellsperside * fccspacing
        latvec(3) = (i - 1) * fccspacing
        a = a + 1; coords(:, a) = latvec
      end if
    end do
    if (.not. periodic) then
      ! Add in the remaining corner
      a = a + 1; coords(:, a) = ncellsperside * fccspacing
    end if

    return

  end subroutine create_fcc_configuration

  recursive subroutine compute_numer_deriv( &
    partnum, dir, model_handle, compute_arguments_handle, DIM, N, coords, &
    cutoff, cutpad, energy, do_update_list, coordsave, neighObject, deriv, &
    deriv_err, ierr)
    use, intrinsic :: iso_c_binding
    use error
    use mod_neighborlist
    implicit none
    integer(c_int), parameter :: cd = c_double ! used for literal constants

    !--Transferred variables
    integer(c_int), intent(in)    :: partnum
    integer(c_int), intent(in)    :: dir
    type(kim_model_handle_type), intent(in) :: model_handle
    type(kim_compute_arguments_handle_type), intent(in) :: &
      compute_arguments_handle
    integer(c_int), intent(in)    :: DIM
    integer(c_int), intent(in)    :: N
    real(c_double), intent(inout) :: coords(dim, n)
    real(c_double), intent(in)    :: cutoff
    real(c_double), intent(in)    :: cutpad
    real(c_double), intent(inout) :: energy
    logical, intent(inout) :: do_update_list
    real(c_double), intent(inout) :: coordsave(dim, n)
    type(neighobject_type), intent(inout) :: neighObject
    real(c_double), intent(out)   :: deriv
    real(c_double), intent(out)   :: deriv_err
    integer(c_int), intent(out)   :: ierr

    !-- Local variables
    real(c_double), parameter :: eps_init = 1.e-6_cd
    integer(c_int), parameter :: number_eps_levels = 15
    real(c_double) eps, deriv_last, deriv_err_last
    integer(c_int) i

    ! Initialize error flag
    ierr = 0

    deriv_last = 0.0_cd ! initialize

    ! Outer loop of Ridders' method for computing numerical derivative
    !
    eps = eps_init
    deriv_err_last = huge(1.0_cd)
    do i = 1, number_eps_levels
      deriv = dfridr(eps, deriv_err)
      if (ierr /= 0) then
        call my_error("compute_numer_deriv")
      end if
      if (deriv_err > deriv_err_last) then
        deriv = deriv_last
        deriv_err = deriv_err_last
        exit
      end if
      eps = eps * 10.0_cd
      deriv_last = deriv
      deriv_err_last = deriv_err
    end do

    return

  contains

    !----------------------------------------------------------------------------
    !
    ! Compute numerical derivative using Ridders' method
    !
    ! Based on code from Numerical Recipes, Press et al., Second Ed., Cambridge,
    ! 1992
    !
    ! Ref: Ridders, C. J. F., "Two algorithms for the calculation of F'(x)=D",
    !      Advances in Engineering Software, Vol. 4, no. 2, pp. 75-76, 1982.
    !
    !
    ! Returns the gradient grad() of a KIM-compliant interatomic model at the
    ! current configuration by Ridders' method of polynomial extrapolation.
    ! An estimate for the error in each component of the gradient is returned in
    ! grad_err().
    !
    !----------------------------------------------------------------------------
    real(c_double) recursive function dfridr(h, err)
      implicit none

      !-- Transferred variables
      real(c_double), intent(inout) :: h
      real(c_double), intent(out)   :: err

      !-- Local variables
      integer(c_int), parameter :: NTAB = 10     ! Maximum size of tableau
      real(c_double), parameter :: CON = 1.4_cd  ! Stepsize incr. by CON at each iter
      real(c_double), parameter :: CON2 = con * con
      real(c_double), parameter :: BIG = huge(1.0_cd)
      real(c_double), parameter :: SAFE = 2.0_cd ! Returns when error is SAFE worse
      ! than the best so far
      integer(c_int) i, j
      real(c_double) errt, fac, hh, a(ntab, ntab), fp, fm, coordorig

      dfridr = 0.0_cd ! initialize
      err = big ! initialize

      if (abs(h) <= tiny(0.0_cd)) then  ! avoid division by zero
        ierr = 1
        return
      end if

      hh = h
      coordorig = coords(dir, partnum)
      coords(dir, partnum) = coordorig + hh
      call update_neighborlist(dim, n, coords, cutoff, cutpad, &
                               do_update_list, coordsave, &
                               neighobject, ierr)
      call kim_compute(model_handle, compute_arguments_handle, ierr)
      if (ierr /= 0) then
        call my_error("kim_api_model_compute")
      end if
      fp = energy
      coords(dir, partnum) = coordorig - hh
      call update_neighborlist(dim, n, coords, cutoff, cutpad, &
                               do_update_list, coordsave, &
                               neighobject, ierr)
      call kim_compute(model_handle, compute_arguments_handle, ierr)
      if (ierr /= 0) then
        call my_error("kim_api_model_compute")
      end if
      fm = energy
      coords(dir, partnum) = coordorig
      call update_neighborlist(dim, n, coords, cutoff, cutpad, &
                               do_update_list, coordsave, &
                               neighobject, ierr)
      a(1, 1) = (fp - fm) / (2.0_cd * hh)
      ! successive columns in the Neville tableau will go to smaller step sizes
      ! and higher orders of extrapolation
      do i = 2, ntab
        ! try new, smaller step size
        hh = hh / con
        coords(dir, partnum) = coordorig + hh
        call update_neighborlist(dim, n, coords, cutoff, cutpad, &
                                 do_update_list, coordsave, &
                                 neighobject, ierr)
        call kim_compute(model_handle, compute_arguments_handle, ierr)
        if (ierr /= 0) then
          call my_error("kim_api_model_compute")
        end if
        fp = energy
        coords(dir, partnum) = coordorig - hh
        call update_neighborlist(dim, n, coords, cutoff, cutpad, &
                                 do_update_list, coordsave, &
                                 neighobject, ierr)
        call kim_compute(model_handle, compute_arguments_handle, ierr)
        if (ierr /= 0) then
          call my_error("kim_api_model_compute")
        end if
        fm = energy
        coords(dir, partnum) = coordorig
        call update_neighborlist(dim, n, coords, cutoff, cutpad, &
                                 do_update_list, coordsave, &
                                 neighobject, ierr)
        a(1, i) = (fp - fm) / (2.0_cd * hh)
        fac = con2
        ! compute extrapolations of various orders, requiring no new function
        ! evaluations
        do j = 2, i
          a(j, i) = (a(j - 1, i) * fac - a(j - 1, i - 1)) / (fac - 1.0_cd)
          fac = con2 * fac
          ! The error strategy is to compute each new extrapolation to one order
          ! lower, both at the present step size and the previous one.
          errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
          if (errt <= err) then ! if error is decreased, save the improved answer
            err = errt
            dfridr = a(j, i)
          end if
        end do
        if (abs(a(i, i) - a(i - 1, i - 1)) >= safe * err) return ! if higher order is worse
        ! by significant factor
        ! `SAFE', then quit early.
      end do
      return
    end function dfridr

  end subroutine compute_numer_deriv

end module mod_utilities

!*******************************************************************************
!**
!**  PROGRAM vc_forces_numer_deriv
!**
!**  KIM compliant program to perform numerical derivative check on a model
!**
!*******************************************************************************

!-------------------------------------------------------------------------------
!
! Main program
!
!-------------------------------------------------------------------------------
program vc_forces_numer_deriv
  use, intrinsic :: iso_c_binding
  use error
  use mod_neighborlist
  use mod_utilities
  implicit none
  integer(c_int), parameter :: cd = c_double ! used for literal constants

  integer(c_int), parameter :: nCellsPerSide = 2
  integer(c_int), parameter :: DIM = 3
  real(c_double), parameter :: cutpad = 0.75_cd
  integer(c_int), parameter :: max_species = 200 ! most species supported
  real(c_double), parameter :: eps_prec = epsilon(1.0_cd)
  real(c_double) FCCspacing

  integer(c_int), parameter :: N = 4 * (ncellsperside)**3 &
                               + 6 * (ncellsperside)**2 &
                               + 3 * (ncellsperside) + 1
  real(c_double), allocatable          :: forces_num(:, :)
  real(c_double), allocatable          :: forces_num_err(:, :)
  type(kim_species_name_type) :: model_species(max_species)
  integer(c_int), target               :: num_species
  character(len=5, kind=c_char)        :: passfail
  real(c_double)                       :: forcediff
  real(c_double)                       :: forcediff_sumsq
  real(c_double)                       :: weight
  real(c_double)                       :: weight_sum
  real(c_double)                       :: alpha
  real(c_double)                       :: term
  real(c_double)                       :: term_max
  real(c_double), allocatable          :: cluster_coords(:, :)
  real(c_double), allocatable          :: cluster_disps(:, :)
  type(kim_species_name_type), allocatable :: cluster_species(:)
  integer(c_int) I, J, Imax, Jmax, species_code
  integer(c_int) seed_size
  integer(c_int), allocatable :: seed(:)

  !
  ! neighbor list
  !
  type(neighobject_type), target :: neighObject
  integer(c_int), allocatable, target :: neighborList(:, :)
  real(c_double), allocatable :: coordsave(:, :)
  logical do_update_list

  !
  ! KIM variables
  !
  character(len=256, kind=c_char) :: testname = "vc_forces_numer_deriv"
  character(len=256, kind=c_char) :: modelname

  type(kim_model_handle_type) :: model_handle
  type(kim_compute_arguments_handle_type) :: compute_arguments_handle
  integer(c_int) ierr, ierr2
  integer(c_int) species_is_supported
  integer(c_int), target :: numberOfParticles
  integer(c_int), target :: particleSpeciesCodes(n)
  integer(c_int), target :: particleContributing(n)
  real(c_double) :: influence_distance
  integer(c_int) :: number_of_neighbor_lists
  real(c_double), allocatable :: cutoffs(:)
  integer(c_int), allocatable :: &
    model_will_not_request_neighbors_of_noncontributing_particles(:)
  real(c_double) :: cutoff
  real(c_double), target :: energy
  real(c_double), target :: particle_energy(n)
  real(c_double), target :: particle_energy_sum
  real(c_double), target :: coords(3, n)
  real(c_double), target :: forces_kim(3, n)
  real(c_double) :: forces(3, n)
  integer(c_int) middleDum
  logical doing_non_contributing
  logical particle_energy_supported
  logical particle_energy_optional
  logical noncontrib_particle_energy_nonzero
  logical forces_optional
  logical model_is_compatible
  integer(c_int) number_of_model_routine_names
  type(kim_model_routine_name_type) model_routine_name
  integer(c_int) present
  integer(c_int) required
  integer(c_int) requested_units_accepted
  real(c_double) rnd, deriv, deriv_err
  real(c_double), pointer :: null_pointer

  nullify (null_pointer)

  doing_non_contributing = .false. ! .true. on 2nd pass through main program

  numberofparticles = n

  term_max = 0.0_cd ! initialize

  ! Initialize error flag
  ierr = 0

  ! Initialize seed for random number generator
  !
  ! NOTE: Here, we set the seed to a fixed value for reproducibility
  call random_seed(size=seed_size) ! Get seed size for this CPU
  allocate (seed(seed_size))
  seed(:) = 13
  call random_seed(put=seed)
  deallocate (seed)

  ! Get KIM Model name to use
  print '("Please enter a valid KIM model name: ")'
  read (*, *) modelname

  ! Print output header
  !
  print *
  print *, 'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES'
  print *
  print '(120(''-''))'
  print '("This is Test          : ",A)', trim(testname)
  print '("Results for KIM Model : ",A)', trim(modelname)

  ! Create empty KIM object
  !
  call kim_model_create(kim_numbering_one_based, &
                        kim_length_unit_a, &
                        kim_energy_unit_ev, &
                        kim_charge_unit_e, &
                        kim_temperature_unit_k, &
                        kim_time_unit_ps, &
                        trim(modelname), &
                        requested_units_accepted, &
                        model_handle, ierr)
  if (ierr /= 0) then
    call my_error("kim_api_create")
  end if

  ! check that we are compatible
  if (requested_units_accepted == 0) then
    call my_error("Must adapt to model units")
  end if

  ! check that we know about all required routines
  call kim_get_number_of_model_routine_names(number_of_model_routine_names)
  do i = 1, number_of_model_routine_names
    call kim_get_model_routine_name(i, model_routine_name, ierr)
    if (ierr /= 0) then
      call my_error("kim_get_model_routine_name")
    end if
    call kim_is_routine_present(model_handle, model_routine_name, present, &
                                required, ierr)
    if (ierr /= 0) then
      call my_error("kim_is_routine_present")
    end if

    if ((present == 1) .and. (required == 1)) then
      if (.not. ((model_routine_name == kim_model_routine_name_create) &
                 .or. (model_routine_name == &
                       kim_model_routine_name_compute_arguments_create) &
                 .or. (model_routine_name == kim_model_routine_name_compute) &
                 .or. (model_routine_name == kim_model_routine_name_refresh) &
                 .or. (model_routine_name == &
                       kim_model_routine_name_compute_arguments_destroy) &
                 .or. (model_routine_name == kim_model_routine_name_destroy))) &
        then
        call my_error("Unknown required ModelRoutineName found.")
      end if
    end if
  end do

  ! create compute_arguments object
  call kim_compute_arguments_create(model_handle, &
                                    compute_arguments_handle, ierr)
  if (ierr /= 0) then
    call my_error("kim_model_compute_arguments_create")
  end if

  call check_model_compatibility(compute_arguments_handle, forces_optional, &
                                 particle_energy_supported, &
                                 particle_energy_optional, &
                                 model_is_compatible, ierr)
  if (ierr /= 0) then
    call my_error("error checking compatibility")
  end if
  if (.not. model_is_compatible) then
    call my_error("incompatibility reported by check_model_compatibility")
  end if

  ! Get list of particle species supported by the model
  !
  call get_model_supported_species(model_handle, max_species, model_species, &
                                   num_species, ierr)
  if (ierr /= 0) then
    call my_error("Get_Model_Supported_Species")
  end if
  ! Setup random cluster
  !
  allocate (cluster_coords(3, n), cluster_disps(3, n), cluster_species(n))
  do i = 1, n
    call random_number(rnd)  ! return random number between 0 and 1
    species_code = 1 + int(rnd * num_species)
    cluster_species(i) = model_species(species_code)
  end do
  fccspacing = 1.0_cd  ! initially generate an FCC cluster with lattice
  ! spacing equal to one. This is scaled below based
  ! on the cutoff radius of the model.
  call create_fcc_configuration(fccspacing, ncellsperside, .false., &
                                cluster_coords, middledum)
  ! Generate random displacements for all particles
  !
  do i = 1, n
    do j = 1, dim
      call random_number(rnd)  ! return random number between 0 and 1
      cluster_disps(j, i) = 0.1_cd * (rnd - 0.5_cd)
    end do
  end do

  ! register memory with the KIM system
  ierr = 0
  call kim_set_argument_pointer( &
    compute_arguments_handle, &
    kim_compute_argument_name_number_of_particles, &
    numberofparticles, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer( &
    compute_arguments_handle, &
    kim_compute_argument_name_particle_species_codes, &
    particlespeciescodes, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer( &
    compute_arguments_handle, &
    kim_compute_argument_name_particle_contributing, &
    particlecontributing, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer( &
    compute_arguments_handle, &
    kim_compute_argument_name_coordinates, coords, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer( &
    compute_arguments_handle, &
    kim_compute_argument_name_partial_energy, energy, ierr2)
  ierr = ierr + ierr2
  if (particle_energy_supported) then
    call kim_set_argument_pointer( &
      compute_arguments_handle, &
      kim_compute_argument_name_partial_particle_energy, particle_energy, ierr2)
    ierr = ierr + ierr2
  end if
  call kim_set_argument_pointer( &
    compute_arguments_handle, &
    kim_compute_argument_name_partial_forces, forces_kim, ierr2)
  ierr = ierr + ierr2
  if (ierr /= 0) then
    call my_error("set_argument_pointer")
  end if

  ! Allocate storage for neighbor lists and
  ! store pointers to neighbor list object and access function
  !
  allocate (coordsave(dim, n))
  allocate (neighborlist(n + 1, n))
  neighobject%neighbor_list_pointer = c_loc(neighborlist)
  neighobject%number_of_particles = n
  allocate (forces_num(dim, n), forces_num_err(dim, n))

  ! Set pointer in KIM object to neighbor list routine and object
  !
  call kim_set_callback_pointer( &
    compute_arguments_handle, &
    kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
    c_funloc(get_neigh), c_loc(neighobject), ierr)
  if (ierr /= 0) then
    call my_error("set_callback_pointer")
  end if

  call kim_get_influence_distance(model_handle, influence_distance)
  call kim_get_number_of_neighbor_lists(model_handle, &
                                        number_of_neighbor_lists)
  allocate (cutoffs(number_of_neighbor_lists), &
            model_will_not_request_neighbors_of_noncontributing_particles( &
            number_of_neighbor_lists))
  call kim_get_neighbor_list_values( &
    model_handle, cutoffs, &
    model_will_not_request_neighbors_of_noncontributing_particles, ierr)
  if (ierr /= 0) then
    call my_error("get_neighbor_list_values")
  end if
  cutoff = maxval(cutoffs)

  ! Scale reference FCC configuration based on cutoff radius.
  fccspacing = 0.75_cd * cutoff ! set the FCC spacing to a fraction
  ! of the cutoff radius
  do i = 1, n
    cluster_coords(:, i) = fccspacing * cluster_coords(:, i)
  end do
  print '("Using FCC lattice parameter: ",f12.5)', fccspacing
  print *

1000 continue ! Start of configuration setup

  if (doing_non_contributing) then
    ! Turn particle energy computation back on, if possible
    if (particle_energy_optional) then
      call kim_set_argument_pointer( &
        compute_arguments_handle, &
        kim_compute_argument_name_partial_particle_energy, &
        particle_energy, ierr)
      if (ierr /= 0) then
        call my_error("set_argument_pointer")
      end if
    end if

    ! Turn force computation back on, if possible
    !
    if (forces_optional) then
      call kim_set_argument_pointer( &
        compute_arguments_handle, kim_compute_argument_name_partial_forces, &
        forces_kim, ierr)
      if (ierr /= 0) then
        call my_error("set_argument_pointer")
      end if
    end if
  end if

  do i = 1, n
    call kim_get_species_support_and_code( &
      model_handle, cluster_species(i), species_is_supported, &
      particlespeciescodes(i), ierr)
  end do
  if (ierr /= 0) then
    call my_error("kim_api_get_species_code")
  end if
  if (.not. doing_non_contributing) then
    do i = 1, n
      particlecontributing(i) = 1  ! all particle contribute
    end do
  else
    do i = 1, n
      ! Random collection of contributing atoms
      call random_number(rnd)  ! return random number between 0 and 1
      if (rnd > 0.5) then
        particlecontributing(i) = 1
      else
        particlecontributing(i) = 0
      end if
    end do
  end if
  do i = 1, n
    coords(:, i) = cluster_coords(:, i) + cluster_disps(:, i)
  end do

  ! Compute neighbor lists
  !
  do_update_list = .true.
  call update_neighborlist(dim, n, coords, cutoff, cutpad, &
                           do_update_list, coordsave, &
                           neighobject, ierr)
  if (ierr /= 0) then
    call my_error("update_neighborlist")
  end if

  ! Call model compute to get forces (gradient)
  !
  call kim_compute(model_handle, compute_arguments_handle, ierr)
  if (ierr /= 0) then
    call my_error("kim_api_model_compute")
  end if

  ! Copy forces in case model will overwrite forces_kim below
  !
  forces = forces_kim

  ! Add up particle_energy to compare with energy
  noncontrib_particle_energy_nonzero = .false.
  if (particle_energy_supported) then
    particle_energy_sum = 0.0_cd
    do i = 1, n
      if (particlecontributing(i) == 0) then
        if (abs(particle_energy(i)) > epsilon(0.0_cd)) then
          noncontrib_particle_energy_nonzero = .true.
        end if
      else
        particle_energy_sum = particle_energy_sum + particle_energy(i)
      end if
    end do
  end if

  ! Print results to screen
  !
  print '(41(''=''))'
  if (.not. doing_non_contributing) then
    print '("Configuration with all contributing particles")'
  else
    print '("Configuration with some non-contributing particles")'
  end if
  if (particle_energy_supported) then
    print '(A25,2X,A25,2X,A15)', "Energy", "Sum Contrib. Energies", "Diff"
    print '(ES25.15,2X,ES25.15,2X,ES15.5)', energy, particle_energy_sum, &
      energy - particle_energy_sum
    if (noncontrib_particle_energy_nonzero) then
      call my_error( &
        "Some non-contributing particles have non-zero &
        &partial_particle_energy")
    end if
  else
    print '("Energy = ",ES25.15)', energy
  end if
  print '(41(''=''))'
  print *

  ! Turn off particle energy computation, if possible
  if (particle_energy_optional) then
    call kim_set_argument_pointer( &
      compute_arguments_handle, &
      kim_compute_argument_name_partial_particle_energy, &
      null_pointer, ierr)
    if (ierr /= 0) then
      call my_error("set_argument_pointer")
    end if
  end if

  ! Turn off force computation, if possible
  !
  if (forces_optional) then
    call kim_set_argument_pointer( &
      compute_arguments_handle, kim_compute_argument_name_partial_forces, &
      null_pointer, ierr)
    if (ierr /= 0) then
      call my_error("set_argument_pointer")
    end if
  end if

  ! Compute gradient using numerical differentiation
  !
  do i = 1, n
    do j = 1, dim
      call compute_numer_deriv(i, j, model_handle, compute_arguments_handle, &
                               dim, n, coords, cutoff, &
                               cutpad, energy, do_update_list, &
                               coordsave, neighobject, deriv, deriv_err, ierr)
      if (ierr /= 0) then
        call my_error("compute_numer_deriv")
      end if
      forces_num(j, i) = -deriv
      forces_num_err(j, i) = deriv_err
    end do
  end do

  ! Continue printing results to screen
  !
  print '(A6,2X,A7,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)', "Part", "Contrib", &
    "Spec", "Dir", "Force_model", "Force_numer", "Force diff", "pred error", &
    "weight", "stat"
  forcediff_sumsq = 0.0_cd
  weight_sum = 0.0_cd
  do i = 1, n
    do j = 1, dim
      forcediff = abs(forces(j, i) - forces_num(j, i))
      if (forcediff < forces_num_err(j, i)) then
        passfail = "ideal"
      else
        passfail = "     "
      end if
      weight = max(abs(forces_num(j, i)), eps_prec) / &
               max(abs(forces_num_err(j, i)), eps_prec)
      term = weight * forcediff**2
      if (term > term_max) then
        term_max = term
        imax = i
        jmax = j
      end if
      forcediff_sumsq = forcediff_sumsq + term
      weight_sum = weight_sum + weight
      if (j == 1) then
        print '(I6,2X,I7,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
          i, particlecontributing(i), &
          particlespeciescodes(i), j, forces(j, i), &
          forces_num(j, i), forcediff, &
          forces_num_err(j, i), weight, passfail
      else
        print '(23X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
          j, forces(j, i), forces_num(j, i), &
          forcediff, forces_num_err(j, i), weight, passfail
      end if
    end do
    print *
  end do
  alpha = sqrt(forcediff_sumsq / weight_sum) / dble(dim * n)
  print *
  print '("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," &
        &(units of force)")', alpha
  print *
  print '(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,'// &
    ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
    imax, jmax, abs(forces(jmax, imax) - forces_num(jmax, imax)), &
    abs(forces(jmax, imax) - forces_num(jmax, imax)) / abs(forces(jmax, imax))

  if (.not. doing_non_contributing) then
    doing_non_contributing = .true.
    print *
    print *
    goto 1000
  end if

  ! Free temporary storage
  !
  deallocate (forces_num)
  deallocate (forces_num_err)
  deallocate (neighborlist)
  deallocate (coordsave)
  deallocate (cutoffs)
  deallocate (model_will_not_request_neighbors_of_noncontributing_particles)

  call kim_compute_arguments_destroy(model_handle, &
                                     compute_arguments_handle, ierr)
  if (ierr /= 0) then
    call my_error("kim_model_compute_arguments_destroy")
  end if
  call kim_model_destroy(model_handle)

  ! Print output footer
  !
  print *
  print '(120(''-''))'

  ! Free cluster storage
  !
  deallocate (cluster_coords, cluster_disps, cluster_species)

  stop

end program vc_forces_numer_deriv