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