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