48 use,
intrinsic :: iso_c_binding
63 integer(c_int),
parameter :: cd = c_double
64 integer(c_int),
parameter :: DIM=3
72 type, bind(c) :: buffer_type
73 real(c_double) :: influence_distance(1)
74 real(c_double) :: pcutoff(1)
75 real(c_double) :: cutsq(1)
76 real(c_double) :: epsilon(1)
77 real(c_double) :: sigma(1)
78 real(c_double) :: shift(1)
96 real(c_double),
intent(in) :: model_epsilon
97 real(c_double),
intent(in) :: model_sigma
98 real(c_double),
intent(in) :: model_shift
99 real(c_double),
intent(in) :: model_cutoff
100 real(c_double),
intent(in) :: r
101 real(c_double),
intent(out) :: phi
104 real(c_double) rsq,sor,sor6,sor12
111 if (r .gt. model_cutoff)
then 115 phi = 4.0_cd*model_epsilon*(sor12-sor6) + model_shift
128 model_cutoff,r,phi,dphi)
132 real(c_double),
intent(in) :: model_epsilon
133 real(c_double),
intent(in) :: model_sigma
134 real(c_double),
intent(in) :: model_shift
135 real(c_double),
intent(in) :: model_cutoff
136 real(c_double),
intent(in) :: r
137 real(c_double),
intent(out) :: phi,dphi
140 real(c_double) rsq,sor,sor6,sor12
147 if (r .gt. model_cutoff)
then 152 phi = 4.0_cd*model_epsilon*(sor12-sor6) + model_shift
153 dphi = 24.0_cd*model_epsilon*(-2.0_cd*sor12+sor6)/r
166 model_cutoff,r,phi,dphi,d2phi)
170 real(c_double),
intent(in) :: model_epsilon
171 real(c_double),
intent(in) :: model_sigma
172 real(c_double),
intent(in) :: model_shift
173 real(c_double),
intent(in) :: model_cutoff
174 real(c_double),
intent(in) :: r
175 real(c_double),
intent(out) :: phi,dphi,d2phi
178 real(c_double) rsq,sor,sor6,sor12
185 if (r .gt. model_cutoff)
then 191 phi = 4.0_cd*model_epsilon*(sor12-sor6) + model_shift
192 dphi = 24.0_cd*model_epsilon*(-2.0_cd*sor12+sor6)/r
193 d2phi = 24.0_cd*model_epsilon*(26.0_cd*sor12-7.0_cd*sor6)/rsq
203 #include "kim_model_compute_log_macros.fd" 212 type(kim_model_compute_handle_type),
intent(in) :: model_compute_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(:)
243 call c_f_pointer(pbuf, buf)
245 model_cutoff = buf%influence_distance(1)
251 call kim_model_compute_is_callback_present(model_compute_handle, &
254 call kim_model_compute_is_callback_present(model_compute_handle, &
266 call kim_model_compute_get_argument_pointer(model_compute_handle, &
271 call kim_model_compute_get_argument_pointer(model_compute_handle, &
273 n, particlespeciescodes, ierr2)
275 call kim_model_compute_get_argument_pointer(model_compute_handle, &
279 call kim_model_compute_get_argument_pointer(model_compute_handle, &
282 call kim_model_compute_get_argument_pointer(model_compute_handle, &
285 call kim_model_compute_get_argument_pointer(model_compute_handle, &
288 call kim_model_compute_get_argument_pointer(model_compute_handle, &
297 if (
associated(energy))
then 302 if (
associated(force))
then 307 if (
associated(enepot))
then 314 if (comp_process_d2edr2.eq.1)
then 315 allocate( r_pairs(2) )
316 allocate( rij_pairs(dim,2) )
317 allocate( i_pairs(2) )
318 allocate( j_pairs(2) )
326 if (particlespeciescodes(i).ne.
speccode)
then 336 if (comp_enepot.eq.1) enepot = 0.0_cd
337 if (comp_energy.eq.1) energy = 0.0_cd
338 if (comp_force.eq.1) force = 0.0_cd
348 if (particlecontributing(i) == 1)
then 352 model_compute_handle, 1, i, numnei, nei1part, ierr)
369 rij(:) = coor(:,j) - coor(:,i)
373 rsqij = dot_product(rij,rij)
374 if ( rsqij .lt. buf%cutsq(1) )
then 377 if (comp_process_d2edr2.eq.1)
then 385 d2eidr = 0.5_cd*d2phi
386 elseif (comp_force.eq.1.or.comp_process_dedr.eq.1)
then 405 if (comp_enepot.eq.1)
then 406 enepot(i) = enepot(i) + 0.5_cd*phi
408 if (comp_energy.eq.1)
then 409 energy = energy + 0.5_cd*phi
414 if (comp_process_dedr.eq.1)
then 415 call kim_model_compute_process_dedr_term( &
416 model_compute_handle, deidr, r, c_loc(rij(1)), i, j, ierr)
420 if (comp_process_d2edr2.eq.1)
then 430 call kim_model_compute_process_d2edr2_term( &
431 model_compute_handle, d2eidr, &
433 c_loc(rij_pairs(1,1)), &
435 c_loc(j_pairs(1)), ierr)
440 if (comp_force.eq.1)
then 441 force(:,i) = force(:,i) + deidr*rij/r
442 force(:,j) = force(:,j) - deidr*rij/r
456 if (comp_process_d2edr2.eq.1)
then 457 deallocate( r_pairs )
458 deallocate( rij_pairs )
459 deallocate( i_pairs )
460 deallocate( j_pairs )
475 subroutine refresh(model_refresh_handle, ierr) bind(c)
480 type(kim_model_refresh_handle_type),
intent(inout) :: model_refresh_handle
481 integer(c_int),
intent(out) :: ierr
484 real(c_double) energy_at_cutoff
485 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
488 call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
489 call c_f_pointer(pbuf, buf)
491 call kim_model_refresh_set_influence_distance_pointer(model_refresh_handle, &
492 buf%influence_distance(1))
493 call kim_model_refresh_set_neighbor_list_cutoffs_pointer(model_refresh_handle, &
494 1, buf%influence_distance(1))
498 buf%influence_distance(1) = buf%Pcutoff(1)
499 buf%cutsq(1) = (buf%Pcutoff(1))**2
505 buf%Pcutoff(1),energy_at_cutoff)
506 buf%shift(1) = -energy_at_cutoff
518 subroutine destroy(model_destroy_handle, ierr) bind(c)
523 type(kim_model_destroy_handle_type),
intent(inout) :: model_destroy_handle
524 integer(c_int),
intent(out) :: ierr
527 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
530 call kim_model_destroy_get_model_buffer_pointer(model_destroy_handle, pbuf)
531 call c_f_pointer(pbuf, buf)
547 #include "kim_model_driver_create_log_macros.fd" 549 requested_length_unit, requested_energy_unit, requested_charge_unit, &
550 requested_temperature_unit, requested_time_unit, ierr) bind(c)
551 use,
intrinsic :: iso_c_binding
563 integer(c_int),
parameter :: cd = c_double
566 type(kim_model_driver_create_handle_type),
intent(inout) &
567 :: model_driver_create_handle
568 type(kim_length_unit_type),
intent(in),
value :: requested_length_unit
569 type(kim_energy_unit_type),
intent(in),
value :: requested_energy_unit
570 type(kim_charge_unit_type),
intent(in),
value :: requested_charge_unit
571 type(kim_temperature_unit_type),
intent(in),
value :: requested_temperature_unit
572 type(kim_time_unit_type),
intent(in),
value :: requested_time_unit
573 integer(c_int),
intent(out) :: ierr
576 integer(c_int) :: number_of_parameter_files
577 character(len=1024) :: parameter_file_name
578 integer(c_int) :: ierr2
579 integer(c_int),
parameter :: one=1
580 type(BUFFER_TYPE),
pointer :: buf;
581 type(kim_species_name_type) species_name
583 real(c_double) factor
584 character(len=100) in_species
585 real(c_double) in_cutoff
586 real(c_double) in_epsilon
587 real(c_double) in_sigma
588 real(c_double) energy_at_cutoff
593 call kim_model_driver_create_set_model_numbering( &
637 call kim_model_driver_create_set_refresh_pointer( &
653 call kim_model_driver_create_get_number_of_parameter_files( &
654 model_driver_create_handle, number_of_parameter_files)
655 if (number_of_parameter_files .ne. 1)
then 664 call kim_model_driver_create_get_parameter_file_name( &
665 model_driver_create_handle, 1, parameter_file_name, ierr)
672 open(10,file=parameter_file_name,status=
"old")
673 read(10,*,iostat=ierr,err=100) in_species
674 read(10,*,iostat=ierr,err=100) in_cutoff
675 read(10,*,iostat=ierr,err=100) in_epsilon
676 read(10,*,iostat=ierr,err=100) in_sigma
691 call kim_species_name_from_string(in_species, species_name)
699 model_driver_create_handle, species_name,
speccode, ierr)
708 model_driver_create_handle, &
710 kim_energy_unit_ev, &
712 kim_temperature_unit_k, &
714 requested_length_unit, &
715 requested_energy_unit, &
716 requested_charge_unit, &
717 requested_temperature_unit, &
718 requested_time_unit, &
719 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
725 in_cutoff = in_cutoff * factor
728 model_driver_create_handle, &
730 kim_energy_unit_ev, &
732 kim_temperature_unit_k, &
734 requested_length_unit, &
735 requested_energy_unit, &
736 requested_charge_unit, &
737 requested_temperature_unit, &
738 requested_time_unit, &
739 0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
745 in_epsilon = in_epsilon * factor
748 model_driver_create_handle, &
750 kim_energy_unit_ev, &
752 kim_temperature_unit_k, &
754 requested_length_unit, &
755 requested_energy_unit, &
756 requested_charge_unit, &
757 requested_temperature_unit, &
758 requested_time_unit, &
759 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
765 in_sigma = in_sigma * factor
771 buf%influence_distance(1) = in_cutoff
772 buf%Pcutoff(1) = in_cutoff
773 buf%cutsq(1) = in_cutoff**2
774 buf%epsilon(1) = in_epsilon
775 buf%sigma(1) = in_sigma
780 in_cutoff, energy_at_cutoff)
781 buf%shift(1) = -energy_at_cutoff
785 model_driver_create_handle, buf%influence_distance(1))
786 call kim_model_driver_create_set_neighbor_list_cutoffs_pointer( &
787 model_driver_create_handle, 1, buf%influence_distance(1))
793 model_driver_create_handle, c_loc(buf))
796 call kim_model_driver_create_set_parameter_pointer( &
797 model_driver_create_handle, buf%pcutoff,
"cutoff", ierr)
804 call kim_model_driver_create_set_parameter_pointer( &
805 model_driver_create_handle, buf%epsilon,
"epsilon", ierr)
812 call kim_model_driver_create_set_parameter_pointer( &
813 model_driver_create_handle, buf%sigma,
"sigma", ierr)
type(kim_numbering_type), public, protected kim_numbering_one_based
type(kim_argument_name_type), public, protected kim_argument_name_number_of_particles
character(len=4096), public kim_log_file
character(len=65536), public kim_log_message
type(kim_language_name_type), public, protected kim_language_name_fortran
type(kim_support_status_type), public, protected kim_support_status_optional
subroutine, public destroy(model_destroy_handle, ierr)
type(kim_callback_name_type), public, protected kim_callback_name_process_dedr_term
subroutine, public calc_phi_dphi_d2phi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi, dphi, d2phi)
type(kim_argument_name_type), public, protected kim_argument_name_partial_forces
type(kim_argument_name_type), public, protected kim_argument_name_particle_contributing
subroutine, public refresh(model_refresh_handle, ierr)
type(kim_argument_name_type), public, protected kim_argument_name_particle_species_codes
type(kim_argument_name_type), public, protected kim_argument_name_partial_particle_energy
type(kim_argument_name_type), public, protected kim_argument_name_partial_energy
integer(c_int), parameter, public speccode
subroutine, public compute_energy_forces(model_compute_handle, ierr)
type(kim_callback_name_type), public, protected kim_callback_name_process_d2edr2_term
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)
type(kim_argument_name_type), public, protected kim_argument_name_coordinates
subroutine, public calc_phi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi)
subroutine, public calc_phi_dphi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi, dphi)