KIM API V2
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--2018, 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  public
35 
36 contains
37  subroutine my_error(message, line, file)
38  implicit none
39  character(len=*, kind=c_char), intent(in) :: message
40  integer, intent(in) :: line
41  character(len=*, kind=c_char), intent(in) :: file
42 
43  print *,"* Error : '", trim(message), "' ", line, ":", &
44  trim(file)
45  stop
46  end subroutine my_error
47 
48  subroutine my_warning(message, line, file)
49  implicit none
50  character(len=*, kind=c_char), intent(in) :: message
51  integer, intent(in) :: line
52  character(len=*, kind=c_char), intent(in) :: file
53 
54  print *,"* Error : '", trim(message), "' ", line, ":", &
55  trim(file)
56  end subroutine my_warning
57 end module error
58 
59 
60 !-------------------------------------------------------------------------------
61 !
62 ! module mod_neighborlist :
63 !
64 ! Module contains type and routines related to neighbor list calculation
65 !
66 !-------------------------------------------------------------------------------
67 
68 module mod_neighborlist
69 
70  use, intrinsic :: iso_c_binding
71 
72  public get_neigh
73 
74  type neighobject_type
75  real(c_double) :: cutoff
76  integer(c_int) :: number_of_particles
77  integer(c_int), pointer :: neighborList(:,:)
78  end type neighobject_type
79 contains
80 
81 !-------------------------------------------------------------------------------
82 !
83 ! get_neigh neighbor list access function
84 !
85 ! This function implements Locator and Iterator mode
86 !
87 !-------------------------------------------------------------------------------
88 subroutine get_neigh(data_object, number_of_neighbor_lists, cutoffs, &
89  neighbor_list_index, request, numnei, pnei1part, ierr) bind(c)
90  use error
91  implicit none
92 
93  !-- Transferred variables
94  type(c_ptr), value, intent(in) :: data_object
95  integer(c_int), value, intent(in) :: number_of_neighbor_lists
96  real(c_double), intent(in) :: cutoffs(number_of_neighbor_lists)
97  integer(c_int), value, intent(in) :: neighbor_list_index
98  integer(c_int), value, intent(in) :: request
99  integer(c_int), intent(out) :: numnei
100  type(c_ptr), intent(out) :: pnei1part
101  integer(c_int), intent(out) :: ierr
102 
103  !-- Local variables
104  integer(c_int), parameter :: DIM = 3
105  integer(c_int) numberOfParticles
106  type(neighObject_type), pointer :: neighObject
107 
108  call c_f_pointer(data_object, neighobject)
109 
110  if (number_of_neighbor_lists /= 1) then
111  call my_warning("invalid number of cutoffs", __line__, __file__)
112  ierr = 1
113  return
114  endif
115 
116  if (cutoffs(1) > neighobject%cutoff) then
117  call my_warning("neighbor list cutoff too small for model cutoff", &
118  __line__, __file__)
119  ierr = 1
120  return
121  endif
122 
123  if (neighbor_list_index /= 1) then
124  call my_warning("wrong list index", __line__, __file__)
125  ierr = 1
126  return
127  endif
128 
129  numberofparticles = neighobject%number_of_particles
130 
131  if ( (request.gt.numberofparticles) .or. (request.lt.1)) then
132  print *, request
133  call my_warning("Invalid part ID in get_neigh", &
134  __line__, __file__)
135  ierr = 1
136  return
137  endif
138 
139  ! set the returned number of neighbors for the returned part
140  numnei = neighobject%neighborList(1,request)
141 
142  ! set the location for the returned neighbor list
143  pnei1part = c_loc(neighobject%neighborList(2,request))
144 
145  ierr = 0
146  return
147 end subroutine get_neigh
148 
149 end module mod_neighborlist
150 
151 
152 module mod_utility
153  implicit none
154  public
155 
156 contains
157 
158 !-------------------------------------------------------------------------------
159 !
160 ! NEIGH_PURE_cluster_neighborlist
161 !
162 !-------------------------------------------------------------------------------
163 subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, &
164  cutoff, neighObject)
165  use, intrinsic :: iso_c_binding
166  use mod_neighborlist
167  implicit none
168 
169  !-- Transferred variables
170  logical, intent(in) :: half
171  integer(c_int), intent(in) :: numberOfParticles
172  real(c_double), dimension(3,numberOfParticles), &
173  intent(in) :: coords
174  real(c_double), intent(in) :: cutoff
175  type(neighObject_type), intent(inout) :: neighObject
176 
177  !-- Local variables
178  integer(c_int) i, j, a
179  real(c_double) dx(3)
180  real(c_double) r2
181  real(c_double) cutoff2
182 
183  neighobject%cutoff = cutoff
184 
185  cutoff2 = cutoff**2
186 
187  do i=1,numberofparticles
188  a = 1
189  do j=1,numberofparticles
190  dx(:) = coords(:, j) - coords(:, i)
191  r2 = dot_product(dx, dx)
192  if (r2.le.cutoff2) then
193  ! part j is a neighbor of part i
194  if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) ) then
195  a = a+1
196  neighobject%neighborList(a,i) = j
197  endif
198  endif
199  enddo
200  ! part i has a-1 neighbors
201  neighobject%neighborList(1,i) = a-1
202  enddo
203 
204  return
205 
206 end subroutine neigh_pure_cluster_neighborlist
207 
208 !-------------------------------------------------------------------------------
209 !
210 ! create_FCC_configuration subroutine
211 !
212 ! creates a cubic configuration of FCC particles with lattice spacing
213 ! `FCCspacing' and `nCellsPerSide' cells along each direction.
214 !
215 ! With periodic==.true. this will result in a total number of particles equal
216 ! to 4*(nCellsPerSide)**3 + 6*(nCellsPerSide)**2 + 3*(nCellsPerSide) + 1
217 !
218 ! With periodic==.false. this will result in a total number of particles equal
219 ! to 4*(nCellsPerSide)**3
220 !
221 ! Returns the Id of the particle situated in the middle of the configuration
222 ! (this particle will have the most neighbors.)
223 !
224 !-------------------------------------------------------------------------------
225 subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, &
226  coords, MiddlePartId)
227  use, intrinsic :: iso_c_binding
228  implicit none
229  integer(c_int), parameter :: cd = c_double ! used for literal constants
230 
231  !-- Transferred variables
232  real(c_double), intent(in) :: FCCspacing
233  integer(c_int), intent(in) :: nCellsPerSide
234  logical, intent(in) :: periodic
235  real(c_double), intent(out) :: coords(3,*)
236  integer(c_int), intent(out) :: MiddlePartId
237  !
238  ! cluster setup variables
239  !
240  real(c_double) FCCshifts(3,4)
241  real(c_double) latVec(3)
242  integer(c_int) a, i, j, k, m
243 
244  ! Create a cubic FCC cluster
245  !
246  fccshifts(1,1) = 0.0_cd
247  fccshifts(2,1) = 0.0_cd
248  fccshifts(3,1) = 0.0_cd
249  fccshifts(1,2) = 0.5_cd*fccspacing
250  fccshifts(2,2) = 0.5_cd*fccspacing
251  fccshifts(3,2) = 0.0_cd
252  fccshifts(1,3) = 0.5_cd*fccspacing
253  fccshifts(2,3) = 0.0_cd
254  fccshifts(3,3) = 0.5_cd*fccspacing
255  fccshifts(1,4) = 0.0_cd
256  fccshifts(2,4) = 0.5_cd*fccspacing
257  fccshifts(3,4) = 0.5_cd*fccspacing
258 
259  middlepartid = 1 ! Always put middle particle as #1
260  a = 1 ! leave space for middle particle as particle #1
261  do i=1,ncellsperside
262  latvec(1) = (i-1)*fccspacing
263  do j=1,ncellsperside
264  latvec(2) = (j-1)*fccspacing
265  do k=1,ncellsperside
266  latvec(3) = (k-1)*fccspacing
267  do m=1,4
268  a = a+1
269  coords(:,a) = latvec + fccshifts(:,m)
270  if ((i.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
271  (k.eq.ncellsperside/2+1) .and. (m.eq.1)) then
272  coords(:,1) = latvec + fccshifts(:,m) ! put middle particle as #1
273  a = a - 1
274  endif
275  enddo
276  enddo
277  if (.not. periodic) then
278  ! Add in the remaining three faces
279  ! pos-x face
280  latvec(1) = ncellsperside*fccspacing
281  latvec(2) = (i-1)*fccspacing
282  latvec(3) = (j-1)*fccspacing
283  a = a+1; coords(:,a) = latvec
284  a = a+1; coords(:,a) = latvec + fccshifts(:,4)
285  ! pos-y face
286  latvec(1) = (i-1)*fccspacing
287  latvec(2) = ncellsperside*fccspacing
288  latvec(3) = (j-1)*fccspacing
289  a = a+1; coords(:,a) = latvec
290  a = a+1; coords(:,a) = latvec + fccshifts(:,3)
291  ! pos-z face
292  latvec(1) = (i-1)*fccspacing
293  latvec(2) = (j-1)*fccspacing
294  latvec(3) = ncellsperside*fccspacing
295  a = a+1; coords(:,a) = latvec
296  a = a+1; coords(:,a) = latvec + fccshifts(:,2)
297  endif
298  enddo
299  if (.not. periodic) then
300  ! Add in the remaining three edges
301  latvec(1) = (i-1)*fccspacing
302  latvec(2) = ncellsperside*fccspacing
303  latvec(3) = ncellsperside*fccspacing
304  a = a+1; coords(:,a) = latvec
305  latvec(1) = ncellsperside*fccspacing
306  latvec(2) = (i-1)*fccspacing
307  latvec(3) = ncellsperside*fccspacing
308  a = a+1; coords(:,a) = latvec
309  latvec(1) = ncellsperside*fccspacing
310  latvec(2) = ncellsperside*fccspacing
311  latvec(3) = (i-1)*fccspacing
312  a = a+1; coords(:,a) = latvec
313  endif
314  enddo
315  if (.not. periodic) then
316  ! Add in the remaining corner
317  a = a+1; coords(:,a) = ncellsperside*fccspacing
318  endif
319 
320  return
321 
322 end subroutine create_fcc_configuration
323 
324 end module mod_utility
325 
326 !*******************************************************************************
327 !**
328 !** PROGRAM vc_forces_numer_deriv
329 !**
330 !** KIM compliant program to perform numerical derivative check on a model
331 !**
332 !*******************************************************************************
333 
334 !-------------------------------------------------------------------------------
335 !
336 ! Main program
337 !
338 !-------------------------------------------------------------------------------
340  use, intrinsic :: iso_c_binding
341  use error
343  use mod_neighborlist
344  use mod_utility
345  implicit none
346  integer(c_int), parameter :: cd = c_double ! used for literal constants
347 
348  integer(c_int), parameter :: ncellsperside = 2
349  integer(c_int), parameter :: dim = 3
350 
351  real(c_double), parameter :: cutpad = 0.75_cd
352  real(c_double), parameter :: fccspacing = 5.260_cd
353  real(c_double), parameter :: min_spacing = 0.8*fccspacing
354  real(c_double), parameter :: max_spacing = 1.2*fccspacing
355  real(c_double), parameter :: spacing_incr = 0.025*fccspacing
356  real(c_double) :: current_spacing
357  real(c_double) :: force_norm
358 
359  character(len=256, kind=c_char) :: modelname
360 
361  integer(c_int), parameter :: &
362  n = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
363 
364  type(neighobject_type), target :: neighobject
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) :: padding_neighbor_hints(1)
373  integer(c_int) :: half_list_hints(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 
385  integer :: middledum
386 
387  ! Initialize error flag
388  ierr = 0
389 
390  ! Get KIM Model name to use
391  print '("Please enter a valid KIM model name: ")'
392  read(*,*) modelname
393 
394  ! Print output header
395  !
396  print *
397  print *,'This is Test : ex_test_Ar_fcc_cluster_fortran'
398  print *
399  print '(80(''-''))'
400  print '("Results for KIM Model : ",A)', trim(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", __line__, __file__)
415  endif
416 
417  ! check that we are compatible
418  if (requested_units_accepted == 0) then
419  call my_error("Must adapt to model units", __line__, __file__)
420  end if
421 
422  ! check that model supports Ar
423  !
424  call kim_model_get_species_support_and_code(model_handle, &
425  kim_species_name_ar, species_is_supported, species_code, ierr)
426  if ((ierr /= 0) .or. (species_is_supported /= 1)) then
427  call my_error("Model does not support Ar", __line__, __file__)
428  endif
429 
430  ! Best-practice is to check that the model is compatible
431  ! but we will skip it here
432 
433  ! create compute_arguments object
434  call kim_model_compute_arguments_create( &
435  model_handle, compute_arguments_handle, ierr)
436  if (ierr /= 0) then
437  call my_error("kim_model_compute_arguments_create", __line__, __file__)
438  endif
439 
440  ! register memory with the KIM system
441  ierr = 0
442  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
443  kim_compute_argument_name_number_of_particles, n, ierr2)
444  ierr = ierr + ierr2
445  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
446  kim_compute_argument_name_particle_species_codes, particle_species_codes, &
447  ierr2)
448  ierr = ierr + ierr2
449  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
450  kim_compute_argument_name_particle_contributing, particle_contributing, &
451  ierr2)
452  ierr = ierr + ierr2
453  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
454  kim_compute_argument_name_coordinates, coords, ierr2)
455  ierr = ierr + ierr2
456  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
457  kim_compute_argument_name_partial_energy, energy, ierr2)
458  ierr = ierr + ierr2
459  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
460  kim_compute_argument_name_partial_forces, forces, ierr2)
461  ierr = ierr + ierr2
462  if (ierr /= 0) then
463  call my_error("set_argument_pointer", __line__, __file__)
464  endif
465 
466  ! Allocate storage for neighbor lists
467  !
468  allocate(neighobject%neighborList(n+1,n))
469  neighobject%number_of_particles = n
470 
471  ! Set pointer in KIM object to neighbor list routine and object
472  !
473  call kim_compute_arguments_set_callback_pointer(compute_arguments_handle, &
474  kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
475  c_funloc(get_neigh), c_loc(neighobject), ierr)
476  if (ierr /= 0) then
477  call my_error("set_callback_pointer", __line__, __file__)
478  end if
479 
480  call kim_model_get_influence_distance(model_handle, influence_distance)
481  call kim_model_get_number_of_neighbor_lists(model_handle, &
482  number_of_neighbor_lists)
483  if (number_of_neighbor_lists /= 1) then
484  call my_error("too many neighbor lists", __line__, __file__)
485  endif
486  call kim_model_get_neighbor_list_values(model_handle, cutoffs, &
487  padding_neighbor_hints, half_list_hints, ierr)
488  if (ierr /= 0) then
489  call my_error("get_neighbor_list_values", __line__, __file__)
490  end if
491  cutoff = cutoffs(1)
492 
493  ! Setup cluster
494  !
495  do i=1,n
496  particle_species_codes(i) = species_code
497  enddo
498 
499  ! setup contributing particles
500  do i=1,n
501  particle_contributing(i) = 1 ! every particle contributes
502  enddo
503 
504  ! print header
505  print '(3A20)', "Energy", "Force Norm", "Lattice Spacing"
506  ! do the computations
507  current_spacing = min_spacing
508  do while (current_spacing < max_spacing)
509 
510  call create_fcc_configuration(current_spacing, ncellsperside, .false., &
511  coords, middledum)
512  ! Compute neighbor lists
513  call neigh_pure_cluster_neighborlist(.false., n, coords, (cutoff+cutpad), &
514  neighobject)
515 
516  ! Call model compute
517  call kim_model_compute(model_handle, compute_arguments_handle, ierr)
518  if (ierr /= 0) then
519  call my_error("kim_api_model_compute", __line__, __file__)
520  endif
521 
522  ! compue force_norm
523  force_norm = 0.0;
524  do i = 1, n
525  do j = 1, dim
526  force_norm = force_norm + forces(j,i)*forces(j,i)
527  end do
528  end do
529  force_norm = sqrt(force_norm)
530 
531  ! Print results to screen
532  !
533  print '(3ES20.10)', energy, force_norm, current_spacing
534 
535  current_spacing = current_spacing + spacing_incr
536  enddo
537 
538  ! Deallocate neighbor list object
539  deallocate( neighobject%neighborList )
540 
541  call kim_model_compute_arguments_destroy(&
542  model_handle, compute_arguments_handle, ierr)
543  if (ierr /= 0) then
544  call my_error("compute_arguments_destroy", __line__, __file__)
545  endif
546  call kim_model_destroy(model_handle)
547 
subroutine my_error(message, line, file)
subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, cutoff, neighObject)
subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, coords, MiddlePartId)
program ex_test_ar_fcc_cluster_fortran
subroutine my_warning(message, line, file)
integer(c_int) function, public get_neigh(pkim, mode, request, part, numnei, pnei1part, pRij)