! ! CDDL HEADER START ! ! The contents of this file are subject to the terms of the Common Development ! and Distribution License Version 1.0 (the "License"). ! ! You can obtain a copy of the license at ! http://www.opensource.org/licenses/CDDL-1.0. See the License for the ! specific language governing permissions and limitations under the License. ! ! When distributing Covered Code, include this CDDL HEADER in each file and ! include the License file in a prominent location with the name LICENSE.CDDL. ! If applicable, add the following below this CDDL HEADER, with the fields ! enclosed by brackets "[]" replaced with your own identifying information: ! ! Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved. ! ! CDDL HEADER END ! ! ! Copyright (c) 2013, Regents of the University of Minnesota. ! All rights reserved. ! ! Contributors: ! Ellad B. Tadmor ! Ryan S. Elliott ! Stephen M. Whalen ! !**************************************************************************** !** !** MODULE Pair_Lennard_Jones_Truncated !** !** Lennard-Jones pair potential KIM Model Driver !** truncated to have zero energy above the cutoff radius !** !** Note: The energy and its derivatives will be discontinuous at the cutoff. !** !** Language: Fortran 2003 !** !**************************************************************************** module Pair_Lennard_Jones_Truncated use, intrinsic :: iso_c_binding use kim_model_driver_headers_module implicit none save private public BUFFER_TYPE, & compute, & compute_arguments_create, & compute_arguments_destroy, & refresh, & destroy, & calc_phi, & calc_phi_dphi, & calc_phi_dphi_d2phi, & speccode ! Below are the definitions and values of all Model parameters integer(c_int), parameter :: cd = c_double ! for literal constants integer(c_int), parameter :: DIM=3 ! dimensionality of space integer(c_int), parameter :: speccode = 1 ! internal species code !------------------------------------------------------------------------------- ! ! Definition of Buffer type ! !------------------------------------------------------------------------------- type, bind(c) :: BUFFER_TYPE real(c_double) :: influence_distance real(c_double) :: cutoff(1) real(c_double) :: cutsq(1) real(c_double) :: eps(1) real(c_double) :: sigma(1) integer(c_int) :: & model_will_not_request_neighbors_of_noncontributing_particles(1) endtype BUFFER_TYPE contains !------------------------------------------------------------------------------- ! ! Calculate pair potential phi(r) ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine calc_phi(model_eps, & model_sigma, & model_cutoff,r,phi) implicit none !-- Transferred variables real(c_double), intent(in) :: model_eps real(c_double), intent(in) :: model_sigma real(c_double), intent(in) :: model_cutoff real(c_double), intent(in) :: r real(c_double), intent(out) :: phi !-- Local variables real(c_double) rsq,sor,sor6,sor12 if (r .gt. model_cutoff) then ! Argument exceeds cutoff radius phi = 0.0_cd return endif rsq = r*r ! r^2 sor = model_sigma/r ! (sig/r) sor6 = sor*sor*sor ! sor6 = sor6*sor6 ! (sig/r)^6 sor12= sor6*sor6 ! (sig/r)^12 phi = 4.0_cd*model_eps*(sor12-sor6) return end subroutine calc_phi !------------------------------------------------------------------------------- ! ! Calculate pair potential phi(r) and its derivative dphi(r) ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine calc_phi_dphi(model_eps, & model_sigma, & model_cutoff,r,phi,dphi) implicit none !-- Transferred variables real(c_double), intent(in) :: model_eps real(c_double), intent(in) :: model_sigma real(c_double), intent(in) :: model_cutoff real(c_double), intent(in) :: r real(c_double), intent(out) :: phi,dphi !-- Local variables real(c_double) rsq,sor,sor6,sor12 if (r .gt. model_cutoff) then ! Argument exceeds cutoff radius phi = 0.0_cd dphi = 0.0_cd return endif rsq = r*r ! r^2 sor = model_sigma/r ! (sig/r) sor6 = sor*sor*sor ! sor6 = sor6*sor6 ! (sig/r)^6 sor12= sor6*sor6 ! (sig/r)^12 phi = 4.0_cd*model_eps*(sor12-sor6) dphi = 24.0_cd*model_eps*(-2.0_cd*sor12+sor6)/r return end subroutine calc_phi_dphi !------------------------------------------------------------------------------- ! ! Calculate pair potential phi(r) and its derivatives dphi(r) and d2phi(r) ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine calc_phi_dphi_d2phi(model_eps, & model_sigma, & model_cutoff,r,phi,dphi,d2phi) implicit none !-- Transferred variables real(c_double), intent(in) :: model_eps real(c_double), intent(in) :: model_sigma real(c_double), intent(in) :: model_cutoff real(c_double), intent(in) :: r real(c_double), intent(out) :: phi,dphi,d2phi !-- Local variables real(c_double) rsq,sor,sor6,sor12 if (r .gt. model_cutoff) then ! Argument exceeds cutoff radius phi = 0.0_cd dphi = 0.0_cd d2phi = 0.0_cd return endif rsq = r*r ! r^2 sor = model_sigma/r ! (sig/r) sor6 = sor*sor*sor ! sor6 = sor6*sor6 ! (sig/r)^6 sor12= sor6*sor6 ! (sig/r)^12 phi = 4.0_cd*model_eps*(sor12-sor6) dphi = 24.0_cd*model_eps*(-2.0_cd*sor12+sor6)/r d2phi = 24.0_cd*model_eps*(26.0_cd*sor12-7.0_cd*sor6)/rsq return end subroutine calc_phi_dphi_d2phi !------------------------------------------------------------------------------- ! ! Compute ! !-------------------------------------------------------------------------------i ! The "recursive" keyword is added below make this routine thread safe recursive subroutine compute(model_compute_handle, & model_compute_arguments_handle, ierr) bind(c) implicit none !-- Transferred variables type(kim_model_compute_handle_type), intent(in) :: model_compute_handle type(kim_model_compute_arguments_handle_type), intent(in) :: & model_compute_arguments_handle integer(c_int), intent(out) :: ierr !-- Local variables real(c_double) :: r,Rsqij,phi,dphi,d2phi,dEidr,d2Eidr real(c_double) :: v real(c_double) :: vir(6) integer(c_int) :: i,j,jj,numnei integer(c_int) :: ierr2 integer(c_int) :: comp_force,comp_energy,comp_particleEnergy,comp_process_dEdr, & comp_process_d2Edr2,comp_virial,comp_particleVirial type(BUFFER_TYPE), pointer :: buf; type(c_ptr) :: pbuf real(c_double), target :: Rij(DIM) ! Quantities for computing d2Edr2 real(c_double), target :: Rij_pairs(DIM,2) real(c_double), target :: r_pairs(2) integer(c_int), target :: i_pairs(2), j_pairs(2) !-- KIM variables integer(c_int), pointer :: N real(c_double), pointer :: energy real(c_double), pointer :: coor(:,:) real(c_double), pointer :: force(:,:) real(c_double), pointer :: particleEnergy(:) real(c_double), pointer :: virial(:) real(c_double), pointer :: particleVirial(:, :) integer(c_int), pointer :: nei1part(:) integer(c_int), pointer :: particleSpeciesCodes(:) integer(c_int), pointer :: particleContributing(:) ! get model buffer from KIM object call kim_get_model_buffer_pointer(model_compute_handle, pbuf) call c_f_pointer(pbuf, buf) ! Unpack data from KIM object ! ierr = 0_c_int call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_NUMBER_OF_PARTICLES, N, ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_SPECIES_CODES, & N, particlespeciesCodes, ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_CONTRIBUTING, N, particlecontributing, & ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_COORDINATES, DIM, N, coor, ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, energy, ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, DIM, N, force, ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_ENERGY, N, particleEnergy, ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_VIRIAL, & 6_c_int, & virial, & ierr2) ierr = ierr + ierr2 call kim_get_argument_pointer( & model_compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_VIRIAL, & 6_c_int, & N, & particleVirial, & ierr2) ierr = ierr + ierr2 if (ierr /= 0_c_int) then call kim_log_entry(model_compute_arguments_handle, & KIM_LOG_VERBOSITY_ERROR, "get_argument_pointer") return endif ! Check to see if we have been asked to compute the energy, forces, energy per particle, ! dEdr, and d2Edr2 ! if (associated(energy)) then comp_energy = 1_c_int else comp_energy = 0_c_int end if if (associated(force)) then comp_force = 1_c_int else comp_force = 0_c_int end if if (associated(particleEnergy)) then comp_particleEnergy = 1_c_int else comp_particleEnergy = 0_c_int end if if (associated(virial)) then comp_virial = 1_c_int else comp_virial = 0_c_int end if if (associated(particleVirial)) then comp_particleVirial = 1_c_int else comp_particleVirial = 0_c_int end if ierr = 0_c_int call kim_is_callback_present( & model_compute_arguments_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_DEDR_TERM, comp_process_dedr, ierr2) ierr = ierr + ierr2 call kim_is_callback_present( & model_compute_arguments_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_D2EDR2_TERM, comp_process_d2edr2, ierr2) ierr = ierr + ierr2 if (ierr /= 0_c_int) then call kim_log_entry(model_compute_arguments_handle, & KIM_LOG_VERBOSITY_ERROR, "get_compute") return endif ! Check to be sure that the species are correct ! ierr = 1_c_int ! assume an error do i = 1_c_int,N if (particleSpeciesCodes(i).ne.speccode) then call kim_log_entry(model_compute_arguments_handle, & KIM_LOG_VERBOSITY_ERROR, "Unexpected species code detected") return endif enddo ierr = 0_c_int ! everything is ok ! Initialize potential energies, forces ! if (comp_particleEnergy.eq.1_c_int) particleEnergy = 0.0_cd if (comp_energy.eq.1_c_int) energy = 0.0_cd if (comp_force.eq.1_c_int) force = 0.0_cd if (comp_virial.eq.1_c_int) virial = 0.0_cd if (comp_particleVirial.eq.1_c_int) particleVirial = 0.0_cd ! ! Compute energy and forces ! ! Loop over particles and compute energy and forces ! do i = 1_c_int,N if (particleContributing(i) == 1_c_int) then ! Set up neighbor list for next particle call kim_get_neighbor_list( & model_compute_arguments_handle, 1_c_int, i, numnei, nei1part, ierr) if (ierr /= 0_c_int) then ! some sort of problem, exit call kim_log_entry(model_compute_arguments_handle, & KIM_LOG_VERBOSITY_ERROR, "GetNeighborList failed") ierr = 1_c_int return endif ! Loop over the neighbors of atom i ! do jj = 1_c_int, numnei j = nei1part(jj) ! get neighbor ID if ( .not.( (particleContributing(j).eq.1_c_int) .and. & j.lt.i) ) then ! effective half list ! compute relative position vector ! Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j ! compute energy and forces ! Rsqij = dot_product(Rij,Rij) ! compute square distance if ( Rsqij .lt. buf%cutsq(1) ) then ! particles are interacting? r = sqrt(Rsqij) ! compute distance if (comp_process_d2Edr2.eq.1_c_int) then call calc_phi_dphi_d2phi(buf%eps(1), & buf%sigma(1), & buf%cutoff(1), & r,phi,dphi,d2phi) ! compute pair potential ! and it derivatives if (particleContributing(j).eq.1_c_int) then dEidr = dphi d2Eidr = d2phi else dEidr = 0.5_cd*dphi d2Eidr = 0.5_cd*d2phi endif elseif ((comp_force.eq.1_c_int).or. & (comp_process_dEdr.eq.1_c_int).or. & (comp_virial.eq.1_c_int).or. & (comp_particleVirial.eq.1_c_int)) then call calc_phi_dphi(buf%eps(1), & buf%sigma(1), & buf%cutoff(1), & r,phi,dphi) ! compute pair potential ! and it derivative if (particleContributing(j).eq.1_c_int) then dEidr = dphi else dEidr = 0.5_cd*dphi endif else call calc_phi(buf%eps(1), & buf%sigma(1), & buf%cutoff(1), & r,phi) ! compute just pair potential endif ! contribution to energy ! if (comp_particleEnergy.eq.1_c_int) then particleEnergy(i) = particleEnergy(i) + 0.5_cd*phi ! accumulate energy if (particleContributing(j).eq.1_c_int) then particleEnergy(j) = particleEnergy(j) + 0.5_cd*phi ! accumulate energy endif endif if (comp_energy.eq.1_c_int) then if (particleContributing(j).eq.1_c_int) then energy = energy + phi ! add half v to total energy else energy = energy + 0.5_cd*phi ! add half v to total energy endif endif ! contribution to process_dEdr ! if (comp_process_dEdr.eq.1_c_int) then call kim_process_dedr_term( & model_compute_arguments_handle, dEidr, r, Rij, i, j, ierr) endif ! contribution to virial and particle virial ! if ((comp_virial.eq.1_c_int).or. & (comp_particleVirial.eq.1_c_int)) then ! Virial has 6 components and is stored as a 6-element ! vector in the following order: xx, yy, zz, yz, xz, xy. v = dEidr / r vir(1) = v * Rij(1) * Rij(1) vir(2) = v * Rij(2) * Rij(2) vir(3) = v * Rij(3) * Rij(3) vir(4) = v * Rij(2) * Rij(3) vir(5) = v * Rij(1) * Rij(3) vir(6) = v * Rij(1) * Rij(2) if (comp_virial.eq.1_c_int) then virial(1) = virial(1) + vir(1) virial(2) = virial(2) + vir(2) virial(3) = virial(3) + vir(3) virial(4) = virial(4) + vir(4) virial(5) = virial(5) + vir(5) virial(6) = virial(6) + vir(6) endif if (comp_particleVirial.eq.1_c_int) then vir = vir * 0.5 particleVirial(1, i) = particleVirial(1, i) + vir(1) particleVirial(2, i) = particleVirial(2, i) + vir(2) particleVirial(3, i) = particleVirial(3, i) + vir(3) particleVirial(4, i) = particleVirial(4, i) + vir(4) particleVirial(5, i) = particleVirial(5, i) + vir(5) particleVirial(6, i) = particleVirial(6, i) + vir(6) particleVirial(1, j) = particleVirial(1, j) + vir(1) particleVirial(2, j) = particleVirial(2, j) + vir(2) particleVirial(3, j) = particleVirial(3, j) + vir(3) particleVirial(4, j) = particleVirial(4, j) + vir(4) particleVirial(5, j) = particleVirial(5, j) + vir(5) particleVirial(6, j) = particleVirial(6, j) + vir(6) endif endif ! contribution to process_d2Edr2 if (comp_process_d2Edr2.eq.1_c_int) then r_pairs(1) = r r_pairs(2) = r Rij_pairs(:,1) = Rij Rij_pairs(:,2) = Rij i_pairs(1) = i i_pairs(2) = i j_pairs(1) = j j_pairs(2) = j call kim_process_d2edr2_term( & model_compute_arguments_handle, d2Eidr, & r_pairs, Rij_pairs, i_pairs, j_pairs, ierr) endif ! contribution to forces ! if (comp_force.eq.1_c_int) then force(:,i) = force(:,i) + dEidr*Rij/r ! accumulate force on atom i force(:,j) = force(:,j) - dEidr*Rij/r ! accumulate force on atom j endif endif ! Check on whether particle jj is interacting with particle i endif ! if ( i.lt.j ) enddo ! loop on jj endif ! Check on whether particle i is contributing enddo ! infinite do loop (terminated by exit statements above) ! No errors ! ierr = 0_c_int return end subroutine compute !------------------------------------------------------------------------------- ! ! Model driver refresh routine ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine refresh(model_refresh_handle, ierr) bind(c) implicit none !-- Transferred variables type(kim_model_refresh_handle_type), intent(inout) :: model_refresh_handle integer(c_int), intent(out) :: ierr !-- Local variables type(BUFFER_TYPE), pointer :: buf; type(c_ptr) :: pbuf ! Get model buffer from KIM object call kim_get_model_buffer_pointer(model_refresh_handle, pbuf) call c_f_pointer(pbuf, buf) call kim_set_influence_distance_pointer(model_refresh_handle, & buf%influence_distance) call kim_set_neighbor_list_pointers(model_refresh_handle, & 1_c_int, buf%cutoff, & buf%model_will_not_request_neighbors_of_noncontributing_particles) ! Set new values in KIM object and buffer ! buf%influence_distance = buf%cutoff(1) buf%cutsq(1) = buf%cutoff(1)*buf%cutoff(1) ierr = 0_c_int return end subroutine refresh !------------------------------------------------------------------------------- ! ! Model driver destroy routine ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine destroy(model_destroy_handle, ierr) bind(c) implicit none !-- Transferred variables type(kim_model_destroy_handle_type), intent(inout) :: model_destroy_handle integer(c_int), intent(out) :: ierr !-- Local variables type(BUFFER_TYPE), pointer :: buf; type(c_ptr) :: pbuf ! get model buffer from KIM object call kim_get_model_buffer_pointer(model_destroy_handle, pbuf) call c_f_pointer(pbuf, buf) deallocate( buf ) ierr = 0_c_int return end subroutine destroy !------------------------------------------------------------------------------- ! ! Model driver compute arguments create routine ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine compute_arguments_create(model_compute_handle, & model_compute_arguments_create_handle, ierr) bind(c) use kim_model_compute_arguments_create_module implicit none !-- Transferred variables type(kim_model_compute_handle_type), intent(in) :: model_compute_handle type(kim_model_compute_arguments_create_handle_type), intent(inout) :: & model_compute_arguments_create_handle integer(c_int), intent(out) :: ierr integer(c_int) ierr2 ! avoid unused dummy argument warnings if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue ierr = 0_c_int ierr2 = 0_c_int ! register arguments call kim_set_argument_support_status( & model_compute_arguments_create_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, & KIM_SUPPORT_STATUS_OPTIONAL, ierr) call kim_set_argument_support_status( & model_compute_arguments_create_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_FORCES, & KIM_SUPPORT_STATUS_OPTIONAL, ierr2) ierr = ierr + ierr2 call kim_set_argument_support_status( & model_compute_arguments_create_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_ENERGY, & KIM_SUPPORT_STATUS_OPTIONAL, ierr2) ierr = ierr + ierr2 call kim_set_argument_support_status( & model_compute_arguments_create_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_VIRIAL, & KIM_SUPPORT_STATUS_OPTIONAL, & ierr2) ierr = ierr + ierr2 call kim_set_argument_support_status( & model_compute_arguments_create_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_PARTICLE_VIRIAL, & KIM_SUPPORT_STATUS_OPTIONAL, & ierr2) ierr = ierr + ierr2 if (ierr /= 0_c_int) then call kim_log_entry(model_compute_arguments_create_handle, & KIM_LOG_VERBOSITY_ERROR, "Unable to register arguments support status") goto 42 end if ! register callbacks call kim_set_callback_support_status( & model_compute_arguments_create_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_DEDR_TERM, & KIM_SUPPORT_STATUS_OPTIONAL, ierr) call kim_set_callback_support_status( & model_compute_arguments_create_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_D2EDR2_TERM, & KIM_SUPPORT_STATUS_OPTIONAL, ierr2) ierr = ierr + ierr2 if (ierr /= 0_c_int) then call kim_log_entry(model_compute_arguments_create_handle, & KIM_LOG_VERBOSITY_ERROR, "Unable to register callback support status") goto 42 end if ierr = 0_c_int 42 continue return end subroutine compute_arguments_create !------------------------------------------------------------------------------- ! ! Model driver compute arguments destroy routine ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine compute_arguments_destroy(model_compute_handle, & model_compute_arguments_destroy_handle, ierr) bind(c) implicit none !-- Transferred variables type(kim_model_compute_handle_type), intent(in) :: model_compute_handle type(kim_model_compute_arguments_destroy_handle_type), intent(in) :: & model_compute_arguments_destroy_handle integer(c_int), intent(out) :: ierr ! avoid unsed dummy argument warnings if (model_compute_handle .eq. KIM_MODEL_COMPUTE_NULL_HANDLE) continue if (model_compute_arguments_destroy_handle .eq. & KIM_MODEL_COMPUTE_ARGUMENTS_DESTROY_NULL_HANDLE) continue ierr = 0_c_int return end subroutine compute_arguments_destroy end module Pair_Lennard_Jones_Truncated !------------------------------------------------------------------------------- ! ! Model driver initialization routine (REQUIRED) ! !------------------------------------------------------------------------------- ! The "recursive" keyword is added below make this routine thread safe recursive subroutine create(model_driver_create_handle, requested_length_unit, & requested_energy_unit, requested_charge_unit, requested_temperature_unit, & requested_time_unit, ierr) bind(c) use, intrinsic :: iso_c_binding use Pair_Lennard_Jones_Truncated use kim_model_driver_headers_module implicit none !-- Transferred variables type(kim_model_driver_create_handle_type), intent(inout) :: model_driver_create_handle type(kim_length_unit_type), intent(in), value :: requested_length_unit type(kim_energy_unit_type), intent(in), value :: requested_energy_unit type(kim_charge_unit_type), intent(in), value :: requested_charge_unit type(kim_temperature_unit_type), intent(in), value :: requested_temperature_unit type(kim_time_unit_type), intent(in), value :: requested_time_unit integer(c_int), intent(out) :: ierr !-- Local variables integer(c_int), parameter :: cd = c_double ! used for literal constants integer(c_int) :: number_of_parameter_files character(len=1024, kind=c_char) :: parameter_file_name integer(c_int) :: ierr2 character(len=100, kind=c_char) :: in_species real(c_double) :: in_cutoff real(c_double) :: in_eps real(c_double) :: in_sigma real(c_double) :: convert_length, convert_energy type(kim_species_name_type) species_name type(BUFFER_TYPE), pointer :: buf ! Register numbering call kim_set_model_numbering( & model_driver_create_handle, KIM_NUMBERING_ONE_BASED, ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, & KIM_LOG_VERBOSITY_ERROR, "Unable to set numbering") goto 42 end if ! Store callback pointers in KIM object call kim_set_routine_pointer( & model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE, & KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(compute), ierr) call kim_set_routine_pointer( & model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_CREATE, & KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(compute_arguments_create), ierr2) ierr = ierr + ierr2 call kim_set_routine_pointer( & model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_DESTROY, & KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(compute_arguments_destroy), ierr2) ierr = ierr + ierr2 call kim_set_routine_pointer( & model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_REFRESH, & KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(refresh), ierr2) ierr = ierr + ierr2 call kim_set_routine_pointer( & model_driver_create_handle, KIM_MODEL_ROUTINE_NAME_DESTROY, & KIM_LANGUAGE_NAME_FORTRAN, 1_c_int, c_funloc(destroy), ierr2) ierr = ierr + ierr2 if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Unable to store callback pointers") goto 42 end if ! Process parameter files call kim_get_number_of_parameter_files( & model_driver_create_handle, number_of_parameter_files) if (number_of_parameter_files .ne. 1_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Wrong number of parameter files") ierr = 1_c_int goto 42 end if ! Read in model parameters from parameter file call kim_get_parameter_file_name( & model_driver_create_handle, 1_c_int, parameter_file_name, ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Unable to get parameter file name") ierr = 1_c_int goto 42 end if allocate(buf) ! Read in model parameters from parameter file ! open(10,file=parameter_file_name,status="old") read(10,*,iostat=ierr,err=100) in_species read(10,*,iostat=ierr,err=100) in_cutoff read(10,*,iostat=ierr,err=100) in_eps read(10,*,iostat=ierr,err=100) in_sigma close(10) goto 200 100 continue ! Reading parameters failed ierr = 1_c_int call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Unable to read model parameters") goto 42 200 continue ! Convert parameters to requested units call kim_convert_unit( & KIM_LENGTH_UNIT_A, & KIM_ENERGY_UNIT_EV, & KIM_CHARGE_UNIT_E, & KIM_TEMPERATURE_UNIT_K, & KIM_TIME_UNIT_PS, & requested_length_unit, & requested_energy_unit, & requested_charge_unit, & requested_temperature_unit, & requested_time_unit, & 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, convert_length, ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Failed to convert units") goto 42 endif in_cutoff = in_cutoff*convert_length in_sigma = in_sigma*convert_length call kim_convert_unit( & KIM_LENGTH_UNIT_A, & KIM_ENERGY_UNIT_EV, & KIM_CHARGE_UNIT_E, & KIM_TEMPERATURE_UNIT_K, & KIM_TIME_UNIT_PS, & requested_length_unit, & requested_energy_unit, & requested_charge_unit, & requested_temperature_unit, & requested_time_unit, & 0.0_cd, 1.0_cd, 0.0_cd, 0.0_cd, 0.0_cd, convert_energy, ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Failed to convert units") goto 42 endif in_eps = in_eps*convert_energy ! Register units call kim_set_units( & model_driver_create_handle, & requested_length_unit, & requested_energy_unit, & KIM_CHARGE_UNIT_UNUSED, & KIM_TEMPERATURE_UNIT_UNUSED, & KIM_TIME_UNIT_UNUSED, ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Unable to set units") goto 42 end if buf%influence_distance = in_cutoff buf%cutoff(1) = in_cutoff buf%cutsq(1) = in_cutoff*in_cutoff buf%eps(1) = in_eps buf%sigma(1) = in_sigma buf%model_will_not_request_neighbors_of_noncontributing_particles(1) = 1_c_int ! Store model cutoff in KIM object call kim_set_influence_distance_pointer( & model_driver_create_handle, buf%influence_distance) call kim_set_neighbor_list_pointers( & model_driver_create_handle, 1_c_int, buf%cutoff, & buf%model_will_not_request_neighbors_of_noncontributing_particles) ! Register buffer in KIM API object call kim_set_model_buffer_pointer( & model_driver_create_handle, c_loc(buf)) ! Register species call kim_from_string(in_species, species_name) call kim_set_species_code( & model_driver_create_handle, species_name, speccode, ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Unable to set species code") goto 42 end if ! Set parameter pointers so they can be changed by the simulator call kim_set_parameter_pointer( & model_driver_create_handle, buf%cutoff, "cutoff", & "Distance beyond which particles do not interact.", ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Failed to set parameter 'cutoff'.") goto 42 endif call kim_set_parameter_pointer( & model_driver_create_handle, buf%eps, "epsilon", & "Energy scaling coefficient (4*epsilon is the depth of the potential & &well).", ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Failed to set parameter 'epsilon'.") goto 42 endif call kim_set_parameter_pointer( & model_driver_create_handle, buf%sigma, "sigma", & "Distance at which energy is zero and force is repulsive.", ierr) if (ierr /= 0_c_int) then call kim_log_entry(model_driver_create_handle, KIM_LOG_VERBOSITY_ERROR, & "Failed to set parameter 'sigma'.") goto 42 endif ierr = 0_c_int 42 continue return end subroutine create