46 use,
intrinsic :: iso_c_binding
64 integer(c_int),
parameter :: cd = c_double
65 integer(c_int),
parameter :: DIM=3
66 integer(c_int),
parameter :: speccode = 1
73 type, bind(c) :: buffer_type
74 real(c_double) :: influence_distance(1)
75 real(c_double) :: pcutoff(1)
76 real(c_double) :: cutsq(1)
77 integer(c_int) :: padding_neighbor_hints(1)
78 integer(c_int) :: half_list_hints(1)
79 real(c_double) :: epsilon(1)
80 real(c_double) :: sigma(1)
81 real(c_double) :: shift(1)
99 real(c_double),
intent(in) :: model_epsilon
100 real(c_double),
intent(in) :: model_sigma
101 real(c_double),
intent(in) :: model_shift
102 real(c_double),
intent(in) :: model_cutoff
103 real(c_double),
intent(in) :: r
104 real(c_double),
intent(out) :: phi
107 real(c_double) rsq,sor,sor6,sor12
114 if (r .gt. model_cutoff)
then 118 phi = 4.0_cd*model_epsilon*(sor12-sor6) + model_shift
131 model_cutoff,r,phi,dphi)
135 real(c_double),
intent(in) :: model_epsilon
136 real(c_double),
intent(in) :: model_sigma
137 real(c_double),
intent(in) :: model_shift
138 real(c_double),
intent(in) :: model_cutoff
139 real(c_double),
intent(in) :: r
140 real(c_double),
intent(out) :: phi,dphi
143 real(c_double) rsq,sor,sor6,sor12
150 if (r .gt. model_cutoff)
then 155 phi = 4.0_cd*model_epsilon*(sor12-sor6) + model_shift
156 dphi = 24.0_cd*model_epsilon*(-2.0_cd*sor12+sor6)/r
169 model_cutoff,r,phi,dphi,d2phi)
173 real(c_double),
intent(in) :: model_epsilon
174 real(c_double),
intent(in) :: model_sigma
175 real(c_double),
intent(in) :: model_shift
176 real(c_double),
intent(in) :: model_cutoff
177 real(c_double),
intent(in) :: r
178 real(c_double),
intent(out) :: phi,dphi,d2phi
181 real(c_double) rsq,sor,sor6,sor12
188 if (r .gt. model_cutoff)
then 194 phi = 4.0_cd*model_epsilon*(sor12-sor6) + model_shift
195 dphi = 24.0_cd*model_epsilon*(-2.0_cd*sor12+sor6)/r
196 d2phi = 24.0_cd*model_epsilon*(26.0_cd*sor12-7.0_cd*sor6)/rsq
206 #include "kim_model_compute_log_macros.fd" 208 model_compute_arguments_handle, ierr) bind(c)
212 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
213 type(kim_model_compute_arguments_handle_type),
intent(in) :: &
214 model_compute_arguments_handle
215 integer(c_int),
intent(out) :: ierr
218 real(c_double) :: r,rsqij,phi,dphi,d2phi,deidr,d2eidr
219 integer(c_int) :: i,j,jj,numnei
220 integer(c_int) :: ierr2
221 integer(c_int) :: comp_force,comp_energy,comp_enepot,comp_process_dedr, &
223 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
225 real(c_double),
pointer :: rij(:)
226 real(c_double),
pointer :: rij_pairs(:,:)
227 real(c_double),
pointer :: r_pairs(:)
228 integer(c_int),
pointer :: i_pairs(:), j_pairs(:)
231 real(c_double) :: model_cutoff
232 integer(c_int),
pointer :: n
233 real(c_double),
pointer :: energy
234 real(c_double),
pointer :: coor(:,:)
235 real(c_double),
pointer :: force(:,:)
236 real(c_double),
pointer :: enepot(:)
237 integer(c_int),
pointer :: nei1part(:)
238 integer(c_int),
pointer :: particlespeciescodes(:)
239 integer(c_int),
pointer :: particlecontributing(:)
241 kim_log_file = __file__
244 call kim_model_compute_get_model_buffer_pointer(model_compute_handle, pbuf)
245 call c_f_pointer(pbuf, buf)
247 model_cutoff = buf%influence_distance(1)
253 call kim_model_compute_arguments_is_callback_present( &
254 model_compute_arguments_handle, &
255 kim_compute_callback_name_process_dedr_term, comp_process_dedr, ierr2)
257 call kim_model_compute_arguments_is_callback_present( &
258 model_compute_arguments_handle, &
259 kim_compute_callback_name_process_d2edr2_term, comp_process_d2edr2, ierr2)
262 kim_log_message =
"get_compute" 270 call kim_model_compute_arguments_get_argument_pointer( &
271 model_compute_arguments_handle, &
272 kim_compute_argument_name_number_of_particles, n, ierr2)
275 call kim_model_compute_arguments_get_argument_pointer( &
276 model_compute_arguments_handle, &
277 kim_compute_argument_name_particle_species_codes, &
278 n, particlespeciescodes, ierr2)
280 call kim_model_compute_arguments_get_argument_pointer( &
281 model_compute_arguments_handle, &
282 kim_compute_argument_name_particle_contributing, n, particlecontributing, &
285 call kim_model_compute_arguments_get_argument_pointer( &
286 model_compute_arguments_handle, &
287 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
289 call kim_model_compute_arguments_get_argument_pointer( &
290 model_compute_arguments_handle, &
291 kim_compute_argument_name_partial_energy, energy, ierr2)
293 call kim_model_compute_arguments_get_argument_pointer( &
294 model_compute_arguments_handle, &
295 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
297 call kim_model_compute_arguments_get_argument_pointer( &
298 model_compute_arguments_handle, &
299 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
302 kim_log_message =
"get_argument_pointer" 307 if (
associated(energy))
then 312 if (
associated(force))
then 317 if (
associated(enepot))
then 324 if (comp_process_d2edr2.eq.1)
then 325 allocate( r_pairs(2) )
326 allocate( rij_pairs(dim,2) )
327 allocate( i_pairs(2) )
328 allocate( j_pairs(2) )
336 if (particlespeciescodes(i).ne.speccode)
then 337 kim_log_message =
"Unexpected species code detected" 346 if (comp_enepot.eq.1) enepot = 0.0_cd
347 if (comp_energy.eq.1) energy = 0.0_cd
348 if (comp_force.eq.1) force = 0.0_cd
358 if (particlecontributing(i) == 1)
then 361 call kim_model_compute_arguments_get_neighbor_list( &
362 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
365 kim_log_message =
"kim_api_get_neigh" 379 rij(:) = coor(:,j) - coor(:,i)
383 rsqij = dot_product(rij,rij)
384 if ( rsqij .lt. buf%cutsq(1) )
then 387 if (comp_process_d2edr2.eq.1)
then 395 d2eidr = 0.5_cd*d2phi
396 elseif (comp_force.eq.1.or.comp_process_dedr.eq.1)
then 415 if (comp_enepot.eq.1)
then 416 enepot(i) = enepot(i) + 0.5_cd*phi
418 if (comp_energy.eq.1)
then 419 energy = energy + 0.5_cd*phi
424 if (comp_process_dedr.eq.1)
then 425 call kim_model_compute_arguments_process_dedr_term( &
426 model_compute_arguments_handle, deidr, r, c_loc(rij(1)), i, j, ierr)
430 if (comp_process_d2edr2.eq.1)
then 440 call kim_model_compute_arguments_process_d2edr2_term( &
441 model_compute_arguments_handle, d2eidr, &
443 c_loc(rij_pairs(1,1)), &
445 c_loc(j_pairs(1)), ierr)
450 if (comp_force.eq.1)
then 451 force(:,i) = force(:,i) + deidr*rij/r
452 force(:,j) = force(:,j) - deidr*rij/r
466 if (comp_process_d2edr2.eq.1)
then 467 deallocate( r_pairs )
468 deallocate( rij_pairs )
469 deallocate( i_pairs )
470 deallocate( j_pairs )
485 subroutine refresh(model_refresh_handle, ierr) bind(c)
489 type(kim_model_refresh_handle_type),
intent(inout) :: model_refresh_handle
490 integer(c_int),
intent(out) :: ierr
493 real(c_double) energy_at_cutoff
494 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
497 call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
498 call c_f_pointer(pbuf, buf)
500 call kim_model_refresh_set_influence_distance_pointer(model_refresh_handle, &
501 buf%influence_distance(1))
502 call kim_model_refresh_set_neighbor_list_pointers(model_refresh_handle, &
503 1, buf%influence_distance, buf%padding_neighbor_hints, buf%half_list_hints)
507 buf%influence_distance(1) = buf%Pcutoff(1)
508 buf%cutsq(1) = (buf%Pcutoff(1))**2
514 buf%Pcutoff(1),energy_at_cutoff)
515 buf%shift(1) = -energy_at_cutoff
527 subroutine destroy(model_destroy_handle, ierr) bind(c)
531 type(kim_model_destroy_handle_type),
intent(inout) :: model_destroy_handle
532 integer(c_int),
intent(out) :: ierr
535 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
538 call kim_model_destroy_get_model_buffer_pointer(model_destroy_handle, pbuf)
539 call c_f_pointer(pbuf, buf)
553 #include "kim_model_compute_arguments_create_log_macros.fd" 555 model_compute_arguments_create_handle, ierr) bind(c)
557 log_entry=>kim_model_compute_arguments_create_log_entry
561 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
562 type(kim_model_compute_arguments_create_handle_type),
intent(inout) :: &
563 model_compute_arguments_create_handle
564 integer(c_int),
intent(out) :: ierr
572 call kim_model_compute_arguments_create_set_argument_support_status( &
573 model_compute_arguments_create_handle, &
574 kim_compute_argument_name_partial_energy, &
575 kim_support_status_optional, ierr)
576 call kim_model_compute_arguments_create_set_argument_support_status( &
577 model_compute_arguments_create_handle, &
578 kim_compute_argument_name_partial_forces, &
579 kim_support_status_optional, ierr2)
581 call kim_model_compute_arguments_create_set_argument_support_status( &
582 model_compute_arguments_create_handle, &
583 kim_compute_argument_name_partial_particle_energy, &
584 kim_support_status_optional, ierr2)
587 kim_log_message =
"Unable to register arguments support_statuss" 594 model_compute_arguments_create_handle, &
595 kim_compute_callback_name_process_dedr_term, &
596 kim_support_status_optional, ierr)
598 model_compute_arguments_create_handle, &
599 kim_compute_callback_name_process_d2edr2_term, &
600 kim_support_status_optional, ierr2)
603 kim_log_message =
"Unable to register callbacks support_statuss" 619 #include "kim_model_compute_arguments_destroy_log_macros.fd" 621 model_compute_arguments_destroy_handle, ierr) bind(c)
623 log_entry=>kim_model_compute_arguments_destroy_log_entry
627 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
628 type(kim_model_compute_arguments_destroy_handle_type),
intent(inout) :: &
629 model_compute_arguments_destroy_handle
630 integer(c_int),
intent(out) :: ierr
646 #include "kim_model_driver_create_log_macros.fd" 648 requested_length_unit, requested_energy_unit, requested_charge_unit, &
649 requested_temperature_unit, requested_time_unit, ierr) bind(c)
650 use,
intrinsic :: iso_c_binding
654 integer(c_int),
parameter :: cd = c_double
657 type(kim_model_driver_create_handle_type),
intent(inout) &
658 :: model_driver_create_handle
659 type(kim_length_unit_type),
intent(in),
value :: requested_length_unit
660 type(kim_energy_unit_type),
intent(in),
value :: requested_energy_unit
661 type(kim_charge_unit_type),
intent(in),
value :: requested_charge_unit
662 type(kim_temperature_unit_type),
intent(in),
value :: requested_temperature_unit
663 type(kim_time_unit_type),
intent(in),
value :: requested_time_unit
664 integer(c_int),
intent(out) :: ierr
667 integer(c_int) :: number_of_parameter_files
668 character(len=1024, kind=c_char) :: parameter_file_name
669 integer(c_int) :: ierr2
670 integer(c_int),
parameter :: one=1
671 type(BUFFER_TYPE),
pointer :: buf;
672 type(kim_species_name_type) species_name
674 real(c_double) factor
675 character(len=100, kind=c_char) in_species
676 real(c_double) in_cutoff
677 real(c_double) in_epsilon
678 real(c_double) in_sigma
679 real(c_double) energy_at_cutoff
681 kim_log_file = __file__
684 call kim_model_driver_create_set_units( &
685 model_driver_create_handle, &
686 requested_length_unit, &
687 requested_energy_unit, &
688 kim_charge_unit_unused, &
689 kim_temperature_unit_unused, &
690 kim_time_unit_unused, ierr)
692 kim_log_message =
"Unable to set units" 698 call kim_model_driver_create_set_model_numbering( &
699 model_driver_create_handle, kim_numbering_one_based, ierr)
701 kim_log_message =
"Unable to set numbering" 708 call kim_model_driver_create_set_compute_pointer( &
709 model_driver_create_handle, kim_language_name_fortran, &
711 call kim_model_driver_create_set_compute_arguments_create_pointer( &
712 model_driver_create_handle, kim_language_name_fortran, &
715 call kim_model_driver_create_set_compute_arguments_destroy_pointer( &
716 model_driver_create_handle, kim_language_name_fortran, &
719 call kim_model_driver_create_set_refresh_pointer( &
720 model_driver_create_handle, kim_language_name_fortran, &
723 call kim_model_driver_create_set_destroy_pointer( &
724 model_driver_create_handle, kim_language_name_fortran, &
728 kim_log_message =
"Unable to store callback pointers" 735 call kim_model_driver_create_get_number_of_parameter_files( &
736 model_driver_create_handle, number_of_parameter_files)
737 if (number_of_parameter_files .ne. 1)
then 738 kim_log_message =
"Wrong number of parameter files" 746 call kim_model_driver_create_get_parameter_file_name( &
747 model_driver_create_handle, 1, parameter_file_name, ierr)
749 kim_log_message =
"Unable to get parameter file name" 754 open(10,file=parameter_file_name,status=
"old")
755 read(10,*,iostat=ierr,err=100) in_species
756 read(10,*,iostat=ierr,err=100) in_cutoff
757 read(10,*,iostat=ierr,err=100) in_epsilon
758 read(10,*,iostat=ierr,err=100) in_sigma
765 kim_log_message =
"Unable to read LJ parameters" 773 call kim_species_name_from_string(in_species, species_name)
775 kim_log_message =
"Unable to set species_name" 780 call kim_model_driver_create_set_species_code( &
781 model_driver_create_handle, species_name, speccode, ierr)
783 kim_log_message =
"Unable to set species code" 789 call kim_model_driver_create_convert_unit( &
790 model_driver_create_handle, &
792 kim_energy_unit_ev, &
794 kim_temperature_unit_k, &
796 requested_length_unit, &
797 requested_energy_unit, &
798 requested_charge_unit, &
799 requested_temperature_unit, &
800 requested_time_unit, &
801 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
803 kim_log_message =
"kim_api_convert_to_act_unit" 807 in_cutoff = in_cutoff * factor
809 call kim_model_driver_create_convert_unit( &
810 model_driver_create_handle, &
812 kim_energy_unit_ev, &
814 kim_temperature_unit_k, &
816 requested_length_unit, &
817 requested_energy_unit, &
818 requested_charge_unit, &
819 requested_temperature_unit, &
820 requested_time_unit, &
821 0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
823 kim_log_message =
"kim_api_convert_to_act_unit" 827 in_epsilon = in_epsilon * factor
829 call kim_model_driver_create_convert_unit( &
830 model_driver_create_handle, &
832 kim_energy_unit_ev, &
834 kim_temperature_unit_k, &
836 requested_length_unit, &
837 requested_energy_unit, &
838 requested_charge_unit, &
839 requested_temperature_unit, &
840 requested_time_unit, &
841 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
843 kim_log_message =
"kim_api_convert_to_act_unit" 847 in_sigma = in_sigma * factor
853 buf%influence_distance(1) = in_cutoff
854 buf%Pcutoff(1) = in_cutoff
855 buf%cutsq(1) = in_cutoff**2
856 buf%padding_neighbor_hints = 1
857 buf%half_list_hints = 0
858 buf%epsilon(1) = in_epsilon
859 buf%sigma(1) = in_sigma
864 in_cutoff, energy_at_cutoff)
865 buf%shift(1) = -energy_at_cutoff
868 call kim_model_driver_create_set_influence_distance_pointer( &
869 model_driver_create_handle, buf%influence_distance(1))
870 call kim_model_driver_create_set_neighbor_list_pointers( &
871 model_driver_create_handle, 1, buf%influence_distance, &
872 buf%padding_neighbor_hints, buf%half_list_hints)
877 call kim_model_driver_create_set_model_buffer_pointer( &
878 model_driver_create_handle, c_loc(buf))
881 call kim_model_driver_create_set_parameter_pointer( &
882 model_driver_create_handle, buf%pcutoff,
"cutoff", ierr)
884 kim_log_message =
"set_parameter" 889 call kim_model_driver_create_set_parameter_pointer( &
890 model_driver_create_handle, buf%epsilon,
"epsilon", ierr)
892 kim_log_message =
"set_parameter" 897 call kim_model_driver_create_set_parameter_pointer( &
898 model_driver_create_handle, buf%sigma,
"sigma", ierr)
900 kim_log_message =
"set_parameter" subroutine, public destroy(model_destroy_handle, ierr)
subroutine, public compute_energy_forces(model_compute_handle, model_compute_arguments_handle, ierr)
subroutine, public calc_phi_dphi_d2phi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi, dphi, d2phi)
static void calc_phi_dphi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi, double *dphi)
subroutine, public refresh(model_refresh_handle, ierr)
static int compute_arguments_create(KIM_ModelCompute const *const modelCompute, KIM_ModelComputeArgumentsCreate *const modelComputeArgumentsCreate)
static void calc_phi(double *epsilon, double *C, double *Rzero, double *shift, double *cutoff, double r, double *phi)
subroutine model_driver_create_routine(model_driver_create_handle, requested_length_unit, requested_energy_unit, requested_charge_unit, requested_temperature_unit, requested_time_unit, ierr)
static int compute_arguments_destroy(KIM_ModelCompute const *const modelCompute, KIM_ModelComputeArgumentsDestroy *const modelComputeArgumentsDestroy)