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

!****************************************************************************
!**
!**  MODULE ex_model_driver_P_LJ
!**
!**  Lennard-Jones pair potential KIM Model Driver
!**  shifted to have zero energy at the cutoff radius
!**
!**  Language: Fortran 2003
!**
!****************************************************************************

module ex_model_driver_p_lj

  use, intrinsic :: iso_c_binding
  use kim_model_driver_headers_module
  implicit none

  save
  private
  public &
    buffer_type, &
    compute_energy_forces, &
    compute_arguments_create, &
    compute_arguments_destroy, &
    refresh, &
    write_model, &
    destroy, &
    calc_phi, &
    calc_phi_dphi, &
    calc_phi_dphi_d2phi, &
    speccode

  ! Below are the definitions and values of all Model parameters
  integer(c_int), parameter          :: cd = c_double  ! for literal constants
  integer(c_int), parameter          :: dim = 3        ! dimensionality of space
  integer(c_int), parameter          :: speccode = 1   ! internal species code

  !-----------------------------------------------------------------------------
  !
  !  Definition of Buffer type
  !
  !-----------------------------------------------------------------------------
  type, bind(c) :: buffer_type
    character(c_char) :: species_name(100)
    real(c_double) :: influence_distance(1)
    real(c_double) :: pcutoff(1)
    real(c_double) :: cutsq(1)
    integer(c_int) :: &
      model_will_not_request_neighbors_of_noncontributing_particles(1)
    real(c_double) :: epsilon(1)
    real(c_double) :: sigma(1)
    real(c_double) :: shift(1)
  end type buffer_type

