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

!
! Release: This file is part of the kim-api.git package.
!


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

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

module mod_neighborlist

  use, intrinsic :: iso_c_binding
  use KIM_API_F03

  public setup_neighborlist_KIM_access, get_neigh

  type neighObject_type
     integer(c_int), pointer :: neighborList(:,:)
     real(c_double), pointer :: RijList(:,:,:)
  end type neighObject_type

contains

!-------------------------------------------------------------------------------
!
! setup_neighborlist_KIM_access :
!
!    Store necessary pointers in KIM API object to access the neighbor list
!    data and methods.
!
!-------------------------------------------------------------------------------
subroutine setup_neighborlist_KIM_access(pkim, neighObject)
  implicit none

  !-- Transferred variables
  type(c_ptr),                    intent(in) :: pkim
  type(neighObject_type), target, intent(in) :: neighObject

  !-- Local variables
  integer(c_int), parameter :: SizeOne = 1
  integer(c_int)            :: ier, idum

  ! store location of neighObject variable
  !
  ier = kim_api_set_data(pkim, "neighObject", SizeOne, c_loc(neighObject))
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_set_data", ier)
     stop
  endif
  ier = kim_api_set_method(pkim, "get_neigh", SizeOne, c_funloc(get_neigh))
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_set_method", ier)
     stop
  endif

  return

end subroutine setup_neighborlist_KIM_access

!-------------------------------------------------------------------------------
!
! get_neigh neighbor list access function
!
! This function implements Locator and Iterator mode
!
!-------------------------------------------------------------------------------
integer(c_int) function get_neigh(pkim,mode,request,part,numnei,pnei1part, &
                                  pRij) bind(c)
  implicit none

  !-- Transferred variables
  type(c_ptr),    intent(in)  :: pkim
  integer(c_int), intent(in)  :: mode
  integer(c_int), intent(in)  :: request
  integer(c_int), intent(out) :: part
  integer(c_int), intent(out) :: numnei
  type(c_ptr),    intent(out) :: pnei1part
  type(c_ptr),    intent(out) :: pRij

  !-- Local variables
  integer(c_int), parameter :: DIM = 3
  integer(c_int), save :: iterVal = 0
  integer(c_int)  N
  integer(c_int)  partToReturn
  integer(c_int), pointer :: numberOfParticles;   type(c_ptr) :: pnParts
  type(neighObject_type), pointer :: neighObject; type(c_ptr) :: pneighObject
  integer(c_int)  ier, idum

  ! unpack number of particles
  pnParts = kim_api_get_data(pkim, "numberOfParticles", ier)
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_get_data", ier)
     stop
  endif
  call c_f_pointer(pnParts, numberOfParticles)

  ! unpack neighbor list object
  pneighObject = kim_api_get_data(pkim, "neighObject", ier)
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "kim_api_get_data", ier)
     stop
  endif
  call c_f_pointer(pneighObject, neighObject)

  N = size(neighObject%neighborList, 2)

  ! check mode and request
  if (mode.eq.0) then ! iterator mode
     if (request.eq.0) then ! reset iterator
        iterVal = 0
        get_neigh = KIM_STATUS_NEIGH_ITER_INIT_OK
        return
     elseif (request.eq.1) then ! increment iterator
        iterVal = iterVal + 1
        if (iterVal.gt.N) then
           get_neigh = KIM_STATUS_NEIGH_ITER_PAST_END
           return
        else
           partToReturn = iterVal
        endif
     else
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "Invalid request in get_neigh", &
                                    KIM_STATUS_NEIGH_INVALID_REQUEST)
        get_neigh = KIM_STATUS_NEIGH_INVALID_REQUEST
        return
     endif
  elseif (mode.eq.1) then ! locator mode
     if ( (request.gt.N) .or. (request.lt.1)) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "Invalid part ID in get_neigh", &
                                    KIM_STATUS_PARTICLE_INVALID_ID)
        get_neigh = KIM_STATUS_PARTICLE_INVALID_ID
        return
     else
        partToReturn = request
     endif
  else ! not iterator or locator mode
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "Invalid mode in get_neigh", &
                                 KIM_STATUS_NEIGH_INVALID_MODE)
     get_neigh = KIM_STATUS_NEIGH_INVALID_MODE
     return
  endif

  ! set the returned part
  part = partToReturn

  ! set the returned number of neighbors for the returned part
  numnei = neighObject%neighborList(1,part)

  ! set the location for the returned neighbor list
  pnei1part = c_loc(neighObject%neighborList(2,part))

  ! set pointer to Rij to appropriate value
  if (associated(neighObject%RijList)) then
    pRij = c_loc(neighObject%RijList(1,1,part))
  else
    pRij = c_null_ptr
  endif

  get_neigh = KIM_STATUS_OK
  return
