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_cutoffs, 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_cutoffs
95  real(c_double), intent(in) :: cutoffs(number_of_cutoffs)
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_cutoffs /= 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  integer(c_int), parameter :: &
360  n = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
361  integer(c_int), parameter :: sizeone = 1
362 
363  type(neighobject_type), target :: neighobject
364 
365  !
366  ! KIM variables
367  !
368  character(len=KIM_KEY_STRING_LENGTH) :: testkimfile = "descriptor.kim"
369  character(len=KIM_KEY_STRING_LENGTH) :: modelname
370  character(len=KIM_KEY_STRING_LENGTH) :: nbc_method
371 
372  type(c_ptr) :: pkim
373  integer(c_int), pointer :: numberofparticles; type(c_ptr) :: pnparts
374  integer(c_int), pointer :: numberofspecies; type(c_ptr) :: pnofspecies
375  integer(c_int), pointer :: particlespecies(:); type(c_ptr) :: pparticlespecies
376  real(c_double), pointer :: cutoff; type(c_ptr) :: pcutoff
377  real(c_double), pointer :: energy; type(c_ptr) :: penergy
378  real(c_double), pointer :: coords(:,:); type(c_ptr) :: pcoor
379  real(c_double), pointer :: forces(:,:); type(c_ptr) :: pforces
380 
381 
382  integer(c_int) ierr
383  integer :: middledum
384  integer :: i, j
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  ! Initialize the KIM object
402  ierr = kim_api_file_init(pkim, testkimfile, modelname)
403  if (ierr.lt.kim_status_ok) then
404  call my_error("kim_api_file_init", __line__, __file__)
405  stop
406  endif
407 
408  ! make sure NEIGH_PURE_F
409  ierr = kim_api_get_nbc_method(pkim, nbc_method)
410  if (ierr.lt.kim_status_ok) then
411  call my_error("kim_api_get_nbc_method", __line__, __file__)
412  stop
413  endif
414  if (index(nbc_method,"NEIGH_PURE_F").ne.1) then
415  ierr = kim_status_fail
416  call my_error("Unknown NBC method", __line__, __file__)
417  stop
418  endif
419 
420  ! Allocate memory via the KIM system
421  call kim_api_allocate(pkim, n, aspecies, ierr)
422  if (ierr.lt.kim_status_ok) then
423  call my_error("kim_api_allocate", __line__, __file__)
424  stop
425  endif
426 
427  ! Allocate storage for neighbor lists
428  !
429  allocate(neighobject%neighborList(n+1,n))
430  neighobject%number_of_particles = n
431 
432  ! Set pointer in KIM object to neighbor list routine and object
433  !
434  ierr = kim_api_set_data(pkim, "neighObject", sizeone, c_loc(neighobject))
435  if (ierr.lt.kim_status_ok) then
436  call my_error("kim_api_set_data", __line__, __file__)
437  stop
438  endif
439 
440  ierr = kim_api_set_method(pkim, "get_neigh", sizeone, c_funloc(get_neigh))
441  if (ierr.lt.kim_status_ok) then
442  call my_error("kim_api_set_method", __line__, __file__)
443  stop
444  endif
445 
446  ! call model's init routine
447  ierr = kim_api_model_init(pkim)
448  if (ierr.lt.kim_status_ok) then
449  call my_error("kim_api_model_init", __line__, __file__)
450  stop
451  endif
452 
453 
454  ! Unpack data from KIM object
455  !
456  call kim_api_getm_data(pkim, ierr, &
457  "numberOfParticles", pnparts, 1, &
458  "numberOfSpecies", pnofspecies, 1, &
459  "particleSpecies", pparticlespecies, 1, &
460  "coordinates", pcoor, 1, &
461  "cutoff", pcutoff, 1, &
462  "energy", penergy, 1, &
463  "forces", pforces, 1)
464  if (ierr.lt.kim_status_ok) then
465  call my_error("kim_api_getm_data", __line__, __file__)
466  stop
467  endif
468  call c_f_pointer(pnparts, numberofparticles)
469  call c_f_pointer(pnofspecies, numberofspecies)
470  call c_f_pointer(pparticlespecies, particlespecies, [n])
471  call c_f_pointer(pcoor, coords, [dim,n])
472  call c_f_pointer(pcutoff, cutoff)
473  call c_f_pointer(penergy, energy)
474  call c_f_pointer(pforces, forces, [dim,n])
475 
476 
477  ! Set values
478  numberofparticles = n
479  numberofspecies = aspecies
480  particlespecies(:) = kim_api_get_species_code(pkim, "Ar", ierr)
481  if (ierr.lt.kim_status_ok) then
482  call my_error("kim_api_get_species_code", __line__, __file__)
483  stop
484  endif
485 
486  ! print header
487  print '(3A20)', "Energy", "Force Norm", "Lattice Spacing"
488  ! do the computations
489  current_spacing = min_spacing
490  do while (current_spacing < max_spacing)
491 
492  call create_fcc_configuration(current_spacing, ncellsperside, .false., &
493  coords, middledum)
494  ! Compute neighbor lists
495  call neigh_pure_cluster_neighborlist(.false., n, coords, (cutoff+cutpad), &
496  neighobject)
497 
498  ! Call model compute
499  ierr = kim_api_model_compute(pkim)
500  if (ierr.lt.kim_status_ok) then
501  call my_error("kim_api_model_compute", __line__, __file__)
502  stop
503  endif
504 
505  ! compue force_norm
506  force_norm = 0.0;
507  do i = 1, n
508  do j = 1, dim
509  force_norm = force_norm + forces(j,i)*forces(j,i)
510  end do
511  end do
512  force_norm = sqrt(force_norm)
513 
514  ! Print results to screen
515  !
516  print '(3ES20.10)', energy, force_norm, current_spacing
517 
518  current_spacing = current_spacing + spacing_incr
519  enddo
520 
521  ! Deallocate neighbor list object
522  deallocate( neighobject%neighborList )
523 
524  ierr = kim_api_model_destroy(pkim)
525  if (ierr.lt.kim_status_ok) then
526  call my_error("kim_api_model_destroy", __line__, __file__)
527  stop
528  endif
529  call kim_api_free(pkim, ierr)
530  if (ierr.lt.kim_status_ok) then
531  call my_error("kim_api_free", __line__, __file__)
532  stop
533  endif
534 
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)