contains

  !-----------------------------------------------------------------------------
  !
  !  Calculate pair potential phi(r)
  !
  !-----------------------------------------------------------------------------
  recursive subroutine calc_phi(model_epsilon, &
                                model_sigma, &
                                model_shift, &
                                model_cutoff, r, phi)
    implicit none

    !-- Transferred variables
    real(c_double), intent(in)  :: model_epsilon
    real(c_double), intent(in)  :: model_sigma
    real(c_double), intent(in)  :: model_shift
    real(c_double), intent(in)  :: model_cutoff
    real(c_double), intent(in)  :: r
    real(c_double), intent(out) :: phi

    !-- Local variables
    real(c_double) rsq, sor, sor6, sor12

    rsq = r * r             !  r^2
    sor = model_sigma / r   !  (sig/r)
    sor6 = sor * sor * sor     !
    sor6 = sor6 * sor6       !  (sig/r)^6
    sor12 = sor6 * sor6       !  (sig/r)^12
    if (r > model_cutoff) then
      ! Argument exceeds cutoff radius
      phi = 0.0_cd
    else
      phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
    end if

  end subroutine calc_phi

  !-----------------------------------------------------------------------------
  !
  !  Calculate pair potential phi(r) and its derivative dphi(r)
  !
  !-----------------------------------------------------------------------------
  recursive subroutine calc_phi_dphi(model_epsilon, &
                                     model_sigma, &
                                     model_shift, &
                                     model_cutoff, r, phi, dphi)
    implicit none

    !-- Transferred variables
    real(c_double), intent(in)  :: model_epsilon
    real(c_double), intent(in)  :: model_sigma
    real(c_double), intent(in)  :: model_shift
    real(c_double), intent(in)  :: model_cutoff
    real(c_double), intent(in)  :: r
    real(c_double), intent(out) :: phi, dphi

    !-- Local variables
    real(c_double) rsq, sor, sor6, sor12

    rsq = r * r             !  r^2
    sor = model_sigma / r   !  (sig/r)
    sor6 = sor * sor * sor     !
    sor6 = sor6 * sor6       !  (sig/r)^6
    sor12 = sor6 * sor6       !  (sig/r)^12
    if (r > model_cutoff) then
      ! Argument exceeds cutoff radius
      phi = 0.0_cd
      dphi = 0.0_cd
    else
      phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
      dphi = 24.0_cd * model_epsilon * (-2.0_cd * sor12 + sor6) / r
    end if

  end subroutine calc_phi_dphi

  !-----------------------------------------------------------------------------
  !
  !  Calculate pair potential phi(r) and its derivatives dphi(r) and d2phi(r)
  !
  !-----------------------------------------------------------------------------
  recursive subroutine calc_phi_dphi_d2phi(model_epsilon, &
                                           model_sigma, &
                                           model_shift, &
                                           model_cutoff, r, phi, dphi, d2phi)
    implicit none

    !-- Transferred variables
    real(c_double), intent(in)  :: model_epsilon
    real(c_double), intent(in)  :: model_sigma
    real(c_double), intent(in)  :: model_shift
    real(c_double), intent(in)  :: model_cutoff
    real(c_double), intent(in)  :: r
    real(c_double), intent(out) :: phi, dphi, d2phi

    !-- Local variables
    real(c_double) rsq, sor, sor6, sor12

    rsq = r * r             !  r^2
    sor = model_sigma / r   !  (sig/r)
    sor6 = sor * sor * sor     !
    sor6 = sor6 * sor6       !  (sig/r)^6
    sor12 = sor6 * sor6       !  (sig/r)^12
    if (r > model_cutoff) then
      ! Argument exceeds cutoff radius
      phi = 0.0_cd
      dphi = 0.0_cd
      d2phi = 0.0_cd
    else
      phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
      dphi = 24.0_cd * model_epsilon * (-2.0_cd * sor12 + sor6) / r
      d2phi = 24.0_cd * model_epsilon * (26.0_cd * sor12 - 7.0_cd * sor6) / rsq
    end if

  end subroutine calc_phi_dphi_d2phi

  !-----------------------------------------------------------------------------
  !
  ! Compute energy and forces on particles from the positions.
  !
  !-----------------------------------------------------------------------------
  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
    real(c_double) :: r, Rsqij, phi, dphi, d2phi, dEidr, d2Eidr
    integer(c_int) :: i, j, jj, numnei
    integer(c_int) :: ierr2
    integer(c_int) :: comp_force, comp_energy, comp_enepot, comp_process_dEdr, &
                      comp_process_d2Edr2
    type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf

    real(c_double), target :: Rij(dim)
    real(c_double), target :: Rij_pairs(dim, 2)
    real(c_double), target :: r_pairs(2)
    integer(c_int), target :: i_pairs(2), j_pairs(2)

    !-- KIM variables
    real(c_double) :: model_cutoff
    integer(c_int), pointer :: N
    real(c_double), pointer :: energy
    real(c_double), pointer :: coor(:, :)
    real(c_double), pointer :: force(:, :)
    real(c_double), pointer :: enepot(:)
    integer(c_int), pointer :: nei1part(:)
    integer(c_int), pointer :: particleSpeciesCodes(:)
    integer(c_int), pointer :: particleContributing(:)

    ! get model buffer from KIM object
    call kim_get_model_buffer_pointer(model_compute_handle, pbuf)
    call c_f_pointer(pbuf, buf)

    model_cutoff = buf%influence_distance(1)

    ! Check to see if we have been asked to compute the forces, energyperpart,
    ! energy and d1Edr
    !
    ierr = 0
    call kim_is_callback_present( &
      model_compute_arguments_handle, &
      kim_compute_callback_name_process_dedr_term, comp_process_dedr, ierr2)
    ierr = ierr + ierr2
    call kim_is_callback_present( &
      model_compute_arguments_handle, &
      kim_compute_callback_name_process_d2edr2_term, comp_process_d2edr2, ierr2)
    ierr = ierr + ierr2
    if (ierr /= 0) then
      call kim_log_entry(model_compute_arguments_handle, &
                         kim_log_verbosity_error, "get_compute")
      ierr = 1
      return
    end if

    ! 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
    call kim_get_argument_pointer( &
      model_compute_arguments_handle, &
      kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
    ierr = ierr + ierr2
    if (ierr /= 0) then
      call kim_log_entry(model_compute_arguments_handle, &
                         kim_log_verbosity_error, "get_argument_pointer")
      ierr = 1
      return
    end if

    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
    if (associated(enepot)) then
      comp_enepot = 1
    else
      comp_enepot = 0
    end if

    ! Check to be sure that the species are correct
    !

    ierr = 1 ! assume an error
    do i = 1, n
      if (particlespeciescodes(i) /= speccode) then
        call kim_log_entry( &
          model_compute_arguments_handle, kim_log_verbosity_error, &
          "Unexpected species code detected")
        ierr = 1
        return
      end if
    end do
    ierr = 0 ! everything is ok

    ! Initialize potential energies, forces
    !
    if (comp_enepot == 1) enepot = 0.0_cd
    if (comp_energy == 1) energy = 0.0_cd
    if (comp_force == 1) force = 0.0_cd

    !
    !  Compute energy and forces
    !

    !  Loop over particles and compute energy and forces
    !
    do i = 1, n

      if (particlecontributing(i) == 1) then
        ! Set up neighbor list for next particle
        !
        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, &
            "kim_api_get_neigh")
          ierr = 1
          return
        end if

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

          j = nei1part(jj)                        ! get neighbor ID

          ! compute relative position vector
          !
          rij(:) = coor(:, j) - coor(:, i)       ! distance vector between i j

          ! compute energy and forces
          !
          rsqij = dot_product(rij, rij)                ! compute square distance
          if (rsqij < buf%cutsq(1)) then         ! particles are interacting?

            r = sqrt(rsqij)                          ! compute distance
            if (comp_process_d2edr2 == 1) then
              call calc_phi_dphi_d2phi(buf%epsilon(1), &
                                       buf%sigma(1), &
                                       buf%shift(1), &
                                       buf%Pcutoff(1), &
                                       r, phi, dphi, d2phi) ! compute pair potential
              !   and it derivatives
              deidr = 0.5_cd * dphi               !      regular contribution
              d2eidr = 0.5_cd * d2phi
            elseif (comp_force == 1 .or. comp_process_dedr == 1) then
              call calc_phi_dphi(buf%epsilon(1), &
                                 buf%sigma(1), &
                                 buf%shift(1), &
                                 buf%Pcutoff(1), &
                                 r, phi, dphi)        ! compute pair potential
              !   and it derivative

              deidr = 0.5_cd * dphi                !      regular contribution
            else
              call calc_phi(buf%epsilon(1), &
                            buf%sigma(1), &
                            buf%shift(1), &
                            buf%Pcutoff(1), &
                            r, phi)                  ! compute just pair potential
            end if

            ! contribution to energy
            !
            if (comp_enepot == 1) then
              enepot(i) = enepot(i) + 0.5_cd * phi    ! accumulate energy
            end if
            if (comp_energy == 1) then
              energy = energy + 0.5_cd * phi       !      add half v to total energy
            end if

            ! contribution to process_dEdr
            !
            if (comp_process_dedr == 1) then
              call kim_process_dedr_term( &
                model_compute_arguments_handle, deidr, r, rij, i, j, ierr)
            end if

            ! contribution to process_d2Edr2
            if (comp_process_d2edr2 == 1) then
              r_pairs(1) = r
              r_pairs(2) = r
              rij_pairs(:, 1) = rij
              rij_pairs(:, 2) = rij
              i_pairs(1) = i
              i_pairs(2) = i
              j_pairs(1) = j
              j_pairs(2) = j

              call kim_process_d2edr2_term( &
                model_compute_arguments_handle, d2eidr, &
                r_pairs, rij_pairs, i_pairs, j_pairs, ierr)
            end if

            ! contribution to forces
            !
            if (comp_force == 1) then
              force(:, i) = force(:, i) + deidr * rij / r ! accumulate force on particle i
              force(:, j) = force(:, j) - deidr * rij / r ! accumulate force on particle j
            end if

          end if

        end do  ! loop on jj

      end if  ! if particleContributing

    end do  ! do i

    ! Everything is great
    !
    ierr = 0
    return

  end subroutine compute_energy_forces

  !-----------------------------------------------------------------------------
  !
  ! Model driver refresh routine
  !
  !-----------------------------------------------------------------------------
  recursive subroutine refresh(model_refresh_handle, ierr) bind(c)
    implicit none

    !-- transferred variables
    type(kim_model_refresh_handle_type), intent(in) :: model_refresh_handle
    integer(c_int), intent(out) :: ierr

    !-- Local variables
    real(c_double) energy_at_cutoff
    type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf

    ! get model buffer from KIM object
    call kim_get_model_buffer_pointer(model_refresh_handle, pbuf)
    call c_f_pointer(pbuf, buf)

    call kim_set_influence_distance_pointer(model_refresh_handle, &
                                            buf%influence_distance(1))
    call kim_set_neighbor_list_pointers( &
      model_refresh_handle, 1, buf%influence_distance, &
      buf%model_will_not_request_neighbors_of_noncontributing_particles)

    ! Set new values in KIM object and buffer
    !
    buf%influence_distance(1) = buf%Pcutoff(1)
    buf%cutsq(1) = (buf%Pcutoff(1))**2
    ! calculate pair potential at r=cutoff with shift=0.0
    call calc_phi(buf%epsilon(1), &
                  buf%sigma(1), &
                  0.0_cd, &
                  buf%Pcutoff(1), &
                  buf%Pcutoff(1), energy_at_cutoff)
    buf%shift(1) = -energy_at_cutoff

    ierr = 0
    return

  end subroutine refresh

  !-----------------------------------------------------------------------------
  !
  ! Model driver write_model routine
  !
  !-----------------------------------------------------------------------------
  recursive subroutine write_model( &
    model_write_parameterized_model_handle, ierr) bind(c)
    implicit none

    !-- transferred variables
    type(kim_model_write_parameterized_model_handle_type), intent(in) &
      :: model_write_parameterized_model_handle
    integer(c_int), intent(out) :: ierr

    !-- Local variables
    integer i
    type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
    character(len=512, kind=c_char) :: path
    character(len=512, kind=c_char) :: model_name
    character(len=512, kind=c_char) :: string_buffer
    character(len=100, kind=c_char) :: species_name

    ! get model buffer from KIM object
    call kim_get_model_buffer_pointer( &
      model_write_parameterized_model_handle, pbuf)
    call c_f_pointer(pbuf, buf)

    call kim_get_path(model_write_parameterized_model_handle, path)
    call kim_get_model_name(model_write_parameterized_model_handle, model_name)

    write (string_buffer, '(A)') trim(model_name)//".params"
    call kim_set_parameter_file_name(model_write_parameterized_model_handle, &
                                     string_buffer)
    write (string_buffer, '(A)') trim(path)//"/"//trim(string_buffer)

    open (42, file=trim(string_buffer), &
          status="REPLACE", action="WRITE", iostat=ierr)
    if (ierr /= 0) then
      call kim_log_entry( &
        model_write_parameterized_model_handle, kim_log_verbosity_error, &
        "Unable to open parameter file for writing.")
      return
    end if

    do i = 1, 100
      species_name(i:i) = buf%species_name(i)
    end do
    write (42, '(A)') trim(species_name)
    write (42, '(ES20.10)') buf%Pcutoff(1)
    write (42, '(ES20.10)') buf%epsilon(1)
    write (42, '(ES20.10)') buf%sigma(1)

    ierr = 0
    return

  end subroutine write_model

  !-----------------------------------------------------------------------------
  !
  ! Model driver destroy routine
  !
  !-----------------------------------------------------------------------------
  recursive subroutine destroy(model_destroy_handle, ierr) bind(c)
    implicit none

    !-- Transferred variables
    type(kim_model_destroy_handle_type), intent(in) :: model_destroy_handle
    integer(c_int), intent(out) :: ierr

    !-- Local variables
    type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf

    ! get model buffer from KIM object
    call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
    call c_f_pointer(pbuf, buf)

    deallocate (buf)

    ierr = 0
    return

  end subroutine destroy

  !-----------------------------------------------------------------------------
  !
  ! Model driver compute arguments create routine
  !
  !-----------------------------------------------------------------------------
  recursive subroutine compute_arguments_create( &
    model_compute_handle, model_compute_arguments_create_handle, ierr) bind(c)
    implicit none

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

    integer(c_int) ierr2

    ! avoid unsed dummy argument warnings
    if (model_compute_handle == 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, ierr)
    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
    call kim_set_argument_support_status( &
      model_compute_arguments_create_handle, &
      kim_compute_argument_name_partial_particle_energy, &
      kim_support_status_optional, ierr2)
    ierr = ierr + ierr2
    if (ierr /= 0) then
      call kim_log_entry( &
        model_compute_arguments_create_handle, kim_log_verbosity_error, &
        "Unable to register arguments support_statuss")
      ierr = 1
      goto 42
    end if

    ! register callbacks
    call kim_set_callback_support_status( &
      model_compute_arguments_create_handle, &
      kim_compute_callback_name_process_dedr_term, &
      kim_support_status_optional, ierr)
    call kim_set_callback_support_status( &
      model_compute_arguments_create_handle, &
      kim_compute_callback_name_process_d2edr2_term, &
      kim_support_status_optional, ierr2)
    ierr = ierr + ierr2
    if (ierr /= 0) then
      call kim_log_entry( &
        model_compute_arguments_create_handle, kim_log_verbosity_error, &
        "Unable to register callbacks support_statuss")
      ierr = 1
      goto 42
    end if

    ierr = 0
