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 !integer(c_int) function Compute_Energy_Forces(pkim) bind(c)
171 !implicit none
172 !
173 !!-- Transferred variables
174 !type(c_ptr), intent(in) :: pkim
175 !
176 !!-- Local variables
177 !real(c_double) :: Rij(DIM)
178 !real(c_double) :: r,Rsqij,phi,dphi,dEidr = 0.0_cd
179 !integer(c_int) :: i,j,jj,numnei,part_ret,comp_force,comp_enepot,comp_virial, &
180 ! comp_energy
181 !character (len=80) :: error_message
182 !
183 !!-- KIM variables
184 !integer(c_int), pointer :: N; type(c_ptr) :: pN
185 !real(c_double), pointer :: energy; type(c_ptr) :: penergy
186 !real(c_double), pointer :: coor(:,:); type(c_ptr) :: pcoor
187 !real(c_double), pointer :: force(:,:); type(c_ptr) :: pforce
188 !real(c_double), pointer :: enepot(:); type(c_ptr) :: penepot
189 !real(c_double), pointer :: Rij_list(:,:); type(c_ptr) :: pRij_list
190 !integer(c_int), pointer :: nei1part(:); type(c_ptr) :: pnei1part
191 !integer(c_int), pointer :: particleSpecies(:);type(c_ptr) :: pparticleSpecies
192 !real(c_double), pointer :: virial(:); type(c_ptr) :: pvirial
193 !integer(c_int) idum
194 !
195 !
196 !! Check to see if we have been asked to compute the forces, energyperpart,
197 !! energy and virial
198 !!
199 !call kim_api_getm_compute(pkim, Compute_Energy_Forces, &
200 ! "energy", comp_energy, 1, &
201 ! "forces", comp_force, 1, &
202 ! "particleEnergy", comp_enepot, 1, &
203 ! "virial", comp_virial, 1)
204 !if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
205 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
206 ! "kim_api_getm_compute", Compute_Energy_Forces)
207 ! return
208 !endif
209 !
210 !! Unpack data from KIM object
211 !!
212 !call kim_api_getm_data(pkim, Compute_Energy_Forces, &
213 ! "numberOfParticles", pN, 1, &
214 ! "particleSpecies", pparticleSpecies,1, &
215 ! "coordinates", pcoor, 1, &
216 ! "energy", penergy, TRUEFALSE(comp_energy.eq.1), &
217 ! "forces", pforce, TRUEFALSE(comp_force.eq.1), &
218 ! "particleEnergy", penepot, TRUEFALSE(comp_enepot.eq.1), &
219 ! "virial", pvirial, TRUEFALSE(comp_virial.eq.1))
220 !if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
221 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
222 ! "kim_api_getm_data_f", Compute_Energy_Forces)
223 ! return
224 !endif
225 !
226 !call c_f_pointer(pN, N)
227 !call c_f_pointer(pparticleSpecies, particleSpecies, [N])
228 !call c_f_pointer(pcoor, coor, [DIM,N])
229 !if (comp_energy.eq.1) call c_f_pointer(penergy, energy)
230 !if (comp_force.eq.1) call c_f_pointer(pforce, force, [DIM,N])
231 !if (comp_enepot.eq.1) call c_f_pointer(penepot, enepot, [N])
232 !if (comp_virial.eq.1) call c_f_pointer(pvirial, virial, [6])
233 !
234 !
235 !! Check to be sure that the species are correct
236 !!
237 !Compute_Energy_Forces = KIM_STATUS_FAIL ! assume an error
238 !do i = 1,N
239 ! if (particleSpecies(i).ne.speccode) then
240 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
241 ! "Unexpected species detected", &
242 ! Compute_Energy_Forces)
243 ! return
244 ! endif
245 !enddo
246 !Compute_Energy_Forces = KIM_STATUS_OK ! everything is ok
247 !
248 !! Initialize potential energies, forces, virial term
249 !!
250 !if (comp_enepot.eq.1) enepot = 0.0_cd
251 !if (comp_energy.eq.1) energy = 0.0_cd
252 !if (comp_force.eq.1) force = 0.0_cd
253 !if (comp_virial.eq.1) virial = 0.0_cd
254 !
255 !
256 !!
257 !! Compute energy and forces
258 !!
259 !
260 !! Loop over particles and compute energy and forces
261 !!
262 !do i=1,N
263 ! Compute_Energy_Forces = kim_api_get_neigh(pkim,1,i,part_ret,numnei, &
264 ! pnei1part,pRij_list)
265 ! if (Compute_Energy_Forces.ne.KIM_STATUS_OK) then
266 ! ! some sort of problem, exit
267 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
268 ! "kim_api_get_neigh", &
269 ! Compute_Energy_Forces)
270 ! Compute_Energy_Forces = KIM_STATUS_FAIL
271 ! return
272 ! endif
273 !
274 ! call c_f_pointer(pnei1part, nei1part, [numnei])
275 !
276 ! ! Loop over the neighbors of particle i
277 ! !
278 ! do jj = 1, numnei
279 !
280 ! j = nei1part(jj) ! get neighbor ID
281 !
282 ! ! compute relative position vector
283 ! !
284 ! Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
285 !
286 ! ! compute energy and forces
287 ! !
288 ! Rsqij = dot_product(Rij,Rij) ! compute square distance
289 ! if ( Rsqij .lt. model_cutsq ) then ! particles are interacting?
290 !
291 ! r = sqrt(Rsqij) ! compute distance
292 ! if (comp_force.eq.1.or.comp_virial.eq.1) then
293 ! call calc_phi_dphi(r,phi,dphi) ! compute pair potential
294 ! ! and it derivative
295 ! dEidr = 0.5_cd*dphi
296 ! else
297 ! call calc_phi(r,phi) ! compute just pair potential
298 ! endif
299 !
300 ! ! contribution to energy
301 ! !
302 ! if (comp_enepot.eq.1) then
303 ! enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
304 ! endif
305 ! if (comp_energy.eq.1) then
306 ! energy = energy + 0.5_cd*phi
307 ! endif
308 !
309 ! ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
310 ! !
311 ! if (comp_virial.eq.1) then
312 ! virial(1) = virial(1) + Rij(1)*Rij(1)*dEidr/r
313 ! virial(2) = virial(2) + Rij(2)*Rij(2)*dEidr/r
314 ! virial(3) = virial(3) + Rij(3)*Rij(3)*dEidr/r
315 ! virial(4) = virial(4) + Rij(2)*Rij(3)*dEidr/r
316 ! virial(5) = virial(5) + Rij(1)*Rij(3)*dEidr/r
317 ! virial(6) = virial(6) + Rij(1)*Rij(2)*dEidr/r
318 ! endif
319 !
320 ! ! contribution to forces
321 ! !
322 ! if (comp_force.eq.1) then
323 ! force(:,i) = force(:,i) + dEidr*Rij/r ! accumulate force on i
324 ! force(:,j) = force(:,j) - dEidr*Rij/r ! accumulate force on j
325 ! endif
326 !
327 ! endif
328 !
329 ! enddo ! loop on jj
330 !
331 !enddo
332 !
333 !! Everything is great
334 !!
335 !Compute_Energy_Forces = KIM_STATUS_OK
336 !return
337 !
338 !end function Compute_Energy_Forces
339 
340 !-------------------------------------------------------------------------------
341 !
342 ! Model destroy routine (REQUIRED)
343 !
344 !-------------------------------------------------------------------------------
345 #include "kim_model_destroy_log_macros.fd"
346 subroutine model_destroy_func(model_destroy_handle, ierr) bind(c)
347  use, intrinsic :: iso_c_binding
348  implicit none
349 
350  !-- Transferred variables
351  type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle
352  integer(c_int), intent(out) :: ierr
353 
354  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
355 
356  kim_log_file = __file__
357 
358  call kim_model_destroy_get_model_buffer_pointer(model_destroy_handle, pbuf)
359  call c_f_pointer(pbuf, buf)
360  kim_log_message = "deallocating model buffer"
361  log_information()
362  deallocate(buf)
363  ierr = 0 ! everything is good
364 end subroutine model_destroy_func
365 
366 !-------------------------------------------------------------------------------
367 !
368 ! Model refresh routine (REQUIRED)
369 !
370 !-------------------------------------------------------------------------------
371 #include "kim_model_refresh_log_macros.fd"
372 subroutine model_refresh_func(model_refresh_handle, ierr) bind(c)
373  use, intrinsic :: iso_c_binding
374  implicit none
375 
376  !-- Transferred variables
377  type(kim_model_refresh_handle_type), intent(inout) :: model_refresh_handle
378  integer(c_int), intent(out) :: ierr
379 
380  type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
381 
382  kim_log_file = __file__
383 
384  call kim_model_refresh_get_model_buffer_pointer(model_refresh_handle, pbuf)
385  call c_f_pointer(pbuf, buf)
386 
387  kim_log_message = "Resettings influence distance and cutoffs"
388  log_information()
389  call kim_model_refresh_set_influence_distance_pointer( &
390  model_refresh_handle, buf%cutoff(1))
391  call kim_model_refresh_set_neighbor_list_pointers( &
392  model_refresh_handle, 1, buf%cutoff, buf%padding_neighbor_hints, &
393  buf%half_list_hints)
394 
395  ierr = 0 ! everything is good
396 end subroutine model_refresh_func
397 
398 !-------------------------------------------------------------------------------
399 !
400 ! Model compute arguments create routine (REQUIRED)
401 !
402 !-------------------------------------------------------------------------------
403 #include "kim_model_compute_arguments_create_log_macros.fd"
404 subroutine model_compute_arguments_create(model_compute_handle, &
405  model_compute_arguments_create_handle, ierr) bind(c)
406  use, intrinsic :: iso_c_binding
408  log_entry=>kim_model_compute_arguments_create_log_entry
409  implicit none
410 
411  !-- Transferred variables
412  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
413  type(kim_model_compute_arguments_create_handle_type), intent(inout) :: &
414  model_compute_arguments_create_handle
415  integer(c_int), intent(out) :: ierr
416 
417  integer(c_int) :: ierr2
418 
419  ierr = 0
420  ierr2 = 0
421 
422  ! register arguments
423  call kim_model_compute_arguments_create_set_argument_support_status( &
424  model_compute_arguments_create_handle, &
425  kim_compute_argument_name_partial_energy, &
426  kim_support_status_optional, ierr2)
427  ierr = ierr + ierr2
428  call kim_model_compute_arguments_create_set_argument_support_status( &
429  model_compute_arguments_create_handle, &
430  kim_compute_argument_name_partial_forces, &
431  kim_support_status_optional, ierr2)
432  ierr = ierr + ierr2
433  call kim_model_compute_arguments_create_set_argument_support_status( &
434  model_compute_arguments_create_handle, &
435  kim_compute_argument_name_partial_particle_energy, &
436  kim_support_status_optional, ierr2)
437  ierr = ierr + ierr2
438  call kim_model_compute_arguments_create_set_argument_support_status( &
439  model_compute_arguments_create_handle, &
440  kim_compute_argument_name_partial_virial, &
441  kim_support_status_optional, ierr2)
442  ierr = ierr + ierr2
443 
444  ! register call backs
445  ! NONE
446 
447  if (ierr /= 0) then
448  ierr = 1
449  kim_log_message = "Unable to successfully create compute_arguments object"
450  log_error()
451  endif
452 
453  return
454 end subroutine model_compute_arguments_create
455 
456 !-------------------------------------------------------------------------------
457 !
458 ! Model compute arguments destroy routine (REQUIRED)
459 !
460 !-------------------------------------------------------------------------------
461 #include "kim_model_compute_arguments_destroy_log_macros.fd"
462 subroutine model_compute_arguments_destroy(model_compute_handle, &
463  model_compute_arguments_destroy_handle, ierr) bind(c)
464  use, intrinsic :: iso_c_binding
466  log_entry=>kim_model_compute_arguments_destroy_log_entry
467  implicit none
468 
469  !-- Transferred variables
470  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
471  type(kim_model_compute_arguments_destroy_handle_type), intent(inout) :: &
472  model_compute_arguments_destroy_handle
473  integer(c_int), intent(out) :: ierr
474 
475  integer(c_int) :: ierr2
476 
477  ierr = 0
478  ierr2 = 0
479 
480  ! nothing to do
481 
482  return
483 end subroutine model_compute_arguments_destroy
484 
485 end module ex_model_ar_p_mlj_f03
486 
487 !-------------------------------------------------------------------------------
488 !
489 ! Model create routine (REQUIRED)
490 !
491 !-------------------------------------------------------------------------------
492 #include "kim_model_create_log_macros.fd"
493 subroutine model_create_routine(model_create_handle, requested_length_unit, &
494  requested_energy_unit, requested_charge_unit, requested_temperature_unit, &
495  requested_time_unit, ierr) bind(c)
496 use, intrinsic :: iso_c_binding
499 implicit none
500 
501 !-- Transferred variables
502 type(kim_model_create_handle_type), intent(inout) :: model_create_handle
503 type(kim_length_unit_type), intent(in) :: requested_length_unit
504 type(kim_energy_unit_type), intent(in) :: requested_energy_unit
505 type(kim_charge_unit_type), intent(in) :: requested_charge_unit
506 type(kim_temperature_unit_type), intent(in) :: requested_temperature_unit
507 type(kim_time_unit_type), intent(in) :: requested_time_unit
508 integer(c_int), intent(out) :: ierr
509 
510 !-- KIM variables
511 integer(c_int) :: ierr2
512 type(buffer_type), pointer :: buf
513 
514 kim_log_file = __file__
515 
516 ierr = 0
517 ierr2 = 0
518 
519 ! set units
520 call kim_model_create_set_units(model_create_handle, &
521  kim_length_unit_a, &
522  kim_energy_unit_ev, &
523  kim_charge_unit_unused, &
524  kim_temperature_unit_unused, &
525  kim_time_unit_unused, &
526  ierr2)
527 ierr = ierr + ierr2
528 
529 ! register species
530 call kim_model_create_set_species_code(model_create_handle, &
531  kim_species_name_ar, speccode, ierr2)
532 ierr = ierr + ierr2
533 
534 ! register numbering
535 call kim_model_create_set_model_numbering(model_create_handle, &
536  kim_numbering_one_based, ierr2);
537 ierr = ierr + ierr2
538 
539 ! register function pointers
540 call kim_model_create_set_compute_pointer(model_create_handle, &
541  kim_language_name_fortran, c_funloc(kim_model_create_string), ierr2)
542 ierr = ierr + ierr2
543 call kim_model_create_set_compute_arguments_create_pointer( &
544  model_create_handle, kim_language_name_fortran, &
545  c_funloc(model_compute_arguments_create), ierr2)
546 ierr = ierr + ierr2
547 call kim_model_create_set_compute_arguments_destroy_pointer( &
548  model_create_handle, kim_language_name_fortran, &
549  c_funloc(model_compute_arguments_destroy), ierr2)
550 ierr = ierr + ierr2
551 call kim_model_create_set_destroy_pointer(model_create_handle, &
552  kim_language_name_fortran, c_funloc(model_destroy_func), ierr2)
553 ierr = ierr + ierr2
554 call kim_model_create_set_refresh_pointer( &
555  model_create_handle, kim_language_name_fortran, &
556  c_funloc(model_refresh_func), ierr2)
557 ierr = ierr + ierr2
558 
559 ! allocate buffer
560 allocate( buf )
561 
562 ! store model buffer in KIM object
563 call kim_model_create_set_model_buffer_pointer(model_create_handle, &
564  c_loc(buf))
565 
566 ! set buffer values
567 buf%influence_distance = model_cutoff
568 buf%cutoff = model_cutoff
569 buf%padding_neighbor_hints = 1
570 buf%half_list_hints = 0
571 
572 ! register influence distance
573 call kim_model_create_set_influence_distance_pointer( &
574  model_create_handle, buf%influence_distance)
575 
576 ! register cutoff
577 call kim_model_create_set_neighbor_list_pointers(model_create_handle, &
578  1, buf%cutoff, buf%padding_neighbor_hints, buf%half_list_hints)
579 
580 if (ierr /= 0) then
581  ierr = 1
582  deallocate( buf )
583  kim_log_message = "Unable to successfully initialize model"
584  log_error()
585 endif
586 
587 return
588 
589 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)