29 use,
intrinsic :: iso_c_binding
34 recursive subroutine my_error(message)
36 character(len=*, kind=c_char),
intent(in) :: message
38 print *,
"* Error : ", trim(message)
44 character(len=*, kind=c_char),
intent(in) :: message
46 print *,
"* Warning : ", trim(message)
60 use,
intrinsic :: iso_c_binding
65 real(c_double) :: cutoff
66 integer(c_int) :: number_of_particles
67 type(c_ptr) :: neighbor_list_pointer
85 recursive subroutine get_neigh(data_object, number_of_neighbor_lists, &
86 cutoffs, neighbor_list_index, request, &
87 numnei, pnei1part, ierr) bind(c)
92 type(c_ptr),
value,
intent(in) :: data_object
93 integer(c_int),
value,
intent(in) :: number_of_neighbor_lists
94 real(c_double),
intent(in) :: cutoffs(*)
95 integer(c_int),
value,
intent(in) :: neighbor_list_index
96 integer(c_int),
value,
intent(in) :: request
97 integer(c_int),
intent(out) :: numnei
98 type(c_ptr),
intent(out) :: pnei1part
99 integer(c_int),
intent(out) :: ierr
102 integer(c_int) numberofparticles
104 integer(c_int),
pointer :: neighborlist(:, :)
106 if (number_of_neighbor_lists > 1)
then 107 call my_warning(
"Model requires too many neighbor lists")
112 call c_f_pointer(data_object, neighobject)
113 numberofparticles = neighobject%number_of_particles
114 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
115 [numberofparticles + 1, numberofparticles])
117 if (cutoffs(neighbor_list_index) > neighobject%cutoff)
then 118 call my_warning(
"neighbor list cutoff too small for model cutoff")
123 if ((request > numberofparticles) .or. (request < 1))
then 125 call my_warning(
"Invalid part ID in get_neigh")
131 numnei = neighborlist(1, request)
134 pnei1part = c_loc(neighborlist(2, request))
155 compute_arguments_handle, forces_optional, particle_energy_supported, &
156 particle_energy_optional, model_is_compatible, ierr)
157 use,
intrinsic :: iso_c_binding
162 type(kim_compute_arguments_handle_type),
intent(in) :: &
163 compute_arguments_handle
164 logical,
intent(out) :: forces_optional
165 logical,
intent(out) :: particle_energy_supported
166 logical,
intent(out) :: particle_energy_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 particle_energy_supported = .false.
181 particle_energy_optional = .false.
182 forces_optional = .false.
186 call kim_get_number_of_compute_argument_names( &
187 number_of_argument_names)
188 do i = 1, number_of_argument_names
189 call kim_get_compute_argument_name(i, argument_name, &
195 call kim_get_argument_support_status( &
196 compute_arguments_handle, argument_name, support_status, ierr)
198 call my_warning(
"can't get argument support_status")
203 if (support_status == kim_support_status_required)
then 205 (argument_name == kim_compute_argument_name_partial_energy) .or. &
207 kim_compute_argument_name_partial_particle_energy) .or. &
208 (argument_name == kim_compute_argument_name_partial_forces)))
then 209 call my_warning(
"unsupported required argument")
216 if ((argument_name == kim_compute_argument_name_partial_energy) .and. &
217 (support_status == kim_support_status_not_supported))
then 218 call my_warning(
"model does not support energy")
222 if (argument_name == kim_compute_argument_name_partial_forces)
then 223 if (support_status == kim_support_status_not_supported)
then 224 call my_warning(
"model does not support forces")
227 else if (support_status == kim_support_status_required)
then 228 forces_optional = .false.
229 else if (support_status == kim_support_status_optional)
then 230 forces_optional = .true.
232 call my_warning(
"unknown support_status for forces")
239 if (argument_name == &
240 kim_compute_argument_name_partial_particle_energy)
then 241 if (support_status == kim_support_status_not_supported)
then 242 call my_warning(
"model does not support partial_particle_energy. & 243 &The associated checks will be disabled.")
244 particle_energy_supported = .false.
245 particle_energy_optional = .false.
246 else if (support_status == kim_support_status_required)
then 247 particle_energy_supported = .true.
248 particle_energy_optional = .false.
249 else if (support_status == kim_support_status_optional)
then 250 particle_energy_supported = .true.
251 particle_energy_optional = .true.
253 call my_warning(
"unknown support_status for particle energy")
261 call kim_get_number_of_compute_callback_names( &
262 number_of_callback_names)
263 do i = 1, number_of_callback_names
264 call kim_get_compute_callback_name(i, callback_name, &
270 call kim_get_callback_support_status( &
271 compute_arguments_handle, callback_name, support_status, ierr)
273 call my_warning(
"can't get call back support_status")
278 if (support_status == kim_support_status_required)
then 279 call my_warning(
"unsupported required call back")
286 model_is_compatible = .true.
298 model_handle, max_species, model_species, num_species, ier)
299 use,
intrinsic :: iso_c_binding
303 type(kim_model_handle_type),
intent(in) :: model_handle
304 integer(c_int),
intent(in) :: max_species
305 type(kim_species_name_type),
intent(out) :: model_species(max_species)
306 integer(c_int),
intent(out) :: num_species
307 integer(c_int),
intent(out) :: ier
311 integer(c_int) total_num_species
312 type(kim_species_name_type) :: species_name
313 integer(c_int) species_is_supported
321 call kim_get_number_of_species_names(total_num_species)
323 if (total_num_species > max_species)
return 326 do i = 1, total_num_species
327 call kim_get_species_name(i, species_name, ier)
328 call kim_get_species_support_and_code(model_handle, species_name, &
329 species_is_supported, code, ier)
330 if ((ier == 0) .and. (species_is_supported /= 0))
then 331 num_species = num_species + 1
332 model_species(num_species) = species_name
342 do_update_list, coordsave, &
344 use,
intrinsic :: iso_c_binding
347 integer(c_int),
parameter :: cd = c_double
350 integer(c_int),
intent(in) :: DIM
351 integer(c_int),
intent(in) :: N
352 real(c_double),
intent(in) :: coords(dim, n)
353 real(c_double),
intent(in) :: cutoff
354 real(c_double),
intent(in) :: cutpad
355 logical,
intent(inout) :: do_update_list
356 real(c_double),
intent(inout) :: coordsave(dim, n)
358 integer(c_int),
intent(out) :: ierr
361 real(c_double) disp, disp1, disp2, cutrange, dispvec(dim)
370 if (.not. do_update_list)
then 377 dispvec(1:dim) = coords(1:dim, i) - coordsave(1:dim, i)
378 disp = sqrt(dot_product(dispvec, dispvec))
379 if (disp >= disp1)
then 382 else if (disp >= disp2)
then 386 do_update_list = (disp1 + disp2 > cutpad)
390 if (do_update_list)
then 393 coordsave(1:dim, 1:n) = coords(1:dim, 1:n)
396 cutrange = cutoff + cutpad
401 do_update_list = .false.
414 half, numberOfParticles, coords, cutoff, neighObject)
415 use,
intrinsic :: iso_c_binding
420 logical,
intent(in) :: half
421 integer(c_int),
intent(in) :: numberOfParticles
422 real(c_double),
dimension(3, numberOfParticles), &
424 real(c_double),
intent(in) :: cutoff
428 integer(c_int),
pointer :: neighborList(:, :)
429 integer(c_int) i, j, a
432 real(c_double) cutoff2
434 call c_f_pointer(neighobject%neighbor_list_pointer, neighborlist, &
435 [numberofparticles + 1, numberofparticles])
437 neighobject%cutoff = cutoff
441 do i = 1, numberofparticles
443 do j = 1, numberofparticles
444 dx(:) = coords(:, j) - coords(:, i)
445 r2 = dot_product(dx, dx)
446 if (r2 <= cutoff2)
then 448 if ((j > i) .OR. ((.not. half) .AND. (i /= j)))
then 450 neighborlist(a, i) = j
455 neighborlist(1, i) = a - 1
480 periodic, coords, MiddlePartId)
481 use,
intrinsic :: iso_c_binding
483 integer(c_int),
parameter :: cd = c_double
486 real(c_double),
intent(in) :: FCCspacing
487 integer(c_int),
intent(in) :: nCellsPerSide
488 logical,
intent(in) :: periodic
489 real(c_double),
intent(out) :: coords(3, *)
490 integer(c_int),
intent(out) :: MiddlePartId
494 real(c_double) FCCshifts(3, 4)
495 real(c_double) latVec(3)
496 integer(c_int) a, i, j, k, m
500 fccshifts(1, 1) = 0.0_cd
501 fccshifts(2, 1) = 0.0_cd
502 fccshifts(3, 1) = 0.0_cd
503 fccshifts(1, 2) = 0.5_cd * fccspacing
504 fccshifts(2, 2) = 0.5_cd * fccspacing
505 fccshifts(3, 2) = 0.0_cd
506 fccshifts(1, 3) = 0.5_cd * fccspacing
507 fccshifts(2, 3) = 0.0_cd
508 fccshifts(3, 3) = 0.5_cd * fccspacing
509 fccshifts(1, 4) = 0.0_cd
510 fccshifts(2, 4) = 0.5_cd * fccspacing
511 fccshifts(3, 4) = 0.5_cd * fccspacing
515 do i = 1, ncellsperside
516 latvec(1) = (i - 1) * fccspacing
517 do j = 1, ncellsperside
518 latvec(2) = (j - 1) * fccspacing
519 do k = 1, ncellsperside
520 latvec(3) = (k - 1) * fccspacing
523 coords(:, a) = latvec + fccshifts(:, m)
524 if ((i == ncellsperside / 2 + 1) &
525 .and. (j == ncellsperside / 2 + 1) &
526 .and. (k == ncellsperside / 2 + 1) &
530 coords(:, 1) = latvec + fccshifts(:, m)
535 if (.not. periodic)
then 538 latvec(1) = ncellsperside * fccspacing
539 latvec(2) = (i - 1) * fccspacing
540 latvec(3) = (j - 1) * fccspacing
541 a = a + 1; coords(:, a) = latvec
542 a = a + 1; coords(:, a) = latvec + fccshifts(:, 4)
544 latvec(1) = (i - 1) * fccspacing
545 latvec(2) = ncellsperside * fccspacing
546 latvec(3) = (j - 1) * fccspacing
547 a = a + 1; coords(:, a) = latvec
548 a = a + 1; coords(:, a) = latvec + fccshifts(:, 3)
550 latvec(1) = (i - 1) * fccspacing
551 latvec(2) = (j - 1) * fccspacing
552 latvec(3) = ncellsperside * fccspacing
553 a = a + 1; coords(:, a) = latvec
554 a = a + 1; coords(:, a) = latvec + fccshifts(:, 2)
557 if (.not. periodic)
then 559 latvec(1) = (i - 1) * fccspacing
560 latvec(2) = ncellsperside * fccspacing
561 latvec(3) = ncellsperside * fccspacing
562 a = a + 1; coords(:, a) = latvec
563 latvec(1) = ncellsperside * fccspacing
564 latvec(2) = (i - 1) * fccspacing
565 latvec(3) = ncellsperside * fccspacing
566 a = a + 1; coords(:, a) = latvec
567 latvec(1) = ncellsperside * fccspacing
568 latvec(2) = ncellsperside * fccspacing
569 latvec(3) = (i - 1) * fccspacing
570 a = a + 1; coords(:, a) = latvec
573 if (.not. periodic)
then 575 a = a + 1; coords(:, a) = ncellsperside * fccspacing
583 partnum, dir, model_handle, compute_arguments_handle, DIM, N, coords, &
584 cutoff, cutpad, energy, do_update_list, coordsave, neighObject, deriv, &
586 use,
intrinsic :: iso_c_binding
590 integer(c_int),
parameter :: cd = c_double
593 integer(c_int),
intent(in) :: partnum
594 integer(c_int),
intent(in) :: dir
595 type(kim_model_handle_type),
intent(in) :: model_handle
596 type(kim_compute_arguments_handle_type),
intent(in) :: &
597 compute_arguments_handle
598 integer(c_int),
intent(in) :: DIM
599 integer(c_int),
intent(in) :: N
600 real(c_double),
intent(inout) :: coords(dim, n)
601 real(c_double),
intent(in) :: cutoff
602 real(c_double),
intent(in) :: cutpad
603 real(c_double),
intent(inout) :: energy
604 logical,
intent(inout) :: do_update_list
605 real(c_double),
intent(inout) :: coordsave(dim, n)
607 real(c_double),
intent(out) :: deriv
608 real(c_double),
intent(out) :: deriv_err
609 integer(c_int),
intent(out) :: ierr
612 real(c_double),
parameter :: eps_init = 1.e-6_cd
613 integer(c_int),
parameter :: number_eps_levels = 15
614 real(c_double) eps, deriv_last, deriv_err_last
625 deriv_err_last = huge(1.0_cd)
626 do i = 1, number_eps_levels
627 deriv =
dfridr(eps, deriv_err)
629 call my_error(
"compute_numer_deriv")
631 if (deriv_err > deriv_err_last)
then 633 deriv_err = deriv_err_last
638 deriv_err_last = deriv_err
662 real(c_double) recursive function dfridr(h, err)
666 real(c_double),
intent(inout) :: h
667 real(c_double),
intent(out) :: err
670 integer(c_int),
parameter :: NTAB = 10
671 real(c_double),
parameter :: CON = 1.4_cd
672 real(c_double),
parameter :: CON2 = con * con
673 real(c_double),
parameter :: BIG = huge(1.0_cd)
674 real(c_double),
parameter :: SAFE = 2.0_cd
677 real(c_double) errt, fac, hh, a(ntab, ntab), fp, fm, coordorig
682 if (abs(h) <= tiny(0.0_cd))
then 688 coordorig = coords(dir, partnum)
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, 1) = (fp - fm) / (2.0_cd * hh)
717 coords(dir, partnum) = coordorig + hh
719 do_update_list, coordsave, &
721 call kim_compute(model_handle, compute_arguments_handle, ierr)
723 call my_error(
"kim_api_model_compute")
726 coords(dir, partnum) = coordorig - hh
728 do_update_list, coordsave, &
730 call kim_compute(model_handle, compute_arguments_handle, ierr)
732 call my_error(
"kim_api_model_compute")
735 coords(dir, partnum) = coordorig
737 do_update_list, coordsave, &
739 a(1, i) = (fp - fm) / (2.0_cd * hh)
744 a(j, i) = (a(j - 1, i) * fac - a(j - 1, i - 1)) / (fac - 1.0_cd)
748 errt = max(abs(a(j, i) - a(j - 1, i)), abs(a(j, i) - a(j - 1, i - 1)))
749 if (errt <= err)
then 754 if (abs(a(i, i) - a(i - 1, i - 1)) >= safe * err)
return 779 use,
intrinsic :: iso_c_binding
784 integer(c_int),
parameter :: cd = c_double
786 integer(c_int),
parameter :: nCellsPerSide = 2
787 integer(c_int),
parameter :: DIM = 3
788 real(c_double),
parameter :: cutpad = 0.75_cd
789 integer(c_int),
parameter :: max_species = 200
790 real(c_double),
parameter :: eps_prec = epsilon(1.0_cd)
791 real(c_double) FCCspacing
793 integer(c_int),
parameter :: N = 4 * (ncellsperside)**3 &
794 + 6 * (ncellsperside)**2 &
795 + 3 * (ncellsperside) + 1
796 real(c_double),
allocatable :: forces_num(:, :)
797 real(c_double),
allocatable :: forces_num_err(:, :)
798 type(kim_species_name_type) :: model_species(max_species)
799 integer(c_int),
target :: num_species
800 character(len=5, kind=c_char) :: passfail
801 real(c_double) :: forcediff
802 real(c_double) :: forcediff_sumsq
803 real(c_double) :: weight
804 real(c_double) :: weight_sum
805 real(c_double) :: alpha
806 real(c_double) :: term
807 real(c_double) :: term_max
808 real(c_double),
allocatable :: cluster_coords(:, :)
809 real(c_double),
allocatable :: cluster_disps(:, :)
810 type(kim_species_name_type),
allocatable :: cluster_species(:)
811 integer(c_int) I, J, Imax, Jmax, species_code
812 integer(c_int) seed_size
813 integer(c_int),
allocatable :: seed(:)
819 integer(c_int),
allocatable,
target :: neighborList(:, :)
820 real(c_double),
allocatable :: coordsave(:, :)
821 logical do_update_list
826 character(len=256, kind=c_char) :: testname =
"vc_forces_numer_deriv" 827 character(len=256, kind=c_char) :: modelname
829 type(kim_model_handle_type) :: model_handle
830 type(kim_compute_arguments_handle_type) :: compute_arguments_handle
831 integer(c_int) ierr, ierr2
832 integer(c_int) species_is_supported
833 integer(c_int),
target :: numberOfParticles
834 integer(c_int),
target :: particleSpeciesCodes(n)
835 integer(c_int),
target :: particleContributing(n)
836 real(c_double) :: influence_distance
837 integer(c_int) :: number_of_neighbor_lists
838 real(c_double),
allocatable :: cutoffs(:)
839 integer(c_int),
allocatable :: &
840 model_will_not_request_neighbors_of_noncontributing_particles(:)
841 real(c_double) :: cutoff
842 real(c_double),
target :: energy
843 real(c_double),
target :: particle_energy(n)
844 real(c_double),
target :: particle_energy_sum
845 real(c_double),
target :: coords(3, n)
846 real(c_double),
target :: forces_kim(3, n)
847 real(c_double) :: forces(3, n)
848 integer(c_int) middleDum
849 logical doing_non_contributing
850 logical particle_energy_supported
851 logical particle_energy_optional
852 logical noncontrib_particle_energy_nonzero
853 logical forces_optional
854 logical model_is_compatible
855 integer(c_int) number_of_model_routine_names
856 type(kim_model_routine_name_type) model_routine_name
857 integer(c_int) present
858 integer(c_int) required
859 integer(c_int) requested_units_accepted
860 real(c_double) rnd, deriv, deriv_err
861 real(c_double),
pointer :: null_pointer
863 nullify (null_pointer)
865 doing_non_contributing = .false.
867 numberofparticles = n
877 call random_seed(size=seed_size)
878 allocate (seed(seed_size))
880 call random_seed(put=seed)
884 print
'("Please enter a valid KIM model name: ")' 885 read (*, *) modelname
890 print *,
'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES' 893 print
'("This is Test : ",A)', trim(testname)
894 print
'("Results for KIM Model : ",A)', trim(modelname)
898 call kim_model_create(kim_numbering_one_based, &
900 kim_energy_unit_ev, &
902 kim_temperature_unit_k, &
905 requested_units_accepted, &
912 if (requested_units_accepted == 0)
then 913 call my_error(
"Must adapt to model units")
917 call kim_get_number_of_model_routine_names(number_of_model_routine_names)
918 do i = 1, number_of_model_routine_names
919 call kim_get_model_routine_name(i, model_routine_name, ierr)
921 call my_error(
"kim_get_model_routine_name")
923 call kim_is_routine_present(model_handle, model_routine_name,
present, &
926 call my_error(
"kim_is_routine_present")
929 if ((
present == 1) .and. (required == 1))
then 930 if (.not. ((model_routine_name == kim_model_routine_name_create) &
931 .or. (model_routine_name == &
932 kim_model_routine_name_compute_arguments_create) &
933 .or. (model_routine_name == kim_model_routine_name_compute) &
934 .or. (model_routine_name == kim_model_routine_name_refresh) &
935 .or. (model_routine_name == &
936 kim_model_routine_name_compute_arguments_destroy) &
937 .or. (model_routine_name == kim_model_routine_name_destroy))) &
939 call my_error(
"Unknown required ModelRoutineName found.")
945 call kim_compute_arguments_create(model_handle, &
946 compute_arguments_handle, ierr)
948 call my_error(
"kim_model_compute_arguments_create")
952 particle_energy_supported, &
953 particle_energy_optional, &
954 model_is_compatible, ierr)
956 call my_error(
"error checking compatibility")
958 if (.not. model_is_compatible)
then 959 call my_error(
"incompatibility reported by check_model_compatibility")
967 call my_error(
"Get_Model_Supported_Species")
971 allocate (cluster_coords(3, n), cluster_disps(3, n), cluster_species(n))
973 call random_number(rnd)
974 species_code = 1 + int(rnd * num_species)
975 cluster_species(i) = model_species(species_code)
981 cluster_coords, middledum)
986 call random_number(rnd)
987 cluster_disps(j, i) = 0.1_cd * (rnd - 0.5_cd)
993 call kim_set_argument_pointer( &
994 compute_arguments_handle, &
995 kim_compute_argument_name_number_of_particles, &
996 numberofparticles, ierr2)
998 call kim_set_argument_pointer( &
999 compute_arguments_handle, &
1000 kim_compute_argument_name_particle_species_codes, &
1001 particlespeciescodes, ierr2)
1003 call kim_set_argument_pointer( &
1004 compute_arguments_handle, &
1005 kim_compute_argument_name_particle_contributing, &
1006 particlecontributing, ierr2)
1008 call kim_set_argument_pointer( &
1009 compute_arguments_handle, &
1010 kim_compute_argument_name_coordinates, coords, ierr2)
1012 call kim_set_argument_pointer( &
1013 compute_arguments_handle, &
1014 kim_compute_argument_name_partial_energy, energy, ierr2)
1016 if (particle_energy_supported)
then 1017 call kim_set_argument_pointer( &
1018 compute_arguments_handle, &
1019 kim_compute_argument_name_partial_particle_energy, particle_energy, ierr2)
1022 call kim_set_argument_pointer( &
1023 compute_arguments_handle, &
1024 kim_compute_argument_name_partial_forces, forces_kim, ierr2)
1027 call my_error(
"set_argument_pointer")
1033 allocate (coordsave(dim, n))
1034 allocate (neighborlist(n + 1, n))
1035 neighobject%neighbor_list_pointer = c_loc(neighborlist)
1036 neighobject%number_of_particles = n
1037 allocate (forces_num(dim, n), forces_num_err(dim, n))
1041 call kim_set_callback_pointer( &
1042 compute_arguments_handle, &
1043 kim_compute_callback_name_get_neighbor_list, kim_language_name_fortran, &
1044 c_funloc(
get_neigh), c_loc(neighobject), ierr)
1046 call my_error(
"set_callback_pointer")
1049 call kim_get_influence_distance(model_handle, influence_distance)
1050 call kim_get_number_of_neighbor_lists(model_handle, &
1051 number_of_neighbor_lists)
1052 allocate (cutoffs(number_of_neighbor_lists), &
1053 model_will_not_request_neighbors_of_noncontributing_particles( &
1054 number_of_neighbor_lists))
1055 call kim_get_neighbor_list_values( &
1056 model_handle, cutoffs, &
1057 model_will_not_request_neighbors_of_noncontributing_particles, ierr)
1059 call my_error(
"get_neighbor_list_values")
1061 cutoff = maxval(cutoffs)
1064 fccspacing = 0.75_cd * cutoff
1067 cluster_coords(:, i) = fccspacing * cluster_coords(:, i)
1069 print
'("Using FCC lattice parameter: ",f12.5)', fccspacing
1074 if (doing_non_contributing)
then 1076 if (particle_energy_optional)
then 1077 call kim_set_argument_pointer( &
1078 compute_arguments_handle, &
1079 kim_compute_argument_name_partial_particle_energy, &
1080 particle_energy, ierr)
1082 call my_error(
"set_argument_pointer")
1088 if (forces_optional)
then 1089 call kim_set_argument_pointer( &
1090 compute_arguments_handle, kim_compute_argument_name_partial_forces, &
1093 call my_error(
"set_argument_pointer")
1099 call kim_get_species_support_and_code( &
1100 model_handle, cluster_species(i), species_is_supported, &
1101 particlespeciescodes(i), ierr)
1104 call my_error(
"kim_api_get_species_code")
1106 if (.not. doing_non_contributing)
then 1108 particlecontributing(i) = 1
1113 call random_number(rnd)
1115 particlecontributing(i) = 1
1117 particlecontributing(i) = 0
1122 coords(:, i) = cluster_coords(:, i) + cluster_disps(:, i)
1127 do_update_list = .true.
1129 do_update_list, coordsave, &
1132 call my_error(
"update_neighborlist")
1137 call kim_compute(model_handle, compute_arguments_handle, ierr)
1139 call my_error(
"kim_api_model_compute")
1147 noncontrib_particle_energy_nonzero = .false.
1148 if (particle_energy_supported)
then 1149 particle_energy_sum = 0.0_cd
1151 if (particlecontributing(i) == 0)
then 1152 if (abs(particle_energy(i)) > epsilon(0.0_cd))
then 1153 noncontrib_particle_energy_nonzero = .true.
1156 particle_energy_sum = particle_energy_sum + particle_energy(i)
1164 if (.not. doing_non_contributing)
then 1165 print
'("Configuration with all contributing particles")' 1167 print
'("Configuration with some non-contributing particles")' 1169 if (particle_energy_supported)
then 1170 print
'(A25,2X,A25,2X,A15)',
"Energy",
"Sum Contrib. Energies",
"Diff" 1171 print
'(ES25.15,2X,ES25.15,2X,ES15.5)', energy, particle_energy_sum, &
1172 energy - particle_energy_sum
1173 if (noncontrib_particle_energy_nonzero)
then 1175 "Some non-contributing particles have non-zero & 1176 &partial_particle_energy")
1179 print
'("Energy = ",ES25.15)', energy
1185 if (particle_energy_optional)
then 1186 call kim_set_argument_pointer( &
1187 compute_arguments_handle, &
1188 kim_compute_argument_name_partial_particle_energy, &
1191 call my_error(
"set_argument_pointer")
1197 if (forces_optional)
then 1198 call kim_set_argument_pointer( &
1199 compute_arguments_handle, kim_compute_argument_name_partial_forces, &
1202 call my_error(
"set_argument_pointer")
1211 dim, n, coords, cutoff, &
1212 cutpad, energy, do_update_list, &
1213 coordsave, neighobject, deriv, deriv_err, ierr)
1215 call my_error(
"compute_numer_deriv")
1217 forces_num(j, i) = -deriv
1218 forces_num_err(j, i) = deriv_err
1224 print
'(A6,2X,A7,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)',
"Part",
"Contrib", &
1225 "Spec",
"Dir",
"Force_model",
"Force_numer",
"Force diff",
"pred error", &
1227 forcediff_sumsq = 0.0_cd
1231 forcediff = abs(forces(j, i) - forces_num(j, i))
1232 if (forcediff < forces_num_err(j, i))
then 1237 weight = max(abs(forces_num(j, i)), eps_prec) / &
1238 max(abs(forces_num_err(j, i)), eps_prec)
1239 term = weight * forcediff**2
1240 if (term > term_max)
then 1245 forcediff_sumsq = forcediff_sumsq + term
1246 weight_sum = weight_sum + weight
1248 print
'(I6,2X,I7,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1249 i, particlecontributing(i), &
1250 particlespeciescodes(i), j, forces(j, i), &
1251 forces_num(j, i), forcediff, &
1252 forces_num_err(j, i), weight, passfail
1254 print
'(23X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
1255 j, forces(j, i), forces_num(j, i), &
1256 forcediff, forces_num_err(j, i), weight, passfail
1261 alpha = sqrt(forcediff_sumsq / weight_sum) / dble(dim * n)
1263 print
'("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," & 1264 &(units of force)")', alpha
1266 print
'(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,'// &
1267 ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
1268 imax, jmax, abs(forces(jmax, imax) - forces_num(jmax, imax)), &
1269 abs(forces(jmax, imax) - forces_num(jmax, imax)) / abs(forces(jmax, imax))
1271 if (.not. doing_non_contributing)
then 1272 doing_non_contributing = .true.
1280 deallocate (forces_num)
1281 deallocate (forces_num_err)
1282 deallocate (neighborlist)
1283 deallocate (coordsave)
1284 deallocate (cutoffs)
1285 deallocate (model_will_not_request_neighbors_of_noncontributing_particles)
1287 call kim_compute_arguments_destroy(model_handle, &
1288 compute_arguments_handle, ierr)
1290 call my_error(
"kim_model_compute_arguments_destroy")
1292 call kim_model_destroy(model_handle)
1297 print
'(120(''-''))' 1301 deallocate (cluster_coords, cluster_disps, cluster_species)
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 check_model_compatibility(compute_arguments_handle, forces_optional, particle_energy_supported, particle_energy_optional, model_is_compatible, ierr)
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)