ex_test_Ar_fcc_cluster_fortran.f90

Go to the documentation of this file.

!
! CDDL HEADER START
!
! The contents of this file are subject to the terms of the Common Development
! and Distribution License Version 1.0 (the "License").
!
! You can obtain a copy of the license at
! http://www.opensource.org/licenses/CDDL-1.0.  See the License for the
! specific language governing permissions and limitations under the License.
!
! When distributing Covered Code, include this CDDL HEADER in each file and
! include the License file in a prominent location with the name LICENSE.CDDL.
! If applicable, add the following below this CDDL HEADER, with the fields
! enclosed by brackets "[]" replaced with your own identifying information:
!
! Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
!
! CDDL HEADER END
!

!
! Copyright (c) 2013--2019, Regents of the University of Minnesota.
! All rights reserved.
!
! Contributors:
!    Ellad B. Tadmor
!    Ryan S. Elliott
!    Stephen M. Whalen
!

module error
  use, intrinsic :: iso_c_binding
  implicit none

  public

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

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

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

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


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

module mod_neighborlist

  use, intrinsic :: iso_c_binding

  public get_neigh

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

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

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

  !-- Local variables
  integer(c_int) numberOfParticles
  type(neighobject_type), pointer :: neighObject
  integer(c_int), pointer :: neighborList(:,:)

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

  if (number_of_neighbor_lists /= 1) then
    call my_warning("invalid number of cutoffs")
    ierr = 1
    return
  endif

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

  if (neighbor_list_index /= 1) then
    call my_warning("wrong list index")
    ierr = 1
    return
  endif

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

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

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

  ierr = 0
  return
end subroutine get_neigh

end module mod_neighborlist


module mod_utility
  implicit none
  public

contains

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

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

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

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

  neighobject%cutoff = cutoff

  cutoff2 = cutoff**2

  do i=1,numberofparticles
     a = 1
     do j=1,numberofparticles
        dx(:) = coords(:, j) - coords(:, i)
        r2 = dot_product(dx, dx)
        if (r2.le.cutoff2) then
           ! part j is a neighbor of part i
           if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) ) then
               a = a+1
               neighborlist(a,i) = j
           endif
        endif
     enddo
     ! part i has a-1 neighbors
     neighborlist(1,i) = a-1
  enddo

  return

end subroutine neigh_pure_cluster_neighborlist

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

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

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

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

  return

end subroutine create_fcc_configuration

end module mod_utility

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

