47 use,
intrinsic :: iso_c_binding
63 integer(c_int),
parameter :: cd = c_double
64 integer(c_int),
parameter :: DIM = 3
65 integer(c_int),
parameter :: speccode = 1
66 real(c_double),
parameter :: model_cutoff = 8.15_cd
68 real(c_double),
parameter :: model_cutsq = model_cutoff**2
85 real(c_double),
parameter :: lj_epsilon = 0.0104_cd
86 real(c_double),
parameter :: lj_sigma = 3.40_cd
87 real(c_double),
parameter :: lj_cutnorm = model_cutoff/lj_sigma
88 real(c_double),
parameter :: lj_A = 12.0_cd*lj_epsilon*(-26.0_cd &
89 + 7.0_cd*lj_cutnorm**6)/(lj_cutnorm**14 &
91 real(c_double),
parameter :: lj_B = 96.0_cd*lj_epsilon*(7.0_cd &
92 - 2.0_cd*lj_cutnorm**6)/(lj_cutnorm**13*lj_sigma)
93 real(c_double),
parameter :: lj_C = 28.0_cd*lj_epsilon*(-13.0_cd &
94 + 4.0_cd*lj_cutnorm**6)/(lj_cutnorm**12)
96 type, bind(c) :: buffer_type
97 real(c_double) :: influence_distance
98 real(c_double) :: cutoff(1)
99 integer(c_int) :: padding_neighbor_hints(1)
100 integer(c_int) :: half_list_hints(1)
110 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
125 if (r .gt. model_cutoff)
then 129 phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
132 end subroutine calc_phi
139 subroutine calc_phi_dphi(r,phi,dphi)
143 real(c_double),
intent(in) :: r
144 real(c_double),
intent(out) :: phi,dphi
147 real(c_double) rsq,sor,sor6,sor12
154 if (r .gt. model_cutoff)
then 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 + 2.0_cd*lj_a*r + lj_b
163 end subroutine calc_phi_dphi
170 #include "kim_model_compute_log_macros.fd" 172 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,comp_virial, comp_energy
185 integer(c_int) :: ierr2
188 integer(c_int),
pointer :: n
189 real(c_double),
pointer :: energy
190 real(c_double),
pointer :: coor(:,:)
191 real(c_double),
pointer :: force(:,:)
192 real(c_double),
pointer :: enepot(:)
193 integer(c_int),
pointer :: nei1part(:)
194 integer(c_int),
pointer :: particlespeciescodes(:)
195 integer(c_int),
pointer :: particlecontributing(:)
196 real(c_double),
pointer :: virial(:)
198 kim_log_file = __file__
203 call kim_model_compute_arguments_get_argument_pointer( &
204 model_compute_arguments_handle, &
205 kim_compute_argument_name_number_of_particles, n, ierr2)
207 call kim_model_compute_arguments_get_argument_pointer( &
208 model_compute_arguments_handle, &
209 kim_compute_argument_name_particle_species_codes, n, particlespeciescodes, &
212 call kim_model_compute_arguments_get_argument_pointer( &
213 model_compute_arguments_handle, &
214 kim_compute_argument_name_particle_contributing, n, particlecontributing, &
217 call kim_model_compute_arguments_get_argument_pointer( &
218 model_compute_arguments_handle, &
219 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
221 call kim_model_compute_arguments_get_argument_pointer( &
222 model_compute_arguments_handle, &
223 kim_compute_argument_name_partial_energy, energy, ierr2)
225 call kim_model_compute_arguments_get_argument_pointer( &
226 model_compute_arguments_handle, &
227 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
229 call kim_model_compute_arguments_get_argument_pointer( &
230 model_compute_arguments_handle, &
231 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
233 call kim_model_compute_arguments_get_argument_pointer( &
234 model_compute_arguments_handle, &
235 kim_compute_argument_name_partial_virial, 6, virial, ierr2)
238 kim_log_message =
"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).ne.speccode)
then 272 kim_log_message =
"Unexpected species code detected" 281 if (comp_enepot.eq.1) enepot = 0.0_cd
282 if (comp_energy.eq.1) energy = 0.0_cd
283 if (comp_force.eq.1) force = 0.0_cd
284 if (comp_virial.eq.1) virial = 0.0_cd
293 if (particlecontributing(i) == 1)
then 295 call kim_model_compute_arguments_get_neighbor_list( &
296 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
299 kim_log_message =
"GetNeighborList failed" 313 rij(:) = coor(:,j) - coor(:,i)
317 rsqij = dot_product(rij,rij)
318 if ( rsqij .lt. model_cutsq )
then 321 if (comp_force.eq.1.or.comp_virial.eq.1)
then 322 call calc_phi_dphi(r,phi,dphi)
331 if (comp_enepot.eq.1)
then 332 enepot(i) = enepot(i) + 0.5_cd*phi
334 if (comp_energy.eq.1)
then 335 energy = energy + 0.5_cd*phi
340 if (comp_virial.eq.1)
then 341 virial(1) = virial(1) + rij(1)*rij(1)*deidr/r
342 virial(2) = virial(2) + rij(2)*rij(2)*deidr/r
343 virial(3) = virial(3) + rij(3)*rij(3)*deidr/r
344 virial(4) = virial(4) + rij(2)*rij(3)*deidr/r
345 virial(5) = virial(5) + rij(1)*rij(3)*deidr/r
346 virial(6) = virial(6) + rij(1)*rij(2)*deidr/r
351 if (comp_force.eq.1)
then 352 force(:,i) = force(:,i) + deidr*rij/r
353 force(:,j) = force(:,j) - deidr*rij/r
376 #include "kim_model_destroy_log_macros.fd" 378 use,
intrinsic :: iso_c_binding
382 type(kim_model_destroy_handle_type),
intent(inout) :: model_destroy_handle
383 integer(c_int),
intent(out) :: ierr
385 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
387 kim_log_file = __file__
389 call kim_model_destroy_get_model_buffer_pointer(model_destroy_handle, pbuf)
390 call c_f_pointer(pbuf, buf)
391 kim_log_message =
"deallocating model buffer" 402 #include "kim_model_refresh_log_macros.fd" 404 use,
intrinsic :: iso_c_binding
408 type(kim_model_refresh_handle_type),
intent(inout) :: model_refresh_handle
409 integer(c_int),
intent(out) :: ierr
411 type(buffer_type),
pointer :: buf;
type(c_ptr) :: pbuf
413 kim_log_file = __file__
415 call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
416 call c_f_pointer(pbuf, buf)
418 kim_log_message =
"Resettings influence distance and cutoffs" 420 call kim_model_refresh_set_influence_distance_pointer( &
421 model_refresh_handle, buf%cutoff(1))
422 call kim_model_refresh_set_neighbor_list_pointers( &
423 model_refresh_handle, 1, buf%cutoff, buf%padding_neighbor_hints, &
434 #include "kim_model_compute_arguments_create_log_macros.fd" 436 model_compute_arguments_create_handle, ierr) bind(c)
437 use,
intrinsic :: iso_c_binding
439 log_entry=>kim_model_compute_arguments_create_log_entry
443 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
444 type(kim_model_compute_arguments_create_handle_type),
intent(inout) :: &
445 model_compute_arguments_create_handle
446 integer(c_int),
intent(out) :: ierr
448 integer(c_int) :: ierr2
454 call kim_model_compute_arguments_create_set_argument_support_status( &
455 model_compute_arguments_create_handle, &
456 kim_compute_argument_name_partial_energy, &
457 kim_support_status_optional, ierr2)
459 call kim_model_compute_arguments_create_set_argument_support_status( &
460 model_compute_arguments_create_handle, &
461 kim_compute_argument_name_partial_forces, &
462 kim_support_status_optional, ierr2)
464 call kim_model_compute_arguments_create_set_argument_support_status( &
465 model_compute_arguments_create_handle, &
466 kim_compute_argument_name_partial_particle_energy, &
467 kim_support_status_optional, ierr2)
469 call kim_model_compute_arguments_create_set_argument_support_status( &
470 model_compute_arguments_create_handle, &
471 kim_compute_argument_name_partial_virial, &
472 kim_support_status_optional, ierr2)
480 kim_log_message =
"Unable to successfully create compute_arguments object" 492 #include "kim_model_compute_arguments_destroy_log_macros.fd" 494 model_compute_arguments_destroy_handle, ierr) bind(c)
495 use,
intrinsic :: iso_c_binding
497 log_entry=>kim_model_compute_arguments_destroy_log_entry
501 type(kim_model_compute_handle_type),
intent(in) :: model_compute_handle
502 type(kim_model_compute_arguments_destroy_handle_type),
intent(inout) :: &
503 model_compute_arguments_destroy_handle
504 integer(c_int),
intent(out) :: ierr
506 integer(c_int) :: ierr2
523 #include "kim_model_create_log_macros.fd" 524 subroutine model_create_routine(model_create_handle, requested_length_unit, &
525 requested_energy_unit, requested_charge_unit, requested_temperature_unit, &
526 requested_time_unit, ierr) bind(c)
527 use,
intrinsic :: iso_c_binding
533 type(kim_model_create_handle_type),
intent(inout) :: model_create_handle
534 type(kim_length_unit_type),
intent(in) :: requested_length_unit
535 type(kim_energy_unit_type),
intent(in) :: requested_energy_unit
536 type(kim_charge_unit_type),
intent(in) :: requested_charge_unit
537 type(kim_temperature_unit_type),
intent(in) :: requested_temperature_unit
538 type(kim_time_unit_type),
intent(in) :: requested_time_unit
539 integer(c_int),
intent(out) :: ierr
542 integer(c_int) :: ierr2
543 type(buffer_type),
pointer :: buf
545 kim_log_file = __file__
551 call kim_model_create_set_units(model_create_handle, &
553 kim_energy_unit_ev, &
554 kim_charge_unit_unused, &
555 kim_temperature_unit_unused, &
556 kim_time_unit_unused, &
561 call kim_model_create_set_species_code(model_create_handle, &
562 kim_species_name_ar, speccode, ierr2)
566 call kim_model_create_set_model_numbering(model_create_handle, &
567 kim_numbering_one_based, ierr2);
571 call kim_model_create_set_compute_pointer(model_create_handle, &
574 call kim_model_create_set_compute_arguments_create_pointer( &
575 model_create_handle, kim_language_name_fortran, &
578 call kim_model_create_set_compute_arguments_destroy_pointer( &
579 model_create_handle, kim_language_name_fortran, &
582 call kim_model_create_set_destroy_pointer(model_create_handle, &
585 call kim_model_create_set_refresh_pointer( &
586 model_create_handle, kim_language_name_fortran, &
594 call kim_model_create_set_model_buffer_pointer(model_create_handle, &
598 buf%influence_distance = model_cutoff
599 buf%cutoff = model_cutoff
600 buf%padding_neighbor_hints = 1
601 buf%half_list_hints = 0
604 call kim_model_create_set_influence_distance_pointer( &
605 model_create_handle, buf%influence_distance)
608 call kim_model_create_set_neighbor_list_pointers(model_create_handle, &
609 1, buf%cutoff, buf%padding_neighbor_hints, buf%half_list_hints)
614 kim_log_message =
"Unable to successfully initialize model" 620 end subroutine model_create_routine
subroutine, public model_destroy_func(model_destroy_handle, ierr)
subroutine, public model_refresh_func(model_refresh_handle, ierr)
subroutine, public model_compute_arguments_destroy(model_compute_handle, model_compute_arguments_destroy_handle, ierr)
subroutine, public model_compute_arguments_create(model_compute_handle, model_compute_arguments_create_handle, ierr)
integer(c_int) function, public compute_energy_forces(pkim)