module error use, intrinsic :: iso_c_binding implicit none public contains recursive subroutine my_error(message) implicit none character(len=*, kind=c_char), intent(in) :: message print *,"* Error : ", trim(message) stop 1 end subroutine my_error recursive subroutine my_warning(message) implicit none character(len=*,kind=c_char),intent(in) :: message print *,"* Warning : ", trim(message) end subroutine my_warning end module error module mod_kim use, intrinsic :: iso_c_binding use error use mod_global use kim_simulator_headers_module use kim_supported_extensions_module implicit none integer(c_int), parameter :: SizeOne = 1 ! ! KIM variables ! integer(c_int) :: ierr, idum,ierr2 integer(c_int) :: requested_units_accepted,number_of_model_routine_names type(kim_model_routine_name_type) model_routine_name integer(c_int) :: present integer(c_int) :: required integer(c_int) :: model_will_not_request_neighbors_of_noncontributing_particles(1) type(kim_compute_arguments_handle_type) :: compute_arguments_handle type(kim_model_handle_type) :: model_handle type(c_ptr) pDummy contains subroutine initialize_kim(modelname) implicit none ! transferred variables character(len=256, kind=c_char) :: modelname ! local variables integer i call kim_model_create(KIM_NUMBERING_ONE_BASED, & KIM_LENGTH_UNIT_A, & KIM_ENERGY_UNIT_EV, & KIM_CHARGE_UNIT_E, & KIM_TEMPERATURE_UNIT_K, & KIM_TIME_UNIT_PS, & trim(modelname), & requested_units_accepted, & model_handle, ierr) if (ierr /= 0) then call my_error("kim_api_create") endif ! check that units of test and model are compatible if (requested_units_accepted == 0) then call my_error("Must adapt to model units") end if ! Get the number of model routine names call kim_get_number_of_model_routine_names(number_of_model_routine_names) do i=1,number_of_model_routine_names call kim_get_model_routine_name(i, model_routine_name, ierr) if (ierr /= 0) then call my_error("kim_get_model_routine_name") endif call kim_is_routine_present(model_handle, model_routine_name, present, & required, ierr) if (ierr /= 0) then call my_error("kim_is_routine_present") endif if ((present == 1) .and. (required == 1)) then if (.not. ((model_routine_name == KIM_MODEL_ROUTINE_NAME_CREATE) & .or. (model_routine_name == & KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_CREATE) & .or. (model_routine_name == KIM_MODEL_ROUTINE_NAME_COMPUTE) & .or. (model_routine_name == KIM_MODEL_ROUTINE_NAME_REFRESH) & .or. (model_routine_name == & KIM_MODEL_ROUTINE_NAME_COMPUTE_ARGUMENTS_DESTROY) & .or. (model_routine_name == KIM_MODEL_ROUTINE_NAME_DESTROY))) then call my_error("Unknown required ModelRoutineName found.") endif endif enddo ! create compute_arguments object call kim_compute_arguments_create(model_handle, compute_arguments_handle, ierr) if (ierr /= 0) then call my_error("kim_model_compute_arguments_create") endif end subroutine initialize_kim function getSpeciesCode_kim(species) implicit none ! transferref variables character(len=2),intent(in) :: species integer(c_int) :: getSpeciesCode_kim ! local variables integer(c_int) :: species_is_supported, speciesCode type(kim_species_name_type) :: kim_species call kim_from_string(species,kim_species) call kim_get_species_support_and_code(model_handle, & kim_species, species_is_supported, speciesCode, ierr) if ((ierr /= 0) .or. (species_is_supported /= 1)) then call my_error("Model does not support " // species) endif getSpeciesCode_kim = speciesCode return end function getSpeciesCode_kim function queryAnalytical_kim() implicit none ! transferred variables integer :: queryAnalytical_kim ! local variables type(kim_support_status_type) :: processdEdr_support_status, processd2Edr2_support_status call kim_get_callback_support_status(compute_arguments_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_DEDR_TERM,& processdEdr_support_status,ierr) if ((ierr /= 0)) then call my_warning("can't get support status of processdEdr") endif call kim_get_callback_support_status(compute_arguments_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_D2EDR2_TERM, & processd2Edr2_support_status,ierr) if ((ierr /= 0)) then call my_warning("can't get support status of processd2Edr2") endif if ( (processdEdr_support_status == KIM_SUPPORT_STATUS_NOT_SUPPORTED) .or. & (processd2Edr2_support_status == KIM_SUPPORT_STATUS_NOT_SUPPORTED) ) then queryAnalytical_kim = 0 return else queryAnalytical_kim = 1 return endif end function queryAnalytical_kim function getCutoff_kim() implicit none ! transferred variables real(c_double) :: getCutoff_kim !-- Local variables real(c_double) :: cutoffs(1),rcutoff ! Get the cutoff radius from the model call kim_get_neighbor_list_values(model_handle, cutoffs, & model_will_not_request_neighbors_of_noncontributing_particles, ierr) if (ierr /= 0) then call my_error("get_neighbor_list_values") end if call kim_get_influence_distance(model_handle, rcutoff) rcutoff = cutoffs(1) getCutoff_kim = rcutoff end function getCutoff_kim subroutine broadcastAtomistic_kim(aS) use mod_atomistic, only : Atomistic, process_dedr,process_d2edr2 implicit none ! Passed variables type(Atomistic),pointer :: aS ! local variables integer i type(kim_species_name_type) species_name integer(c_int) :: number_of_neighbor_lists ierr = 0 ! broadcast number of particles call kim_set_argument_pointer(compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_NUMBER_OF_PARTICLES, as%numParticles, ierr2) ierr = ierr + ierr2 ! broadcast particleSpecies call kim_set_argument_pointer(compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_SPECIES_CODES, aS%particleSpecies, & ierr2) ierr = ierr + ierr2 ! broadcast particle_contributing flags call kim_set_argument_pointer(compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTICLE_CONTRIBUTING, aS%particle_contributing,ierr2) ierr = ierr + ierr2 ! broadcast coordinates call kim_set_argument_pointer(compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_COORDINATES, aS%coords, ierr2) ierr = ierr + ierr2 ! broadcase partial energy call kim_set_argument_pointer(compute_arguments_handle, & KIM_COMPUTE_ARGUMENT_NAME_PARTIAL_ENERGY, aS%energy, ierr2) ierr = ierr + ierr2 ! Set pointer in KIM object to neighbor list routine and object ! call kim_set_callback_pointer(compute_arguments_handle, & KIM_COMPUTE_CALLBACK_NAME_GET_NEIGHBOR_LIST, KIM_LANGUAGE_NAME_FORTRAN, & c_funloc(get_neigh), c_loc(aS%neighObject), ierr) if (ierr /= 0) then call my_error("set_callback_pointer get_neigh") end if ! process_dedr if (aS%processdEdrFlag .eq. 1) then call kim_set_callback_pointer(compute_arguments_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_DEDR_TERM, KIM_LANGUAGE_NAME_FORTRAN, & c_funloc(process_dedr), pDummy, ierr) if (ierr /= 0) then call my_error("set_callback_pointer process_dedr") end if endif if (aS%processd2Edr2Flag .eq. 1) then ! process_d2edr2 call kim_set_callback_pointer(compute_arguments_handle, & KIM_COMPUTE_CALLBACK_NAME_PROCESS_D2EDR2_TERM, KIM_LANGUAGE_NAME_FORTRAN, & c_funloc(process_d2edr2), pDummy, ierr) if (ierr /= 0) then call my_error("set_callback_pointer process_d2edr2") end if end if call kim_get_number_of_neighbor_lists(model_handle, & number_of_neighbor_lists) if (number_of_neighbor_lists /= 1) then call my_error("too many neighbor lists") endif return end subroutine broadcastAtomistic_kim subroutine get_neigh(data_object, number_of_neighbor_lists, cutoffs, & neighbor_list_index, request, numnei, pnei1part, ier) bind(c) use mod_atomistic, only : NeighObject_type implicit none !-- Transferred variables type(c_ptr), value, intent(in) :: data_object integer(c_int), value, intent(in) :: number_of_neighbor_lists real(c_double), intent(in) :: cutoffs(number_of_neighbor_lists) integer(c_int), value, intent(in) :: neighbor_list_index integer(c_int), value, intent(in) :: request integer(c_int), intent(out) :: numnei type(c_ptr), intent(out) :: pnei1part integer(c_int), intent(out) :: ier !-- Local variables type(NeighObject_type), pointer :: neighObject !integer(c_int), pointer :: nat; type(c_ptr) :: pnat call c_f_pointer(data_object, neighObject) ! unpack number of particles !call kim_get_argument_pointer(compute_arguments_handle,KIM_COMPUTE_ARGUMENT_NAME_NUMBER_OF_PARTICLES,pnat,ierr2) !call c_f_pointer(pnat, nat) if (number_of_neighbor_lists /= 1) then call my_error("invalid number of cutoffs") ier = 1 return endif if (neighbor_list_index /= 1) then call my_error("wrong list index") ier = 1 return endif !if ( (request.gt.nat) .or. (request.lt.1)) then ! call my_error("Invalid part ID in get_neigh") ! ier = 1 ! return !endif ! set the returned number of neighbors for the returned part numnei = neighObject%neighborList(1,request) ! set the location for the returned neighbor list pnei1part = c_loc(neighObject%neighborList(2,request)) ier = 0 return end subroutine get_neigh integer function compute_kim() implicit none call kim_compute(model_handle, compute_arguments_handle, ierr) compute_kim = ierr end function compute_kim subroutine free_KIM_API_object !------------------------------------------------------------------------------- ! ! free_KIM_API_object : Deallocate storage and destroy KIM API object ! !------------------------------------------------------------------------------- implicit none call kim_compute_arguments_destroy(& model_handle, compute_arguments_handle, ierr) if (ierr /= 0) then call my_error("compute_arguments_destroy") endif call kim_model_destroy(model_handle) return end subroutine free_KIM_API_object subroutine write_results_kim(species,cS,cT,fsgeType) use mod_fsge, only : voigtStrainGradient,voigtStrain,Fsge use mod_crystal, only : Crystal implicit none ! Passed variables character(len=*),intent(in) :: cS type(Crystal),intent(in) :: cT type(Fsge),intent(in) :: fsgeType character(len=2),intent(in) :: species ! Local variables integer(c_int) :: runit,m,n,i integer,external :: GetFreeUnit character(len=1) :: cr character(len=30) :: cvalstr,dvalstr character(len=2) :: mstr,nstr real(c_double) :: a,c integer(c_int) :: cpntsD(6),cpntsC(4) character(len=2000):: basisAtomsCoordsFractionString, & latticeConstantString, & crystalString, & speciesString, & CString, & DString,& propertyString,& tempstr cr = char(10) crystalString = ' "' // trim(cS) //'"' // cr if (trim(cS) .eq. 'fcc') then a = cT%LatticeType%latVec(2,1)*2._cd write(basisAtomsCoordsFractionString,*) & ' [' // cr // & ' 0' // cr // & ' 0' // cr // & ' 0' // cr // & ' ]' // cr // & ' [' // cr // & ' 0' // cr // & ' 0.5' // cr // & ' 0.5' // cr // & ' ]' // cr // & ' [' // cr // & ' 0.5 ' // cr // & ' 0 ' // cr // & ' 0.5' // cr // & ' ]' // cr // & ' [' // cr // & ' 0.5' // cr // & ' 0.5' // cr // & ' 0' // cr // & ' ]' // cr write(speciesString,*) & ' "' // species //'"' // cr // & ' "' // species //'"' // cr // & ' "' // species //'"' // cr // & ' "' // species //'"' // cr write(propertyString,*) & ' "property-id" "tag:staff@noreply.openkim.org,2016-05-24:property/' propertyString = trim(adjustl(propertyString)) // & 'elastic-constants-first-strain-gradient-isothermal-cubic-crystal-npt"' // cr elseif(trim(cS) .eq. 'bcc') then a = cT%LatticeType%latVec(2,1)*2._cd write(basisAtomsCoordsFractionString,*) & ' [' // cr // & ' 0' // cr // & ' 0' // cr // & ' 0' // cr // & ' ]' // cr // & ' [' // cr // & ' 0.5' // cr // & ' 0.5' // cr // & ' 0.5' // cr // & ' ]' // cr write(speciesString,*) & ' "' // species //'"' // cr // & ' "' // species //'"' // cr write(propertyString,*) & ' "property-id" "tag:staff@noreply.openkim.org,2016-05-24:property/' propertyString = trim(adjustl(propertyString)) // & 'elastic-constants-first-strain-gradient-isothermal-cubic-crystal-npt"' // cr elseif(trim(cS) .eq. 'hex') then c = cT%LatticeType%latVec(3,3) write(basisAtomsCoordsFractionString,*) & ' [' // cr // & ' 0' // cr // & ' 0' // cr // & ' 0' // cr // & ' ]' // cr write(speciesString,*) & ' "' // species //'"' // cr write(propertyString,*) & ' "property-id" "tag:staff@noreply.openkim.org,2016-05-24:property/' propertyString = trim(adjustl(propertyString)) // & 'elastic-constants-first-strain-gradient-isothermal-monoatomic-hexagonal-crystal-npt"' // cr endif write(tempstr,*) a write(latticeConstantString,*) & ' "a" {' // cr // & ' "source-value" '//trim(adjustl(tempstr)) // cr // & ' "source-unit" "angstrom"' // cr // & ' }' // cr if (trim(cS) .eq. 'hex' .or. trim(cS) .eq. 'hcp') then write(tempstr,*) c latticeConstantString = trim(latticeConstantString) // & ' "c" {' // cr // & ' "source-value" '//trim(adjustl(tempstr)) // cr // & ' "source-unit" "angstrom"' // cr // & ' }' // cr endif CString='' do i=1,size(fsgeType%IndComponentsC,1) m = fsgeType%IndComponentsC(i,1); write(mstr,'(i2)') m n = fsgeType%IndComponentsC(i,2); write(nstr,'(i2)') n cpntsC(1:2) = voigtStrain(m) cpntsC(3:4) = voigtStrain(n) write(cvalstr,*) fsgeType%Cijkl(cpntsC(1),cpntsC(2),cpntsC(3),cpntsC(4)) write(tempstr,*)& ' "c' //trim(adjustL(mstr)) // trim(adjustl(nstr)) // '" {' // cr // & ' "source-value" ' // trim(cvalstr) // cr // & ' "source-unit" "eV/angstrom^3"' // cr // & ' }' // cr CString = trim(CString) // trim(tempstr) enddo DString='' do i=1,size(fsgeType%IndComponentsD,1) m = fsgeType%IndComponentsD(i,1); write(mstr,'(i2)') m n = fsgeType%IndComponentsD(i,2); write(nstr,'(i2)') n cpntsD(1:3) = voigtStrainGradient(m) cpntsD(4:6) = voigtStrainGradient(n) write(dvalstr,'(F20.8)') fsgeType%Dijmkln(cpntsD(1),cpntsD(2),cpntsD(3),cpntsD(4),cpntsD(5),cpntsD(6)) write(tempstr,*) & ' "d-' // trim(adjustl(mstr)) // '-' // trim(adjustl(nstr)) // '" {' // cr // & ' "source-value" ' // trim(adjustl(dvalstr)) // cr // & ' "source-unit" "eV/angstrom"' // cr // & ' }' // cr DString = trim(Dstring) // trim(tempstr) enddo ! Now write the results.edn file runit = GetFreeUnit() open(unit=runit,file='output/results.edn',status='replace') write(runit,*) & '{' // cr // & !!!!!!!!!!!!!!!! ' "short-name" {' // cr // & ' "source-value" [' // cr // & trim(crystalString) // & ' ]' // cr // & ' }' // cr // & !!!!!!!!!!!!!!!! trim(latticeConstantString) // & !!!!!!!!!!!!!!!! trim(propertyString) // & !!!!!!!!!!!!!!!! ' "basis-atom-coordinates" {' // cr // & ' "source-value" [' // cr // & trim(basisAtomsCoordsFractionString) // & ' ]' // cr // & ' }' // cr // & !!!!!!!!!!!!!!!! ' "species" {' // cr // & ' "source-value" [' // cr // & trim(speciesString) // & ' ]' // cr // & ' }' // cr // & !!!!!!!!!!!!!!!! ' "temperature" {' // cr // & ' "source-value" 0' // cr // & ' "source-unit" "K"' // cr // & ' }' // cr // & !!!!!!!!!!!!!!!! trim(CString) // & !!!!!!!!!!!!!!!! trim(DString) // & !!!!!!!!!!!!!!!! ' "instance-id" 1' // cr // & '}' end subroutine write_results_kim end module mod_kim integer function GetFreeUnit() implicit none logical InUse do GetFreeUnit=7,98 inquire(unit=GetFreeUnit,opened=InUse) if (.not. InUse) return enddo write(*,*) "Could not obtain a free unit handle in function GetFreeUnit" stop end function GetFreeUnit