end function get_neigh

end module mod_neighborlist

!*******************************************************************************
!**
!**  PROGRAM utility_config_ener_forces
!**
!**  KIM compliant program to compute the energy and forces for an
!**  isolated cluster configuration read from file.
!**
!**  Works with the following NBC methods:
!**        NEIGH_RVEC_H
!**        NEIGH_PURE_H
!**        NEIGH_RVEC_F
!**        NEIGH_PURE_F
!**        MI_OPBC_H
!**        MI_OPBC_F
!**        CLUSTER
!**
!*******************************************************************************

!-------------------------------------------------------------------------------
!
! Main program
!
!-------------------------------------------------------------------------------
program utility_config_ener_forces
  use, intrinsic :: iso_c_binding
  use KIM_API_F03
  use mod_neighborlist
  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  ! max species allowed
  integer(c_int), parameter :: max_NBCs       = 20  ! maximum number of NBCs
  real(c_double), parameter :: eps_prec       = epsilon(1.0_cd)
  integer(c_int)            :: in
  integer(c_int)            :: N
  real(c_double)            :: max_force_component
  real(c_double)            :: scaled_eps_prec
  integer(c_int)            :: SizeOne = 1
  integer(c_int)            :: num_species_in_config
  integer(c_int)            :: num_NBCs
  logical                   :: found
  integer(c_int)            :: i,j
  real(c_double)            :: force_err(DIM),ave_force_error
  character(len=KIM_KEY_STRING_LENGTH) :: species_in_config(max_species)
  character(len=KIM_KEY_STRING_LENGTH) :: model_NBCs(max_NBCs)

  !
  ! neighbor list
  !
  type(neighObject_type), target :: neighObject
  real(c_double), allocatable    :: coordsave(:,:)
  logical do_update_list

  !
  ! KIM variables
  !
  character(len=KIM_KEY_STRING_LENGTH) :: testname     = "utility_config_ener_forces"
  character(len=KIM_KEY_STRING_LENGTH) :: modelname
  character(len=KIM_KEY_STRING_LENGTH) :: configfile
  ! configuration species (element symbols)
  character(len=KIM_KEY_STRING_LENGTH), pointer :: conf_species(:)
  real(c_double), pointer :: conf_coors(:,:)  ! configuration coordinates
  real(c_double), pointer :: conf_forces(:,:) ! configuration forces
  real(c_double)          :: conf_energy      ! configuration energy

  character(len=KIM_KEY_STRING_LENGTH) :: NBC_Method
  ! 0- NEIGH_RVEC_H, 1- NEIGH_PURE_H, 2- NEIGH_RVEC_F, 3- NEIGH_PURE_F,
  ! 4- MI_OPBC_H,    5- MI_OPBC_F,    6- CLUSTER
  integer(c_int) nbc
  type(c_ptr)             :: pkim
  integer(c_int)          :: ier, idum, inbc
  character(len=10000)    :: test_descriptor_string
  integer(c_int), pointer :: numberOfParticles;    type(c_ptr) pnParts
  integer(c_int), pointer :: numContrib;           type(c_ptr) pnumContrib
  integer(c_int), pointer :: numberOfSpecies;      type(c_ptr) pnOfSpecies
  integer(c_int), pointer :: particleSpecies(:);   type(c_ptr) pparticleSpecies
  real(c_double), pointer :: cutoff;               type(c_ptr) pcutoff
  real(c_double), pointer :: energy;               type(c_ptr) penergy
  real(c_double), pointer :: coords(:,:);          type(c_ptr) pcoor
  real(c_double), pointer :: forces(:,:);          type(c_ptr) pforces
  real(c_double), pointer :: boxSideLengths(:);    type(c_ptr) pboxSideLengths

  ! Initialize error flag
  ier = KIM_STATUS_OK

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

  ! Get filename of file containing configuration
  print '("Please enter the configuration file name: ")'
  read(*,*) configfile

  ! Read in configuration file
  in = 20
  open(unit=in,file=trim(configfile))
  read(in,*,err=10) N
  if (N<1) then
     ier = KIM_STATUS_FAIL
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
           "Error: Number of particles is less than 1", ier)
     stop
  endif
  ! dynamically allocate memory
  allocate(conf_species(N),conf_coors(DIM,N),conf_forces(DIM,N))
