kim-api  2.1.4-git+v2.1.3-git-1-g7847914a.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
ex_model_Ar_P_MLJ_Fortran.f90
Go to the documentation of this file.
1 !
2 ! CDDL HEADER START
3 !
4 ! The contents of this file are subject to the terms of the Common Development
5 ! and Distribution License Version 1.0 (the "License").
6 !
7 ! You can obtain a copy of the license at
8 ! http://www.opensource.org/licenses/CDDL-1.0. See the License for the
9 ! specific language governing permissions and limitations under the License.
10 !
11 ! When distributing Covered Code, include this CDDL HEADER in each file and
12 ! include the License file in a prominent location with the name LICENSE.CDDL.
13 ! If applicable, add the following below this CDDL HEADER, with the fields
14 ! enclosed by brackets "[]" replaced with your own identifying information:
15 !
16 ! Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved.
17 !
18 ! CDDL HEADER END
19 !
20 
21 !
22 ! Copyright (c) 2013--2019, Regents of the University of Minnesota.
23 ! All rights reserved.
24 !
25 ! Contributors:
26 ! Ryan S. Elliott
27 ! Ellad B. Tadmor
28 ! Valeriu Smirichinski
29 ! Stephen M. Whalen
30 !
31 
32 !****************************************************************************
33 !**
34 !** MODULE ex_model_Ar_P_MLJ_F03
35 !**
36 !** Modified Lennard-Jones pair potential (with smooth cutoff) model for Ar
37 !**
38 !** Reference: Ashcroft and Mermin
39 !**
40 !** Language: Fortran 2003
41 !**
42 !****************************************************************************
43 
44 
46 
47 use, intrinsic :: iso_c_binding
49 implicit none
50 
51 save
52 private
53 public compute_energy_forces, &
57  model_cutoff, &
58  speccode, &
59  buffer_type
60 
61 ! Below are the definitions and values of all Model parameters
62 integer(c_int), parameter :: cd = c_double ! used for literal constants
63 integer(c_int), parameter :: dim = 3 ! dimensionality of space
64 integer(c_int), parameter :: speccode = 1 ! internal species code
65 real(c_double), parameter :: model_cutoff = 8.15_cd ! cutoff radius
66  ! in angstroms
67 real(c_double), parameter :: model_cutsq = model_cutoff**2
68 
69 !-------------------------------------------------------------------------------
70 ! Below are the definitions and values of all additional model parameters
71 !
72 ! Recall that the Fortran 2003 format for declaring parameters is as follows:
73 !
74 ! integer(c_int), parameter :: parname = value ! This defines an integer
75 ! ! parameter called `parname'
76 ! ! with a value equal to
77 ! ! `value' (a number)
78 !
79 ! real(c_double), parameter :: parname = value ! This defines a real(c_double)
80 ! ! parameter called `parname'
81 ! ! with a value equal to
82 ! ! `value' (a number)
83 !-------------------------------------------------------------------------------
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 &
89  *lj_sigma**2)
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)
94 
95 type, bind(c) :: buffer_type
96  real(c_double) :: influence_distance
97  real(c_double) :: cutoff(1)
98  integer(c_int) :: &
99  model_will_not_request_neighbors_of_noncontributing_particles(1)
100 end type buffer_type
101 
102 contains
103 
104 !-------------------------------------------------------------------------------
105 !
106 ! Calculate pair potential phi(r)
107 !
108 !-------------------------------------------------------------------------------
109 recursive subroutine calc_phi(r,phi)
110 implicit none
111 
112 !-- Transferred variables
113 real(c_double), intent(in) :: r
114 real(c_double), intent(out) :: phi
115 
116 !-- Local variables
117 real(c_double) rsq,sor,sor6,sor12
118 
119 rsq = r*r ! r^2
120 sor = lj_sigma/r ! (sig/r)
121 sor6 = sor*sor*sor !
122 sor6 = sor6*sor6 ! (sig/r)^6
123 sor12= sor6*sor6 ! (sig/r)^12
124 if (r .gt. model_cutoff) then
125  ! Argument exceeds cutoff radius
126  phi = 0.0_cd
127 else
128  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
129 endif
130 
131 end subroutine calc_phi
132 
133 !-------------------------------------------------------------------------------
134 !
135 ! Calculate pair potential phi(r) and its derivative dphi(r)
136 !
137 !-------------------------------------------------------------------------------
138 recursive subroutine calc_phi_dphi(r,phi,dphi)
139 implicit none
140 
141 !-- Transferred variables
142 real(c_double), intent(in) :: r
143 real(c_double), intent(out) :: phi,dphi
144 
145 !-- Local variables
146 real(c_double) rsq,sor,sor6,sor12
147 
148 rsq = r*r ! r^2
149 sor = lj_sigma/r ! (sig/r)
150 sor6 = sor*sor*sor !
151 sor6 = sor6*sor6 ! (sig/r)^6
152 sor12= sor6*sor6 ! (sig/r)^12
153 if (r .gt. model_cutoff) then
154  ! Argument exceeds cutoff radius
155  phi = 0.0_cd
156  dphi = 0.0_cd
157 else
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
160 endif
161 
162 end subroutine calc_phi_dphi
163 
164 !-------------------------------------------------------------------------------
165 !
166 ! Compute energy and forces on particles from the positions.
167 !
168 !-------------------------------------------------------------------------------
169 recursive subroutine compute_energy_forces(model_compute_handle, &
170  model_compute_arguments_handle, ierr) bind(c)
171 implicit none
172 
173 !-- Transferred variables
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
178 
179 !-- Local variables
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
184 
185 !-- KIM variables
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(:)
195 
196 ! Unpack data from KIM object
197 !
198 ierr = 0
199 call kim_get_argument_pointer( &
200  model_compute_arguments_handle, &
201  kim_compute_argument_name_number_of_particles, n, ierr2)
202 ierr = ierr + ierr2
203 call kim_get_argument_pointer( &
204  model_compute_arguments_handle, &
205  kim_compute_argument_name_particle_species_codes, n, particlespeciescodes, &
206  ierr2)
207 ierr = ierr + ierr2
208 call kim_get_argument_pointer( &
209  model_compute_arguments_handle, &
210  kim_compute_argument_name_particle_contributing, n, particlecontributing, &
211  ierr2)
212 ierr = ierr + ierr2
213 call kim_get_argument_pointer( &
214  model_compute_arguments_handle, &
215  kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
216 ierr = ierr + ierr2
217 call kim_get_argument_pointer( &
218  model_compute_arguments_handle, &
219  kim_compute_argument_name_partial_energy, energy, ierr2)
220 ierr = ierr + ierr2
221 call kim_get_argument_pointer( &
222  model_compute_arguments_handle, &
223  kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
224 ierr = ierr + ierr2
225 call kim_get_argument_pointer( &
226  model_compute_arguments_handle, &
227  kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
228 ierr = ierr + ierr2
229 call kim_get_argument_pointer( &
230  model_compute_arguments_handle, &
231  kim_compute_argument_name_partial_virial, 6, virial, ierr2)
232 ierr = ierr + ierr2
233 if (ierr /= 0) then
234  call kim_log_entry(model_compute_arguments_handle, &
235  kim_log_verbosity_error, "get data")
236  ierr = 1
237  return
238 endif
239 
240 ! Check to see if we have been asked to compute the forces, energyperpart,
241 ! energy and virial
242 !
243 if (associated(energy)) then
244  comp_energy = 1
245 else
246  comp_energy = 0
247 end if
248 if (associated(force)) then
249  comp_force = 1
250 else
251  comp_force = 0
252 end if
253 if (associated(enepot)) then
254  comp_enepot = 1
255 else
256  comp_enepot = 0
257 end if
258 if (associated(virial)) then
259  comp_virial = 1
260 else
261  comp_virial = 0
262 end if
263 
264 ! Check to be sure that the species are correct
265 !
266 ierr = 1 ! assume an error
267 do i = 1,n
268  if (particlespeciescodes(i).ne.speccode) then
269  call kim_log_entry(model_compute_handle, &
270  kim_log_verbosity_error, "Unexpected species code detected")
271  ierr = 1
272  return
273  endif
274 enddo
275 ierr = 0 ! everything is ok
276 
277 ! Initialize potential energies, forces, virial term
278 !
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
283 
284 !
285 ! Compute energy and forces
286 !
287 
288 ! Loop over particles and compute energy and forces
289 !
290 do i = 1, n
291  if (particlecontributing(i) == 1) then
292  ! Set up neighbor list for next particle
293  call kim_get_neighbor_list( &
294  model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
295  if (ierr /= 0) then
296  ! some sort of problem, exit
297  call kim_log_entry(&
298  model_compute_arguments_handle, kim_log_verbosity_error, &
299  "GetNeighborList failed")
300  ierr = 1
301  return
302  endif
303 
304  ! Loop over the neighbors of particle i
305  !
306  do jj = 1, numnei
307 
308  j = nei1part(jj) ! get neighbor ID
309 
310  ! compute relative position vector
311  !
312  rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
313 
314  ! compute energy and forces
315  !
316  rsqij = dot_product(rij,rij) ! compute square distance
317  if ( rsqij .lt. model_cutsq ) then ! particles are interacting?
318 
319  r = sqrt(rsqij) ! compute distance
320  if (comp_force.eq.1.or.comp_virial.eq.1) then
321  call calc_phi_dphi(r,phi,dphi) ! compute pair potential
322  ! and it derivative
323  deidr = 0.5_cd*dphi
324  else
325  call calc_phi(r,phi) ! compute just pair potential
326  endif
327 
328  ! contribution to energy
329  !
330  if (comp_enepot.eq.1) then
331  enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
332  endif
333  if (comp_energy.eq.1) then
334  energy = energy + 0.5_cd*phi
335  endif
336 
337  ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
338  !
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
346  endif
347 
348  ! contribution to forces
349  !
350  if (comp_force.eq.1) then
351  force(:,i) = force(:,i) + deidr*rij/r ! accumulate force on i
352  force(:,j) = force(:,j) - deidr*rij/r ! accumulate force on j
353  endif
354 
355  endif
356 
357  enddo ! loop on jj
358 
359  endif ! if particleContributing
360 
361 enddo ! do i
362 
363 ! Everything is great
364 !
365 ierr = 0
366 return
367 
368 end subroutine compute_energy_forces
369 
370 !-------------------------------------------------------------------------------
371 !
372 ! Model destroy routine (REQUIRED)
373 !
374 !-------------------------------------------------------------------------------
375 recursive subroutine model_destroy_func(model_destroy_handle, ierr) bind(c)
376  use, intrinsic :: iso_c_binding
377  implicit none
378 
379  !-- Transferred variables
380  type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle
381  integer(c_int), intent(out) :: ierr
382 
383  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
384 
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")
389  deallocate(buf)
390  ierr = 0 ! everything is good
391 end subroutine model_destroy_func
392 
393 !-------------------------------------------------------------------------------
394 !
395 ! Model compute arguments create routine (REQUIRED)
396 !
397 !-------------------------------------------------------------------------------
398 recursive subroutine model_compute_arguments_create(model_compute_handle, &
399  model_compute_arguments_create_handle, ierr) bind(c)
400  use, intrinsic :: iso_c_binding
401  implicit none
402 
403  !-- Transferred variables
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
408 
409  integer(c_int) :: ierr2
410 
411  ! avoid unused dummy argument warnings
412  if (model_compute_handle .eq. kim_model_compute_null_handle) continue
413 
414  ierr = 0
415  ierr2 = 0
416 
417  ! register arguments
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)
422  ierr = ierr + 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)
427  ierr = ierr + 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)
432  ierr = ierr + 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)
437  ierr = ierr + ierr2
438 
439  ! register call backs
440  ! NONE
441 
442  if (ierr /= 0) then
443  ierr = 1
444  call kim_log_entry( &
445  model_compute_arguments_create_handle, kim_log_verbosity_error, &
446  "Unable to successfully create compute_arguments object")
447  endif
448 
449  ierr = 0
450  return
451 end subroutine model_compute_arguments_create
452 
453 !-------------------------------------------------------------------------------
454 !
455 ! Model compute arguments destroy routine (REQUIRED)
456 !
457 !-------------------------------------------------------------------------------
458 recursive subroutine model_compute_arguments_destroy(model_compute_handle, &
459  model_compute_arguments_destroy_handle, ierr) bind(c)
460  use, intrinsic :: iso_c_binding
461  implicit none
462 
463  !-- Transferred variables
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
468 
469  ! avoid unused dummy argument warnings
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
473 
474  ierr = 0
475  return
476 end subroutine model_compute_arguments_destroy
477 
478 end module ex_model_ar_p_mlj_f03
479 
480 !-------------------------------------------------------------------------------
481 !
482 ! Model create routine (REQUIRED)
483 !
484 !-------------------------------------------------------------------------------
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
491 implicit none
492 
493 !-- Transferred variables
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
501 
502 !-- KIM variables
503 integer(c_int) :: ierr2
504 type(buffer_type), pointer :: buf
505 
506 ierr = 0
507 ierr2 = 0
508 
509 ! avoid unsed dummy argument warnings
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
515 
516 ! set units
517 call kim_set_units(model_create_handle, &
518  kim_length_unit_a, &
519  kim_energy_unit_ev, &
520  kim_charge_unit_unused, &
521  kim_temperature_unit_unused, &
522  kim_time_unit_unused, &
523  ierr2)
524 ierr = ierr + ierr2
525 
526 ! register species
527 call kim_set_species_code(model_create_handle, &
528  kim_species_name_ar, speccode, ierr2)
529 ierr = ierr + ierr2
530 
531 ! register numbering
532 call kim_set_model_numbering(model_create_handle, &
533  kim_numbering_one_based, ierr2);
534 ierr = ierr + ierr2
535 
536 ! register function pointers
537 call kim_set_routine_pointer(model_create_handle, &
538  kim_model_routine_name_compute, kim_language_name_fortran, &
539  1, c_funloc(compute_energy_forces), ierr2)
540 ierr = ierr + ierr2
541 call kim_set_routine_pointer( &
542  model_create_handle, kim_model_routine_name_compute_arguments_create, &
543  kim_language_name_fortran, 1, c_funloc(model_compute_arguments_create), ierr2)
544 ierr = ierr + ierr2
545 call kim_set_routine_pointer( &
546  model_create_handle, kim_model_routine_name_compute_arguments_destroy, &
547  kim_language_name_fortran, 1, c_funloc(model_compute_arguments_destroy), &
548  ierr2)
549 ierr = ierr + ierr2
550 call kim_set_routine_pointer(model_create_handle, &
551  kim_model_routine_name_destroy, kim_language_name_fortran, &
552  1, c_funloc(model_destroy_func), ierr2)
553 ierr = ierr + ierr2
554 
555 ! allocate buffer
556 allocate( buf )
557 
558 ! store model buffer in KIM object
559 call kim_set_model_buffer_pointer(model_create_handle, &
560  c_loc(buf))
561 
562 ! set buffer values
563 buf%influence_distance = model_cutoff
564 buf%cutoff = model_cutoff
565 buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
566 
567 ! register influence distance
568 call kim_set_influence_distance_pointer( &
569  model_create_handle, buf%influence_distance)
570 
571 ! register cutoff
572 call kim_set_neighbor_list_pointers(model_create_handle, &
573  1, buf%cutoff, &
574  buf%model_will_not_request_neighbors_of_noncontributing_particles)
575 
576 if (ierr /= 0) then
577  ierr = 1
578  deallocate( buf )
579  call kim_log_entry(model_create_handle, &
580  kim_log_verbosity_error, "Unable to successfully initialize model")
581 endif
582 
583 ierr = 0
584 return
585 
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