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 real(c_double) :: epsilon(1)
78 real(c_double) :: sigma(1)
79 real(c_double) :: shift(1)
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
112 if (r .gt. 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
148 if (r .gt. 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
186 if (r .gt. 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
204 #include "kim_model_compute_log_macros.fd" 206 model_compute_arguments_handle, ierr) bind(c)
210 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
211 type(kim_model_compute_arguments_handle_type),
intent(in) :: &
212 model_compute_arguments_handle
213 integer(c_int),
intent(out) :: ierr
216 real(c_double) :: r,rsqij,phi,dphi,d2phi,deidr,d2eidr
217 integer(c_int) :: i,j,jj,numnei
218 integer(c_int) :: ierr2
219 integer(c_int) :: comp_force,comp_energy,comp_enepot,comp_process_dedr, &
221 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
223 real(c_double),
pointer :: rij(:)
224 real(c_double),
pointer :: rij_pairs(:,:)
225 real(c_double),
pointer :: r_pairs(:)
226 integer(c_int),
pointer :: i_pairs(:), j_pairs(:)
229 real(c_double) :: model_cutoff
230 integer(c_int),
pointer :: n
231 real(c_double),
pointer :: energy
232 real(c_double),
pointer :: coor(:,:)
233 real(c_double),
pointer :: force(:,:)
234 real(c_double),
pointer :: enepot(:)
235 integer(c_int),
pointer :: nei1part(:)
236 integer(c_int),
pointer :: particlespeciescodes(:)
237 integer(c_int),
pointer :: particlecontributing(:)
239 kim_log_file = __file__
242 call kim_model_compute_get_model_buffer_pointer(model_compute_handle, pbuf)
243 call c_f_pointer(pbuf, buf)
245 model_cutoff = buf%influence_distance(1)
251 call kim_model_compute_arguments_is_callback_present( &
252 model_compute_arguments_handle, &
253 kim_compute_callback_name_process_dedr_term, comp_process_dedr, ierr2)
255 call kim_model_compute_arguments_is_callback_present( &
256 model_compute_arguments_handle, &
257 kim_compute_callback_name_process_d2edr2_term, comp_process_d2edr2, ierr2)
260 kim_log_message =
"get_compute" 268 call kim_model_compute_arguments_get_argument_pointer( &
269 model_compute_arguments_handle, &
270 kim_compute_argument_name_number_of_particles, n, ierr2)
273 call kim_model_compute_arguments_get_argument_pointer( &
274 model_compute_arguments_handle, &
275 kim_compute_argument_name_particle_species_codes, &
276 n, particlespeciescodes, ierr2)
278 call kim_model_compute_arguments_get_argument_pointer( &
279 model_compute_arguments_handle, &
280 kim_compute_argument_name_particle_contributing, n, particlecontributing, &
283 call kim_model_compute_arguments_get_argument_pointer( &
284 model_compute_arguments_handle, &
285 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
287 call kim_model_compute_arguments_get_argument_pointer( &
288 model_compute_arguments_handle, &
289 kim_compute_argument_name_partial_energy, energy, ierr2)
291 call kim_model_compute_arguments_get_argument_pointer( &
292 model_compute_arguments_handle, &
293 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
295 call kim_model_compute_arguments_get_argument_pointer( &
296 model_compute_arguments_handle, &
297 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
300 kim_log_message =
"get_argument_pointer" 305 if (
associated(energy))
then 310 if (
associated(force))
then 315 if (
associated(enepot))
then 322 if (comp_process_d2edr2.eq.1)
then 323 allocate( r_pairs(2) )
324 allocate( rij_pairs(dim,2) )
325 allocate( i_pairs(2) )
326 allocate( j_pairs(2) )
334 if (particlespeciescodes(i).ne.speccode)
then 335 kim_log_message =
"Unexpected species code detected" 344 if (comp_enepot.eq.1) enepot = 0.0_cd
345 if (comp_energy.eq.1) energy = 0.0_cd
346 if (comp_force.eq.1) force = 0.0_cd
356 if (particlecontributing(i) == 1)
then 359 call kim_model_compute_arguments_get_neighbor_list( &
360 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
363 kim_log_message =
"kim_api_get_neigh" 377 rij(:) = coor(:,j) - coor(:,i)
381 rsqij = dot_product(rij,rij)
382 if ( rsqij .lt. buf%cutsq(1) )
then 385 if (comp_process_d2edr2.eq.1)
then 393 d2eidr = 0.5_cd*d2phi
394 elseif (comp_force.eq.1.or.comp_process_dedr.eq.1)
then 413 if (comp_enepot.eq.1)
then 414 enepot(i) = enepot(i) + 0.5_cd*phi
416 if (comp_energy.eq.1)
then 417 energy = energy + 0.5_cd*phi
422 if (comp_process_dedr.eq.1)
then 423 call kim_model_compute_arguments_process_dedr_term( &
424 model_compute_arguments_handle, deidr, r, c_loc(rij(1)), i, j, ierr)
428 if (comp_process_d2edr2.eq.1)
then 438 call kim_model_compute_arguments_process_d2edr2_term( &
439 model_compute_arguments_handle, d2eidr, &
441 c_loc(rij_pairs(1,1)), &
443 c_loc(j_pairs(1)), ierr)
448 if (comp_force.eq.1)
then 449 force(:,i) = force(:,i) + deidr*rij/r
450 force(:,j) = force(:,j) - deidr*rij/r
464 if (comp_process_d2edr2.eq.1)
then 465 deallocate( r_pairs )
466 deallocate( rij_pairs )
467 deallocate( i_pairs )
468 deallocate( j_pairs )
483 subroutine refresh(model_refresh_handle, ierr) bind(c)
487 type(kim_model_refresh_handle_type),
intent(inout) :: model_refresh_handle
488 integer(c_int),
intent(out) :: ierr
491 real(c_double) energy_at_cutoff
492 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
495 call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
496 call c_f_pointer(pbuf, buf)
498 call kim_model_refresh_set_influence_distance_pointer(model_refresh_handle, &
499 buf%influence_distance(1))
500 call kim_model_refresh_set_neighbor_list_cutoffs_pointer(model_refresh_handle, &
501 1, buf%influence_distance(1))
505 buf%influence_distance(1) = buf%Pcutoff(1)
506 buf%cutsq(1) = (buf%Pcutoff(1))**2
512 buf%Pcutoff(1),energy_at_cutoff)
513 buf%shift(1) = -energy_at_cutoff
525 subroutine destroy(model_destroy_handle, ierr) bind(c)
529 type(kim_model_destroy_handle_type),
intent(inout) :: model_destroy_handle
530 integer(c_int),
intent(out) :: ierr
533 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
536 call kim_model_destroy_get_model_buffer_pointer(model_destroy_handle, pbuf)
537 call c_f_pointer(pbuf, buf)
551 #include "kim_model_compute_arguments_create_log_macros.fd" 553 model_compute_arguments_create_handle, ierr) bind(c)
555 log_entry=>kim_model_compute_arguments_create_log_entry
559 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
560 type(kim_model_compute_arguments_create_handle_type),
intent(inout) :: &
561 model_compute_arguments_create_handle
562 integer(c_int),
intent(out) :: ierr
570 call kim_model_compute_arguments_create_set_argument_support_status( &
571 model_compute_arguments_create_handle, &
572 kim_compute_argument_name_partial_energy, &
573 kim_support_status_optional, ierr)
574 call kim_model_compute_arguments_create_set_argument_support_status( &
575 model_compute_arguments_create_handle, &
576 kim_compute_argument_name_partial_forces, &
577 kim_support_status_optional, ierr2)
579 call kim_model_compute_arguments_create_set_argument_support_status( &
580 model_compute_arguments_create_handle, &
581 kim_compute_argument_name_partial_particle_energy, &
582 kim_support_status_optional, ierr2)
585 kim_log_message =
"Unable to register arguments support_statuss" 592 model_compute_arguments_create_handle, &
593 kim_compute_callback_name_process_dedr_term, &
594 kim_support_status_optional, ierr)
596 model_compute_arguments_create_handle, &
597 kim_compute_callback_name_process_d2edr2_term, &
598 kim_support_status_optional, ierr2)
601 kim_log_message =
"Unable to register callbacks support_statuss" 617 #include "kim_model_compute_arguments_destroy_log_macros.fd" 619 model_compute_arguments_destroy_handle, ierr) bind(c)
621 log_entry=>kim_model_compute_arguments_destroy_log_entry
625 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
626 type(kim_model_compute_arguments_destroy_handle_type),
intent(inout) :: &
627 model_compute_arguments_destroy_handle
628 integer(c_int),
intent(out) :: ierr
644 #include "kim_model_driver_create_log_macros.fd" 646 requested_length_unit, requested_energy_unit, requested_charge_unit, &
647 requested_temperature_unit, requested_time_unit, ierr) bind(c)
648 use,
intrinsic :: iso_c_binding
652 integer(c_int),
parameter :: cd = c_double
655 type(kim_model_driver_create_handle_type),
intent(inout) &
656 :: model_driver_create_handle
657 type(kim_length_unit_type),
intent(in),
value :: requested_length_unit
658 type(kim_energy_unit_type),
intent(in),
value :: requested_energy_unit
659 type(kim_charge_unit_type),
intent(in),
value :: requested_charge_unit
660 type(kim_temperature_unit_type),
intent(in),
value :: requested_temperature_unit
661 type(kim_time_unit_type),
intent(in),
value :: requested_time_unit
662 integer(c_int),
intent(out) :: ierr
665 integer(c_int) :: number_of_parameter_files
666 character(len=1024, kind=c_char) :: parameter_file_name
667 integer(c_int) :: ierr2
668 integer(c_int),
parameter :: one=1
669 type(BUFFER_TYPE),
pointer :: buf;
670 type(kim_species_name_type) species_name
672 real(c_double) factor
673 character(len=100, kind=c_char) in_species
674 real(c_double) in_cutoff
675 real(c_double) in_epsilon
676 real(c_double) in_sigma
677 real(c_double) energy_at_cutoff
679 kim_log_file = __file__
682 call kim_model_driver_create_set_units( &
683 model_driver_create_handle, &
684 requested_length_unit, &
685 requested_energy_unit, &
686 kim_charge_unit_unused, &
687 kim_temperature_unit_unused, &
688 kim_time_unit_unused, ierr)
690 kim_log_message =
"Unable to set units" 696 call kim_model_driver_create_set_model_numbering( &
697 model_driver_create_handle, kim_numbering_one_based, ierr)
699 kim_log_message =
"Unable to set numbering" 706 call kim_model_driver_create_set_compute_pointer( &
707 model_driver_create_handle, kim_language_name_fortran, &
709 call kim_model_driver_create_set_compute_arguments_create_pointer( &
710 model_driver_create_handle, kim_language_name_fortran, &
713 call kim_model_driver_create_set_compute_arguments_destroy_pointer( &
714 model_driver_create_handle, kim_language_name_fortran, &
717 call kim_model_driver_create_set_refresh_pointer( &
718 model_driver_create_handle, kim_language_name_fortran, &
721 call kim_model_driver_create_set_destroy_pointer( &
722 model_driver_create_handle, kim_language_name_fortran, &
726 kim_log_message =
"Unable to store callback pointers" 733 call kim_model_driver_create_get_number_of_parameter_files( &
734 model_driver_create_handle, number_of_parameter_files)
735 if (number_of_parameter_files .ne. 1)
then 736 kim_log_message =
"Wrong number of parameter files" 744 call kim_model_driver_create_get_parameter_file_name( &
745 model_driver_create_handle, 1, parameter_file_name, ierr)
747 kim_log_message =
"Unable to get parameter file name" 752 open(10,file=parameter_file_name,status=
"old")
753 read(10,*,iostat=ierr,err=100) in_species
754 read(10,*,iostat=ierr,err=100) in_cutoff
755 read(10,*,iostat=ierr,err=100) in_epsilon
756 read(10,*,iostat=ierr,err=100) in_sigma
763 kim_log_message =
"Unable to read LJ parameters" 771 call kim_species_name_from_string(in_species, species_name)
773 kim_log_message =
"Unable to set species_name" 778 call kim_model_driver_create_set_species_code( &
779 model_driver_create_handle, species_name, speccode, ierr)
781 kim_log_message =
"Unable to set species code" 787 call kim_model_driver_create_convert_unit( &
788 model_driver_create_handle, &
790 kim_energy_unit_ev, &
792 kim_temperature_unit_k, &
794 requested_length_unit, &
795 requested_energy_unit, &
796 requested_charge_unit, &
797 requested_temperature_unit, &
798 requested_time_unit, &
799 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
801 kim_log_message =
"kim_api_convert_to_act_unit" 805 in_cutoff = in_cutoff * factor
807 call kim_model_driver_create_convert_unit( &
808 model_driver_create_handle, &
810 kim_energy_unit_ev, &
812 kim_temperature_unit_k, &
814 requested_length_unit, &
815 requested_energy_unit, &
816 requested_charge_unit, &
817 requested_temperature_unit, &
818 requested_time_unit, &
819 0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
821 kim_log_message =
"kim_api_convert_to_act_unit" 825 in_epsilon = in_epsilon * factor
827 call kim_model_driver_create_convert_unit( &
828 model_driver_create_handle, &
830 kim_energy_unit_ev, &
832 kim_temperature_unit_k, &
834 requested_length_unit, &
835 requested_energy_unit, &
836 requested_charge_unit, &
837 requested_temperature_unit, &
838 requested_time_unit, &
839 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
841 kim_log_message =
"kim_api_convert_to_act_unit" 845 in_sigma = in_sigma * factor
851 buf%influence_distance(1) = in_cutoff
852 buf%Pcutoff(1) = in_cutoff
853 buf%cutsq(1) = in_cutoff**2
854 buf%epsilon(1) = in_epsilon
855 buf%sigma(1) = in_sigma
860 in_cutoff, energy_at_cutoff)
861 buf%shift(1) = -energy_at_cutoff
864 call kim_model_driver_create_set_influence_distance_pointer( &
865 model_driver_create_handle, buf%influence_distance(1))
866 call kim_model_driver_create_set_neighbor_list_cutoffs_pointer( &
867 model_driver_create_handle, 1, buf%influence_distance(1))
872 call kim_model_driver_create_set_model_buffer_pointer( &
873 model_driver_create_handle, c_loc(buf))
876 call kim_model_driver_create_set_parameter_pointer( &
877 model_driver_create_handle, buf%pcutoff,
"cutoff", ierr)
879 kim_log_message =
"set_parameter" 884 call kim_model_driver_create_set_parameter_pointer( &
885 model_driver_create_handle, buf%epsilon,
"epsilon", ierr)
887 kim_log_message =
"set_parameter" 892 call kim_model_driver_create_set_parameter_pointer( &
893 model_driver_create_handle, buf%sigma,
"sigma", ierr)
895 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)