!@--- IN FUTURE WILL BE CHANGED TO READ IN PERIODIC BOX INFO
!@--- read(in,*,err=20) conf_boxsize(:)
  do i=1,N
     read(in,*,err=30) conf_species(i), conf_coors(:,i)
  enddo
  read(in,*,err=40) conf_energy
  max_force_component = 0.0_cd
  do i=1,N
     read(in,*,err=50) conf_forces(:,i)
     max_force_component = max(max_force_component, &
                               abs(conf_forces(1,i)), &
                               abs(conf_forces(2,i)), &
                               abs(conf_forces(3,i)))
  enddo
  scaled_eps_prec = max_force_component*eps_prec
  goto 100

  ! Error handing on input
  !

10  continue
    ier = KIM_STATUS_FAIL
    idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
           "Error reading in the number of particles N in the configuration", ier)
    stop

!20 continue
   ier = KIM_STATUS_FAIL
   idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
          "Error reading in the periodic box size in the configuration", ier)
   stop

30 continue
   ier = KIM_STATUS_FAIL
   idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
          "Error reading in a data line (species + coors) in configuration", &
          ier)
   stop

40 continue
   ier = KIM_STATUS_FAIL
   idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
          "Error reading in the configuration energy", ier)
   stop

50 continue
   ier = KIM_STATUS_FAIL
   idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
          "Error reading in forces in configuration", ier)
   stop



  ! No errors -- continue from here

