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