kim-api  2.1.4-git+v2.1.3-git-1-g7847914a.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 ! CDDL HEADER START
3 !
4 ! The contents of this file are subject to the terms of the Common Development
5 ! and Distribution License Version 1.0 (the "License").
6 !
7 ! You can obtain a copy of the license at
8 ! http://www.opensource.org/licenses/CDDL-1.0. See the License for the
9 ! specific language governing permissions and limitations under the License.
10 !
11 ! When distributing Covered Code, include this CDDL HEADER in each file and
12 ! include the License file in a prominent location with the name LICENSE.CDDL.
13 ! If applicable, add the following below this CDDL HEADER, with the fields
14 ! enclosed by brackets "[]" replaced with your own identifying information:
15 !
16 ! Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
17 !
18 ! CDDL HEADER END
19 !
20 
21 !
22 ! Copyright (c) 2013--2019, Regents of the University of Minnesota.
23 ! All rights reserved.
24 !
25 ! Contributors:
26 ! Ellad B. Tadmor
27 ! Ryan S. Elliott
28 ! Stephen M. Whalen
29 !
30 
31 module error
32  use, intrinsic :: iso_c_binding
33  implicit none
34 
35  public
36 
37 contains
38  recursive subroutine my_error(message)
39  implicit none
40  character(len=*, kind=c_char), intent(in) :: message
41 
42  print *,"* Error : ", trim(message)
43  stop 1
44  end subroutine my_error
45 
46  recursive subroutine my_warning(message)
47  implicit none
48  character(len=*, kind=c_char), intent(in) :: message
49 
50  print *,"* Warning : ", trim(message)
51  end subroutine my_warning
52 end module error
53 
54 
55 !-------------------------------------------------------------------------------
56 !
57 ! module mod_neighborlist :
58 !
59 ! Module contains type and routines related to neighbor list calculation
60 !
61 !-------------------------------------------------------------------------------
62 
64 
65  use, intrinsic :: iso_c_binding
66 
67  public get_neigh
68 
69  type, bind(c) :: neighobject_type
70  real(c_double) :: cutoff
71  integer(c_int) :: number_of_particles
72  type(c_ptr) :: neighbor_list_pointer
73  end type neighobject_type
74 contains
75 
76 !-------------------------------------------------------------------------------
77 !
78 ! get_neigh neighbor list access function
79 !
80 ! This function implements Locator and Iterator mode
81 !
82 !-------------------------------------------------------------------------------
83 recursive subroutine get_neigh(data_object, number_of_neighbor_lists, cutoffs, &
84  neighbor_list_index, request, numnei, pnei1part, ierr) bind(c)
85  use error
86  implicit none
87 
88  !-- Transferred variables
89  type(c_ptr), value, intent(in) :: data_object
90  integer(c_int), value, intent(in) :: number_of_neighbor_lists
91  real(c_double), intent(in) :: cutoffs(*)
92  integer(c_int), value, intent(in) :: neighbor_list_index
93  integer(c_int), value, intent(in) :: request
94  integer(c_int), intent(out) :: numnei
95  type(c_ptr), intent(out) :: pnei1part
96  integer(c_int), intent(out) :: ierr
97 
98  !-- Local variables
99  integer(c_int) numberOfParticles
100  type(neighobject_type), pointer :: neighObject
101  integer(c_int), pointer :: neighborList(:,:)
102 
103  call c_f_pointer(data_object, neighobject)
104  numberofparticles = neighobject%number_of_particles
105  call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
106  [numberofparticles+1, numberofparticles])
107 
108  if (number_of_neighbor_lists /= 1) then
109  call my_warning("invalid number of cutoffs")
110  ierr = 1
111  return
112  endif
113 
114  if (cutoffs(1) > neighobject%cutoff) then
115  call my_warning("neighbor list cutoff too small for model cutoff")
116  ierr = 1
117  return
118  endif
119 
120  if (neighbor_list_index /= 1) then
121  call my_warning("wrong list index")
122  ierr = 1
123  return
124  endif
125 
126  if ( (request.gt.numberofparticles) .or. (request.lt.1)) then
127  print *, request
128  call my_warning("Invalid part ID in get_neigh")
129  ierr = 1
130  return
131  endif
132 
133  ! set the returned number of neighbors for the returned part
134  numnei = neighborlist(1,request)
135 
136  ! set the location for the returned neighbor list
137  pnei1part = c_loc(neighborlist(2,request))
138 
139  ierr = 0
140  return
141 end subroutine get_neigh
142 
143 end module mod_neighborlist
144 
145 
147  implicit none
148  public
149 
150 contains
151 
152 !-------------------------------------------------------------------------------
153 !
154 ! NEIGH_PURE_cluster_neighborlist
155 !
156 !-------------------------------------------------------------------------------
157 recursive subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, &
158  coords, cutoff, neighObject)
159  use, intrinsic :: iso_c_binding
160  use mod_neighborlist
161  implicit none
162 
163  !-- Transferred variables
164  logical, intent(in) :: half
165  integer(c_int), intent(in) :: numberOfParticles
166  real(c_double), dimension(3,numberOfParticles), &
167  intent(in) :: coords
168  real(c_double), intent(in) :: cutoff
169  type(neighobject_type), intent(inout) :: neighObject
170 
171  !-- Local variables
172  integer(c_int) i, j, a
173  real(c_double) dx(3)
174  real(c_double) r2
175  real(c_double) cutoff2
176  integer(c_int), pointer :: neighborList(:,:)
177 
178  call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
179  [numberofparticles+1, numberofparticles])
180 
181  neighobject%cutoff = cutoff
182 
183  cutoff2 = cutoff**2
184 
185  do i=1,numberofparticles
186  a = 1
187  do j=1,numberofparticles
188  dx(:) = coords(:, j) - coords(:, i)
189  r2 = dot_product(dx, dx)
190  if (r2.le.cutoff2) then
191  ! part j is a neighbor of part i
192  if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) ) then
193  a = a+1
194  neighborlist(a,i) = j
195  endif
196  endif
197  enddo
198  ! part i has a-1 neighbors
199  neighborlist(1,i) = a-1
200  enddo
201 
202  return
203 
204 end subroutine neigh_pure_cluster_neighborlist
205 
206 !-------------------------------------------------------------------------------
207 !
208 ! create_FCC_configuration subroutine
209 !
210 ! creates a cubic configuration of FCC particles with lattice spacing
211 ! `FCCspacing' and `nCellsPerSide' cells along each direction.
212 !
213 ! With periodic==.true. this will result in a total number of particles equal
214 ! to 4*(nCellsPerSide)**3 + 6*(nCellsPerSide)**2 + 3*(nCellsPerSide) + 1
215 !
216 ! With periodic==.false. this will result in a total number of particles equal
217 ! to 4*(nCellsPerSide)**3
218 !
219 ! Returns the Id of the particle situated in the middle of the configuration
220 ! (this particle will have the most neighbors.)
221 !
222 !-------------------------------------------------------------------------------
223 recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, &
224  periodic, coords, MiddlePartId)
225  use, intrinsic :: iso_c_binding
226  implicit none
227  integer(c_int), parameter :: cd = c_double ! used for literal constants
228 
229  !-- Transferred variables
230  real(c_double), intent(in) :: FCCspacing
231  integer(c_int), intent(in) :: nCellsPerSide
232  logical, intent(in) :: periodic
233  real(c_double), intent(out) :: coords(3,*)
234  integer(c_int), intent(out) :: MiddlePartId
235  !
236  ! cluster setup variables
237  !
238  real(c_double) FCCshifts(3,4)
239  real(c_double) latVec(3)
240  integer(c_int) a, i, j, k, m
241 
242  ! Create a cubic FCC cluster
243  !
244  fccshifts(1,1) = 0.0_cd
245  fccshifts(2,1) = 0.0_cd
246  fccshifts(3,1) = 0.0_cd
247  fccshifts(1,2) = 0.5_cd*fccspacing
248  fccshifts(2,2) = 0.5_cd*fccspacing
249  fccshifts(3,2) = 0.0_cd
250  fccshifts(1,3) = 0.5_cd*fccspacing
251  fccshifts(2,3) = 0.0_cd
252  fccshifts(3,3) = 0.5_cd*fccspacing
253  fccshifts(1,4) = 0.0_cd
254  fccshifts(2,4) = 0.5_cd*fccspacing
255  fccshifts(3,4) = 0.5_cd*fccspacing
256 
257  middlepartid = 1 ! Always put middle particle as #1
258  a = 1 ! leave space for middle particle as particle #1
259  do i=1,ncellsperside
260  latvec(1) = (i-1)*fccspacing
261  do j=1,ncellsperside
262  latvec(2) = (j-1)*fccspacing
263  do k=1,ncellsperside
264  latvec(3) = (k-1)*fccspacing
265  do m=1,4
266  a = a+1
267  coords(:,a) = latvec + fccshifts(:,m)
268  if ((i.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
269  (k.eq.ncellsperside/2+1) .and. (m.eq.1)) then
270  coords(:,1) = latvec + fccshifts(:,m) ! put middle particle as #1
271  a = a - 1
272  endif
273  enddo
274  enddo
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  endif
296  enddo
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  endif
312  enddo
313  if (.not. periodic) then
314  ! Add in the remaining corner
315  a = a+1; coords(:,a) = ncellsperside*fccspacing
316  endif
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 :: &
361  N = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
362 
363  type(neighobject_type), target :: neighObject
364  integer(c_int), allocatable, target :: neighborList(:,:)
365 
366  type(kim_model_handle_type) :: model_handle
367  type(kim_compute_arguments_handle_type) :: compute_arguments_handle
368  real(c_double) :: influence_distance
369  integer(c_int) :: number_of_neighbor_lists
370  real(c_double) :: cutoff
371  real(c_double) :: cutoffs(1)
372  integer(c_int) :: &
373  model_will_not_request_neighbors_of_noncontributing_particles(1)
374  integer(c_int), target :: particle_species_codes(n)
375  integer(c_int), target :: particle_contributing(n)
376  real(c_double), target :: energy
377  real(c_double), target :: coords(dim, n)
378  real(c_double), target :: forces(dim, n)
379  integer(c_int) i, j, ierr, ierr2
380 
381  integer(c_int) species_is_supported
382  integer(c_int) species_code
383  integer(c_int) requested_units_accepted
384  integer(c_int) number_of_model_routine_names
385  type(kim_model_routine_name_type) model_routine_name
386  integer(c_int) present
387  integer(c_int) required
388  type(kim_supported_extensions_type), target :: supported_extensions
389  character(len=KIM_MAX_EXTENSION_ID_LENGTH, kind=c_char) :: id_string
390 
391  integer :: middledum
392 
393  ! Initialize
394  ierr = 0
395  supported_extensions%number_of_supported_extensions = 0
396 
397  ! Get KIM Model name to use
398  print '("Please enter a valid KIM model name: ")'
399  read(*,*) modelname
400 
401  ! Create empty KIM object
402  !
403  call kim_model_create(kim_numbering_one_based, &
404  kim_length_unit_a, &
405  kim_energy_unit_ev, &
406  kim_charge_unit_e, &
407  kim_temperature_unit_k, &
408  kim_time_unit_ps, &
409  trim(modelname), &
410  requested_units_accepted, &
411  model_handle, ierr)
412  if (ierr /= 0) then
413  call my_error("kim_api_create")
414  endif
415 
416  ! check that we are compatible
417  if (requested_units_accepted == 0) then
418  call my_error("Must adapt to model units")
419  end if
420 
421  ! check that we know about all required routines
422  call kim_get_number_of_model_routine_names(number_of_model_routine_names)
423  do i=1,number_of_model_routine_names
424  call kim_get_model_routine_name(i, model_routine_name, ierr)
425  if (ierr /= 0) then
426  call my_error("kim_get_model_routine_name")
427  endif
428  call kim_is_routine_present(model_handle, model_routine_name, present, &
429  required, ierr)
430  if (ierr /= 0) then
431  call my_error("kim_is_routine_present")
432  endif
433 
434  if ((present == 1) .and. (required == 1)) then
435  if (.not. ((model_routine_name == kim_model_routine_name_create) &
436  .or. (model_routine_name == &
437  kim_model_routine_name_compute_arguments_create) &
438  .or. (model_routine_name == kim_model_routine_name_compute) &
439  .or. (model_routine_name == kim_model_routine_name_refresh) &
440  .or. (model_routine_name == &
441  kim_model_routine_name_compute_arguments_destroy) &
442  .or. (model_routine_name == kim_model_routine_name_destroy))) then
443 
444  call my_error("Unknown required ModelRoutineName found.")
445  endif
446  endif
447  enddo
448 
449  ! check that model supports Ar
450  !
451  call kim_get_species_support_and_code(model_handle, &
452  kim_species_name_ar, species_is_supported, species_code, ierr)
453  if ((ierr /= 0) .or. (species_is_supported /= 1)) then
454  call my_error("Model does not support Ar")
455  endif
456 
457  ! Check supported extensions, if any
458  call kim_is_routine_present(model_handle, &
459  kim_model_routine_name_extension, present, required, ierr)
460  if (ierr /= 0) then
461  call my_error("Unable to get Extension present/required.")
462  endif
463  if (present /= 0) then
464  call kim_extension(model_handle, kim_supported_extensions_id, &
465  c_loc(supported_extensions), ierr)
466  if (ierr /= 0) then
467  call my_error("Error returned from kim_model_extension().")
468  endif
469  write (*,'(A,I2,A)') "Model Supports ", &
470  supported_extensions%number_of_supported_extensions, &
471  " Extensions:"
472  do i = 1, supported_extensions%number_of_supported_extensions
473  call kim_c_char_array_to_string(&
474  supported_extensions%supported_extension_id(:,i), id_string)
475  write (*,'(A,I2,A,A,A,A,I2)') " supportedExtensionID[", i, '] = "', &
476  trim(id_string), '" ', &
477  "which has required = ", &
478  supported_extensions%supported_extension_required(i)
479  end do
480  endif
481 
482  ! Best-practice is to check that the model is compatible
483  ! but we will skip it here
484 
485  ! create compute_arguments object
486  call kim_compute_arguments_create( &
487  model_handle, compute_arguments_handle, ierr)
488  if (ierr /= 0) then
489  call my_error("kim_model_compute_arguments_create")
490  endif
491 
492  ! register memory with the KIM system
493  ierr = 0
494  call kim_set_argument_pointer(compute_arguments_handle, &
495  kim_compute_argument_name_number_of_particles, n, ierr2)
496  ierr = ierr + ierr2
497  call kim_set_argument_pointer(compute_arguments_handle, &
498  kim_compute_argument_name_particle_species_codes, particle_species_codes, &
499  ierr2)
500  ierr = ierr + ierr2
501  call kim_set_argument_pointer(compute_arguments_handle, &
502  kim_compute_argument_name_particle_contributing, particle_contributing, &
503  ierr2)
504  ierr = ierr + ierr2
505  call kim_set_argument_pointer(compute_arguments_handle, &
506  kim_compute_argument_name_coordinates, coords, ierr2)
507  ierr = ierr + ierr2
508  call kim_set_argument_pointer(compute_arguments_handle, &
509  kim_compute_argument_name_partial_energy, energy, ierr2)
510  ierr = ierr + ierr2
511  call kim_set_argument_pointer(compute_arguments_handle, &
512  kim_compute_argument_name_partial_forces, forces, ierr2)
513  ierr = ierr + ierr2
514  if (ierr /= 0) then
515  call my_error("set_argument_pointer")
516  endif
517 
518  ! Allocate storage for neighbor lists
519  !
520  allocate(neighborlist(n+1,n))
521  neighobject%neighbor_list_pointer = c_loc(neighborlist)
522  neighobject%number_of_particles = n
523 
524  ! Set pointer in KIM object to neighbor list routine and object
525  !
526  call kim_set_callback_pointer(compute_arguments_handle, &
527  kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
528  c_funloc(get_neigh), c_loc(neighobject), ierr)
529  if (ierr /= 0) then
530  call my_error("set_callback_pointer")
531  end if
532 
533  call kim_get_influence_distance(model_handle, influence_distance)
534  call kim_get_number_of_neighbor_lists(model_handle, &
535  number_of_neighbor_lists)
536  if (number_of_neighbor_lists /= 1) then
537  call my_error("too many neighbor lists")
538  endif
539  call kim_get_neighbor_list_values(model_handle, cutoffs, &
540  model_will_not_request_neighbors_of_noncontributing_particles, ierr)
541  if (ierr /= 0) then
542  call my_error("get_neighbor_list_values")
543  end if
544  cutoff = cutoffs(1)
545 
546  ! Setup cluster
547  !
548  do i=1,n
549  particle_species_codes(i) = species_code
550  enddo
551 
552  ! setup contributing particles
553  do i=1,n
554  particle_contributing(i) = 1 ! every particle contributes
555  enddo
556 
557  ! Print output header
558  !
559  print *
560  print *,'This is Test : ex_test_Ar_fcc_cluster_fortran'
561  print *
562  print '(80(''-''))'
563  print '("Results for KIM Model : ",A)', trim(modelname)
564 
565  ! print header
566  print '(3A20)', "Energy", "Force Norm", "Lattice Spacing"
567  ! do the computations
568  current_spacing = min_spacing
569  do while (current_spacing < max_spacing)
570 
571  call create_fcc_configuration(current_spacing, ncellsperside, .false., &
572  coords, middledum)
573  ! Compute neighbor lists
574  call neigh_pure_cluster_neighborlist(.false., n, coords, (cutoff+cutpad), &
575  neighobject)
576 
577  ! Call model compute
578  call kim_compute(model_handle, compute_arguments_handle, ierr)
579  if (ierr /= 0) then
580  call my_error("kim_api_model_compute")
581  endif
582 
583  ! compue force_norm
584  force_norm = 0.0;
585  do i = 1, n
586  do j = 1, dim
587  force_norm = force_norm + forces(j,i)*forces(j,i)
588  end do
589  end do
590  force_norm = sqrt(force_norm)
591 
592  ! Print results to screen
593  !
594  print '(3ES20.10)', energy, force_norm, current_spacing
595 
596  current_spacing = current_spacing + spacing_incr
597  enddo
598 
599  ! Deallocate neighbor list object
600  deallocate( neighborlist )
601 
602  call kim_compute_arguments_destroy(&
603  model_handle, compute_arguments_handle, ierr)
604  if (ierr /= 0) then
605  call my_error("compute_arguments_destroy")
606  endif
607  call kim_model_destroy(model_handle)
608 
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.