100 continue

  ! Generate list of species appearing in config
  !
  num_species_in_config = 1
  species_in_config(1) = conf_species(1)
  do i=2,N
     found = .false.
     do j=1,num_species_in_config
        found = (trim(conf_species(i))==trim(species_in_config(j)))
        if (found) exit
     enddo
     if (.not.found) then
        num_species_in_config = num_species_in_config + 1
        species_in_config(num_species_in_config) = conf_species(i)
     endif
  enddo

  ! Get list of NBCs supported by the model
  !
  call Get_Model_NBC_methods(modelname, max_NBCs, model_NBCs, num_NBCs, ier)
  if (ier.lt.KIM_STATUS_OK) then
     idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                 "Get_Model_NBC_methods", ier)
     stop
  endif

  ! Print output header
  !
  print *
  print *,'VERIFICATION CHECK: COMPUTE ENERGY AND FORCES FOR CONFIGURATION'
  print *
  print '(120(''-''))'
  print '("This is Test          : ",A)', trim(testname)
  print '("Results for KIM Model : ",A)', trim(modelname)
  print *
  print *,'*** Read-in Configuration ***'
  print *
  print '(a,f20.10)', "Energy = ",conf_energy
  print *
  print '(A6,2X,A4,2X,A)',"Part","Spec","Coordinates"
  do i=1,N
     print '(I6,2X,A4,2X,3ES25.15)',i,conf_species(i), conf_coors(:,i)
  enddo
  print *
  print '(A6,2X,A4,2X,A)',"Part","Spec","Force"
  do i=1,N
     print '(I6,2X,A4,2X,3ES25.15)',i,conf_species(i),conf_forces(:,i)
  enddo
  print *

  ! Loop over all NBCs and compute energy and forces for each one
  !
  do inbc = 1, num_NBCs

     ! Write out KIM descriptor string for Test for current NBC
     !
     call Write_KIM_descriptor(model_NBCs(inbc), max_species, species_in_config, &
                               num_species_in_config, test_descriptor_string, ier)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "Write_KIM_descriptor", ier)
        stop
     endif

     ! Create empty KIM object conforming to fields in the KIM descriptor files
     ! of the Test and Model
     !
     ier = kim_api_string_init(pkim,trim(test_descriptor_string)//char(0), &
                               modelname)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_string_init", ier)
        stop
     endif

     ! Double check that the NBC method being used is what we think it is
     !
     ier = kim_api_get_nbc_method(pkim, NBC_Method)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_get_nbc_method", ier)
        stop
     endif
     if (index(NBC_Method,trim(model_NBCs(inbc))).ne.1) then
        ier = KIM_STATUS_FAIL
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
         "Internal Error: Selected NBC method different from requested value", &
         ier)
        stop
     endif

     ! Set NBC code based on selected NBC method
     !
     if (index(NBC_Method,"NEIGH_RVEC_H").eq.1) then
        nbc = 0
     elseif (index(NBC_Method,"NEIGH_PURE_H").eq.1) then
        nbc = 1
     elseif (index(NBC_Method,"NEIGH_RVEC_F").eq.1) then
        nbc = 2
     elseif (index(NBC_Method,"NEIGH_PURE_F").eq.1) then
        nbc = 3
     elseif (index(NBC_Method,"MI_OPBC_H").eq.1) then
        nbc = 4
     elseif (index(NBC_Method,"MI_OPBC_F").eq.1) then
        nbc = 5
     elseif (index(NBC_Method,"CLUSTER").eq.1) then
        nbc = 6
     else
        ier = KIM_STATUS_FAIL
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "Unknown NBC method", ier)
        stop
     endif

     ! Allocate memory via the KIM system
     !
     call kim_api_allocate(pkim, N, num_species_in_config, ier)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_allocate", ier)
        stop
     endif

     ! Allocate storage for neighbor lists and
     ! store pointers to neighbor list object and access function
     !
     if (nbc.le.5) then
        allocate(neighObject%neighborList(N+1,N))
        if (nbc.eq.0 .or. nbc.eq.2) then
          allocate(neighObject%RijList(DIM,N+1,N))
        endif
        ier = kim_api_set_data(pkim, "neighObject", SizeOne, &
                               c_loc(neighObject))
        if (ier.lt.KIM_STATUS_OK) then
          idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                      "kim_api_set_data", ier)
          stop
        endif
     endif

     ! Set pointer in KIM object to neighbor list routine
     !
     if (nbc.ne.6) then
        ier = kim_api_set_method(pkim, "get_neigh", SizeOne, &
                                 c_funloc(get_neigh))
        if (ier.lt.KIM_STATUS_OK) then
           idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                       "kim_api_set_method", ier)
           stop
        endif
     endif

     ! Initialize Model
     !
     ier = kim_api_model_init(pkim)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_model_init", ier)
        stop
     endif

     ! Unpack data from KIM object
     !
     call kim_api_getm_data(pkim, ier, &
          "numberOfParticles",           pnParts,           1,                               &
          "numberContributingParticles", pnumContrib,       TRUEFALSE(nbc.eq.0.or.nbc.eq.1.or.nbc.eq.4), &
          "numberOfSpecies",             pnOfSpecies,       1,                               &
          "particleSpecies",             pparticleSpecies,  1,                               &
          "coordinates",                 pcoor,             1,                               &
          "cutoff",                      pcutoff,           1,                               &
          "boxSideLengths",              pboxSideLengths,   TRUEFALSE(nbc.eq.4.or.nbc.eq.5), &
          "energy",                      penergy,           1,                               &
          "forces",                      pforces,           1)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_getm_data", ier)
        stop
     endif
     call c_f_pointer(pnParts, numberOfParticles)
     if ((nbc.eq.0).or.(nbc.eq.1).or.(nbc.eq.4)) call c_f_pointer(pnumContrib, &
                                                                  numContrib)
     call c_f_pointer(pnOfSpecies,      numberOfSpecies)
     call c_f_pointer(pparticleSpecies, particleSpecies, [N])
     call c_f_pointer(pcoor, coords, [DIM,N])
     call c_f_pointer(pcutoff, cutoff)
     if ((nbc.eq.4).or.(nbc.eq.5)) call c_f_pointer(pboxSideLengths, &
                                                    boxSideLengths, [DIM])
     call c_f_pointer(penergy, energy)
     call c_f_pointer(pforces, forces, [DIM,N])

     ! Set values in KIM object
     !
     numberOfParticles   = N
     if (nbc.eq.0.or.nbc.eq.1.or.nbc.eq.4) numContrib = N
     numberOfSpecies = num_species_in_config
     do i=1,N
        particleSpecies(i) = kim_api_get_species_code(pkim,&
                                                      trim(conf_species(i)),ier)
     enddo
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_get_species_code", ier)
        stop
     endif
     do i=1,N
        coords(:,i) = conf_coors(:,i)
     enddo
     ! set boxSideLengths large enough to make the cluster isolated
     if (nbc.eq.4.or.nbc.eq.5) boxSideLengths(:) = 600.0_cd

     ! Compute neighbor lists
     !
     if (nbc.le.5) then
        do_update_list = .true.
        allocate(coordsave(DIM,N))
        call update_neighborlist(DIM,N,coords,cutoff,cutpad,boxSideLengths, &
                                 NBC_Method,do_update_list,coordsave, &
                                 neighObject,ier)
        if (ier.lt.KIM_STATUS_OK) then
           idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                       "update_neighborlist", ier)
           stop
        endif
     endif

     ! Call model compute to get energy and forces
     !
     ier = kim_api_model_compute(pkim)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_model_compute", ier)
        stop
     endif

     ! print results to screen
     !
     print '(41(''=''))'
     print '("NBC Method = ",A28)', trim(NBC_Method)
     print '(41(''=''))'
     print *
     print '(a,f20.10)', "Energy = ",energy
     print *
     print '(A6,2X,A4,2X,A)',"Part","Spec","Computed Force"
     do i=1,N
        print '(I6,2X,A4,2X,3ES25.15)',i,conf_species(i),forces(:,i)
     enddo
     print *
     print *,'*** Energy and Forces Agreement ***'
     print *
     print '(A6,2X,A4,2X,A)',"Part","Spec","Force Error"
     ave_force_error = 0.0_cd
     do i=1,N
        do j=1,DIM
           force_err(j) = abs(forces(j,i)-conf_forces(j,i))/ &
                          max(abs(conf_forces(j,i)),scaled_eps_prec)
        enddo
        print '(I6,2X,A4,2X,3ES25.15)',i,conf_species(i),force_err(:)
        ave_force_error = ave_force_error + dot_product(force_err,force_err)
     enddo
     ave_force_error = sqrt(ave_force_error)/dble(DIM*N)
     print *
     print '(a,ES25.15)', "Average force error = ",ave_force_error
     print *
     print '(a,ES25.15)', "Energy error        = ",abs(energy-conf_energy)/ &
                                                   max(abs(conf_energy), &
                                                   eps_prec)
     print *

     ! Free temporary storage
     !
     if (nbc.le.5) then ! deallocate neighbor list storage
        deallocate(neighObject%neighborList)
        deallocate(coordsave)
        if (nbc.eq.0.or.nbc.eq.2) then
           deallocate(neighObject%RijList)
        endif
     endif
     ier = kim_api_model_destroy(pkim)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_model_destroy", ier)
        stop
     endif
     call kim_api_free(pkim, ier)
     if (ier.lt.KIM_STATUS_OK) then
        idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                                    "kim_api_free", ier)
        stop
     endif

  enddo ! loop over NBC methods

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

  ! Free configuration storage
  !
  deallocate(conf_species,conf_coors,conf_forces)
  stop

