46 use,
intrinsic :: iso_c_binding
61 integer(c_int),
parameter :: cd = c_double
62 integer(c_int),
parameter :: dim = 3
66 real(c_double),
parameter :: model_cutsq =
model_cutoff**2
83 real(c_double),
parameter :: lj_epsilon = 0.0104_cd
84 real(c_double),
parameter :: lj_sigma = 3.40_cd
85 real(c_double),
parameter :: lj_cutnorm =
model_cutoff / lj_sigma
86 real(c_double),
parameter :: &
87 lj_a = 12.0_cd * lj_epsilon * (-26.0_cd + 7.0_cd * lj_cutnorm**6) &
88 / (lj_cutnorm**14 * lj_sigma**2)
89 real(c_double),
parameter :: &
90 lj_b = 96.0_cd * lj_epsilon * (7.0_cd - 2.0_cd * lj_cutnorm**6) &
91 / (lj_cutnorm**13 * lj_sigma)
92 real(c_double),
parameter :: &
93 lj_c = 28.0_cd * lj_epsilon * (-13.0_cd + 4.0_cd * lj_cutnorm**6) &
96 type, bind(c) :: buffer_type
97 real(c_double) :: influence_distance
98 real(c_double) :: cutoff(1)
100 model_will_not_request_neighbors_of_noncontributing_particles(1)
110 recursive subroutine calc_phi(r, phi)
114 real(c_double),
intent(in) :: r
115 real(c_double),
intent(out) :: phi
118 real(c_double) rsq, sor, sor6, sor12
122 sor6 = sor * sor * sor
129 phi = 4.0_cd * lj_epsilon * (sor12 - sor6) + lj_a * rsq + lj_b * r + lj_c
143 real(c_double),
intent(in) :: r
144 real(c_double),
intent(out) :: phi, dphi
147 real(c_double) rsq, sor, sor6, sor12
151 sor6 = sor * sor * sor
159 phi = 4.0_cd * lj_epsilon * (sor12 - sor6) + lj_a * rsq + lj_b * r + lj_c
160 dphi = 24.0_cd * lj_epsilon * (-2.0_cd * sor12 + sor6) / r &
161 + 2.0_cd * lj_a * r + lj_b
172 model_compute_handle, model_compute_arguments_handle, ierr) bind(c)
176 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
177 type(kim_model_compute_arguments_handle_type),
intent(in) :: &
178 model_compute_arguments_handle
179 integer(c_int),
intent(out) :: ierr
182 real(c_double) :: Rij(dim)
183 real(c_double) :: r, Rsqij, phi, dphi, dEidr = 0.0_cd
184 integer(c_int) :: i, j, jj, numnei, comp_force, comp_enepot, &
185 comp_virial, comp_energy
186 integer(c_int) :: ierr2
189 integer(c_int),
pointer :: N
190 real(c_double),
pointer :: energy
191 real(c_double),
pointer :: coor(:, :)
192 real(c_double),
pointer :: force(:, :)
193 real(c_double),
pointer :: enepot(:)
194 integer(c_int),
pointer :: nei1part(:)
195 integer(c_int),
pointer :: particleSpeciesCodes(:)
196 integer(c_int),
pointer :: particleContributing(:)
197 real(c_double),
pointer :: virial(:)
202 call kim_get_argument_pointer( &
203 model_compute_arguments_handle, &
204 kim_compute_argument_name_number_of_particles, n, ierr2)
206 call kim_get_argument_pointer( &
207 model_compute_arguments_handle, &
208 kim_compute_argument_name_particle_species_codes, n, &
209 particlespeciescodes, ierr2)
211 call kim_get_argument_pointer( &
212 model_compute_arguments_handle, &
213 kim_compute_argument_name_particle_contributing, n, &
214 particlecontributing, ierr2)
216 call kim_get_argument_pointer( &
217 model_compute_arguments_handle, &
218 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
220 call kim_get_argument_pointer( &
221 model_compute_arguments_handle, &
222 kim_compute_argument_name_partial_energy, energy, ierr2)
224 call kim_get_argument_pointer( &
225 model_compute_arguments_handle, &
226 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
228 call kim_get_argument_pointer( &
229 model_compute_arguments_handle, &
230 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
232 call kim_get_argument_pointer( &
233 model_compute_arguments_handle, &
234 kim_compute_argument_name_partial_virial, 6, virial, ierr2)
237 call kim_log_entry(model_compute_arguments_handle, &
238 kim_log_verbosity_error,
"get data")
246 if (
associated(energy))
then 251 if (
associated(force))
then 256 if (
associated(enepot))
then 261 if (
associated(virial))
then 271 if (particlespeciescodes(i) /=
speccode)
then 272 call kim_log_entry( &
273 model_compute_handle, kim_log_verbosity_error, &
274 "Unexpected species code detected")
283 if (comp_enepot == 1) enepot = 0.0_cd
284 if (comp_energy == 1) energy = 0.0_cd
285 if (comp_force == 1) force = 0.0_cd
286 if (comp_virial == 1) virial = 0.0_cd
295 if (particlecontributing(i) == 1)
then 297 call kim_get_neighbor_list( &
298 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
301 call kim_log_entry( &
302 model_compute_arguments_handle, kim_log_verbosity_error, &
303 "GetNeighborList failed")
316 rij(:) = coor(:, j) - coor(:, i)
320 rsqij = dot_product(rij, rij)
321 if (rsqij < model_cutsq)
then 324 if (comp_force == 1 .or. comp_virial == 1)
then 327 deidr = 0.5_cd * dphi
334 if (comp_enepot == 1)
then 335 enepot(i) = enepot(i) + 0.5_cd * phi
337 if (comp_energy == 1)
then 338 energy = energy + 0.5_cd * phi
343 if (comp_virial == 1)
then 344 virial(1) = virial(1) + rij(1) * rij(1) * deidr / r
345 virial(2) = virial(2) + rij(2) * rij(2) * deidr / r
346 virial(3) = virial(3) + rij(3) * rij(3) * deidr / r
347 virial(4) = virial(4) + rij(2) * rij(3) * deidr / r
348 virial(5) = virial(5) + rij(1) * rij(3) * deidr / r
349 virial(6) = virial(6) + rij(1) * rij(2) * deidr / r
354 if (comp_force == 1)
then 355 force(:, i) = force(:, i) + deidr * rij / r
356 force(:, j) = force(:, j) - deidr * rij / r
380 use,
intrinsic :: iso_c_binding
384 type(kim_model_destroy_handle_type),
intent(inout) :: model_destroy_handle
385 integer(c_int),
intent(out) :: ierr
387 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
389 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
390 call c_f_pointer(pbuf, buf)
391 call kim_log_entry(model_destroy_handle, &
392 kim_log_verbosity_error,
"deallocating model buffer")
403 model_compute_handle, model_compute_arguments_create_handle, ierr) bind(c)
404 use,
intrinsic :: iso_c_binding
408 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
409 type(kim_model_compute_arguments_create_handle_type),
intent(inout) :: &
410 model_compute_arguments_create_handle
411 integer(c_int),
intent(out) :: ierr
413 integer(c_int) :: ierr2
416 if (model_compute_handle == kim_model_compute_null_handle)
continue 422 call kim_set_argument_support_status( &
423 model_compute_arguments_create_handle, &
424 kim_compute_argument_name_partial_energy, &
425 kim_support_status_optional, ierr2)
427 call kim_set_argument_support_status( &
428 model_compute_arguments_create_handle, &
429 kim_compute_argument_name_partial_forces, &
430 kim_support_status_optional, ierr2)
432 call kim_set_argument_support_status( &
433 model_compute_arguments_create_handle, &
434 kim_compute_argument_name_partial_particle_energy, &
435 kim_support_status_optional, ierr2)
437 call kim_set_argument_support_status( &
438 model_compute_arguments_create_handle, &
439 kim_compute_argument_name_partial_virial, &
440 kim_support_status_optional, ierr2)
448 call kim_log_entry( &
449 model_compute_arguments_create_handle, kim_log_verbosity_error, &
450 "Unable to successfully create compute_arguments object")
463 model_compute_handle, model_compute_arguments_destroy_handle, ierr) bind(c)
464 use,
intrinsic :: iso_c_binding
468 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
469 type(kim_model_compute_arguments_destroy_handle_type),
intent(inout) :: &
470 model_compute_arguments_destroy_handle
471 integer(c_int),
intent(out) :: ierr
474 if (model_compute_handle == kim_model_compute_null_handle)
continue 475 if (model_compute_arguments_destroy_handle == &
476 kim_model_compute_arguments_destroy_null_handle)
continue 489 recursive subroutine model_create_routine( &
490 model_create_handle, requested_length_unit, requested_energy_unit, &
491 requested_charge_unit, requested_temperature_unit, requested_time_unit, &
493 use,
intrinsic :: iso_c_binding
499 type(kim_model_create_handle_type),
intent(inout) :: model_create_handle
500 type(kim_length_unit_type),
intent(in),
value :: requested_length_unit
501 type(kim_energy_unit_type),
intent(in),
value :: requested_energy_unit
502 type(kim_charge_unit_type),
intent(in),
value :: requested_charge_unit
503 type(kim_temperature_unit_type),
intent(in),
value :: &
504 requested_temperature_unit
505 type(kim_time_unit_type),
intent(in),
value :: requested_time_unit
506 integer(c_int),
intent(out) :: ierr
509 integer(c_int) :: ierr2
510 type(buffer_type),
pointer :: buf
516 if (requested_length_unit == kim_length_unit_unused)
continue 517 if (requested_energy_unit == kim_energy_unit_unused)
continue 518 if (requested_charge_unit == kim_charge_unit_unused)
continue 519 if (requested_temperature_unit == kim_temperature_unit_unused)
continue 520 if (requested_time_unit == kim_time_unit_unused)
continue 523 call kim_set_units(model_create_handle, &
525 kim_energy_unit_ev, &
526 kim_charge_unit_unused, &
527 kim_temperature_unit_unused, &
528 kim_time_unit_unused, &
533 call kim_set_species_code(model_create_handle, &
534 kim_species_name_ar,
speccode, ierr2)
538 call kim_set_model_numbering(model_create_handle, &
539 kim_numbering_one_based, ierr2)
543 call kim_set_routine_pointer( &
544 model_create_handle, &
545 kim_model_routine_name_compute, kim_language_name_fortran, &
548 call kim_set_routine_pointer( &
549 model_create_handle, kim_model_routine_name_compute_arguments_create, &
553 call kim_set_routine_pointer( &
554 model_create_handle, kim_model_routine_name_compute_arguments_destroy, &
558 call kim_set_routine_pointer( &
559 model_create_handle, &
560 kim_model_routine_name_destroy, kim_language_name_fortran, &
568 call kim_set_model_buffer_pointer(model_create_handle, &
574 buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
577 call kim_set_influence_distance_pointer( &
578 model_create_handle, buf%influence_distance)
581 call kim_set_neighbor_list_pointers( &
582 model_create_handle, 1, buf%cutoff, &
583 buf%model_will_not_request_neighbors_of_noncontributing_particles)
588 call kim_log_entry(model_create_handle, kim_log_verbosity_error, &
589 "Unable to successfully initialize model")
595 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