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=*), intent(in) :: message
39  integer, intent(in) :: line
40  character(len=*), 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=*), intent(in) :: message
50  integer, intent(in) :: line
51  character(len=*), 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 
68 
69  use, intrinsic :: iso_c_binding
70 
71  public get_neigh
72 
74  integer(c_int) :: number_of_particles
75  integer(c_int), pointer :: neighborlist(:,:)
76  real(c_double), pointer :: rijlist(:,:,:)
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, neighbor_list_index, request, numnei, &
88  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) :: neighbor_list_index
95  integer(c_int), value, intent(in) :: request
96  integer(c_int), intent(out) :: numnei
97  type(c_ptr), intent(out) :: pnei1part
98  integer(c_int), intent(out) :: ierr
99 
100  !-- Local variables
101  integer(c_int), parameter :: dim = 3
102  integer(c_int) numberofparticles
103  type(neighobject_type), pointer :: neighobject
104 
105  if (neighbor_list_index /= 1) then
106  call my_warning("wrong list index", __line__, __file__)
107  ierr = 1
108  return
109  endif
110 
111  call c_f_pointer(data_object, neighobject)
112 
113  numberofparticles = neighobject%number_of_particles
114 
115  if ( (request.gt.numberofparticles) .or. (request.lt.1)) then
116  call my_warning("Invalid part ID in get_neigh", &
117  __line__, __file__)
118  ierr = 1
119  return
120  endif
121 
122  ! set the returned number of neighbors for the returned part
123  numnei = neighobject%neighborList(1,request)
124 
125  ! set the location for the returned neighbor list
126  pnei1part = c_loc(neighobject%neighborList(2,request))
127 
128  ierr = 0
129  return
130 end subroutine get_neigh
131 
132 end module mod_neighborlist
133 
134 !-------------------------------------------------------------------------------
135 !
136 ! NEIGH_PURE_cluster_neighborlist
137 !
138 !-------------------------------------------------------------------------------
139 subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, &
140  cutoff, neighObject)
141  use, intrinsic :: iso_c_binding
142  use mod_neighborlist
143  implicit none
144 
145  !-- Transferred variables
146  logical, intent(in) :: half
147  integer(c_int), intent(in) :: numberOfParticles
148  real(c_double), dimension(3,numberOfParticles), &
149  intent(in) :: coords
150  real(c_double), intent(in) :: cutoff
151  type(neighObject_type), intent(inout) :: neighObject
152 
153  !-- Local variables
154  integer(c_int) i, j, a
155  real(c_double) dx(3)
156  real(c_double) r2
157  real(c_double) cutoff2
158 
159  cutoff2 = cutoff**2
160 
161  do i=1,numberofparticles
162  a = 1
163  do j=1,numberofparticles
164  dx(:) = coords(:, j) - coords(:, i)
165  r2 = dot_product(dx, dx)
166  if (r2.le.cutoff2) then
167  ! part j is a neighbor of part i
168  if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) ) then
169  a = a+1
170  neighobject%neighborList(a,i) = j
171  endif
172  endif
173  enddo
174  ! part i has a-1 neighbors
175  neighobject%neighborList(1,i) = a-1
176  enddo
177 
178  return
179 
180 end subroutine neigh_pure_cluster_neighborlist
181 
182 
183 !-------------------------------------------------------------------------------
184 !
185 ! create_FCC_configuration subroutine
186 !
187 ! creates a cubic configuration of FCC particles with lattice spacing
188 ! `FCCspacing' and `nCellsPerSide' cells along each direction.
189 !
190 ! With periodic==.true. this will result in a total number of particles equal
191 ! to 4*(nCellsPerSide)**3 + 6*(nCellsPerSide)**2 + 3*(nCellsPerSide) + 1
192 !
193 ! With periodic==.false. this will result in a total number of particles equal
194 ! to 4*(nCellsPerSide)**3
195 !
196 ! Returns the Id of the particle situated in the middle of the configuration
197 ! (this particle will have the most neighbors.)
198 !
199 !-------------------------------------------------------------------------------
200 subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, &
201  coords, MiddlePartId)
202  use, intrinsic :: iso_c_binding
203  implicit none
204  integer(c_int), parameter :: cd = c_double ! used for literal constants
205 
206  !-- Transferred variables
207  real(c_double), intent(in) :: FCCspacing
208  integer(c_int), intent(in) :: nCellsPerSide
209  logical, intent(in) :: periodic
210  real(c_double), intent(out) :: coords(3,*)
211  integer(c_int), intent(out) :: MiddlePartId
212  !
213  ! cluster setup variables
214  !
215  real(c_double) FCCshifts(3,4)
216  real(c_double) latVec(3)
217  integer(c_int) a, i, j, k, m
218 
219  ! Create a cubic FCC cluster
220  !
221  fccshifts(1,1) = 0.0_cd
222  fccshifts(2,1) = 0.0_cd
223  fccshifts(3,1) = 0.0_cd
224  fccshifts(1,2) = 0.5_cd*fccspacing
225  fccshifts(2,2) = 0.5_cd*fccspacing
226  fccshifts(3,2) = 0.0_cd
227  fccshifts(1,3) = 0.5_cd*fccspacing
228  fccshifts(2,3) = 0.0_cd
229  fccshifts(3,3) = 0.5_cd*fccspacing
230  fccshifts(1,4) = 0.0_cd
231  fccshifts(2,4) = 0.5_cd*fccspacing
232  fccshifts(3,4) = 0.5_cd*fccspacing
233 
234  middlepartid = 1 ! Always put middle particle as #1
235  a = 1 ! leave space for middle particle as particle #1
236  do i=1,ncellsperside
237  latvec(1) = (i-1)*fccspacing
238  do j=1,ncellsperside
239  latvec(2) = (j-1)*fccspacing
240  do k=1,ncellsperside
241  latvec(3) = (k-1)*fccspacing
242  do m=1,4
243  a = a+1
244  coords(:,a) = latvec + fccshifts(:,m)
245  if ((i.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
246  (k.eq.ncellsperside/2+1) .and. (m.eq.1)) then
247  coords(:,1) = latvec + fccshifts(:,m) ! put middle particle as #1
248  a = a - 1
249  endif
250  enddo
251  enddo
252  if (.not. periodic) then
253  ! Add in the remaining three faces
254  ! pos-x face
255  latvec(1) = ncellsperside*fccspacing
256  latvec(2) = (i-1)*fccspacing
257  latvec(3) = (j-1)*fccspacing
258  a = a+1; coords(:,a) = latvec
259  a = a+1; coords(:,a) = latvec + fccshifts(:,4)
260  ! pos-y face
261  latvec(1) = (i-1)*fccspacing
262  latvec(2) = ncellsperside*fccspacing
263  latvec(3) = (j-1)*fccspacing
264  a = a+1; coords(:,a) = latvec
265  a = a+1; coords(:,a) = latvec + fccshifts(:,3)
266  ! pos-z face
267  latvec(1) = (i-1)*fccspacing
268  latvec(2) = (j-1)*fccspacing
269  latvec(3) = ncellsperside*fccspacing
270  a = a+1; coords(:,a) = latvec
271  a = a+1; coords(:,a) = latvec + fccshifts(:,2)
272  endif
273  enddo
274  if (.not. periodic) then
275  ! Add in the remaining three edges
276  latvec(1) = (i-1)*fccspacing
277  latvec(2) = ncellsperside*fccspacing
278  latvec(3) = ncellsperside*fccspacing
279  a = a+1; coords(:,a) = latvec
280  latvec(1) = ncellsperside*fccspacing
281  latvec(2) = (i-1)*fccspacing
282  latvec(3) = ncellsperside*fccspacing
283  a = a+1; coords(:,a) = latvec
284  latvec(1) = ncellsperside*fccspacing
285  latvec(2) = ncellsperside*fccspacing
286  latvec(3) = (i-1)*fccspacing
287  a = a+1; coords(:,a) = latvec
288  endif
289  enddo
290  if (.not. periodic) then
291  ! Add in the remaining corner
292  a = a+1; coords(:,a) = ncellsperside*fccspacing
293  endif
294 
295  return
296 
297 end subroutine create_fcc_configuration
298 
299 
300 
301 !*******************************************************************************
302 !**
303 !** PROGRAM vc_forces_numer_deriv
304 !**
305 !** KIM compliant program to perform numerical derivative check on a model
306 !**
307 !*******************************************************************************
308 
309 !-------------------------------------------------------------------------------
310 !
311 ! Main program
312 !
313 !-------------------------------------------------------------------------------
315  use, intrinsic :: iso_c_binding
316  use error
320  use kim_model_module
324  use mod_neighborlist
325  implicit none
326  integer(c_int), parameter :: cd = c_double ! used for literal constants
327 
328  integer(c_int), parameter :: ncellsperside = 2
329  integer(c_int), parameter :: dim = 3
330 
331  real(c_double), parameter :: cutpad = 0.75_cd
332  real(c_double), parameter :: fccspacing = 5.260_cd
333  real(c_double), parameter :: min_spacing = 0.8*fccspacing
334  real(c_double), parameter :: max_spacing = 1.2*fccspacing
335  real(c_double), parameter :: spacing_incr = 0.025*fccspacing
336  real(c_double) :: current_spacing
337 
338  character(len=256) :: modelname
339 
340  integer(c_int), parameter :: &
341  n = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
342 
343  type(neighobject_type), target :: neighobject
344 
345  type(kim_model_handle_type) :: model_handle
346  real(c_double) :: influence_distance
347  integer(c_int) :: number_of_cutoffs
348  real(c_double) :: cutoff
349  real(c_double) :: cutoffs(1)
350  integer(c_int), target :: particle_species_codes(n)
351  integer(c_int), target :: particle_contributing(n)
352  real(c_double), target :: energy
353  real(c_double), target :: coords(dim, n)
354  integer(c_int) i, ierr, ierr2
355 
356  integer(c_int) species_is_supported
357  integer(c_int) species_code
358  integer(c_int) requested_units_accepted
359  real(c_double), pointer :: null_pointer
360 
361  integer :: middledum
362 
363  nullify(null_pointer)
364 
365  ! Initialize error flag
366  ierr = 0
367 
368  ! Get KIM Model name to use
369  print '("Please enter a valid KIM model name: ")'
370  read(*,*) modelname
371 
372  ! Print output header
373  !
374  print *
375  print *,'This is Test : ex_test_Ar_fcc_cluster.'
376  print *
377  print '(120(''-''))'
378  print '("Results for KIM Model : ",A)', trim(modelname)
379 
380  ! Create empty KIM object
381  !
383  kim_length_unit_a, &
384  kim_energy_unit_ev, &
385  kim_charge_unit_e, &
386  kim_temperature_unit_k, &
387  kim_time_unit_ps, &
388  trim(modelname), &
389  requested_units_accepted, &
390  model_handle, ierr)
391  if (ierr /= 0) then
392  call my_error("kim_api_create", __line__, __file__)
393  endif
394 
395  ! check that we are compatible
396  if (requested_units_accepted == 0) then
397  call my_error("Must adapt to model units", __line__, __file__)
398  end if
399 
400  ! check that model supports Ar
401  !
402  call kim_model_get_species_support_and_code(model_handle, &
403  kim_species_name_ar, species_is_supported, species_code, ierr)
404  if ((ierr /= 0) .or. (species_is_supported /= 1)) then
405  call my_error("Model does not support Ar", __line__, __file__)
406  endif
407 
408  ! it would be good to check that the model is compatible
409  ! but we will skip it here
410 
411  ! register memory with the KIM system
412  ierr = 0
413  call kim_model_set_argument_pointer(model_handle, &
415  ierr = ierr + ierr2
416  call kim_model_set_argument_pointer(model_handle, &
417  kim_argument_name_particle_species_codes, particle_species_codes, ierr2)
418  ierr = ierr + ierr2
419  call kim_model_set_argument_pointer(model_handle, &
420  kim_argument_name_particle_contributing, particle_contributing, ierr2)
421  ierr = ierr + ierr2
422  call kim_model_set_argument_pointer(model_handle, &
423  kim_argument_name_coordinates, coords, ierr2)
424  ierr = ierr + ierr2
425  call kim_model_set_argument_pointer(model_handle, &
426  kim_argument_name_partial_energy, energy, ierr2)
427  ierr = ierr + ierr2
428  if (ierr /= 0) then
429  call my_error("set_argument_pointer", __line__, __file__)
430  endif
431 
432  ! Set pointer in KIM object to neighbor list routine and object
433  !
434  call kim_model_set_callback_pointer(model_handle, &
436  c_funloc(get_neigh), c_loc(neighobject), ierr)
437  if (ierr /= 0) then
438  call my_error("set_callback_pointer", __line__, __file__)
439  end if
440 
441  call kim_model_get_influence_distance(model_handle, influence_distance)
442  call kim_model_get_number_of_neighbor_list_cutoffs(model_handle, &
443  number_of_cutoffs)
444  if (number_of_cutoffs /= 1) then
445  call my_error("too many cutoffs", __line__, __file__)
446  endif
447  call kim_model_get_neighbor_list_cutoffs(model_handle, cutoffs, ierr)
448  if (ierr /= 0) then
449  call my_error("get_cutoffs", __line__, __file__)
450  end if
451  cutoff = cutoffs(1)
452 
453  ! Setup cluster
454  !
455  do i=1,n
456  particle_species_codes(i) = species_code
457  enddo
458 
459  ! setup contributing particles
460  do i=1,n
461  particle_contributing(i) = 1 ! every particle contributes
462  enddo
463 
464  ! Allocate storage for neighbor lists and
465  ! store pointers to neighbor list object and access function
466  !
467  allocate(neighobject%neighborList(n+1,n))
468  neighobject%number_of_particles = n
469 
470  ! do the computations
471  current_spacing = min_spacing
472  do while (current_spacing < max_spacing)
473 
474  call create_fcc_configuration(current_spacing, ncellsperside, .false., &
475  coords, middledum)
476  ! Compute neighbor lists
477  !
478  call neigh_pure_cluster_neighborlist(.false., n, coords, cutoff+cutpad, &
479  neighobject)
480  if (ierr /= 0) then
481  call my_error("update_neighborlist", __line__, __file__)
482  endif
483 
484  ! Call model compute to get forces (gradient)
485  !
486  call kim_model_compute(model_handle, ierr)
487  if (ierr /= 0) then
488  call my_error("kim_api_model_compute", __line__, __file__)
489  endif
490 
491  ! Print results to screen
492  !
493  print '(2ES20.10)', energy, current_spacing
494 
495  current_spacing = current_spacing + spacing_incr
496  enddo
497 
498  call kim_model_destroy(model_handle)
499 
500 end program ex_test_ar_fcc_cluster
type(kim_numbering_type), public, protected kim_numbering_one_based
type(kim_argument_name_type), public, protected kim_argument_name_number_of_particles
type(kim_language_name_type), public, protected kim_language_name_fortran
program ex_test_ar_fcc_cluster
subroutine my_error(message, line, file)
type(kim_argument_name_type), public, protected kim_argument_name_particle_contributing
type(kim_callback_name_type), public, protected kim_callback_name_get_neighbor_list
type(kim_argument_name_type), public, protected kim_argument_name_particle_species_codes
subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, cutoff, neighObject)
type(kim_argument_name_type), public, protected kim_argument_name_partial_energy
subroutine my_warning(message, line, file)
type(kim_argument_name_type), public, protected kim_argument_name_coordinates
subroutine, public get_neigh(data_object, neighbor_list_index, request, numnei, pnei1part, ierr)
subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, coords, MiddlePartId)
type(kim_species_name_type), public, protected kim_species_name_ar