53 use,
intrinsic :: iso_c_binding
69 integer(c_int),
parameter :: cd = c_double
70 integer(c_int),
parameter :: dim = 3
88 real(c_double),
parameter :: lj_spring = 0.00051226_cd
89 real(c_double),
parameter :: lj_sigma = 5.26_cd / 1.74724_cd
98 type, bind(c) :: buffer_type
99 real(c_double) :: influence_distance
100 real(c_double) :: cutoff(2)
102 model_will_not_request_neighbors_of_noncontributing_particles(2)
113 recursive subroutine calc_phi(r, phi, dphi, calc_deriv)
117 real(c_double),
intent(in) :: r
118 real(c_double),
intent(out) :: phi
119 real(c_double),
intent(out) :: dphi
120 logical,
intent(in) :: calc_deriv
123 real(c_double) rsq, sor, sor6, sor12
127 sor6 = sor * sor * sor
133 if (calc_deriv) dphi = 0.0_cd
135 phi = 4.0_cd * (sor12 - sor6)
136 if (calc_deriv) dphi = 24.0_cd * (-2.0_cd * sor12 + sor6) / r
146 recursive subroutine calc_spring_energyamp(model_compute_arguments_handle, &
147 atom, coor, eps, ierr)
151 type(kim_model_compute_arguments_handle_type), &
152 intent(in) :: model_compute_arguments_handle
153 integer(c_int),
intent(in) :: atom
154 real(c_double),
intent(in) :: coor(:, :)
155 real(c_double),
intent(out) :: eps
156 integer(c_int),
intent(out) :: ierr
160 real(c_double) rrel(dim)
161 real(c_double) rsqrel
162 integer(c_int) numneishort
163 integer(c_int),
pointer :: neishort(:)
166 call kim_get_neighbor_list( &
167 model_compute_arguments_handle, 1, atom, numneishort, neishort, ierr)
168 if (ierr /= 0)
return 171 do kk = 1, numneishort
173 rrel(:) = coor(:, k) - coor(:, atom)
174 rsqrel = dot_product(rrel, rrel)
175 if (rsqrel < model_cutsq1) eps = eps + rsqrel
177 eps = 0.5_cd * lj_spring * eps
179 end subroutine calc_spring_energyamp
186 recursive subroutine calc_spring_force(model_compute_arguments_handle, atom, &
187 coor, eps, phi, force, ierr)
191 type(kim_model_compute_arguments_handle_type), &
192 intent(in) :: model_compute_arguments_handle
193 integer(c_int),
intent(in) :: atom
194 real(c_double),
intent(in) :: coor(:, :)
195 real(c_double),
intent(in) :: eps
196 real(c_double),
intent(in) :: phi
197 real(c_double),
intent(inout) :: force(:, :)
198 integer(c_int),
intent(out) :: ierr
202 real(c_double) rrel(dim), dforce(dim)
203 real(c_double) rsqrel
204 integer(c_int) numneishort
205 integer(c_int),
pointer :: neishort(:)
208 call kim_get_neighbor_list( &
209 model_compute_arguments_handle, 1, atom, numneishort, neishort, ierr)
210 if (ierr /= 0)
return 214 do kk = 1, numneishort
216 rrel(:) = coor(:, k) - coor(:, atom)
217 rsqrel = dot_product(rrel, rrel)
218 if (rsqrel < model_cutsq1)
then 219 dforce(:) = 0.5_cd * eps * lj_spring * rrel(:) * phi
220 force(:, atom) = force(:, atom) + dforce(:)
221 force(:, k) = force(:, k) - dforce(:)
225 end subroutine calc_spring_force
233 model_compute_handle, model_compute_arguments_handle, ierr) bind(c)
237 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
238 type(kim_model_compute_arguments_handle_type),
intent(in) :: &
239 model_compute_arguments_handle
240 integer(c_int),
intent(out) :: ierr
243 real(c_double) :: Rij(dim), dforce(dim)
244 real(c_double) :: r, Rsqij, phi, dphi, dEidr, epsi, epsj
245 integer(c_int) :: i, j, jj, comp_force, comp_enepot, comp_energy
247 integer(c_int) :: numnei
248 integer(c_int) :: ierr2
249 logical :: calc_deriv
252 integer(c_int),
pointer :: N
253 real(c_double),
pointer :: energy
254 real(c_double),
pointer :: coor(:, :)
255 real(c_double),
pointer :: force(:, :)
256 real(c_double),
pointer :: enepot(:)
257 integer(c_int),
pointer :: nei1part(:)
258 integer(c_int),
pointer :: particleSpeciesCodes(:)
259 integer(c_int),
pointer :: particleContributing(:)
265 call kim_get_argument_pointer( &
266 model_compute_arguments_handle, &
267 kim_compute_argument_name_number_of_particles, n, ierr2)
269 call kim_get_argument_pointer( &
270 model_compute_arguments_handle, &
271 kim_compute_argument_name_particle_species_codes, n, &
272 particlespeciescodes, ierr2)
274 call kim_get_argument_pointer( &
275 model_compute_arguments_handle, &
276 kim_compute_argument_name_particle_contributing, n, &
277 particlecontributing, ierr2)
279 call kim_get_argument_pointer( &
280 model_compute_arguments_handle, &
281 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
283 call kim_get_argument_pointer( &
284 model_compute_arguments_handle, &
285 kim_compute_argument_name_partial_energy, energy, ierr2)
287 call kim_get_argument_pointer( &
288 model_compute_arguments_handle, &
289 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
291 call kim_get_argument_pointer( &
292 model_compute_arguments_handle, &
293 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
300 call kim_log_entry(model_compute_arguments_handle, &
301 kim_log_verbosity_error,
"get data")
308 if (
associated(energy))
then 313 if (
associated(force))
then 318 if (
associated(enepot))
then 328 calc_deriv = comp_force == 1
334 if (particlespeciescodes(i) /=
speccode)
then 335 call kim_log_entry( &
336 model_compute_handle, kim_log_verbosity_error, &
337 "Unexpected species code detected")
345 if (comp_enepot == 1) enepot = 0.0_cd
346 if (comp_energy == 1) energy = 0.0_cd
347 if (comp_force == 1) force = 0.0_cd
349 if (calc_deriv) deidr = 0.0_cd
358 if (particlecontributing(i) == 1)
then 360 call kim_get_neighbor_list( &
361 model_compute_arguments_handle, 2, i, numnei, nei1part, ierr)
364 call kim_log_entry( &
365 model_compute_arguments_handle, kim_log_verbosity_error, &
366 "GetNeighborList failed")
372 call calc_spring_energyamp(model_compute_arguments_handle, &
376 call kim_log_entry( &
377 model_compute_handle, kim_log_verbosity_error, &
378 "GetNeighborList failed")
390 call calc_spring_energyamp(model_compute_arguments_handle, j, coor, &
394 call kim_log_entry( &
395 model_compute_handle, kim_log_verbosity_error, &
396 "GetNeighborList failed")
403 rij(:) = coor(:, j) - coor(:, i)
407 rsqij = dot_product(rij, rij)
408 if (rsqij < model_cutsq2)
then 411 call calc_phi(r, phi, dphi, calc_deriv)
412 if (calc_deriv) deidr = 0.5_cd * dphi
416 if (comp_enepot == 1)
then 417 enepot(i) = enepot(i) + 0.5_cd * epsi * epsj * phi
419 if (comp_energy == 1)
then 420 energy = energy + 0.5_cd * epsi * epsj * phi
439 if (comp_force == 1)
then 441 call calc_spring_force(model_compute_arguments_handle, i, coor, &
442 epsj, phi, force, ierr)
445 call kim_log_entry( &
446 model_compute_handle, kim_log_verbosity_error, &
447 "GetNeighborList failed")
452 call calc_spring_force(model_compute_arguments_handle, j, coor, &
453 epsi, phi, force, ierr)
456 call kim_log_entry( &
457 model_compute_handle, kim_log_verbosity_error, &
458 "GetNeighborList failed")
463 dforce(:) = epsi * epsj * deidr * rij(:) / r
464 force(:, i) = force(:, i) + dforce(:)
465 force(:, j) = force(:, j) - dforce(:)
489 use,
intrinsic :: iso_c_binding
493 type(kim_model_destroy_handle_type),
intent(inout) :: model_destroy_handle
494 integer(c_int),
intent(out) :: ierr
496 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
498 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
499 call c_f_pointer(pbuf, buf)
500 call kim_log_entry(model_destroy_handle, kim_log_verbosity_information, &
501 "deallocating model buffer")
512 model_compute_handle, model_compute_arguments_create_handle, ierr) bind(c)
513 use,
intrinsic :: iso_c_binding
517 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
518 type(kim_model_compute_arguments_create_handle_type),
intent(inout) :: &
519 model_compute_arguments_create_handle
520 integer(c_int),
intent(out) :: ierr
522 integer(c_int) :: ierr2
525 if (model_compute_handle == kim_model_compute_null_handle)
continue 531 call kim_set_argument_support_status( &
532 model_compute_arguments_create_handle, &
533 kim_compute_argument_name_partial_energy, &
534 kim_support_status_optional, ierr2)
536 call kim_set_argument_support_status( &
537 model_compute_arguments_create_handle, &
538 kim_compute_argument_name_partial_forces, &
539 kim_support_status_optional, ierr2)
541 call kim_set_argument_support_status( &
542 model_compute_arguments_create_handle, &
543 kim_compute_argument_name_partial_particle_energy, &
544 kim_support_status_optional, ierr2)
557 call kim_log_entry( &
558 model_compute_arguments_create_handle, &
559 kim_log_verbosity_error, &
560 "Unable to successfully create compute_arguments object")
572 model_compute_handle, model_compute_arguments_destroy_handle, ierr) bind(c)
573 use,
intrinsic :: iso_c_binding
577 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
578 type(kim_model_compute_arguments_destroy_handle_type),
intent(inout) :: &
579 model_compute_arguments_destroy_handle
580 integer(c_int),
intent(out) :: ierr
583 if (model_compute_handle == kim_model_compute_null_handle)
continue 584 if (model_compute_arguments_destroy_handle == &
585 kim_model_compute_arguments_destroy_null_handle)
continue 598 recursive subroutine model_create_routine( &
599 model_create_handle, requested_length_unit, requested_energy_unit, &
600 requested_charge_unit, requested_temperature_unit, requested_time_unit, &
602 use,
intrinsic :: iso_c_binding
608 type(kim_model_create_handle_type),
intent(inout) :: model_create_handle
609 type(kim_length_unit_type),
intent(in),
value :: requested_length_unit
610 type(kim_energy_unit_type),
intent(in),
value :: requested_energy_unit
611 type(kim_charge_unit_type),
intent(in),
value :: requested_charge_unit
612 type(kim_temperature_unit_type),
intent(in),
value :: &
613 requested_temperature_unit
614 type(kim_time_unit_type),
intent(in),
value :: requested_time_unit
615 integer(c_int),
intent(out) :: ierr
618 integer(c_int) :: ierr2
619 type(buffer_type),
pointer :: buf
625 if (requested_length_unit == kim_length_unit_unused)
continue 626 if (requested_energy_unit == kim_energy_unit_unused)
continue 627 if (requested_charge_unit == kim_charge_unit_unused)
continue 628 if (requested_temperature_unit == kim_temperature_unit_unused)
continue 629 if (requested_time_unit == kim_time_unit_unused)
continue 632 call kim_set_units(model_create_handle, &
634 kim_energy_unit_ev, &
635 kim_charge_unit_unused, &
636 kim_temperature_unit_unused, &
637 kim_time_unit_unused, &
642 call kim_set_species_code(model_create_handle, &
643 kim_species_name_ar,
speccode, ierr2)
647 call kim_set_model_numbering(model_create_handle, &
648 kim_numbering_one_based, ierr2)
652 call kim_set_routine_pointer( &
653 model_create_handle, &
654 kim_model_routine_name_compute, kim_language_name_fortran, &
657 call kim_set_routine_pointer( &
658 model_create_handle, &
659 kim_model_routine_name_compute_arguments_create, &
660 kim_language_name_fortran, &
663 call kim_set_routine_pointer( &
664 model_create_handle, &
665 kim_model_routine_name_compute_arguments_destroy, &
666 kim_language_name_fortran, &
669 call kim_set_routine_pointer( &
670 model_create_handle, &
671 kim_model_routine_name_destroy, kim_language_name_fortran, 1, &
679 call kim_set_model_buffer_pointer(model_create_handle, &
686 buf%model_will_not_request_neighbors_of_noncontributing_particles(1) = 0
687 buf%model_will_not_request_neighbors_of_noncontributing_particles(2) = 1
690 call kim_set_influence_distance_pointer( &
691 model_create_handle, buf%influence_distance)
694 call kim_set_neighbor_list_pointers( &
695 model_create_handle, 2, buf%cutoff, &
696 buf%model_will_not_request_neighbors_of_noncontributing_particles)
701 call kim_log_entry( &
702 model_create_handle, kim_log_verbosity_error, &
703 "Unable to successfully initialize model")
708 end subroutine model_create_routine
recursive subroutine, public model_compute_arguments_create(model_compute_handle, model_compute_arguments_create_handle, ierr)
recursive subroutine, public model_compute_arguments_destroy(model_compute_handle, model_compute_arguments_destroy_handle, ierr)
real(c_double), parameter, public model_cutoff1
real(c_double), parameter, public model_cutoff2
integer(c_int), parameter, public speccode
static void calc_phi(double const *epsilon, double const *C, double const *Rzero, double const *shift, double const cutoff, double const r, double *phi)
recursive subroutine, public compute_energy_forces(model_compute_handle, model_compute_arguments_handle, ierr)
recursive subroutine, public model_destroy_func(model_destroy_handle, ierr)