kim-api  2.3.1-git+v2.3.0-git-2-g378406f9.GNU.GNU.
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
ex_test_Ar_fcc_cluster_fortran.f90
Go to the documentation of this file.
1 !
2 ! KIM-API: An API for interatomic models
3 ! Copyright (c) 2013--2022, Regents of the University of Minnesota.
4 ! All rights reserved.
5 !
6 ! Contributors:
7 ! Ellad B. Tadmor
8 ! Ryan S. Elliott
9 ! Stephen M. Whalen
10 !
11 ! SPDX-License-Identifier: LGPL-2.1-or-later
12 !
13 ! This library is free software; you can redistribute it and/or
14 ! modify it under the terms of the GNU Lesser General Public
15 ! License as published by the Free Software Foundation; either
16 ! version 2.1 of the License, or (at your option) any later version.
17 !
18 ! This library is distributed in the hope that it will be useful,
19 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
20 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21 ! Lesser General Public License for more details.
22 !
23 ! You should have received a copy of the GNU Lesser General Public License
24 ! along with this library; if not, write to the Free Software Foundation,
25 ! Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26 !
27 
28 module error
29  use, intrinsic :: iso_c_binding
30  implicit none
31 
32  public
33 
34 contains
35  recursive subroutine my_error(message)
36  implicit none
37  character(len=*, kind=c_char), intent(in) :: message
38 
39  print *, "* Error : ", trim(message)
40  stop 1
41  end subroutine my_error
42 
43  recursive subroutine my_warning(message)
44  implicit none
45  character(len=*, kind=c_char), intent(in) :: message
46 
47  print *, "* Warning : ", trim(message)
48  end subroutine my_warning
49 end module error
50 
51 !-------------------------------------------------------------------------------
52 !
53 ! module mod_neighborlist :
54 !
55 ! Module contains type and routines related to neighbor list calculation
56 !
57 !-------------------------------------------------------------------------------
58 
60 
61  use, intrinsic :: iso_c_binding
62 
63  public get_neigh
64 
65  type, bind(c) :: neighobject_type
66  real(c_double) :: cutoff
67  integer(c_int) :: number_of_particles
68  type(c_ptr) :: neighbor_list_pointer
69  end type neighobject_type
70 contains
71 
72  !-----------------------------------------------------------------------------
73  !
74  ! get_neigh neighbor list access function
75  !
76  ! This function implements Locator and Iterator mode
77  !
78  !-----------------------------------------------------------------------------
79  recursive subroutine get_neigh(data_object, number_of_neighbor_lists, &
80  cutoffs, neighbor_list_index, request, &
81  numnei, pnei1part, ierr) bind(c)
82  use error
83  implicit none
84 
85  !-- Transferred variables
86  type(c_ptr), value, intent(in) :: data_object
87  integer(c_int), value, intent(in) :: number_of_neighbor_lists
88  real(c_double), intent(in) :: cutoffs(*)
89  integer(c_int), value, intent(in) :: neighbor_list_index
90  integer(c_int), value, intent(in) :: request
91  integer(c_int), intent(out) :: numnei
92  type(c_ptr), intent(out) :: pnei1part
93  integer(c_int), intent(out) :: ierr
94 
95  !-- Local variables
96  integer(c_int) numberOfParticles
97  type(neighobject_type), pointer :: neighObject
98  integer(c_int), pointer :: neighborList(:, :)
99 
100  call c_f_pointer(data_object, neighobject)
101  numberofparticles = neighobject%number_of_particles
102  call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
103  [numberofparticles + 1, numberofparticles])
104 
105  if (number_of_neighbor_lists /= 1) then
106  call my_warning("invalid number of cutoffs")
107  ierr = 1
108  return
109  end if
110 
111  if (cutoffs(1) > neighobject%cutoff) then
112  call my_warning("neighbor list cutoff too small for model cutoff")
113  ierr = 1
114  return
115  end if
116 
117  if (neighbor_list_index /= 1) then
118  call my_warning("wrong list index")
119  ierr = 1
120  return
121  end if
122 
123  if ((request > numberofparticles) .or. (request < 1)) then
124  print *, request
125  call my_warning("Invalid part ID in get_neigh")
126  ierr = 1
127  return
128  end if
129 
130  ! set the returned number of neighbors for the returned part
131  numnei = neighborlist(1, request)
132 
133  ! set the location for the returned neighbor list
134  pnei1part = c_loc(neighborlist(2, request))
135 
136  ierr = 0
137  return
138  end subroutine get_neigh
139 
140 end module mod_neighborlist
141 
143  implicit none
144  public
145 
146 contains
147 
148  !-----------------------------------------------------------------------------
149  !
150  ! NEIGH_PURE_cluster_neighborlist
151  !
152  !-----------------------------------------------------------------------------
153  recursive subroutine neigh_pure_cluster_neighborlist( &
154  half, numberOfParticles, coords, cutoff, neighObject)
155  use, intrinsic :: iso_c_binding
156  use mod_neighborlist
157  implicit none
158 
159  !-- Transferred variables
160  logical, intent(in) :: half
161  integer(c_int), intent(in) :: numberOfParticles
162  real(c_double), dimension(3, numberOfParticles), &
163  intent(in) :: coords
164  real(c_double), intent(in) :: cutoff
165  type(neighobject_type), intent(inout) :: neighObject
166 
167  !-- Local variables
168  integer(c_int) i, j, a
169  real(c_double) dx(3)
170  real(c_double) r2
171  real(c_double) cutoff2
172  integer(c_int), pointer :: neighborList(:, :)
173 
174  call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
175  [numberofparticles + 1, numberofparticles])
176 
177  neighobject%cutoff = cutoff
178 
179  cutoff2 = cutoff**2
180 
181  do i = 1, numberofparticles
182  a = 1
183  do j = 1, numberofparticles
184  dx(:) = coords(:, j) - coords(:, i)
185  r2 = dot_product(dx, dx)
186  if (r2 <= cutoff2) then
187  ! part j is a neighbor of part i
188  if ((j > i) .OR. ((.not. half) .AND. (i /= j))) then
189  a = a + 1
190  neighborlist(a, i) = j
191  end if
192  end if
193  end do
194  ! part i has a-1 neighbors
195  neighborlist(1, i) = a - 1
196  end do
197 
198  return
199 
200  end subroutine neigh_pure_cluster_neighborlist
201 
202  !-----------------------------------------------------------------------------
203  !
204  ! create_FCC_configuration subroutine
205  !
206  ! creates a cubic configuration of FCC particles with lattice spacing
207  ! `FCCspacing' and `nCellsPerSide' cells along each direction.
208  !
209  ! With periodic==.true. this will result in a total number of particles equal
210  ! to 4*(nCellsPerSide)**3 + 6*(nCellsPerSide)**2 + 3*(nCellsPerSide) + 1
211  !
212  ! With periodic==.false. this will result in a total number of particles equal
213  ! to 4*(nCellsPerSide)**3
214  !
215  ! Returns the Id of the particle situated in the middle of the configuration
216  ! (this particle will have the most neighbors.)
217  !
218  !-----------------------------------------------------------------------------
219  recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, &
220  periodic, coords, MiddlePartId)
221  use, intrinsic :: iso_c_binding
222  implicit none
223  integer(c_int), parameter :: cd = c_double ! used for literal constants
224 
225  !-- Transferred variables
226  real(c_double), intent(in) :: FCCspacing
227  integer(c_int), intent(in) :: nCellsPerSide
228  logical, intent(in) :: periodic
229  real(c_double), intent(out) :: coords(3, *)
230  integer(c_int), intent(out) :: MiddlePartId
231  !
232  ! cluster setup variables
233  !
234  real(c_double) FCCshifts(3, 4)
235  real(c_double) latVec(3)
236  integer(c_int) a, i, j, k, m
237 
238  ! Create a cubic FCC cluster
239  !
240  fccshifts(1, 1) = 0.0_cd
241  fccshifts(2, 1) = 0.0_cd
242  fccshifts(3, 1) = 0.0_cd
243  fccshifts(1, 2) = 0.5_cd * fccspacing
244  fccshifts(2, 2) = 0.5_cd * fccspacing
245  fccshifts(3, 2) = 0.0_cd
246  fccshifts(1, 3) = 0.5_cd * fccspacing
247  fccshifts(2, 3) = 0.0_cd
248  fccshifts(3, 3) = 0.5_cd * fccspacing
249  fccshifts(1, 4) = 0.0_cd
250  fccshifts(2, 4) = 0.5_cd * fccspacing
251  fccshifts(3, 4) = 0.5_cd * fccspacing
252 
253  middlepartid = 1 ! Always put middle particle as #1
254  a = 1 ! leave space for middle particle as particle #1
255  do i = 1, ncellsperside
256  latvec(1) = (i - 1) * fccspacing
257  do j = 1, ncellsperside
258  latvec(2) = (j - 1) * fccspacing
259  do k = 1, ncellsperside
260  latvec(3) = (k - 1) * fccspacing
261  do m = 1, 4
262  a = a + 1
263  coords(:, a) = latvec + fccshifts(:, m)
264  if ((i == ncellsperside / 2 + 1) &
265  .and. (j == ncellsperside / 2 + 1) &
266  .and. (k == ncellsperside / 2 + 1) &
267  .and. (m == 1)) &
268  then
269  ! put middle particle as #1
270  coords(:, 1) = latvec + fccshifts(:, m)
271  a = a - 1
272  end if
273  end do
274  end do
275  if (.not. periodic) then
276  ! Add in the remaining three faces
277  ! pos-x face
278  latvec(1) = ncellsperside * fccspacing
279  latvec(2) = (i - 1) * fccspacing
280  latvec(3) = (j - 1) * fccspacing
281  a = a + 1; coords(:, a) = latvec
282  a = a + 1; coords(:, a) = latvec + fccshifts(:, 4)
283  ! pos-y face
284  latvec(1) = (i - 1) * fccspacing
285  latvec(2) = ncellsperside * fccspacing
286  latvec(3) = (j - 1) * fccspacing
287  a = a + 1; coords(:, a) = latvec
288  a = a + 1; coords(:, a) = latvec + fccshifts(:, 3)
289  ! pos-z face
290  latvec(1) = (i - 1) * fccspacing
291  latvec(2) = (j - 1) * fccspacing
292  latvec(3) = ncellsperside * fccspacing
293  a = a + 1; coords(:, a) = latvec
294  a = a + 1; coords(:, a) = latvec + fccshifts(:, 2)
295  end if
296  end do
297  if (.not. periodic) then
298  ! Add in the remaining three edges
299  latvec(1) = (i - 1) * fccspacing
300  latvec(2) = ncellsperside * fccspacing
301  latvec(3) = ncellsperside * fccspacing
302  a = a + 1; coords(:, a) = latvec
303  latvec(1) = ncellsperside * fccspacing
304  latvec(2) = (i - 1) * fccspacing
305  latvec(3) = ncellsperside * fccspacing
306  a = a + 1; coords(:, a) = latvec
307  latvec(1) = ncellsperside * fccspacing
308  latvec(2) = ncellsperside * fccspacing
309  latvec(3) = (i - 1) * fccspacing
310  a = a + 1; coords(:, a) = latvec
311  end if
312  end do
313  if (.not. periodic) then
314  ! Add in the remaining corner
315  a = a + 1; coords(:, a) = ncellsperside * fccspacing
316  end if
317 
318  return
319 
320  end subroutine create_fcc_configuration
321 
322 end module mod_utility
323 
324 !*******************************************************************************
325 !**
326 !** PROGRAM vc_forces_numer_deriv
327 !**
328 !** KIM compliant program to perform numerical derivative check on a model
329 !**
330 !*******************************************************************************
331 
332 !-------------------------------------------------------------------------------
333 !
334 ! Main program
335 !
336 !-------------------------------------------------------------------------------
338  use, intrinsic :: iso_c_binding
339  use error
342  use mod_neighborlist
343  use mod_utility
344  implicit none
345  integer(c_int), parameter :: cd = c_double ! used for literal constants
346 
347  integer(c_int), parameter :: nCellsPerSide = 2
348  integer(c_int), parameter :: DIM = 3
349 
350  real(c_double), parameter :: cutpad = 0.75_cd
351  real(c_double), parameter :: FCCspacing = 5.260_cd
352  real(c_double), parameter :: min_spacing = 0.8 * fccspacing
353  real(c_double), parameter :: max_spacing = 1.2 * fccspacing
354  real(c_double), parameter :: spacing_incr = 0.025 * fccspacing
355  real(c_double) :: current_spacing
356  real(c_double) :: force_norm
357 
358  character(len=256, kind=c_char) :: modelname
359 
360  integer(c_int), parameter :: N = 4 * (ncellsperside)**3 &
361  + 6 * (ncellsperside)**2 &
362  + 3 * (ncellsperside) + 1
363 
364  type(neighobject_type), target :: neighObject
365  integer(c_int), allocatable, target :: neighborList(:, :)
366 
367  type(kim_model_handle_type) :: model_handle
368  type(kim_compute_arguments_handle_type) :: compute_arguments_handle
369  real(c_double) :: influence_distance
370  integer(c_int) :: number_of_neighbor_lists
371  real(c_double) :: cutoff
372  real(c_double) :: cutoffs(1)
373  integer(c_int) :: &
374  model_will_not_request_neighbors_of_noncontributing_particles(1)
375  integer(c_int), target :: particle_species_codes(n)
376  integer(c_int), target :: particle_contributing(n)
377  real(c_double), target :: energy
378  real(c_double), target :: coords(dim, n)
379  real(c_double), target :: forces(dim, n)
380  integer(c_int) i, j, ierr, ierr2
381 
382  integer(c_int) species_is_supported
383  integer(c_int) species_code
384  integer(c_int) requested_units_accepted
385  integer(c_int) number_of_model_routine_names
386  type(kim_model_routine_name_type) model_routine_name
387  integer(c_int) present
388  integer(c_int) required
389  type(kim_supported_extensions_type), target :: supported_extensions
390  character(len=KIM_MAX_EXTENSION_ID_LENGTH, kind=c_char) :: id_string
391 
392  integer :: middledum
393 
394  ! Initialize
395  ierr = 0
396  supported_extensions%number_of_supported_extensions = 0
397 
398  ! Get KIM Model name to use
399  print '("Please enter a valid KIM model name: ")'
400  read (*, *) modelname
401 
402  ! Create empty KIM object
403  !
404  call kim_model_create(kim_numbering_one_based, &
405  kim_length_unit_a, &
406  kim_energy_unit_ev, &
407  kim_charge_unit_e, &
408  kim_temperature_unit_k, &
409  kim_time_unit_ps, &
410  trim(modelname), &
411  requested_units_accepted, &
412  model_handle, ierr)
413  if (ierr /= 0) then
414  call my_error("kim_api_create")
415  end if
416 
417  ! check that we are compatible
418  if (requested_units_accepted == 0) then
419  call my_error("Must adapt to model units")
420  end if
421 
422  ! check that we know about all required routines
423  call kim_get_number_of_model_routine_names(number_of_model_routine_names)
424  do i = 1, number_of_model_routine_names
425  call kim_get_model_routine_name(i, model_routine_name, ierr)
426  if (ierr /= 0) then
427  call my_error("kim_get_model_routine_name")
428  end if
429  call kim_is_routine_present(model_handle, model_routine_name, present, &
430  required, ierr)
431  if (ierr /= 0) then
432  call my_error("kim_is_routine_present")
433  end if
434 
435  if ((present == 1) .and. (required == 1)) then
436  if (.not. ((model_routine_name == kim_model_routine_name_create) &
437  .or. (model_routine_name == &
438  kim_model_routine_name_compute_arguments_create) &
439  .or. (model_routine_name == kim_model_routine_name_compute) &
440  .or. (model_routine_name == kim_model_routine_name_refresh) &
441  .or. (model_routine_name == &
442  kim_model_routine_name_compute_arguments_destroy) &
443  .or. (model_routine_name == kim_model_routine_name_destroy))) &
444  then
445  call my_error("Unknown required ModelRoutineName found.")
446  end if
447  end if
448  end do
449 
450  ! check that model supports Ar
451  !
452  call kim_get_species_support_and_code( &
453  model_handle, kim_species_name_ar, species_is_supported, species_code, ierr)
454  if ((ierr /= 0) .or. (species_is_supported /= 1)) then
455  call my_error("Model does not support Ar")
456  end if
457 
458  ! Check supported extensions, if any
459  call kim_is_routine_present( &
460  model_handle, kim_model_routine_name_extension, present, required, ierr)
461  if (ierr /= 0) then
462  call my_error("Unable to get Extension present/required.")
463  end if
464  if (present /= 0) then
465  call kim_extension(model_handle, kim_supported_extensions_id, &
466  c_loc(supported_extensions), ierr)
467  if (ierr /= 0) then
468  call my_error("Error returned from kim_model_extension().")
469  end if
470  write (*, '(A,I2,A)') "Model Supports ", &
471  supported_extensions%number_of_supported_extensions, &
472  " Extensions:"
473  do i = 1, supported_extensions%number_of_supported_extensions
474  call kim_c_char_array_to_string( &
475  supported_extensions%supported_extension_id(:, i), id_string)
476  write (*, '(A,I2,A,A,A,A,I2)') " supportedExtensionID[", i, '] = "', &
477  trim(id_string), '" ', &
478  "which has required = ", &
479  supported_extensions%supported_extension_required(i)
480  end do
481  end if
482 
483  ! Best-practice is to check that the model is compatible
484  ! but we will skip it here
485 
486  ! create compute_arguments object
487  call kim_compute_arguments_create( &
488  model_handle, compute_arguments_handle, ierr)
489  if (ierr /= 0) then
490  call my_error("kim_model_compute_arguments_create")
491  end if
492 
493  ! register memory with the KIM system
494  ierr = 0
495  call kim_set_argument_pointer( &
496  compute_arguments_handle, kim_compute_argument_name_number_of_particles, &
497  n, ierr2)
498  ierr = ierr + ierr2
499  call kim_set_argument_pointer( &
500  compute_arguments_handle, &
501  kim_compute_argument_name_particle_species_codes, &
502  particle_species_codes, ierr2)
503  ierr = ierr + ierr2
504  call kim_set_argument_pointer( &
505  compute_arguments_handle, &
506  kim_compute_argument_name_particle_contributing, particle_contributing, &
507  ierr2)
508  ierr = ierr + ierr2
509  call kim_set_argument_pointer( &
510  compute_arguments_handle, kim_compute_argument_name_coordinates, coords, &
511  ierr2)
512  ierr = ierr + ierr2
513  call kim_set_argument_pointer( &
514  compute_arguments_handle, kim_compute_argument_name_partial_energy, &
515  energy, ierr2)
516  ierr = ierr + ierr2
517  call kim_set_argument_pointer( &
518  compute_arguments_handle, kim_compute_argument_name_partial_forces, &
519  forces, ierr2)
520  ierr = ierr + ierr2
521  if (ierr /= 0) then
522  call my_error("set_argument_pointer")
523  end if
524 
525  ! Allocate storage for neighbor lists
526  !
527  allocate (neighborlist(n + 1, n))
528  neighobject%neighbor_list_pointer = c_loc(neighborlist)
529  neighobject%number_of_particles = n
530 
531  ! Set pointer in KIM object to neighbor list routine and object
532  !
533  call kim_set_callback_pointer( &
534  compute_arguments_handle, kim_compute_callback_name_get_neighbor_list, &
535  kim_language_name_fortran, c_funloc(get_neigh), c_loc(neighobject), ierr)
536  if (ierr /= 0) then
537  call my_error("set_callback_pointer")
538  end if
539 
540  call kim_get_influence_distance(model_handle, influence_distance)
541  call kim_get_number_of_neighbor_lists(model_handle, &
542  number_of_neighbor_lists)
543  if (number_of_neighbor_lists /= 1) then
544  call my_error("too many neighbor lists")
545  end if
546  call kim_get_neighbor_list_values( &
547  model_handle, cutoffs, &
548  model_will_not_request_neighbors_of_noncontributing_particles, ierr)
549  if (ierr /= 0) then
550  call my_error("get_neighbor_list_values")
551  end if
552  cutoff = cutoffs(1)
553 
554  ! Setup cluster
555  !
556  do i = 1, n
557  particle_species_codes(i) = species_code
558  end do
559 
560  ! setup contributing particles
561  do i = 1, n
562  particle_contributing(i) = 1 ! every particle contributes
563  end do
564 
565  ! Print output header
566  !
567  print *
568  print *, 'This is Test : ex_test_Ar_fcc_cluster_fortran'
569  print *
570  print '(80(''-''))'
571  print '("Results for KIM Model : ",A)', trim(modelname)
572 
573  ! print header
574  print '(3A20)', "Energy", "Force Norm", "Lattice Spacing"
575  ! do the computations
576  current_spacing = min_spacing
577  do while (current_spacing < max_spacing)
578 
579  call create_fcc_configuration(current_spacing, ncellsperside, .false., &
580  coords, middledum)
581  ! Compute neighbor lists
582  call neigh_pure_cluster_neighborlist(.false., n, coords, &
583  (cutoff + cutpad), neighobject)
584 
585  ! Call model compute
586  call kim_compute(model_handle, compute_arguments_handle, ierr)
587  if (ierr /= 0) then
588  call my_error("kim_api_model_compute")
589  end if
590 
591  ! compue force_norm
592  force_norm = 0.0
593  do i = 1, n
594  do j = 1, dim
595  force_norm = force_norm + forces(j, i) * forces(j, i)
596  end do
597  end do
598  force_norm = sqrt(force_norm)
599 
600  ! Print results to screen
601  !
602  print '(3ES20.10)', energy, force_norm, current_spacing
603 
604  current_spacing = current_spacing + spacing_incr
605  end do
606 
607  ! Deallocate neighbor list object
608  deallocate (neighborlist)
609 
610  call kim_compute_arguments_destroy( &
611  model_handle, compute_arguments_handle, ierr)
612  if (ierr /= 0) then
613  call my_error("compute_arguments_destroy")
614  end if
615  call kim_model_destroy(model_handle)
616 
recursive subroutine, public get_neigh(data_object, number_of_neighbor_lists, cutoffs, neighbor_list_index, request, numnei, pnei1part, ierr)
recursive subroutine my_warning(message)
recursive subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, cutoff, neighObject)
recursive subroutine my_error(message)
recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, coords, MiddlePartId)
character(len= *, kind=c_char), parameter, public kim_supported_extensions_id
program ex_test_ar_fcc_cluster_fortran
The only standard extension defined by the KIM API.