32 use,
intrinsic :: iso_c_binding
37 recursive subroutine my_error(message)
39 character(len=*, kind=c_char),
intent(in) :: message
41 print *,
"* Error : ", trim(message)
47 character(len=*, kind=c_char),
intent(in) :: message
49 print *,
"* Warning : ", trim(message)
63 use,
intrinsic :: iso_c_binding
68 real(c_double) :: cutoff
69 integer(c_int) :: number_of_particles
70 type(c_ptr) :: neighbor_list_pointer
88 recursive subroutine get_neigh(data_object, number_of_neighbor_lists, &
89 cutoffs, neighbor_list_index, request, &
90 numnei, pnei1part, ierr) bind(c)
95 type(c_ptr),
value,
intent(in) :: data_object
96 integer(c_int),
value,
intent(in) :: number_of_neighbor_lists
97 real(c_double),
intent(in) :: cutoffs(*)
98 integer(c_int),
value,
intent(in) :: neighbor_list_index
99 integer(c_int),
value,
intent(in) :: request
100 integer(c_int),
intent(out) :: numnei
101 type(c_ptr),
intent(out) :: pnei1part
102 integer(c_int),
intent(out) :: ierr
105 integer(c_int) numberofparticles
107 integer(c_int),
pointer :: neighborlist(:, :)
109 if (number_of_neighbor_lists > 1)
then 110 call my_warning(
"Model requires too many neighbor lists")
115 call c_f_pointer(data_object, neighobject)
116 numberofparticles = neighobject%number_of_particles
117 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
118 [numberofparticles + 1, numberofparticles])
120 if (cutoffs(neighbor_list_index) > neighobject%cutoff)
then 121 call my_warning(
"neighbor list cutoff too small for model cutoff")
126 if ((request > numberofparticles) .or. (request < 1))
then 128 call my_warning(
"Invalid part ID in get_neigh")
134 numnei = neighborlist(1, request)
137 pnei1part = c_loc(neighborlist(2, request))
158 compute_arguments_handle, forces_optional, model_is_compatible, ierr)
159 use,
intrinsic :: iso_c_binding
164 type(kim_compute_arguments_handle_type),
intent(in) :: &
165 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
179 model_is_compatible = .false.
180 forces_optional = .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_handle, max_species, 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
296 call kim_get_number_of_species_names(total_num_species)
298 if (total_num_species > max_species)
return 301 do i = 1, total_num_species
302 call kim_get_species_name(i, species_name, ier)
303 call kim_get_species_support_and_code(model_handle, species_name, &
304 species_is_supported, code, ier)
305 if ((ier == 0) .and. (species_is_supported /= 0))
then 306 num_species = num_species + 1
307 model_species(num_species) = species_name
317 do_update_list, coordsave, &
319 use,
intrinsic :: iso_c_binding
322 integer(c_int),
parameter :: cd = c_double
325 integer(c_int),
intent(in) :: DIM
326 integer(c_int),
intent(in) :: N
327 real(c_double),
intent(in) :: coords(dim, n)
328 real(c_double),
intent(in) :: cutoff
329 real(c_double),
intent(in) :: cutpad
330 logical,
intent(inout) :: do_update_list
331 real(c_double),
intent(inout) :: coordsave(dim, n)
333 integer(c_int),
intent(out) :: ierr
336 real(c_double) disp, disp1, disp2, cutrange, dispvec(dim)
345 if (.not. do_update_list)
then 352 dispvec(1:dim) = coords(1:dim, i) - coordsave(1:dim, i)
353 disp = sqrt(dot_product(dispvec, dispvec))
354 if (disp >= disp1)
then 357 else if (disp >= disp2)
then 361 do_update_list = (disp1 + disp2 > cutpad)
365 if (do_update_list)
then 368 coordsave(1:dim, 1:n) = coords(1:dim, 1:n)
371 cutrange = cutoff + cutpad
376 do_update_list = .false.
389 half, numberOfParticles, 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 <= cutoff2)
then 423 if ((j > i) .OR. ((.not. half) .AND. (i /= 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
490 do i = 1, ncellsperside
491 latvec(1) = (i - 1) * fccspacing
492 do j = 1, ncellsperside
493 latvec(2) = (j - 1) * fccspacing
494 do k = 1, ncellsperside
495 latvec(3) = (k - 1) * fccspacing
498 coords(:, a) = latvec + fccshifts(:, m)
499 if ((i == ncellsperside / 2 + 1) &
500 .and. (j == ncellsperside / 2 + 1) &
501 .and. (k == ncellsperside / 2 + 1) &
505 coords(:, 1) = latvec + fccshifts(:, m)
510 if (.not. periodic)
then 513 latvec(1) = ncellsperside * fccspacing
514 latvec(2) = (i - 1) * fccspacing
515 latvec(3) = (j - 1) * fccspacing
516 a = a + 1; coords(:, a) = latvec
517 a = a + 1; coords(:, a) = latvec + fccshifts(:, 4)
519 latvec(1) = (i - 1) * fccspacing
520 latvec(2) = ncellsperside * fccspacing
521 latvec(3) = (j - 1) * fccspacing
522 a = a + 1; coords(:, a) = latvec
523 a = a + 1; coords(:, a) = latvec + fccshifts(:, 3)
525 latvec(1) = (i - 1) * fccspacing
526 latvec(2) = (j - 1) * fccspacing
527 latvec(3) = ncellsperside * fccspacing
528 a = a + 1; coords(:, a) = latvec
529 a = a + 1; coords(:, a) = latvec + fccshifts(:, 2)
532 if (.not. periodic)
then 534 latvec(1) = (i - 1) * fccspacing
535 latvec(2) = ncellsperside * fccspacing
536 latvec(3) = ncellsperside * fccspacing
537 a = a + 1; coords(:, a) = latvec
538 latvec(1) = ncellsperside * fccspacing
539 latvec(2) = (i - 1) * fccspacing
540 latvec(3) = ncellsperside * fccspacing
541 a = a + 1; coords(:, a) = latvec
542 latvec(1) = ncellsperside * fccspacing
543 latvec(2) = ncellsperside * fccspacing
544 latvec(3) = (i - 1) * fccspacing
545 a = a + 1; coords(:, a) = latvec
548 if (.not. periodic)
then 550 a = a + 1; coords(:, a) = ncellsperside * fccspacing
558 partnum, dir, model_handle, compute_arguments_handle, DIM, N, coords, &
559 cutoff, cutpad, energy, do_update_list, coordsave, neighObject, deriv, &
561 use,
intrinsic :: iso_c_binding
565 integer(c_int),
parameter :: cd = c_double
568 integer(c_int),
intent(in) :: partnum
569 integer(c_int),
intent(in) :: dir
570 type(kim_model_handle_type),
intent(in) :: model_handle
571 type(kim_compute_arguments_handle_type),
intent(in) :: &
572 compute_arguments_handle
573 integer(c_int),
intent(in) :: DIM
574 integer(c_int),
intent(in) :: N
575 real(c_double),
intent(inout) :: coords(dim, n)
576 real(c_double),
intent(in) :: cutoff
577 real(c_double),
intent(in) :: cutpad
578 real(c_double),
intent(inout) :: energy
579 logical,
intent(inout) :: do_update_list
580 real(c_double),
intent(inout) :: coordsave(dim, n)
582 real(c_double),
intent(out) :: deriv
583 real(c_double),
intent(out) :: deriv_err
584 integer(c_int),
intent(out) :: ierr
587 real(c_double),
parameter :: eps_init = 1.e-6_cd
588 integer(c_int),
parameter :: number_eps_levels = 15
589 real(c_double) eps, deriv_last, deriv_err_last
600 deriv_err_last = huge(1.0_cd)
601 do i = 1, number_eps_levels
602 deriv =
dfridr(eps, deriv_err)
604 call my_error(
"compute_numer_deriv")
606 if (deriv_err > deriv_err_last)
then 608 deriv_err = deriv_err_last
613 deriv_err_last = deriv_err
637 real(c_double) recursive function dfridr(h, err)
641 real(c_double),
intent(inout) :: h
642 real(c_double),
intent(out) :: err
645 integer(c_int),
parameter :: NTAB = 10
646 real(c_double),
parameter :: CON = 1.4_cd
647 real(c_double),
parameter :: CON2 = con * con
648 real(c_double),
parameter :: BIG = huge(1.0_cd)
649 real(c_double),
parameter :: SAFE = 2.0_cd
652 real(c_double) errt, fac, hh, a(ntab, ntab), fp, fm, coordorig
657 if (abs(h) <= tiny(0.0_cd))
then 663 coordorig = coords(dir, partnum)
664 coords(dir, partnum) = coordorig + hh
666 do_update_list, coordsave, &
668 call kim_compute(model_handle, compute_arguments_handle, ierr)
670 call my_error(
"kim_api_model_compute")
673 coords(dir, partnum) = coordorig - hh
675 do_update_list, coordsave, &
677 call kim_compute(model_handle, compute_arguments_handle, ierr)
679 call my_error(
"kim_api_model_compute")
682 coords(dir, partnum) = coordorig
684 do_update_list, coordsave, &
686 a(1, 1) = (fp - fm) / (2.0_cd * hh)
692 coords(dir, partnum) = coordorig + hh
694 do_update_list, coordsave, &
696 call kim_compute(model_handle, compute_arguments_handle, ierr)
698 call my_error(
"kim_api_model_compute")
701 coords(dir, partnum) = coordorig - hh
703 do_update_list, coordsave, &
705 call kim_compute(model_handle, compute_arguments_handle, ierr)
707 call my_error(
"kim_api_model_compute")
710 coords(dir, partnum) = coordorig
712 do_update_list, coordsave, &
714 a(1, i) = (fp - fm) / (2.0_cd * hh)
719 a(j, i) = (a(j - 1, i) * fac - a(j - 1, i - 1)) / (fac - 1.0_cd)
723 errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
724 if (errt <= err)
then 729 if (abs(a(i, i) - a(i - 1, i - 1)) >= safe * err)
return 754 use,
intrinsic :: iso_c_binding
759 integer(c_int),
parameter :: cd = c_double
761 integer(c_int),
parameter :: nCellsPerSide = 2
762 integer(c_int),
parameter :: DIM = 3
763 real(c_double),
parameter :: cutpad = 0.75_cd
764 integer(c_int),
parameter :: max_species = 200
765 real(c_double),
parameter :: eps_prec = epsilon(1.0_cd)
766 real(c_double) FCCspacing
768 integer(c_int),
parameter :: N = 4 * (ncellsperside)**3 &
769 + 6 * (ncellsperside)**2 &
770 + 3 * (ncellsperside) + 1
771 real(c_double),
allocatable :: forces_num(:, :)
772 real(c_double),
allocatable :: forces_num_err(:, :)
773 type(kim_species_name_type) :: model_species(max_species)
774 integer(c_int),
target :: num_species
775 character(len=5, kind=c_char) :: passfail
776 real(c_double) :: forcediff
777 real(c_double) :: forcediff_sumsq
778 real(c_double) :: weight
779 real(c_double) :: weight_sum
780 real(c_double) :: alpha
781 real(c_double) :: term
782 real(c_double) :: term_max
783 real(c_double),
allocatable :: cluster_coords(:, :)
784 real(c_double),
allocatable :: cluster_disps(:, :)
785 type(kim_species_name_type),
allocatable :: cluster_species(:)
786 integer(c_int) I, J, Imax, Jmax, species_code
787 integer(c_int) seed_size
788 integer(c_int),
allocatable :: seed(:)
794 integer(c_int),
allocatable,
target :: neighborList(:, :)
795 real(c_double),
allocatable :: coordsave(:, :)
796 logical do_update_list
801 character(len=256, kind=c_char) :: testname =
"vc_forces_numer_deriv" 802 character(len=256, kind=c_char) :: modelname
804 type(kim_model_handle_type) :: model_handle
805 type(kim_compute_arguments_handle_type) :: compute_arguments_handle
806 integer(c_int) ierr, ierr2
807 integer(c_int) species_is_supported
808 integer(c_int),
target :: numberOfParticles
809 integer(c_int),
target :: particleSpeciesCodes(n)
810 integer(c_int),
target :: particleContributing(n)
811 real(c_double) :: influence_distance
812 integer(c_int) :: number_of_neighbor_lists
813 real(c_double),
allocatable :: cutoffs(:)
814 integer(c_int),
allocatable :: &
815 model_will_not_request_neighbors_of_noncontributing_particles(:)
816 real(c_double) :: cutoff
817 real(c_double),
target :: energy
818 real(c_double),
target :: coords(3, n)
819 real(c_double),
target :: forces_kim(3, n)
820 real(c_double) :: forces(3, n)
821 integer(c_int) middleDum
822 logical forces_optional
823 logical model_is_compatible
824 integer(c_int) number_of_model_routine_names
825 type(kim_model_routine_name_type) model_routine_name
826 integer(c_int) present
827 integer(c_int) required
828 integer(c_int) requested_units_accepted
829 real(c_double) rnd, deriv, deriv_err
830 real(c_double),
pointer :: null_pointer
832 nullify (null_pointer)
834 numberofparticles = n
844 call random_seed(size=seed_size)
845 allocate (seed(seed_size))
847 call random_seed(put=seed)
851 print
'("Please enter a valid KIM model name: ")' 852 read (*, *) modelname
857 print *,
'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES' 860 print
'("This is Test : ",A)', trim(testname)
861 print
'("Results for KIM Model : ",A)', trim(modelname)
865 call kim_model_create(kim_numbering_one_based, &
867 kim_energy_unit_ev, &
869 kim_temperature_unit_k, &
872 requested_units_accepted, &
879 if (requested_units_accepted == 0)
then 880 call my_error(
"Must adapt to model units")
884 call kim_get_number_of_model_routine_names(number_of_model_routine_names)
885 do i = 1, number_of_model_routine_names
886 call kim_get_model_routine_name(i, model_routine_name, ierr)
888 call my_error(
"kim_get_model_routine_name")
890 call kim_is_routine_present(model_handle, model_routine_name,
present, &
893 call my_error(
"kim_is_routine_present")
896 if ((
present == 1) .and. (required == 1))
then 897 if (.not. ((model_routine_name == kim_model_routine_name_create) &
898 .or. (model_routine_name == &
899 kim_model_routine_name_compute_arguments_create) &
900 .or. (model_routine_name == kim_model_routine_name_compute) &
901 .or. (model_routine_name == kim_model_routine_name_refresh) &
902 .or. (model_routine_name == &
903 kim_model_routine_name_compute_arguments_destroy) &
904 .or. (model_routine_name == kim_model_routine_name_destroy))) &
906 call my_error(
"Unknown required ModelRoutineName found.")
912 call kim_compute_arguments_create(model_handle, &
913 compute_arguments_handle, ierr)
915 call my_error(
"kim_model_compute_arguments_create")
919 model_is_compatible, ierr)
921 call my_error(
"error checking compatibility")
923 if (.not. model_is_compatible)
then 924 call my_error(
"incompatibility reported by check_model_compatibility")
932 call my_error(
"Get_Model_Supported_Species")
936 allocate (cluster_coords(3, n), cluster_disps(3, n), cluster_species(n))
938 call random_number(rnd)
939 species_code = 1 + int(rnd * num_species)
940 cluster_species(i) = model_species(species_code)
946 cluster_coords, middledum)
951 call random_number(rnd)
952 cluster_disps(j, i) = 0.1_cd * (rnd - 0.5_cd)
958 call kim_set_argument_pointer( &
959 compute_arguments_handle, &
960 kim_compute_argument_name_number_of_particles, &
961 numberofparticles, ierr2)
963 call kim_set_argument_pointer( &
964 compute_arguments_handle, &
965 kim_compute_argument_name_particle_species_codes, &
966 particlespeciescodes, ierr2)
968 call kim_set_argument_pointer( &
969 compute_arguments_handle, &
970 kim_compute_argument_name_particle_contributing, &
971 particlecontributing, ierr2)
973 call kim_set_argument_pointer( &
974 compute_arguments_handle, &
975 kim_compute_argument_name_coordinates, coords, ierr2)
977 call kim_set_argument_pointer( &
978 compute_arguments_handle, &
979 kim_compute_argument_name_partial_energy, energy, ierr2)
981 call kim_set_argument_pointer( &
982 compute_arguments_handle, &
983 kim_compute_argument_name_partial_forces, forces_kim, ierr2)
986 call my_error(
"set_argument_pointer")
992 allocate (neighborlist(n + 1, n))
993 neighobject%neighbor_list_pointer = c_loc(neighborlist)
994 neighobject%number_of_particles = n
998 call kim_set_callback_pointer( &
999 compute_arguments_handle, &
1000 kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
1001 c_funloc(
get_neigh), c_loc(neighobject), ierr)
1003 call my_error(
"set_callback_pointer")
1006 call kim_get_influence_distance(model_handle, influence_distance)
1007 call kim_get_number_of_neighbor_lists(model_handle, &
1008 number_of_neighbor_lists)
1009 allocate (cutoffs(number_of_neighbor_lists), &
1010 model_will_not_request_neighbors_of_noncontributing_particles( &
1011 number_of_neighbor_lists))
1012 call kim_get_neighbor_list_values( &
1013 model_handle, cutoffs, &
1014 model_will_not_request_neighbors_of_noncontributing_particles, ierr)
1016 call my_error(
"get_neighbor_list_values")
1018 cutoff = maxval(cutoffs)
1021 fccspacing = 0.75_cd * cutoff
1024 cluster_coords(:, i) = fccspacing * cluster_coords(:, i)
1026 print
'("Using FCC lattice parameter: ",f12.5)', fccspacing
1029 call kim_get_species_support_and_code( &
1030 model_handle, cluster_species(i), species_is_supported, &
1031 particlespeciescodes(i), ierr)
1034 call my_error(
"kim_api_get_species_code")
1037 particlecontributing(i) = 1
1040 coords(:, i) = cluster_coords(:, i) + cluster_disps(:, i)
1045 do_update_list = .true.
1046 allocate (coordsave(dim, n))
1048 do_update_list, coordsave, &
1051 call my_error(
"update_neighborlist")
1056 call kim_compute(model_handle, compute_arguments_handle, ierr)
1058 call my_error(
"kim_api_model_compute")
1068 print
'("Energy = ",ES25.15)', energy
1074 if (forces_optional)
then 1075 call kim_set_argument_pointer( &
1076 compute_arguments_handle, kim_compute_argument_name_partial_forces, &
1079 call my_error(
"set_argument_pointer")
1085 allocate (forces_num(dim, n), forces_num_err(dim, n))
1089 dim, n, coords, cutoff, &
1090 cutpad, energy, do_update_list, &
1091 coordsave, neighobject, deriv, deriv_err, ierr)
1093 call my_error(
"compute_numer_deriv")
1095 forces_num(j, i) = -deriv
1096 forces_num_err(j, i) = deriv_err
1102 print
'(A6,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)',
"Part",
"Spec",
"Dir", &
1103 "Force_model",
"Force_numer",
"Force diff",
"pred error",
"weight", &
1105 forcediff_sumsq = 0.0_cd
1109 forcediff = abs(forces(j, i) - forces_num(j, i))
1110 if (forcediff < forces_num_err(j, i))
then 1115 weight = max(abs(forces_num(j, i)), eps_prec) / &
1116 max(abs(forces_num_err(j, i)), eps_prec)
1117 term = weight * forcediff**2
1118 if (term > term_max)
then 1123 forcediff_sumsq = forcediff_sumsq + term
1124 weight_sum = weight_sum + weight
1126 print
'(I6,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1127 i, particlespeciescodes(i), j, forces(j, i), forces_num(j, i), &
1128 forcediff, forces_num_err(j, i), weight, passfail
1130 print
'(14X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1131 j, forces(j, i), forces_num(j, i), &
1132 forcediff, forces_num_err(j, i), weight, passfail
1137 alpha = sqrt(forcediff_sumsq / weight_sum) / dble(dim * n)
1139 print
'("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," & 1140 &(units of force)")', alpha
1142 print
'(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,'// &
1143 ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
1144 imax, jmax, abs(forces(jmax, imax) - forces_num(jmax, imax)), &
1145 abs(forces(jmax, imax) - forces_num(jmax, imax)) / abs(forces(jmax, imax))
1149 deallocate (forces_num)
1150 deallocate (forces_num_err)
1151 deallocate (neighborlist)
1152 deallocate (coordsave)
1153 deallocate (cutoffs)
1154 deallocate (model_will_not_request_neighbors_of_noncontributing_particles)
1156 call kim_compute_arguments_destroy(model_handle, &
1157 compute_arguments_handle, ierr)
1159 call my_error(
"kim_model_compute_arguments_destroy")
1161 call kim_model_destroy(model_handle)
1166 print
'(120(''-''))' 1170 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)