kim-api  2.1.4-git+v2.1.3-git-1-g7847914a.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_name
711 integer(c_int) :: ierr2
712 type(buffer_type), pointer :: buf;
713 type(kim_species_name_type) species_name
714 ! define variables for all model parameters to be read in
715 real(c_double) factor
716 character(len=100, kind=c_char) in_species
717 real(c_double) in_cutoff
718 real(c_double) in_epsilon
719 real(c_double) in_sigma
720 real(c_double) energy_at_cutoff
721 
722 ! register units
723 call kim_set_units( &
724  model_driver_create_handle, &
725  requested_length_unit, &
726  requested_energy_unit, &
727  kim_charge_unit_unused, &
728  kim_temperature_unit_unused, &
729  kim_time_unit_unused, ierr)
730 if (ierr /= 0) then
731  call kim_log_entry(model_driver_create_handle, &
732  kim_log_verbosity_error, "Unable to set units")
733  ierr = 1
734  goto 42
735 end if
736 
737 ! register numbering
738 call kim_set_model_numbering( &
739  model_driver_create_handle, kim_numbering_one_based, ierr)
740 if (ierr /= 0) then
741  call kim_log_entry(model_driver_create_handle, &
742  kim_log_verbosity_error, "Unable to set numbering")
743  ierr = 1
744  goto 42
745 end if
746 
747 
748 ! store callback pointers in KIM object
749 call kim_set_routine_pointer( &
750  model_driver_create_handle, kim_model_routine_name_compute, &
751  kim_language_name_fortran, 1, c_funloc(compute_energy_forces), ierr)
752 call kim_set_routine_pointer( &
753  model_driver_create_handle, kim_model_routine_name_compute_arguments_create, &
754  kim_language_name_fortran, 1, c_funloc(compute_arguments_create), ierr2)
755 ierr = ierr + ierr2
756 call kim_set_routine_pointer( &
757  model_driver_create_handle, kim_model_routine_name_refresh, &
758  kim_language_name_fortran, 1, c_funloc(refresh), ierr2)
759 ierr = ierr + ierr2
760 call kim_set_routine_pointer( &
761  model_driver_create_handle, &
762  kim_model_routine_name_write_parameterized_model, &
763  kim_language_name_fortran, 0, c_funloc(write_model), ierr2)
764 ierr = ierr + ierr2
765 call kim_set_routine_pointer( &
766  model_driver_create_handle, &
767  kim_model_routine_name_compute_arguments_destroy, &
768  kim_language_name_fortran, 1, c_funloc(compute_arguments_destroy), ierr2)
769 ierr = ierr + ierr2
770 call kim_set_routine_pointer( &
771  model_driver_create_handle, kim_model_routine_name_destroy, &
772  kim_language_name_fortran, 1, c_funloc(destroy), ierr2)
773 ierr = ierr + ierr2
774 if (ierr /= 0) then
775  call kim_log_entry(model_driver_create_handle, &
776  kim_log_verbosity_error, "Unable to store callback pointers")
777  ierr = 1
778  goto 42
779 end if
780 
781 
782 ! process parameter files
783 call kim_get_number_of_parameter_files( &
784  model_driver_create_handle, number_of_parameter_files)
785 if (number_of_parameter_files .ne. 1) then
786  call kim_log_entry(model_driver_create_handle, &
787  kim_log_verbosity_error, "Wrong number of parameter files")
788  ierr = 1
789  goto 42
790 end if
791 
792 ! Read in model parameters from parameter file
793 !
794 call kim_get_parameter_file_name( &
795  model_driver_create_handle, 1, parameter_file_name, ierr)
796 if (ierr /= 0) then
797  call kim_log_entry(model_driver_create_handle, &
798  kim_log_verbosity_error, "Unable to get parameter file name")
799  ierr = 1
800  goto 42
801 end if
802 open(10,file=parameter_file_name,status="old")
803 read(10,*,iostat=ierr,err=100) in_species
804 read(10,*,iostat=ierr,err=100) in_cutoff
805 read(10,*,iostat=ierr,err=100) in_epsilon
806 read(10,*,iostat=ierr,err=100) in_sigma
807 close(10)
808 
809 goto 200
810 100 continue
811 ! reading parameters failed
812 call kim_log_entry(model_driver_create_handle, &
813  kim_log_verbosity_error, "Unable to read LJ parameters")
814 ierr = 1
815 goto 42
816 
817 200 continue
818 
819 
820 ! register species
821 call kim_from_string(in_species, species_name)
822 call kim_set_species_code( &
823  model_driver_create_handle, species_name, speccode, ierr)
824 if (ierr /= 0) then
825  call kim_log_entry(model_driver_create_handle, &
826  kim_log_verbosity_error, "Unable to set species code")
827  ierr = 1
828  goto 42
829 end if
830 
831 ! convert to appropriate units
832 call kim_convert_unit( &
833  kim_length_unit_a, &
834  kim_energy_unit_ev, &
835  kim_charge_unit_e, &
836  kim_temperature_unit_k, &
837  kim_time_unit_ps, &
838  requested_length_unit, &
839  requested_energy_unit, &
840  requested_charge_unit, &
841  requested_temperature_unit, &
842  requested_time_unit, &
843  1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
844 if (ierr /= 0) then
845  call kim_log_entry(model_driver_create_handle, &
846  kim_log_verbosity_error, "kim_api_convert_to_act_unit")
847  ierr = 1
848  goto 42
849 endif
850 in_cutoff = in_cutoff * factor
851 
852 call kim_convert_unit( &
853  kim_length_unit_a, &
854  kim_energy_unit_ev, &
855  kim_charge_unit_e, &
856  kim_temperature_unit_k, &
857  kim_time_unit_ps, &
858  requested_length_unit, &
859  requested_energy_unit, &
860  requested_charge_unit, &
861  requested_temperature_unit, &
862  requested_time_unit, &
863  0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
864 if (ierr /= 0) then
865  call kim_log_entry(model_driver_create_handle, &
866  kim_log_verbosity_error, "kim_api_convert_to_act_unit")
867  ierr = 1
868  goto 42
869 endif
870 in_epsilon = in_epsilon * factor
871 
872 call kim_convert_unit( &
873  kim_length_unit_a, &
874  kim_energy_unit_ev, &
875  kim_charge_unit_e, &
876  kim_temperature_unit_k, &
877  kim_time_unit_ps, &
878  requested_length_unit, &
879  requested_energy_unit, &
880  requested_charge_unit, &
881  requested_temperature_unit, &
882  requested_time_unit, &
883  1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, factor, ierr)
884 if (ierr /= 0) then
885  call kim_log_entry(model_driver_create_handle, &
886  kim_log_verbosity_error, "kim_api_convert_to_act_unit")
887  ierr = 1
888  goto 42
889 endif
890 in_sigma = in_sigma * factor
891 
892 allocate( buf )
893 
894 ! setup buffer
895 ! set value of parameters
896 do i=1,100
897  buf%species_name(i) = in_species(i:i)
898 end do
899 buf%influence_distance(1) = in_cutoff
900 buf%Pcutoff(1) = in_cutoff
901 buf%cutsq(1) = in_cutoff**2
902 buf%model_will_not_request_neighbors_of_noncontributing_particles = 1
903 buf%epsilon(1) = in_epsilon
904 buf%sigma(1) = in_sigma
905 call calc_phi(in_epsilon, &
906  in_sigma, &
907  0.0_cd, &
908  in_cutoff, &
909  in_cutoff, energy_at_cutoff)
910 buf%shift(1) = -energy_at_cutoff
911 
912 ! store model cutoff in KIM object
913 call kim_set_influence_distance_pointer( &
914  model_driver_create_handle, buf%influence_distance(1))
915 call kim_set_neighbor_list_pointers( &
916  model_driver_create_handle, 1, buf%influence_distance, &
917  buf%model_will_not_request_neighbors_of_noncontributing_particles)
918 
919 ! end setup buffer
920 
921 ! store in model buffer
922 call kim_set_model_buffer_pointer( &
923  model_driver_create_handle, c_loc(buf))
924 
925 ! set pointers to parameters in KIM object
926 call kim_set_parameter_pointer( &
927  model_driver_create_handle, buf%pcutoff, "cutoff", &
928  "Distance beyond which particles do not interact with one another.", ierr)
929 if (ierr /= 0) then
930  call kim_log_entry(model_driver_create_handle, &
931  kim_log_verbosity_error, "set_parameter")
932  ierr = 1
933  goto 42
934 endif
935 
936 call kim_set_parameter_pointer( &
937  model_driver_create_handle, buf%epsilon, "epsilon", &
938  "Maximum depth of the potential well.", ierr)
939 if (ierr /= 0) then
940  call kim_log_entry(model_driver_create_handle, &
941  kim_log_verbosity_error, "set_parameter")
942  ierr = 1
943  goto 42
944 endif
945 
946 call kim_set_parameter_pointer( &
947  model_driver_create_handle, buf%sigma, "sigma", &
948  "Distance at which energy is zero and force is repulsive.", ierr)
949 if (ierr /= 0) then
950  call kim_log_entry(model_driver_create_handle, &
951  kim_log_verbosity_error, "set_parameter")
952  ierr = 1
953  goto 42
954 endif
955 
956 ierr = 0
957 42 continue
958 return
959 
960 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)