kim-api 2.3.0+AppleClang.AppleClang.GNU
An Application Programming Interface (API) for the Knowledgebase of Interatomic Models (KIM).
Loading...
Searching...
No Matches
ex_model_Ar_P_MLJ_Fortran.f90
Go to the documentation of this file.
1!
2! KIM-API: An API for interatomic models
3! Copyright (c) 2013--2022, Regents of the University of Minnesota.
4! All rights reserved.
5!
6! Contributors:
7! Ryan S. Elliott
8! Ellad B. Tadmor
9! Valeriu Smirichinski
10! Stephen M. Whalen
11!
12! SPDX-License-Identifier: LGPL-2.1-or-later
13!
14! This library is free software; you can redistribute it and/or
15! modify it under the terms of the GNU Lesser General Public
16! License as published by the Free Software Foundation; either
17! version 2.1 of the License, or (at your option) any later version.
18!
19! This library is distributed in the hope that it will be useful,
20! but WITHOUT ANY WARRANTY; without even the implied warranty of
21! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
22! Lesser General Public License for more details.
23!
24! You should have received a copy of the GNU Lesser General Public License
25! along with this library; if not, write to the Free Software Foundation,
26! Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
27!
28!
29! Copyright (c) 2013--2022, Regents of the University of Minnesota.
30! All rights reserved.
31!
32
33!****************************************************************************
34!**
35!** MODULE ex_model_Ar_P_MLJ_F03
36!**
37!** Modified Lennard-Jones pair potential (with smooth cutoff) model for Ar
38!**
39!** Reference: Ashcroft and Mermin
40!**
41!** Language: Fortran 2003
42!**
43!****************************************************************************
44
46
47 use, intrinsic :: iso_c_binding
49 implicit none
50
51 save
52 private
53 public compute_energy_forces, &
58 speccode, &
60
61! Below are the definitions and values of all Model parameters
62 integer(c_int), parameter :: cd = c_double ! used for literal constants
63 integer(c_int), parameter :: DIM = 3 ! dimensionality of space
64 integer(c_int), parameter :: speccode = 1 ! internal species code
65 real(c_double), parameter :: model_cutoff = 8.15_cd ! cutoff radius
66 ! in angstroms
67 real(c_double), parameter :: model_cutsq = model_cutoff**2
68
69 !-----------------------------------------------------------------------------
70 ! Below are the definitions and values of all additional model parameters
71 !
72 ! Recall that the Fortran 2003 format for declaring parameters is as follows:
73 !
74 ! integer(c_int), parameter :: parname = value ! This defines an integer
75 ! ! parameter called `parname'
76 ! ! with a value equal to
77 ! ! `value' (a number)
78 !
79 ! real(c_double), parameter :: parname = value ! This defines a real(c_double)
80 ! ! parameter called `parname'
81 ! ! with a value equal to
82 ! ! `value' (a number)
83 !-----------------------------------------------------------------------------
84 real(c_double), parameter :: lj_epsilon = 0.0104_cd
85 real(c_double), parameter :: lj_sigma = 3.40_cd
86 real(c_double), parameter :: lj_cutnorm = model_cutoff / lj_sigma
87 real(c_double), parameter :: &
88 lj_a = 12.0_cd * lj_epsilon * (-26.0_cd + 7.0_cd * lj_cutnorm**6) &
89 / (lj_cutnorm**14 * lj_sigma**2)
90 real(c_double), parameter :: &
91 lj_b = 96.0_cd * lj_epsilon * (7.0_cd - 2.0_cd * lj_cutnorm**6) &
92 / (lj_cutnorm**13 * lj_sigma)
93 real(c_double), parameter :: &
94 lj_c = 28.0_cd * lj_epsilon * (-13.0_cd + 4.0_cd * lj_cutnorm**6) &
95 / (lj_cutnorm**12)
96
97 type, bind(c) :: buffer_type
98 real(c_double) :: influence_distance
99 real(c_double) :: cutoff(1)
100 integer(c_int) :: &
101 model_will_not_request_neighbors_of_noncontributing_particles(1)
102 end type buffer_type
103
104contains
105
106 !-----------------------------------------------------------------------------
107 !
108 ! Calculate pair potential phi(r)
109 !
110 !-----------------------------------------------------------------------------
111 recursive subroutine calc_phi(r, phi)
112 implicit none
113
114 !-- Transferred variables
115 real(c_double), intent(in) :: r
116 real(c_double), intent(out) :: phi
117
118 !-- Local variables
119 real(c_double) rsq, sor, sor6, sor12
120
121 rsq = r * r ! r^2
122 sor = lj_sigma / r ! (sig/r)
123 sor6 = sor * sor * sor !
124 sor6 = sor6 * sor6 ! (sig/r)^6
125 sor12 = sor6 * sor6 ! (sig/r)^12
126 if (r > model_cutoff) then
127 ! Argument exceeds cutoff radius
128 phi = 0.0_cd
129 else
130 phi = 4.0_cd * lj_epsilon * (sor12 - sor6) + lj_a * rsq + lj_b * r + lj_c
131 end if
132
133 end subroutine calc_phi
134
135 !-----------------------------------------------------------------------------
136 !
137 ! Calculate pair potential phi(r) and its derivative dphi(r)
138 !
139 !-----------------------------------------------------------------------------
140 recursive subroutine calc_phi_dphi(r, phi, dphi)
141 implicit none
142
143 !-- Transferred variables
144 real(c_double), intent(in) :: r
145 real(c_double), intent(out) :: phi, dphi
146
147 !-- Local variables
148 real(c_double) rsq, sor, sor6, sor12
149
150 rsq = r * r ! r^2
151 sor = lj_sigma / r ! (sig/r)
152 sor6 = sor * sor * sor !
153 sor6 = sor6 * sor6 ! (sig/r)^6
154 sor12 = sor6 * sor6 ! (sig/r)^12
155 if (r > model_cutoff) then
156 ! Argument exceeds cutoff radius
157 phi = 0.0_cd
158 dphi = 0.0_cd
159 else
160 phi = 4.0_cd * lj_epsilon * (sor12 - sor6) + lj_a * rsq + lj_b * r + lj_c
161 dphi = 24.0_cd * lj_epsilon * (-2.0_cd * sor12 + sor6) / r &
162 + 2.0_cd * lj_a * r + lj_b
163 end if
164
165 end subroutine calc_phi_dphi
166
167 !-----------------------------------------------------------------------------
168 !
169 ! Compute energy and forces on particles from the positions.
170 !
171 !-----------------------------------------------------------------------------
172 recursive subroutine compute_energy_forces( &
173 model_compute_handle, model_compute_arguments_handle, ierr) bind(c)
174 implicit none
175
176 !-- Transferred variables
177 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
178 type(kim_model_compute_arguments_handle_type), intent(in) :: &
179 model_compute_arguments_handle
180 integer(c_int), intent(out) :: ierr
181
182 !-- Local variables
183 real(c_double) :: rij(dim)
184 real(c_double) :: r, rsqij, phi, dphi, deidr = 0.0_cd
185 integer(c_int) :: i, j, jj, numnei, comp_force, comp_enepot, &
186 comp_virial, comp_energy
187 integer(c_int) :: ierr2
188
189 !-- KIM variables
190 integer(c_int), pointer :: n
191 real(c_double), pointer :: energy
192 real(c_double), pointer :: coor(:, :)
193 real(c_double), pointer :: force(:, :)
194 real(c_double), pointer :: enepot(:)
195 integer(c_int), pointer :: nei1part(:)
196 integer(c_int), pointer :: particlespeciescodes(:)
197 integer(c_int), pointer :: particlecontributing(:)
198 real(c_double), pointer :: virial(:)
199
200 ! Unpack data from KIM object
201 !
202 ierr = 0
203 call kim_get_argument_pointer( &
204 model_compute_arguments_handle, &
205 kim_compute_argument_name_number_of_particles, n, ierr2)
206 ierr = ierr + ierr2
207 call kim_get_argument_pointer( &
208 model_compute_arguments_handle, &
209 kim_compute_argument_name_particle_species_codes, n, &
210 particlespeciescodes, ierr2)
211 ierr = ierr + ierr2
212 call kim_get_argument_pointer( &
213 model_compute_arguments_handle, &
214 kim_compute_argument_name_particle_contributing, n, &
215 particlecontributing, ierr2)
216 ierr = ierr + ierr2
217 call kim_get_argument_pointer( &
218 model_compute_arguments_handle, &
219 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
220 ierr = ierr + ierr2
221 call kim_get_argument_pointer( &
222 model_compute_arguments_handle, &
223 kim_compute_argument_name_partial_energy, energy, ierr2)
224 ierr = ierr + ierr2
225 call kim_get_argument_pointer( &
226 model_compute_arguments_handle, &
227 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
228 ierr = ierr + ierr2
229 call kim_get_argument_pointer( &
230 model_compute_arguments_handle, &
231 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
232 ierr = ierr + ierr2
233 call kim_get_argument_pointer( &
234 model_compute_arguments_handle, &
235 kim_compute_argument_name_partial_virial, 6, virial, ierr2)
236 ierr = ierr + ierr2
237 if (ierr /= 0) then
238 call kim_log_entry(model_compute_arguments_handle, &
239 kim_log_verbosity_error, "get data")
240 ierr = 1
241 return
242 end if
243
244 ! Check to see if we have been asked to compute the forces, energyperpart,
245 ! energy and virial
246 !
247 if (associated(energy)) then
248 comp_energy = 1
249 else
250 comp_energy = 0
251 end if
252 if (associated(force)) then
253 comp_force = 1
254 else
255 comp_force = 0
256 end if
257 if (associated(enepot)) then
258 comp_enepot = 1
259 else
260 comp_enepot = 0
261 end if
262 if (associated(virial)) then
263 comp_virial = 1
264 else
265 comp_virial = 0
266 end if
267
268 ! Check to be sure that the species are correct
269 !
270 ierr = 1 ! assume an error
271 do i = 1, n
272 if (particlespeciescodes(i) /= speccode) then
273 call kim_log_entry( &
274 model_compute_handle, kim_log_verbosity_error, &
275 "Unexpected species code detected")
276 ierr = 1
277 return
278 end if
279 end do
280 ierr = 0 ! everything is ok
281
282 ! Initialize potential energies, forces, virial term
283 !
284 if (comp_enepot == 1) enepot = 0.0_cd
285 if (comp_energy == 1) energy = 0.0_cd
286 if (comp_force == 1) force = 0.0_cd
287 if (comp_virial == 1) virial = 0.0_cd
288
289 !
290 ! Compute energy and forces
291 !
292
293 ! Loop over particles and compute energy and forces
294 !
295 do i = 1, n
296 if (particlecontributing(i) == 1) then
297 ! Set up neighbor list for next particle
298 call kim_get_neighbor_list( &
299 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
300 if (ierr /= 0) then
301 ! some sort of problem, exit
302 call kim_log_entry( &
303 model_compute_arguments_handle, kim_log_verbosity_error, &
304 "GetNeighborList failed")
305 ierr = 1
306 return
307 end if
308
309 ! Loop over the neighbors of particle i
310 !
311 do jj = 1, numnei
312
313 j = nei1part(jj) ! get neighbor ID
314
315 ! compute relative position vector
316 !
317 rij(:) = coor(:, j) - coor(:, i) ! distance vector between i j
318
319 ! compute energy and forces
320 !
321 rsqij = dot_product(rij, rij) ! compute square distance
322 if (rsqij < model_cutsq) then ! particles are interacting?
323
324 r = sqrt(rsqij) ! compute distance
325 if (comp_force == 1 .or. comp_virial == 1) then
326 call calc_phi_dphi(r, phi, dphi) ! compute pair potential
327 ! and it derivative
328 deidr = 0.5_cd * dphi
329 else
330 call calc_phi(r, phi) ! compute just pair potential
331 end if
332
333 ! contribution to energy
334 !
335 if (comp_enepot == 1) then
336 enepot(i) = enepot(i) + 0.5_cd * phi ! accumulate energy
337 end if
338 if (comp_energy == 1) then
339 energy = energy + 0.5_cd * phi
340 end if
341
342 ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r
343 !
344 if (comp_virial == 1) then
345 virial(1) = virial(1) + rij(1) * rij(1) * deidr / r
346 virial(2) = virial(2) + rij(2) * rij(2) * deidr / r
347 virial(3) = virial(3) + rij(3) * rij(3) * deidr / r
348 virial(4) = virial(4) + rij(2) * rij(3) * deidr / r
349 virial(5) = virial(5) + rij(1) * rij(3) * deidr / r
350 virial(6) = virial(6) + rij(1) * rij(2) * deidr / r
351 end if
352
353 ! contribution to forces
354 !
355 if (comp_force == 1) then
356 force(:, i) = force(:, i) + deidr * rij / r ! accumulate force on i
357 force(:, j) = force(:, j) - deidr * rij / r ! accumulate force on j
358 end if
359
360 end if
361
362 end do ! loop on jj
363
364 end if ! if particleContributing
365
366 end do ! do i
367
368 ! Everything is great
369 !
370 ierr = 0
371 return
372
373 end subroutine compute_energy_forces
374
375 !-----------------------------------------------------------------------------
376 !
377 ! Model destroy routine (REQUIRED)
378 !
379 !-----------------------------------------------------------------------------
380 recursive subroutine model_destroy_func(model_destroy_handle, ierr) bind(c)
381 use, intrinsic :: iso_c_binding
382 implicit none
383
384 !-- Transferred variables
385 type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle
386 integer(c_int), intent(out) :: ierr
387
388 type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
389
390 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
391 call c_f_pointer(pbuf, buf)
392 call kim_log_entry(model_destroy_handle, &
393 kim_log_verbosity_information, &
394 "deallocating model buffer")
395 deallocate (buf)
396 ierr = 0 ! everything is good
397 end subroutine model_destroy_func
398
399 !-----------------------------------------------------------------------------
400 !
401 ! Model compute arguments create routine (REQUIRED)
402 !
403 !-----------------------------------------------------------------------------
404 recursive subroutine model_compute_arguments_create( &
405 model_compute_handle, model_compute_arguments_create_handle, ierr) bind(c)
406 use, intrinsic :: iso_c_binding
407 implicit none
408
409 !-- Transferred variables
410 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
411 type(kim_model_compute_arguments_create_handle_type), intent(inout) :: &
412 model_compute_arguments_create_handle
413 integer(c_int), intent(out) :: ierr
414
415 integer(c_int) :: ierr2
416
417 ! avoid unused dummy argument warnings
418 if (model_compute_handle == kim_model_compute_null_handle) continue
419
420 ierr = 0
421 ierr2 = 0
422
423 ! register arguments
424 call kim_set_argument_support_status( &
425 model_compute_arguments_create_handle, &
426 kim_compute_argument_name_partial_energy, &
427 kim_support_status_optional, ierr2)
428 ierr = ierr + ierr2
429 call kim_set_argument_support_status( &
430 model_compute_arguments_create_handle, &
431 kim_compute_argument_name_partial_forces, &
432 kim_support_status_optional, ierr2)
433 ierr = ierr + ierr2
434 call kim_set_argument_support_status( &
435 model_compute_arguments_create_handle, &
436 kim_compute_argument_name_partial_particle_energy, &
437 kim_support_status_optional, ierr2)
438 ierr = ierr + ierr2
439 call kim_set_argument_support_status( &
440 model_compute_arguments_create_handle, &
441 kim_compute_argument_name_partial_virial, &
442 kim_support_status_optional, ierr2)
443 ierr = ierr + ierr2
444
445 ! register call backs
446 ! NONE
447
448 if (ierr /= 0) then
449 ierr = 1
450 call kim_log_entry( &
451 model_compute_arguments_create_handle, kim_log_verbosity_error, &
452 "Unable to successfully create compute_arguments object")
453 end if
454
455 ierr = 0
456 return
457 end subroutine model_compute_arguments_create
458
459 !-----------------------------------------------------------------------------
460 !
461 ! Model compute arguments destroy routine (REQUIRED)
462 !
463 !-----------------------------------------------------------------------------
464 recursive subroutine model_compute_arguments_destroy( &
465 model_compute_handle, model_compute_arguments_destroy_handle, ierr) bind(c)
466 use, intrinsic :: iso_c_binding
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 ! avoid unused dummy argument warnings
476 if (model_compute_handle == kim_model_compute_null_handle) continue
477 if (model_compute_arguments_destroy_handle == &
478 kim_model_compute_arguments_destroy_null_handle) continue
479
480 ierr = 0
481 return
483
484end module ex_model_ar_p_mlj_f03
485
486!-------------------------------------------------------------------------------
487!
488! Model create routine (REQUIRED)
489!
490!-------------------------------------------------------------------------------
491recursive subroutine model_create_routine( &
492 model_create_handle, requested_length_unit, requested_energy_unit, &
493 requested_charge_unit, requested_temperature_unit, requested_time_unit, &
494 ierr) bind(c)
495 use, intrinsic :: iso_c_binding
498 implicit none
499
500 !-- Transferred variables
501 type(kim_model_create_handle_type), intent(inout) :: model_create_handle
502 type(kim_length_unit_type), intent(in), value :: requested_length_unit
503 type(kim_energy_unit_type), intent(in), value :: requested_energy_unit
504 type(kim_charge_unit_type), intent(in), value :: requested_charge_unit
505 type(kim_temperature_unit_type), intent(in), value :: &
506 requested_temperature_unit
507 type(kim_time_unit_type), intent(in), value :: 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 ierr = 0
515 ierr2 = 0
516
517 ! avoid unsed dummy argument warnings
518 if (requested_length_unit == kim_length_unit_unused) continue
519 if (requested_energy_unit == kim_energy_unit_unused) continue
520 if (requested_charge_unit == kim_charge_unit_unused) continue
521 if (requested_temperature_unit == kim_temperature_unit_unused) continue
522 if (requested_time_unit == kim_time_unit_unused) continue
523
524 ! set units
525 call kim_set_units(model_create_handle, &
526 kim_length_unit_a, &
527 kim_energy_unit_ev, &
528 kim_charge_unit_unused, &
529 kim_temperature_unit_unused, &
530 kim_time_unit_unused, &
531 ierr2)
532 ierr = ierr + ierr2
533
534 ! register species
535 call kim_set_species_code(model_create_handle, &
536 kim_species_name_ar, speccode, ierr2)
537 ierr = ierr + ierr2
538
539 ! register numbering
540 call kim_set_model_numbering(model_create_handle, &
541 kim_numbering_one_based, ierr2)
542 ierr = ierr + ierr2
543
544 ! register function pointers
545 call kim_set_routine_pointer( &
546 model_create_handle, &
547 kim_model_routine_name_compute, kim_language_name_fortran, &
548 1, c_funloc(compute_energy_forces), ierr2)
549 ierr = ierr + ierr2
550 call kim_set_routine_pointer( &
551 model_create_handle, kim_model_routine_name_compute_arguments_create, &
552 kim_language_name_fortran, 1, c_funloc(model_compute_arguments_create), &
553 ierr2)
554 ierr = ierr + ierr2
555 call kim_set_routine_pointer( &
556 model_create_handle, kim_model_routine_name_compute_arguments_destroy, &
557 kim_language_name_fortran, 1, c_funloc(model_compute_arguments_destroy), &
558 ierr2)
559 ierr = ierr + ierr2
560 call kim_set_routine_pointer( &
561 model_create_handle, &
562 kim_model_routine_name_destroy, kim_language_name_fortran, &
563 1, c_funloc(model_destroy_func), ierr2)
564 ierr = ierr + ierr2
565
566 ! allocate buffer
567 allocate (buf)
568
569 ! store model buffer in KIM object
570 call kim_set_model_buffer_pointer(model_create_handle, &
571 c_loc(buf))
572
573 ! set buffer values
574 buf%influence_distance = model_cutoff
575 buf%cutoff = model_cutoff
576 buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
577
578 ! register influence distance
579 call kim_set_influence_distance_pointer( &
580 model_create_handle, buf%influence_distance)
581
582 ! register cutoff
583 call kim_set_neighbor_list_pointers( &
584 model_create_handle, 1, buf%cutoff, &
585 buf%model_will_not_request_neighbors_of_noncontributing_particles)
586
587 if (ierr /= 0) then
588 ierr = 1
589 deallocate (buf)
590 call kim_log_entry(model_create_handle, kim_log_verbosity_error, &
591 "Unable to successfully initialize model")
592 end if
593
594 ierr = 0
595 return
596
597end subroutine model_create_routine
static void calc_phi(double const *epsilon, double const *C, double const *Rzero, double const *shift, double const cutoff, double const r, double *phi)
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)
recursive subroutine, public model_compute_arguments_destroy(model_compute_handle, model_compute_arguments_destroy_handle, ierr)
recursive subroutine, public model_compute_arguments_create(model_compute_handle, model_compute_arguments_create_handle, ierr)
integer(c_int), parameter, public speccode
recursive subroutine, public compute_energy_forces(model_compute_handle, model_compute_arguments_handle, ierr)
recursive subroutine, public model_destroy_func(model_destroy_handle, ierr)
real(c_double), parameter, public model_cutoff