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