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