end program utility_config_ener_forces

!-------------------------------------------------------------------------------
!
!  Get the number and identities of NBCs supported by the Model.
!
!-------------------------------------------------------------------------------
subroutine Get_Model_NBC_methods(modelname, max_NBCs, model_NBCs, num_NBCs, ier)
use, intrinsic :: iso_c_binding
use KIM_API_F03
implicit none

!-- Transferred variables
character(len=KIM_KEY_STRING_LENGTH), intent(in)   :: modelname
integer(c_int),                       intent(in)   :: max_NBCs
character(len=KIM_KEY_STRING_LENGTH), intent(out)  :: model_NBCs(max_NBCs)
integer(c_int),                       intent(out)  :: num_NBCs
integer(c_int),                       intent(out)  :: ier

!-- Local variables
integer(c_int) :: index

!-- KIM variables
type(c_ptr) pkim

! Initialize error flag
ier = KIM_STATUS_OK

! Generate empty KIM object containing generic model info
ier = kim_api_model_info(pkim, modelname)
if (ier.lt.KIM_STATUS_OK) return

! Identify supported NBCs by seeking index to each of the NBCs in the KIM model
num_NBCs = 0

! NEIGH_RVEC_H
index = kim_api_get_index(pkim, "NEIGH_RVEC_H", ier)
if (index.ge.0) then
   num_NBCs = num_NBCs + 1
   if (num_NBCs>max_NBCs) then
      ier = KIM_STATUS_FAIL
      return
   endif
   model_NBCs(num_NBCs) = "NEIGH_RVEC_H"
endif

! NEIGH_PURE_H
index = kim_api_get_index(pkim, "NEIGH_PURE_H", ier)
if (index.ge.0) then
   num_NBCs = num_NBCs + 1
   if (num_NBCs>max_NBCs) then
      ier = KIM_STATUS_FAIL
      return
   endif
   model_NBCs(num_NBCs) = "NEIGH_PURE_H"
endif

! NEIGH_RVEC_F
index = kim_api_get_index(pkim, "NEIGH_RVEC_F", ier)
if (index.ge.0) then
   num_NBCs = num_NBCs + 1
   if (num_NBCs>max_NBCs) then
      ier = KIM_STATUS_FAIL
      return
   endif
   model_NBCs(num_NBCs) = "NEIGH_RVEC_F"
endif

! NEIGH_PURE_F
index = kim_api_get_index(pkim, "NEIGH_PURE_F", ier)
if (index.ge.0) then
   num_NBCs = num_NBCs + 1
   if (num_NBCs>max_NBCs) then
      ier = KIM_STATUS_FAIL
      return
   endif
   model_NBCs(num_NBCs) = "NEIGH_PURE_F"
endif

! MI_OPBC_H
index = kim_api_get_index(pkim, "MI_OPBC_H", ier)
if (index.ge.0) then
   num_NBCs = num_NBCs + 1
   if (num_NBCs>max_NBCs) then
      ier = KIM_STATUS_FAIL
      return
   endif
   model_NBCs(num_NBCs) = "MI_OPBC_H"
endif

