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.
181 forces_optional = .false.
185 call kim_get_number_of_compute_argument_names(&
186 number_of_argument_names)
187 do i=1,number_of_argument_names
188 call kim_get_compute_argument_name(i, argument_name, &
194 call kim_get_argument_support_status( &
195 compute_arguments_handle, argument_name, support_status, ierr)
197 call my_warning(
"can't get argument support_status")
202 if (support_status == kim_support_status_required)
then 204 (argument_name == kim_compute_argument_name_partial_energy) .or. &
205 (argument_name == kim_compute_argument_name_partial_forces)))
then 206 call my_warning(
"unsupported required argument")
213 if ((argument_name == kim_compute_argument_name_partial_energy) .and. &
214 (support_status == kim_support_status_not_supported))
then 215 call my_warning(
"model does not support energy")
219 if (argument_name == kim_compute_argument_name_partial_forces)
then 220 if (support_status == kim_support_status_not_supported)
then 221 call my_warning(
"model does not support forces")
224 else if (support_status == kim_support_status_required)
then 225 forces_optional = .false.
226 else if (support_status == kim_support_status_optional)
then 227 forces_optional = .true.
229 call my_warning(
"unknown support_status for forces")
237 call kim_get_number_of_compute_callback_names( &
238 number_of_callback_names)
239 do i=1,number_of_callback_names
240 call kim_get_compute_callback_name(i, callback_name, &
246 call kim_get_callback_support_status( &
247 compute_arguments_handle, callback_name, support_status, ierr)
249 call my_warning(
"can't get call back support_status")
254 if (support_status == kim_support_status_required)
then 255 call my_warning(
"unsupported required call back")
262 model_is_compatible = .true.
274 model_species, num_species, ier)
275 use,
intrinsic :: iso_c_binding
279 type(kim_model_handle_type),
intent(in) :: model_handle
280 integer(c_int),
intent(in) :: max_species
281 type(kim_species_name_type),
intent(out) :: model_species(max_species)
282 integer(c_int),
intent(out) :: num_species
283 integer(c_int),
intent(out) :: ier
287 integer(c_int) total_num_species
288 type(kim_species_name_type) :: species_name
289 integer(c_int) species_is_supported
297 call kim_get_number_of_species_names(total_num_species)
299 if (total_num_species .gt. max_species)
return 302 do i=1,total_num_species
303 call kim_get_species_name(i, species_name, ier)
304 call kim_get_species_support_and_code(model_handle, species_name, &
305 species_is_supported, code, ier)
306 if ((ier == 0) .and. (species_is_supported .ne. 0))
then 307 num_species = num_species + 1
308 model_species(num_species) = species_name
318 do_update_list,coordsave, &
320 use,
intrinsic :: iso_c_binding
323 integer(c_int),
parameter :: cd = c_double
326 integer(c_int),
intent(in) :: DIM
327 integer(c_int),
intent(in) :: N
328 real(c_double),
intent(in) :: coords(dim,n)
329 real(c_double),
intent(in) :: cutoff
330 real(c_double),
intent(in) :: cutpad
331 logical,
intent(inout) :: do_update_list
332 real(c_double),
intent(inout) :: coordsave(dim,n)
334 integer(c_int),
intent(out) :: ierr
337 real(c_double) disp, disp1, disp2, cutrange, dispvec(dim)
346 if (.not.do_update_list)
then 353 dispvec(1:dim) = coords(1:dim,i) - coordsave(1:dim,i)
354 disp = sqrt( dot_product(dispvec,dispvec) )
355 if (disp >= disp1)
then 358 else if (disp >= disp2)
then 362 do_update_list = ( disp1 + disp2 > cutpad )
366 if (do_update_list)
then 369 coordsave(1:dim,1:n) = coords(1:dim,1:n)
372 cutrange = cutoff + cutpad
377 do_update_list = .false.
392 coords, cutoff, neighObject)
393 use,
intrinsic :: iso_c_binding
398 logical,
intent(in) :: half
399 integer(c_int),
intent(in) :: numberOfParticles
400 real(c_double),
dimension(3,numberOfParticles), &
402 real(c_double),
intent(in) :: cutoff
406 integer(c_int),
pointer :: neighborList(:,:)
407 integer(c_int) i, j, a
410 real(c_double) cutoff2
412 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
413 [numberofparticles+1, numberofparticles])
415 neighobject%cutoff = cutoff
419 do i=1,numberofparticles
421 do j=1,numberofparticles
422 dx(:) = coords(:, j) - coords(:, i)
423 r2 = dot_product(dx, dx)
424 if (r2.le.cutoff2)
then 426 if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) )
then 428 neighborlist(a,i) = j
433 neighborlist(1,i) = a-1
458 periodic, coords, MiddlePartId)
459 use,
intrinsic :: iso_c_binding
461 integer(c_int),
parameter :: cd = c_double
464 real(c_double),
intent(in) :: FCCspacing
465 integer(c_int),
intent(in) :: nCellsPerSide
466 logical,
intent(in) :: periodic
467 real(c_double),
intent(out) :: coords(3,*)
468 integer(c_int),
intent(out) :: MiddlePartId
472 real(c_double) FCCshifts(3,4)
473 real(c_double) latVec(3)
474 integer(c_int) a, i, j, k, m
478 fccshifts(1,1) = 0.0_cd
479 fccshifts(2,1) = 0.0_cd
480 fccshifts(3,1) = 0.0_cd
481 fccshifts(1,2) = 0.5_cd*fccspacing
482 fccshifts(2,2) = 0.5_cd*fccspacing
483 fccshifts(3,2) = 0.0_cd
484 fccshifts(1,3) = 0.5_cd*fccspacing
485 fccshifts(2,3) = 0.0_cd
486 fccshifts(3,3) = 0.5_cd*fccspacing
487 fccshifts(1,4) = 0.0_cd
488 fccshifts(2,4) = 0.5_cd*fccspacing
489 fccshifts(3,4) = 0.5_cd*fccspacing
494 latvec(1) = (i-1)*fccspacing
496 latvec(2) = (j-1)*fccspacing
498 latvec(3) = (k-1)*fccspacing
501 coords(:,a) = latvec + fccshifts(:,m)
502 if ((i.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
503 (k.eq.ncellsperside/2+1) .and. (m.eq.1))
then 504 coords(:,1) = latvec + fccshifts(:,m)
509 if (.not. periodic)
then 512 latvec(1) = ncellsperside*fccspacing
513 latvec(2) = (i-1)*fccspacing
514 latvec(3) = (j-1)*fccspacing
515 a = a+1; coords(:,a) = latvec
516 a = a+1; coords(:,a) = latvec + fccshifts(:,4)
518 latvec(1) = (i-1)*fccspacing
519 latvec(2) = ncellsperside*fccspacing
520 latvec(3) = (j-1)*fccspacing
521 a = a+1; coords(:,a) = latvec
522 a = a+1; coords(:,a) = latvec + fccshifts(:,3)
524 latvec(1) = (i-1)*fccspacing
525 latvec(2) = (j-1)*fccspacing
526 latvec(3) = ncellsperside*fccspacing
527 a = a+1; coords(:,a) = latvec
528 a = a+1; coords(:,a) = latvec + fccshifts(:,2)
531 if (.not. periodic)
then 533 latvec(1) = (i-1)*fccspacing
534 latvec(2) = ncellsperside*fccspacing
535 latvec(3) = ncellsperside*fccspacing
536 a = a+1; coords(:,a) = latvec
537 latvec(1) = ncellsperside*fccspacing
538 latvec(2) = (i-1)*fccspacing
539 latvec(3) = ncellsperside*fccspacing
540 a = a+1; coords(:,a) = latvec
541 latvec(1) = ncellsperside*fccspacing
542 latvec(2) = ncellsperside*fccspacing
543 latvec(3) = (i-1)*fccspacing
544 a = a+1; coords(:,a) = latvec
547 if (.not. periodic)
then 549 a = a+1; coords(:,a) = ncellsperside*fccspacing
557 compute_arguments_handle,DIM,N,coords,cutoff,cutpad,energy,do_update_list, &
558 coordsave,neighObject,deriv,deriv_err,ierr)
559 use,
intrinsic :: iso_c_binding
563 integer(c_int),
parameter :: cd = c_double
566 integer(c_int),
intent(in) :: partnum
567 integer(c_int),
intent(in) :: dir
568 type(kim_model_handle_type),
intent(in) :: model_handle
569 type(kim_compute_arguments_handle_type),
intent(in) :: compute_arguments_handle
570 integer(c_int),
intent(in) :: DIM
571 integer(c_int),
intent(in) :: N
572 real(c_double),
intent(inout) :: coords(dim,n)
573 real(c_double),
intent(in) :: cutoff
574 real(c_double),
intent(in) :: cutpad
575 real(c_double),
intent(inout) :: energy
576 logical,
intent(inout) :: do_update_list
577 real(c_double),
intent(inout) :: coordsave(dim,n)
579 real(c_double),
intent(out) :: deriv
580 real(c_double),
intent(out) :: deriv_err
581 integer(c_int),
intent(out) :: ierr
584 real(c_double),
parameter :: eps_init = 1.e-6_cd
585 integer(c_int),
parameter :: number_eps_levels = 15
586 real(c_double) eps, deriv_last, deriv_err_last
597 deriv_err_last = huge(1.0_cd)
598 do i=1,number_eps_levels
599 deriv =
dfridr(eps,deriv_err)
601 call my_error(
"compute_numer_deriv")
603 if (deriv_err>deriv_err_last)
then 605 deriv_err = deriv_err_last
610 deriv_err_last = deriv_err
634 real(c_double) recursive function dfridr(h,err)
638 real(c_double),
intent(inout) :: h
639 real(c_double),
intent(out) :: err
642 integer(c_int),
parameter :: NTAB=10
643 real(c_double),
parameter :: CON=1.4_cd
644 real(c_double),
parameter :: CON2=con*con
645 real(c_double),
parameter :: BIG=huge(1.0_cd)
646 real(c_double),
parameter :: SAFE=2.0_cd
649 real(c_double) errt,fac,hh,a(ntab,ntab),fp,fm,coordorig
654 if (abs(h).le.tiny(0.0_cd))
then 660 coordorig = coords(dir,partnum)
661 coords(dir,partnum) = coordorig + hh
663 do_update_list,coordsave, &
665 call kim_compute(model_handle, compute_arguments_handle, ierr)
667 call my_error(
"kim_api_model_compute")
670 coords(dir,partnum) = coordorig - hh
672 do_update_list,coordsave, &
674 call kim_compute(model_handle, compute_arguments_handle, ierr)
676 call my_error(
"kim_api_model_compute")
679 coords(dir,partnum) = coordorig
681 do_update_list,coordsave, &
683 a(1,1)=(fp-fm)/(2.0_cd*hh)
689 coords(dir,partnum) = coordorig + hh
691 do_update_list,coordsave, &
693 call kim_compute(model_handle, compute_arguments_handle, ierr)
695 call my_error(
"kim_api_model_compute")
698 coords(dir,partnum) = coordorig - hh
700 do_update_list,coordsave, &
702 call kim_compute(model_handle, compute_arguments_handle, ierr)
704 call my_error(
"kim_api_model_compute")
707 coords(dir,partnum) = coordorig
709 do_update_list,coordsave, &
711 a(1,i)=(fp-fm)/(2.0_cd*hh)
716 a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.0_cd)
720 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
721 if (errt.le.err)
then 726 if (abs(a(i,i)-a(i-1,i-1)).ge.safe*err)
return 751 use,
intrinsic :: iso_c_binding
756 integer(c_int),
parameter :: cd = c_double
758 integer(c_int),
parameter :: nCellsPerSide = 2
759 integer(c_int),
parameter :: DIM = 3
760 real(c_double),
parameter :: cutpad = 0.75_cd
761 integer(c_int),
parameter :: max_species = 200
762 real(c_double),
parameter :: eps_prec = epsilon(1.0_cd)
763 real(c_double) FCCspacing
765 integer(c_int),
parameter :: &
766 N = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
767 real(c_double),
allocatable :: forces_num(:,:)
768 real(c_double),
allocatable :: forces_num_err(:,:)
769 type(kim_species_name_type) :: model_species(max_species)
770 integer(c_int),
target :: num_species
771 character(len=5, kind=c_char) :: passfail
772 real(c_double) :: forcediff
773 real(c_double) :: forcediff_sumsq
774 real(c_double) :: weight
775 real(c_double) :: weight_sum
776 real(c_double) :: alpha
777 real(c_double) :: term
778 real(c_double) :: term_max
779 real(c_double),
allocatable :: cluster_coords(:,:)
780 real(c_double),
allocatable :: cluster_disps(:,:)
781 type(kim_species_name_type),
allocatable :: cluster_species(:)
782 integer(c_int) I,J,Imax,Jmax,species_code
783 integer(c_int) seed_size
784 integer(c_int),
allocatable :: seed(:)
790 integer(c_int),
allocatable,
target :: neighborList(:,:)
791 real(c_double),
allocatable :: coordsave(:,:)
792 logical do_update_list
797 character(len=256, kind=c_char) :: testname =
"vc_forces_numer_deriv" 798 character(len=256, kind=c_char) :: modelname
800 type(kim_model_handle_type) :: model_handle
801 type(kim_compute_arguments_handle_type) :: compute_arguments_handle
802 integer(c_int) ierr, ierr2
803 integer(c_int) species_is_supported
804 integer(c_int),
target :: numberOfParticles
805 integer(c_int),
target :: particleSpeciesCodes(n)
806 integer(c_int),
target :: particleContributing(n)
807 real(c_double) :: influence_distance
808 integer(c_int) :: number_of_neighbor_lists
809 real(c_double),
allocatable :: cutoffs(:)
810 integer(c_int),
allocatable :: &
811 model_will_not_request_neighbors_of_noncontributing_particles(:)
812 real(c_double) :: cutoff
813 real(c_double),
target :: energy
814 real(c_double),
target :: coords(3,n)
815 real(c_double),
target :: forces_kim(3,n)
816 real(c_double) :: forces(3,n)
817 integer(c_int) middleDum
818 logical forces_optional
819 logical model_is_compatible
820 integer(c_int) number_of_model_routine_names
821 type(kim_model_routine_name_type) model_routine_name
822 integer(c_int) present
823 integer(c_int) required
824 integer(c_int) requested_units_accepted
825 real(c_double) rnd, deriv, deriv_err
826 real(c_double),
pointer :: null_pointer
828 nullify(null_pointer)
830 numberofparticles = n
840 call random_seed(size=seed_size)
841 allocate(seed(seed_size))
843 call random_seed(put=seed)
847 print
'("Please enter a valid KIM model name: ")' 853 print *,
'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES' 856 print
'("This is Test : ",A)', trim(testname)
857 print
'("Results for KIM Model : ",A)', trim(modelname)
861 call kim_model_create(kim_numbering_one_based, &
863 kim_energy_unit_ev, &
865 kim_temperature_unit_k, &
868 requested_units_accepted, &
875 if (requested_units_accepted == 0)
then 876 call my_error(
"Must adapt to model units")
880 call kim_get_number_of_model_routine_names(number_of_model_routine_names)
881 do i=1,number_of_model_routine_names
882 call kim_get_model_routine_name(i, model_routine_name, ierr)
884 call my_error(
"kim_get_model_routine_name")
886 call kim_is_routine_present(model_handle, model_routine_name,
present, &
889 call my_error(
"kim_is_routine_present")
892 if ((
present == 1) .and. (required == 1))
then 893 if (.not. ((model_routine_name == kim_model_routine_name_create) &
894 .or. (model_routine_name == &
895 kim_model_routine_name_compute_arguments_create) &
896 .or. (model_routine_name == kim_model_routine_name_compute) &
897 .or. (model_routine_name == kim_model_routine_name_refresh) &
898 .or. (model_routine_name == &
899 kim_model_routine_name_compute_arguments_destroy) &
900 .or. (model_routine_name == kim_model_routine_name_destroy)))
then 902 call my_error(
"Unknown required ModelRoutineName found.")
908 call kim_compute_arguments_create(model_handle, &
909 compute_arguments_handle, ierr)
911 call my_error(
"kim_model_compute_arguments_create")
915 model_is_compatible, ierr)
917 call my_error(
"error checking compatibility")
919 if (.not. model_is_compatible)
then 920 call my_error(
"incompatibility reported by check_model_compatibility")
928 call my_error(
"Get_Model_Supported_Species")
932 allocate(cluster_coords(3,n),cluster_disps(3,n),cluster_species(n))
934 call random_number(rnd)
935 species_code = 1 + int(rnd*num_species)
936 cluster_species(i) = model_species(species_code)
942 cluster_coords, middledum)
947 call random_number(rnd)
948 cluster_disps(j,i) = 0.1_cd*(rnd-0.5_cd)
954 call kim_set_argument_pointer(compute_arguments_handle, &
955 kim_compute_argument_name_number_of_particles, numberofparticles, ierr2)
957 call kim_set_argument_pointer(compute_arguments_handle, &
958 kim_compute_argument_name_particle_species_codes, particlespeciescodes, &
961 call kim_set_argument_pointer(compute_arguments_handle, &
962 kim_compute_argument_name_particle_contributing, particlecontributing, &
965 call kim_set_argument_pointer(compute_arguments_handle, &
966 kim_compute_argument_name_coordinates, coords, ierr2)
968 call kim_set_argument_pointer(compute_arguments_handle, &
969 kim_compute_argument_name_partial_energy, energy, ierr2)
971 call kim_set_argument_pointer(compute_arguments_handle, &
972 kim_compute_argument_name_partial_forces, forces_kim, ierr2)
975 call my_error(
"set_argument_pointer")
981 allocate(neighborlist(n+1,n))
982 neighobject%neighbor_list_pointer = c_loc(neighborlist)
983 neighobject%number_of_particles = n
987 call kim_set_callback_pointer(compute_arguments_handle, &
988 kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
989 c_funloc(
get_neigh), c_loc(neighobject), ierr)
991 call my_error(
"set_callback_pointer")
994 call kim_get_influence_distance(model_handle, influence_distance)
995 call kim_get_number_of_neighbor_lists(model_handle, &
996 number_of_neighbor_lists)
997 allocate(cutoffs(number_of_neighbor_lists), &
998 model_will_not_request_neighbors_of_noncontributing_particles( &
999 number_of_neighbor_lists))
1000 call kim_get_neighbor_list_values(model_handle, cutoffs, &
1001 model_will_not_request_neighbors_of_noncontributing_particles, ierr)
1003 call my_error(
"get_neighbor_list_values")
1005 cutoff = maxval(cutoffs)
1008 fccspacing = 0.75_cd*cutoff
1011 cluster_coords(:,i) = fccspacing*cluster_coords(:,i)
1013 print
'("Using FCC lattice parameter: ",f12.5)', fccspacing
1016 call kim_get_species_support_and_code(model_handle, &
1017 cluster_species(i), species_is_supported, particlespeciescodes(i), ierr)
1020 call my_error(
"kim_api_get_species_code")
1023 particlecontributing(i) = 1
1026 coords(:,i) = cluster_coords(:,i) + cluster_disps(:,i)
1031 do_update_list = .true.
1032 allocate(coordsave(dim,n))
1034 do_update_list,coordsave, &
1037 call my_error(
"update_neighborlist")
1042 call kim_compute(model_handle, compute_arguments_handle, ierr)
1044 call my_error(
"kim_api_model_compute")
1054 print
'("Energy = ",ES25.15)', energy
1060 if (forces_optional)
then 1061 call kim_set_argument_pointer(compute_arguments_handle, &
1062 kim_compute_argument_name_partial_forces, null_pointer, ierr)
1064 call my_error(
"set_argument_pointer")
1070 allocate(forces_num(dim,n),forces_num_err(dim,n))
1074 dim,n,coords,cutoff, &
1075 cutpad, energy, do_update_list, &
1076 coordsave,neighobject,deriv,deriv_err,ierr)
1078 call my_error(
"compute_numer_deriv")
1080 forces_num(j,i) = -deriv
1081 forces_num_err(j,i) = deriv_err
1087 print
'(A6,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)',
"Part",
"Spec",
"Dir", &
1088 "Force_model",
"Force_numer",
"Force diff",
"pred error",
"weight", &
1090 forcediff_sumsq = 0.0_cd
1094 forcediff = abs(forces(j,i)-forces_num(j,i))
1095 if (forcediff<forces_num_err(j,i))
then 1100 weight = max(abs(forces_num(j,i)),eps_prec)/ &
1101 max(abs(forces_num_err(j,i)),eps_prec)
1102 term = weight*forcediff**2
1103 if (term.gt.term_max)
then 1108 forcediff_sumsq = forcediff_sumsq + term
1109 weight_sum = weight_sum + weight
1111 print
'(I6,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1112 i,particlespeciescodes(i),j,forces(j,i),forces_num(j,i), &
1113 forcediff,forces_num_err(j,i),weight,passfail
1115 print
'(14X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1116 j,forces(j,i),forces_num(j,i), &
1117 forcediff,forces_num_err(j,i),weight,passfail
1122 alpha = sqrt(forcediff_sumsq/weight_sum)/dble(dim*n)
1124 print
'("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," (units of force)")', &
1127 print
'(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,' // &
1128 ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
1129 imax,jmax,abs(forces(jmax,imax)-forces_num(jmax,imax)), &
1130 abs(forces(jmax,imax)-forces_num(jmax,imax))/abs(forces(jmax,imax))
1134 deallocate(forces_num)
1135 deallocate(forces_num_err)
1136 deallocate(neighborlist)
1137 deallocate(coordsave)
1138 deallocate(cutoffs, model_will_not_request_neighbors_of_noncontributing_particles)
1140 call kim_compute_arguments_destroy(model_handle, &
1141 compute_arguments_handle, ierr)
1143 call my_error(
"kim_model_compute_arguments_destroy")
1145 call kim_model_destroy(model_handle)
1150 print
'(120(''-''))' 1154 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)