41 subroutine my_error(message, line, file)
43 character(len=*),
intent(in) :: message
44 integer,
intent(in) :: line
45 character(len=*),
intent(in) :: file
47 print *,
"* Error : '", trim(message),
"' ", line,
":", &
54 character(len=*),
intent(in) :: message
55 integer,
intent(in) :: line
56 character(len=*),
intent(in) :: file
58 print *,
"* Warning : '", trim(message),
"' ", line,
":", &
74 use,
intrinsic :: iso_c_binding
79 integer(c_int) :: number_of_particles
80 integer(c_int),
pointer :: neighborList(:,:)
81 real(c_double),
pointer :: RijList(:,:,:)
92 subroutine get_neigh(data_object, neighbor_list_index, request, numnei, &
93 pnei1part, ierr) bind(c)
98 type(c_ptr),
value,
intent(in) :: data_object
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),
parameter :: DIM = 3
107 integer(c_int) numberOfParticles
108 type(neighObject_type),
pointer :: neighObject
110 if (neighbor_list_index /= 1)
then 111 call my_warning(
"wrong list index", __line__, __file__)
116 call c_f_pointer(data_object, neighobject)
118 numberofparticles = neighobject%number_of_particles
120 if ( (request.gt.numberofparticles) .or. (request.lt.1))
then 121 call my_warning(
"Invalid part ID in get_neigh", &
128 numnei = neighobject%neighborList(1,request)
131 pnei1part = c_loc(neighobject%neighborList(2,request))
153 use,
intrinsic :: iso_c_binding
164 integer(c_int),
parameter :: cd = c_double
166 integer(c_int),
parameter :: ncellsperside = 2
167 integer(c_int),
parameter :: dim = 3
168 real(c_double),
parameter :: cutpad = 0.75_cd
169 integer(c_int),
parameter :: max_species = 200
170 real(c_double),
parameter :: eps_prec = epsilon(1.0_cd)
171 real(c_double) fccspacing
173 integer(c_int),
parameter :: &
174 n = 4*(ncellsperside)**3 + 6*(ncellsperside)**2 + 3*(ncellsperside) + 1
175 real(c_double),
allocatable :: forces_num(:,:)
176 real(c_double),
allocatable :: forces_num_err(:,:)
177 type(kim_species_name_type) :: model_species(max_species)
178 integer(c_int),
target :: num_species
179 character(len=5) :: passfail
180 real(c_double) :: forcediff
181 real(c_double) :: forcediff_sumsq
182 real(c_double) :: weight
183 real(c_double) :: weight_sum
184 real(c_double) :: alpha
185 real(c_double) :: term
186 real(c_double) :: term_max
187 real(c_double),
allocatable :: cluster_coords(:,:)
188 real(c_double),
allocatable :: cluster_disps(:,:)
189 type(kim_species_name_type),
allocatable :: cluster_species(:)
190 integer(c_int) i,j,imax,jmax,species_code
196 real(c_double),
allocatable :: coordsave(:,:)
197 logical do_update_list
202 character(len=256) :: testname =
"vc_forces_numer_deriv" 203 character(len=256) :: modelname
205 type(kim_model_handle_type) :: model_handle
206 integer(c_int) ierr, ierr2
207 integer(c_int) species_is_supported
208 integer(c_int),
target :: numberofparticles
209 integer(c_int),
target :: particlespeciescodes(n)
210 integer(c_int),
target :: particlecontributing(n)
211 real(c_double) :: influence_distance
212 integer(c_int) :: number_of_cutoffs
213 real(c_double) :: cutoffs(1)
214 real(c_double) :: cutoff
215 real(c_double),
target :: energy
216 real(c_double),
target :: coords(3,n)
217 real(c_double),
target :: forces(3,n)
218 integer(c_int) middledum
219 logical forces_optional
220 logical model_is_compatible
221 integer(c_int) requested_units_accepted
222 real(c_double) rnd, deriv, deriv_err
223 real(c_double),
pointer :: null_pointer
225 nullify(null_pointer)
227 numberofparticles = n
235 print
'("Please enter a valid KIM model name: ")' 241 print *,
'VERIFICATION CHECK: NUMERICAL DERIVATIVE VERIFICATION OF FORCES' 244 print
'("This is Test : ",A)', trim(testname)
245 print
'("Results for KIM Model : ",A)', trim(modelname)
251 kim_energy_unit_ev, &
253 kim_temperature_unit_k, &
256 requested_units_accepted, &
259 call my_error(
"kim_api_create", __line__, __file__)
263 if (requested_units_accepted == 0)
then 264 call my_error(
"Must adapt to model units", __line__, __file__)
270 call my_error(
"error checking compatibility", __line__, __file__)
272 if (.not. model_is_compatible)
then 273 call my_error(
"incompatibility reported by check_model_compatibility", &
282 call my_error(
"Get_Model_Supported_Species", __line__, __file__)
286 allocate(cluster_coords(3,n),cluster_disps(3,n),cluster_species(n))
288 call random_number(rnd)
289 species_code = 1 + int(rnd*num_species)
290 cluster_species(i) = model_species(species_code)
296 cluster_coords, middledum)
301 call random_number(rnd)
302 cluster_disps(j,i) = 0.1_cd*(rnd-0.5_cd)
308 call kim_model_set_argument_pointer(model_handle, &
311 call kim_model_set_argument_pointer(model_handle, &
314 call kim_model_set_argument_pointer(model_handle, &
317 call kim_model_set_argument_pointer(model_handle, &
320 call kim_model_set_argument_pointer(model_handle, &
323 call kim_model_set_argument_pointer(model_handle, &
327 call my_error(
"set_argument_pointer", __line__, __file__)
333 allocate(neighobject%neighborList(n+1,n))
334 neighobject%number_of_particles = n
340 c_funloc(
get_neigh), c_loc(neighobject), ierr)
342 call my_error(
"set_callback_pointer", __line__, __file__)
345 call kim_model_get_influence_distance(model_handle, influence_distance)
346 call kim_model_get_number_of_neighbor_list_cutoffs(model_handle, number_of_cutoffs)
347 if (number_of_cutoffs /= 1)
then 348 call my_error(
"too many cutoffs", __line__, __file__)
350 call kim_model_get_neighbor_list_cutoffs(model_handle, cutoffs, ierr)
352 call my_error(
"get_cutoffs", __line__, __file__)
357 fccspacing = 0.75_cd*cutoff
360 cluster_coords(:,i) = fccspacing*cluster_coords(:,i)
362 print
'("Using FCC lattice parameter: ",f12.5)', fccspacing
365 call kim_model_get_species_support_and_code(model_handle, cluster_species(i), &
366 species_is_supported, particlespeciescodes(i), ierr)
369 call my_error(
"kim_api_get_species_code", __line__, __file__)
372 particlecontributing(i) = 1
375 coords(:,i) = cluster_coords(:,i) + cluster_disps(:,i)
380 do_update_list = .true.
381 allocate(coordsave(dim,n))
383 do_update_list,coordsave, &
386 call my_error(
"update_neighborlist", __line__, __file__)
393 call my_error(
"kim_api_model_compute", __line__, __file__)
399 print
'("Energy = ",ES25.15)', energy
405 if (forces_optional)
then 406 call kim_model_set_argument_pointer(model_handle, &
409 call my_error(
"set_argument_pointer", __line__, __file__)
415 allocate(forces_num(dim,n),forces_num_err(dim,n))
419 cutpad, energy, do_update_list, &
420 coordsave,neighobject,deriv,deriv_err,ierr)
422 call my_error(
"compute_numer_deriv", __line__, __file__)
424 forces_num(j,i) = -deriv
425 forces_num_err(j,i) = deriv_err
431 print
'(A6,2X,A4,2X,A3,2X,2A25,3A15,2X,A4)',
"Part",
"Spec",
"Dir", &
432 "Force_model",
"Force_numer",
"Force diff",
"pred error",
"weight", &
434 forcediff_sumsq = 0.0_cd
438 forcediff = abs(forces(j,i)-forces_num(j,i))
439 if (forcediff<forces_num_err(j,i))
then 444 weight = max(abs(forces_num(j,i)),eps_prec)/ &
445 max(abs(forces_num_err(j,i)),eps_prec)
446 term = weight*forcediff**2
447 if (term.gt.term_max)
then 452 forcediff_sumsq = forcediff_sumsq + term
453 weight_sum = weight_sum + weight
455 print
'(I6,2X,I4,2X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
456 i,particlespeciescodes(i),j,forces(j,i),forces_num(j,i), &
457 forcediff,forces_num_err(j,i),weight,passfail
459 print
'(14X,I3,2X,2ES25.15,3ES15.5,2X,A5)', &
460 j,forces(j,i),forces_num(j,i), &
461 forcediff,forces_num_err(j,i),weight,passfail
466 alpha = sqrt(forcediff_sumsq/weight_sum)/dble(dim*n)
468 print
'("alpha = |Force_model - Force_numer|_w/(DIM*N) = ",ES15.5," (units of force)")', &
471 print
'(''Maximum term obtained for Part = '',I6,'', Dir = '',I1,' // &
472 ''', forcediff = '',ES15.5, '', forcediff/force_model = '',ES15.5)', &
473 imax,jmax,abs(forces(jmax,imax)-forces_num(jmax,imax)), &
474 abs(forces(jmax,imax)-forces_num(jmax,imax))/abs(forces(jmax,imax))
478 deallocate(forces_num)
479 deallocate(forces_num_err)
480 deallocate(neighobject%neighborList)
481 deallocate(coordsave)
492 deallocate(cluster_coords,cluster_disps,cluster_species)
505 model_is_compatible, ierr)
506 use,
intrinsic :: iso_c_binding
515 type(kim_model_handle_type),
intent(in) :: model_handle
516 logical,
intent(out) :: forces_optional
517 logical,
intent(out) :: model_is_compatible
518 integer(c_int),
intent(out) :: ierr
522 integer(c_int) number_of_arguments
523 integer(c_int) number_of_callbacks
524 type(kim_argument_name_type) argument_name
525 type(kim_support_status_type) support_status
526 type(kim_callback_name_type) callback_name
530 model_is_compatible = .false.
534 call kim_argument_name_get_number_of_arguments(number_of_arguments)
535 do i=1,number_of_arguments
536 call kim_argument_name_get_argument_name(i, argument_name, ierr)
542 call kim_model_get_argument_support_status(model_handle, argument_name, &
543 support_status, ierr)
545 call my_warning(
"can't get argument support_status", &
554 call my_warning(
"unsupported required argument", &
564 call my_warning(
"model does not support energy", &
571 call my_warning(
"model does not support forces", &
576 forces_optional = .false.
578 forces_optional = .true.
580 call my_warning(
"unknown support_status for forces", &
589 call kim_callback_name_get_number_of_callbacks(number_of_callbacks)
590 do i=1,number_of_callbacks
591 call kim_callback_name_get_callback_name(i, callback_name, ierr)
598 support_status, ierr)
600 call my_warning(
"can't get call back support_status", &
607 call my_warning(
"unsupported required call back", &
615 model_is_compatible = .true.
628 use,
intrinsic :: iso_c_binding
635 type(kim_model_handle_type),
intent(in) :: model_handle
636 integer(c_int),
intent(in) :: max_species
637 type(kim_species_name_type),
intent(out) :: model_species(max_species)
638 integer(c_int),
intent(out) :: num_species
639 integer(c_int),
intent(out) :: ier
643 integer(c_int) total_num_species
644 type(kim_species_name_type) :: species_name
645 integer(c_int) species_is_supported
651 call kim_species_name_get_number_of_species(total_num_species)
653 if (total_num_species .gt. max_species)
return 656 do i=1,total_num_species
657 call kim_species_name_get_species_name(i, species_name, ier)
658 call kim_model_get_species_support_and_code(model_handle, species_name, &
659 species_is_supported, code, ier)
660 if ((ier == 0) .and. (species_is_supported .ne. 0))
then 661 num_species = num_species + 1
662 model_species(num_species) = species_name
672 do_update_list,coordsave, &
674 use,
intrinsic :: iso_c_binding
677 integer(c_int),
parameter :: cd = c_double
680 integer(c_int),
intent(in) :: DIM
681 integer(c_int),
intent(in) :: N
682 real(c_double),
intent(in) :: coords(DIM,N)
683 real(c_double),
intent(in) :: cutoff
684 real(c_double),
intent(in) :: cutpad
685 logical,
intent(inout) :: do_update_list
686 real(c_double),
intent(inout) :: coordsave(DIM,N)
687 type(neighObject_type),
intent(inout) :: neighObject
688 integer(c_int),
intent(out) :: ierr
693 real(c_double) disp, disp1, disp2, cutrange, dispvec(DIM)
702 if (.not.do_update_list)
then 709 dispvec(1:dim) = coords(1:dim,i) - coordsave(1:dim,i)
710 disp = sqrt( dot_product(dispvec,dispvec) )
711 if (disp >= disp1)
then 714 else if (disp >= disp2)
then 718 do_update_list = ( disp1 + disp2 > cutpad )
722 if (do_update_list)
then 725 coordsave(1:dim,1:n) = coords(1:dim,1:n)
728 cutrange = cutoff + cutpad
733 do_update_list = .false.
749 use,
intrinsic :: iso_c_binding
754 logical,
intent(in) :: half
755 integer(c_int),
intent(in) :: numberOfParticles
756 real(c_double),
dimension(3,numberOfParticles), &
758 real(c_double),
intent(in) :: cutoff
759 type(neighObject_type),
intent(inout) :: neighObject
762 integer(c_int) i, j, a
765 real(c_double) cutoff2
769 do i=1,numberofparticles
771 do j=1,numberofparticles
772 dx(:) = coords(:, j) - coords(:, i)
773 r2 = dot_product(dx, dx)
774 if (r2.le.cutoff2)
then 776 if ( (j .gt. i) .OR. ((.not. half) .AND. (i.ne.j)) )
then 778 neighobject%neighborList(a,i) = j
783 neighobject%neighborList(1,i) = a-1
808 coords, MiddlePartId)
809 use,
intrinsic :: iso_c_binding
811 integer(c_int),
parameter :: cd = c_double
814 real(c_double),
intent(in) :: FCCspacing
815 integer(c_int),
intent(in) :: nCellsPerSide
816 logical,
intent(in) :: periodic
817 real(c_double),
intent(out) :: coords(3,*)
818 integer(c_int),
intent(out) :: MiddlePartId
822 real(c_double) FCCshifts(3,4)
823 real(c_double) latVec(3)
824 integer(c_int) a, i, j, k, m
828 fccshifts(1,1) = 0.0_cd
829 fccshifts(2,1) = 0.0_cd
830 fccshifts(3,1) = 0.0_cd
831 fccshifts(1,2) = 0.5_cd*fccspacing
832 fccshifts(2,2) = 0.5_cd*fccspacing
833 fccshifts(3,2) = 0.0_cd
834 fccshifts(1,3) = 0.5_cd*fccspacing
835 fccshifts(2,3) = 0.0_cd
836 fccshifts(3,3) = 0.5_cd*fccspacing
837 fccshifts(1,4) = 0.0_cd
838 fccshifts(2,4) = 0.5_cd*fccspacing
839 fccshifts(3,4) = 0.5_cd*fccspacing
844 latvec(1) = (i-1)*fccspacing
846 latvec(2) = (j-1)*fccspacing
848 latvec(3) = (k-1)*fccspacing
851 coords(:,a) = latvec + fccshifts(:,m)
852 if ((i.eq.ncellsperside/2+1).and.(j.eq.ncellsperside/2+1) .and. &
853 (k.eq.ncellsperside/2+1) .and. (m.eq.1))
then 854 coords(:,1) = latvec + fccshifts(:,m)
859 if (.not. periodic)
then 862 latvec(1) = ncellsperside*fccspacing
863 latvec(2) = (i-1)*fccspacing
864 latvec(3) = (j-1)*fccspacing
865 a = a+1; coords(:,a) = latvec
866 a = a+1; coords(:,a) = latvec + fccshifts(:,4)
868 latvec(1) = (i-1)*fccspacing
869 latvec(2) = ncellsperside*fccspacing
870 latvec(3) = (j-1)*fccspacing
871 a = a+1; coords(:,a) = latvec
872 a = a+1; coords(:,a) = latvec + fccshifts(:,3)
874 latvec(1) = (i-1)*fccspacing
875 latvec(2) = (j-1)*fccspacing
876 latvec(3) = ncellsperside*fccspacing
877 a = a+1; coords(:,a) = latvec
878 a = a+1; coords(:,a) = latvec + fccshifts(:,2)
881 if (.not. periodic)
then 883 latvec(1) = (i-1)*fccspacing
884 latvec(2) = ncellsperside*fccspacing
885 latvec(3) = ncellsperside*fccspacing
886 a = a+1; coords(:,a) = latvec
887 latvec(1) = ncellsperside*fccspacing
888 latvec(2) = (i-1)*fccspacing
889 latvec(3) = ncellsperside*fccspacing
890 a = a+1; coords(:,a) = latvec
891 latvec(1) = ncellsperside*fccspacing
892 latvec(2) = ncellsperside*fccspacing
893 latvec(3) = (i-1)*fccspacing
894 a = a+1; coords(:,a) = latvec
897 if (.not. periodic)
then 899 a = a+1; coords(:,a) = ncellsperside*fccspacing
907 cutoff,cutpad,energy, do_update_list, &
908 coordsave,neighObject,deriv,deriv_err,ierr)
909 use,
intrinsic :: iso_c_binding
915 integer(c_int),
parameter :: cd = c_double
918 integer(c_int),
intent(in) :: partnum
919 integer(c_int),
intent(in) :: dir
920 type(kim_model_handle_type),
intent(in) :: model_handle
921 integer(c_int),
intent(in) :: DIM
922 integer(c_int),
intent(in) :: N
923 real(c_double),
intent(inout) :: coords(DIM,N)
924 real(c_double),
intent(in) :: cutoff
925 real(c_double),
intent(in) :: cutpad
926 real(c_double),
intent(inout) :: energy
927 logical,
intent(inout) :: do_update_list
928 real(c_double),
intent(inout) :: coordsave(DIM,N)
929 type(neighObject_type),
intent(inout) :: neighObject
930 real(c_double),
intent(out) :: deriv
931 real(c_double),
intent(out) :: deriv_err
932 integer(c_int),
intent(out) :: ierr
935 real(c_double),
parameter :: eps_init = 1.e-6_cd
936 integer(c_int),
parameter :: number_eps_levels = 15
937 real(c_double) eps, deriv_last, deriv_err_last
948 deriv_err_last = huge(1.0_cd)
949 do i=1,number_eps_levels
950 deriv =
dfridr(eps,deriv_err)
952 call my_error(
"compute_numer_deriv", &
955 if (deriv_err>deriv_err_last)
then 957 deriv_err = deriv_err_last
962 deriv_err_last = deriv_err
986 real(c_double) function dfridr(h,err)
990 real(c_double),
intent(inout) :: h
991 real(c_double),
intent(out) :: err
994 integer(c_int),
parameter :: ntab=10
995 real(c_double),
parameter :: con=1.4_cd
996 real(c_double),
parameter :: con2=con*con
997 real(c_double),
parameter :: big=huge(1.0_cd)
998 real(c_double),
parameter :: safe=2.0_cd
1001 real(c_double) errt,fac,hh,a(ntab,ntab),fp,fm,coordorig
1005 if (h.eq.0.0_cd)
then 1011 coordorig = coords(dir,partnum)
1012 coords(dir,partnum) = coordorig + hh
1014 do_update_list,coordsave, &
1018 call my_error(
"kim_api_model_compute", &
1022 coords(dir,partnum) = coordorig - hh
1024 do_update_list,coordsave, &
1028 call my_error(
"kim_api_model_compute", &
1032 coords(dir,partnum) = coordorig
1034 do_update_list,coordsave, &
1036 a(1,1)=(fp-fm)/(2.0_cd*hh)
1043 coords(dir,partnum) = coordorig + hh
1045 do_update_list,coordsave, &
1049 call my_error(
"kim_api_model_compute", &
1053 coords(dir,partnum) = coordorig - hh
1055 do_update_list,coordsave, &
1059 call my_error(
"kim_api_model_compute", &
1063 coords(dir,partnum) = coordorig
1065 do_update_list,coordsave, &
1067 a(1,i)=(fp-fm)/(2.0_cd*hh)
1072 a(j,i)=(a(j-1,i)*fac-a(j-1,i-1))/(fac-1.0_cd)
1076 errt = max(abs(a(j,i)-a(j-1,i)),abs(a(j,i)-a(j-1,i-1)))
1077 if (errt.le.err)
then 1082 if (abs(a(i,i)-a(i-1,i-1)).ge.safe*err)
return type(kim_numbering_type), public, protected kim_numbering_one_based
type(kim_argument_name_type), public, protected kim_argument_name_number_of_particles
type(kim_language_name_type), public, protected kim_language_name_fortran
type(kim_support_status_type), public, protected kim_support_status_optional
subroutine my_error(message, line, file)
subroutine get_model_supported_species(model_handle, max_species, model_species, num_species, ier)
subroutine update_neighborlist(DIM, N, coords, cutoff, cutpad, do_update_list, coordsave, neighObject, ierr)
type(kim_argument_name_type), public, protected kim_argument_name_partial_forces
type(kim_argument_name_type), public, protected kim_argument_name_particle_contributing
type(kim_callback_name_type), public, protected kim_callback_name_get_neighbor_list
type(kim_support_status_type), public, protected kim_support_status_required
subroutine check_model_compatibility(model_handle, forces_optional, model_is_compatible, ierr)
type(kim_argument_name_type), public, protected kim_argument_name_particle_species_codes
subroutine neigh_pure_cluster_neighborlist(half, numberOfParticles, coords, cutoff, neighObject)
program vc_forces_numer_deriv
type(kim_argument_name_type), public, protected kim_argument_name_partial_energy
type(kim_support_status_type), public, protected kim_support_status_not_supported
subroutine my_warning(message, line, file)
real(c_double) function dfridr(h, err)
type(kim_argument_name_type), public, protected kim_argument_name_coordinates
subroutine compute_numer_deriv(partnum, dir, model_handle, DIM, N, coords, cutoff, cutpad, energy, do_update_list, coordsave, neighObject, deriv, deriv_err, ierr)
subroutine, public get_neigh(data_object, neighbor_list_index, request, numnei, pnei1part, ierr)
subroutine create_fcc_configuration(FCCspacing, nCellsPerSide, periodic, coords, MiddlePartId)