!****************************************************************************
!
! 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) 2012, Regents of the University of Minnesota.  All rights reserved.
!
!****************************************************************************
!
! MODULE Exp6_KongChakrabarty_1973_ArNe
!
!
! Compute energy and forces on an isolated cluster of atoms using
! the exp-6 potential shifted to have zero energy at cutoff.
! Parameters due to Hogervorst [1] and mixing rule between species
! due to Kong and Chakrabarty [2].
!
! Author : E. B. Tadmor (10-MAR-12), updated to kim api v1.6 4-AUG-14,
!          updated to kim api v2.0 30-JUN-18 by D. S. Karls.
!
! Support atom species:
!
!  1 = "Ar"
!  2 = "Ne"
!
!  References:
!
!  [1] W. Hogervorst, Physica, vol 51, 77-89 (1971). (Like interactions.)
!
!  [2] C. L. Kong and M. R. Chakrabarty, J. Phys. Chem., Vol. 77, 2668-2670
!      (1973). (Mixing rule.)
!
!****************************************************************************

module Exp6_KongChakrabarty_1973_ArNe

use, intrinsic :: iso_c_binding
use kim_model_headers_module
implicit none

save
private
public Compute_Energy_Forces, &
       compute_arguments_create, &
       compute_arguments_destroy, &
       destroy, &
       speccodeAr, &
       speccodeNe, &
       cutoff, &
       buffer_type

integer(c_int), parameter :: cd = c_double ! used for literal constants

!
! Define global potential parameters
!
real(c_double), parameter :: kb     = 8.6173e-5_cd ! Boltzmann's constant [eV/K]
real(c_double), parameter :: cutoff = 8.15_cd      ! cutoff radius [A] (arbitrary value)
real(c_double), parameter :: eps11  = 138.0_cd*kb  ! Ar-Ar epsilon parameter [eV]
real(c_double), parameter :: rm11   = 3.77_cd      ! Ar-Ar rm parameter [A]
real(c_double), parameter :: alf11  = 14.8_cd      ! Ar-Ar alpha parameter [-]
real(c_double), parameter :: eps22  = 43.0_cd*kb   ! Ne-Ne epsilon parameter [eV]
real(c_double), parameter :: rm22   = 3.03_cd      ! Ne-Ne rm parameter [A]
real(c_double), parameter :: alf22  = 16.0_cd      ! Ne-Ne alpha parameter [-]
real(c_double), parameter :: eps12  = 68.89_cd*kb  ! Ar-Ne epsilon parameter [eV]
real(c_double), parameter :: rm12   = 3.447_cd     ! Ar-Ne rm parameter [A]
real(c_double), parameter :: alf12  = 15.52_cd     ! Ar-Ne alpha parameter [-]
integer(c_int), parameter :: speccodeAr = 1, speccodeNe = 2

type, bind(c) :: buffer_type
  real(c_double) :: influence_distance
  real(c_double) :: cutoff(1)
  real(c_double) :: cutsq(1)
  integer(c_int) :: &
    model_will_not_request_neighbors_of_noncontributing_particles(1)
end type buffer_type

contains

!-------------------------------------------------------------------------------
!
!  Calculate exp-6 pair potential phi(r) and its derivative dphi(r)
!  Note: potential is shifted to have zero energy at r=cutoff
!
!-------------------------------------------------------------------------------
! The "recursive" keyword is added below make this routine thread safe
recursive subroutine calc_phi_dphi(r,phi,dphi,eps,rm,alf,cutoff)
  implicit none

  !-- Transferred variables
  real(c_double), intent(in)  :: r
  real(c_double), intent(in)  :: eps,rm,alf,cutoff
  real(c_double), intent(out) :: phi,dphi

  !-- Local variables
  real(c_double) soa,amp,ror,cor,phicut

  if (r .gt. cutoff) then
    ! Argument exceeds cutoff radius
    phi = 0.0_cd
    dphi = 0.0_cd
  else
    soa = 6.0_cd/alf
    amp = eps/(1.0_cd-soa)
    cor = cutoff/rm
    ror = r/rm
    phicut = amp*(soa*exp(alf*(1.0_cd-cor))-cor**(-6))
    phi = amp*(soa*exp(alf*(1.0_cd-ror))-ror**(-6)) - phicut
    dphi = -6.0_cd*amp/rm*(exp(alf*(1.0_cd-ror))-ror**(-7))
  endif

  return