! MI_OPBC_F
index = kim_api_get_index(pkim, "MI_OPBC_F", ier)
if (index.ge.0) then
   num_NBCs = num_NBCs + 1
   if (num_NBCs>max_NBCs) then
      ier = KIM_STATUS_FAIL
      return
   endif
   model_NBCs(num_NBCs) = "MI_OPBC_F"
endif

! CLUSTER
index = kim_api_get_index(pkim, "CLUSTER", ier)
if (index.ge.0) then
   num_NBCs = num_NBCs + 1
   if (num_NBCs>max_NBCs) then
      ier = KIM_STATUS_FAIL
      return
   endif
   model_NBCs(num_NBCs) = "CLUSTER"
endif

! free temporary storage
call kim_api_free(pkim,ier)

return

end subroutine Get_Model_NBC_methods

!-------------------------------------------------------------------------------
!
!  Write KIM descriptor file for MiniMol for given NBC and set of
!  supported species
!
!-------------------------------------------------------------------------------
subroutine Write_KIM_descriptor(NBC_method, max_species, model_species, num_species, &
                                kim_descriptor, ier)
use, intrinsic :: iso_c_binding
use KIM_API_F03
implicit none

!-- Transferred variables
character(len=KIM_KEY_STRING_LENGTH), intent(in)   :: NBC_method
integer(c_int),                       intent(in)   :: max_species
character(len=KIM_KEY_STRING_LENGTH), intent(in)   :: model_species(max_species)
integer(c_int),                       intent(in)   :: num_species
character(len=10000),                 intent(out)  :: kim_descriptor
integer(c_int),                       intent(out)  :: ier

!-- Local variables
integer(c_int) :: i
character(len=103)  :: divider
character(len=1)    :: cr
character(len=1024) :: species_line
character(len=32)   :: nbcline

! Initialize error flag
ier = KIM_STATUS_OK

! Define frequently used variables
!
cr = char(10)
divider = '#######################################################################################################'

! Write Minimol descriptor file into string kim_descriptor
!
kim_descriptor = &
   divider                                                                      // cr // &
   '#'                                                                          // cr // &
   '# CDDL HEADER START'                                                        // cr // &
   '#'                                                                          // cr // &
   '# The contents of this file are subject to the terms of the Common Development' // cr // &
   '# and Distribution License Version 1.0 (the "License").'                    // cr // &
   '#'                                                                          // cr // &
   '# You can obtain a copy of the license at'                                  // cr // &
   '# http://www.opensource.org/licenses/CDDL-1.0.  See the License for the'    // cr // &
   '# specific language governing permissions and limitations under the License.' // cr // &
   '#'                                                                          // cr // &
   '# When distributing Covered Code, include this CDDL HEADER in each file and'// cr // &
   '# include the License file in a prominent location with the name LICENSE.CDDL.' // cr // &
   '# If applicable, add the following below this CDDL HEADER, with the fields' // cr // &
   '# enclosed by brackets "[]" replaced with your own identifying information:'// cr // &
   '#'                                                                          // cr // &
   '# Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.' // cr // &
   '#'                                                                          // cr // &
   '# CDDL HEADER END'                                                          // cr // &
   '#'                                                                          // cr // &
                                                                                   cr // &
   '#'                                                                          // cr // &
   '# Copyright (c) 2013--2018, Regents of the University of Minnesota.'        // cr // &
   '# All rights reserved.'                                                     // cr // &
   '#'                                                                          // cr // &
   '# Contributors:'                                                            // cr // &
   '#    Automatically generated by calling Test'                               // cr // &
   '#'                                                                          // cr // &
                                                                                   cr // &
   '#'                                                                          // cr // &
   '# See src/standard.kim for documentation about this file'                   // cr // &
   '#'                                                                          // cr // &
   divider                                                                      // cr // &
                                                                                   cr // &
   'KIM_API_Version := 1.6.0'                                                   // cr // &
                                                                                   cr // &
   'Unit_length      := A'                                                      // cr // &
   'Unit_energy      := eV'                                                     // cr // &
   'Unit_charge      := e'                                                      // cr // &
   'Unit_temperature := K'                                                      // cr // &
   'Unit_time        := ps'                                                     // cr // &
                                                                                   cr // &
                                                                                   cr // &
   divider                                                                      // cr // &
   'PARTICLE_SPECIES:'                                                          // cr // &
   '# Symbol/name               Type                    code'                   // cr

   do i = 1,num_species
      write(species_line,'(a,'' '',''spec'',20x,i4)') trim(model_species(i)),0
      kim_descriptor = trim(kim_descriptor) // trim(species_line) // cr
   enddo

   nbcline = NBC_method(1:28) // 'flag'

   kim_descriptor = trim(kim_descriptor) // &
                                                                                   cr // &
                                                                                   cr // &
   divider                                                                      // cr // &
   'CONVENTIONS:'                                                               // cr // &
   '# Name                      Type'                                           // cr // &
                                                                                   cr // &
   'OneBasedLists               flag'                                           // cr // &
                                                                                   cr // &
   'Neigh_BothAccess            flag'                                           // cr // &
                                                                                   cr // &
   nbcline                                                                      // cr // &
                                                                                   cr // &
   divider                                                                      // cr // &
   'MODEL_INPUT:'                                                               // cr // &
   '# Name                      Type         Unit       Shape              requirements' // cr // &
   'numberOfParticles           integer      none       []'                     // cr // &
                                                                                   cr // &
   'numberOfSpecies             integer      none       []'                     // cr // &
                                                                                   cr // &
   'particleSpecies             integer      none       [numberOfParticles]'    // cr // &
                                                                                   cr // &
   'coordinates                 double       length     [numberOfParticles,3]'  // cr // &
                                                                                   cr // &
   'get_neigh                   method       none       []'                     // cr // &
                                                                                   cr // &
   'neighObject                 pointer      none       []'                     // cr // &
                                                                                   cr // &
   'numberContributingParticles integer      none       []'                     // cr // &
                                                                                   cr // &
   'boxSideLengths              double       length     [3]'                    // cr // &
                                                                                   cr // &
                                                                                   cr // &
   divider                                                                      // cr // &
   'MODEL_OUTPUT:'                                                              // cr // &
   '# Name                      Type         Unit       Shape              requirements' // cr // &
                                                                                   cr // &
   'destroy                     method       none       []'                     // cr // &
                                                                                   cr // &
   'compute                     method       none       []'                     // cr // &
                                                                                   cr // &
   'cutoff                      double       length     []'                     // cr // &
                                                                                   cr // &
   'energy                      double       energy     []'                     // cr // &
                                                                                   cr // &
   'forces                      double       force      [numberOfParticles,3]'  // cr // &
                                                                                   cr // &
   divider                                                                            // cr

