47 use,
intrinsic :: iso_c_binding
62 integer(c_int),
parameter :: cd = c_double
63 integer(c_int),
parameter :: dim = 3
67 real(c_double),
parameter :: model_cutsq =
model_cutoff**2
84 real(c_double),
parameter :: lj_epsilon = 0.0104_cd
85 real(c_double),
parameter :: lj_sigma = 3.40_cd
86 real(c_double),
parameter :: lj_cutnorm =
model_cutoff/lj_sigma
87 real(c_double),
parameter :: lj_a = 12.0_cd*lj_epsilon*(-26.0_cd &
88 + 7.0_cd*lj_cutnorm**6)/(lj_cutnorm**14 &
90 real(c_double),
parameter :: lj_b = 96.0_cd*lj_epsilon*(7.0_cd &
91 - 2.0_cd*lj_cutnorm**6)/(lj_cutnorm**13*lj_sigma)
92 real(c_double),
parameter :: lj_c = 28.0_cd*lj_epsilon*(-13.0_cd &
93 + 4.0_cd*lj_cutnorm**6)/(lj_cutnorm**12)
95 type, bind(c) :: buffer_type
96 real(c_double) :: influence_distance
97 real(c_double) :: cutoff(1)
99 model_will_not_request_neighbors_of_noncontributing_particles(1)
109 recursive subroutine calc_phi(r,phi)
113 real(c_double),
intent(in) :: r
114 real(c_double),
intent(out) :: phi
117 real(c_double) rsq,sor,sor6,sor12
128 phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
142 real(c_double),
intent(in) :: r
143 real(c_double),
intent(out) :: phi,dphi
146 real(c_double) rsq,sor,sor6,sor12
158 phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
159 dphi = 24.0_cd*lj_epsilon*(-2.0_cd*sor12+sor6)/r + 2.0_cd*lj_a*r + lj_b
170 model_compute_arguments_handle, ierr) bind(c)
174 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
175 type(kim_model_compute_arguments_handle_type),
intent(in) :: &
176 model_compute_arguments_handle
177 integer(c_int),
intent(out) :: ierr
180 real(c_double) :: Rij(dim)
181 real(c_double) :: r,Rsqij,phi,dphi,dEidr = 0.0_cd
182 integer(c_int) :: i,j,jj,numnei,comp_force,comp_enepot,comp_virial, comp_energy
183 integer(c_int) :: ierr2
186 integer(c_int),
pointer :: N
187 real(c_double),
pointer :: energy
188 real(c_double),
pointer :: coor(:,:)
189 real(c_double),
pointer :: force(:,:)
190 real(c_double),
pointer :: enepot(:)
191 integer(c_int),
pointer :: nei1part(:)
192 integer(c_int),
pointer :: particleSpeciesCodes(:)
193 integer(c_int),
pointer :: particleContributing(:)
194 real(c_double),
pointer :: virial(:)
199 call kim_get_argument_pointer( &
200 model_compute_arguments_handle, &
201 kim_compute_argument_name_number_of_particles, n, ierr2)
203 call kim_get_argument_pointer( &
204 model_compute_arguments_handle, &
205 kim_compute_argument_name_particle_species_codes, n, particlespeciescodes, &
208 call kim_get_argument_pointer( &
209 model_compute_arguments_handle, &
210 kim_compute_argument_name_particle_contributing, n, particlecontributing, &
213 call kim_get_argument_pointer( &
214 model_compute_arguments_handle, &
215 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
217 call kim_get_argument_pointer( &
218 model_compute_arguments_handle, &
219 kim_compute_argument_name_partial_energy, energy, ierr2)
221 call kim_get_argument_pointer( &
222 model_compute_arguments_handle, &
223 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
225 call kim_get_argument_pointer( &
226 model_compute_arguments_handle, &
227 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
229 call kim_get_argument_pointer( &
230 model_compute_arguments_handle, &
231 kim_compute_argument_name_partial_virial, 6, virial, ierr2)
234 call kim_log_entry(model_compute_arguments_handle, &
235 kim_log_verbosity_error,
"get data")
243 if (
associated(energy))
then 248 if (
associated(force))
then 253 if (
associated(enepot))
then 258 if (
associated(virial))
then 268 if (particlespeciescodes(i).ne.
speccode)
then 269 call kim_log_entry(model_compute_handle, &
270 kim_log_verbosity_error,
"Unexpected species code detected")
279 if (comp_enepot.eq.1) enepot = 0.0_cd
280 if (comp_energy.eq.1) energy = 0.0_cd
281 if (comp_force.eq.1) force = 0.0_cd
282 if (comp_virial.eq.1) virial = 0.0_cd
291 if (particlecontributing(i) == 1)
then 293 call kim_get_neighbor_list( &
294 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
298 model_compute_arguments_handle, kim_log_verbosity_error, &
299 "GetNeighborList failed")
312 rij(:) = coor(:,j) - coor(:,i)
316 rsqij = dot_product(rij,rij)
317 if ( rsqij .lt. model_cutsq )
then 320 if (comp_force.eq.1.or.comp_virial.eq.1)
then 330 if (comp_enepot.eq.1)
then 331 enepot(i) = enepot(i) + 0.5_cd*phi
333 if (comp_energy.eq.1)
then 334 energy = energy + 0.5_cd*phi
339 if (comp_virial.eq.1)
then 340 virial(1) = virial(1) + rij(1)*rij(1)*deidr/r
341 virial(2) = virial(2) + rij(2)*rij(2)*deidr/r
342 virial(3) = virial(3) + rij(3)*rij(3)*deidr/r
343 virial(4) = virial(4) + rij(2)*rij(3)*deidr/r
344 virial(5) = virial(5) + rij(1)*rij(3)*deidr/r
345 virial(6) = virial(6) + rij(1)*rij(2)*deidr/r
350 if (comp_force.eq.1)
then 351 force(:,i) = force(:,i) + deidr*rij/r
352 force(:,j) = force(:,j) - deidr*rij/r
376 use,
intrinsic :: iso_c_binding
380 type(kim_model_destroy_handle_type),
intent(inout) :: model_destroy_handle
381 integer(c_int),
intent(out) :: ierr
383 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
385 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
386 call c_f_pointer(pbuf, buf)
387 call kim_log_entry(model_destroy_handle, &
388 kim_log_verbosity_error,
"deallocating model buffer")
399 model_compute_arguments_create_handle, ierr) bind(c)
400 use,
intrinsic :: iso_c_binding
404 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
405 type(kim_model_compute_arguments_create_handle_type),
intent(inout) :: &
406 model_compute_arguments_create_handle
407 integer(c_int),
intent(out) :: ierr
409 integer(c_int) :: ierr2
412 if (model_compute_handle .eq. kim_model_compute_null_handle)
continue 418 call kim_set_argument_support_status( &
419 model_compute_arguments_create_handle, &
420 kim_compute_argument_name_partial_energy, &
421 kim_support_status_optional, ierr2)
423 call kim_set_argument_support_status( &
424 model_compute_arguments_create_handle, &
425 kim_compute_argument_name_partial_forces, &
426 kim_support_status_optional, ierr2)
428 call kim_set_argument_support_status( &
429 model_compute_arguments_create_handle, &
430 kim_compute_argument_name_partial_particle_energy, &
431 kim_support_status_optional, ierr2)
433 call kim_set_argument_support_status( &
434 model_compute_arguments_create_handle, &
435 kim_compute_argument_name_partial_virial, &
436 kim_support_status_optional, ierr2)
444 call kim_log_entry( &
445 model_compute_arguments_create_handle, kim_log_verbosity_error, &
446 "Unable to successfully create compute_arguments object")
459 model_compute_arguments_destroy_handle, ierr) bind(c)
460 use,
intrinsic :: iso_c_binding
464 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
465 type(kim_model_compute_arguments_destroy_handle_type),
intent(inout) :: &
466 model_compute_arguments_destroy_handle
467 integer(c_int),
intent(out) :: ierr
470 if (model_compute_handle .eq. kim_model_compute_null_handle)
continue 471 if (model_compute_arguments_destroy_handle .eq. &
472 kim_model_compute_arguments_destroy_null_handle)
continue 485 recursive subroutine model_create_routine(model_create_handle, &
486 requested_length_unit, requested_energy_unit, requested_charge_unit, &
487 requested_temperature_unit, requested_time_unit, ierr) bind(c)
488 use,
intrinsic :: iso_c_binding
494 type(kim_model_create_handle_type),
intent(inout) :: model_create_handle
495 type(kim_length_unit_type),
intent(in),
value :: requested_length_unit
496 type(kim_energy_unit_type),
intent(in),
value :: requested_energy_unit
497 type(kim_charge_unit_type),
intent(in),
value :: requested_charge_unit
498 type(kim_temperature_unit_type),
intent(in),
value :: requested_temperature_unit
499 type(kim_time_unit_type),
intent(in),
value :: requested_time_unit
500 integer(c_int),
intent(out) :: ierr
503 integer(c_int) :: ierr2
504 type(buffer_type),
pointer :: buf
510 if (requested_length_unit .eq. kim_length_unit_unused)
continue 511 if (requested_energy_unit .eq. kim_energy_unit_unused)
continue 512 if (requested_charge_unit .eq. kim_charge_unit_unused)
continue 513 if (requested_temperature_unit .eq. kim_temperature_unit_unused)
continue 514 if (requested_time_unit .eq. kim_time_unit_unused)
continue 517 call kim_set_units(model_create_handle, &
519 kim_energy_unit_ev, &
520 kim_charge_unit_unused, &
521 kim_temperature_unit_unused, &
522 kim_time_unit_unused, &
527 call kim_set_species_code(model_create_handle, &
528 kim_species_name_ar,
speccode, ierr2)
532 call kim_set_model_numbering(model_create_handle, &
533 kim_numbering_one_based, ierr2);
537 call kim_set_routine_pointer(model_create_handle, &
538 kim_model_routine_name_compute, kim_language_name_fortran, &
541 call kim_set_routine_pointer( &
542 model_create_handle, kim_model_routine_name_compute_arguments_create, &
545 call kim_set_routine_pointer( &
546 model_create_handle, kim_model_routine_name_compute_arguments_destroy, &
550 call kim_set_routine_pointer(model_create_handle, &
551 kim_model_routine_name_destroy, kim_language_name_fortran, &
559 call kim_set_model_buffer_pointer(model_create_handle, &
565 buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
568 call kim_set_influence_distance_pointer( &
569 model_create_handle, buf%influence_distance)
572 call kim_set_neighbor_list_pointers(model_create_handle, &
574 buf%model_will_not_request_neighbors_of_noncontributing_particles)
579 call kim_log_entry(model_create_handle, &
580 kim_log_verbosity_error,
"Unable to successfully initialize model")
586 end subroutine model_create_routine
recursive subroutine, public compute_energy_forces(model_compute_handle, model_compute_arguments_handle, ierr)
static void calc_phi_dphi(double const *epsilon, double const *C, double const *Rzero, double const *shift, double const cutoff, double const r, double *phi, double *dphi)
recursive subroutine, public model_destroy_func(model_destroy_handle, ierr)
recursive subroutine, public model_compute_arguments_create(model_compute_handle, model_compute_arguments_create_handle, ierr)
static void calc_phi(double const *epsilon, double const *C, double const *Rzero, double const *shift, double const cutoff, double const r, double *phi)
real(c_double), parameter, public model_cutoff
recursive subroutine, public model_compute_arguments_destroy(model_compute_handle, model_compute_arguments_destroy_handle, ierr)
integer(c_int), parameter, public speccode