33 use,
intrinsic :: iso_c_binding
38 recursive subroutine my_error(message)
40 character(len=*, kind=c_char),
intent(in) :: message
42 print *,
"* Error : ", trim(message)
48 character(len=*, kind=c_char),
intent(in) :: message
50 print *,
"* Warning : ", trim(message)
65 use,
intrinsic :: iso_c_binding
70 real(c_double) :: cutoff
71 integer(c_int) :: number_of_particles
72 type(c_ptr) :: neighbor_list_pointer
90 recursive subroutine get_neigh(data_object, number_of_neighbor_lists, cutoffs, &
91 neighbor_list_index, request, numnei, pnei1part, ierr) bind(c)
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
106 integer(c_int) numberofparticles
108 integer(c_int),
pointer :: neighborlist(:,:)
110 if (number_of_neighbor_lists > 1)
then 111 call my_warning(
"Model requires too many neighbor lists")
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])
121 if (cutoffs(neighbor_list_index) > neighobject%cutoff)
then 122 call my_warning(
"neighbor list cutoff too small for model cutoff")
127 if ( (request.gt.numberofparticles) .or. (request.lt.1))
then 129 call my_warning(
"Invalid part ID in get_neigh")
135 numnei = neighborlist(1,request)
138 pnei1part = c_loc(neighborlist(2,request))
159 forces_optional, model_is_compatible, ierr)
160 use,
intrinsic :: iso_c_binding
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
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
180 model_is_compatible = .false.
184 call kim_get_number_of_compute_argument_names(&
185 number_of_argument_names)
186 do i=1,number_of_argument_names
187 call kim_get_compute_argument_name(i, argument_name, &
193 call kim_get_argument_support_status( &
194 compute_arguments_handle, argument_name, support_status, ierr)
196 call my_warning(
"can't get argument support_status")
201 if (support_status == kim_support_status_required)
then 203 (argument_name == kim_compute_argument_name_partial_energy) .or. &
204 (argument_name == kim_compute_argument_name_partial_forces)))
then 205 call my_warning(
"unsupported required argument")
212 if ((argument_name == kim_compute_argument_name_partial_energy) .and. &
213 (support_status == kim_support_status_not_supported))
then 214 call my_warning(
"model does not support energy")
218 if (argument_name == kim_compute_argument_name_partial_forces)
then 219 if (support_status == kim_support_status_not_supported)
then 220 call my_warning(
"model does not support forces")
223 else if (support_status == kim_support_status_required)
then 224 forces_optional = .false.
225 else if (support_status == kim_support_status_optional)
then 226 forces_optional = .true.
228 call my_warning(
"unknown support_status for forces")
236 call kim_get_number_of_compute_callback_names( &
237 number_of_callback_names)
238 do i=1,number_of_callback_names
239 call kim_get_compute_callback_name(i, callback_name, &
245 call kim_get_callback_support_status( &
246 compute_arguments_handle, callback_name, support_status, ierr)
248 call my_warning(
"can't get call back support_status")
253 if (support_status == kim_support_status_required)
then 254 call my_warning(
"unsupported required call back")
261 model_is_compatible = .true.
273 model_species, num_species, ier)
274 use,
intrinsic :: iso_c_binding
278 type(kim_model_handle_type),
intent(in) :: model_handle
279 integer(c_int),
intent(in) :: max_species
280 type(kim_species_name_type),
intent(out) :: model_species(max_species)
281 integer(c_int),
intent(out) :: num_species
282 integer(c_int),
intent(out) :: ier
286 integer(c_int) total_num_species
287 type(kim_species_name_type) :: species_name
288 integer(c_int) species_is_supported
294 call kim_get_number_of_species_names(total_num_species)
296 if (total_num_species .gt. max_species)
return 299 do i=1,total_num_species
300 call kim_get_species_name(i, species_name, ier)
301 call kim_get_species_support_and_code(model_handle, species_name, &
302 species_is_supported, code, ier)
303 if ((ier == 0) .and. (species_is_supported .ne. 0))
then 304 num_species = num_species + 1
305 model_species(num_species) = species_name
315 do_update_list,coordsave, &
317 use,
intrinsic :: iso_c_binding
320 integer(c_int),
parameter :: cd = c_double
323 integer(c_int),
intent(in) :: DIM
324 integer(c_int),
intent(in) :: N
325 real(c_double),
intent(in) :: coords(dim,n)
326 real(c_double),
intent(in) :: cutoff
327 real(c_double),
intent(in) :: cutpad
328 logical,
intent(inout) :: do_update_list
329 real(c_double),
intent(inout) :: coordsave(dim,n)
331 integer(c_int),
intent(out) :: ierr
334 real(c_double) disp, disp1, disp2, cutrange, dispvec(dim)
343 if (.not.do_update_list)
then 350 dispvec(1:dim) = coords(1:dim,i) - coordsave(1:dim,i)
351 disp = sqrt( dot_product(dispvec,dispvec) )
352 if (disp >= disp1)
then 355 else if (disp >= disp2)
then 359 do_update_list = ( disp1 + disp2 > cutpad )
363 if (do_update_list)
then 366 coordsave(1:dim,1:n) = coords(1:dim,1:n)
369 cutrange = cutoff + cutpad
374 do_update_list = .false.
389 coords, cutoff, neighObject)
390 use,
intrinsic :: iso_c_binding
395 logical,
intent(in) :: half
396 integer(c_int),
intent(in) :: numberOfParticles
397 real(c_double),
dimension(3,numberOfParticles), &
399 real(c_double),
intent(in) :: cutoff
403 integer(c_int),
pointer :: neighborList(:,:)
404 integer(c_int) i, j, a
407 real(c_double) cutoff2
409 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
410 [numberofparticles+1, numberofparticles])
412 neighobject%cutoff = cutoff
416 do i=1,numberofparticles
418 do j=1,numberofparticles
419 dx(:) = coords(:, j) - coords(:, i)
420 r2 = dot_product(dx, dx)
421 if (r2.le.cutoff2)
then 423 if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) )
then 425 neighborlist(a,i) = j
430 neighborlist(1,i) = a-1
455 periodic, coords, MiddlePartId)
456 use,
intrinsic :: iso_c_binding
458 integer(c_int),
parameter :: cd = c_double
461 real(c_double),
intent(in) :: FCCspacing
462 integer(c_int),
intent(in) :: nCellsPerSide
463 logical,
intent(in) :: periodic
464 real(c_double),
intent(out) :: coords(3,*)
465 integer(c_int),
intent(out) :: MiddlePartId
469 real(c_double) FCCshifts(3,4)
470 real(c_double) latVec(3)
471 integer(c_int) a, i, j, k, m
475 fccshifts(1,1) = 0.0_cd
476 fccshifts(2,1) = 0.0_cd
477 fccshifts(3,1) = 0.0_cd
478 fccshifts(1,2) = 0.5_cd*fccspacing
479 fccshifts(2,2) = 0.5_cd*fccspacing
480 fccshifts(3,2) = 0.0_cd
481 fccshifts(1,3) = 0.5_cd*fccspacing
482 fccshifts(2,3) = 0.0_cd
483 fccshifts(3,3) = 0.5_cd*fccspacing
484 fccshifts(1,4) = 0.0_cd
485 fccshifts(2,4) = 0.5_cd*fccspacing
486 fccshifts(3,4) = 0.5_cd*fccspacing
491 latvec(1) = (i-1)*fccspacing
493 latvec(2) = (j-1)*fccspacing
495 latvec(3) = (k-1)*fccspacing
498 coords(:,a) = latvec + fccshifts(:,m)
499 if ((i.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
500 (k.eq.ncellsperside/2+1) .and. (m.eq.1))
then 501 coords(:,1) = latvec + fccshifts(:,m)
506 if (.not. periodic)
then 509 latvec(1) = ncellsperside*fccspacing
510 latvec(2) = (i-1)*fccspacing
511 latvec(3) = (j-1)*fccspacing
512 a = a+1; coords(:,a) = latvec
513 a = a+1; coords(:,a) = latvec + fccshifts(:,4)
515 latvec(1) = (i-1)*fccspacing
516 latvec(2) = ncellsperside*fccspacing
517 latvec(3) = (j-1)*fccspacing
518 a = a+1; coords(:,a) = latvec
519 a = a+1; coords(:,a) = latvec + fccshifts(:,3)
521 latvec(1) = (i-1)*fccspacing
522 latvec(2) = (j-1)*fccspacing
523 latvec(3) = ncellsperside*fccspacing
524 a = a+1; coords(:,a) = latvec
525 a = a+1; coords(:,a) = latvec + fccshifts(:,2)
528 if (.not. periodic)
then 530 latvec(1) = (i-1)*fccspacing
531 latvec(2) = ncellsperside*fccspacing
532 latvec(3) = ncellsperside*fccspacing
533 a = a+1; coords(:,a) = latvec
534 latvec(1) = ncellsperside*fccspacing
535 latvec(2) = (i-1)*fccspacing
536 latvec(3) = ncellsperside*fccspacing
537 a = a+1; coords(:,a) = latvec
538 latvec(1) = ncellsperside*fccspacing
539 latvec(2) = ncellsperside*fccspacing
540 latvec(3) = (i-1)*fccspacing
541 a = a+1; coords(:,a) = latvec
544 if (.not. periodic)
then 546 a = a+1; coords(:,a) = ncellsperside*fccspacing
554 compute_arguments_handle,DIM,N,coords,cutoff,cutpad,energy,do_update_list, &
555 coordsave,neighObject,deriv,deriv_err,ierr)
556 use,
intrinsic :: iso_c_binding
560 integer(c_int),
parameter :: cd = c_double
563 integer(c_int),
intent(in) :: partnum
564 integer(c_int),
intent(in) :: dir
565 type(kim_model_handle_type),
intent(in) :: model_handle
566 type(kim_compute_arguments_handle_type),
intent(in) :: compute_arguments_handle
567 integer(c_int),
intent(in) :: DIM
568 integer(c_int),
intent(in) :: N
569 real(c_double),
intent(inout) :: coords(dim,n)
570 real(c_double),
intent(in) :: cutoff
571 real(c_double),
intent(in) :: cutpad
572 real(c_double),
intent(inout) :: energy
573 logical,
intent(inout) :: do_update_list
574 real(c_double),
intent(inout) :: coordsave(dim,n)
576 real(c_double),
intent(out) :: deriv
577 real(c_double),
intent(out) :: deriv_err
578 integer(c_int),
intent(out) :: ierr
581 real(c_double),
parameter :: eps_init = 1.e-6_cd
582 integer(c_int),
parameter :: number_eps_levels = 15
583 real(c_double) eps, deriv_last, deriv_err_last
594 deriv_err_last = huge(1.0_cd)
595 do i=1,number_eps_levels
596 deriv =
dfridr(eps,deriv_err)
598 call my_error(
"compute_numer_deriv")
600 if (deriv_err>deriv_err_last)
then 602 deriv_err = deriv_err_last
607 deriv_err_last = deriv_err
631 real(c_double) recursive function dfridr(h,err)
635 real(c_double),
intent(inout) :: h
636 real(c_double),
intent(out) :: err
639 integer(c_int),
parameter :: NTAB=10
640 real(c_double),
parameter :: CON=1.4_cd
641 real(c_double),
parameter :: CON2=con*con
642 real(c_double),
parameter :: BIG=huge(1.0_cd)
643 real(c_double),
parameter :: SAFE=2.0_cd
646 real(c_double) errt,fac,hh,a(ntab,ntab),fp,fm,coordorig
650 if (abs(h).le.tiny(0.0_cd))
then 656 coordorig = coords(dir,partnum)
657 coords(dir,partnum) = coordorig + hh
659 do_update_list,coordsave, &
661 call kim_compute(model_handle, compute_arguments_handle, ierr)
663 call my_error(
"kim_api_model_compute")
666 coords(dir,partnum) = coordorig - hh
668 do_update_list,coordsave, &
670 call kim_compute(model_handle, compute_arguments_handle, ierr)
672 call my_error(
"kim_api_model_compute")
675 coords(dir,partnum) = coordorig
677 do_update_list,coordsave, &
679 a(1,1)=(fp-fm)/(2.0_cd*hh)
686 coords(dir,partnum) = coordorig + hh
688 do_update_list,coordsave, &
690 call kim_compute(model_handle, compute_arguments_handle, ierr)
692 call my_error(
"kim_api_model_compute")
695 coords(dir,partnum) = coordorig - hh
697 do_update_list,coordsave, &
699 call kim_compute(model_handle, compute_arguments_handle, ierr)
701 call my_error(
"kim_api_model_compute")
704 coords(dir,partnum) = coordorig
706 do_update_list,coordsave, &
708 a(1,i)=(fp-fm)/(2.0_cd*hh)
713 a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.0_cd)
717 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
718 if (errt.le.err)
then 723 if (abs(a(i,i)-a(i-1,i-1)).ge.safe*err)
return 748 use,
intrinsic :: iso_c_binding
753 integer(c_int),
parameter :: cd = c_double
755 integer(c_int),
parameter :: nCellsPerSide = 2
756 integer(c_int),
parameter :: DIM = 3
757 real(c_double),
parameter :: cutpad = 0.75_cd
758 integer(c_int),
parameter :: max_species = 200
759 real(c_double),
parameter :: eps_prec = epsilon(1.0_cd)
760 real(c_double) FCCspacing
762 integer(c_int),
parameter :: &
763 N = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
764 real(c_double),
allocatable :: forces_num(:,:)
765 real(c_double),
allocatable :: forces_num_err(:,:)
766 type(kim_species_name_type) :: model_species(max_species)
767 integer(c_int),
target :: num_species
768 character(len=5, kind=c_char) :: passfail
769 real(c_double) :: forcediff
770 real(c_double) :: forcediff_sumsq
771 real(c_double) :: weight
772 real(c_double) :: weight_sum
773 real(c_double) :: alpha
774 real(c_double) :: term
775 real(c_double) :: term_max
776 real(c_double),
allocatable :: cluster_coords(:,:)
777 real(c_double),
allocatable :: cluster_disps(:,:)
778 type(kim_species_name_type),
allocatable :: cluster_species(:)
779 integer(c_int) I,J,Imax,Jmax,species_code
780 integer(c_int) seed_size
781 integer(c_int),
allocatable :: seed(:)
787 integer(c_int),
allocatable,
target :: neighborList(:,:)
788 real(c_double),
allocatable :: coordsave(:,:)
789 logical do_update_list
794 character(len=256, kind=c_char) :: testname =
"vc_forces_numer_deriv" 795 character(len=256, kind=c_char) :: modelname
797 type(kim_model_handle_type) :: model_handle
798 type(kim_compute_arguments_handle_type) :: compute_arguments_handle
799 integer(c_int) ierr, ierr2
800 integer(c_int) species_is_supported
801 integer(c_int),
target :: numberOfParticles
802 integer(c_int),
target :: particleSpeciesCodes(n)
803 integer(c_int),
target :: particleContributing(n)
804 real(c_double) :: influence_distance
805 integer(c_int) :: number_of_neighbor_lists
806 real(c_double),
allocatable :: cutoffs(:)
807 integer(c_int),
allocatable :: &
808 model_will_not_request_neighbors_of_noncontributing_particles(:)
809 real(c_double) :: cutoff
810 real(c_double),
target :: energy
811 real(c_double),
target :: coords(3,n)
812 real(c_double),
target :: forces_kim(3,n)
813 real(c_double) :: forces(3,n)
814 integer(c_int) middleDum
815 logical forces_optional
816 logical model_is_compatible
817 integer(c_int) number_of_model_routine_names
818 type(kim_model_routine_name_type) model_routine_name
819 integer(c_int) present
820 integer(c_int) required
821 integer(c_int) requested_units_accepted
822 real(c_double) rnd, deriv, deriv_err
823 real(c_double),
pointer :: null_pointer
825 nullify(null_pointer)
827 numberofparticles = n
837 call random_seed(size=seed_size)
838 allocate(seed(seed_size))
840 call random_seed(put=seed)
844 print
'("Please enter a valid KIM model name: ")' 850 print *,
'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES' 853 print
'("This is Test : ",A)', trim(testname)
854 print
'("Results for KIM Model : ",A)', trim(modelname)
858 call kim_model_create(kim_numbering_one_based, &
860 kim_energy_unit_ev, &
862 kim_temperature_unit_k, &
865 requested_units_accepted, &
872 if (requested_units_accepted == 0)
then 873 call my_error(
"Must adapt to model units")
877 call kim_get_number_of_model_routine_names(number_of_model_routine_names)
878 do i=1,number_of_model_routine_names
879 call kim_get_model_routine_name(i, model_routine_name, ierr)
881 call my_error(
"kim_get_model_routine_name")
883 call kim_is_routine_present(model_handle, model_routine_name,
present, &
886 call my_error(
"kim_is_routine_present")
889 if ((
present == 1) .and. (required == 1))
then 890 if (.not. ((model_routine_name == kim_model_routine_name_create) &
891 .or. (model_routine_name == &
892 kim_model_routine_name_compute_arguments_create) &
893 .or. (model_routine_name == kim_model_routine_name_compute) &
894 .or. (model_routine_name == kim_model_routine_name_refresh) &
895 .or. (model_routine_name == &
896 kim_model_routine_name_compute_arguments_destroy) &
897 .or. (model_routine_name == kim_model_routine_name_destroy)))
then 899 call my_error(
"Unknown required ModelRoutineName found.")
905 call kim_compute_arguments_create(model_handle, &
906 compute_arguments_handle, ierr)
908 call my_error(
"kim_model_compute_arguments_create")
912 model_is_compatible, ierr)
914 call my_error(
"error checking compatibility")
916 if (.not. model_is_compatible)
then 917 call my_error(
"incompatibility reported by check_model_compatibility")
925 call my_error(
"Get_Model_Supported_Species")
929 allocate(cluster_coords(3,n),cluster_disps(3,n),cluster_species(n))
931 call random_number(rnd)
932 species_code = 1 + int(rnd*num_species)
933 cluster_species(i) = model_species(species_code)
939 cluster_coords, middledum)
944 call random_number(rnd)
945 cluster_disps(j,i) = 0.1_cd*(rnd-0.5_cd)
951 call kim_set_argument_pointer(compute_arguments_handle, &
952 kim_compute_argument_name_number_of_particles, numberofparticles, ierr2)
954 call kim_set_argument_pointer(compute_arguments_handle, &
955 kim_compute_argument_name_particle_species_codes, particlespeciescodes, &
958 call kim_set_argument_pointer(compute_arguments_handle, &
959 kim_compute_argument_name_particle_contributing, particlecontributing, &
962 call kim_set_argument_pointer(compute_arguments_handle, &
963 kim_compute_argument_name_coordinates, coords, ierr2)
965 call kim_set_argument_pointer(compute_arguments_handle, &
966 kim_compute_argument_name_partial_energy, energy, ierr2)
968 call kim_set_argument_pointer(compute_arguments_handle, &
969 kim_compute_argument_name_partial_forces, forces_kim, ierr2)
972 call my_error(
"set_argument_pointer")
978 allocate(neighborlist(n+1,n))
979 neighobject%neighbor_list_pointer = c_loc(neighborlist)
980 neighobject%number_of_particles = n
984 call kim_set_callback_pointer(compute_arguments_handle, &
985 kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
986 c_funloc(
get_neigh), c_loc(neighobject), ierr)
988 call my_error(
"set_callback_pointer")
991 call kim_get_influence_distance(model_handle, influence_distance)
992 call kim_get_number_of_neighbor_lists(model_handle, &
993 number_of_neighbor_lists)
994 allocate(cutoffs(number_of_neighbor_lists), &
995 model_will_not_request_neighbors_of_noncontributing_particles( &
996 number_of_neighbor_lists))
997 call kim_get_neighbor_list_values(model_handle, cutoffs, &
998 model_will_not_request_neighbors_of_noncontributing_particles, ierr)
1000 call my_error(
"get_neighbor_list_values")
1002 cutoff = maxval(cutoffs)
1005 fccspacing = 0.75_cd*cutoff
1008 cluster_coords(:,i) = fccspacing*cluster_coords(:,i)
1010 print
'("Using FCC lattice parameter: ",f12.5)', fccspacing
1013 call kim_get_species_support_and_code(model_handle, &
1014 cluster_species(i), species_is_supported, particlespeciescodes(i), ierr)
1017 call my_error(
"kim_api_get_species_code")
1020 particlecontributing(i) = 1
1023 coords(:,i) = cluster_coords(:,i) + cluster_disps(:,i)
1028 do_update_list = .true.
1029 allocate(coordsave(dim,n))
1031 do_update_list,coordsave, &
1034 call my_error(
"update_neighborlist")
1039 call kim_compute(model_handle, compute_arguments_handle, ierr)
1041 call my_error(
"kim_api_model_compute")
1051 print
'("Energy = ",ES25.15)', energy
1057 if (forces_optional)
then 1058 call kim_set_argument_pointer(compute_arguments_handle, &
1059 kim_compute_argument_name_partial_forces, null_pointer, ierr)
1061 call my_error(
"set_argument_pointer")
1067 allocate(forces_num(dim,n),forces_num_err(dim,n))
1071 dim,n,coords,cutoff, &
1072 cutpad, energy, do_update_list, &
1073 coordsave,neighobject,deriv,deriv_err,ierr)
1075 call my_error(
"compute_numer_deriv")
1077 forces_num(j,i) = -deriv
1078 forces_num_err(j,i) = deriv_err
1084 print
'(A6,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)',
"Part",
"Spec",
"Dir", &
1085 "Force_model",
"Force_numer",
"Force diff",
"pred error",
"weight", &
1087 forcediff_sumsq = 0.0_cd
1091 forcediff = abs(forces(j,i)-forces_num(j,i))
1092 if (forcediff<forces_num_err(j,i))
then 1097 weight = max(abs(forces_num(j,i)),eps_prec)/ &
1098 max(abs(forces_num_err(j,i)),eps_prec)
1099 term = weight*forcediff**2
1100 if (term.gt.term_max)
then 1105 forcediff_sumsq = forcediff_sumsq + term
1106 weight_sum = weight_sum + weight
1108 print
'(I6,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1109 i,particlespeciescodes(i),j,forces(j,i),forces_num(j,i), &
1110 forcediff,forces_num_err(j,i),weight,passfail
1112 print
'(14X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1113 j,forces(j,i),forces_num(j,i), &
1114 forcediff,forces_num_err(j,i),weight,passfail
1119 alpha = sqrt(forcediff_sumsq/weight_sum)/dble(dim*n)
1121 print
'("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," (units of force)")', &
1124 print
'(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,' // &
1125 ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
1126 imax,jmax,abs(forces(jmax,imax)-forces_num(jmax,imax)), &
1127 abs(forces(jmax,imax)-forces_num(jmax,imax))/abs(forces(jmax,imax))
1131 deallocate(forces_num)
1132 deallocate(forces_num_err)
1133 deallocate(neighborlist)
1134 deallocate(coordsave)
1135 deallocate(cutoffs, model_will_not_request_neighbors_of_noncontributing_particles)
1137 call kim_compute_arguments_destroy(model_handle, &
1138 compute_arguments_handle, ierr)
1140 call my_error(
"kim_model_compute_arguments_destroy")
1142 call kim_model_destroy(model_handle)
1147 print
'(120(''-''))' 1151 deallocate(cluster_coords,cluster_disps,cluster_species)
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)