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