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, &
56  model_cutoff, &
57  speccode, &
58  buffer_type
59 
60 ! Below are the definitions and values of all Model parameters
61 integer(c_int), parameter :: cd = c_double ! used for literal constants
62 integer(c_int), parameter :: DIM = 3 ! dimensionality of space
63 integer(c_int), parameter :: speccode = 1 ! internal species code
64 real(c_double), parameter :: model_cutoff = 8.15_cd ! cutoff radius
65  ! in angstroms
66 real(c_double), parameter :: model_cutsq = model_cutoff**2
67 
68 !-------------------------------------------------------------------------------
69 ! Below are the definitions and values of all additional model parameters
70 !
71 ! Recall that the Fortran 2003 format for declaring parameters is as follows:
72 !
73 ! integer(c_int), parameter :: parname = value ! This defines an integer
74 ! ! parameter called `parname'
75 ! ! with a value equal to
76 ! ! `value' (a number)
77 !
78 ! real(c_double), parameter :: parname = value ! This defines a real(c_double)
79 ! ! parameter called `parname'
80 ! ! with a value equal to
81 ! ! `value' (a number)
82 !-------------------------------------------------------------------------------
83 real(c_double), parameter :: lj_epsilon = 0.0104_cd
84 real(c_double), parameter :: lj_sigma = 3.40_cd
85 real(c_double), parameter :: lj_cutnorm = model_cutoff/lj_sigma
86 real(c_double), parameter :: lj_A = 12.0_cd*lj_epsilon*(-26.0_cd &
87  + 7.0_cd*lj_cutnorm**6)/(lj_cutnorm**14 &
88  *lj_sigma**2)
89 real(c_double), parameter :: lj_B = 96.0_cd*lj_epsilon*(7.0_cd &
90  - 2.0_cd*lj_cutnorm**6)/(lj_cutnorm**13*lj_sigma)
91 real(c_double), parameter :: lj_C = 28.0_cd*lj_epsilon*(-13.0_cd &
92  + 4.0_cd*lj_cutnorm**6)/(lj_cutnorm**12)
93 
94 type, bind(c) :: buffer_type
95  real(c_double) :: influence_distance
96  real(c_double) :: cutoff(1)
97  integer(c_int) :: padding_neighbor_hints(1)
98  integer(c_int) :: half_list_hints(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 !integer(c_int) function Compute_Energy_Forces(pkim) bind(c)
169 !implicit none
170 !
171 !!-- Transferred variables
172 !type(c_ptr), intent(in) :: pkim
173 !
174 !!-- Local variables
175 !real(c_double) :: Rij(DIM)
176 !real(c_double) :: r,Rsqij,phi,dphi,dEidr = 0.0_cd
177 !integer(c_int) :: i,j,jj,numnei,part_ret,comp_force,comp_enepot,comp_virial, &
178 ! comp_energy
179 !character (len=80) :: error_message
180 !
181 !!-- KIM variables
182 !integer(c_int), pointer :: N; type(c_ptr) :: pN
183 !real(c_double), pointer :: energy; type(c_ptr) :: penergy
184 !real(c_double), pointer :: coor(:,:); type(c_ptr) :: pcoor
185 !real(c_double), pointer :: force(:,:); type(c_ptr) :: pforce
186 !real(c_double), pointer :: enepot(:); type(c_ptr) :: penepot
187 !real(c_double), pointer :: Rij_list(:,:); type(c_ptr) :: pRij_list
188 !integer(c_int), pointer :: nei1part(:); type(c_ptr) :: pnei1part
189 !integer(c_int), pointer :: particleSpecies(:);type(c_ptr) :: pparticleSpecies
190 !real(c_double), pointer :: virial(:); type(c_ptr) :: pvirial
191 !integer(c_int) idum
192 !
193 !
194 !! Check to see if we have been asked to compute the forces, energyperpart,
195 !! energy and virial
196 !!
197 !call kim_api_getm_compute(pkim, Compute_Energy_Forces, &
198 ! "energy", comp_energy, 1, &
199 ! "forces", comp_force, 1, &
200 ! "particleEnergy", comp_enepot, 1, &
201 ! "virial", comp_virial, 1)
202 !if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
203 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
204 ! "kim_api_getm_compute", Compute_Energy_Forces)
205 ! return
206 !endif
207 !
208 !! Unpack data from KIM object
209 !!
210 !call kim_api_getm_data(pkim, Compute_Energy_Forces, &
211 ! "numberOfParticles", pN, 1, &
212 ! "particleSpecies", pparticleSpecies,1, &
213 ! "coordinates", pcoor, 1, &
214 ! "energy", penergy, TRUEFALSE(comp_energy.eq.1), &
215 ! "forces", pforce, TRUEFALSE(comp_force.eq.1), &
216 ! "particleEnergy", penepot, TRUEFALSE(comp_enepot.eq.1), &
217 ! "virial", pvirial, TRUEFALSE(comp_virial.eq.1))
218 !if (Compute_Energy_Forces.lt.KIM_STATUS_OK) then
219 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
220 ! "kim_api_getm_data_f", Compute_Energy_Forces)
221 ! return
222 !endif
223 !
224 !call c_f_pointer(pN, N)
225 !call c_f_pointer(pparticleSpecies, particleSpecies, [N])
226 !call c_f_pointer(pcoor, coor, [DIM,N])
227 !if (comp_energy.eq.1) call c_f_pointer(penergy, energy)
228 !if (comp_force.eq.1) call c_f_pointer(pforce, force, [DIM,N])
229 !if (comp_enepot.eq.1) call c_f_pointer(penepot, enepot, [N])
230 !if (comp_virial.eq.1) call c_f_pointer(pvirial, virial, [6])
231 !
232 !
233 !! Check to be sure that the species are correct
234 !!
235 !Compute_Energy_Forces = KIM_STATUS_FAIL ! assume an error
236 !do i = 1,N
237 ! if (particleSpecies(i).ne.speccode) then
238 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
239 ! "Unexpected species detected", &
240 ! Compute_Energy_Forces)
241 ! return
242 ! endif
243 !enddo
244 !Compute_Energy_Forces = KIM_STATUS_OK ! everything is ok
245 !
246 !! Initialize potential energies, forces, virial term
247 !!
248 !if (comp_enepot.eq.1) enepot = 0.0_cd
249 !if (comp_energy.eq.1) energy = 0.0_cd
250 !if (comp_force.eq.1) force = 0.0_cd
251 !if (comp_virial.eq.1) virial = 0.0_cd
252 !
253 !
254 !!
255 !! Compute energy and forces
256 !!
257 !
258 !! Loop over particles and compute energy and forces
259 !!
260 !do i=1,N
261 ! Compute_Energy_Forces = kim_api_get_neigh(pkim,1,i,part_ret,numnei, &
262 ! pnei1part,pRij_list)
263 ! if (Compute_Energy_Forces.ne.KIM_STATUS_OK) then
264 ! ! some sort of problem, exit
265 ! idum = kim_api_report_error(__LINE__, THIS_FILE_NAME, &
266 ! "kim_api_get_neigh", &
267 ! Compute_Energy_Forces)
268 ! Compute_Energy_Forces = KIM_STATUS_FAIL
269 ! return
270 ! endif
271 !
272 ! call c_f_pointer(pnei1part, nei1part, [numnei])
273 !
274 ! ! Loop over the neighbors of particle i
275 ! !
276 ! do jj = 1, numnei
277 !
278 ! j = nei1part(jj) ! get neighbor ID
279 !
280 ! ! compute relative position vector
281 ! !
282 ! Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
283 !
284 ! ! compute energy and forces
285 ! !
286 ! Rsqij = dot_product(Rij,Rij) ! compute square distance
287 ! if ( Rsqij .lt. model_cutsq ) then ! particles are interacting?
288 !
289 ! r = sqrt(Rsqij) ! compute distance
290 ! if (comp_force.eq.1.or.comp_virial.eq.1) then
291 ! call calc_phi_dphi(r,phi,dphi) ! compute pair potential
292 ! ! and it derivative
293 ! dEidr = 0.5_cd*dphi
294 ! else
295 ! call calc_phi(r,phi) ! compute just pair potential
296 ! endif
297 !
298 ! ! contribution to energy
299 ! !
300 ! if (comp_enepot.eq.1) then
301 ! enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
302 ! endif
303 ! if (comp_energy.eq.1) then
304 ! energy = energy + 0.5_cd*phi
305 ! endif
306 !
307 ! ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
308 ! !
309 ! if (comp_virial.eq.1) then
310 ! virial(1) = virial(1) + Rij(1)*Rij(1)*dEidr/r
311 ! virial(2) = virial(2) + Rij(2)*Rij(2)*dEidr/r
312 ! virial(3) = virial(3) + Rij(3)*Rij(3)*dEidr/r
313 ! virial(4) = virial(4) + Rij(2)*Rij(3)*dEidr/r
314 ! virial(5) = virial(5) + Rij(1)*Rij(3)*dEidr/r
315 ! virial(6) = virial(6) + Rij(1)*Rij(2)*dEidr/r
316 ! endif
317 !
318 ! ! contribution to forces
319 ! !
320 ! if (comp_force.eq.1) then
321 ! force(:,i) = force(:,i) + dEidr*Rij/r ! accumulate force on i
322 ! force(:,j) = force(:,j) - dEidr*Rij/r ! accumulate force on j
323 ! endif
324 !
325 ! endif
326 !
327 ! enddo ! loop on jj
328 !
329 !enddo
330 !
331 !! Everything is great
332 !!
333 !Compute_Energy_Forces = KIM_STATUS_OK
334 !return
335 !
336 !end function Compute_Energy_Forces
337 
338 !-------------------------------------------------------------------------------
339 !
340 ! Model compute arguments create routine (REQUIRED)
341 !
342 !-------------------------------------------------------------------------------
343 #include "kim_model_compute_arguments_create_log_macros.fd"
344 subroutine model_compute_arguments_create(model_compute_handle, &
345  model_compute_arguments_create_handle, ierr) bind(c)
346  use, intrinsic :: iso_c_binding
348  log_entry=>kim_model_compute_arguments_create_log_entry
349  implicit none
350 
351  !-- Transferred variables
352  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
353  type(kim_model_compute_arguments_create_handle_type), intent(inout) :: &
354  model_compute_arguments_create_handle
355  integer(c_int), intent(out) :: ierr
356 
357  integer(c_int) :: ierr2
358  character(kind=c_char, len=100000) compute_arguments_string
359 
360  ierr = 0
361  ierr2 = 0
362 
363  ! register arguments
364  call kim_model_compute_arguments_create_set_argument_support_status( &
365  model_compute_arguments_create_handle, &
366  kim_compute_argument_name_partial_energy, &
367  kim_support_status_optional, ierr2)
368  ierr = ierr + ierr2
369  call kim_model_compute_arguments_create_set_argument_support_status( &
370  model_compute_arguments_create_handle, &
371  kim_compute_argument_name_partial_forces, &
372  kim_support_status_optional, ierr2)
373  ierr = ierr + ierr2
374  call kim_model_compute_arguments_create_set_argument_support_status( &
375  model_compute_arguments_create_handle, &
376  kim_compute_argument_name_partial_particle_energy, &
377  kim_support_status_optional, ierr2)
378  ierr = ierr + ierr2
379  call kim_model_compute_arguments_create_set_argument_support_status( &
380  model_compute_arguments_create_handle, &
381  kim_compute_argument_name_partial_virial, &
382  kim_support_status_optional, ierr2)
383  ierr = ierr + ierr2
384 
385  ! register call backs
386  ! NONE
387 
388  if (ierr /= 0) then
389  ierr = 1
390  kim_log_message = "Unable to successfully create compute_arguments object"
391  log_error()
392  else
394  model_compute_arguments_create_handle, compute_arguments_string)
395  print '(A)', trim(compute_arguments_string)
396  endif
397 
398  return
399 end subroutine model_compute_arguments_create
400 
401 !-------------------------------------------------------------------------------
402 !
403 ! Model compute arguments destroy routine (REQUIRED)
404 !
405 !-------------------------------------------------------------------------------
406 #include "kim_model_compute_arguments_destroy_log_macros.fd"
407 subroutine model_compute_arguments_destroy(model_compute_handle, &
408  model_compute_arguments_destroy_handle, ierr) bind(c)
409  use, intrinsic :: iso_c_binding
411  log_entry=>kim_model_compute_arguments_destroy_log_entry
412  implicit none
413 
414  !-- Transferred variables
415  type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
416  type(kim_model_compute_arguments_destroy_handle_type), intent(inout) :: &
417  model_compute_arguments_destroy_handle
418  integer(c_int), intent(out) :: ierr
419 
420  integer(c_int) :: ierr2
421 
422  ierr = 0
423  ierr2 = 0
424 
425  ! nothing to do
426 
427  return
428 end subroutine model_compute_arguments_destroy
429 
430 end module ex_model_ar_p_mlj_f03
431 
432 !-------------------------------------------------------------------------------
433 !
434 ! Model create routine (REQUIRED)
435 !
436 !-------------------------------------------------------------------------------
437 #include "kim_model_create_log_macros.fd"
438 subroutine model_create_routine(model_create_handle, requested_length_unit, &
439  requested_energy_unit, requested_charge_unit, requested_temperature_unit, &
440  requested_time_unit, ierr) bind(c)
441 use, intrinsic :: iso_c_binding
444 implicit none
445 
446 !-- Transferred variables
447 type(kim_model_create_handle_type), intent(inout) :: model_create_handle
448 type(kim_length_unit_type), intent(in) :: requested_length_unit
449 type(kim_energy_unit_type), intent(in) :: requested_energy_unit
450 type(kim_charge_unit_type), intent(in) :: requested_charge_unit
451 type(kim_temperature_unit_type), intent(in) :: requested_temperature_unit
452 type(kim_time_unit_type), intent(in) :: requested_time_unit
453 integer(c_int), intent(out) :: ierr
454 
455 !-- KIM variables
456 integer(c_int) :: ierr2
457 type(buffer_type), pointer :: buf
458 character(kind=c_char, len=100000) model_create_string
459 
460 kim_log_file = __file__
461 
462 ierr = 0
463 ierr2 = 0
464 
465 ! set units
466 call kim_model_create_set_units(model_create_handle, &
467  kim_length_unit_a, &
468  kim_energy_unit_ev, &
469  kim_charge_unit_unused, &
470  kim_temperature_unit_unused, &
471  kim_time_unit_unused, &
472  ierr2)
473 ierr = ierr + ierr2
474 
475 ! register species
476 call kim_model_create_set_species_code(model_create_handle, &
477  kim_species_name_ar, speccode, ierr2)
478 ierr = ierr + ierr2
479 
480 ! register numbering
481 call kim_model_create_set_model_numbering(model_create_handle, &
482  kim_numbering_one_based, ierr2);
483 ierr = ierr + ierr2
484 
485 ! register function pointers
486 call kim_model_create_set_compute_pointer(model_create_handle, &
487  kim_language_name_fortran, c_funloc(kim_model_create_string), ierr2)
488 ierr = ierr + ierr2
489 call kim_model_create_set_compute_arguments_create_pointer( &
490  model_create_handle, kim_language_name_fortran, &
491  c_funloc(model_compute_arguments_create), ierr2)
492 ierr = ierr + ierr2
493 call kim_model_create_set_compute_arguments_destroy_pointer( &
494  model_create_handle, kim_language_name_fortran, &
495  c_funloc(model_compute_arguments_destroy), ierr2)
496 ierr = ierr + ierr2
497 call kim_model_create_set_destroy_pointer(model_create_handle, &
498  kim_language_name_fortran, c_funloc(kim_model_create_string), ierr2)
499 ierr = ierr + ierr2
500 call kim_model_create_set_refresh_pointer( &
501  model_create_handle, kim_language_name_fortran, &
502  c_funloc(kim_model_create_string), ierr2)
503 ierr = ierr + ierr2
504 
505 ! allocate buffer
506 allocate( buf )
507 
508 ! store model buffer in KIM object
509 call kim_model_create_set_model_buffer_pointer(model_create_handle, &
510  c_loc(buf))
511 
512 ! set buffer values
513 buf%influence_distance = model_cutoff
514 buf%cutoff = model_cutoff
515 buf%padding_neighbor_hints = 1
516 buf%half_list_hints = 0
517 
518 ! register influence distance
519 call kim_model_create_set_influence_distance_pointer( &
520  model_create_handle, buf%influence_distance)
521 
522 ! register cutoff
523 call kim_model_create_set_neighbor_list_pointers(model_create_handle, &
524  1, buf%cutoff, buf%padding_neighbor_hints, buf%half_list_hints)
525 
526 if (ierr /= 0) then
527  ierr = 1
528  deallocate( buf )
529  kim_log_message = "Unable to successfully initialize model"
530  log_error()
531 endif
532 
533 call kim_model_create_string(model_create_handle, model_create_string)
534 print '(A)', trim(model_create_string)
535 stop
536 
537 return
538 
539 end subroutine model_create_routine
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)