42 use,
intrinsic :: iso_c_binding
62 integer(c_int),
parameter :: cd = c_double
63 integer(c_int),
parameter :: dim = 3
71 type, bind(c) :: buffer_type
72 character(c_char) :: species_name(100)
73 real(c_double) :: influence_distance(1)
74 real(c_double) :: pcutoff(1)
75 real(c_double) :: cutsq(1)
77 model_will_not_request_neighbors_of_noncontributing_particles(1)
78 real(c_double) :: epsilon(1)
79 real(c_double) :: sigma(1)
80 real(c_double) :: shift(1)
90 recursive subroutine calc_phi(model_epsilon, &
97 real(c_double),
intent(in) :: model_epsilon
98 real(c_double),
intent(in) :: model_sigma
99 real(c_double),
intent(in) :: model_shift
100 real(c_double),
intent(in) :: model_cutoff
101 real(c_double),
intent(in) :: r
102 real(c_double),
intent(out) :: phi
105 real(c_double) rsq, sor, sor6, sor12
108 sor = model_sigma / r
109 sor6 = sor * sor * sor
112 if (r > model_cutoff)
then 116 phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
129 model_cutoff, r, phi, dphi)
133 real(c_double),
intent(in) :: model_epsilon
134 real(c_double),
intent(in) :: model_sigma
135 real(c_double),
intent(in) :: model_shift
136 real(c_double),
intent(in) :: model_cutoff
137 real(c_double),
intent(in) :: r
138 real(c_double),
intent(out) :: phi, dphi
141 real(c_double) rsq, sor, sor6, sor12
144 sor = model_sigma / r
145 sor6 = sor * sor * sor
148 if (r > model_cutoff)
then 153 phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
154 dphi = 24.0_cd * model_epsilon * (-2.0_cd * sor12 + sor6) / r
167 model_cutoff, r, phi, dphi, d2phi)
171 real(c_double),
intent(in) :: model_epsilon
172 real(c_double),
intent(in) :: model_sigma
173 real(c_double),
intent(in) :: model_shift
174 real(c_double),
intent(in) :: model_cutoff
175 real(c_double),
intent(in) :: r
176 real(c_double),
intent(out) :: phi, dphi, d2phi
179 real(c_double) rsq, sor, sor6, sor12
182 sor = model_sigma / r
183 sor6 = sor * sor * sor
186 if (r > model_cutoff)
then 192 phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
193 dphi = 24.0_cd * model_epsilon * (-2.0_cd * sor12 + sor6) / r
194 d2phi = 24.0_cd * model_epsilon * (26.0_cd * sor12 - 7.0_cd * sor6) / rsq
205 model_compute_handle, model_compute_arguments_handle, ierr) bind(c)
209 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
210 type(kim_model_compute_arguments_handle_type),
intent(in) :: &
211 model_compute_arguments_handle
212 integer(c_int),
intent(out) :: ierr
215 real(c_double) :: r, Rsqij, phi, dphi, d2phi, dEidr, d2Eidr
216 integer(c_int) :: i, j, jj, numnei
217 integer(c_int) :: ierr2
218 integer(c_int) :: comp_force, comp_energy, comp_enepot, comp_process_dEdr, &
220 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
222 real(c_double),
target :: Rij(dim)
223 real(c_double),
target :: Rij_pairs(dim, 2)
224 real(c_double),
target :: r_pairs(2)
225 integer(c_int),
target :: i_pairs(2), j_pairs(2)
228 real(c_double) :: model_cutoff
229 integer(c_int),
pointer :: N
230 real(c_double),
pointer :: energy
231 real(c_double),
pointer :: coor(:, :)
232 real(c_double),
pointer :: force(:, :)
233 real(c_double),
pointer :: enepot(:)
234 integer(c_int),
pointer :: nei1part(:)
235 integer(c_int),
pointer :: particleSpeciesCodes(:)
236 integer(c_int),
pointer :: particleContributing(:)
239 call kim_get_model_buffer_pointer(model_compute_handle, pbuf)
240 call c_f_pointer(pbuf, buf)
242 model_cutoff = buf%influence_distance(1)
248 call kim_is_callback_present( &
249 model_compute_arguments_handle, &
250 kim_compute_callback_name_process_dedr_term, comp_process_dedr, ierr2)
252 call kim_is_callback_present( &
253 model_compute_arguments_handle, &
254 kim_compute_callback_name_process_d2edr2_term, comp_process_d2edr2, ierr2)
257 call kim_log_entry(model_compute_arguments_handle, &
258 kim_log_verbosity_error,
"get_compute")
266 call kim_get_argument_pointer( &
267 model_compute_arguments_handle, &
268 kim_compute_argument_name_number_of_particles, n, ierr2)
271 call kim_get_argument_pointer( &
272 model_compute_arguments_handle, &
273 kim_compute_argument_name_particle_species_codes, &
274 n, particlespeciescodes, ierr2)
276 call kim_get_argument_pointer( &
277 model_compute_arguments_handle, &
278 kim_compute_argument_name_particle_contributing, &
279 n, particlecontributing, ierr2)
281 call kim_get_argument_pointer( &
282 model_compute_arguments_handle, &
283 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
285 call kim_get_argument_pointer( &
286 model_compute_arguments_handle, &
287 kim_compute_argument_name_partial_energy, energy, ierr2)
289 call kim_get_argument_pointer( &
290 model_compute_arguments_handle, &
291 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
293 call kim_get_argument_pointer( &
294 model_compute_arguments_handle, &
295 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
298 call kim_log_entry(model_compute_arguments_handle, &
299 kim_log_verbosity_error,
"get_argument_pointer")
304 if (
associated(energy))
then 309 if (
associated(force))
then 314 if (
associated(enepot))
then 325 if (particlespeciescodes(i) /=
speccode)
then 326 call kim_log_entry( &
327 model_compute_arguments_handle, kim_log_verbosity_error, &
328 "Unexpected species code detected")
337 if (comp_enepot == 1) enepot = 0.0_cd
338 if (comp_energy == 1) energy = 0.0_cd
339 if (comp_force == 1) force = 0.0_cd
349 if (particlecontributing(i) == 1)
then 352 call kim_get_neighbor_list( &
353 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
356 call kim_log_entry( &
357 model_compute_arguments_handle, kim_log_verbosity_error, &
371 rij(:) = coor(:, j) - coor(:, i)
375 rsqij = dot_product(rij, rij)
376 if (rsqij < buf%cutsq(1))
then 379 if (comp_process_d2edr2 == 1)
then 386 deidr = 0.5_cd * dphi
387 d2eidr = 0.5_cd * d2phi
388 elseif (comp_force == 1 .or. comp_process_dedr == 1)
then 396 deidr = 0.5_cd * dphi
407 if (comp_enepot == 1)
then 408 enepot(i) = enepot(i) + 0.5_cd * phi
410 if (comp_energy == 1)
then 411 energy = energy + 0.5_cd * phi
416 if (comp_process_dedr == 1)
then 417 call kim_process_dedr_term( &
418 model_compute_arguments_handle, deidr, r, rij, i, j, ierr)
422 if (comp_process_d2edr2 == 1)
then 425 rij_pairs(:, 1) = rij
426 rij_pairs(:, 2) = rij
432 call kim_process_d2edr2_term( &
433 model_compute_arguments_handle, d2eidr, &
434 r_pairs, rij_pairs, i_pairs, j_pairs, ierr)
439 if (comp_force == 1)
then 440 force(:, i) = force(:, i) + deidr * rij / r
441 force(:, j) = force(:, j) - deidr * rij / r
464 recursive subroutine refresh(model_refresh_handle, ierr) bind(c)
468 type(kim_model_refresh_handle_type),
intent(in) :: model_refresh_handle
469 integer(c_int),
intent(out) :: ierr
472 real(c_double) energy_at_cutoff
473 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
476 call kim_get_model_buffer_pointer(model_refresh_handle, pbuf)
477 call c_f_pointer(pbuf, buf)
479 call kim_set_influence_distance_pointer(model_refresh_handle, &
480 buf%influence_distance(1))
481 call kim_set_neighbor_list_pointers( &
482 model_refresh_handle, 1, buf%influence_distance, &
483 buf%model_will_not_request_neighbors_of_noncontributing_particles)
487 buf%influence_distance(1) = buf%Pcutoff(1)
488 buf%cutsq(1) = (buf%Pcutoff(1))**2
494 buf%Pcutoff(1), energy_at_cutoff)
495 buf%shift(1) = -energy_at_cutoff
508 model_write_parameterized_model_handle, ierr) bind(c)
512 type(kim_model_write_parameterized_model_handle_type),
intent(in) &
513 :: model_write_parameterized_model_handle
514 integer(c_int),
intent(out) :: ierr
518 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
519 character(len=512, kind=c_char) :: path
520 character(len=512, kind=c_char) :: model_name
521 character(len=512, kind=c_char) :: string_buffer
522 character(len=100, kind=c_char) :: species_name
525 call kim_get_model_buffer_pointer( &
526 model_write_parameterized_model_handle, pbuf)
527 call c_f_pointer(pbuf, buf)
529 call kim_get_path(model_write_parameterized_model_handle, path)
530 call kim_get_model_name(model_write_parameterized_model_handle, model_name)
532 write (string_buffer,
'(A)') trim(model_name)//
".params" 533 call kim_set_parameter_file_name(model_write_parameterized_model_handle, &
535 write (string_buffer,
'(A)') trim(path)//
"/"//trim(string_buffer)
537 open (42, file=trim(string_buffer), &
538 status=
"REPLACE", action=
"WRITE", iostat=ierr)
540 call kim_log_entry( &
541 model_write_parameterized_model_handle, kim_log_verbosity_error, &
542 "Unable to open parameter file for writing.")
547 species_name(i:i) = buf%species_name(i)
549 write (42,
'(A)') trim(species_name)
550 write (42,
'(ES20.10)') buf%Pcutoff(1)
551 write (42,
'(ES20.10)') buf%epsilon(1)
552 write (42,
'(ES20.10)') buf%sigma(1)
564 recursive subroutine destroy(model_destroy_handle, ierr) bind(c)
568 type(kim_model_destroy_handle_type),
intent(in) :: model_destroy_handle
569 integer(c_int),
intent(out) :: ierr
572 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
575 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
576 call c_f_pointer(pbuf, buf)
591 model_compute_handle, model_compute_arguments_create_handle, ierr) bind(c)
595 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
596 type(kim_model_compute_arguments_create_handle_type),
intent(in) :: &
597 model_compute_arguments_create_handle
598 integer(c_int),
intent(out) :: ierr
603 if (model_compute_handle == kim_model_compute_null_handle)
continue 609 call kim_set_argument_support_status( &
610 model_compute_arguments_create_handle, &
611 kim_compute_argument_name_partial_energy, &
612 kim_support_status_optional, ierr)
613 call kim_set_argument_support_status( &
614 model_compute_arguments_create_handle, &
615 kim_compute_argument_name_partial_forces, &
616 kim_support_status_optional, ierr2)
618 call kim_set_argument_support_status( &
619 model_compute_arguments_create_handle, &
620 kim_compute_argument_name_partial_particle_energy, &
621 kim_support_status_optional, ierr2)
624 call kim_log_entry( &
625 model_compute_arguments_create_handle, kim_log_verbosity_error, &
626 "Unable to register arguments support_statuss")
632 call kim_set_callback_support_status( &
633 model_compute_arguments_create_handle, &
634 kim_compute_callback_name_process_dedr_term, &
635 kim_support_status_optional, ierr)
636 call kim_set_callback_support_status( &
637 model_compute_arguments_create_handle, &
638 kim_compute_callback_name_process_d2edr2_term, &
639 kim_support_status_optional, ierr2)
642 call kim_log_entry( &
643 model_compute_arguments_create_handle, kim_log_verbosity_error, &
644 "Unable to register callbacks support_statuss")
661 model_compute_handle, model_compute_arguments_destroy_handle, ierr) bind(c)
665 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
666 type(kim_model_compute_arguments_destroy_handle_type),
intent(in) :: &
667 model_compute_arguments_destroy_handle
668 integer(c_int),
intent(out) :: ierr
671 if (model_compute_handle == kim_model_compute_null_handle)
continue 672 if (model_compute_arguments_destroy_handle == &
673 kim_model_compute_arguments_destroy_null_handle)
continue 686 recursive subroutine model_driver_create_routine( &
687 model_driver_create_handle, requested_length_unit, requested_energy_unit, &
688 requested_charge_unit, requested_temperature_unit, requested_time_unit, &
690 use,
intrinsic :: iso_c_binding
694 integer(c_int),
parameter :: cd = c_double
697 type(kim_model_driver_create_handle_type),
intent(in) &
698 :: model_driver_create_handle
699 type(kim_length_unit_type),
intent(in),
value :: requested_length_unit
700 type(kim_energy_unit_type),
intent(in),
value :: requested_energy_unit
701 type(kim_charge_unit_type),
intent(in),
value :: requested_charge_unit
702 type(kim_temperature_unit_type),
intent(in),
value :: &
703 requested_temperature_unit
704 type(kim_time_unit_type),
intent(in),
value :: requested_time_unit
705 integer(c_int),
intent(out) :: ierr
709 integer(c_int) :: number_of_parameter_files
710 character(len=1024, kind=c_char) :: parameter_file_directory_name
711 character(len=1024, kind=c_char) :: parameter_file_basename
712 character(len=1024, kind=c_char) :: parameter_file_name
713 integer(c_int) :: ierr2
714 type(buffer_type),
pointer :: buf
715 type(kim_species_name_type) species_name
717 real(c_double) factor
718 character(len=100, kind=c_char) in_species
719 real(c_double) in_cutoff
720 real(c_double) in_epsilon
721 real(c_double) in_sigma
722 real(c_double) energy_at_cutoff
725 call kim_set_units( &
726 model_driver_create_handle, &
727 requested_length_unit, &
728 requested_energy_unit, &
729 kim_charge_unit_unused, &
730 kim_temperature_unit_unused, &
731 kim_time_unit_unused, ierr)
733 call kim_log_entry(model_driver_create_handle, &
734 kim_log_verbosity_error,
"Unable to set units")
740 call kim_set_model_numbering( &
741 model_driver_create_handle, kim_numbering_one_based, ierr)
743 call kim_log_entry(model_driver_create_handle, &
744 kim_log_verbosity_error,
"Unable to set numbering")
750 call kim_set_routine_pointer( &
751 model_driver_create_handle, kim_model_routine_name_compute, &
753 call kim_set_routine_pointer( &
754 model_driver_create_handle, &
755 kim_model_routine_name_compute_arguments_create, &
758 call kim_set_routine_pointer( &
759 model_driver_create_handle, kim_model_routine_name_refresh, &
760 kim_language_name_fortran, 1, c_funloc(
refresh), ierr2)
762 call kim_set_routine_pointer( &
763 model_driver_create_handle, &
764 kim_model_routine_name_write_parameterized_model, &
765 kim_language_name_fortran, 0, c_funloc(
write_model), ierr2)
767 call kim_set_routine_pointer( &
768 model_driver_create_handle, &
769 kim_model_routine_name_compute_arguments_destroy, &
772 call kim_set_routine_pointer( &
773 model_driver_create_handle, kim_model_routine_name_destroy, &
774 kim_language_name_fortran, 1, c_funloc(
destroy), ierr2)
777 call kim_log_entry( &
778 model_driver_create_handle, kim_log_verbosity_error, &
779 "Unable to store callback pointers")
785 call kim_get_number_of_parameter_files( &
786 model_driver_create_handle, number_of_parameter_files)
787 if (number_of_parameter_files /= 1)
then 788 call kim_log_entry( &
789 model_driver_create_handle, kim_log_verbosity_error, &
790 "Wrong number of parameter files")
797 call kim_get_parameter_file_directory_name( &
798 model_driver_create_handle, parameter_file_directory_name)
799 call kim_get_parameter_file_basename( &
800 model_driver_create_handle, 1, parameter_file_basename, ierr)
802 call kim_log_entry( &
803 model_driver_create_handle, kim_log_verbosity_error, &
804 "Unable to get parameter file basename")
808 parameter_file_name = trim(parameter_file_directory_name)//
"/"// &
809 trim(parameter_file_basename)
810 open (10, file=parameter_file_name, status=
"old")
811 read (10, *, iostat=ierr, err=100) in_species
812 read (10, *, iostat=ierr, err=100) in_cutoff
813 read (10, *, iostat=ierr, err=100) in_epsilon
814 read (10, *, iostat=ierr, err=100) in_sigma
820 call kim_log_entry(model_driver_create_handle, &
821 kim_log_verbosity_error,
"Unable to read LJ parameters")
828 call kim_from_string(in_species, species_name)
829 call kim_set_species_code( &
830 model_driver_create_handle, species_name,
speccode, ierr)
832 call kim_log_entry(model_driver_create_handle, &
833 kim_log_verbosity_error,
"Unable to set species code")
839 call kim_convert_unit( &
841 kim_energy_unit_ev, &
843 kim_temperature_unit_k, &
845 requested_length_unit, &
846 requested_energy_unit, &
847 requested_charge_unit, &
848 requested_temperature_unit, &
849 requested_time_unit, &
850 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
852 call kim_log_entry(model_driver_create_handle, &
853 kim_log_verbosity_error,
"kim_api_convert_to_act_unit")
857 in_cutoff = in_cutoff * factor
859 call kim_convert_unit( &
861 kim_energy_unit_ev, &
863 kim_temperature_unit_k, &
865 requested_length_unit, &
866 requested_energy_unit, &
867 requested_charge_unit, &
868 requested_temperature_unit, &
869 requested_time_unit, &
870 0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
872 call kim_log_entry(model_driver_create_handle, &
873 kim_log_verbosity_error,
"kim_api_convert_to_act_unit")
877 in_epsilon = in_epsilon * factor
879 call kim_convert_unit( &
881 kim_energy_unit_ev, &
883 kim_temperature_unit_k, &
885 requested_length_unit, &
886 requested_energy_unit, &
887 requested_charge_unit, &
888 requested_temperature_unit, &
889 requested_time_unit, &
890 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
892 call kim_log_entry(model_driver_create_handle, &
893 kim_log_verbosity_error,
"kim_api_convert_to_act_unit")
897 in_sigma = in_sigma * factor
904 buf%species_name(i) = in_species(i:i)
906 buf%influence_distance(1) = in_cutoff
907 buf%Pcutoff(1) = in_cutoff
908 buf%cutsq(1) = in_cutoff**2
909 buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
910 buf%epsilon(1) = in_epsilon
911 buf%sigma(1) = in_sigma
916 in_cutoff, energy_at_cutoff)
917 buf%shift(1) = -energy_at_cutoff
920 call kim_set_influence_distance_pointer( &
921 model_driver_create_handle, buf%influence_distance(1))
922 call kim_set_neighbor_list_pointers( &
923 model_driver_create_handle, 1, buf%influence_distance, &
924 buf%model_will_not_request_neighbors_of_noncontributing_particles)
929 call kim_set_model_buffer_pointer( &
930 model_driver_create_handle, c_loc(buf))
933 call kim_set_parameter_pointer( &
934 model_driver_create_handle, buf%pcutoff,
"cutoff", &
935 "Distance beyond which particles do not interact with one another.", ierr)
937 call kim_log_entry(model_driver_create_handle, &
938 kim_log_verbosity_error,
"set_parameter")
943 call kim_set_parameter_pointer( &
944 model_driver_create_handle, buf%epsilon,
"epsilon", &
945 "Maximum depth of the potential well.", ierr)
947 call kim_log_entry(model_driver_create_handle, &
948 kim_log_verbosity_error,
"set_parameter")
953 call kim_set_parameter_pointer( &
954 model_driver_create_handle, buf%sigma,
"sigma", &
955 "Distance at which energy is zero and force is repulsive.", ierr)
957 call kim_log_entry(model_driver_create_handle, &
958 kim_log_verbosity_error,
"set_parameter")
967 end subroutine model_driver_create_routine
recursive subroutine, public compute_energy_forces(model_compute_handle, model_compute_arguments_handle, ierr)
recursive subroutine, public destroy(model_destroy_handle, ierr)
recursive subroutine, public compute_arguments_destroy(model_compute_handle, model_compute_arguments_destroy_handle, ierr)
recursive subroutine, public refresh(model_refresh_handle, ierr)
recursive subroutine, public calc_phi_dphi_d2phi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi, dphi, d2phi)
recursive subroutine, public calc_phi_dphi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi, dphi)
recursive subroutine, public compute_arguments_create(model_compute_handle, model_compute_arguments_create_handle, ierr)
integer(c_int), parameter, public speccode
recursive subroutine, public calc_phi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi)
recursive subroutine, public write_model(model_write_parameterized_model_handle, ierr)