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
48 implicit none
49 
50 save
51 private
52 public compute_energy_forces, &
55  model_cutoff, &
56  speccode, &
57  buffer_type
58 
59 ! Below are the definitions and values of all Model parameters
60 integer(c_int), parameter :: cd = c_double ! used for literal constants
61 integer(c_int), parameter :: DIM = 3 ! dimensionality of space
62 integer(c_int), parameter :: speccode = 1 ! internal species code
63 real(c_double), parameter :: model_cutoff = 8.15_cd ! cutoff radius
64  ! in angstroms
65 real(c_double), parameter :: model_cutsq = model_cutoff**2
66 
67 !-------------------------------------------------------------------------------
68 ! Below are the definitions and values of all additional model parameters
69 !
70 ! Recall that the Fortran 2003 format for declaring parameters is as follows:
71 !
72 ! integer(c_int), parameter :: parname = value ! This defines an integer
73 ! ! parameter called `parname'
74 ! ! with a value equal to
75 ! ! `value' (a number)
76 !
77 ! real(c_double), parameter :: parname = value ! This defines a real(c_double)
78 ! ! parameter called `parname'
79 ! ! with a value equal to
80 ! ! `value' (a number)
81 !-------------------------------------------------------------------------------
82 real(c_double), parameter :: lj_epsilon = 0.0104_cd
83 real(c_double), parameter :: lj_sigma = 3.40_cd
84 real(c_double), parameter :: lj_cutnorm = model_cutoff/lj_sigma
85 real(c_double), parameter :: lj_a = 12.0_cd*lj_epsilon*(-26.0_cd &
86  + 7.0_cd*lj_cutnorm**6)/(lj_cutnorm**14 &
87  *lj_sigma**2)
88 real(c_double), parameter :: lj_b = 96.0_cd*lj_epsilon*(7.0_cd &
89  - 2.0_cd*lj_cutnorm**6)/(lj_cutnorm**13*lj_sigma)
90 real(c_double), parameter :: lj_c = 28.0_cd*lj_epsilon*(-13.0_cd &
91  + 4.0_cd*lj_cutnorm**6)/(lj_cutnorm**12)
92 
93 type, bind(c) :: buffer_type
94  real(c_double) :: cutoff(1)
95 end type buffer_type
96 
97 contains
98 
99 !-------------------------------------------------------------------------------
100 !
101 ! Calculate pair potential phi(r)
102 !
103 !-------------------------------------------------------------------------------
104 subroutine calc_phi(r,phi)
105 implicit none
106 
107 !-- Transferred variables
108 real(c_double), intent(in) :: r
109 real(c_double), intent(out) :: phi
110 
111 !-- Local variables
112 real(c_double) rsq,sor,sor6,sor12
113 
114 rsq = r*r ! r^2
115 sor = lj_sigma/r ! (sig/r)
116 sor6 = sor*sor*sor !
117 sor6 = sor6*sor6 ! (sig/r)^6
118 sor12= sor6*sor6 ! (sig/r)^12
119 if (r .gt. model_cutoff) then
120  ! Argument exceeds cutoff radius
121  phi = 0.0_cd
122 else
123  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
124 endif
125 
126 end subroutine calc_phi
127 
128 !-------------------------------------------------------------------------------
129 !
130 ! Calculate pair potential phi(r) and its derivative dphi(r)
131 !
132 !-------------------------------------------------------------------------------
133 subroutine calc_phi_dphi(r,phi,dphi)
134 implicit none
135 
136 !-- Transferred variables
137 real(c_double), intent(in) :: r
138 real(c_double), intent(out) :: phi,dphi
139 
140 !-- Local variables
141 real(c_double) rsq,sor,sor6,sor12
142 
143 rsq = r*r ! r^2
144 sor = lj_sigma/r ! (sig/r)
145 sor6 = sor*sor*sor !
146 sor6 = sor6*sor6 ! (sig/r)^6
147 sor12= sor6*sor6 ! (sig/r)^12
148 if (r .gt. model_cutoff) then
149  ! Argument exceeds cutoff radius
150  phi = 0.0_cd
151  dphi = 0.0_cd
152 else
153  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
154  dphi = 24.0_cd*lj_epsilon*(-2.0_cd*sor12+sor6)/r + 2.0_cd*lj_a*r + lj_b
155 endif
156 
157 end subroutine calc_phi_dphi
158 
159 !-------------------------------------------------------------------------------
160 !
161 ! Compute energy and forces on particles from the positions.
162 !
163 !-------------------------------------------------------------------------------
164 #include "kim_model_compute_log_macros.fd"
165 subroutine compute_energy_forces(model_compute_handle, ierr) bind(c)
169 implicit none
170 
171 !-- Transferred variables
172 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
173 integer(c_int), intent(out) :: ierr
174 
175 !-- Local variables
176 real(c_double) :: rij(dim)
177 real(c_double) :: r,rsqij,phi,dphi,deidr = 0.0_cd
178 integer(c_int) :: i,j,jj,numnei,comp_force,comp_enepot,comp_virial, comp_energy
179 integer(c_int) :: ierr2
180 
181 !-- KIM variables
182 integer(c_int), pointer :: n
183 real(c_double), pointer :: energy
184 real(c_double), pointer :: coor(:,:)
185 real(c_double), pointer :: force(:,:)
186 real(c_double), pointer :: enepot(:)
187 integer(c_int), pointer :: nei1part(:)
188 integer(c_int), pointer :: particlespeciescodes(:)
189 integer(c_int), pointer :: particlecontributing(:)
190 real(c_double), pointer :: virial(:)
191 
192 kim_log_file = __file__
193 
194 ! Unpack data from KIM object
195 !
196 ierr = 0
197 call kim_model_compute_get_argument_pointer(model_compute_handle, &
199 ierr = ierr + ierr2
200 call kim_model_compute_get_argument_pointer(model_compute_handle, &
201  kim_argument_name_particle_species_codes, n, particlespeciescodes, ierr2)
202 ierr = ierr + ierr2
203 call kim_model_compute_get_argument_pointer(model_compute_handle, &
204  kim_argument_name_particle_contributing, n, particlecontributing, ierr2)
205 ierr = ierr + ierr2
206 call kim_model_compute_get_argument_pointer(model_compute_handle, &
207  kim_argument_name_coordinates, dim, n, coor, ierr2)
208 ierr = ierr + ierr2
209 call kim_model_compute_get_argument_pointer(model_compute_handle, &
210  kim_argument_name_partial_energy, energy, ierr2)
211 ierr = ierr + ierr2
212 call kim_model_compute_get_argument_pointer(model_compute_handle, &
213  kim_argument_name_partial_forces, dim, n, force, ierr2)
214 ierr = ierr + ierr2
215 call kim_model_compute_get_argument_pointer(model_compute_handle, &
217 ierr = ierr + ierr2
218 call kim_model_compute_get_argument_pointer(model_compute_handle, &
219  kim_argument_name_partial_virial, 6, virial, ierr2)
220 ierr = ierr + ierr2
221 if (ierr /= 0) then
222  kim_log_message = "get data"
223  log_error()
224  return
225 endif
226 
227 ! Check to see if we have been asked to compute the forces, energyperpart,
228 ! energy and virial
229 !
230 if (associated(energy)) then
231  comp_energy = 1
232 else
233  comp_energy = 0
234 end if
235 if (associated(force)) then
236  comp_force = 1
237 else
238  comp_force = 0
239 end if
240 if (associated(enepot)) then
241  comp_enepot = 1
242 else
243  comp_enepot = 0
244 end if
245 if (associated(virial)) then
246  comp_virial = 1
247 else
248  comp_virial = 0
249 end if
250 
251 ! Check to be sure that the species are correct
252 !
253 ierr = 1 ! assume an error
254 do i = 1,n
255  if (particlespeciescodes(i).ne.speccode) then
256  kim_log_message = "Unexpected species code detected"
257  log_error()
258  return
259  endif
260 enddo
261 ierr = 0 ! everything is ok
262 
263 ! Initialize potential energies, forces, virial term
264 !
265 if (comp_enepot.eq.1) enepot = 0.0_cd
266 if (comp_energy.eq.1) energy = 0.0_cd
267 if (comp_force.eq.1) force = 0.0_cd
268 if (comp_virial.eq.1) virial = 0.0_cd
269 
270 !
271 ! Compute energy and forces
272 !
273 
274 ! Loop over particles and compute energy and forces
275 !
276 do i = 1, n
277  if (particlecontributing(i) == 1) then
278  ! Set up neighbor list for next particle
280  model_compute_handle, 1, i, numnei, nei1part, ierr)
281  if (ierr /= 0) then
282  ! some sort of problem, exit
283  kim_log_message = "GetNeighborList failed"
284  log_error()
285  ierr = 1
286  return
287  endif
288 
289  ! Loop over the neighbors of particle i
290  !
291  do jj = 1, numnei
292 
293  j = nei1part(jj) ! get neighbor ID
294 
295  ! compute relative position vector
296  !
297  rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
298 
299  ! compute energy and forces
300  !
301  rsqij = dot_product(rij,rij) ! compute square distance
302  if ( rsqij .lt. model_cutsq ) then ! particles are interacting?
303 
304  r = sqrt(rsqij) ! compute distance
305  if (comp_force.eq.1.or.comp_virial.eq.1) then
306  call calc_phi_dphi(r,phi,dphi) ! compute pair potential
307  ! and it derivative
308  deidr = 0.5_cd*dphi ! regular contribution
309  else
310  call calc_phi(r,phi) ! compute just pair potential
311  endif
312 
313  ! contribution to energy
314  !
315  if (comp_enepot.eq.1) then
316  enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
317  endif
318  if (comp_energy.eq.1) then
319  energy = energy + 0.5_cd*phi ! add half v to total energy
320  endif
321 
322  ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
323  !
324  if (comp_virial.eq.1) then
325  virial(1) = virial(1) + rij(1)*rij(1)*deidr/r
326  virial(2) = virial(2) + rij(2)*rij(2)*deidr/r
327  virial(3) = virial(3) + rij(3)*rij(3)*deidr/r
328  virial(4) = virial(4) + rij(2)*rij(3)*deidr/r
329  virial(5) = virial(5) + rij(1)*rij(3)*deidr/r
330  virial(6) = virial(6) + rij(1)*rij(2)*deidr/r
331  endif
332 
333  ! contribution to forces
334  !
335  if (comp_force.eq.1) then
336  force(:,i) = force(:,i) + deidr*rij/r ! accumulate force on i
337  force(:,j) = force(:,j) - deidr*rij/r ! accumulate force on j
338  endif
339 
340  endif
341 
342  enddo ! loop on jj
343 
344  endif ! if particleContributing
345 
346 enddo ! do i
347 
348 ! Everything is great
349 !
350 ierr = 0
351 return
352 
353 end subroutine compute_energy_forces
354 
355 !-------------------------------------------------------------------------------
356 !
357 ! Model destroy routine (REQUIRED)
358 !
359 !-------------------------------------------------------------------------------
360 #include "kim_model_destroy_log_macros.fd"
361 subroutine model_destroy_func(model_destroy_handle, ierr) bind(c)
362  use, intrinsic :: iso_c_binding
365  implicit none
366 
367  !-- Transferred variables
368  type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle
369  integer(c_int), intent(out) :: ierr
370 
371  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
372 
373  kim_log_file = __file__
374 
375  call kim_model_destroy_get_model_buffer_pointer(model_destroy_handle, pbuf)
376  call c_f_pointer(pbuf, buf)
377  kim_log_message = "deallocating model buffer"
378  log_information()
379  deallocate(buf)
380  ierr = 0 ! everything is good
381 end subroutine model_destroy_func
382 
383 !-------------------------------------------------------------------------------
384 !
385 ! Model refresh routine (REQUIRED)
386 !
387 !-------------------------------------------------------------------------------
388 #include "kim_model_refresh_log_macros.fd"
389 subroutine model_refresh_func(model_refresh_handle, ierr) bind(c)
390  use, intrinsic :: iso_c_binding
393  implicit none
394 
395  !-- Transferred variables
396  type(kim_model_refresh_handle_type), intent(inout) :: model_refresh_handle
397  integer(c_int), intent(out) :: ierr
398 
399  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
400 
401  kim_log_file = __file__
402 
403  call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
404  call c_f_pointer(pbuf, buf)
405 
406  kim_log_message = "Resettings influence distance and cutoffs"
407  log_information()
408  call kim_model_refresh_set_influence_distance_pointer( &
409  model_refresh_handle, buf%cutoff(1))
410  call kim_model_refresh_set_neighbor_list_cutoffs_pointer( &
411  model_refresh_handle, 1, buf%cutoff)
412 
413  ierr = 0 ! everything is good
414 end subroutine model_refresh_func
415 
416 end module ex_model_ar_p_mlj_f03
417 
418 !-------------------------------------------------------------------------------
419 !
420 ! Model create routine (REQUIRED)
421 !
422 !-------------------------------------------------------------------------------
423 #include "kim_model_create_log_macros.fd"
424 subroutine model_create_routine(model_create_handle, requested_length_unit, &
425  requested_energy_unit, requested_charge_unit, requested_temperature_unit, &
426  requested_time_unit, ierr) bind(c)
427 use, intrinsic :: iso_c_binding
437 implicit none
438 
439 !-- Transferred variables
440 type(kim_model_create_handle_type), intent(inout) :: model_create_handle
441 type(kim_length_unit_type), intent(in) :: requested_length_unit
442 type(kim_energy_unit_type), intent(in) :: requested_energy_unit
443 type(kim_charge_unit_type), intent(in) :: requested_charge_unit
444 type(kim_temperature_unit_type), intent(in) :: requested_temperature_unit
445 type(kim_time_unit_type), intent(in) :: requested_time_unit
446 integer(c_int), intent(out) :: ierr
447 
448 !-- KIM variables
449 integer(c_int) :: ierr2
450 type(buffer_type), pointer :: buf
451 
452 kim_log_file = __file__
453 
454 ierr2 = 0
455 
456 ! register numbering
457 call kim_model_create_set_model_numbering(model_create_handle, &
458  kim_numbering_one_based, ierr2);
459 ierr = ierr + ierr2
460 
461 ! register species
462 call kim_model_create_set_species_code(model_create_handle, &
464 ierr = ierr + ierr2
465 
466 ! register arguments
468  model_create_handle, kim_argument_name_partial_energy, &
470 ierr = ierr + ierr2
472  model_create_handle, kim_argument_name_partial_forces, &
474 ierr = ierr + ierr2
476  model_create_handle, kim_argument_name_partial_particle_energy, &
478 ierr = ierr + ierr2
480  model_create_handle, kim_argument_name_partial_virial, &
482 ierr = ierr + ierr2
483 
484 ! register call backs
485 ! NONE
486 
487 ! register function pointers
488 call kim_model_create_set_compute_pointer(model_create_handle, &
490 ierr = ierr + ierr2
491 call kim_model_create_set_destroy_pointer(model_create_handle, &
493 ierr = ierr + ierr2
494 call kim_model_create_set_refresh_pointer( &
495  model_create_handle, kim_language_name_fortran, &
496  c_funloc(model_refresh_func), ierr2)
497 ierr = ierr + ierr2
498 
499 ! set units
500 call kim_model_create_set_units(model_create_handle, &
501  kim_length_unit_a, &
502  kim_energy_unit_ev, &
503  kim_charge_unit_unused, &
504  kim_temperature_unit_unused, &
505  kim_time_unit_unused, &
506  ierr2)
507 ierr = ierr + ierr2
508 
509 ! store model cutoff in KIM object
510 allocate( buf )
511 buf%cutoff = model_cutoff
512 call kim_model_create_set_model_buffer_pointer(model_create_handle, &
513  c_loc(buf))
514 
515 ! register influence distance and cutoffs
517  model_create_handle, buf%cutoff(1))
518 
519 call kim_model_create_set_neighbor_list_cutoffs_pointer(model_create_handle, &
520  1, buf%cutoff)
521 
522 if (ierr /= 0) then
523  ierr = 1
524  kim_log_message = "Unable to successfully initialize model"
525  log_error()
526 endif
527 
528 return
529 
530 end subroutine model_create_routine
type(kim_numbering_type), public, protected kim_numbering_one_based
type(kim_argument_name_type), public, protected kim_argument_name_number_of_particles
character(len=4096), public kim_log_file
character(len=65536), public kim_log_message
type(kim_language_name_type), public, protected kim_language_name_fortran
type(kim_support_status_type), public, protected kim_support_status_optional
type(kim_argument_name_type), public, protected kim_argument_name_partial_forces
type(kim_argument_name_type), public, protected kim_argument_name_particle_contributing
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)
subroutine model_create_routine(model_create_handle, requested_length_unit, requested_energy_unit, requested_charge_unit, requested_temperature_unit, requested_time_unit, ierr)
type(kim_argument_name_type), public, protected kim_argument_name_particle_species_codes
subroutine, public model_destroy_func(model_destroy_handle, ierr)
type(kim_argument_name_type), public, protected kim_argument_name_partial_particle_energy
type(kim_argument_name_type), public, protected kim_argument_name_partial_virial
type(kim_argument_name_type), public, protected kim_argument_name_partial_energy
subroutine, public compute_energy_forces(model_compute_handle, ierr)
type(kim_argument_name_type), public, protected kim_argument_name_coordinates
static void calc_phi(double const *epsilon, double const *C, double const *Rzero, double const *shift, double const cutoff, double const r, double *phi)
subroutine, public model_refresh_func(model_refresh_handle, ierr)
real(c_double), parameter, public model_cutoff
integer(c_int), parameter, public speccode
type(kim_species_name_type), public, protected kim_species_name_ar