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_cutoffs, 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_cutoffs
96  real(c_double), intent(in) :: cutoffs(number_of_cutoffs)
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_cutoffs /= 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_cutoffs
370  real(c_double) :: cutoff
371  real(c_double) :: cutoffs(1)
372  integer(c_int), target :: particle_species_codes(n)
373  integer(c_int), target :: particle_contributing(n)
374  real(c_double), target :: energy
375  real(c_double), target :: coords(dim, n)
376  real(c_double), target :: forces(dim, n)
377  integer(c_int) i, j, ierr, ierr2
378 
379  integer(c_int) species_is_supported
380  integer(c_int) species_code
381  integer(c_int) requested_units_accepted
382 
383  integer :: middledum
384 
385  ! Initialize error flag
386  ierr = 0
387 
388  ! Get KIM Model name to use
389  print '("Please enter a valid KIM model name: ")'
390  read(*,*) modelname
391 
392  ! Print output header
393  !
394  print *
395  print *,'This is Test : ex_test_Ar_fcc_cluster_fortran'
396  print *
397  print '(80(''-''))'
398  print '("Results for KIM Model : ",A)', trim(modelname)
399 
400  ! Create empty KIM object
401  !
402  call kim_model_create(kim_numbering_one_based, &
403  kim_length_unit_a, &
404  kim_energy_unit_ev, &
405  kim_charge_unit_e, &
406  kim_temperature_unit_k, &
407  kim_time_unit_ps, &
408  trim(modelname), &
409  requested_units_accepted, &
410  model_handle, ierr)
411  if (ierr /= 0) then
412  call my_error("kim_api_create", __line__, __file__)
413  endif
414 
415  ! check that we are compatible
416  if (requested_units_accepted == 0) then
417  call my_error("Must adapt to model units", __line__, __file__)
418  end if
419 
420  ! check that model supports Ar
421  !
422  call kim_model_get_species_support_and_code(model_handle, &
423  kim_species_name_ar, species_is_supported, species_code, ierr)
424  if ((ierr /= 0) .or. (species_is_supported /= 1)) then
425  call my_error("Model does not support Ar", __line__, __file__)
426  endif
427 
428  ! Best-practice is to check that the model is compatible
429  ! but we will skip it here
430 
431  ! create compute_arguments object
432  call kim_model_compute_arguments_create( &
433  model_handle, compute_arguments_handle, ierr)
434  if (ierr /= 0) then
435  call my_error("kim_model_compute_arguments_create", __line__, __file__)
436  endif
437 
438  ! register memory with the KIM system
439  ierr = 0
440  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
441  kim_compute_argument_name_number_of_particles, n, ierr2)
442  ierr = ierr + ierr2
443  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
444  kim_compute_argument_name_particle_species_codes, particle_species_codes, &
445  ierr2)
446  ierr = ierr + ierr2
447  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
448  kim_compute_argument_name_particle_contributing, particle_contributing, &
449  ierr2)
450  ierr = ierr + ierr2
451  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
452  kim_compute_argument_name_coordinates, coords, ierr2)
453  ierr = ierr + ierr2
454  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
455  kim_compute_argument_name_partial_energy, energy, ierr2)
456  ierr = ierr + ierr2
457  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
458  kim_compute_argument_name_partial_forces, forces, ierr2)
459  ierr = ierr + ierr2
460  if (ierr /= 0) then
461  call my_error("set_argument_pointer", __line__, __file__)
462  endif
463 
464  ! Allocate storage for neighbor lists
465  !
466  allocate(neighobject%neighborList(n+1,n))
467  neighobject%number_of_particles = n
468 
469  ! Set pointer in KIM object to neighbor list routine and object
470  !
471  call kim_compute_arguments_set_callback_pointer(compute_arguments_handle, &
472  kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
473  c_funloc(get_neigh), c_loc(neighobject), ierr)
474  if (ierr /= 0) then
475  call my_error("set_callback_pointer", __line__, __file__)
476  end if
477 
478  call kim_model_get_influence_distance(model_handle, influence_distance)
479  call kim_model_get_number_of_neighbor_list_cutoffs(model_handle, &
480  number_of_cutoffs)
481  if (number_of_cutoffs /= 1) then
482  call my_error("too many cutoffs", __line__, __file__)
483  endif
484  call kim_model_get_neighbor_list_cutoffs(model_handle, cutoffs, ierr)
485  if (ierr /= 0) then
486  call my_error("get_cutoffs", __line__, __file__)
487  end if
488  cutoff = cutoffs(1)
489 
490  ! Setup cluster
491  !
492  do i=1,n
493  particle_species_codes(i) = species_code
494  enddo
495 
496  ! setup contributing particles
497  do i=1,n
498  particle_contributing(i) = 1 ! every particle contributes
499  enddo
500 
501  ! print header
502  print '(3A20)', "Energy", "Force Norm", "Lattice Spacing"
503  ! do the computations
504  current_spacing = min_spacing
505  do while (current_spacing < max_spacing)
506 
507  call create_fcc_configuration(current_spacing, ncellsperside, .false., &
508  coords, middledum)
509  ! Compute neighbor lists
510  call neigh_pure_cluster_neighborlist(.false., n, coords, (cutoff+cutpad), &
511  neighobject)
512 
513  ! Call model compute
514  call kim_model_compute(model_handle, compute_arguments_handle, ierr)
515  if (ierr /= 0) then
516  call my_error("kim_api_model_compute", __line__, __file__)
517  endif
518 
519  ! compue force_norm
520  force_norm = 0.0;
521  do i = 1, n
522  do j = 1, dim
523  force_norm = force_norm + forces(j,i)*forces(j,i)
524  end do
525  end do
526  force_norm = sqrt(force_norm)
527 
528  ! Print results to screen
529  !
530  print '(3ES20.10)', energy, force_norm, current_spacing
531 
532  current_spacing = current_spacing + spacing_incr
533  enddo
534 
535  ! Deallocate neighbor list object
536  deallocate( neighobject%neighborList )
537 
538  call kim_model_compute_arguments_destroy(&
539  model_handle, compute_arguments_handle, ierr)
540  if (ierr /= 0) then
541  call my_error("compute_arguments_destroy", __line__, __file__)
542  endif
543  call kim_model_destroy(model_handle)
544 
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)