kim-api 2.3.0+AppleClang.AppleClang.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
Loading...
Searching...
No Matches
utility_forces_numer_deriv.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 public
32
33contains
34 recursive subroutine my_error(message)
35 implicit none
36 character(len=*, kind=c_char), intent(in) :: message
37
38 print *, "* Error : ", trim(message)
39 stop
40 end subroutine my_error
41
42 recursive subroutine my_warning(message)
43 implicit none
44 character(len=*, kind=c_char), intent(in) :: message
45
46 print *, "* Warning : ", trim(message)
47 end subroutine my_warning
48end module error
49
50!-------------------------------------------------------------------------------
51!
52! module mod_neighborlist :
53!
54! Module contains type and routines related to neighbor list calculation
55!
56!-------------------------------------------------------------------------------
57
59
60 use, intrinsic :: iso_c_binding
61
62 public get_neigh
63
64 type, bind(c) :: neighObject_type
65 real(c_double) :: cutoff
66 integer(c_int) :: number_of_particles
67 type(c_ptr) :: neighbor_list_pointer
68 end type neighobject_type
69contains
70
71!-------------------------------------------------------------------------------
72!
73! get_neigh neighbor list access function
74!
75! This function implements Locator and Iterator mode
76!
77!-------------------------------------------------------------------------------
78!-------------------------------------------------------------------------------
79!
80! get_neigh neighbor list access function
81!
82! This function implements Locator and Iterator mode
83!
84!-------------------------------------------------------------------------------
85 recursive subroutine get_neigh(data_object, number_of_neighbor_lists, &
86 cutoffs, neighbor_list_index, request, &
87 numnei, pnei1part, ierr) bind(c)
88 use error
89 implicit none
90
91 !-- Transferred variables
92 type(c_ptr), value, intent(in) :: data_object
93 integer(c_int), value, intent(in) :: number_of_neighbor_lists
94 real(c_double), intent(in) :: cutoffs(*)
95 integer(c_int), value, intent(in) :: neighbor_list_index
96 integer(c_int), value, intent(in) :: request
97 integer(c_int), intent(out) :: numnei
98 type(c_ptr), intent(out) :: pnei1part
99 integer(c_int), intent(out) :: ierr
100
101 !-- Local variables
102 integer(c_int) numberOfParticles
103 type(neighObject_type), pointer :: neighObject
104 integer(c_int), pointer :: neighborList(:, :)
105
106 if (number_of_neighbor_lists > 1) then
107 call my_warning("Model requires too many neighbor lists")
108 ierr = 1
109 return
110 end if
111
112 call c_f_pointer(data_object, neighobject)
113 numberofparticles = neighobject%number_of_particles
114 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
115 [numberofparticles + 1, numberofparticles])
116
117 if (cutoffs(neighbor_list_index) > neighobject%cutoff) then
118 call my_warning("neighbor list cutoff too small for model cutoff")
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
144 implicit none
145 public
146
147contains
148
149!-------------------------------------------------------------------------------
150!
151! Check if we are compatible with the model
152!
153!-------------------------------------------------------------------------------
154 recursive subroutine check_model_compatibility( &
155 compute_arguments_handle, forces_optional, particle_energy_supported, &
156 particle_energy_optional, model_is_compatible, ierr)
157 use, intrinsic :: iso_c_binding
158 use error
159 implicit none
160
161 !-- Transferred variables
162 type(kim_compute_arguments_handle_type), intent(in) :: &
163 compute_arguments_handle
164 logical, intent(out) :: forces_optional
165 logical, intent(out) :: particle_energy_supported
166 logical, intent(out) :: particle_energy_optional
167 logical, intent(out) :: model_is_compatible
168 integer(c_int), intent(out) :: ierr
169
170 !-- Local variables
171 integer(c_int) i
172 integer(c_int) number_of_argument_names
173 integer(c_int) number_of_callback_names
174 type(kim_compute_argument_name_type) argument_name
175 type(kim_support_status_type) support_status
176 type(kim_compute_callback_name_type) callback_name
177
178 ! assume fail
179 model_is_compatible = .false.
180 particle_energy_supported = .false.
181 particle_energy_optional = .false.
182 forces_optional = .false.
183 ierr = 0
184
185 ! check arguments
186 call kim_get_number_of_compute_argument_names( &
187 number_of_argument_names)
188 do i = 1, number_of_argument_names
189 call kim_get_compute_argument_name(i, argument_name, &
190 ierr)
191 if (ierr /= 0) then
192 call my_warning("can't get argument name")
193 return
194 end if
195 call kim_get_argument_support_status( &
196 compute_arguments_handle, argument_name, support_status, ierr)
197 if (ierr /= 0) then
198 call my_warning("can't get argument support_status")
199 return
200 end if
201
202 ! can only handle energy, particle_energy and forces as required args
203 if (support_status == kim_support_status_required) then
204 if (.not. ( &
205 (argument_name == kim_compute_argument_name_partial_energy) .or. &
206 (argument_name == &
207 kim_compute_argument_name_partial_particle_energy) .or. &
208 (argument_name == kim_compute_argument_name_partial_forces))) then
209 call my_warning("unsupported required argument")
210 ierr = 0
211 return
212 end if
213 end if
214
215 ! need both energy and forces not "notSupported"
216 if ((argument_name == kim_compute_argument_name_partial_energy) .and. &
217 (support_status == kim_support_status_not_supported)) then
218 call my_warning("model does not support energy")
219 ierr = 0
220 return
221 end if
222 if (argument_name == kim_compute_argument_name_partial_forces) then
223 if (support_status == kim_support_status_not_supported) then
224 call my_warning("model does not support forces")
225 ierr = 0
226 return
227 else if (support_status == kim_support_status_required) then
228 forces_optional = .false.
229 else if (support_status == kim_support_status_optional) then
230 forces_optional = .true.
231 else
232 call my_warning("unknown support_status for forces")
233 ierr = 0
234 return
235 end if
236 end if
237
238 ! check support for particle_energy
239 if (argument_name == &
240 kim_compute_argument_name_partial_particle_energy) then
241 if (support_status == kim_support_status_not_supported) then
242 call my_warning("model does not support partial_particle_energy. &
243 &The associated checks will be disabled.")
244 particle_energy_supported = .false.
245 particle_energy_optional = .false.
246 else if (support_status == kim_support_status_required) then
247 particle_energy_supported = .true.
248 particle_energy_optional = .false.
249 else if (support_status == kim_support_status_optional) then
250 particle_energy_supported = .true.
251 particle_energy_optional = .true.
252 else
253 call my_warning("unknown support_status for particle energy")
254 ierr = 0
255 return
256 end if
257 end if
258 end do
259
260 ! check call backs
261 call kim_get_number_of_compute_callback_names( &
262 number_of_callback_names)
263 do i = 1, number_of_callback_names
264 call kim_get_compute_callback_name(i, callback_name, &
265 ierr)
266 if (ierr /= 0) then
267 call my_warning("can't get call back name")
268 return
269 end if
270 call kim_get_callback_support_status( &
271 compute_arguments_handle, callback_name, support_status, ierr)
272 if (ierr /= 0) then
273 call my_warning("can't get call back support_status")
274 return
275 end if
276
277 ! cannot handle any "required" call backs
278 if (support_status == kim_support_status_required) then
279 call my_warning("unsupported required call back")
280 ierr = 0
281 return
282 end if
283 end do
284
285 ! got to here, then everything must be OK
286 model_is_compatible = .true.
287 ierr = 0
288 return
289 end subroutine check_model_compatibility
290
291 !-----------------------------------------------------------------------------
292 !
293 ! Get number and identities of particle species supported by
294 ! KIM Model `modelname'
295 !
296 !-----------------------------------------------------------------------------
297 recursive subroutine get_model_supported_species( &
298 model_handle, max_species, model_species, num_species, ier)
299 use, intrinsic :: iso_c_binding
300 implicit none
301
302 !-- Transferred variables
303 type(kim_model_handle_type), intent(in) :: model_handle
304 integer(c_int), intent(in) :: max_species
305 type(kim_species_name_type), intent(out) :: model_species(max_species)
306 integer(c_int), intent(out) :: num_species
307 integer(c_int), intent(out) :: ier
308
309 !-- Local variables
310 integer(c_int) i
311 integer(c_int) total_num_species
312 type(kim_species_name_type) :: species_name
313 integer(c_int) species_is_supported
314 integer(c_int) code
315
316 ! Initialize error flag
317 ier = 1
318
319 num_species = 0 ! initialize
320
321 call kim_get_number_of_species_names(total_num_species)
322
323 if (total_num_species > max_species) return
324
325 num_species = 0
326 do i = 1, total_num_species
327 call kim_get_species_name(i, species_name, ier)
328 call kim_get_species_support_and_code(model_handle, species_name, &
329 species_is_supported, code, ier)
330 if ((ier == 0) .and. (species_is_supported /= 0)) then
331 num_species = num_species + 1
332 model_species(num_species) = species_name
333 end if
334 end do
335
336 ier = 0
337 return
338
339 end subroutine get_model_supported_species
340
341 recursive subroutine update_neighborlist(DIM, N, coords, cutoff, cutpad, &
342 do_update_list, coordsave, &
343 neighObject, ierr)
344 use, intrinsic :: iso_c_binding
346 implicit none
347 integer(c_int), parameter :: cd = c_double ! used for literal constants
348
349 !-- Transferred variables
350 integer(c_int), intent(in) :: dim
351 integer(c_int), intent(in) :: n
352 real(c_double), intent(in) :: coords(dim, n)
353 real(c_double), intent(in) :: cutoff
354 real(c_double), intent(in) :: cutpad
355 logical, intent(inout) :: do_update_list
356 real(c_double), intent(inout) :: coordsave(dim, n)
357 type(neighobject_type), intent(inout) :: neighobject
358 integer(c_int), intent(out) :: ierr
359
360 !-- Local variables
361 real(c_double) disp, disp1, disp2, cutrange, dispvec(dim)
362 integer(c_int) i
363
364 ! Initialize error code
365 !
366 ierr = 0
367
368 ! Update neighbor lists if necessary
369 !
370 if (.not. do_update_list) then ! if update not requested
371
372 ! check whether a neighbor list update is necessary even if it hasn't been
373 ! requested using the "two max sum" criterion
374 disp1 = 0.0_cd
375 disp2 = 0.0_cd
376 do i = 1, n
377 dispvec(1:dim) = coords(1:dim, i) - coordsave(1:dim, i)
378 disp = sqrt(dot_product(dispvec, dispvec))
379 if (disp >= disp1) then ! 1st position taken
380 disp2 = disp1 ! push current 1st into 2nd place
381 disp1 = disp ! and put this one into current 1st
382 else if (disp >= disp2) then ! 2nd position taken
383 disp2 = disp
384 end if
385 end do
386 do_update_list = (disp1 + disp2 > cutpad)
387
388 end if
389
390 if (do_update_list) then
391
392 ! save current coordinates
393 coordsave(1:dim, 1:n) = coords(1:dim, 1:n)
394
395 ! compute neighbor lists
396 cutrange = cutoff + cutpad
397 call neigh_pure_cluster_neighborlist(.false., n, coords, cutrange, &
398 neighobject)
399
400 ! neighbor list uptodate, no need to compute again for now
401 do_update_list = .false.
402 end if
403
404 return
405
406 end subroutine update_neighborlist
407
408 !-----------------------------------------------------------------------------
409 !
410 ! NEIGH_PURE_cluster_neighborlist
411 !
412 !-----------------------------------------------------------------------------
413 recursive subroutine neigh_pure_cluster_neighborlist( &
414 half, numberOfParticles, coords, cutoff, neighObject)
415 use, intrinsic :: iso_c_binding
417 implicit none
418
419 !-- Transferred variables
420 logical, intent(in) :: half
421 integer(c_int), intent(in) :: numberofparticles
422 real(c_double), dimension(3, numberOfParticles), &
423 intent(in) :: coords
424 real(c_double), intent(in) :: cutoff
425 type(neighobject_type), intent(inout) :: neighobject
426
427 !-- Local variables
428 integer(c_int), pointer :: neighborlist(:, :)
429 integer(c_int) i, j, a
430 real(c_double) dx(3)
431 real(c_double) r2
432 real(c_double) cutoff2
433
434 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
435 [numberofparticles + 1, numberofparticles])
436
437 neighobject%cutoff = cutoff
438
439 cutoff2 = cutoff**2
440
441 do i = 1, numberofparticles
442 a = 1
443 do j = 1, numberofparticles
444 dx(:) = coords(:, j) - coords(:, i)
445 r2 = dot_product(dx, dx)
446 if (r2 <= cutoff2) then
447 ! part j is a neighbor of part i
448 if ((j > i) .OR. ((.not. half) .AND. (i /= j))) then
449 a = a + 1
450 neighborlist(a, i) = j
451 end if
452 end if
453 end do
454 ! part i has a-1 neighbors
455 neighborlist(1, i) = a - 1
456 end do
457
458 return
459
461
462 !-----------------------------------------------------------------------------
463 !
464 ! create_FCC_configuration subroutine
465 !
466 ! creates a cubic configuration of FCC particles with lattice spacing
467 ! `FCCspacing' and `nCellsPerSide' cells along each direction.
468 !
469 ! With periodic==.true. this will result in a total number of particles equal
470 ! to 4*(nCellsPerSide)**3 + 6*(nCellsPerSide)**2 + 3*(nCellsPerSide) + 1
471 !
472 ! With periodic==.false. this will result in a total number of particles equal
473 ! to 4*(nCellsPerSide)**3
474 !
475 ! Returns the Id of the particle situated in the middle of the configuration
476 ! (this particle will have the most neighbors.)
477 !
478 !-----------------------------------------------------------------------------
479 recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, &
480 periodic, coords, MiddlePartId)
481 use, intrinsic :: iso_c_binding
482 implicit none
483 integer(c_int), parameter :: cd = c_double ! used for literal constants
484
485 !-- Transferred variables
486 real(c_double), intent(in) :: fccspacing
487 integer(c_int), intent(in) :: ncellsperside
488 logical, intent(in) :: periodic
489 real(c_double), intent(out) :: coords(3, *)
490 integer(c_int), intent(out) :: middlepartid
491 !
492 ! cluster setup variables
493 !
494 real(c_double) fccshifts(3, 4)
495 real(c_double) latvec(3)
496 integer(c_int) a, i, j, k, m
497
498 ! Create a cubic FCC cluster
499 !
500 fccshifts(1, 1) = 0.0_cd
501 fccshifts(2, 1) = 0.0_cd
502 fccshifts(3, 1) = 0.0_cd
503 fccshifts(1, 2) = 0.5_cd * fccspacing
504 fccshifts(2, 2) = 0.5_cd * fccspacing
505 fccshifts(3, 2) = 0.0_cd
506 fccshifts(1, 3) = 0.5_cd * fccspacing
507 fccshifts(2, 3) = 0.0_cd
508 fccshifts(3, 3) = 0.5_cd * fccspacing
509 fccshifts(1, 4) = 0.0_cd
510 fccshifts(2, 4) = 0.5_cd * fccspacing
511 fccshifts(3, 4) = 0.5_cd * fccspacing
512
513 middlepartid = 1 ! Always put middle particle as #1
514 a = 1 ! leave space for middle particle as particle #1
515 do i = 1, ncellsperside
516 latvec(1) = (i - 1) * fccspacing
517 do j = 1, ncellsperside
518 latvec(2) = (j - 1) * fccspacing
519 do k = 1, ncellsperside
520 latvec(3) = (k - 1) * fccspacing
521 do m = 1, 4
522 a = a + 1
523 coords(:, a) = latvec + fccshifts(:, m)
524 if ((i == ncellsperside / 2 + 1) &
525 .and. (j == ncellsperside / 2 + 1) &
526 .and. (k == ncellsperside / 2 + 1) &
527 .and. (m == 1)) &
528 then
529 ! put middle particle as #1
530 coords(:, 1) = latvec + fccshifts(:, m)
531 a = a - 1
532 end if
533 end do
534 end do
535 if (.not. periodic) then
536 ! Add in the remaining three faces
537 ! pos-x face
538 latvec(1) = ncellsperside * fccspacing
539 latvec(2) = (i - 1) * fccspacing
540 latvec(3) = (j - 1) * fccspacing
541 a = a + 1; coords(:, a) = latvec
542 a = a + 1; coords(:, a) = latvec + fccshifts(:, 4)
543 ! pos-y face
544 latvec(1) = (i - 1) * fccspacing
545 latvec(2) = ncellsperside * fccspacing
546 latvec(3) = (j - 1) * fccspacing
547 a = a + 1; coords(:, a) = latvec
548 a = a + 1; coords(:, a) = latvec + fccshifts(:, 3)
549 ! pos-z face
550 latvec(1) = (i - 1) * fccspacing
551 latvec(2) = (j - 1) * fccspacing
552 latvec(3) = ncellsperside * fccspacing
553 a = a + 1; coords(:, a) = latvec
554 a = a + 1; coords(:, a) = latvec + fccshifts(:, 2)
555 end if
556 end do
557 if (.not. periodic) then
558 ! Add in the remaining three edges
559 latvec(1) = (i - 1) * fccspacing
560 latvec(2) = ncellsperside * fccspacing
561 latvec(3) = ncellsperside * fccspacing
562 a = a + 1; coords(:, a) = latvec
563 latvec(1) = ncellsperside * fccspacing
564 latvec(2) = (i - 1) * fccspacing
565 latvec(3) = ncellsperside * fccspacing
566 a = a + 1; coords(:, a) = latvec
567 latvec(1) = ncellsperside * fccspacing
568 latvec(2) = ncellsperside * fccspacing
569 latvec(3) = (i - 1) * fccspacing
570 a = a + 1; coords(:, a) = latvec
571 end if
572 end do
573 if (.not. periodic) then
574 ! Add in the remaining corner
575 a = a + 1; coords(:, a) = ncellsperside * fccspacing
576 end if
577
578 return
579
580 end subroutine create_fcc_configuration
581
582 recursive subroutine compute_numer_deriv( &
583 partnum, dir, model_handle, compute_arguments_handle, DIM, N, coords, &
584 cutoff, cutpad, energy, do_update_list, coordsave, neighObject, deriv, &
585 deriv_err, ierr)
586 use, intrinsic :: iso_c_binding
587 use error
589 implicit none
590 integer(c_int), parameter :: cd = c_double ! used for literal constants
591
592 !--Transferred variables
593 integer(c_int), intent(in) :: partnum
594 integer(c_int), intent(in) :: dir
595 type(kim_model_handle_type), intent(in) :: model_handle
596 type(kim_compute_arguments_handle_type), intent(in) :: &
597 compute_arguments_handle
598 integer(c_int), intent(in) :: dim
599 integer(c_int), intent(in) :: n
600 real(c_double), intent(inout) :: coords(dim, n)
601 real(c_double), intent(in) :: cutoff
602 real(c_double), intent(in) :: cutpad
603 real(c_double), intent(inout) :: energy
604 logical, intent(inout) :: do_update_list
605 real(c_double), intent(inout) :: coordsave(dim, n)
606 type(neighobject_type), intent(inout) :: neighobject
607 real(c_double), intent(out) :: deriv
608 real(c_double), intent(out) :: deriv_err
609 integer(c_int), intent(out) :: ierr
610
611 !-- Local variables
612 real(c_double), parameter :: eps_init = 1.e-6_cd
613 integer(c_int), parameter :: number_eps_levels = 15
614 real(c_double) eps, deriv_last, deriv_err_last
615 integer(c_int) i
616
617 ! Initialize error flag
618 ierr = 0
619
620 deriv_last = 0.0_cd ! initialize
621
622 ! Outer loop of Ridders' method for computing numerical derivative
623 !
624 eps = eps_init
625 deriv_err_last = huge(1.0_cd)
626 do i = 1, number_eps_levels
627 deriv = dfridr(eps, deriv_err)
628 if (ierr /= 0) then
629 call my_error("compute_numer_deriv")
630 end if
631 if (deriv_err > deriv_err_last) then
632 deriv = deriv_last
633 deriv_err = deriv_err_last
634 exit
635 end if
636 eps = eps * 10.0_cd
637 deriv_last = deriv
638 deriv_err_last = deriv_err
639 end do
640
641 return
642
643 contains
644
645 !----------------------------------------------------------------------------
646 !
647 ! Compute numerical derivative using Ridders' method
648 !
649 ! Based on code from Numerical Recipes, Press et al., Second Ed., Cambridge,
650 ! 1992
651 !
652 ! Ref: Ridders, C. J. F., "Two algorithms for the calculation of F'(x)=D",
653 ! Advances in Engineering Software, Vol. 4, no. 2, pp. 75-76, 1982.
654 !
655 !
656 ! Returns the gradient grad() of a KIM-compliant interatomic model at the
657 ! current configuration by Ridders' method of polynomial extrapolation.
658 ! An estimate for the error in each component of the gradient is returned in
659 ! grad_err().
660 !
661 !----------------------------------------------------------------------------
662 real(c_double) recursive function dfridr(h, err)
663 implicit none
664
665 !-- Transferred variables
666 real(c_double), intent(inout) :: h
667 real(c_double), intent(out) :: err
668
669 !-- Local variables
670 integer(c_int), parameter :: ntab = 10 ! Maximum size of tableau
671 real(c_double), parameter :: con = 1.4_cd ! Stepsize incr. by CON at each iter
672 real(c_double), parameter :: con2 = con * con
673 real(c_double), parameter :: big = huge(1.0_cd)
674 real(c_double), parameter :: safe = 2.0_cd ! Returns when error is SAFE worse
675 ! than the best so far
676 integer(c_int) i, j
677 real(c_double) errt, fac, hh, a(ntab, ntab), fp, fm, coordorig
678
679 dfridr = 0.0_cd ! initialize
680 err = big ! initialize
681
682 if (abs(h) <= tiny(0.0_cd)) then ! avoid division by zero
683 ierr = 1
684 return
685 end if
686
687 hh = h
688 coordorig = coords(dir, partnum)
689 coords(dir, partnum) = coordorig + hh
690 call update_neighborlist(dim, n, coords, cutoff, cutpad, &
691 do_update_list, coordsave, &
692 neighobject, ierr)
693 call kim_compute(model_handle, compute_arguments_handle, ierr)
694 if (ierr /= 0) then
695 call my_error("kim_api_model_compute")
696 end if
697 fp = energy
698 coords(dir, partnum) = coordorig - hh
699 call update_neighborlist(dim, n, coords, cutoff, cutpad, &
700 do_update_list, coordsave, &
701 neighobject, ierr)
702 call kim_compute(model_handle, compute_arguments_handle, ierr)
703 if (ierr /= 0) then
704 call my_error("kim_api_model_compute")
705 end if
706 fm = energy
707 coords(dir, partnum) = coordorig
708 call update_neighborlist(dim, n, coords, cutoff, cutpad, &
709 do_update_list, coordsave, &
710 neighobject, ierr)
711 a(1, 1) = (fp - fm) / (2.0_cd * hh)
712 ! successive columns in the Neville tableau will go to smaller step sizes
713 ! and higher orders of extrapolation
714 do i = 2, ntab
715 ! try new, smaller step size
716 hh = hh / con
717 coords(dir, partnum) = coordorig + hh
718 call update_neighborlist(dim, n, coords, cutoff, cutpad, &
719 do_update_list, coordsave, &
720 neighobject, ierr)
721 call kim_compute(model_handle, compute_arguments_handle, ierr)
722 if (ierr /= 0) then
723 call my_error("kim_api_model_compute")
724 end if
725 fp = energy
726 coords(dir, partnum) = coordorig - hh
727 call update_neighborlist(dim, n, coords, cutoff, cutpad, &
728 do_update_list, coordsave, &
729 neighobject, ierr)
730 call kim_compute(model_handle, compute_arguments_handle, ierr)
731 if (ierr /= 0) then
732 call my_error("kim_api_model_compute")
733 end if
734 fm = energy
735 coords(dir, partnum) = coordorig
736 call update_neighborlist(dim, n, coords, cutoff, cutpad, &
737 do_update_list, coordsave, &
738 neighobject, ierr)
739 a(1, i) = (fp - fm) / (2.0_cd * hh)
740 fac = con2
741 ! compute extrapolations of various orders, requiring no new function
742 ! evaluations
743 do j = 2, i
744 a(j, i) = (a(j - 1, i) * fac - a(j - 1, i - 1)) / (fac - 1.0_cd)
745 fac = con2 * fac
746 ! The error strategy is to compute each new extrapolation to one order
747 ! lower, both at the present step size and the previous one.
748 errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
749 if (errt <= err) then ! if error is decreased, save the improved answer
750 err = errt
751 dfridr = a(j, i)
752 end if
753 end do
754 if (abs(a(i, i) - a(i - 1, i - 1)) >= safe * err) return ! if higher order is worse
755 ! by significant factor
756 ! `SAFE', then quit early.
757 end do
758 return
759 end function dfridr
760
761 end subroutine compute_numer_deriv
762
763end module mod_utilities
764
765!*******************************************************************************
766!**
767!** PROGRAM vc_forces_numer_deriv
768!**
769!** KIM compliant program to perform numerical derivative check on a model
770!**
771!*******************************************************************************
772
773!-------------------------------------------------------------------------------
774!
775! Main program
776!
777!-------------------------------------------------------------------------------
779 use, intrinsic :: iso_c_binding
780 use error
782 use mod_utilities
783 implicit none
784 integer(c_int), parameter :: cd = c_double ! used for literal constants
785
786 integer(c_int), parameter :: ncellsperside = 2
787 integer(c_int), parameter :: dim = 3
788 real(c_double), parameter :: cutpad = 0.75_cd
789 integer(c_int), parameter :: max_species = 200 ! most species supported
790 real(c_double), parameter :: eps_prec = epsilon(1.0_cd)
791 real(c_double) fccspacing
792
793 integer(c_int), parameter :: n = 4 * (ncellsperside)**3 &
794 + 6 * (ncellsperside)**2 &
795 + 3 * (ncellsperside) + 1
796 real(c_double), allocatable :: forces_num(:, :)
797 real(c_double), allocatable :: forces_num_err(:, :)
798 type(kim_species_name_type) :: model_species(max_species)
799 integer(c_int), target :: num_species
800 character(len=5, kind=c_char) :: passfail
801 real(c_double) :: forcediff
802 real(c_double) :: forcediff_sumsq
803 real(c_double) :: weight
804 real(c_double) :: weight_sum
805 real(c_double) :: alpha
806 real(c_double) :: term
807 real(c_double) :: term_max
808 real(c_double), allocatable :: cluster_coords(:, :)
809 real(c_double), allocatable :: cluster_disps(:, :)
810 type(kim_species_name_type), allocatable :: cluster_species(:)
811 integer(c_int) i, j, imax, jmax, species_code
812 integer(c_int) seed_size
813 integer(c_int), allocatable :: seed(:)
814
815 !
816 ! neighbor list
817 !
818 type(neighobject_type), target :: neighobject
819 integer(c_int), allocatable, target :: neighborlist(:, :)
820 real(c_double), allocatable :: coordsave(:, :)
821 logical do_update_list
822
823 !
824 ! KIM variables
825 !
826 character(len=256, kind=c_char) :: testname = "vc_forces_numer_deriv"
827 character(len=256, kind=c_char) :: modelname
828
829 type(kim_model_handle_type) :: model_handle
830 type(kim_compute_arguments_handle_type) :: compute_arguments_handle
831 integer(c_int) ierr, ierr2
832 integer(c_int) species_is_supported
833 integer(c_int), target :: numberofparticles
834 integer(c_int), target :: particlespeciescodes(n)
835 integer(c_int), target :: particlecontributing(n)
836 real(c_double) :: influence_distance
837 integer(c_int) :: number_of_neighbor_lists
838 real(c_double), allocatable :: cutoffs(:)
839 integer(c_int), allocatable :: &
840 model_will_not_request_neighbors_of_noncontributing_particles(:)
841 real(c_double) :: cutoff
842 real(c_double), target :: energy
843 real(c_double), target :: particle_energy(n)
844 real(c_double), target :: particle_energy_sum
845 real(c_double), target :: coords(3, n)
846 real(c_double), target :: forces_kim(3, n)
847 real(c_double) :: forces(3, n)
848 integer(c_int) middledum
849 logical doing_non_contributing
850 logical particle_energy_supported
851 logical particle_energy_optional
852 logical noncontrib_particle_energy_nonzero
853 logical forces_optional
854 logical model_is_compatible
855 integer(c_int) number_of_model_routine_names
856 type(kim_model_routine_name_type) model_routine_name
857 integer(c_int) present
858 integer(c_int) required
859 integer(c_int) requested_units_accepted
860 real(c_double) rnd, deriv, deriv_err
861 real(c_double), pointer :: null_pointer
862
863 nullify (null_pointer)
864
865 doing_non_contributing = .false. ! .true. on 2nd pass through main program
866
867 numberofparticles = n
868
869 term_max = 0.0_cd ! initialize
870
871 ! Initialize error flag
872 ierr = 0
873
874 ! Initialize seed for random number generator
875 !
876 ! NOTE: Here, we set the seed to a fixed value for reproducibility
877 call random_seed(size=seed_size) ! Get seed size for this CPU
878 allocate (seed(seed_size))
879 seed(:) = 13
880 call random_seed(put=seed)
881 deallocate (seed)
882
883 ! Get KIM Model name to use
884 print '("Please enter a valid KIM model name: ")'
885 read (*, *) modelname
886
887 ! Print output header
888 !
889 print *
890 print *, 'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES'
891 print *
892 print '(120(''-''))'
893 print '("This is Test : ",A)', trim(testname)
894 print '("Results for KIM Model : ",A)', trim(modelname)
895
896 ! Create empty KIM object
897 !
898 call kim_model_create(kim_numbering_one_based, &
899 kim_length_unit_a, &
900 kim_energy_unit_ev, &
901 kim_charge_unit_e, &
902 kim_temperature_unit_k, &
903 kim_time_unit_ps, &
904 trim(modelname), &
905 requested_units_accepted, &
906 model_handle, ierr)
907 if (ierr /= 0) then
908 call my_error("kim_api_create")
909 end if
910
911 ! check that we are compatible
912 if (requested_units_accepted == 0) then
913 call my_error("Must adapt to model units")
914 end if
915
916 ! check that we know about all required routines
917 call kim_get_number_of_model_routine_names(number_of_model_routine_names)
918 do i = 1, number_of_model_routine_names
919 call kim_get_model_routine_name(i, model_routine_name, ierr)
920 if (ierr /= 0) then
921 call my_error("kim_get_model_routine_name")
922 end if
923 call kim_is_routine_present(model_handle, model_routine_name, present, &
924 required, ierr)
925 if (ierr /= 0) then
926 call my_error("kim_is_routine_present")
927 end if
928
929 if ((present == 1) .and. (required == 1)) then
930 if (.not. ((model_routine_name == kim_model_routine_name_create) &
931 .or. (model_routine_name == &
932 kim_model_routine_name_compute_arguments_create) &
933 .or. (model_routine_name == kim_model_routine_name_compute) &
934 .or. (model_routine_name == kim_model_routine_name_refresh) &
935 .or. (model_routine_name == &
936 kim_model_routine_name_compute_arguments_destroy) &
937 .or. (model_routine_name == kim_model_routine_name_destroy))) &
938 then
939 call my_error("Unknown required ModelRoutineName found.")
940 end if
941 end if
942 end do
943
944 ! create compute_arguments object
945 call kim_compute_arguments_create(model_handle, &
946 compute_arguments_handle, ierr)
947 if (ierr /= 0) then
948 call my_error("kim_model_compute_arguments_create")
949 end if
950
951 call check_model_compatibility(compute_arguments_handle, forces_optional, &
952 particle_energy_supported, &
953 particle_energy_optional, &
954 model_is_compatible, ierr)
955 if (ierr /= 0) then
956 call my_error("error checking compatibility")
957 end if
958 if (.not. model_is_compatible) then
959 call my_error("incompatibility reported by check_model_compatibility")
960 end if
961
962 ! Get list of particle species supported by the model
963 !
964 call get_model_supported_species(model_handle, max_species, model_species, &
965 num_species, ierr)
966 if (ierr /= 0) then
967 call my_error("Get_Model_Supported_Species")
968 end if
969 ! Setup random cluster
970 !
971 allocate (cluster_coords(3, n), cluster_disps(3, n), cluster_species(n))
972 do i = 1, n
973 call random_number(rnd) ! return random number between 0 and 1
974 species_code = 1 + int(rnd * num_species)
975 cluster_species(i) = model_species(species_code)
976 end do
977 fccspacing = 1.0_cd ! initially generate an FCC cluster with lattice
978 ! spacing equal to one. This is scaled below based
979 ! on the cutoff radius of the model.
980 call create_fcc_configuration(fccspacing, ncellsperside, .false., &
981 cluster_coords, middledum)
982 ! Generate random displacements for all particles
983 !
984 do i = 1, n
985 do j = 1, dim
986 call random_number(rnd) ! return random number between 0 and 1
987 cluster_disps(j, i) = 0.1_cd * (rnd - 0.5_cd)
988 end do
989 end do
990
991 ! register memory with the KIM system
992 ierr = 0
993 call kim_set_argument_pointer( &
994 compute_arguments_handle, &
995 kim_compute_argument_name_number_of_particles, &
996 numberofparticles, ierr2)
997 ierr = ierr + ierr2
998 call kim_set_argument_pointer( &
999 compute_arguments_handle, &
1000 kim_compute_argument_name_particle_species_codes, &
1001 particlespeciescodes, ierr2)
1002 ierr = ierr + ierr2
1003 call kim_set_argument_pointer( &
1004 compute_arguments_handle, &
1005 kim_compute_argument_name_particle_contributing, &
1006 particlecontributing, ierr2)
1007 ierr = ierr + ierr2
1008 call kim_set_argument_pointer( &
1009 compute_arguments_handle, &
1010 kim_compute_argument_name_coordinates, coords, ierr2)
1011 ierr = ierr + ierr2
1012 call kim_set_argument_pointer( &
1013 compute_arguments_handle, &
1014 kim_compute_argument_name_partial_energy, energy, ierr2)
1015 ierr = ierr + ierr2
1016 if (particle_energy_supported) then
1017 call kim_set_argument_pointer( &
1018 compute_arguments_handle, &
1019 kim_compute_argument_name_partial_particle_energy, particle_energy, ierr2)
1020 ierr = ierr + ierr2
1021 end if
1022 call kim_set_argument_pointer( &
1023 compute_arguments_handle, &
1024 kim_compute_argument_name_partial_forces, forces_kim, ierr2)
1025 ierr = ierr + ierr2
1026 if (ierr /= 0) then
1027 call my_error("set_argument_pointer")
1028 end if
1029
1030 ! Allocate storage for neighbor lists and
1031 ! store pointers to neighbor list object and access function
1032 !
1033 allocate (coordsave(dim, n))
1034 allocate (neighborlist(n + 1, n))
1035 neighobject%neighbor_list_pointer = c_loc(neighborlist)
1036 neighobject%number_of_particles = n
1037 allocate (forces_num(dim, n), forces_num_err(dim, n))
1038
1039 ! Set pointer in KIM object to neighbor list routine and object
1040 !
1041 call kim_set_callback_pointer( &
1042 compute_arguments_handle, &
1043 kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
1044 c_funloc(get_neigh), c_loc(neighobject), ierr)
1045 if (ierr /= 0) then
1046 call my_error("set_callback_pointer")
1047 end if
1048
1049 call kim_get_influence_distance(model_handle, influence_distance)
1050 call kim_get_number_of_neighbor_lists(model_handle, &
1051 number_of_neighbor_lists)
1052 allocate (cutoffs(number_of_neighbor_lists), &
1053 model_will_not_request_neighbors_of_noncontributing_particles( &
1054 number_of_neighbor_lists))
1055 call kim_get_neighbor_list_values( &
1056 model_handle, cutoffs, &
1057 model_will_not_request_neighbors_of_noncontributing_particles, ierr)
1058 if (ierr /= 0) then
1059 call my_error("get_neighbor_list_values")
1060 end if
1061 cutoff = maxval(cutoffs)
1062
1063 ! Scale reference FCC configuration based on cutoff radius.
1064 fccspacing = 0.75_cd * cutoff ! set the FCC spacing to a fraction
1065 ! of the cutoff radius
1066 do i = 1, n
1067 cluster_coords(:, i) = fccspacing * cluster_coords(:, i)
1068 end do
1069 print '("Using FCC lattice parameter: ",f12.5)', fccspacing
1070 print *
1071
10721000 continue ! Start of configuration setup
1073
1074 if (doing_non_contributing) then
1075 ! Turn particle energy computation back on, if possible
1076 if (particle_energy_optional) then
1077 call kim_set_argument_pointer( &
1078 compute_arguments_handle, &
1079 kim_compute_argument_name_partial_particle_energy, &
1080 particle_energy, ierr)
1081 if (ierr /= 0) then
1082 call my_error("set_argument_pointer")
1083 end if
1084 end if
1085
1086 ! Turn force computation back on, if possible
1087 !
1088 if (forces_optional) then
1089 call kim_set_argument_pointer( &
1090 compute_arguments_handle, kim_compute_argument_name_partial_forces, &
1091 forces_kim, ierr)
1092 if (ierr /= 0) then
1093 call my_error("set_argument_pointer")
1094 end if
1095 end if
1096 end if
1097
1098 do i = 1, n
1099 call kim_get_species_support_and_code( &
1100 model_handle, cluster_species(i), species_is_supported, &
1101 particlespeciescodes(i), ierr)
1102 end do
1103 if (ierr /= 0) then
1104 call my_error("kim_api_get_species_code")
1105 end if
1106 if (.not. doing_non_contributing) then
1107 do i = 1, n
1108 particlecontributing(i) = 1 ! all particle contribute
1109 end do
1110 else
1111 do i = 1, n
1112 ! Random collection of contributing atoms
1113 call random_number(rnd) ! return random number between 0 and 1
1114 if (rnd > 0.5) then
1115 particlecontributing(i) = 1
1116 else
1117 particlecontributing(i) = 0
1118 end if
1119 end do
1120 end if
1121 do i = 1, n
1122 coords(:, i) = cluster_coords(:, i) + cluster_disps(:, i)
1123 end do
1124
1125 ! Compute neighbor lists
1126 !
1127 do_update_list = .true.
1128 call update_neighborlist(dim, n, coords, cutoff, cutpad, &
1129 do_update_list, coordsave, &
1130 neighobject, ierr)
1131 if (ierr /= 0) then
1132 call my_error("update_neighborlist")
1133 end if
1134
1135 ! Call model compute to get forces (gradient)
1136 !
1137 call kim_compute(model_handle, compute_arguments_handle, ierr)
1138 if (ierr /= 0) then
1139 call my_error("kim_api_model_compute")
1140 end if
1141
1142 ! Copy forces in case model will overwrite forces_kim below
1143 !
1144 forces = forces_kim
1145
1146 ! Add up particle_energy to compare with energy
1147 noncontrib_particle_energy_nonzero = .false.
1148 if (particle_energy_supported) then
1149 particle_energy_sum = 0.0_cd
1150 do i = 1, n
1151 if (particlecontributing(i) == 0) then
1152 if (abs(particle_energy(i)) > epsilon(0.0_cd)) then
1153 noncontrib_particle_energy_nonzero = .true.
1154 end if
1155 else
1156 particle_energy_sum = particle_energy_sum + particle_energy(i)
1157 end if
1158 end do
1159 end if
1160
1161 ! Print results to screen
1162 !
1163 print '(41(''=''))'
1164 if (.not. doing_non_contributing) then
1165 print '("Configuration with all contributing particles")'
1166 else
1167 print '("Configuration with some non-contributing particles")'
1168 end if
1169 if (particle_energy_supported) then
1170 print '(A25,2X,A25,2X,A15)', "Energy", "Sum Contrib. Energies", "Diff"
1171 print '(ES25.15,2X,ES25.15,2X,ES15.5)', energy, particle_energy_sum, &
1172 energy - particle_energy_sum
1173 if (noncontrib_particle_energy_nonzero) then
1174 call my_error( &
1175 "Some non-contributing particles have non-zero &
1176 &partial_particle_energy")
1177 end if
1178 else
1179 print '("Energy = ",ES25.15)', energy
1180 end if
1181 print '(41(''=''))'
1182 print *
1183
1184 ! Turn off particle energy computation, if possible
1185 if (particle_energy_optional) then
1186 call kim_set_argument_pointer( &
1187 compute_arguments_handle, &
1188 kim_compute_argument_name_partial_particle_energy, &
1189 null_pointer, ierr)
1190 if (ierr /= 0) then
1191 call my_error("set_argument_pointer")
1192 end if
1193 end if
1194
1195 ! Turn off force computation, if possible
1196 !
1197 if (forces_optional) then
1198 call kim_set_argument_pointer( &
1199 compute_arguments_handle, kim_compute_argument_name_partial_forces, &
1200 null_pointer, ierr)
1201 if (ierr /= 0) then
1202 call my_error("set_argument_pointer")
1203 end if
1204 end if
1205
1206 ! Compute gradient using numerical differentiation
1207 !
1208 do i = 1, n
1209 do j = 1, dim
1210 call compute_numer_deriv(i, j, model_handle, compute_arguments_handle, &
1211 dim, n, coords, cutoff, &
1212 cutpad, energy, do_update_list, &
1213 coordsave, neighobject, deriv, deriv_err, ierr)
1214 if (ierr /= 0) then
1215 call my_error("compute_numer_deriv")
1216 end if
1217 forces_num(j, i) = -deriv
1218 forces_num_err(j, i) = deriv_err
1219 end do
1220 end do
1221
1222 ! Continue printing results to screen
1223 !
1224 print '(A6,2X,A7,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)', "Part", "Contrib", &
1225 "Spec", "Dir", "Force_model", "Force_numer", "Force diff", "pred error", &
1226 "weight", "stat"
1227 forcediff_sumsq = 0.0_cd
1228 weight_sum = 0.0_cd
1229 do i = 1, n
1230 do j = 1, dim
1231 forcediff = abs(forces(j, i) - forces_num(j, i))
1232 if (forcediff < forces_num_err(j, i)) then
1233 passfail = "ideal"
1234 else
1235 passfail = " "
1236 end if
1237 weight = max(abs(forces_num(j, i)), eps_prec) / &
1238 max(abs(forces_num_err(j, i)), eps_prec)
1239 term = weight * forcediff**2
1240 if (term > term_max) then
1241 term_max = term
1242 imax = i
1243 jmax = j
1244 end if
1245 forcediff_sumsq = forcediff_sumsq + term
1246 weight_sum = weight_sum + weight
1247 if (j == 1) then
1248 print '(I6,2X,I7,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1249 i, particlecontributing(i), &
1250 particlespeciescodes(i), j, forces(j, i), &
1251 forces_num(j, i), forcediff, &
1252 forces_num_err(j, i), weight, passfail
1253 else
1254 print '(23X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1255 j, forces(j, i), forces_num(j, i), &
1256 forcediff, forces_num_err(j, i), weight, passfail
1257 end if
1258 end do
1259 print *
1260 end do
1261 alpha = sqrt(forcediff_sumsq / weight_sum) / dble(dim * n)
1262 print *
1263 print '("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," &
1264 &(units of force)")', alpha
1265 print *
1266 print '(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,'// &
1267 ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
1268 imax, jmax, abs(forces(jmax, imax) - forces_num(jmax, imax)), &
1269 abs(forces(jmax, imax) - forces_num(jmax, imax)) / abs(forces(jmax, imax))
1270
1271 if (.not. doing_non_contributing) then
1272 doing_non_contributing = .true.
1273 print *
1274 print *
1275 goto 1000
1276 end if
1277
1278 ! Free temporary storage
1279 !
1280 deallocate (forces_num)
1281 deallocate (forces_num_err)
1282 deallocate (neighborlist)
1283 deallocate (coordsave)
1284 deallocate (cutoffs)
1285 deallocate (model_will_not_request_neighbors_of_noncontributing_particles)
1286
1287 call kim_compute_arguments_destroy(model_handle, &
1288 compute_arguments_handle, ierr)
1289 if (ierr /= 0) then
1290 call my_error("kim_model_compute_arguments_destroy")
1291 end if
1292 call kim_model_destroy(model_handle)
1293
1294 ! Print output footer
1295 !
1296 print *
1297 print '(120(''-''))'
1298
1299 ! Free cluster storage
1300 !
1301 deallocate (cluster_coords, cluster_disps, cluster_species)
1302
1303 stop
1304
1305end program vc_forces_numer_deriv
recursive subroutine my_error(message)
recursive subroutine my_warning(message)
recursive subroutine, public get_neigh(data_object, number_of_neighbor_lists, cutoffs, neighbor_list_index, request, numnei, pnei1part, ierr)
recursive subroutine get_model_supported_species(model_handle, max_species, model_species, num_species, ier)
recursive subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, cutoff, neighObject)
recursive subroutine check_model_compatibility(compute_arguments_handle, forces_optional, particle_energy_supported, particle_energy_optional, model_is_compatible, ierr)
recursive subroutine compute_numer_deriv(partnum, dir, model_handle, compute_arguments_handle, DIM, N, coords, cutoff, cutpad, energy, do_update_list, coordsave, neighObject, deriv, deriv_err, ierr)
recursive subroutine update_neighborlist(DIM, N, coords, cutoff, cutpad, do_update_list, coordsave, neighObject, ierr)
recursive subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, coords, MiddlePartId)
program vc_forces_numer_deriv
real(c_double) recursive function dfridr(h, err)