utility_forces_numer_deriv.f90

Go to the documentation of this file.

!
! 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--2019, Regents of the University of Minnesota.
! All rights reserved.
!
! Contributors:
!    Ellad B. Tadmor
!    Ryan S. Elliott
!    Stephen M. Whalen
!


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(number_of_neighbor_lists)
  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(:,:)

  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
  endif

  if ( (request.gt.numberofparticles) .or. (request.lt.1)) then
    print *, request
    call my_warning("Invalid part ID in get_neigh")
    ierr = 1
    return
  endif

  ! 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, 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) :: 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.
  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 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_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
  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

call kim_get_number_of_species_names(total_num_species)

if (total_num_species .gt. 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 .ne. 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
      endif
   enddo
   do_update_list = ( disp1 + disp2 > cutpad )

endif

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.
endif

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.le.cutoff2) then
           ! part j is a neighbor of part i
           if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) ) then
               a = a+1
               neighborlist(a,i) = j
           endif
        endif
     enddo
     ! part i has a-1 neighbors
     neighborlist(1,i) = a-1
  enddo

  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.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
                   (k.eq.ncellsperside/2+1) .and. (m.eq.1)) then
                 coords(:,1) = latvec + fccshifts(:,m) ! put middle particle as #1
                 a = a - 1
              endif
           enddo
        enddo
        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)
         endif
     enddo
     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
      endif
  enddo
  if (.not. periodic) then
      ! Add in the remaining corner
      a = a+1; coords(:,a) = ncellsperside*fccspacing
  endif

  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")
   endif
   if (deriv_err>deriv_err_last) then
      deriv  = deriv_last
      deriv_err = deriv_err_last
      exit
   endif
   eps = eps*10.0_cd
   deriv_last  = deriv
   deriv_err_last = deriv_err
enddo

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

   if (abs(h).le.tiny(0.0_cd)) then  ! avoid division by zero
      ierr = 1
      return
   endif

   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")
   endif
   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")
   endif
   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)
   err=big
   ! 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")
      endif
      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")
      endif
      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.le.err) then ! if error is decreased, save the improved answer
            err = errt
            dfridr=a(j,i)
         endif
      enddo
      if (abs(a(i,i)-a(i-1,i-1)).ge.safe*err) return ! if higher order is worse
                                                     ! by significant factor
                                                     ! `SAFE', then quit early.
   enddo
   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 :: coords(3,n)
  real(c_double), target :: forces_kim(3,n)
  real(c_double) :: forces(3,n)
  integer(c_int) middleDum
  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)

  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")
  endif

  ! 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")
    endif
    call kim_is_routine_present(model_handle, model_routine_name, present, &
      required, ierr)
    if (ierr /= 0) then
      call my_error("kim_is_routine_present")
    endif

    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.")
      endif
    endif
  enddo

  ! 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")
  endif

  call check_model_compatibility(compute_arguments_handle, forces_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")
  endif
  ! 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)
  enddo
  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)
     enddo
  enddo

  ! 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
  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")
  endif

  ! Allocate storage for neighbor lists and
  ! store pointers to neighbor list object and access function
  !
  allocate(neighborlist(n+1,n))
  neighobject%neighbor_list_pointer = c_loc(neighborlist)
  neighobject%number_of_particles = 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)
  enddo
  print '("Using FCC lattice parameter: ",f12.5)', fccspacing

  do i=1,n
    call kim_get_species_support_and_code(model_handle, &
      cluster_species(i), species_is_supported, particlespeciescodes(i), ierr)
  enddo
  if (ierr /= 0) then
    call my_error("kim_api_get_species_code")
  endif
  do i=1,n
    particlecontributing(i) = 1  ! every particle contributes
  enddo
  do i=1,n
     coords(:,i) = cluster_coords(:,i) + cluster_disps(:,i)
  enddo

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

  ! 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")
  endif

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

  ! Print results to screen
  !
  print '(41(''=''))'
  print '("Energy = ",ES25.15)', energy
  print '(41(''=''))'
  print *

  ! 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")
    endif
  endif

  ! Compute gradient using numerical differentiation
  !
  allocate(forces_num(dim,n),forces_num_err(dim,n))
  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")
        endif
        forces_num(j,i) = -deriv
        forces_num_err(j,i) = deriv_err
     enddo
  enddo

  ! Continue printing results to screen
  !
  print '(A6,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)',"Part","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 = "     "
        endif
        weight = max(abs(forces_num(j,i)),eps_prec)/ &
                 max(abs(forces_num_err(j,i)),eps_prec)
        term = weight*forcediff**2
        if (term.gt.term_max) then
           term_max = term
           imax = i
           jmax = j
        endif
        forcediff_sumsq = forcediff_sumsq + term
        weight_sum = weight_sum + weight
        if (j.eq.1) then
           print '(I6,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
                  i,particlespeciescodes(i),j,forces(j,i),forces_num(j,i), &
                  forcediff,forces_num_err(j,i),weight,passfail
        else
           print '(14X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
                  j,forces(j,i),forces_num(j,i), &
                  forcediff,forces_num_err(j,i),weight,passfail
        endif
     enddo
     print *
  enddo
  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))

  ! Free temporary storage
  !
  deallocate(forces_num)
  deallocate(forces_num_err)
  deallocate(neighborlist)
  deallocate(coordsave)
  deallocate(cutoffs, 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")
  endif
  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