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_neighbor_lists, 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_neighbor_lists
105  real(c_double), intent(in) :: cutoffs(number_of_neighbor_lists)
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_neighbor_lists /= 1) then
120  call my_warning("invalid number of neighbor_lists", __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  integer(c_int) seed_size
808  integer(c_int), allocatable :: seed(:)
809 
810  !
811  ! neighbor list
812  !
813  type(neighobject_type), target :: neighobject
814  real(c_double), allocatable :: coordsave(:,:)
815  logical do_update_list
816 
817  !
818  ! KIM variables
819  !
820  character(len=256, kind=c_char) :: testname = "vc_forces_numer_deriv"
821  character(len=256, kind=c_char) :: modelname
822 
823  type(kim_model_handle_type) :: model_handle
824  type(kim_compute_arguments_handle_type) :: compute_arguments_handle
825  integer(c_int) ierr, ierr2
826  integer(c_int) species_is_supported
827  integer(c_int), target :: numberofparticles
828  integer(c_int), target :: particlespeciescodes(n)
829  integer(c_int), target :: particlecontributing(n)
830  real(c_double) :: influence_distance
831  integer(c_int) :: number_of_neighbor_lists
832  real(c_double) :: cutoffs(1)
833  integer(c_int) :: padding_neighbor_hints(1)
834  integer(c_int) :: half_list_hints(1)
835  real(c_double) :: cutoff
836  real(c_double), target :: energy
837  real(c_double), target :: coords(3,n)
838  real(c_double), target :: forces_kim(3,n)
839  real(c_double) :: forces(3,n)
840  integer(c_int) middledum
841  logical forces_optional
842  logical model_is_compatible
843  integer(c_int) requested_units_accepted
844  real(c_double) rnd, deriv, deriv_err
845  real(c_double), pointer :: null_pointer
846 
847  nullify(null_pointer)
848 
849  numberofparticles = n
850 
851  term_max = 0.0_cd ! initialize
852 
853  ! Initialize error flag
854  ierr = 0
855 
856  ! Initialize seed for random number generator
857  !
858  ! NOTE: Here, we set the seed to a fixed value for reproducibility
859  call random_seed(size=seed_size) ! Get seed size for this CPU
860  allocate(seed(seed_size))
861  seed(:) = 13
862  call random_seed(put=seed)
863  deallocate(seed)
864 
865  ! Get KIM Model name to use
866  print '("Please enter a valid KIM model name: ")'
867  read(*,*) modelname
868 
869  ! Print output header
870  !
871  print *
872  print *,'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES'
873  print *
874  print '(120(''-''))'
875  print '("This is Test : ",A)', trim(testname)
876  print '("Results for KIM Model : ",A)', trim(modelname)
877 
878  ! Create empty KIM object
879  !
880  call kim_model_create(kim_numbering_one_based, &
881  kim_length_unit_a, &
882  kim_energy_unit_ev, &
883  kim_charge_unit_e, &
884  kim_temperature_unit_k, &
885  kim_time_unit_ps, &
886  trim(modelname), &
887  requested_units_accepted, &
888  model_handle, ierr)
889  if (ierr /= 0) then
890  call my_error("kim_api_create", __line__, __file__)
891  endif
892 
893  ! check that we are compatible
894  if (requested_units_accepted == 0) then
895  call my_error("Must adapt to model units", __line__, __file__)
896  end if
897 
898  ! create compute_arguments object
899  call kim_model_compute_arguments_create(model_handle, &
900  compute_arguments_handle, ierr)
901  if (ierr /= 0) then
902  call my_error("kim_model_compute_arguments_create", __line__, __file__)
903  endif
904 
905  call check_model_compatibility(compute_arguments_handle, forces_optional, &
906  model_is_compatible, ierr)
907  if (ierr /= 0) then
908  call my_error("error checking compatibility", __line__, __file__)
909  end if
910  if (.not. model_is_compatible) then
911  call my_error("incompatibility reported by check_model_compatibility", &
912  __line__, __file__)
913  end if
914 
915  ! Get list of particle species supported by the model
916  !
917  call get_model_supported_species(model_handle, max_species, model_species, &
918  num_species, ierr)
919  if (ierr /= 0) then
920  call my_error("Get_Model_Supported_Species", __line__, __file__)
921  endif
922  ! Setup random cluster
923  !
924  allocate(cluster_coords(3,n),cluster_disps(3,n),cluster_species(n))
925  do i=1,n
926  call random_number(rnd) ! return random number between 0 and 1
927  species_code = 1 + int(rnd*num_species)
928  cluster_species(i) = model_species(species_code)
929  enddo
930  fccspacing = 1.0_cd ! initially generate an FCC cluster with lattice
931  ! spacing equal to one. This is scaled below based
932  ! on the cutoff radius of the model.
933  call create_fcc_configuration(fccspacing, ncellsperside, .false., &
934  cluster_coords, middledum)
935  ! Generate random displacements for all particles
936  !
937  do i=1,n
938  do j=1,dim
939  call random_number(rnd) ! return random number between 0 and 1
940  cluster_disps(j,i) = 0.1_cd*(rnd-0.5_cd)
941  enddo
942  enddo
943 
944  ! register memory with the KIM system
945  ierr = 0
946  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
947  kim_compute_argument_name_number_of_particles, numberofparticles, ierr2)
948  ierr = ierr + ierr2
949  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
950  kim_compute_argument_name_particle_species_codes, particlespeciescodes, &
951  ierr2)
952  ierr = ierr + ierr2
953  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
954  kim_compute_argument_name_particle_contributing, particlecontributing, &
955  ierr2)
956  ierr = ierr + ierr2
957  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
958  kim_compute_argument_name_coordinates, coords, ierr2)
959  ierr = ierr + ierr2
960  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
961  kim_compute_argument_name_partial_energy, energy, ierr2)
962  ierr = ierr + ierr2
963  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
964  kim_compute_argument_name_partial_forces, forces_kim, ierr2)
965  ierr = ierr + ierr2
966  if (ierr /= 0) then
967  call my_error("set_argument_pointer", __line__, __file__)
968  endif
969 
970  ! Allocate storage for neighbor lists and
971  ! store pointers to neighbor list object and access function
972  !
973  allocate(neighobject%neighborList(n+1,n))
974  neighobject%number_of_particles = n
975 
976  ! Set pointer in KIM object to neighbor list routine and object
977  !
978  call kim_compute_arguments_set_callback_pointer(compute_arguments_handle, &
979  kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
980  c_funloc(get_neigh), c_loc(neighobject), ierr)
981  if (ierr /= 0) then
982  call my_error("set_callback_pointer", __line__, __file__)
983  end if
984 
985  call kim_model_get_influence_distance(model_handle, influence_distance)
986  call kim_model_get_number_of_neighbor_lists(model_handle, &
987  number_of_neighbor_lists)
988  if (number_of_neighbor_lists /= 1) then
989  call my_error("too many neighbor lists", __line__, __file__)
990  endif
991  call kim_model_get_neighbor_list_values(model_handle, cutoffs, &
992  padding_neighbor_hints, half_list_hints, ierr)
993  if (ierr /= 0) then
994  call my_error("get_neighbor_list_values", __line__, __file__)
995  end if
996  cutoff = cutoffs(1)
997 
998  ! Scale reference FCC configuration based on cutoff radius.
999  fccspacing = 0.75_cd*cutoff ! set the FCC spacing to a fraction
1000  ! of the cutoff radius
1001  do i=1,n
1002  cluster_coords(:,i) = fccspacing*cluster_coords(:,i)
1003  enddo
1004  print '("Using FCC lattice parameter: ",f12.5)', fccspacing
1005 
1006  do i=1,n
1007  call kim_model_get_species_support_and_code(model_handle, &
1008  cluster_species(i), species_is_supported, particlespeciescodes(i), ierr)
1009  enddo
1010  if (ierr /= 0) then
1011  call my_error("kim_api_get_species_code", __line__, __file__)
1012  endif
1013  do i=1,n
1014  particlecontributing(i) = 1 ! every particle contributes
1015  enddo
1016  do i=1,n
1017  coords(:,i) = cluster_coords(:,i) + cluster_disps(:,i)
1018  enddo
1019 
1020  ! Compute neighbor lists
1021  !
1022  do_update_list = .true.
1023  allocate(coordsave(dim,n))
1024  call update_neighborlist(dim,n,coords,cutoff,cutpad, &
1025  do_update_list,coordsave, &
1026  neighobject,ierr)
1027  if (ierr /= 0) then
1028  call my_error("update_neighborlist", __line__, __file__)
1029  endif
1030 
1031  ! Call model compute to get forces (gradient)
1032  !
1033  call kim_model_compute(model_handle, compute_arguments_handle, ierr)
1034  if (ierr /= 0) then
1035  call my_error("kim_api_model_compute", __line__, __file__)
1036  endif
1037 
1038  ! Copy forces in case model will overwrite forces_kim below
1039  !
1040  forces = forces_kim
1041 
1042  ! Print results to screen
1043  !
1044  print '(41(''=''))'
1045  print '("Energy = ",ES25.15)', energy
1046  print '(41(''=''))'
1047  print *
1048 
1049  ! Turn off force computation, if possible
1050  !
1051  if (forces_optional) then
1052  call kim_compute_arguments_set_argument_pointer(compute_arguments_handle, &
1053  kim_compute_argument_name_partial_forces, null_pointer, ierr)
1054  if (ierr /= 0) then
1055  call my_error("set_argument_pointer", __line__, __file__)
1056  endif
1057  endif
1058 
1059  ! Compute gradient using numerical differentiation
1060  !
1061  allocate(forces_num(dim,n),forces_num_err(dim,n))
1062  do i=1,n
1063  do j=1,dim
1064  call compute_numer_deriv(i,j,model_handle,compute_arguments_handle,&
1065  dim,n,coords,cutoff, &
1066  cutpad, energy, do_update_list, &
1067  coordsave,neighobject,deriv,deriv_err,ierr)
1068  if (ierr /= 0) then
1069  call my_error("compute_numer_deriv", __line__, __file__)
1070  endif
1071  forces_num(j,i) = -deriv
1072  forces_num_err(j,i) = deriv_err
1073  enddo
1074  enddo
1075 
1076  ! Continue printing results to screen
1077  !
1078  print '(A6,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)',"Part","Spec","Dir", &
1079  "Force_model", "Force_numer", "Force diff", "pred error", "weight", &
1080  "stat"
1081  forcediff_sumsq = 0.0_cd
1082  weight_sum = 0.0_cd
1083  do i=1,n
1084  do j=1,dim
1085  forcediff = abs(forces(j,i)-forces_num(j,i))
1086  if (forcediff<forces_num_err(j,i)) then
1087  passfail = "ideal"
1088  else
1089  passfail = " "
1090  endif
1091  weight = max(abs(forces_num(j,i)),eps_prec)/ &
1092  max(abs(forces_num_err(j,i)),eps_prec)
1093  term = weight*forcediff**2
1094  if (term.gt.term_max) then
1095  term_max = term
1096  imax = i
1097  jmax = j
1098  endif
1099  forcediff_sumsq = forcediff_sumsq + term
1100  weight_sum = weight_sum + weight
1101  if (j.eq.1) then
1102  print '(I6,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1103  i,particlespeciescodes(i),j,forces(j,i),forces_num(j,i), &
1104  forcediff,forces_num_err(j,i),weight,passfail
1105  else
1106  print '(14X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1107  j,forces(j,i),forces_num(j,i), &
1108  forcediff,forces_num_err(j,i),weight,passfail
1109  endif
1110  enddo
1111  print *
1112  enddo
1113  alpha = sqrt(forcediff_sumsq/weight_sum)/dble(dim*n)
1114  print *
1115  print '("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," (units of force)")', &
1116  alpha
1117  print *
1118  print '(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,' // &
1119  ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
1120  imax,jmax,abs(forces(jmax,imax)-forces_num(jmax,imax)), &
1121  abs(forces(jmax,imax)-forces_num(jmax,imax))/abs(forces(jmax,imax))
1122 
1123  ! Free temporary storage
1124  !
1125  deallocate(forces_num)
1126  deallocate(forces_num_err)
1127  deallocate(neighobject%neighborList)
1128  deallocate(coordsave)
1129 
1130  call kim_model_compute_arguments_destroy(model_handle, &
1131  compute_arguments_handle, ierr)
1132  if (ierr /= 0) then
1133  call my_error("kim_model_compute_arguments_destroy", __line__, __file__)
1134  endif
1135  call kim_model_destroy(model_handle)
1136 
1137  ! Print output footer
1138  !
1139  print *
1140  print '(120(''-''))'
1141 
1142  ! Free cluster storage
1143  !
1144  deallocate(cluster_coords,cluster_disps,cluster_species)
1145 
1146  stop
1147 
1148 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)