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