!-------------------------------------------------------------------------------
!
! Main program
!
!-------------------------------------------------------------------------------
program ex_test_ar_fcc_cluster_fortran
  use, intrinsic :: iso_c_binding
  use error
  use kim_simulator_headers_module
  use kim_supported_extensions_module
  use mod_neighborlist
  use mod_utility
  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
  real(c_double), parameter :: FCCspacing     = 5.260_cd
  real(c_double), parameter :: min_spacing    = 0.8*fccspacing
  real(c_double), parameter :: max_spacing    = 1.2*fccspacing
  real(c_double), parameter :: spacing_incr   = 0.025*fccspacing
  real(c_double) :: current_spacing
  real(c_double) :: force_norm

  character(len=256, kind=c_char) :: modelname

  integer(c_int), parameter :: &
    N = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1

  type(neighobject_type), target :: neighObject
  integer(c_int), allocatable, target :: neighborList(:,:)

  type(kim_model_handle_type) :: model_handle
  type(kim_compute_arguments_handle_type) :: compute_arguments_handle
  real(c_double) :: influence_distance
  integer(c_int) :: number_of_neighbor_lists
  real(c_double) :: cutoff
  real(c_double) :: cutoffs(1)
  integer(c_int) :: &
    model_will_not_request_neighbors_of_noncontributing_particles(1)
  integer(c_int), target :: particle_species_codes(n)
  integer(c_int), target :: particle_contributing(n)
  real(c_double), target :: energy
  real(c_double), target :: coords(dim, n)
  real(c_double), target :: forces(dim, n)
  integer(c_int) i, j, ierr, ierr2

  integer(c_int) species_is_supported
  integer(c_int) species_code
  integer(c_int) requested_units_accepted
  integer(c_int) number_of_model_routine_names
  type(kim_model_routine_name_type) model_routine_name
  integer(c_int) present
  integer(c_int) required
  type(kim_supported_extensions_type), target :: supported_extensions
  character(len=KIM_MAX_EXTENSION_ID_LENGTH, kind=c_char) :: id_string

  integer :: middledum

  ! Initialize
  ierr = 0
  supported_extensions%number_of_supported_extensions = 0

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

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

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

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

    if ((present == 1) .and. (required == 1)) then
      if (.not. ((model_routine_name == kim_model_routine_name_create) &
        .or. (model_routine_name == &
        kim_model_routine_name_compute_arguments_create) &
        .or. (model_routine_name == kim_model_routine_name_compute) &
        .or. (model_routine_name == kim_model_routine_name_refresh) &
        .or. (model_routine_name == &
        kim_model_routine_name_compute_arguments_destroy) &
        .or. (model_routine_name == kim_model_routine_name_destroy))) then

        call my_error("Unknown required ModelRoutineName found.")
      endif
    endif
  enddo

  ! check that model supports Ar
  !
  call kim_get_species_support_and_code(model_handle, &
    kim_species_name_ar, species_is_supported, species_code, ierr)
  if ((ierr /= 0) .or. (species_is_supported /= 1)) then
    call my_error("Model does not support Ar")
  endif

  ! Check supported extensions, if any
  call kim_is_routine_present(model_handle, &
    kim_model_routine_name_extension, present, required, ierr)
  if (ierr /= 0) then
    call my_error("Unable to get Extension present/required.")
  endif
  if (present /= 0) then
    call kim_extension(model_handle, kim_supported_extensions_id, &
      c_loc(supported_extensions), ierr)
    if (ierr /= 0) then
      call my_error("Error returned from kim_model_extension().")
    endif
    write (*,'(A,I2,A)') "Model Supports ", &
      supported_extensions%number_of_supported_extensions, &
      " Extensions:"
    do i = 1, supported_extensions%number_of_supported_extensions
      call kim_c_char_array_to_string(&
        supported_extensions%supported_extension_id(:,i), id_string)
      write (*,'(A,I2,A,A,A,A,I2)') " supportedExtensionID[", i, '] = "', &
        trim(id_string), '" ', &
        "which has required = ", &
        supported_extensions%supported_extension_required(i)
    end do
  endif

  ! Best-practice is to check that the model is compatible
  ! but we will skip it here

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

  ! register memory with the KIM system
  ierr = 0
  call kim_set_argument_pointer(compute_arguments_handle, &
    kim_compute_argument_name_number_of_particles, n, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer(compute_arguments_handle, &
    kim_compute_argument_name_particle_species_codes, particle_species_codes, &
    ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer(compute_arguments_handle, &
    kim_compute_argument_name_particle_contributing, particle_contributing, &
    ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer(compute_arguments_handle, &
    kim_compute_argument_name_coordinates, coords, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer(compute_arguments_handle, &
    kim_compute_argument_name_partial_energy, energy, ierr2)
  ierr = ierr + ierr2
  call kim_set_argument_pointer(compute_arguments_handle, &
    kim_compute_argument_name_partial_forces, forces, ierr2)
  ierr = ierr + ierr2
  if (ierr /= 0) then
     call my_error("set_argument_pointer")
  endif

  ! Allocate storage for neighbor lists
  !
  allocate(neighborlist(n+1,n))
  neighobject%neighbor_list_pointer = c_loc(neighborlist)
  neighobject%number_of_particles = n

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

  call kim_get_influence_distance(model_handle, influence_distance)
  call kim_get_number_of_neighbor_lists(model_handle, &
    number_of_neighbor_lists)
  if (number_of_neighbor_lists /= 1) then
    call my_error("too many neighbor lists")
  endif
  call kim_get_neighbor_list_values(model_handle, cutoffs, &
    model_will_not_request_neighbors_of_noncontributing_particles, ierr)
  if (ierr /= 0) then
    call my_error("get_neighbor_list_values")
  end if
  cutoff = cutoffs(1)

  ! Setup cluster
  !
  do i=1,n
    particle_species_codes(i) = species_code
  enddo

  ! setup contributing particles
  do i=1,n
    particle_contributing(i) = 1  ! every particle contributes
  enddo

  ! Print output header
  !
  print *
  print *,'This is Test : ex_test_Ar_fcc_cluster_fortran'
  print *
  print '(80(''-''))'
  print '("Results for KIM Model : ",A)', trim(modelname)

  ! print header
  print '(3A20)', "Energy", "Force Norm", "Lattice Spacing"
  ! do the computations
  current_spacing = min_spacing
  do while (current_spacing < max_spacing)

    call create_fcc_configuration(current_spacing, ncellsperside, .false., &
      coords, middledum)
    ! Compute neighbor lists
    call neigh_pure_cluster_neighborlist(.false., n, coords, (cutoff+cutpad), &
                                         neighobject)

    ! Call model compute
    call kim_compute(model_handle, compute_arguments_handle, ierr)
    if (ierr /= 0) then
      call my_error("kim_api_model_compute")
    endif

    ! compue force_norm
    force_norm = 0.0;
    do i = 1, n
      do j = 1, dim
        force_norm = force_norm + forces(j,i)*forces(j,i)
      end do
    end do
    force_norm = sqrt(force_norm)

    ! Print results to screen
    !
    print '(3ES20.10)', energy, force_norm, current_spacing

    current_spacing = current_spacing + spacing_incr
  enddo

  ! Deallocate neighbor list object
  deallocate( neighborlist )

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

end program ex_test_ar_fcc_cluster_fortran