return

end subroutine Write_KIM_descriptor

subroutine update_neighborlist(DIM,N,coords,cutoff,cutpad,boxSideLengths, &
                               NBC_Method,do_update_list,coordsave, &
                               neighObject,ier)
use, intrinsic :: iso_c_binding
use KIM_API_F03
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
real(c_double),         intent(in)    :: boxSideLengths(DIM)
character(len=64),      intent(in)    :: NBC_Method
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)   :: ier

!-- Local variables
! 0- NEIGH_RVEC_H, 1- NEIGH_PURE_H, 2- NEIGH_RVEC_F, 3- NEIGH_PURE_F,
! 4- MI_OPBC_H, 5- MI_OPBC_F 
integer(c_int) nbc
real(c_double) disp, disp1, disp2, cutrange, dispvec(DIM)
integer(c_int) i, idum

! Initialize error code
!
ier = KIM_STATUS_OK

! Determine which NBC scenario to use
!
if (index(NBC_Method,"NEIGH_RVEC_H").eq.1) then
   nbc = 0
elseif (index(NBC_Method,"NEIGH_PURE_H").eq.1) then
   nbc = 1
elseif (index(NBC_Method,"NEIGH_RVEC_F").eq.1) then
   nbc = 2
elseif (index(NBC_Method,"NEIGH_PURE_F").eq.1) then
   nbc = 3
elseif (index(NBC_Method,"MI_OPBC_H").eq.1) then
   nbc = 4
elseif (index(NBC_Method,"MI_OPBC_F").eq.1) then
   nbc = 5
else
   ier = KIM_STATUS_FAIL
   idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
                               "Unknown NBC method", ier)
   stop
endif

! Update neighbor lists if necessary, if not just update Rij vectors if these
! are used
!
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
   if (nbc.eq.0) then
      call NEIGH_RVEC_cluster_neighborlist(.true., N, coords, cutrange, &
                                           neighObject)
   elseif (nbc.eq.1) then
      call NEIGH_PURE_cluster_neighborlist(.true., N, coords, cutrange, &
                                           neighObject)
   elseif (nbc.eq.2) then
      call NEIGH_RVEC_cluster_neighborlist(.false., N, coords, cutrange, &
                                           neighObject)
   elseif (nbc.eq.3) then
      call NEIGH_PURE_cluster_neighborlist(.false., N, coords, cutrange, &
                                           neighObject)
   elseif (nbc.eq.4) then
      call MI_OPBC_cluster_neighborlist(.true., N, coords, cutrange, &
                                        boxSideLengths, neighObject)
   elseif (nbc.eq.5) then
      call MI_OPBC_cluster_neighborlist(.false., N, coords, cutrange, &
                                        boxSideLengths, neighObject)
   endif

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

