kim-api 2.3.0+AppleClang.AppleClang.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
Loading...
Searching...
No Matches
ex_test_Ar_fcc_cluster_fortran.f90
Go to the documentation of this file.
1!
2! KIM-API: An API for interatomic models
3! Copyright (c) 2013--2022, Regents of the University of Minnesota.
4! All rights reserved.
5!
6! Contributors:
7! Ellad B. Tadmor
8! Ryan S. Elliott
9! Stephen M. Whalen
10!
11! SPDX-License-Identifier: LGPL-2.1-or-later
12!
13! This library is free software; you can redistribute it and/or
14! modify it under the terms of the GNU Lesser General Public
15! License as published by the Free Software Foundation; either
16! version 2.1 of the License, or (at your option) any later version.
17!
18! This library is distributed in the hope that it will be useful,
19! but WITHOUT ANY WARRANTY; without even the implied warranty of
20! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
21! Lesser General Public License for more details.
22!
23! You should have received a copy of the GNU Lesser General Public License
24! along with this library; if not, write to the Free Software Foundation,
25! Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
26!
27
28module error
29 use, intrinsic :: iso_c_binding
30 implicit none
31
32 public
33
34contains
35 recursive subroutine my_error(message)
36 implicit none
37 character(len=*, kind=c_char), intent(in) :: message
38
39 print *, "* Error : ", trim(message)
40 stop 1
41 end subroutine my_error
42
43 recursive subroutine my_warning(message)
44 implicit none
45 character(len=*, kind=c_char), intent(in) :: message
46
47 print *, "* Warning : ", trim(message)
48 end subroutine my_warning
49end module error
50
51!-------------------------------------------------------------------------------
52!
53! module mod_neighborlist :
54!
55! Module contains type and routines related to neighbor list calculation
56!
57!-------------------------------------------------------------------------------
58
60
61 use, intrinsic :: iso_c_binding
62
63 public get_neigh
64
65 type, bind(c) :: neighObject_type
66 real(c_double) :: cutoff
67 integer(c_int) :: number_of_particles
68 type(c_ptr) :: neighbor_list_pointer
69 end type neighobject_type
70contains
71
72 !-----------------------------------------------------------------------------
73 !
74 ! get_neigh neighbor list access function
75 !
76 ! This function implements Locator and Iterator mode
77 !
78 !-----------------------------------------------------------------------------
79 recursive subroutine get_neigh(data_object, number_of_neighbor_lists, &
80 cutoffs, neighbor_list_index, request, &
81 numnei, pnei1part, ierr) bind(c)
82 use error
83 implicit none
84
85 !-- Transferred variables
86 type(c_ptr), value, intent(in) :: data_object
87 integer(c_int), value, intent(in) :: number_of_neighbor_lists
88 real(c_double), intent(in) :: cutoffs(*)
89 integer(c_int), value, intent(in) :: neighbor_list_index
90 integer(c_int), value, intent(in) :: request
91 integer(c_int), intent(out) :: numnei
92 type(c_ptr), intent(out) :: pnei1part
93 integer(c_int), intent(out) :: ierr
94
95 !-- Local variables
96 integer(c_int) numberofparticles
97 type(neighobject_type), pointer :: neighobject
98 integer(c_int), pointer :: neighborlist(:, :)
99
100 call c_f_pointer(data_object, neighobject)
101 numberofparticles = neighobject%number_of_particles
102 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
103 [numberofparticles + 1, numberofparticles])
104
105 if (number_of_neighbor_lists /= 1) then
106 call my_warning("invalid number of cutoffs")
107 ierr = 1
108 return
109 end if
110
111 if (cutoffs(1) > neighobject%cutoff) then
112 call my_warning("neighbor list cutoff too small for model cutoff")
113 ierr = 1
114 return
115 end if
116
117 if (neighbor_list_index /= 1) then
118 call my_warning("wrong list index")
119 ierr = 1
120 return
121 end if
122
123 if ((request > numberofparticles) .or. (request < 1)) then
124 print *, request
125 call my_warning("Invalid part ID in get_neigh")
126 ierr = 1
127 return
128 end if
129
130 ! set the returned number of neighbors for the returned part
131 numnei = neighborlist(1, request)
132
133 ! set the location for the returned neighbor list
134 pnei1part = c_loc(neighborlist(2, request))
135
136 ierr = 0
137 return
138 end subroutine get_neigh
139
140end module mod_neighborlist
141
143 implicit none
144 public
145
146contains
147
148 !-----------------------------------------------------------------------------
149 !
150 ! NEIGH_PURE_cluster_neighborlist
151 !
152 !-----------------------------------------------------------------------------
153 recursive subroutine neigh_pure_cluster_neighborlist( &
154 half, numberOfParticles, coords, cutoff, neighObject)
155 use, intrinsic :: iso_c_binding
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 integer(c_int), pointer :: neighborlist(:, :)
173
174 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
175 [numberofparticles + 1, numberofparticles])
176
177 neighobject%cutoff = cutoff
178
179 cutoff2 = cutoff**2
180
181 do i = 1, numberofparticles
182 a = 1
183 do j = 1, numberofparticles
184 dx(:) = coords(:, j) - coords(:, i)
185 r2 = dot_product(dx, dx)
186 if (r2 <= cutoff2) then
187 ! part j is a neighbor of part i
188 if ((j > i) .OR. ((.not. half) .AND. (i /= j))) then
189 a = a + 1
190 neighborlist(a, i) = j
191 end if
192 end if
193 end do
194 ! part i has a-1 neighbors
195 neighborlist(1, i) = a - 1
196 end do
197
198 return
199
201
202 !-----------------------------------------------------------------------------
203 !
204 ! create_FCC_configuration subroutine
205 !
206 ! creates a cubic configuration of FCC particles with lattice spacing
207 ! `FCCspacing' and `nCellsPerSide' cells along each direction.
208 !
209 ! With periodic==.true. this will result in a total number of particles equal
210 ! to 4*(nCellsPerSide)**3 + 6*(nCellsPerSide)**2 + 3*(nCellsPerSide) + 1
211 !
212 ! With periodic==.false. this will result in a total number of particles equal
213 ! to 4*(nCellsPerSide)**3
214 !
215 ! Returns the Id of the particle situated in the middle of the configuration
216 ! (this particle will have the most neighbors.)
217 !
218 !-----------------------------------------------------------------------------
219 recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, &
220 periodic, coords, MiddlePartId)
221 use, intrinsic :: iso_c_binding
222 implicit none
223 integer(c_int), parameter :: cd = c_double ! used for literal constants
224
225 !-- Transferred variables
226 real(c_double), intent(in) :: fccspacing
227 integer(c_int), intent(in) :: ncellsperside
228 logical, intent(in) :: periodic
229 real(c_double), intent(out) :: coords(3, *)
230 integer(c_int), intent(out) :: middlepartid
231 !
232 ! cluster setup variables
233 !
234 real(c_double) fccshifts(3, 4)
235 real(c_double) latvec(3)
236 integer(c_int) a, i, j, k, m
237
238 ! Create a cubic FCC cluster
239 !
240 fccshifts(1, 1) = 0.0_cd
241 fccshifts(2, 1) = 0.0_cd
242 fccshifts(3, 1) = 0.0_cd
243 fccshifts(1, 2) = 0.5_cd * fccspacing
244 fccshifts(2, 2) = 0.5_cd * fccspacing
245 fccshifts(3, 2) = 0.0_cd
246 fccshifts(1, 3) = 0.5_cd * fccspacing
247 fccshifts(2, 3) = 0.0_cd
248 fccshifts(3, 3) = 0.5_cd * fccspacing
249 fccshifts(1, 4) = 0.0_cd
250 fccshifts(2, 4) = 0.5_cd * fccspacing
251 fccshifts(3, 4) = 0.5_cd * fccspacing
252
253 middlepartid = 1 ! Always put middle particle as #1
254 a = 1 ! leave space for middle particle as particle #1
255 do i = 1, ncellsperside
256 latvec(1) = (i - 1) * fccspacing
257 do j = 1, ncellsperside
258 latvec(2) = (j - 1) * fccspacing
259 do k = 1, ncellsperside
260 latvec(3) = (k - 1) * fccspacing
261 do m = 1, 4
262 a = a + 1
263 coords(:, a) = latvec + fccshifts(:, m)
264 if ((i == ncellsperside / 2 + 1) &
265 .and. (j == ncellsperside / 2 + 1) &
266 .and. (k == ncellsperside / 2 + 1) &
267 .and. (m == 1)) &
268 then
269 ! put middle particle as #1
270 coords(:, 1) = latvec + fccshifts(:, m)
271 a = a - 1
272 end if
273 end do
274 end do
275 if (.not. periodic) then
276 ! Add in the remaining three faces
277 ! pos-x face
278 latvec(1) = ncellsperside * fccspacing
279 latvec(2) = (i - 1) * fccspacing
280 latvec(3) = (j - 1) * fccspacing
281 a = a + 1; coords(:, a) = latvec
282 a = a + 1; coords(:, a) = latvec + fccshifts(:, 4)
283 ! pos-y face
284 latvec(1) = (i - 1) * fccspacing
285 latvec(2) = ncellsperside * fccspacing
286 latvec(3) = (j - 1) * fccspacing
287 a = a + 1; coords(:, a) = latvec
288 a = a + 1; coords(:, a) = latvec + fccshifts(:, 3)
289 ! pos-z face
290 latvec(1) = (i - 1) * fccspacing
291 latvec(2) = (j - 1) * fccspacing
292 latvec(3) = ncellsperside * fccspacing
293 a = a + 1; coords(:, a) = latvec
294 a = a + 1; coords(:, a) = latvec + fccshifts(:, 2)
295 end if
296 end do
297 if (.not. periodic) then
298 ! Add in the remaining three edges
299 latvec(1) = (i - 1) * fccspacing
300 latvec(2) = ncellsperside * fccspacing
301 latvec(3) = ncellsperside * fccspacing
302 a = a + 1; coords(:, a) = latvec
303 latvec(1) = ncellsperside * fccspacing
304 latvec(2) = (i - 1) * fccspacing
305 latvec(3) = ncellsperside * fccspacing
306 a = a + 1; coords(:, a) = latvec
307 latvec(1) = ncellsperside * fccspacing
308 latvec(2) = ncellsperside * fccspacing
309 latvec(3) = (i - 1) * fccspacing
310 a = a + 1; coords(:, a) = latvec
311 end if
312 end do
313 if (.not. periodic) then
314 ! Add in the remaining corner
315 a = a + 1; coords(:, a) = ncellsperside * fccspacing
316 end if
317
318 return
319
320 end subroutine create_fcc_configuration
321
322end module mod_utility
323
324!*******************************************************************************
325!**
326!** PROGRAM vc_forces_numer_deriv
327!**
328!** KIM compliant program to perform numerical derivative check on a model
329!**
330!*******************************************************************************
331
332!-------------------------------------------------------------------------------
333!
334! Main program
335!
336!-------------------------------------------------------------------------------
338 use, intrinsic :: iso_c_binding
339 use error
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
350 real(c_double), parameter :: cutpad = 0.75_cd
351 real(c_double), parameter :: fccspacing = 5.260_cd
352 real(c_double), parameter :: min_spacing = 0.8 * fccspacing
353 real(c_double), parameter :: max_spacing = 1.2 * fccspacing
354 real(c_double), parameter :: spacing_incr = 0.025 * fccspacing
355 real(c_double) :: current_spacing
356 real(c_double) :: force_norm
357
358 character(len=256, kind=c_char) :: modelname
359
360 integer(c_int), parameter :: n = 4 * (ncellsperside)**3 &
361 + 6 * (ncellsperside)**2 &
362 + 3 * (ncellsperside) + 1
363
364 type(neighobject_type), target :: neighobject
365 integer(c_int), allocatable, target :: neighborlist(:, :)
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) :: &
374 model_will_not_request_neighbors_of_noncontributing_particles(1)
375 integer(c_int), target :: particle_species_codes(n)
376 integer(c_int), target :: particle_contributing(n)
377 real(c_double), target :: energy
378 real(c_double), target :: coords(dim, n)
379 real(c_double), target :: forces(dim, n)
380 integer(c_int) i, j, ierr, ierr2
381
382 integer(c_int) species_is_supported
383 integer(c_int) species_code
384 integer(c_int) requested_units_accepted
385 integer(c_int) number_of_model_routine_names
386 type(kim_model_routine_name_type) model_routine_name
387 integer(c_int) present
388 integer(c_int) required
389 type(kim_supported_extensions_type), target :: supported_extensions
390 character(len=KIM_MAX_EXTENSION_ID_LENGTH, kind=c_char) :: id_string
391
392 integer :: middledum
393
394 ! Initialize
395 ierr = 0
396 supported_extensions%number_of_supported_extensions = 0
397
398 ! Get KIM Model name to use
399 print '("Please enter a valid KIM model name: ")'
400 read (*, *) modelname
401
402 ! Create empty KIM object
403 !
404 call kim_model_create(kim_numbering_one_based, &
405 kim_length_unit_a, &
406 kim_energy_unit_ev, &
407 kim_charge_unit_e, &
408 kim_temperature_unit_k, &
409 kim_time_unit_ps, &
410 trim(modelname), &
411 requested_units_accepted, &
412 model_handle, ierr)
413 if (ierr /= 0) then
414 call my_error("kim_api_create")
415 end if
416
417 ! check that we are compatible
418 if (requested_units_accepted == 0) then
419 call my_error("Must adapt to model units")
420 end if
421
422 ! check that we know about all required routines
423 call kim_get_number_of_model_routine_names(number_of_model_routine_names)
424 do i = 1, number_of_model_routine_names
425 call kim_get_model_routine_name(i, model_routine_name, ierr)
426 if (ierr /= 0) then
427 call my_error("kim_get_model_routine_name")
428 end if
429 call kim_is_routine_present(model_handle, model_routine_name, present, &
430 required, ierr)
431 if (ierr /= 0) then
432 call my_error("kim_is_routine_present")
433 end if
434
435 if ((present == 1) .and. (required == 1)) then
436 if (.not. ((model_routine_name == kim_model_routine_name_create) &
437 .or. (model_routine_name == &
438 kim_model_routine_name_compute_arguments_create) &
439 .or. (model_routine_name == kim_model_routine_name_compute) &
440 .or. (model_routine_name == kim_model_routine_name_refresh) &
441 .or. (model_routine_name == &
442 kim_model_routine_name_compute_arguments_destroy) &
443 .or. (model_routine_name == kim_model_routine_name_destroy))) &
444 then
445 call my_error("Unknown required ModelRoutineName found.")
446 end if
447 end if
448 end do
449
450 ! check that model supports Ar
451 !
452 call kim_get_species_support_and_code( &
453 model_handle, kim_species_name_ar, species_is_supported, species_code, ierr)
454 if ((ierr /= 0) .or. (species_is_supported /= 1)) then
455 call my_error("Model does not support Ar")
456 end if
457
458 ! Check supported extensions, if any
459 call kim_is_routine_present( &
460 model_handle, kim_model_routine_name_extension, present, required, ierr)
461 if (ierr /= 0) then
462 call my_error("Unable to get Extension present/required.")
463 end if
464 if (present /= 0) then
465 call kim_extension(model_handle, kim_supported_extensions_id, &
466 c_loc(supported_extensions), ierr)
467 if (ierr /= 0) then
468 call my_error("Error returned from kim_model_extension().")
469 end if
470 write (*, '(A,I2,A)') "Model Supports ", &
471 supported_extensions%number_of_supported_extensions, &
472 " Extensions:"
473 do i = 1, supported_extensions%number_of_supported_extensions
474 call kim_c_char_array_to_string( &
475 supported_extensions%supported_extension_id(:, i), id_string)
476 write (*, '(A,I2,A,A,A,A,I2)') " supportedExtensionID[", i, '] = "', &
477 trim(id_string), '" ', &
478 "which has required = ", &
479 supported_extensions%supported_extension_required(i)
480 end do
481 end if
482
483 ! Best-practice is to check that the model is compatible
484 ! but we will skip it here
485
486 ! create compute_arguments object
487 call kim_compute_arguments_create( &
488 model_handle, compute_arguments_handle, ierr)
489 if (ierr /= 0) then
490 call my_error("kim_model_compute_arguments_create")
491 end if
492
493 ! register memory with the KIM system
494 ierr = 0
495 call kim_set_argument_pointer( &
496 compute_arguments_handle, kim_compute_argument_name_number_of_particles, &
497 n, ierr2)
498 ierr = ierr + ierr2
499 call kim_set_argument_pointer( &
500 compute_arguments_handle, &
501 kim_compute_argument_name_particle_species_codes, &
502 particle_species_codes, ierr2)
503 ierr = ierr + ierr2
504 call kim_set_argument_pointer( &
505 compute_arguments_handle, &
506 kim_compute_argument_name_particle_contributing, particle_contributing, &
507 ierr2)
508 ierr = ierr + ierr2
509 call kim_set_argument_pointer( &
510 compute_arguments_handle, kim_compute_argument_name_coordinates, coords, &
511 ierr2)
512 ierr = ierr + ierr2
513 call kim_set_argument_pointer( &
514 compute_arguments_handle, kim_compute_argument_name_partial_energy, &
515 energy, ierr2)
516 ierr = ierr + ierr2
517 call kim_set_argument_pointer( &
518 compute_arguments_handle, kim_compute_argument_name_partial_forces, &
519 forces, ierr2)
520 ierr = ierr + ierr2
521 if (ierr /= 0) then
522 call my_error("set_argument_pointer")
523 end if
524
525 ! Allocate storage for neighbor lists
526 !
527 allocate (neighborlist(n + 1, n))
528 neighobject%neighbor_list_pointer = c_loc(neighborlist)
529 neighobject%number_of_particles = n
530
531 ! Set pointer in KIM object to neighbor list routine and object
532 !
533 call kim_set_callback_pointer( &
534 compute_arguments_handle, kim_compute_callback_name_get_neighbor_list, &
535 kim_language_name_fortran, c_funloc(get_neigh), c_loc(neighobject), ierr)
536 if (ierr /= 0) then
537 call my_error("set_callback_pointer")
538 end if
539
540 call kim_get_influence_distance(model_handle, influence_distance)
541 call kim_get_number_of_neighbor_lists(model_handle, &
542 number_of_neighbor_lists)
543 if (number_of_neighbor_lists /= 1) then
544 call my_error("too many neighbor lists")
545 end if
546 call kim_get_neighbor_list_values( &
547 model_handle, cutoffs, &
548 model_will_not_request_neighbors_of_noncontributing_particles, ierr)
549 if (ierr /= 0) then
550 call my_error("get_neighbor_list_values")
551 end if
552 cutoff = cutoffs(1)
553
554 ! Setup cluster
555 !
556 do i = 1, n
557 particle_species_codes(i) = species_code
558 end do
559
560 ! setup contributing particles
561 do i = 1, n
562 particle_contributing(i) = 1 ! every particle contributes
563 end do
564
565 ! Print output header
566 !
567 print *
568 print *, 'This is Test : ex_test_Ar_fcc_cluster_fortran'
569 print *
570 print '(80(''-''))'
571 print '("Results for KIM Model : ",A)', trim(modelname)
572
573 ! print header
574 print '(3A20)', "Energy", "Force Norm", "Lattice Spacing"
575 ! do the computations
576 current_spacing = min_spacing
577 do while (current_spacing < max_spacing)
578
579 call create_fcc_configuration(current_spacing, ncellsperside, .false., &
580 coords, middledum)
581 ! Compute neighbor lists
582 call neigh_pure_cluster_neighborlist(.false., n, coords, &
583 (cutoff + cutpad), neighobject)
584
585 ! Call model compute
586 call kim_compute(model_handle, compute_arguments_handle, ierr)
587 if (ierr /= 0) then
588 call my_error("kim_api_model_compute")
589 end if
590
591 ! compue force_norm
592 force_norm = 0.0
593 do i = 1, n
594 do j = 1, dim
595 force_norm = force_norm + forces(j, i) * forces(j, i)
596 end do
597 end do
598 force_norm = sqrt(force_norm)
599
600 ! Print results to screen
601 !
602 print '(3ES20.10)', energy, force_norm, current_spacing
603
604 current_spacing = current_spacing + spacing_incr
605 end do
606
607 ! Deallocate neighbor list object
608 deallocate (neighborlist)
609
610 call kim_compute_arguments_destroy( &
611 model_handle, compute_arguments_handle, ierr)
612 if (ierr /= 0) then
613 call my_error("compute_arguments_destroy")
614 end if
615 call kim_model_destroy(model_handle)
616
program ex_test_ar_fcc_cluster_fortran
recursive subroutine my_error(message)
recursive subroutine my_warning(message)
The only standard extension defined by the KIM API.
character(len= *, kind=c_char), parameter, public kim_supported_extensions_id
recursive subroutine, public get_neigh(data_object, number_of_neighbor_lists, cutoffs, neighbor_list_index, request, numnei, pnei1part, ierr)
recursive subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, cutoff, neighObject)
recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, coords, MiddlePartId)