42  continue
    return

  end subroutine compute_arguments_create

  !-----------------------------------------------------------------------------
  !
  ! Model driver compute arguments destroy routine
  !
  !-----------------------------------------------------------------------------
  recursive subroutine compute_arguments_destroy( &
    model_compute_handle, model_compute_arguments_destroy_handle, ierr) bind(c)
    implicit none

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

    ! avoid unsed dummy argument warnings
    if (model_compute_handle == kim_model_compute_null_handle) continue
    if (model_compute_arguments_destroy_handle == &
        kim_model_compute_arguments_destroy_null_handle) continue

    ierr = 0
    return
  end subroutine compute_arguments_destroy

end module ex_model_driver_p_lj

!-----------------------------------------------------------------------------
!
! Model driver create routine (REQUIRED)
!
!-----------------------------------------------------------------------------
recursive subroutine model_driver_create_routine( &
  model_driver_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 ex_model_driver_p_lj
  use kim_model_driver_headers_module
  implicit none
  integer(c_int), parameter :: cd = c_double ! used for literal constants

  !-- Transferred variables
  type(kim_model_driver_create_handle_type), intent(in) &
    :: model_driver_create_handle
  type(kim_length_unit_type), intent(in), value :: requested_length_unit
  type(kim_energy_unit_type), intent(in), value :: requested_energy_unit
  type(kim_charge_unit_type), intent(in), value :: requested_charge_unit
  type(kim_temperature_unit_type), intent(in), value :: &
    requested_temperature_unit
  type(kim_time_unit_type), intent(in), value :: requested_time_unit
  integer(c_int), intent(out) :: ierr

  !-- Local variables
  integer i
  integer(c_int) :: number_of_parameter_files
  character(len=1024, kind=c_char) :: parameter_file_directory_name
  character(len=1024, kind=c_char) :: parameter_file_basename
  character(len=1024, kind=c_char) :: parameter_file_name
  integer(c_int) :: ierr2
  type(buffer_type), pointer :: buf
  type(kim_species_name_type) species_name
  ! define variables for all model parameters to be read in
  real(c_double) factor
  character(len=100, kind=c_char) in_species
  real(c_double) in_cutoff
  real(c_double) in_epsilon
  real(c_double) in_sigma
  real(c_double) energy_at_cutoff

  ! register units
  call kim_set_units( &
    model_driver_create_handle, &
    requested_length_unit, &
    requested_energy_unit, &
    kim_charge_unit_unused, &
    kim_temperature_unit_unused, &
    kim_time_unit_unused, ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "Unable to set units")
    ierr = 1
    goto 42
  end if

  ! register numbering
  call kim_set_model_numbering( &
    model_driver_create_handle, kim_numbering_one_based, ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "Unable to set numbering")
    ierr = 1
    goto 42
  end if

  ! store callback pointers in KIM object
  call kim_set_routine_pointer( &
    model_driver_create_handle, kim_model_routine_name_compute, &
    kim_language_name_fortran, 1, c_funloc(compute_energy_forces), ierr)
  call kim_set_routine_pointer( &
    model_driver_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_driver_create_handle, kim_model_routine_name_refresh, &
    kim_language_name_fortran, 1, c_funloc(refresh), ierr2)
  ierr = ierr + ierr2
  call kim_set_routine_pointer( &
    model_driver_create_handle, &
    kim_model_routine_name_write_parameterized_model, &
    kim_language_name_fortran, 0, c_funloc(write_model), ierr2)
  ierr = ierr + ierr2
  call kim_set_routine_pointer( &
    model_driver_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_driver_create_handle, kim_model_routine_name_destroy, &
    kim_language_name_fortran, 1, c_funloc(destroy), ierr2)
  ierr = ierr + ierr2
  if (ierr /= 0) then
    call kim_log_entry( &
      model_driver_create_handle, kim_log_verbosity_error, &
      "Unable to store callback pointers")
    ierr = 1
    goto 42
  end if

  ! process parameter files
  call kim_get_number_of_parameter_files( &
    model_driver_create_handle, number_of_parameter_files)
  if (number_of_parameter_files /= 1) then
    call kim_log_entry( &
      model_driver_create_handle, kim_log_verbosity_error, &
      "Wrong number of parameter files")
    ierr = 1
    goto 42
  end if

  ! Read in model parameters from parameter file
  !
  call kim_get_parameter_file_directory_name( &
    model_driver_create_handle, parameter_file_directory_name)
  call kim_get_parameter_file_basename( &
    model_driver_create_handle, 1, parameter_file_basename, ierr)
  if (ierr /= 0) then
    call kim_log_entry( &
      model_driver_create_handle, kim_log_verbosity_error, &
      "Unable to get parameter file basename")
    ierr = 1
    goto 42
  end if
  parameter_file_name = trim(parameter_file_directory_name)//"/"// &
                        trim(parameter_file_basename)
  open (10, file=parameter_file_name, status="old")
  read (10, *, iostat=ierr, err=100) in_species
  read (10, *, iostat=ierr, err=100) in_cutoff
  read (10, *, iostat=ierr, err=100) in_epsilon
  read (10, *, iostat=ierr, err=100) in_sigma
  close (10)

  goto 200
100 continue
  ! reading parameters failed
  call kim_log_entry(model_driver_create_handle, &
                     kim_log_verbosity_error, "Unable to read LJ parameters")
  ierr = 1
  goto 42

200 continue

  ! register species
  call kim_from_string(in_species, species_name)
  call kim_set_species_code( &
    model_driver_create_handle, species_name, speccode, ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "Unable to set species code")
    ierr = 1
    goto 42
  end if

  ! convert to appropriate units
  call kim_convert_unit( &
    kim_length_unit_a, &
    kim_energy_unit_ev, &
    kim_charge_unit_e, &
    kim_temperature_unit_k, &
    kim_time_unit_ps, &
    requested_length_unit, &
    requested_energy_unit, &
    requested_charge_unit, &
    requested_temperature_unit, &
    requested_time_unit, &
    1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "kim_api_convert_to_act_unit")
    ierr = 1
    goto 42
  end if
  in_cutoff = in_cutoff * factor

  call kim_convert_unit( &
    kim_length_unit_a, &
    kim_energy_unit_ev, &
    kim_charge_unit_e, &
    kim_temperature_unit_k, &
    kim_time_unit_ps, &
    requested_length_unit, &
    requested_energy_unit, &
    requested_charge_unit, &
    requested_temperature_unit, &
    requested_time_unit, &
    0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "kim_api_convert_to_act_unit")
    ierr = 1
    goto 42
  end if
  in_epsilon = in_epsilon * factor

  call kim_convert_unit( &
    kim_length_unit_a, &
    kim_energy_unit_ev, &
    kim_charge_unit_e, &
    kim_temperature_unit_k, &
    kim_time_unit_ps, &
    requested_length_unit, &
    requested_energy_unit, &
    requested_charge_unit, &
    requested_temperature_unit, &
    requested_time_unit, &
    1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "kim_api_convert_to_act_unit")
    ierr = 1
    goto 42
  end if
  in_sigma = in_sigma * factor

  allocate (buf)

  ! setup buffer
  ! set value of parameters
  do i = 1, 100
    buf%species_name(i) = in_species(i:i)
  end do
  buf%influence_distance(1) = in_cutoff
  buf%Pcutoff(1) = in_cutoff
  buf%cutsq(1) = in_cutoff**2
  buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
  buf%epsilon(1) = in_epsilon
  buf%sigma(1) = in_sigma
  call calc_phi(in_epsilon, &
                in_sigma, &
                0.0_cd, &
                in_cutoff, &
                in_cutoff, energy_at_cutoff)
  buf%shift(1) = -energy_at_cutoff

  ! store model cutoff in KIM object
  call kim_set_influence_distance_pointer( &
    model_driver_create_handle, buf%influence_distance(1))
  call kim_set_neighbor_list_pointers( &
    model_driver_create_handle, 1, buf%influence_distance, &
    buf%model_will_not_request_neighbors_of_noncontributing_particles)

  ! end setup buffer

  ! store in model buffer
  call kim_set_model_buffer_pointer( &
    model_driver_create_handle, c_loc(buf))

  ! set pointers to parameters in KIM object
  call kim_set_parameter_pointer( &
    model_driver_create_handle, buf%pcutoff, "cutoff", &
    "Distance beyond which particles do not interact with one another.", ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "set_parameter")
    ierr = 1
    goto 42
  end if

  call kim_set_parameter_pointer( &
    model_driver_create_handle, buf%epsilon, "epsilon", &
    "Maximum depth of the potential well.", ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "set_parameter")
    ierr = 1
    goto 42
  end if

  call kim_set_parameter_pointer( &
    model_driver_create_handle, buf%sigma, "sigma", &
    "Distance at which energy is zero and force is repulsive.", ierr)
  if (ierr /= 0) then
    call kim_log_entry(model_driver_create_handle, &
                       kim_log_verbosity_error, "set_parameter")
    ierr = 1
    goto 42
  end if

  ierr = 0
42 continue
  return

end subroutine model_driver_create_routine