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