else

   ! Even though neighbor lists are uptodate Rij vectors still need to be
   ! computed
   if (nbc.eq.0.or.nbc.eq.2) &
      call NEIGH_RVEC_update_Rij_vectors(DIM, N, coords, neighObject)

endif

return

end subroutine update_neighborlist



!-------------------------------------------------------------------------------
!
! MI_OPBC_cluster_neighborlist : construct a half or full neighbor list using
!                                the particle coordinates in coords()
!
!-------------------------------------------------------------------------------
subroutine MI_OPBC_cluster_neighborlist(half, numberOfParticles, coords, rcut, &
                                        boxSideLengths, neighObject)
  use, intrinsic :: iso_c_binding
  use KIM_API_F03
  use mod_neighborlist
  implicit none
  integer(c_int), parameter :: cd = c_double ! used for literal constants

  !-- 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)                  :: rcut
  real(c_double), dimension(3), intent(in)    :: boxSideLengths
  type(neighObject_type),       intent(inout) :: neighObject

  !-- Local variables
  integer(c_int) i, j, a
  real(c_double) dx(3)
  real(c_double) r2
  real(c_double) rcut2

  rcut2 = rcut**2

  do i=1,numberOfParticles
     a = 1
     do j=1,numberOfParticles
        dx(:) = coords(:, j) - coords(:, i)
        where (abs(dx) > 0.5_cd*boxSideLengths)  ! apply PBC
           dx = dx - sign(boxSideLengths,dx)
        endwhere
        r2 = dot_product(dx, dx)
        if (r2.le.rcut2) then
           if (i.ne.j) then
              if ( (j .gt. i) .or. ((.not. half) .AND. (i.ne.j)) ) then
                  ! part j is a neighbor of part i
                  a = a+1
                  neighObject%neighborList(a,i) = j
              endif
           endif
        endif
     enddo
     ! part i has a-1 neighbors
     neighObject%neighborList(1,i) = a-1
  enddo

  return

end subroutine MI_OPBC_cluster_neighborlist

!-------------------------------------------------------------------------------
!
! NEIGH_PURE_cluster_neighborlist
!
!-------------------------------------------------------------------------------
subroutine NEIGH_PURE_cluster_neighborlist(half, numberOfParticles, coords, &
                                           cutoff, neighObject)
  use, intrinsic :: iso_c_binding
  use KIM_API_F03
  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) i, j, a
  real(c_double) dx(3)
  real(c_double) r2
  real(c_double) cutoff2

  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
               neighObject%neighborList(a,i) = j
           endif
        endif
     enddo
     ! part i has a-1 neighbors
     neighObject%neighborList(1,i) = a-1
  enddo

  return

end subroutine NEIGH_PURE_cluster_neighborlist

!-------------------------------------------------------------------------------
!
! NEIGH_RVEC_cluster_neighborlist
!
!-------------------------------------------------------------------------------
subroutine NEIGH_RVEC_cluster_neighborlist(half, numberOfParticles, coords, &
                                           cutoff, neighObject)
  use, intrinsic :: iso_c_binding
  use KIM_API_F03
  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) i, j, a
  real(c_double) dx(3)
  real(c_double) r2
  real(c_double) cutoff2

  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
           if ((half .and. i.lt.j) .or. (.not.half .and. i.ne.j)) then
              ! part j is a neighbor of part i
              a = a+1
              neighObject%neighborList(a,i) = j
              neighObject%RijList(:,a-1,i) = dx
           endif
        endif
     enddo
     ! part i has a-1 neighbors
     neighObject%neighborList(1,i) = a-1
  enddo

  return

end subroutine NEIGH_RVEC_cluster_neighborlist

!-------------------------------------------------------------------------------
!
! NEIGH_RVEC_update_Rij_vectors
!
!-------------------------------------------------------------------------------
subroutine NEIGH_RVEC_update_Rij_vectors(DIM, N, coords, neighObject)
  use, intrinsic :: iso_c_binding
  use mod_neighborlist
  implicit none

  !-- Transferred variables
  integer(c_int), intent(in)  :: DIM
  integer(c_int), intent(in)  :: N
  real(c_double), intent(in)  :: coords(DIM,N)
  type(neighObject_type), intent(inout) :: neighObject

  !-- Local variables
  integer(c_int) i, j, jj, NN
  real(c_double) dx(DIM)

  do i=1,N
     NN = neighObject%neighborList(1,i)
     do jj=1,NN
        j = neighObject%neighborList(jj+1,i)
        dx(:) = coords(:, j) - coords(:, i)
        neighObject%RijList(:,jj,i) = dx(:)
     enddo
  enddo

  return

end subroutine NEIGH_RVEC_update_Rij_vectors
