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 
45 #include "KIM_API_status.h"
46 #define THIS_FILE_NAME __FILE__
47 #define TRUEFALSE(TRUTH) merge(1,0,(TRUTH))
48 
50 
51 use, intrinsic :: iso_c_binding
52 use kim_api_f03
53 implicit none
54 
55 save
56 private
57 public compute_energy_forces, &
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 contains
95 
96 !-------------------------------------------------------------------------------
97 !
98 ! Calculate pair potential phi(r)
99 !
100 !-------------------------------------------------------------------------------
101 subroutine calc_phi(r,phi)
102 implicit none
103 
104 !-- Transferred variables
105 real(c_double), intent(in) :: r
106 real(c_double), intent(out) :: phi
107 
108 !-- Local variables
109 real(c_double) rsq,sor,sor6,sor12
110 
111 rsq = r*r ! r^2
112 sor = lj_sigma/r ! (sig/r)
113 sor6 = sor*sor*sor !
114 sor6 = sor6*sor6 ! (sig/r)^6
115 sor12= sor6*sor6 ! (sig/r)^12
116 if (r .gt. model_cutoff) then
117  ! Argument exceeds cutoff radius
118  phi = 0.0_cd
119 else
120  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
121 endif
122 
123 end subroutine calc_phi
124 
125 !-------------------------------------------------------------------------------
126 !
127 ! Calculate pair potential phi(r) and its derivative dphi(r)
128 !
129 !-------------------------------------------------------------------------------
130 subroutine calc_phi_dphi(r,phi,dphi)
131 implicit none
132 
133 !-- Transferred variables
134 real(c_double), intent(in) :: r
135 real(c_double), intent(out) :: phi,dphi
136 
137 !-- Local variables
138 real(c_double) rsq,sor,sor6,sor12
139 
140 rsq = r*r ! r^2
141 sor = lj_sigma/r ! (sig/r)
142 sor6 = sor*sor*sor !
143 sor6 = sor6*sor6 ! (sig/r)^6
144 sor12= sor6*sor6 ! (sig/r)^12
145 if (r .gt. model_cutoff) then
146  ! Argument exceeds cutoff radius
147  phi = 0.0_cd
148  dphi = 0.0_cd
149 else
150  phi = 4.0_cd*lj_epsilon*(sor12-sor6) + lj_a*rsq + lj_b*r + lj_c
151  dphi = 24.0_cd*lj_epsilon*(-2.0_cd*sor12+sor6)/r + 2.0_cd*lj_a*r + lj_b
152 endif
153 
154 end subroutine calc_phi_dphi
155 
156 !-------------------------------------------------------------------------------
157 !
158 ! Compute energy and forces on particles from the positions.
159 !
160 !-------------------------------------------------------------------------------
161 integer(c_int) function compute_energy_forces(pkim) bind(c)
162 implicit none
163 
164 !-- Transferred variables
165 type(c_ptr), intent(in) :: pkim
166 
167 !-- Local variables
168 real(c_double) :: rij(dim)
169 real(c_double) :: r,rsqij,phi,dphi,deidr = 0.0_cd
170 integer(c_int) :: i,j,jj,numnei,part_ret,comp_force,comp_enepot,comp_virial, &
171  comp_energy
172 character (len=80) :: error_message
173 
174 !-- KIM variables
175 integer(c_int), pointer :: n; type(c_ptr) :: pn
176 real(c_double), pointer :: energy; type(c_ptr) :: penergy
177 real(c_double), pointer :: coor(:,:); type(c_ptr) :: pcoor
178 real(c_double), pointer :: force(:,:); type(c_ptr) :: pforce
179 real(c_double), pointer :: enepot(:); type(c_ptr) :: penepot
180 real(c_double), pointer :: boxsidelengths(:); type(c_ptr) :: pboxsidelengths
181 real(c_double), pointer :: rij_list(:,:); type(c_ptr) :: prij_list
182 integer(c_int), pointer :: numcontrib; type(c_ptr) :: pnumcontrib
183 integer(c_int), pointer :: nei1part(:); type(c_ptr) :: pnei1part
184 integer(c_int), pointer :: particlespecies(:);type(c_ptr) :: pparticlespecies
185 real(c_double), pointer :: virial(:); type(c_ptr) :: pvirial
186 character(len=KIM_KEY_STRING_LENGTH) nbc_method
187 integer(c_int) iterorloca
188 integer(c_int) halforfull
189 integer(c_int) nbc
190 integer(c_int) numbercontrib
191 integer(c_int) idum
192 
193 numbercontrib = 0 ! initialize
194 
195 ! Determine neighbor list boundary condition (NBC)
196 ! and half versus full mode:
197 ! *****************************
198 ! * HalfOrFull = 1 -- Half
199 ! * = 2 -- Full
200 ! *****************************
201 !
202 !
203 compute_energy_forces = kim_api_get_nbc_method(pkim, nbc_method)
204 if (compute_energy_forces.lt.kim_status_ok) then
205  idum = kim_api_report_error(__line__, this_file_name, &
206  "kim_api_get_nbc_method", &
208  return
209 endif
210 if (index(nbc_method,"NEIGH_RVEC_H").eq.1) then
211  nbc = 0
212  halforfull = 1
213 elseif (index(nbc_method,"NEIGH_RVEC_F").eq.1) then
214  nbc = 0
215  halforfull = 2
216 elseif (index(nbc_method,"NEIGH_PURE_H").eq.1) then
217  nbc = 1
218  halforfull = 1
219 elseif (index(nbc_method,"NEIGH_PURE_F").eq.1) then
220  nbc = 1
221  halforfull = 2
222 elseif (index(nbc_method,"MI_OPBC_H").eq.1) then
223  nbc = 2
224  halforfull = 1
225 elseif (index(nbc_method,"MI_OPBC_F").eq.1) then
226  nbc = 2
227  halforfull = 2
228 elseif (index(nbc_method,"CLUSTER").eq.1) then
229  nbc = 3
230  halforfull = 1
231 else
232  compute_energy_forces = kim_status_fail
233  idum = kim_api_report_error(__line__, this_file_name, &
234  "Unknown NBC method", compute_energy_forces)
235  return
236 endif
237 
238 ! Determine neighbor list handling mode
239 !
240 if (nbc.ne.3) then
241  !*****************************
242  !* IterOrLoca = 1 -- Iterator
243  !* = 2 -- Locator
244  !*****************************
245  iterorloca = kim_api_get_neigh_mode(pkim, compute_energy_forces)
246  if (compute_energy_forces.lt.kim_status_ok) then
247  idum = kim_api_report_error(__line__, this_file_name, &
248  "kim_api_get_neigh_mode", &
250  return
251  endif
252  if (iterorloca.ne.1 .and. iterorloca.ne.2) then
253  compute_energy_forces = kim_status_fail
254  write(error_message,'(a,i1)') &
255  'Unsupported IterOrLoca mode = ',iterorloca
256  idum = kim_api_report_error(__line__, this_file_name, &
257  error_message, compute_energy_forces)
258  return
259  endif
260 else
261  iterorloca = 2 ! for CLUSTER NBC
262 endif
263 
264 ! Check to see if we have been asked to compute the forces, energyperpart,
265 ! energy and virial
266 !
267 call kim_api_getm_compute(pkim, compute_energy_forces, &
268  "energy", comp_energy, 1, &
269  "forces", comp_force, 1, &
270  "particleEnergy", comp_enepot, 1, &
271  "virial", comp_virial, 1)
272 if (compute_energy_forces.lt.kim_status_ok) then
273  idum = kim_api_report_error(__line__, this_file_name, &
274  "kim_api_getm_compute", compute_energy_forces)
275  return
276 endif
277 
278 ! Unpack data from KIM object
279 !
280 call kim_api_getm_data(pkim, compute_energy_forces, &
281  "numberOfParticles", pn, 1, &
282  "particleSpecies", pparticlespecies,1, &
283  "coordinates", pcoor, 1, &
284  "numberContributingParticles", pnumcontrib, truefalse(halforfull.eq.1), &
285  "boxSideLengths", pboxsidelengths, truefalse(nbc.eq.2), &
286  "energy", penergy, truefalse(comp_energy.eq.1), &
287  "forces", pforce, truefalse(comp_force.eq.1), &
288  "particleEnergy", penepot, truefalse(comp_enepot.eq.1), &
289  "virial", pvirial, truefalse(comp_virial.eq.1))
290 if (compute_energy_forces.lt.kim_status_ok) then
291  idum = kim_api_report_error(__line__, this_file_name, &
292  "kim_api_getm_data_f", compute_energy_forces)
293  return
294 endif
295 
296 call c_f_pointer(pn, n)
297 call c_f_pointer(pparticlespecies, particlespecies, [n])
298 call c_f_pointer(pcoor, coor, [dim,n])
299 if (halforfull.eq.1) call c_f_pointer(pnumcontrib, numcontrib)
300 if (nbc.eq.2) call c_f_pointer(pboxsidelengths, boxsidelengths, [dim])
301 if (comp_energy.eq.1) call c_f_pointer(penergy, energy)
302 if (comp_force.eq.1) call c_f_pointer(pforce, force, [dim,n])
303 if (comp_enepot.eq.1) call c_f_pointer(penepot, enepot, [n])
304 if (comp_virial.eq.1) call c_f_pointer(pvirial, virial, [6])
305 
306 if (halforfull.eq.1) then
307  if (nbc.ne.3) then ! non-CLUSTER cases
308  numbercontrib = numcontrib
309  else ! CLUSTER case
310  numbercontrib = n
311  endif
312 endif
313 
314 ! Check to be sure that the species are correct
315 !
316 compute_energy_forces = kim_status_fail ! assume an error
317 do i = 1,n
318  if (particlespecies(i).ne.speccode) then
319  idum = kim_api_report_error(__line__, this_file_name, &
320  "Unexpected species detected", &
322  return
323  endif
324 enddo
325 compute_energy_forces = kim_status_ok ! everything is ok
326 
327 ! Initialize potential energies, forces, virial term
328 !
329 if (comp_enepot.eq.1) enepot = 0.0_cd
330 if (comp_energy.eq.1) energy = 0.0_cd
331 if (comp_force.eq.1) force = 0.0_cd
332 if (comp_virial.eq.1) virial = 0.0_cd
333 
334 ! Initialize neighbor handling for CLUSTER NBC
335 !
336 if (nbc.eq.3) then
337  allocate( nei1part(n) )
338 endif
339 
340 ! Initialize neighbor handling for Iterator mode
341 !
342 if (iterorloca.eq.1) then
343  compute_energy_forces = kim_api_get_neigh(pkim,0,0,part_ret,numnei, &
344  pnei1part,prij_list)
345  ! check for successful initialization
346  if (compute_energy_forces.ne.kim_status_neigh_iter_init_ok) then
347  idum = kim_api_report_error(__line__, this_file_name, &
348  "kim_api_get_neigh", compute_energy_forces)
349  compute_energy_forces = kim_status_fail
350  return
351  endif
352 endif
353 
354 !
355 ! Compute energy and forces
356 !
357 
358 ! Loop over particles and compute energy and forces
359 !
360 i = 0
361 do
362 
363  ! Set up neighbor list for next particle for all NBC methods
364  !
365  if (iterorloca.eq.1) then ! ITERATOR mode
366  compute_energy_forces = kim_api_get_neigh(pkim,0,1,part_ret,numnei, &
367  pnei1part,prij_list)
368  if (compute_energy_forces.eq.kim_status_neigh_iter_past_end) exit
369  ! incremented past the end of the list,
370  ! terminate loop
371  if (compute_energy_forces.lt.kim_status_ok) then
372  ! some sort of problem, exit
373  idum = kim_api_report_error(__line__, this_file_name, &
374  "kim_api_get_neigh", &
376  return
377  endif
378 
379  i = part_ret
380 
381  else ! LOCATOR mode
382  i = i + 1
383  if (i.gt.n) exit ! incremented past end of list,
384  ! terminate loop
385  if (nbc.eq.3) then ! CLUSTER NBC method
386  numnei = n - i ! number of neighbors in list i+1, ..., N
387  nei1part(1:numnei) = (/ (i+jj, jj = 1,numnei) /)
388  compute_energy_forces = kim_status_ok
389  else
390  compute_energy_forces = kim_api_get_neigh(pkim,1,i,part_ret,numnei, &
391  pnei1part,prij_list)
392  if (compute_energy_forces.ne.kim_status_ok) then
393  ! some sort of problem, exit
394  idum = kim_api_report_error(__line__, this_file_name, &
395  "kim_api_get_neigh", &
397  compute_energy_forces = kim_status_fail
398  return
399  endif
400  endif
401  endif
402 
403  if (nbc.ne.3) call c_f_pointer(pnei1part, nei1part, [numnei])
404  if (nbc.eq.0) call c_f_pointer(prij_list, rij_list, [dim,numnei])
405 
406  ! Loop over the neighbors of particle i
407  !
408  do jj = 1, numnei
409 
410  j = nei1part(jj) ! get neighbor ID
411 
412  ! compute relative position vector
413  !
414  if (nbc.ne.0) then ! all methods except NEIGH_RVEC
415  rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j
416  else
417  rij(:) = rij_list(:,jj)
418  endif
419 
420  ! apply periodic boundary conditions if required
421  !
422  if (nbc.eq.2) then
423  where ( abs(rij) .gt. 0.5_cd*boxsidelengths )
424  ! periodic boundary conditions
425  rij = rij - sign(boxsidelengths,rij) ! applied where needed.
426  end where
427  endif
428 
429  ! compute energy and forces
430  !
431  rsqij = dot_product(rij,rij) ! compute square distance
432  if ( rsqij .lt. model_cutsq ) then ! particles are interacting?
433 
434  r = sqrt(rsqij) ! compute distance
435  if (comp_force.eq.1.or.comp_virial.eq.1) then
436  call calc_phi_dphi(r,phi,dphi) ! compute pair potential
437  ! and it derivative
438  if ((halforfull.eq.1) .and. &
439  (j .le. numbercontrib)) then ! HALF mode
440  deidr = dphi ! double contribution
441  else ! FULL mode
442  deidr = 0.5_cd*dphi ! regular contribution
443  endif
444  else
445  call calc_phi(r,phi) ! compute just pair potential
446  endif
447 
448  ! contribution to energy
449  !
450  if (comp_enepot.eq.1) then
451  enepot(i) = enepot(i) + 0.5_cd*phi ! accumulate energy
452  if ((halforfull.eq.1) .and. &
453  (j .le. numbercontrib)) & ! HALF mode
454  enepot(j) = enepot(j) + 0.5_cd*phi! (i and j share it)
455  endif
456  if (comp_energy.eq.1) then
457  if ((halforfull.eq.1) .and. &
458  (j .le. numbercontrib)) then ! HALF mode
459  energy = energy + phi ! add v to total energy
460  else ! FULL mode
461  energy = energy + 0.5_cd*phi ! add half v to total energy
462  endif
463  endif
464 
465  ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
466  !
467  if (comp_virial.eq.1) then
468  virial(1) = virial(1) + rij(1)*rij(1)*deidr/r
469  virial(2) = virial(2) + rij(2)*rij(2)*deidr/r
470  virial(3) = virial(3) + rij(3)*rij(3)*deidr/r
471  virial(4) = virial(4) + rij(2)*rij(3)*deidr/r
472  virial(5) = virial(5) + rij(1)*rij(3)*deidr/r
473  virial(6) = virial(6) + rij(1)*rij(2)*deidr/r
474  endif
475 
476  ! contribution to forces
477  !
478  if (comp_force.eq.1) then
479  force(:,i) = force(:,i) + deidr*rij/r ! accumulate force on i
480  force(:,j) = force(:,j) - deidr*rij/r ! accumulate force on j
481  endif
482 
483  endif
484 
485  enddo ! loop on jj
486 
487 enddo ! infinite do loop (terminated by exit statements above)
488 
489 ! Free temporary storage
490 !
491 if (nbc.eq.3) deallocate( nei1part )
492 
493 ! Everything is great
494 !
495 compute_energy_forces = kim_status_ok
496 return
497 
498 end function compute_energy_forces
499 
500 end module ex_model_ar_p_mlj_f03
501 
502 !-------------------------------------------------------------------------------
503 !
504 ! Model initialization routine (REQUIRED)
505 !
506 !-------------------------------------------------------------------------------
507 integer(c_int) function model_init(pkim) bind(c)
508 use, intrinsic :: iso_c_binding
510 use kim_api_f03
511 implicit none
512 
513 !-- Transferred variables
514 type(c_ptr), intent(in) :: pkim
515 
516 !-- Local variables
517 integer(c_int), parameter :: one=1
518 integer(c_int) ier, idum
519 
520 !-- KIM variables
521 real(c_double), pointer :: cutoff; type(c_ptr) :: pcutoff
522 
523 ! store pointer to compute function in KIM object
524 ier = kim_api_set_method(pkim,"compute",one,c_funloc(compute_energy_forces))
525 if (ier.lt.kim_status_ok) then
526  idum = kim_api_report_error(__line__, this_file_name, &
527  "kim_api_set_method", ier)
528  goto 42
529 endif
530 
531 ! store model cutoff in KIM object
532 pcutoff = kim_api_get_data(pkim,"cutoff",ier)
533 if (ier.lt.kim_status_ok) then
534  idum = kim_api_report_error(__line__, this_file_name, &
535  "kim_api_get_data", ier)
536  goto 42
537 endif
538 call c_f_pointer(pcutoff, cutoff)
539 cutoff = model_cutoff
540 
541 ier = kim_status_ok
542 42 continue
543 model_init = ier
544 return
545 
546 end function model_init
real(c_double), parameter, public model_cutoff
integer(c_int) function, public compute_energy_forces(pkim)