end subroutine calc_phi_dphi

!-------------------------------------------------------------------------------
!
! Compute energy and forces on atoms from the positions.
!
!-------------------------------------------------------------------------------
! The "recursive" keyword is added below make this routine thread safe
recursive subroutine Compute_Energy_Forces(model_compute_handle, &
  model_compute_arguments_handle, ierr) bind(c)
  implicit none

  !-- Transferred variables
  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
  type(kim_model_compute_arguments_handle_type), intent(in) :: &
    model_compute_arguments_handle
  integer(c_int), intent(out) :: ierr

  !--Local Variables
  integer(c_int), parameter :: DIM=3
  integer(c_int) :: i, j, jj, numnei, ierr2
  integer(c_int) :: comp_energy, comp_force
  real(c_double) Rij(DIM), Rsqij, r, phi, dphi

  !-- KIM variables
  integer(c_int), pointer :: N
  real(c_double), pointer :: energy
  real(c_double), pointer :: coor(:,:)
  real(c_double), pointer :: force(:,:)
  integer(c_int), pointer :: nei1part(:)
  integer(c_int), pointer :: particleSpeciesCodes(:)
  integer(c_int), pointer :: particleContributing(:)
  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf

  ! Avoid unused argument warning
  if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue

  ! Get buffer address from KIM API
  call kim_get_model_buffer_pointer(model_compute_handle, pbuf)
  call c_f_pointer(pbuf, buf)

  ierr = 0
  ierr2 = 0

  ! Unpack data from KIM object
  !
  ierr = 0
  call kim_get_argument_pointer( &
    model_compute_arguments_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_NUMBER_OF_PARTICLES, N, ierr2)
  ierr = ierr + ierr2
  call kim_get_argument_pointer( &
    model_compute_arguments_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_SPECIES_CODES, N, particleSpeciesCodes, &
    ierr2)
  ierr = ierr + ierr2
  call kim_get_argument_pointer( &
    model_compute_arguments_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_CONTRIBUTING, N, particleContributing, &
    ierr2)
  ierr = ierr + ierr2
  call kim_get_argument_pointer( &
    model_compute_arguments_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_COORDINATES, dim, N, coor, ierr2)
  ierr = ierr + ierr2
  call kim_get_argument_pointer( &
    model_compute_arguments_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, energy, ierr2)
  ierr = ierr + ierr2
  call kim_get_argument_pointer( &
    model_compute_arguments_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, dim, N, force, ierr2)
  ierr = ierr + ierr2
  if (ierr /= 0) then
    call kim_log_entry(model_compute_arguments_handle, &
      KIM_LOG_VERBOSITY_ERROR, &
      "Failed to retrieve data from KIM API compute-arguments object")
    return
  endif

  ! Check to see if we have been asked to compute the energy and/or forces
  !
  if (associated(energy)) then
    comp_energy = 1
  else
    comp_energy = 0
  end if
  if (associated(force)) then
    comp_force = 1
  else
    comp_force = 0
  end if

  ! Check to be sure that the species are correct
  !
  ierr = 1 ! assume an error
  do i = 1,N
    if (particleSpeciesCodes(i).ne.speccodeAr .and. particleSpeciesCodes(i).ne.speccodeNe) then
      call kim_log_entry(model_compute_arguments_handle, &
        KIM_LOG_VERBOSITY_ERROR, "Unexpected species code detected")
      return
    endif
  enddo
  ierr = 0 ! everything is ok

  ! Compute energy and forces
  if (comp_energy.eq.1) energy = 0.0_cd
  if (comp_force.eq.1)  force  = 0.0_cd

  do i=1,N
    if (particleContributing(i) == 1) then
      ! Get neighbor list of current atom
      call kim_get_neighbor_list( &
         model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
      if (ierr /= 0) then
         ! some sort of problem, exit
         call kim_log_entry(model_compute_arguments_handle, &
           KIM_LOG_VERBOSITY_ERROR, "GetNeighborList failed")
         ierr = 1
         return
      endif

      do jj = 1, numnei

        j = nei1part(jj) ! Get index of neighbor

        if ( .not.( (particleContributing(j) == 1) .and. &
             j.lt.i) ) then ! effective half list
          Rij(:) = coor(:,j) - coor(:,i)
          Rsqij = dot_product(Rij,Rij)

          if ( Rsqij.lt.buf%cutsq(1) ) then
             r = sqrt(Rsqij)
             if (particleSpeciesCodes(i).eq.1 .and. particleSpeciesCodes(j).eq.1) then
                call calc_phi_dphi(r,phi,dphi,eps11,rm11,alf11,buf%cutoff(1))
             elseif (particleSpeciesCodes(i).eq.2 .and. particleSpeciesCodes(j).eq.2) then
                call calc_phi_dphi(r,phi,dphi,eps22,rm22,alf22,buf%cutoff(1))
             else
                call calc_phi_dphi(r,phi,dphi,eps12,rm12,alf12,buf%cutoff(1))
             endif
             if (comp_energy.eq.1) then
               if (particleContributing(j) == 1) then
                 energy = energy + phi
               else
                 energy = energy + 0.5_cd*phi
               endif
             endif

             if (comp_force.eq.1) then
               if (particleContributing(j).eq.1) then
                 force(:,i) = force(:,i) + dphi*Rij/r
                 force(:,j) = force(:,j) - dphi*Rij/r
               else
                 force(:,i) = force(:,i) + 0.5*dphi*Rij/r
                 force(:,j) = force(:,j) - 0.5*dphi*Rij/r
               endif
             endif

          endif
        endif
     enddo ! End loop over neighbors of atom i
    endif ! Check on whether particle is contributing
enddo ! End primary loop over particles

! No errors
ierr = 0

return
end subroutine Compute_Energy_Forces

!-------------------------------------------------------------------------------
!
! Model destroy routine
!
!-------------------------------------------------------------------------------
! The "recursive" keyword is added below make this routine thread safe
recursive subroutine destroy(model_destroy_handle, ierr) bind(c)
  use, intrinsic :: iso_c_binding
  implicit none

  !-- Transferred variables
  type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle
  integer(c_int), intent(out) :: ierr
  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf

  call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
  call c_f_pointer(pbuf, buf)

  call kim_log_entry(model_destroy_handle, KIM_LOG_VERBOSITY_INFORMATION, &
    "deallocating model buffer")
  deallocate(buf)
  ierr = 0  ! everything is good
end subroutine destroy

!-------------------------------------------------------------------------------
!
! Model compute arguments create routine (REQUIRED)
!
!-------------------------------------------------------------------------------
! The "recursive" keyword is added below make this routine thread safe
recursive subroutine compute_arguments_create(model_compute_handle, &
  model_compute_arguments_create_handle, ierr) bind(c)
  use, intrinsic :: iso_c_binding
  use kim_model_compute_arguments_create_module
  implicit none

  !-- Transferred variables
  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
  type(kim_model_compute_arguments_create_handle_type), intent(inout) :: &
    model_compute_arguments_create_handle
  integer(c_int), intent(out) :: ierr

  integer(c_int) :: ierr2

  ! Avoid unused argument warning
  if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue

  ierr = 0
  ierr2 = 0

  ! register arguments
  call kim_set_argument_support_status( &
    model_compute_arguments_create_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, &
    KIM_SUPPORT_STATUS_OPTIONAL, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_support_status( &
    model_compute_arguments_create_handle, &
    KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, &
    KIM_SUPPORT_STATUS_OPTIONAL, ierr2)
  ierr = ierr + ierr2

  ! register call backs (ProcessDEDrTerm, ProcessD2EDr2Term
  ! NONE

  if (ierr /= 0) then
    ierr = 1
    call kim_log_entry(model_compute_arguments_create_handle, &
      KIM_LOG_VERBOSITY_ERROR, &
      "Unable to successfully create compute-arguments object")
  endif

  return
end subroutine compute_arguments_create

!-------------------------------------------------------------------------------
!
! Model compute arguments destroy routine (REQUIRED)
!
!-------------------------------------------------------------------------------
! The "recursive" keyword is added below make this routine thread safe
recursive subroutine compute_arguments_destroy(model_compute_handle, &
  model_compute_arguments_destroy_handle, ierr) bind(c)
  use, intrinsic :: iso_c_binding
  implicit none

  !-- Transferred variables
  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
  type(kim_model_compute_arguments_destroy_handle_type), intent(inout) :: &
    model_compute_arguments_destroy_handle
  integer(c_int), intent(out) :: ierr

  ! avoid unused dummy argument warnings
  if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue
  if (model_compute_arguments_destroy_handle .eq. &
    KIM_MODEL_COMPUTE_ARGUMENTS_DESTROY_NULL_HANDLE) continue

  ierr = 0
  return
end subroutine compute_arguments_destroy

end module Exp6_KongChakrabarty_1973_ArNe

!-------------------------------------------------------------------------------
!
! Model initialization routine (REQUIRED)
!
!-------------------------------------------------------------------------------
! The "recursive" keyword is added below make this routine thread safe
recursive subroutine create(model_create_handle, requested_length_unit, &
  requested_energy_unit, requested_charge_unit, requested_temperature_unit, &
  requested_time_unit, ierr) bind(c)
  use, intrinsic :: iso_c_binding
  use Exp6_KongChakrabarty_1973_ArNe
  use kim_model_headers_module
  implicit none

  !-- Transferred variables
  type(kim_model_create_handle_type), intent(inout) :: model_create_handle
  type(kim_length_unit_type), intent(in) :: requested_length_unit
  type(kim_energy_unit_type), intent(in) :: requested_energy_unit
  type(kim_charge_unit_type), intent(in) :: requested_charge_unit
  type(kim_temperature_unit_type), intent(in) :: requested_temperature_unit
  type(kim_time_unit_type), intent(in) :: requested_time_unit
  integer(c_int), intent(out) :: ierr

  !-- KIM variables
  integer(c_int) :: ierr2
  type(buffer_type), pointer :: buf

  ierr = 0
  ierr2 = 0

  call kim_set_units(model_create_handle, &
    KIM_LENGTH_UNIT_A, &
    KIM_ENERGY_UNIT_EV, &
    KIM_CHARGE_UNIT_UNUSED, &
    KIM_TEMPERATURE_UNIT_UNUSED, &
    KIM_TIME_UNIT_UNUSED, &
    ierr2)
  ierr = ierr + ierr2

  ! register species
  call kim_set_species_code(model_create_handle, &
    KIM_SPECIES_NAME_AR, speccodeAr, ierr2)
  ierr = ierr + ierr2
  call kim_set_species_code(model_create_handle, &
    KIM_SPECIES_NAME_NE, speccodeNe, ierr2)
  ierr = ierr + ierr2

  ! register numbering
  call kim_set_model_numbering(model_create_handle, &
    KIM_NUMBERING_ONE_BASED, ierr2);
  ierr = ierr + ierr2

  ! register function pointers
  call kim_set_routine_pointer(model_create_handle, &
    KIM_MODEL_ROUTINE_NAME_COMPUTE, KIM_LANGUAGE_NAME_FORTRAN, &
    1, c_funloc(Compute_Energy_Forces), ierr2)
  ierr = ierr + ierr2
  call kim_set_routine_pointer( &
    model_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_CREATE, &
    KIM_LANGUAGE_NAME_FORTRAN, 1, c_funloc(compute_arguments_create), ierr2)
  ierr = ierr + ierr2
  call kim_set_routine_pointer( &
    model_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_DESTROY, &
    KIM_LANGUAGE_NAME_FORTRAN, 1, c_funloc(compute_arguments_destroy), ierr2)
  ierr = ierr + ierr2
  call kim_set_routine_pointer(model_create_handle, &
    KIM_MODEL_ROUTINE_NAME_DESTROY, KIM_LANGUAGE_NAME_FORTRAN, &
    1, c_funloc(destroy), ierr2)
  ierr = ierr + ierr2

  ! allocate buffer
  allocate(buf)

  ! store model buffer in KIM object
  call kim_set_model_buffer_pointer(model_create_handle, &
    c_loc(buf))

  ! set buffer values
  buf%influence_distance = cutoff
  buf%cutoff(1) = cutoff
  buf%cutsq(1) = cutoff*cutoff
  buf%model_will_not_request_neighbors_of_noncontributing_particles(1) = 1

  ! register influence distance
  call kim_set_influence_distance_pointer( &
    model_create_handle, buf%influence_distance)

  ! register cutoff
  call kim_set_neighbor_list_pointers(model_create_handle, &
    1, buf%cutoff, &
    buf%model_will_not_request_neighbors_of_noncontributing_particles)

  if (ierr /= 0) then
    ierr = 1
    call kim_log_entry(model_create_handle, KIM_LOG_VERBOSITY_ERROR, &
      "Unable to successfully initialize model")
  endif
end subroutine create