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 end type buffer_type
100 
101 contains
102 
103 !-------------------------------------------------------------------------------
104 !
105 ! Calculate pair potential phi(r)
106 !
107 !-------------------------------------------------------------------------------
108 subroutine calc_phi(r,phi)
109 implicit none
110 
111 !-- Transferred variables
112 real(c_double), intent(in) :: r
113 real(c_double), intent(out) :: phi
114 
115 !-- Local variables
116 real(c_double) rsq,sor,sor6,sor12
117 
118 rsq = r*r ! r^2
119 sor = lj_sigma/r ! (sig/r)
120 sor6 = sor*sor*sor !
121 sor6 = sor6*sor6 ! (sig/r)^6
122 sor12= sor6*sor6 ! (sig/r)^12
123 if (r .gt. model_cutoff) then
124  ! Argument exceeds cutoff radius
125  phi = 0.0_cd
126 else
127  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
128 endif
129 
130 end subroutine calc_phi
131 
132 !-------------------------------------------------------------------------------
133 !
134 ! Calculate pair potential phi(r) and its derivative dphi(r)
135 !
136 !-------------------------------------------------------------------------------
137 subroutine calc_phi_dphi(r,phi,dphi)
138 implicit none
139 
140 !-- Transferred variables
141 real(c_double), intent(in) :: r
142 real(c_double), intent(out) :: phi,dphi
143 
144 !-- Local variables
145 real(c_double) rsq,sor,sor6,sor12
146 
147 rsq = r*r ! r^2
148 sor = lj_sigma/r ! (sig/r)
149 sor6 = sor*sor*sor !
150 sor6 = sor6*sor6 ! (sig/r)^6
151 sor12= sor6*sor6 ! (sig/r)^12
152 if (r .gt. model_cutoff) then
153  ! Argument exceeds cutoff radius
154  phi = 0.0_cd
155  dphi = 0.0_cd
156 else
157  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
158  dphi = 24.0_cd*lj_epsilon*(-2.0_cd*sor12+sor6)/r + 2.0_cd*lj_a*r + lj_b
159 endif
160 
161 end subroutine calc_phi_dphi
162 
163 !-------------------------------------------------------------------------------
164 !
165 ! Compute energy and forces on particles from the positions.
166 !
167 !-------------------------------------------------------------------------------
168 #include "kim_model_compute_log_macros.fd"
169 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 kim_log_file = __file__
197 
198 ! Unpack data from KIM object
199 !
200 ierr = 0
201 call kim_model_compute_arguments_get_argument_pointer( &
202  model_compute_arguments_handle, &
203  kim_compute_argument_name_number_of_particles, n, ierr2)
204 ierr = ierr + ierr2
205 call kim_model_compute_arguments_get_argument_pointer( &
206  model_compute_arguments_handle, &
207  kim_compute_argument_name_particle_species_codes, n, particlespeciescodes, &
208  ierr2)
209 ierr = ierr + ierr2
210 call kim_model_compute_arguments_get_argument_pointer( &
211  model_compute_arguments_handle, &
212  kim_compute_argument_name_particle_contributing, n, particlecontributing, &
213  ierr2)
214 ierr = ierr + ierr2
215 call kim_model_compute_arguments_get_argument_pointer( &
216  model_compute_arguments_handle, &
217  kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
218 ierr = ierr + ierr2
219 call kim_model_compute_arguments_get_argument_pointer( &
220  model_compute_arguments_handle, &
221  kim_compute_argument_name_partial_energy, energy, ierr2)
222 ierr = ierr + ierr2
223 call kim_model_compute_arguments_get_argument_pointer( &
224  model_compute_arguments_handle, &
225  kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
226 ierr = ierr + ierr2
227 call kim_model_compute_arguments_get_argument_pointer( &
228  model_compute_arguments_handle, &
229  kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
230 ierr = ierr + ierr2
231 call kim_model_compute_arguments_get_argument_pointer( &
232  model_compute_arguments_handle, &
233  kim_compute_argument_name_partial_virial, 6, virial, ierr2)
234 ierr = ierr + ierr2
235 if (ierr /= 0) then
236  kim_log_message = "get data"
237  log_error()
238  return
239 endif
240 
241 ! Check to see if we have been asked to compute the forces, energyperpart,
242 ! energy and virial
243 !
244 if (associated(energy)) then
245  comp_energy = 1
246 else
247  comp_energy = 0
248 end if
249 if (associated(force)) then
250  comp_force = 1
251 else
252  comp_force = 0
253 end if
254 if (associated(enepot)) then
255  comp_enepot = 1
256 else
257  comp_enepot = 0
258 end if
259 if (associated(virial)) then
260  comp_virial = 1
261 else
262  comp_virial = 0
263 end if
264 
265 ! Check to be sure that the species are correct
266 !
267 ierr = 1 ! assume an error
268 do i = 1,n
269  if (particlespeciescodes(i).ne.speccode) then
270  kim_log_message = "Unexpected species code detected"
271  log_error()
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_model_compute_arguments_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  kim_log_message = "GetNeighborList failed"
298  log_error()
299  ierr = 1
300  return
301  endif
302 
303  ! Loop over the neighbors of particle i
304  !
305  do jj = 1, numnei
306 
307  j = nei1part(jj) ! get neighbor ID
308 
309  ! compute relative position vector
310  !
311  rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
312 
313  ! compute energy and forces
314  !
315  rsqij = dot_product(rij,rij) ! compute square distance
316  if ( rsqij .lt. model_cutsq ) then ! particles are interacting?
317 
318  r = sqrt(rsqij) ! compute distance
319  if (comp_force.eq.1.or.comp_virial.eq.1) then
320  call calc_phi_dphi(r,phi,dphi) ! compute pair potential
321  ! and it derivative
322  deidr = 0.5_cd*dphi
323  else
324  call calc_phi(r,phi) ! compute just pair potential
325  endif
326 
327  ! contribution to energy
328  !
329  if (comp_enepot.eq.1) then
330  enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
331  endif
332  if (comp_energy.eq.1) then
333  energy = energy + 0.5_cd*phi
334  endif
335 
336  ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
337  !
338  if (comp_virial.eq.1) then
339  virial(1) = virial(1) + rij(1)*rij(1)*deidr/r
340  virial(2) = virial(2) + rij(2)*rij(2)*deidr/r
341  virial(3) = virial(3) + rij(3)*rij(3)*deidr/r
342  virial(4) = virial(4) + rij(2)*rij(3)*deidr/r
343  virial(5) = virial(5) + rij(1)*rij(3)*deidr/r
344  virial(6) = virial(6) + rij(1)*rij(2)*deidr/r
345  endif
346 
347  ! contribution to forces
348  !
349  if (comp_force.eq.1) then
350  force(:,i) = force(:,i) + deidr*rij/r ! accumulate force on i
351  force(:,j) = force(:,j) - deidr*rij/r ! accumulate force on j
352  endif
353 
354  endif
355 
356  enddo ! loop on jj
357 
358  endif ! if particleContributing
359 
360 enddo ! do i
361 
362 ! Everything is great
363 !
364 ierr = 0
365 return
366 
367 end subroutine compute_energy_forces
368 
369 !-------------------------------------------------------------------------------
370 !
371 ! Model destroy routine (REQUIRED)
372 !
373 !-------------------------------------------------------------------------------
374 #include "kim_model_destroy_log_macros.fd"
375 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  kim_log_file = __file__
386 
387  call kim_model_destroy_get_model_buffer_pointer(model_destroy_handle, pbuf)
388  call c_f_pointer(pbuf, buf)
389  kim_log_message = "deallocating model buffer"
390  log_information()
391  deallocate(buf)
392  ierr = 0 ! everything is good
393 end subroutine model_destroy_func
394 
395 !-------------------------------------------------------------------------------
396 !
397 ! Model refresh routine (REQUIRED)
398 !
399 !-------------------------------------------------------------------------------
400 #include "kim_model_refresh_log_macros.fd"
401 subroutine model_refresh_func(model_refresh_handle, ierr) bind(c)
402  use, intrinsic :: iso_c_binding
403  implicit none
404 
405  !-- Transferred variables
406  type(kim_model_refresh_handle_type), intent(inout) :: model_refresh_handle
407  integer(c_int), intent(out) :: ierr
408 
409  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
410 
411  kim_log_file = __file__
412 
413  call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
414  call c_f_pointer(pbuf, buf)
415 
416  kim_log_message = "Resettings influence distance and cutoffs"
417  log_information()
418  call kim_model_refresh_set_influence_distance_pointer( &
419  model_refresh_handle, buf%cutoff(1))
420  call kim_model_refresh_set_neighbor_list_cutoffs_pointer( &
421  model_refresh_handle, 1, buf%cutoff)
422 
423  ierr = 0 ! everything is good
424 end subroutine model_refresh_func
425 
426 !-------------------------------------------------------------------------------
427 !
428 ! Model compute arguments create routine (REQUIRED)
429 !
430 !-------------------------------------------------------------------------------
431 #include "kim_model_compute_arguments_create_log_macros.fd"
432 subroutine model_compute_arguments_create(model_compute_handle, &
433  model_compute_arguments_create_handle, ierr) bind(c)
434  use, intrinsic :: iso_c_binding
436  log_entry=>kim_model_compute_arguments_create_log_entry
437  implicit none
438 
439  !-- Transferred variables
440  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
441  type(kim_model_compute_arguments_create_handle_type), intent(inout) :: &
442  model_compute_arguments_create_handle
443  integer(c_int), intent(out) :: ierr
444 
445  integer(c_int) :: ierr2
446 
447  ierr = 0
448  ierr2 = 0
449 
450  ! register arguments
451  call kim_model_compute_arguments_create_set_argument_support_status( &
452  model_compute_arguments_create_handle, &
453  kim_compute_argument_name_partial_energy, &
454  kim_support_status_optional, ierr2)
455  ierr = ierr + ierr2
456  call kim_model_compute_arguments_create_set_argument_support_status( &
457  model_compute_arguments_create_handle, &
458  kim_compute_argument_name_partial_forces, &
459  kim_support_status_optional, ierr2)
460  ierr = ierr + ierr2
461  call kim_model_compute_arguments_create_set_argument_support_status( &
462  model_compute_arguments_create_handle, &
463  kim_compute_argument_name_partial_particle_energy, &
464  kim_support_status_optional, ierr2)
465  ierr = ierr + ierr2
466  call kim_model_compute_arguments_create_set_argument_support_status( &
467  model_compute_arguments_create_handle, &
468  kim_compute_argument_name_partial_virial, &
469  kim_support_status_optional, ierr2)
470  ierr = ierr + ierr2
471 
472  ! register call backs
473  ! NONE
474 
475  if (ierr /= 0) then
476  ierr = 1
477  kim_log_message = "Unable to successfully create compute_arguments object"
478  log_error()
479  endif
480 
481  return
482 end subroutine model_compute_arguments_create
483 
484 !-------------------------------------------------------------------------------
485 !
486 ! Model compute arguments destroy routine (REQUIRED)
487 !
488 !-------------------------------------------------------------------------------
489 #include "kim_model_compute_arguments_destroy_log_macros.fd"
490 subroutine model_compute_arguments_destroy(model_compute_handle, &
491  model_compute_arguments_destroy_handle, ierr) bind(c)
492  use, intrinsic :: iso_c_binding
494  log_entry=>kim_model_compute_arguments_destroy_log_entry
495  implicit none
496 
497  !-- Transferred variables
498  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
499  type(kim_model_compute_arguments_destroy_handle_type), intent(inout) :: &
500  model_compute_arguments_destroy_handle
501  integer(c_int), intent(out) :: ierr
502 
503  integer(c_int) :: ierr2
504 
505  ierr = 0
506  ierr2 = 0
507 
508  ! nothing to do
509 
510  return
511 end subroutine model_compute_arguments_destroy
512 
513 end module ex_model_ar_p_mlj_f03
514 
515 !-------------------------------------------------------------------------------
516 !
517 ! Model create routine (REQUIRED)
518 !
519 !-------------------------------------------------------------------------------
520 #include "kim_model_create_log_macros.fd"
521 subroutine model_create_routine(model_create_handle, requested_length_unit, &
522  requested_energy_unit, requested_charge_unit, requested_temperature_unit, &
523  requested_time_unit, ierr) bind(c)
524 use, intrinsic :: iso_c_binding
527 implicit none
528 
529 !-- Transferred variables
530 type(kim_model_create_handle_type), intent(inout) :: model_create_handle
531 type(kim_length_unit_type), intent(in) :: requested_length_unit
532 type(kim_energy_unit_type), intent(in) :: requested_energy_unit
533 type(kim_charge_unit_type), intent(in) :: requested_charge_unit
534 type(kim_temperature_unit_type), intent(in) :: requested_temperature_unit
535 type(kim_time_unit_type), intent(in) :: requested_time_unit
536 integer(c_int), intent(out) :: ierr
537 
538 !-- KIM variables
539 integer(c_int) :: ierr2
540 type(buffer_type), pointer :: buf
541 
542 kim_log_file = __file__
543 
544 ierr = 0
545 ierr2 = 0
546 
547 ! set units
548 call kim_model_create_set_units(model_create_handle, &
549  kim_length_unit_a, &
550  kim_energy_unit_ev, &
551  kim_charge_unit_unused, &
552  kim_temperature_unit_unused, &
553  kim_time_unit_unused, &
554  ierr2)
555 ierr = ierr + ierr2
556 
557 ! register species
558 call kim_model_create_set_species_code(model_create_handle, &
559  kim_species_name_ar, speccode, ierr2)
560 ierr = ierr + ierr2
561 
562 ! register numbering
563 call kim_model_create_set_model_numbering(model_create_handle, &
564  kim_numbering_one_based, ierr2);
565 ierr = ierr + ierr2
566 
567 ! register function pointers
568 call kim_model_create_set_compute_pointer(model_create_handle, &
569  kim_language_name_fortran, c_funloc(compute_energy_forces), ierr2)
570 ierr = ierr + ierr2
571 call kim_model_create_set_compute_arguments_create_pointer( &
572  model_create_handle, kim_language_name_fortran, &
573  c_funloc(model_compute_arguments_create), ierr2)
574 ierr = ierr + ierr2
575 call kim_model_create_set_compute_arguments_destroy_pointer( &
576  model_create_handle, kim_language_name_fortran, &
577  c_funloc(model_compute_arguments_destroy), ierr2)
578 ierr = ierr + ierr2
579 call kim_model_create_set_destroy_pointer(model_create_handle, &
580  kim_language_name_fortran, c_funloc(model_destroy_func), ierr2)
581 ierr = ierr + ierr2
582 call kim_model_create_set_refresh_pointer( &
583  model_create_handle, kim_language_name_fortran, &
584  c_funloc(model_refresh_func), ierr2)
585 ierr = ierr + ierr2
586 
587 ! allocate buffer
588 allocate( buf )
589 
590 ! store model buffer in KIM object
591 call kim_model_create_set_model_buffer_pointer(model_create_handle, &
592  c_loc(buf))
593 
594 ! set buffer values
595 buf%influence_distance = model_cutoff
596 buf%cutoff = model_cutoff
597 
598 ! register influence distance
599 call kim_model_create_set_influence_distance_pointer( &
600  model_create_handle, buf%influence_distance)
601 
602 ! register cutoff
603 call kim_model_create_set_neighbor_list_cutoffs_pointer(model_create_handle, &
604  1, buf%cutoff)
605 
606 if (ierr /= 0) then
607  ierr = 1
608  deallocate( buf )
609  kim_log_message = "Unable to successfully initialize model"
610  log_error()
611 endif
612 
613 return
614 
615 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)