KIM API V2
ex_model_Ar_P_MLJ_F03.F03
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--2018, 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, &
58  model_cutoff, &
59  speccode, &
60  buffer_type
61 
62 ! Below are the definitions and values of all Model parameters
63 integer(c_int), parameter :: cd = c_double ! used for literal constants
64 integer(c_int), parameter :: DIM = 3 ! dimensionality of space
65 integer(c_int), parameter :: speccode = 1 ! internal species code
66 real(c_double), parameter :: model_cutoff = 8.15_cd ! cutoff radius
67  ! in angstroms
68 real(c_double), parameter :: model_cutsq = model_cutoff**2
69 
70 !-------------------------------------------------------------------------------
71 ! Below are the definitions and values of all additional model parameters
72 !
73 ! Recall that the Fortran 2003 format for declaring parameters is as follows:
74 !
75 ! integer(c_int), parameter :: parname = value ! This defines an integer
76 ! ! parameter called `parname'
77 ! ! with a value equal to
78 ! ! `value' (a number)
79 !
80 ! real(c_double), parameter :: parname = value ! This defines a real(c_double)
81 ! ! parameter called `parname'
82 ! ! with a value equal to
83 ! ! `value' (a number)
84 !-------------------------------------------------------------------------------
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 &
90  *lj_sigma**2)
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)
95 
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)
101 end type buffer_type
102 
103 contains
104 
105 !-------------------------------------------------------------------------------
106 !
107 ! Calculate pair potential phi(r)
108 !
109 !-------------------------------------------------------------------------------
110 subroutine calc_phi(r,phi)
111 implicit none
112 
113 !-- Transferred variables
114 real(c_double), intent(in) :: r
115 real(c_double), intent(out) :: phi
116 
117 !-- Local variables
118 real(c_double) rsq,sor,sor6,sor12
119 
120 rsq = r*r ! r^2
121 sor = lj_sigma/r ! (sig/r)
122 sor6 = sor*sor*sor !
123 sor6 = sor6*sor6 ! (sig/r)^6
124 sor12= sor6*sor6 ! (sig/r)^12
125 if (r .gt. model_cutoff) then
126  ! Argument exceeds cutoff radius
127  phi = 0.0_cd
128 else
129  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
130 endif
131 
132 end subroutine calc_phi
133 
134 !-------------------------------------------------------------------------------
135 !
136 ! Calculate pair potential phi(r) and its derivative dphi(r)
137 !
138 !-------------------------------------------------------------------------------
139 subroutine calc_phi_dphi(r,phi,dphi)
140 implicit none
141 
142 !-- Transferred variables
143 real(c_double), intent(in) :: r
144 real(c_double), intent(out) :: phi,dphi
145 
146 !-- Local variables
147 real(c_double) rsq,sor,sor6,sor12
148 
149 rsq = r*r ! r^2
150 sor = lj_sigma/r ! (sig/r)
151 sor6 = sor*sor*sor !
152 sor6 = sor6*sor6 ! (sig/r)^6
153 sor12= sor6*sor6 ! (sig/r)^12
154 if (r .gt. model_cutoff) then
155  ! Argument exceeds cutoff radius
156  phi = 0.0_cd
157  dphi = 0.0_cd
158 else
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
161 endif
162 
163 end subroutine calc_phi_dphi
164 
165 !-------------------------------------------------------------------------------
166 !
167 ! Compute energy and forces on particles from the positions.
168 !
169 !-------------------------------------------------------------------------------
170 #include "kim_model_compute_log_macros.fd"
171 subroutine compute_energy_forces(model_compute_handle, &
172  model_compute_arguments_handle, ierr) bind(c)
173 implicit none
174 
175 !-- Transferred variables
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
180 
181 !-- Local variables
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
186 
187 !-- KIM variables
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(:)
197 
198 kim_log_file = __file__
199 
200 ! Unpack data from KIM object
201 !
202 ierr = 0
203 call kim_model_compute_arguments_get_argument_pointer( &
204  model_compute_arguments_handle, &
205  kim_compute_argument_name_number_of_particles, n, ierr2)
206 ierr = ierr + 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, &
210  ierr2)
211 ierr = ierr + ierr2
212 call kim_model_compute_arguments_get_argument_pointer( &
213  model_compute_arguments_handle, &
214  kim_compute_argument_name_particle_contributing, n, particlecontributing, &
215  ierr2)
216 ierr = ierr + ierr2
217 call kim_model_compute_arguments_get_argument_pointer( &
218  model_compute_arguments_handle, &
219  kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
220 ierr = ierr + ierr2
221 call kim_model_compute_arguments_get_argument_pointer( &
222  model_compute_arguments_handle, &
223  kim_compute_argument_name_partial_energy, energy, ierr2)
224 ierr = ierr + 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)
228 ierr = ierr + 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)
232 ierr = ierr + 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)
236 ierr = ierr + ierr2
237 if (ierr /= 0) then
238  kim_log_message = "get data"
239  log_error()
240  return
241 endif
242 
243 ! Check to see if we have been asked to compute the forces, energyperpart,
244 ! energy and virial
245 !
246 if (associated(energy)) then
247  comp_energy = 1
248 else
249  comp_energy = 0
250 end if
251 if (associated(force)) then
252  comp_force = 1
253 else
254  comp_force = 0
255 end if
256 if (associated(enepot)) then
257  comp_enepot = 1
258 else
259  comp_enepot = 0
260 end if
261 if (associated(virial)) then
262  comp_virial = 1
263 else
264  comp_virial = 0
265 end if
266 
267 ! Check to be sure that the species are correct
268 !
269 ierr = 1 ! assume an error
270 do i = 1,n
271  if (particlespeciescodes(i).ne.speccode) then
272  kim_log_message = "Unexpected species code detected"
273  log_error()
274  return
275  endif
276 enddo
277 ierr = 0 ! everything is ok
278 
279 ! Initialize potential energies, forces, virial term
280 !
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
285 
286 !
287 ! Compute energy and forces
288 !
289 
290 ! Loop over particles and compute energy and forces
291 !
292 do i = 1, n
293  if (particlecontributing(i) == 1) then
294  ! Set up neighbor list for next particle
295  call kim_model_compute_arguments_get_neighbor_list( &
296  model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
297  if (ierr /= 0) then
298  ! some sort of problem, exit
299  kim_log_message = "GetNeighborList failed"
300  log_error()
301  ierr = 1
302  return
303  endif
304 
305  ! Loop over the neighbors of particle i
306  !
307  do jj = 1, numnei
308 
309  j = nei1part(jj) ! get neighbor ID
310 
311  ! compute relative position vector
312  !
313  rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
314 
315  ! compute energy and forces
316  !
317  rsqij = dot_product(rij,rij) ! compute square distance
318  if ( rsqij .lt. model_cutsq ) then ! particles are interacting?
319 
320  r = sqrt(rsqij) ! compute distance
321  if (comp_force.eq.1.or.comp_virial.eq.1) then
322  call calc_phi_dphi(r,phi,dphi) ! compute pair potential
323  ! and it derivative
324  deidr = 0.5_cd*dphi
325  else
326  call calc_phi(r,phi) ! compute just pair potential
327  endif
328 
329  ! contribution to energy
330  !
331  if (comp_enepot.eq.1) then
332  enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
333  endif
334  if (comp_energy.eq.1) then
335  energy = energy + 0.5_cd*phi
336  endif
337 
338  ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
339  !
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
347  endif
348 
349  ! contribution to forces
350  !
351  if (comp_force.eq.1) then
352  force(:,i) = force(:,i) + deidr*rij/r ! accumulate force on i
353  force(:,j) = force(:,j) - deidr*rij/r ! accumulate force on j
354  endif
355 
356  endif
357 
358  enddo ! loop on jj
359 
360  endif ! if particleContributing
361 
362 enddo ! do i
363 
364 ! Everything is great
365 !
366 ierr = 0
367 return
368 
369 end subroutine compute_energy_forces
370 
371 !-------------------------------------------------------------------------------
372 !
373 ! Model destroy routine (REQUIRED)
374 !
375 !-------------------------------------------------------------------------------
376 #include "kim_model_destroy_log_macros.fd"
377 subroutine model_destroy_func(model_destroy_handle, ierr) bind(c)
378  use, intrinsic :: iso_c_binding
379  implicit none
380 
381  !-- Transferred variables
382  type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle
383  integer(c_int), intent(out) :: ierr
384 
385  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
386 
387  kim_log_file = __file__
388 
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"
392  log_information()
393  deallocate(buf)
394  ierr = 0 ! everything is good
395 end subroutine model_destroy_func
396 
397 !-------------------------------------------------------------------------------
398 !
399 ! Model refresh routine (REQUIRED)
400 !
401 !-------------------------------------------------------------------------------
402 #include "kim_model_refresh_log_macros.fd"
403 subroutine model_refresh_func(model_refresh_handle, ierr) bind(c)
404  use, intrinsic :: iso_c_binding
405  implicit none
406 
407  !-- Transferred variables
408  type(kim_model_refresh_handle_type), intent(inout) :: model_refresh_handle
409  integer(c_int), intent(out) :: ierr
410 
411  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
412 
413  kim_log_file = __file__
414 
415  call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
416  call c_f_pointer(pbuf, buf)
417 
418  kim_log_message = "Resettings influence distance and cutoffs"
419  log_information()
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, &
424  buf%half_list_hints)
425 
426  ierr = 0 ! everything is good
427 end subroutine model_refresh_func
428 
429 !-------------------------------------------------------------------------------
430 !
431 ! Model compute arguments create routine (REQUIRED)
432 !
433 !-------------------------------------------------------------------------------
434 #include "kim_model_compute_arguments_create_log_macros.fd"
435 subroutine model_compute_arguments_create(model_compute_handle, &
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
440  implicit none
441 
442  !-- Transferred variables
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
447 
448  integer(c_int) :: ierr2
449 
450  ierr = 0
451  ierr2 = 0
452 
453  ! register arguments
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)
458  ierr = ierr + 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)
463  ierr = ierr + 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)
468  ierr = ierr + 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)
473  ierr = ierr + ierr2
474 
475  ! register call backs
476  ! NONE
477 
478  if (ierr /= 0) then
479  ierr = 1
480  kim_log_message = "Unable to successfully create compute_arguments object"
481  log_error()
482  endif
483 
484  return
485 end subroutine model_compute_arguments_create
486 
487 !-------------------------------------------------------------------------------
488 !
489 ! Model compute arguments destroy routine (REQUIRED)
490 !
491 !-------------------------------------------------------------------------------
492 #include "kim_model_compute_arguments_destroy_log_macros.fd"
493 subroutine model_compute_arguments_destroy(model_compute_handle, &
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
498  implicit none
499 
500  !-- Transferred variables
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
505 
506  integer(c_int) :: ierr2
507 
508  ierr = 0
509  ierr2 = 0
510 
511  ! nothing to do
512 
513  return
514 end subroutine model_compute_arguments_destroy
515 
516 end module ex_model_ar_p_mlj_f03
517 
518 !-------------------------------------------------------------------------------
519 !
520 ! Model create routine (REQUIRED)
521 !
522 !-------------------------------------------------------------------------------
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
530 implicit none
531 
532 !-- Transferred variables
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
540 
541 !-- KIM variables
542 integer(c_int) :: ierr2
543 type(buffer_type), pointer :: buf
544 
545 kim_log_file = __file__
546 
547 ierr = 0
548 ierr2 = 0
549 
550 ! set units
551 call kim_model_create_set_units(model_create_handle, &
552  kim_length_unit_a, &
553  kim_energy_unit_ev, &
554  kim_charge_unit_unused, &
555  kim_temperature_unit_unused, &
556  kim_time_unit_unused, &
557  ierr2)
558 ierr = ierr + ierr2
559 
560 ! register species
561 call kim_model_create_set_species_code(model_create_handle, &
562  kim_species_name_ar, speccode, ierr2)
563 ierr = ierr + ierr2
564 
565 ! register numbering
566 call kim_model_create_set_model_numbering(model_create_handle, &
567  kim_numbering_one_based, ierr2);
568 ierr = ierr + ierr2
569 
570 ! register function pointers
571 call kim_model_create_set_compute_pointer(model_create_handle, &
572  kim_language_name_fortran, c_funloc(compute_energy_forces), ierr2)
573 ierr = ierr + ierr2
574 call kim_model_create_set_compute_arguments_create_pointer( &
575  model_create_handle, kim_language_name_fortran, &
576  c_funloc(model_compute_arguments_create), ierr2)
577 ierr = ierr + ierr2
578 call kim_model_create_set_compute_arguments_destroy_pointer( &
579  model_create_handle, kim_language_name_fortran, &
580  c_funloc(model_compute_arguments_destroy), ierr2)
581 ierr = ierr + ierr2
582 call kim_model_create_set_destroy_pointer(model_create_handle, &
583  kim_language_name_fortran, c_funloc(model_destroy_func), ierr2)
584 ierr = ierr + ierr2
585 call kim_model_create_set_refresh_pointer( &
586  model_create_handle, kim_language_name_fortran, &
587  c_funloc(model_refresh_func), ierr2)
588 ierr = ierr + ierr2
589 
590 ! allocate buffer
591 allocate( buf )
592 
593 ! store model buffer in KIM object
594 call kim_model_create_set_model_buffer_pointer(model_create_handle, &
595  c_loc(buf))
596 
597 ! set buffer values
598 buf%influence_distance = model_cutoff
599 buf%cutoff = model_cutoff
600 buf%padding_neighbor_hints = 1
601 buf%half_list_hints = 0
602 
603 ! register influence distance
604 call kim_model_create_set_influence_distance_pointer( &
605  model_create_handle, buf%influence_distance)
606 
607 ! register cutoff
608 call kim_model_create_set_neighbor_list_pointers(model_create_handle, &
609  1, buf%cutoff, buf%padding_neighbor_hints, buf%half_list_hints)
610 
611 if (ierr /= 0) then
612  ierr = 1
613  deallocate( buf )
614  kim_log_message = "Unable to successfully initialize model"
615  log_error()
616 endif
617 
618 return
619 
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)