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_driver_P_LJ.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!****************************************************************************
30!**
31!** MODULE ex_model_driver_P_LJ
32!**
33!** Lennard-Jones pair potential KIM Model Driver
34!** shifted to have zero energy at the cutoff radius
35!**
36!** Language: Fortran 2003
37!**
38!****************************************************************************
39
41
42 use, intrinsic :: iso_c_binding
44 implicit none
45
46 save
47 private
48 public &
53 refresh, &
55 destroy, &
56 calc_phi, &
60
61 ! Below are the definitions and values of all Model parameters
62 integer(c_int), parameter :: cd = c_double ! for literal constants
63 integer(c_int), parameter :: DIM = 3 ! dimensionality of space
64 integer(c_int), parameter :: speccode = 1 ! internal species code
65
66 !-----------------------------------------------------------------------------
67 !
68 ! Definition of Buffer type
69 !
70 !-----------------------------------------------------------------------------
71 type, bind(c) :: BUFFER_TYPE
72 character(c_char) :: species_name(100)
73 real(c_double) :: influence_distance(1)
74 real(c_double) :: pcutoff(1)
75 real(c_double) :: cutsq(1)
76 integer(c_int) :: &
77 model_will_not_request_neighbors_of_noncontributing_particles(1)
78 real(c_double) :: epsilon(1)
79 real(c_double) :: sigma(1)
80 real(c_double) :: shift(1)
81 end type buffer_type
82
83contains
84
85 !-----------------------------------------------------------------------------
86 !
87 ! Calculate pair potential phi(r)
88 !
89 !-----------------------------------------------------------------------------
90 recursive subroutine calc_phi(model_epsilon, &
91 model_sigma, &
92 model_shift, &
93 model_cutoff, r, phi)
94 implicit none
95
96 !-- Transferred variables
97 real(c_double), intent(in) :: model_epsilon
98 real(c_double), intent(in) :: model_sigma
99 real(c_double), intent(in) :: model_shift
100 real(c_double), intent(in) :: model_cutoff
101 real(c_double), intent(in) :: r
102 real(c_double), intent(out) :: phi
103
104 !-- Local variables
105 real(c_double) rsq, sor, sor6, sor12
106
107 rsq = r * r ! r^2
108 sor = model_sigma / r ! (sig/r)
109 sor6 = sor * sor * sor !
110 sor6 = sor6 * sor6 ! (sig/r)^6
111 sor12 = sor6 * sor6 ! (sig/r)^12
112 if (r > model_cutoff) then
113 ! Argument exceeds cutoff radius
114 phi = 0.0_cd
115 else
116 phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
117 end if
118
119 end subroutine calc_phi
120
121 !-----------------------------------------------------------------------------
122 !
123 ! Calculate pair potential phi(r) and its derivative dphi(r)
124 !
125 !-----------------------------------------------------------------------------
126 recursive subroutine calc_phi_dphi(model_epsilon, &
127 model_sigma, &
128 model_shift, &
129 model_cutoff, r, phi, dphi)
130 implicit none
131
132 !-- Transferred variables
133 real(c_double), intent(in) :: model_epsilon
134 real(c_double), intent(in) :: model_sigma
135 real(c_double), intent(in) :: model_shift
136 real(c_double), intent(in) :: model_cutoff
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 = model_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 > 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 * model_epsilon * (sor12 - sor6) + model_shift
154 dphi = 24.0_cd * model_epsilon * (-2.0_cd * sor12 + sor6) / r
155 end if
156
157 end subroutine calc_phi_dphi
158
159 !-----------------------------------------------------------------------------
160 !
161 ! Calculate pair potential phi(r) and its derivatives dphi(r) and d2phi(r)
162 !
163 !-----------------------------------------------------------------------------
164 recursive subroutine calc_phi_dphi_d2phi(model_epsilon, &
165 model_sigma, &
166 model_shift, &
167 model_cutoff, r, phi, dphi, d2phi)
168 implicit none
169
170 !-- Transferred variables
171 real(c_double), intent(in) :: model_epsilon
172 real(c_double), intent(in) :: model_sigma
173 real(c_double), intent(in) :: model_shift
174 real(c_double), intent(in) :: model_cutoff
175 real(c_double), intent(in) :: r
176 real(c_double), intent(out) :: phi, dphi, d2phi
177
178 !-- Local variables
179 real(c_double) rsq, sor, sor6, sor12
180
181 rsq = r * r ! r^2
182 sor = model_sigma / r ! (sig/r)
183 sor6 = sor * sor * sor !
184 sor6 = sor6 * sor6 ! (sig/r)^6
185 sor12 = sor6 * sor6 ! (sig/r)^12
186 if (r > model_cutoff) then
187 ! Argument exceeds cutoff radius
188 phi = 0.0_cd
189 dphi = 0.0_cd
190 d2phi = 0.0_cd
191 else
192 phi = 4.0_cd * model_epsilon * (sor12 - sor6) + model_shift
193 dphi = 24.0_cd * model_epsilon * (-2.0_cd * sor12 + sor6) / r
194 d2phi = 24.0_cd * model_epsilon * (26.0_cd * sor12 - 7.0_cd * sor6) / rsq
195 end if
196
197 end subroutine calc_phi_dphi_d2phi
198
199 !-----------------------------------------------------------------------------
200 !
201 ! Compute energy and forces on particles from the positions.
202 !
203 !-----------------------------------------------------------------------------
204 recursive subroutine compute_energy_forces( &
205 model_compute_handle, model_compute_arguments_handle, ierr) bind(c)
206 implicit none
207
208 !-- Transferred variables
209 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
210 type(kim_model_compute_arguments_handle_type), intent(in) :: &
211 model_compute_arguments_handle
212 integer(c_int), intent(out) :: ierr
213
214 !-- Local variables
215 real(c_double) :: r, rsqij, phi, dphi, d2phi, deidr, d2eidr
216 integer(c_int) :: i, j, jj, numnei
217 integer(c_int) :: ierr2
218 integer(c_int) :: comp_force, comp_energy, comp_enepot, comp_process_dedr, &
219 comp_process_d2edr2
220 type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
221
222 real(c_double), target :: rij(dim)
223 real(c_double), target :: rij_pairs(dim, 2)
224 real(c_double), target :: r_pairs(2)
225 integer(c_int), target :: i_pairs(2), j_pairs(2)
226
227 !-- KIM variables
228 real(c_double) :: model_cutoff
229 integer(c_int), pointer :: n
230 real(c_double), pointer :: energy
231 real(c_double), pointer :: coor(:, :)
232 real(c_double), pointer :: force(:, :)
233 real(c_double), pointer :: enepot(:)
234 integer(c_int), pointer :: nei1part(:)
235 integer(c_int), pointer :: particlespeciescodes(:)
236 integer(c_int), pointer :: particlecontributing(:)
237
238 ! get model buffer from KIM object
239 call kim_get_model_buffer_pointer(model_compute_handle, pbuf)
240 call c_f_pointer(pbuf, buf)
241
242 model_cutoff = buf%influence_distance(1)
243
244 ! Check to see if we have been asked to compute the forces, energyperpart,
245 ! energy and d1Edr
246 !
247 ierr = 0
248 call kim_is_callback_present( &
249 model_compute_arguments_handle, &
250 kim_compute_callback_name_process_dedr_term, comp_process_dedr, ierr2)
251 ierr = ierr + ierr2
252 call kim_is_callback_present( &
253 model_compute_arguments_handle, &
254 kim_compute_callback_name_process_d2edr2_term, comp_process_d2edr2, ierr2)
255 ierr = ierr + ierr2
256 if (ierr /= 0) then
257 call kim_log_entry(model_compute_arguments_handle, &
258 kim_log_verbosity_error, "get_compute")
259 ierr = 1
260 return
261 end if
262
263 ! Unpack data from KIM object
264 !
265 ierr = 0
266 call kim_get_argument_pointer( &
267 model_compute_arguments_handle, &
268 kim_compute_argument_name_number_of_particles, n, ierr2)
269 ierr = ierr + ierr2
270
271 call kim_get_argument_pointer( &
272 model_compute_arguments_handle, &
273 kim_compute_argument_name_particle_species_codes, &
274 n, particlespeciescodes, ierr2)
275 ierr = ierr + ierr2
276 call kim_get_argument_pointer( &
277 model_compute_arguments_handle, &
278 kim_compute_argument_name_particle_contributing, &
279 n, particlecontributing, ierr2)
280 ierr = ierr + ierr2
281 call kim_get_argument_pointer( &
282 model_compute_arguments_handle, &
283 kim_compute_argument_name_coordinates, dim, n, coor, ierr2)
284 ierr = ierr + ierr2
285 call kim_get_argument_pointer( &
286 model_compute_arguments_handle, &
287 kim_compute_argument_name_partial_energy, energy, ierr2)
288 ierr = ierr + ierr2
289 call kim_get_argument_pointer( &
290 model_compute_arguments_handle, &
291 kim_compute_argument_name_partial_forces, dim, n, force, ierr2)
292 ierr = ierr + ierr2
293 call kim_get_argument_pointer( &
294 model_compute_arguments_handle, &
295 kim_compute_argument_name_partial_particle_energy, n, enepot, ierr2)
296 ierr = ierr + ierr2
297 if (ierr /= 0) then
298 call kim_log_entry(model_compute_arguments_handle, &
299 kim_log_verbosity_error, "get_argument_pointer")
300 ierr = 1
301 return
302 end if
303
304 if (associated(energy)) then
305 comp_energy = 1
306 else
307 comp_energy = 0
308 end if
309 if (associated(force)) then
310 comp_force = 1
311 else
312 comp_force = 0
313 end if
314 if (associated(enepot)) then
315 comp_enepot = 1
316 else
317 comp_enepot = 0
318 end if
319
320 ! Check to be sure that the species are correct
321 !
322
323 ierr = 1 ! assume an error
324 do i = 1, n
325 if (particlespeciescodes(i) /= speccode) then
326 call kim_log_entry( &
327 model_compute_arguments_handle, kim_log_verbosity_error, &
328 "Unexpected species code detected")
329 ierr = 1
330 return
331 end if
332 end do
333 ierr = 0 ! everything is ok
334
335 ! Initialize potential energies, forces
336 !
337 if (comp_enepot == 1) enepot = 0.0_cd
338 if (comp_energy == 1) energy = 0.0_cd
339 if (comp_force == 1) force = 0.0_cd
340
341 !
342 ! Compute energy and forces
343 !
344
345 ! Loop over particles and compute energy and forces
346 !
347 do i = 1, n
348
349 if (particlecontributing(i) == 1) then
350 ! Set up neighbor list for next particle
351 !
352 call kim_get_neighbor_list( &
353 model_compute_arguments_handle, 1, i, numnei, nei1part, ierr)
354 if (ierr /= 0) then
355 ! some sort of problem, exit
356 call kim_log_entry( &
357 model_compute_arguments_handle, kim_log_verbosity_error, &
358 "kim_api_get_neigh")
359 ierr = 1
360 return
361 end if
362
363 ! Loop over the neighbors of particle i
364 !
365 do jj = 1, numnei
366
367 j = nei1part(jj) ! get neighbor ID
368
369 ! compute relative position vector
370 !
371 rij(:) = coor(:, j) - coor(:, i) ! distance vector between i j
372
373 ! compute energy and forces
374 !
375 rsqij = dot_product(rij, rij) ! compute square distance
376 if (rsqij < buf%cutsq(1)) then ! particles are interacting?
377
378 r = sqrt(rsqij) ! compute distance
379 if (comp_process_d2edr2 == 1) then
380 call calc_phi_dphi_d2phi(buf%epsilon(1), &
381 buf%sigma(1), &
382 buf%shift(1), &
383 buf%Pcutoff(1), &
384 r, phi, dphi, d2phi) ! compute pair potential
385 ! and it derivatives
386 deidr = 0.5_cd * dphi ! regular contribution
387 d2eidr = 0.5_cd * d2phi
388 elseif (comp_force == 1 .or. comp_process_dedr == 1) then
389 call calc_phi_dphi(buf%epsilon(1), &
390 buf%sigma(1), &
391 buf%shift(1), &
392 buf%Pcutoff(1), &
393 r, phi, dphi) ! compute pair potential
394 ! and it derivative
395
396 deidr = 0.5_cd * dphi ! regular contribution
397 else
398 call calc_phi(buf%epsilon(1), &
399 buf%sigma(1), &
400 buf%shift(1), &
401 buf%Pcutoff(1), &
402 r, phi) ! compute just pair potential
403 end if
404
405 ! contribution to energy
406 !
407 if (comp_enepot == 1) then
408 enepot(i) = enepot(i) + 0.5_cd * phi ! accumulate energy
409 end if
410 if (comp_energy == 1) then
411 energy = energy + 0.5_cd * phi ! add half v to total energy
412 end if
413
414 ! contribution to process_dEdr
415 !
416 if (comp_process_dedr == 1) then
417 call kim_process_dedr_term( &
418 model_compute_arguments_handle, deidr, r, rij, i, j, ierr)
419 end if
420
421 ! contribution to process_d2Edr2
422 if (comp_process_d2edr2 == 1) then
423 r_pairs(1) = r
424 r_pairs(2) = r
425 rij_pairs(:, 1) = rij
426 rij_pairs(:, 2) = rij
427 i_pairs(1) = i
428 i_pairs(2) = i
429 j_pairs(1) = j
430 j_pairs(2) = j
431
432 call kim_process_d2edr2_term( &
433 model_compute_arguments_handle, d2eidr, &
434 r_pairs, rij_pairs, i_pairs, j_pairs, ierr)
435 end if
436
437 ! contribution to forces
438 !
439 if (comp_force == 1) then
440 force(:, i) = force(:, i) + deidr * rij / r ! accumulate force on particle i
441 force(:, j) = force(:, j) - deidr * rij / r ! accumulate force on particle j
442 end if
443
444 end if
445
446 end do ! loop on jj
447
448 end if ! if particleContributing
449
450 end do ! do i
451
452 ! Everything is great
453 !
454 ierr = 0
455 return
456
457 end subroutine compute_energy_forces
458
459 !-----------------------------------------------------------------------------
460 !
461 ! Model driver refresh routine
462 !
463 !-----------------------------------------------------------------------------
464 recursive subroutine refresh(model_refresh_handle, ierr) bind(c)
465 implicit none
466
467 !-- transferred variables
468 type(kim_model_refresh_handle_type), intent(in) :: model_refresh_handle
469 integer(c_int), intent(out) :: ierr
470
471 !-- Local variables
472 real(c_double) energy_at_cutoff
473 type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
474
475 ! get model buffer from KIM object
476 call kim_get_model_buffer_pointer(model_refresh_handle, pbuf)
477 call c_f_pointer(pbuf, buf)
478
479 call kim_set_influence_distance_pointer(model_refresh_handle, &
480 buf%influence_distance(1))
481 call kim_set_neighbor_list_pointers( &
482 model_refresh_handle, 1, buf%influence_distance, &
483 buf%model_will_not_request_neighbors_of_noncontributing_particles)
484
485 ! Set new values in KIM object and buffer
486 !
487 buf%influence_distance(1) = buf%Pcutoff(1)
488 buf%cutsq(1) = (buf%Pcutoff(1))**2
489 ! calculate pair potential at r=cutoff with shift=0.0
490 call calc_phi(buf%epsilon(1), &
491 buf%sigma(1), &
492 0.0_cd, &
493 buf%Pcutoff(1), &
494 buf%Pcutoff(1), energy_at_cutoff)
495 buf%shift(1) = -energy_at_cutoff
496
497 ierr = 0
498 return
499
500 end subroutine refresh
501
502 !-----------------------------------------------------------------------------
503 !
504 ! Model driver write_model routine
505 !
506 !-----------------------------------------------------------------------------
507 recursive subroutine write_model( &
508 model_write_parameterized_model_handle, ierr) bind(c)
509 implicit none
510
511 !-- transferred variables
512 type(kim_model_write_parameterized_model_handle_type), intent(in) &
513 :: model_write_parameterized_model_handle
514 integer(c_int), intent(out) :: ierr
515
516 !-- Local variables
517 integer i
518 type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
519 character(len=512, kind=c_char) :: path
520 character(len=512, kind=c_char) :: model_name
521 character(len=512, kind=c_char) :: string_buffer
522 character(len=100, kind=c_char) :: species_name
523
524 ! get model buffer from KIM object
525 call kim_get_model_buffer_pointer( &
526 model_write_parameterized_model_handle, pbuf)
527 call c_f_pointer(pbuf, buf)
528
529 call kim_get_path(model_write_parameterized_model_handle, path)
530 call kim_get_model_name(model_write_parameterized_model_handle, model_name)
531
532 write (string_buffer, '(A)') trim(model_name)//".params"
533 call kim_set_parameter_file_name(model_write_parameterized_model_handle, &
534 string_buffer)
535 write (string_buffer, '(A)') trim(path)//"/"//trim(string_buffer)
536
537 open (42, file=trim(string_buffer), &
538 status="REPLACE", action="WRITE", iostat=ierr)
539 if (ierr /= 0) then
540 call kim_log_entry( &
541 model_write_parameterized_model_handle, kim_log_verbosity_error, &
542 "Unable to open parameter file for writing.")
543 return
544 end if
545
546 do i = 1, 100
547 species_name(i:i) = buf%species_name(i)
548 end do
549 write (42, '(A)') trim(species_name)
550 write (42, '(ES20.10)') buf%Pcutoff(1)
551 write (42, '(ES20.10)') buf%epsilon(1)
552 write (42, '(ES20.10)') buf%sigma(1)
553
554 ierr = 0
555 return
556
557 end subroutine write_model
558
559 !-----------------------------------------------------------------------------
560 !
561 ! Model driver destroy routine
562 !
563 !-----------------------------------------------------------------------------
564 recursive subroutine destroy(model_destroy_handle, ierr) bind(c)
565 implicit none
566
567 !-- Transferred variables
568 type(kim_model_destroy_handle_type), intent(in) :: model_destroy_handle
569 integer(c_int), intent(out) :: ierr
570
571 !-- Local variables
572 type(buffer_type), pointer :: buf; type(c_ptr) :: pbuf
573
574 ! get model buffer from KIM object
575 call kim_get_model_buffer_pointer(model_destroy_handle, pbuf)
576 call c_f_pointer(pbuf, buf)
577
578 deallocate (buf)
579
580 ierr = 0
581 return
582
583 end subroutine destroy
584
585 !-----------------------------------------------------------------------------
586 !
587 ! Model driver compute arguments create routine
588 !
589 !-----------------------------------------------------------------------------
590 recursive subroutine compute_arguments_create( &
591 model_compute_handle, model_compute_arguments_create_handle, ierr) bind(c)
592 implicit none
593
594 !-- Transferred variables
595 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
596 type(kim_model_compute_arguments_create_handle_type), intent(in) :: &
597 model_compute_arguments_create_handle
598 integer(c_int), intent(out) :: ierr
599
600 integer(c_int) ierr2
601
602 ! avoid unsed dummy argument warnings
603 if (model_compute_handle == kim_model_compute_null_handle) continue
604
605 ierr = 0
606 ierr2 = 0
607
608 ! register arguments
609 call kim_set_argument_support_status( &
610 model_compute_arguments_create_handle, &
611 kim_compute_argument_name_partial_energy, &
612 kim_support_status_optional, ierr)
613 call kim_set_argument_support_status( &
614 model_compute_arguments_create_handle, &
615 kim_compute_argument_name_partial_forces, &
616 kim_support_status_optional, ierr2)
617 ierr = ierr + ierr2
618 call kim_set_argument_support_status( &
619 model_compute_arguments_create_handle, &
620 kim_compute_argument_name_partial_particle_energy, &
621 kim_support_status_optional, ierr2)
622 ierr = ierr + ierr2
623 if (ierr /= 0) then
624 call kim_log_entry( &
625 model_compute_arguments_create_handle, kim_log_verbosity_error, &
626 "Unable to register arguments support_statuss")
627 ierr = 1
628 goto 42
629 end if
630
631 ! register callbacks
632 call kim_set_callback_support_status( &
633 model_compute_arguments_create_handle, &
634 kim_compute_callback_name_process_dedr_term, &
635 kim_support_status_optional, ierr)
636 call kim_set_callback_support_status( &
637 model_compute_arguments_create_handle, &
638 kim_compute_callback_name_process_d2edr2_term, &
639 kim_support_status_optional, ierr2)
640 ierr = ierr + ierr2
641 if (ierr /= 0) then
642 call kim_log_entry( &
643 model_compute_arguments_create_handle, kim_log_verbosity_error, &
644 "Unable to register callbacks support_statuss")
645 ierr = 1
646 goto 42
647 end if
648
649 ierr = 0
65042 continue
651 return
652
653 end subroutine compute_arguments_create
654
655 !-----------------------------------------------------------------------------
656 !
657 ! Model driver compute arguments destroy routine
658 !
659 !-----------------------------------------------------------------------------
660 recursive subroutine compute_arguments_destroy( &
661 model_compute_handle, model_compute_arguments_destroy_handle, ierr) bind(c)
662 implicit none
663
664 !-- Transferred variables
665 type(kim_model_compute_handle_type), intent(in) :: model_compute_handle
666 type(kim_model_compute_arguments_destroy_handle_type), intent(in) :: &
667 model_compute_arguments_destroy_handle
668 integer(c_int), intent(out) :: ierr
669
670 ! avoid unsed dummy argument warnings
671 if (model_compute_handle == kim_model_compute_null_handle) continue
672 if (model_compute_arguments_destroy_handle == &
673 kim_model_compute_arguments_destroy_null_handle) continue
674
675 ierr = 0
676 return
677 end subroutine compute_arguments_destroy
678
679end module ex_model_driver_p_lj
680
681!-----------------------------------------------------------------------------
682!
683! Model driver create routine (REQUIRED)
684!
685!-----------------------------------------------------------------------------
686recursive subroutine model_driver_create_routine( &
687 model_driver_create_handle, requested_length_unit, requested_energy_unit, &
688 requested_charge_unit, requested_temperature_unit, requested_time_unit, &
689 ierr) bind(c)
690 use, intrinsic :: iso_c_binding
693 implicit none
694 integer(c_int), parameter :: cd = c_double ! used for literal constants
695
696 !-- Transferred variables
697 type(kim_model_driver_create_handle_type), intent(in) &
698 :: model_driver_create_handle
699 type(kim_length_unit_type), intent(in), value :: requested_length_unit
700 type(kim_energy_unit_type), intent(in), value :: requested_energy_unit
701 type(kim_charge_unit_type), intent(in), value :: requested_charge_unit
702 type(kim_temperature_unit_type), intent(in), value :: &
703 requested_temperature_unit
704 type(kim_time_unit_type), intent(in), value :: requested_time_unit
705 integer(c_int), intent(out) :: ierr
706
707 !-- Local variables
708 integer i
709 integer(c_int) :: number_of_parameter_files
710 character(len=1024, kind=c_char) :: parameter_file_directory_name
711 character(len=1024, kind=c_char) :: parameter_file_basename
712 character(len=1024, kind=c_char) :: parameter_file_name
713 integer(c_int) :: ierr2
714 type(buffer_type), pointer :: buf
715 type(kim_species_name_type) species_name
716 ! define variables for all model parameters to be read in
717 real(c_double) factor
718 character(len=100, kind=c_char) in_species
719 real(c_double) in_cutoff
720 real(c_double) in_epsilon
721 real(c_double) in_sigma
722 real(c_double) energy_at_cutoff
723
724 ! register units
725 call kim_set_units( &
726 model_driver_create_handle, &
727 requested_length_unit, &
728 requested_energy_unit, &
729 kim_charge_unit_unused, &
730 kim_temperature_unit_unused, &
731 kim_time_unit_unused, ierr)
732 if (ierr /= 0) then
733 call kim_log_entry(model_driver_create_handle, &
734 kim_log_verbosity_error, "Unable to set units")
735 ierr = 1
736 goto 42
737 end if
738
739 ! register numbering
740 call kim_set_model_numbering( &
741 model_driver_create_handle, kim_numbering_one_based, ierr)
742 if (ierr /= 0) then
743 call kim_log_entry(model_driver_create_handle, &
744 kim_log_verbosity_error, "Unable to set numbering")
745 ierr = 1
746 goto 42
747 end if
748
749 ! store callback pointers in KIM object
750 call kim_set_routine_pointer( &
751 model_driver_create_handle, kim_model_routine_name_compute, &
752 kim_language_name_fortran, 1, c_funloc(compute_energy_forces), ierr)
753 call kim_set_routine_pointer( &
754 model_driver_create_handle, &
755 kim_model_routine_name_compute_arguments_create, &
756 kim_language_name_fortran, 1, c_funloc(compute_arguments_create), ierr2)
757 ierr = ierr + ierr2
758 call kim_set_routine_pointer( &
759 model_driver_create_handle, kim_model_routine_name_refresh, &
760 kim_language_name_fortran, 1, c_funloc(refresh), ierr2)
761 ierr = ierr + ierr2
762 call kim_set_routine_pointer( &
763 model_driver_create_handle, &
764 kim_model_routine_name_write_parameterized_model, &
765 kim_language_name_fortran, 0, c_funloc(write_model), ierr2)
766 ierr = ierr + ierr2
767 call kim_set_routine_pointer( &
768 model_driver_create_handle, &
769 kim_model_routine_name_compute_arguments_destroy, &
770 kim_language_name_fortran, 1, c_funloc(compute_arguments_destroy), ierr2)
771 ierr = ierr + ierr2
772 call kim_set_routine_pointer( &
773 model_driver_create_handle, kim_model_routine_name_destroy, &
774 kim_language_name_fortran, 1, c_funloc(destroy), ierr2)
775 ierr = ierr + ierr2
776 if (ierr /= 0) then
777 call kim_log_entry( &
778 model_driver_create_handle, kim_log_verbosity_error, &
779 "Unable to store callback pointers")
780 ierr = 1
781 goto 42
782 end if
783
784 ! process parameter files
785 call kim_get_number_of_parameter_files( &
786 model_driver_create_handle, number_of_parameter_files)
787 if (number_of_parameter_files /= 1) then
788 call kim_log_entry( &
789 model_driver_create_handle, kim_log_verbosity_error, &
790 "Wrong number of parameter files")
791 ierr = 1
792 goto 42
793 end if
794
795 ! Read in model parameters from parameter file
796 !
797 call kim_get_parameter_file_directory_name( &
798 model_driver_create_handle, parameter_file_directory_name)
799 call kim_get_parameter_file_basename( &
800 model_driver_create_handle, 1, parameter_file_basename, ierr)
801 if (ierr /= 0) then
802 call kim_log_entry( &
803 model_driver_create_handle, kim_log_verbosity_error, &
804 "Unable to get parameter file basename")
805 ierr = 1
806 goto 42
807 end if
808 parameter_file_name = trim(parameter_file_directory_name)//"/"// &
809 trim(parameter_file_basename)
810 open (10, file=parameter_file_name, status="old")
811 read (10, *, iostat=ierr, err=100) in_species
812 read (10, *, iostat=ierr, err=100) in_cutoff
813 read (10, *, iostat=ierr, err=100) in_epsilon
814 read (10, *, iostat=ierr, err=100) in_sigma
815 close (10)
816
817 goto 200
818100 continue
819 ! reading parameters failed
820 call kim_log_entry(model_driver_create_handle, &
821 kim_log_verbosity_error, "Unable to read LJ parameters")
822 ierr = 1
823 goto 42
824
825200 continue
826
827 ! register species
828 call kim_from_string(in_species, species_name)
829 call kim_set_species_code( &
830 model_driver_create_handle, species_name, speccode, ierr)
831 if (ierr /= 0) then
832 call kim_log_entry(model_driver_create_handle, &
833 kim_log_verbosity_error, "Unable to set species code")
834 ierr = 1
835 goto 42
836 end if
837
838 ! convert to appropriate units
839 call kim_convert_unit( &
840 kim_length_unit_a, &
841 kim_energy_unit_ev, &
842 kim_charge_unit_e, &
843 kim_temperature_unit_k, &
844 kim_time_unit_ps, &
845 requested_length_unit, &
846 requested_energy_unit, &
847 requested_charge_unit, &
848 requested_temperature_unit, &
849 requested_time_unit, &
850 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
851 if (ierr /= 0) then
852 call kim_log_entry(model_driver_create_handle, &
853 kim_log_verbosity_error, "kim_api_convert_to_act_unit")
854 ierr = 1
855 goto 42
856 end if
857 in_cutoff = in_cutoff * factor
858
859 call kim_convert_unit( &
860 kim_length_unit_a, &
861 kim_energy_unit_ev, &
862 kim_charge_unit_e, &
863 kim_temperature_unit_k, &
864 kim_time_unit_ps, &
865 requested_length_unit, &
866 requested_energy_unit, &
867 requested_charge_unit, &
868 requested_temperature_unit, &
869 requested_time_unit, &
870 0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
871 if (ierr /= 0) then
872 call kim_log_entry(model_driver_create_handle, &
873 kim_log_verbosity_error, "kim_api_convert_to_act_unit")
874 ierr = 1
875 goto 42
876 end if
877 in_epsilon = in_epsilon * factor
878
879 call kim_convert_unit( &
880 kim_length_unit_a, &
881 kim_energy_unit_ev, &
882 kim_charge_unit_e, &
883 kim_temperature_unit_k, &
884 kim_time_unit_ps, &
885 requested_length_unit, &
886 requested_energy_unit, &
887 requested_charge_unit, &
888 requested_temperature_unit, &
889 requested_time_unit, &
890 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
891 if (ierr /= 0) then
892 call kim_log_entry(model_driver_create_handle, &
893 kim_log_verbosity_error, "kim_api_convert_to_act_unit")
894 ierr = 1
895 goto 42
896 end if
897 in_sigma = in_sigma * factor
898
899 allocate (buf)
900
901 ! setup buffer
902 ! set value of parameters
903 do i = 1, 100
904 buf%species_name(i) = in_species(i:i)
905 end do
906 buf%influence_distance(1) = in_cutoff
907 buf%Pcutoff(1) = in_cutoff
908 buf%cutsq(1) = in_cutoff**2
909 buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
910 buf%epsilon(1) = in_epsilon
911 buf%sigma(1) = in_sigma
912 call calc_phi(in_epsilon, &
913 in_sigma, &
914 0.0_cd, &
915 in_cutoff, &
916 in_cutoff, energy_at_cutoff)
917 buf%shift(1) = -energy_at_cutoff
918
919 ! store model cutoff in KIM object
920 call kim_set_influence_distance_pointer( &
921 model_driver_create_handle, buf%influence_distance(1))
922 call kim_set_neighbor_list_pointers( &
923 model_driver_create_handle, 1, buf%influence_distance, &
924 buf%model_will_not_request_neighbors_of_noncontributing_particles)
925
926 ! end setup buffer
927
928 ! store in model buffer
929 call kim_set_model_buffer_pointer( &
930 model_driver_create_handle, c_loc(buf))
931
932 ! set pointers to parameters in KIM object
933 call kim_set_parameter_pointer( &
934 model_driver_create_handle, buf%pcutoff, "cutoff", &
935 "Distance beyond which particles do not interact with one another.", ierr)
936 if (ierr /= 0) then
937 call kim_log_entry(model_driver_create_handle, &
938 kim_log_verbosity_error, "set_parameter")
939 ierr = 1
940 goto 42
941 end if
942
943 call kim_set_parameter_pointer( &
944 model_driver_create_handle, buf%epsilon, "epsilon", &
945 "Maximum depth of the potential well.", ierr)
946 if (ierr /= 0) then
947 call kim_log_entry(model_driver_create_handle, &
948 kim_log_verbosity_error, "set_parameter")
949 ierr = 1
950 goto 42
951 end if
952
953 call kim_set_parameter_pointer( &
954 model_driver_create_handle, buf%sigma, "sigma", &
955 "Distance at which energy is zero and force is repulsive.", ierr)
956 if (ierr /= 0) then
957 call kim_log_entry(model_driver_create_handle, &
958 kim_log_verbosity_error, "set_parameter")
959 ierr = 1
960 goto 42
961 end if
962
963 ierr = 0
96442 continue
965 return
966
967end subroutine model_driver_create_routine
recursive subroutine, public refresh(model_refresh_handle, ierr)
recursive subroutine, public compute_energy_forces(model_compute_handle, model_compute_arguments_handle, ierr)
recursive subroutine, public calc_phi_dphi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi, dphi)
integer(c_int), parameter, public speccode
recursive subroutine, public calc_phi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi)
recursive subroutine, public calc_phi_dphi_d2phi(model_epsilon, model_sigma, model_shift, model_cutoff, r, phi, dphi, d2phi)
recursive subroutine, public destroy(model_destroy_handle, ierr)
recursive subroutine, public write_model(model_write_parameterized_model_handle, ierr)