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(number_of_neighbor_lists)
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 call c_f_pointer(data_object, neighobject)
111 numberofparticles = neighobject%number_of_particles
112 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
113 [numberofparticles+1, numberofparticles])
115 if (cutoffs(neighbor_list_index) > neighobject%cutoff)
then 116 call my_warning(
"neighbor list cutoff too small for model cutoff")
121 if ( (request.gt.numberofparticles) .or. (request.lt.1))
then 123 call my_warning(
"Invalid part ID in get_neigh")
129 numnei = neighborlist(1,request)
132 pnei1part = c_loc(neighborlist(2,request))
153 forces_optional, model_is_compatible, ierr)
154 use,
intrinsic :: iso_c_binding
159 type(kim_compute_arguments_handle_type),
intent(in) :: compute_arguments_handle
160 logical,
intent(out) :: forces_optional
161 logical,
intent(out) :: model_is_compatible
162 integer(c_int),
intent(out) :: ierr
166 integer(c_int) number_of_argument_names
167 integer(c_int) number_of_callback_names
168 type(kim_compute_argument_name_type) argument_name
169 type(kim_support_status_type) support_status
170 type(kim_compute_callback_name_type) callback_name
174 model_is_compatible = .false.
178 call kim_get_number_of_compute_argument_names(&
179 number_of_argument_names)
180 do i=1,number_of_argument_names
181 call kim_get_compute_argument_name(i, argument_name, &
187 call kim_get_argument_support_status( &
188 compute_arguments_handle, argument_name, support_status, ierr)
190 call my_warning(
"can't get argument support_status")
195 if (support_status == kim_support_status_required)
then 197 (argument_name == kim_compute_argument_name_partial_energy) .or. &
198 (argument_name == kim_compute_argument_name_partial_forces)))
then 199 call my_warning(
"unsupported required argument")
206 if ((argument_name == kim_compute_argument_name_partial_energy) .and. &
207 (support_status == kim_support_status_not_supported))
then 208 call my_warning(
"model does not support energy")
212 if (argument_name == kim_compute_argument_name_partial_forces)
then 213 if (support_status == kim_support_status_not_supported)
then 214 call my_warning(
"model does not support forces")
217 else if (support_status == kim_support_status_required)
then 218 forces_optional = .false.
219 else if (support_status == kim_support_status_optional)
then 220 forces_optional = .true.
222 call my_warning(
"unknown support_status for forces")
230 call kim_get_number_of_compute_callback_names( &
231 number_of_callback_names)
232 do i=1,number_of_callback_names
233 call kim_get_compute_callback_name(i, callback_name, &
239 call kim_get_callback_support_status( &
240 compute_arguments_handle, callback_name, support_status, ierr)
242 call my_warning(
"can't get call back support_status")
247 if (support_status == kim_support_status_required)
then 248 call my_warning(
"unsupported required call back")
255 model_is_compatible = .true.
267 model_species, num_species, ier)
268 use,
intrinsic :: iso_c_binding
272 type(kim_model_handle_type),
intent(in) :: model_handle
273 integer(c_int),
intent(in) :: max_species
274 type(kim_species_name_type),
intent(out) :: model_species(max_species)
275 integer(c_int),
intent(out) :: num_species
276 integer(c_int),
intent(out) :: ier
280 integer(c_int) total_num_species
281 type(kim_species_name_type) :: species_name
282 integer(c_int) species_is_supported
288 call kim_get_number_of_species_names(total_num_species)
290 if (total_num_species .gt. max_species)
return 293 do i=1,total_num_species
294 call kim_get_species_name(i, species_name, ier)
295 call kim_get_species_support_and_code(model_handle, species_name, &
296 species_is_supported, code, ier)
297 if ((ier == 0) .and. (species_is_supported .ne. 0))
then 298 num_species = num_species + 1
299 model_species(num_species) = species_name
309 do_update_list,coordsave, &
311 use,
intrinsic :: iso_c_binding
314 integer(c_int),
parameter :: cd = c_double
317 integer(c_int),
intent(in) :: DIM
318 integer(c_int),
intent(in) :: N
319 real(c_double),
intent(in) :: coords(dim,n)
320 real(c_double),
intent(in) :: cutoff
321 real(c_double),
intent(in) :: cutpad
322 logical,
intent(inout) :: do_update_list
323 real(c_double),
intent(inout) :: coordsave(dim,n)
325 integer(c_int),
intent(out) :: ierr
328 real(c_double) disp, disp1, disp2, cutrange, dispvec(dim)
337 if (.not.do_update_list)
then 344 dispvec(1:dim) = coords(1:dim,i) - coordsave(1:dim,i)
345 disp = sqrt( dot_product(dispvec,dispvec) )
346 if (disp >= disp1)
then 349 else if (disp >= disp2)
then 353 do_update_list = ( disp1 + disp2 > cutpad )
357 if (do_update_list)
then 360 coordsave(1:dim,1:n) = coords(1:dim,1:n)
363 cutrange = cutoff + cutpad
368 do_update_list = .false.
383 coords, cutoff, neighObject)
384 use,
intrinsic :: iso_c_binding
389 logical,
intent(in) :: half
390 integer(c_int),
intent(in) :: numberOfParticles
391 real(c_double),
dimension(3,numberOfParticles), &
393 real(c_double),
intent(in) :: cutoff
397 integer(c_int),
pointer :: neighborList(:,:)
398 integer(c_int) i, j, a
401 real(c_double) cutoff2
403 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
404 [numberofparticles+1, numberofparticles])
406 neighobject%cutoff = cutoff
410 do i=1,numberofparticles
412 do j=1,numberofparticles
413 dx(:) = coords(:, j) - coords(:, i)
414 r2 = dot_product(dx, dx)
415 if (r2.le.cutoff2)
then 417 if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) )
then 419 neighborlist(a,i) = j
424 neighborlist(1,i) = a-1
449 periodic, coords, MiddlePartId)
450 use,
intrinsic :: iso_c_binding
452 integer(c_int),
parameter :: cd = c_double
455 real(c_double),
intent(in) :: FCCspacing
456 integer(c_int),
intent(in) :: nCellsPerSide
457 logical,
intent(in) :: periodic
458 real(c_double),
intent(out) :: coords(3,*)
459 integer(c_int),
intent(out) :: MiddlePartId
463 real(c_double) FCCshifts(3,4)
464 real(c_double) latVec(3)
465 integer(c_int) a, i, j, k, m
469 fccshifts(1,1) = 0.0_cd
470 fccshifts(2,1) = 0.0_cd
471 fccshifts(3,1) = 0.0_cd
472 fccshifts(1,2) = 0.5_cd*fccspacing
473 fccshifts(2,2) = 0.5_cd*fccspacing
474 fccshifts(3,2) = 0.0_cd
475 fccshifts(1,3) = 0.5_cd*fccspacing
476 fccshifts(2,3) = 0.0_cd
477 fccshifts(3,3) = 0.5_cd*fccspacing
478 fccshifts(1,4) = 0.0_cd
479 fccshifts(2,4) = 0.5_cd*fccspacing
480 fccshifts(3,4) = 0.5_cd*fccspacing
485 latvec(1) = (i-1)*fccspacing
487 latvec(2) = (j-1)*fccspacing
489 latvec(3) = (k-1)*fccspacing
492 coords(:,a) = latvec + fccshifts(:,m)
493 if ((i.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
494 (k.eq.ncellsperside/2+1) .and. (m.eq.1))
then 495 coords(:,1) = latvec + fccshifts(:,m)
500 if (.not. periodic)
then 503 latvec(1) = ncellsperside*fccspacing
504 latvec(2) = (i-1)*fccspacing
505 latvec(3) = (j-1)*fccspacing
506 a = a+1; coords(:,a) = latvec
507 a = a+1; coords(:,a) = latvec + fccshifts(:,4)
509 latvec(1) = (i-1)*fccspacing
510 latvec(2) = ncellsperside*fccspacing
511 latvec(3) = (j-1)*fccspacing
512 a = a+1; coords(:,a) = latvec
513 a = a+1; coords(:,a) = latvec + fccshifts(:,3)
515 latvec(1) = (i-1)*fccspacing
516 latvec(2) = (j-1)*fccspacing
517 latvec(3) = ncellsperside*fccspacing
518 a = a+1; coords(:,a) = latvec
519 a = a+1; coords(:,a) = latvec + fccshifts(:,2)
522 if (.not. periodic)
then 524 latvec(1) = (i-1)*fccspacing
525 latvec(2) = ncellsperside*fccspacing
526 latvec(3) = ncellsperside*fccspacing
527 a = a+1; coords(:,a) = latvec
528 latvec(1) = ncellsperside*fccspacing
529 latvec(2) = (i-1)*fccspacing
530 latvec(3) = ncellsperside*fccspacing
531 a = a+1; coords(:,a) = latvec
532 latvec(1) = ncellsperside*fccspacing
533 latvec(2) = ncellsperside*fccspacing
534 latvec(3) = (i-1)*fccspacing
535 a = a+1; coords(:,a) = latvec
538 if (.not. periodic)
then 540 a = a+1; coords(:,a) = ncellsperside*fccspacing
548 compute_arguments_handle,DIM,N,coords,cutoff,cutpad,energy,do_update_list, &
549 coordsave,neighObject,deriv,deriv_err,ierr)
550 use,
intrinsic :: iso_c_binding
554 integer(c_int),
parameter :: cd = c_double
557 integer(c_int),
intent(in) :: partnum
558 integer(c_int),
intent(in) :: dir
559 type(kim_model_handle_type),
intent(in) :: model_handle
560 type(kim_compute_arguments_handle_type),
intent(in) :: compute_arguments_handle
561 integer(c_int),
intent(in) :: DIM
562 integer(c_int),
intent(in) :: N
563 real(c_double),
intent(inout) :: coords(dim,n)
564 real(c_double),
intent(in) :: cutoff
565 real(c_double),
intent(in) :: cutpad
566 real(c_double),
intent(inout) :: energy
567 logical,
intent(inout) :: do_update_list
568 real(c_double),
intent(inout) :: coordsave(dim,n)
570 real(c_double),
intent(out) :: deriv
571 real(c_double),
intent(out) :: deriv_err
572 integer(c_int),
intent(out) :: ierr
575 real(c_double),
parameter :: eps_init = 1.e-6_cd
576 integer(c_int),
parameter :: number_eps_levels = 15
577 real(c_double) eps, deriv_last, deriv_err_last
588 deriv_err_last = huge(1.0_cd)
589 do i=1,number_eps_levels
590 deriv =
dfridr(eps,deriv_err)
592 call my_error(
"compute_numer_deriv")
594 if (deriv_err>deriv_err_last)
then 596 deriv_err = deriv_err_last
601 deriv_err_last = deriv_err
625 real(c_double) recursive function dfridr(h,err)
629 real(c_double),
intent(inout) :: h
630 real(c_double),
intent(out) :: err
633 integer(c_int),
parameter :: NTAB=10
634 real(c_double),
parameter :: CON=1.4_cd
635 real(c_double),
parameter :: CON2=con*con
636 real(c_double),
parameter :: BIG=huge(1.0_cd)
637 real(c_double),
parameter :: SAFE=2.0_cd
640 real(c_double) errt,fac,hh,a(ntab,ntab),fp,fm,coordorig
644 if (abs(h).le.tiny(0.0_cd))
then 650 coordorig = coords(dir,partnum)
651 coords(dir,partnum) = coordorig + hh
653 do_update_list,coordsave, &
655 call kim_compute(model_handle, compute_arguments_handle, ierr)
657 call my_error(
"kim_api_model_compute")
660 coords(dir,partnum) = coordorig - hh
662 do_update_list,coordsave, &
664 call kim_compute(model_handle, compute_arguments_handle, ierr)
666 call my_error(
"kim_api_model_compute")
669 coords(dir,partnum) = coordorig
671 do_update_list,coordsave, &
673 a(1,1)=(fp-fm)/(2.0_cd*hh)
680 coords(dir,partnum) = coordorig + hh
682 do_update_list,coordsave, &
684 call kim_compute(model_handle, compute_arguments_handle, ierr)
686 call my_error(
"kim_api_model_compute")
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
700 do_update_list,coordsave, &
702 a(1,i)=(fp-fm)/(2.0_cd*hh)
707 a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.0_cd)
711 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
712 if (errt.le.err)
then 717 if (abs(a(i,i)-a(i-1,i-1)).ge.safe*err)
return 742 use,
intrinsic :: iso_c_binding
747 integer(c_int),
parameter :: cd = c_double
749 integer(c_int),
parameter :: nCellsPerSide = 2
750 integer(c_int),
parameter :: DIM = 3
751 real(c_double),
parameter :: cutpad = 0.75_cd
752 integer(c_int),
parameter :: max_species = 200
753 real(c_double),
parameter :: eps_prec = epsilon(1.0_cd)
754 real(c_double) FCCspacing
756 integer(c_int),
parameter :: &
757 N = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
758 real(c_double),
allocatable :: forces_num(:,:)
759 real(c_double),
allocatable :: forces_num_err(:,:)
760 type(kim_species_name_type) :: model_species(max_species)
761 integer(c_int),
target :: num_species
762 character(len=5, kind=c_char) :: passfail
763 real(c_double) :: forcediff
764 real(c_double) :: forcediff_sumsq
765 real(c_double) :: weight
766 real(c_double) :: weight_sum
767 real(c_double) :: alpha
768 real(c_double) :: term
769 real(c_double) :: term_max
770 real(c_double),
allocatable :: cluster_coords(:,:)
771 real(c_double),
allocatable :: cluster_disps(:,:)
772 type(kim_species_name_type),
allocatable :: cluster_species(:)
773 integer(c_int) I,J,Imax,Jmax,species_code
774 integer(c_int) seed_size
775 integer(c_int),
allocatable :: seed(:)
781 integer(c_int),
allocatable,
target :: neighborList(:,:)
782 real(c_double),
allocatable :: coordsave(:,:)
783 logical do_update_list
788 character(len=256, kind=c_char) :: testname =
"vc_forces_numer_deriv" 789 character(len=256, kind=c_char) :: modelname
791 type(kim_model_handle_type) :: model_handle
792 type(kim_compute_arguments_handle_type) :: compute_arguments_handle
793 integer(c_int) ierr, ierr2
794 integer(c_int) species_is_supported
795 integer(c_int),
target :: numberOfParticles
796 integer(c_int),
target :: particleSpeciesCodes(n)
797 integer(c_int),
target :: particleContributing(n)
798 real(c_double) :: influence_distance
799 integer(c_int) :: number_of_neighbor_lists
800 real(c_double),
allocatable :: cutoffs(:)
801 integer(c_int),
allocatable :: &
802 model_will_not_request_neighbors_of_noncontributing_particles(:)
803 real(c_double) :: cutoff
804 real(c_double),
target :: energy
805 real(c_double),
target :: coords(3,n)
806 real(c_double),
target :: forces_kim(3,n)
807 real(c_double) :: forces(3,n)
808 integer(c_int) middleDum
809 logical forces_optional
810 logical model_is_compatible
811 integer(c_int) number_of_model_routine_names
812 type(kim_model_routine_name_type) model_routine_name
813 integer(c_int) present
814 integer(c_int) required
815 integer(c_int) requested_units_accepted
816 real(c_double) rnd, deriv, deriv_err
817 real(c_double),
pointer :: null_pointer
819 nullify(null_pointer)
821 numberofparticles = n
831 call random_seed(size=seed_size)
832 allocate(seed(seed_size))
834 call random_seed(put=seed)
838 print
'("Please enter a valid KIM model name: ")' 844 print *,
'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES' 847 print
'("This is Test : ",A)', trim(testname)
848 print
'("Results for KIM Model : ",A)', trim(modelname)
852 call kim_model_create(kim_numbering_one_based, &
854 kim_energy_unit_ev, &
856 kim_temperature_unit_k, &
859 requested_units_accepted, &
866 if (requested_units_accepted == 0)
then 867 call my_error(
"Must adapt to model units")
871 call kim_get_number_of_model_routine_names(number_of_model_routine_names)
872 do i=1,number_of_model_routine_names
873 call kim_get_model_routine_name(i, model_routine_name, ierr)
875 call my_error(
"kim_get_model_routine_name")
877 call kim_is_routine_present(model_handle, model_routine_name,
present, &
880 call my_error(
"kim_is_routine_present")
883 if ((
present == 1) .and. (required == 1))
then 884 if (.not. ((model_routine_name == kim_model_routine_name_create) &
885 .or. (model_routine_name == &
886 kim_model_routine_name_compute_arguments_create) &
887 .or. (model_routine_name == kim_model_routine_name_compute) &
888 .or. (model_routine_name == kim_model_routine_name_refresh) &
889 .or. (model_routine_name == &
890 kim_model_routine_name_compute_arguments_destroy) &
891 .or. (model_routine_name == kim_model_routine_name_destroy)))
then 893 call my_error(
"Unknown required ModelRoutineName found.")
899 call kim_compute_arguments_create(model_handle, &
900 compute_arguments_handle, ierr)
902 call my_error(
"kim_model_compute_arguments_create")
906 model_is_compatible, ierr)
908 call my_error(
"error checking compatibility")
910 if (.not. model_is_compatible)
then 911 call my_error(
"incompatibility reported by check_model_compatibility")
919 call my_error(
"Get_Model_Supported_Species")
923 allocate(cluster_coords(3,n),cluster_disps(3,n),cluster_species(n))
925 call random_number(rnd)
926 species_code = 1 + int(rnd*num_species)
927 cluster_species(i) = model_species(species_code)
933 cluster_coords, middledum)
938 call random_number(rnd)
939 cluster_disps(j,i) = 0.1_cd*(rnd-0.5_cd)
945 call kim_set_argument_pointer(compute_arguments_handle, &
946 kim_compute_argument_name_number_of_particles, numberofparticles, ierr2)
948 call kim_set_argument_pointer(compute_arguments_handle, &
949 kim_compute_argument_name_particle_species_codes, particlespeciescodes, &
952 call kim_set_argument_pointer(compute_arguments_handle, &
953 kim_compute_argument_name_particle_contributing, particlecontributing, &
956 call kim_set_argument_pointer(compute_arguments_handle, &
957 kim_compute_argument_name_coordinates, coords, ierr2)
959 call kim_set_argument_pointer(compute_arguments_handle, &
960 kim_compute_argument_name_partial_energy, energy, ierr2)
962 call kim_set_argument_pointer(compute_arguments_handle, &
963 kim_compute_argument_name_partial_forces, forces_kim, ierr2)
966 call my_error(
"set_argument_pointer")
972 allocate(neighborlist(n+1,n))
973 neighobject%neighbor_list_pointer = c_loc(neighborlist)
974 neighobject%number_of_particles = n
978 call kim_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)
982 call my_error(
"set_callback_pointer")
985 call kim_get_influence_distance(model_handle, influence_distance)
986 call kim_get_number_of_neighbor_lists(model_handle, &
987 number_of_neighbor_lists)
988 allocate(cutoffs(number_of_neighbor_lists), &
989 model_will_not_request_neighbors_of_noncontributing_particles( &
990 number_of_neighbor_lists))
991 call kim_get_neighbor_list_values(model_handle, cutoffs, &
992 model_will_not_request_neighbors_of_noncontributing_particles, ierr)
994 call my_error(
"get_neighbor_list_values")
996 cutoff = maxval(cutoffs)
999 fccspacing = 0.75_cd*cutoff
1002 cluster_coords(:,i) = fccspacing*cluster_coords(:,i)
1004 print
'("Using FCC lattice parameter: ",f12.5)', fccspacing
1007 call kim_get_species_support_and_code(model_handle, &
1008 cluster_species(i), species_is_supported, particlespeciescodes(i), ierr)
1011 call my_error(
"kim_api_get_species_code")
1014 particlecontributing(i) = 1
1017 coords(:,i) = cluster_coords(:,i) + cluster_disps(:,i)
1022 do_update_list = .true.
1023 allocate(coordsave(dim,n))
1025 do_update_list,coordsave, &
1028 call my_error(
"update_neighborlist")
1033 call kim_compute(model_handle, compute_arguments_handle, ierr)
1035 call my_error(
"kim_api_model_compute")
1045 print
'("Energy = ",ES25.15)', energy
1051 if (forces_optional)
then 1052 call kim_set_argument_pointer(compute_arguments_handle, &
1053 kim_compute_argument_name_partial_forces, null_pointer, ierr)
1055 call my_error(
"set_argument_pointer")
1061 allocate(forces_num(dim,n),forces_num_err(dim,n))
1065 dim,n,coords,cutoff, &
1066 cutpad, energy, do_update_list, &
1067 coordsave,neighobject,deriv,deriv_err,ierr)
1069 call my_error(
"compute_numer_deriv")
1071 forces_num(j,i) = -deriv
1072 forces_num_err(j,i) = deriv_err
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", &
1081 forcediff_sumsq = 0.0_cd
1085 forcediff = abs(forces(j,i)-forces_num(j,i))
1086 if (forcediff<forces_num_err(j,i))
then 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 1099 forcediff_sumsq = forcediff_sumsq + term
1100 weight_sum = weight_sum + weight
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
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
1113 alpha = sqrt(forcediff_sumsq/weight_sum)/dble(dim*n)
1115 print
'("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," (units of force)")', &
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))
1125 deallocate(forces_num)
1126 deallocate(forces_num_err)
1127 deallocate(neighborlist)
1128 deallocate(coordsave)
1129 deallocate(cutoffs, model_will_not_request_neighbors_of_noncontributing_particles)
1131 call kim_compute_arguments_destroy(model_handle, &
1132 compute_arguments_handle, ierr)
1134 call my_error(
"kim_model_compute_arguments_destroy")
1136 call kim_model_destroy(model_handle)
1141 print
'(120(''-''))' 1145 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)