! 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) 2012, Regents of the University of Minnesota. All rights reserved. ! ! Contributors: ! Ryan S. Elliott ! Ellad B. Tadmor ! Valeriu Smirichinski ! Feng Gao ! Eunkoo Lee ! Bohumir Jelinek ! !**************************************************************************** !** !** MODULE kcc_meam !** !** kcc_meam pair potential model for !** !** Reference: !** !** Language: Fortran 90 !** !** Release: !** !**************************************************************************** #include "KIM_API_status.h" #define THIS_FILE_NAME __FILE__ #define TRUEFALSE(TRUTH) merge(1,0,(TRUTH)) ! ! All MEAM parameters stored in this module ! module meam_parameters implicit none save public type meam_params integer :: max_comp integer :: LineTable ! lj_1 double precision, pointer :: PhiTab(:,:,:) double precision, pointer :: DPhiTab(:,:,:) ! lj_2 double precision :: Rmin double precision :: Rsqmin double precision :: Rcutoff double precision :: OverDeltaRsq ! eam1 double precision, pointer :: Ec(:,:) double precision, pointer :: Re(:) double precision, pointer :: B(:,:) double precision, pointer :: A(:) double precision, pointer :: beta(:,:) double precision, pointer :: t(:,:) ! eam2 double precision, pointer :: RhoZ(:) double precision, pointer :: Rho0_Bar(:) double precision, pointer :: Z(:) double precision, pointer :: Z2(:) double precision, pointer :: A2(:) ! eam3 double precision, pointer :: Cmin(:,:,:) double precision, pointer :: Cmax(:,:,:) double precision :: Dcutoff double precision :: OverDeltaR double precision, pointer :: Eu_d(:) ! eam4 double precision, pointer :: Re_Alloy(:,:) double precision, pointer :: alpha(:,:) double precision, pointer :: omega(:,:) ! ener (only amass) double precision, pointer :: Amass(:) ! scrn double precision :: scr_2nn end type meam_params type (meam_params) mp; pointer (a_mp,mp) contains subroutine allocate_meam_parameters(max_comp,LineTable) implicit none !--Transferred variables integer, intent(in) :: max_comp integer, intent(in) :: LineTable mp%max_comp = max_comp mp%LineTable = LineTable allocate(mp%PhiTab(LineTable+1,max_comp,max_comp)) allocate(mp%DPhiTab(LineTable+1,max_comp,max_comp)) allocate(mp%Ec(max_comp,max_comp)) allocate(mp%Re(max_comp)) allocate(mp%B(max_comp,max_comp)) allocate(mp%A(max_comp)) allocate(mp%beta(max_comp,4)) allocate(mp%t(max_comp,3)) allocate(mp%RhoZ(max_comp)) allocate(mp%Rho0_Bar(max_comp)) allocate(mp%Z(max_comp)) allocate(mp%Z2(max_comp)) allocate(mp%A2(max_comp)) allocate(mp%Cmin(max_comp,max_comp,max_comp)) allocate(mp%Cmax(max_comp,max_comp,max_comp)) allocate(mp%Eu_d(max_comp)) allocate(mp%Re_Alloy(max_comp,max_comp)) allocate(mp%alpha(max_comp,max_comp)) allocate(mp%omega(max_comp,max_comp)) allocate(mp%Amass(max_comp)) end subroutine allocate_meam_parameters subroutine deallocate_meam_parameters implicit none deallocate(mp%PhiTab) deallocate(mp%DPhiTab) deallocate(mp%Ec) deallocate(mp%Re) deallocate(mp%B) deallocate(mp%A) deallocate(mp%beta) deallocate(mp%t) deallocate(mp%RhoZ) deallocate(mp%Rho0_Bar) deallocate(mp%Z) deallocate(mp%Z2) deallocate(mp%A2) deallocate(mp%Cmin) deallocate(mp%Cmax) deallocate(mp%Eu_d) deallocate(mp%Re_Alloy) deallocate(mp%alpha) deallocate(mp%omega) deallocate(mp%Amass) end subroutine deallocate_meam_parameters end module meam_parameters module kcc_meam use KIM_API implicit none save private public Compute_Energy_Forces, & ReInit, & Destroy, & Ncomp, & Icomp ! Below are the definitions and values of all Model parameters integer, parameter :: DIM=3 ! dimensionality of space integer Ncomp character*2, allocatable :: Icomp(:) contains !------------------------------------------------------------------------------- ! ! Compute energy and forces on atoms from the positions. ! !------------------------------------------------------------------------------- integer function Compute_Energy_Forces(pkim) use meam_parameters implicit none !-- Transferred variables integer(kind=kim_intptr), intent(in) :: pkim !-- Local variables integer ier integer i,idum !-- KIM variables character*64 NBC_Method; pointer(pNBC_Method,NBC_Method) integer HalfOrFull ! 1 - half, 2 - full integer numberOfParticles; pointer(pnAtoms,numberOfParticles) integer nparticleTypes; pointer(pnparticleTypes,nparticleTypes) integer ptypesdum(1); pointer(pptypes,ptypesdum) double precision cutoff; pointer(pcutoff,cutoff) double precision energy; pointer(penergy,energy) double precision coordum(DIM,1); pointer(pcoor,coordum) double precision forcedum(DIM,1); pointer(pforce,forcedum) double precision enepotdum(1); pointer(penepot,enepotdum) double precision virialdum(1); pointer(pvirial,virialdum) double precision, pointer :: coor(:,:),force(:,:),ene_pot(:),virial_global(:) integer, pointer :: particleTypes(:) integer comp_energy, comp_force, comp_enepot, comp_virial double precision boxSideLengths(3); pointer(pboxSideLengths,boxSideLengths) integer numContrib; pointer(pnumContrib,numContrib) integer nei1atom(1); pointer(pnei1atom,nei1atom) double precision Rij_dummy(3,1); pointer(pRij_dummy,Rij_dummy) ! Get model buffer from KIM object ! a_mp = kim_api_get_model_buffer_f(pkim, ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_model_buffer_f", ier) goto 42 endif ! Determine whether half or full lists are being used ! pNBC_Method = kim_api_get_nbc_method_f(pkim, ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_nbc_method_f", ier) goto 42 endif HalfOrFull = 1 if (index(NBC_Method,"MI_OPBC_F").eq.1) then HalfOrFull = 2 else ier = KIM_STATUS_FAIL idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "Unsupported NBC type", ier) goto 42 endif call free(pNBC_Method) ! don't forget to release the memory... ! Check to see if we have been asked to compute the energy, forces, energyperatom, ! and virial ! call kim_api_getm_compute_f(pkim, ier, & "energy", comp_energy, 1, & "forces", comp_force, 1, & "particleEnergy", comp_enepot, 1, & "virial", comp_virial, 1) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_getm_compute_f", ier) goto 42 endif if (comp_enepot.eq.1) then ier = KIM_STATUS_FAIL idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "ERROR: Energy per atom not supported by MEAM yet.", ier) goto 42 endif ! Unpack data from KIM object ! call kim_api_getm_data_f(pkim, ier, & "numberOfParticles", pnAtoms, 1, & "numberContributingParticles", pnumContrib, TRUEFALSE(HalfOrFull.eq.1), & "numberParticleTypes", pnparticleTypes, 1, & "particleTypes", pptypes, 1, & "cutoff", pcutoff, 1, & "coordinates", pcoor, 1, & "boxSideLengths", pboxSideLengths, 1, & "energy", penergy, comp_energy, & "forces", pforce, comp_force, & "particleEnergy", penepot, comp_enepot, & "virial", pvirial, comp_virial) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_getm_data_f", ier) goto 42 endif call KIM_to_F90_int_array_1d(ptypesdum,particleTypes,numberOfParticles) call KIM_to_F90_real_array_2d(coordum,coor,DIM,numberOfParticles) if (comp_force.eq.1) & call KIM_to_F90_real_array_2d(forcedum,force,DIM,numberOfParticles) if (comp_enepot.eq.1) & call KIM_to_F90_real_array_1d(enepotdum,ene_pot,numberOfParticles) if (comp_virial.eq.1) & call KIM_to_F90_real_array_1d(virialdum,virial_global,6) ! Check to be sure that the atom types are correct by comparing ! the provided species codes to the value given here (which should ! be the same as that given in the .kim file). ! ier = KIM_STATUS_FAIL ! assume an error do i = 1,numberOfParticles if (particleTypes(i) .gt. Ncomp) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "Wrong Atom Type", ier) goto 42 endif enddo ier = KIM_STATUS_OK ! everything is ok ! Initialize potential energies, forces, virial term ! if (comp_enepot.eq.1) ene_pot(1:numberOfParticles) = 0.d0 if (comp_energy.eq.1) energy = 0.d0 if (comp_force.eq.1) force(1:3,1:numberOfParticles) = 0.d0 if (comp_virial.eq.1) virial_global = 0.d0 ! Call MEAM potential ! !@ NEED TO ADD ENERGY PER ATOM !@ call MEAM_potential(pkim,numberOfParticles,coor,particleTypes,boxSideLengths, & energy, force, virial_global, & comp_energy, comp_force, comp_virial, ier) ier = KIM_STATUS_OK 42 continue Compute_Energy_Forces = ier return end function Compute_Energy_Forces !------------------------------------------------------------------------------- ! ! Model reinitialization routine ! !------------------------------------------------------------------------------- integer function ReInit(pkim) implicit none !-- Transferred variables integer(kind=kim_intptr), intent(in) :: pkim ! assume all is well ! ReInit = KIM_STATUS_OK end function ReInit !------------------------------------------------------------------------------- ! ! Model destroy routine ! !------------------------------------------------------------------------------- integer function Destroy(pkim) use KIM_API use meam_parameters implicit none !-- Transferred variables integer(kind=kim_intptr), intent(in) :: pkim !-- Local variables integer ier integer idum ! get model buffer from KIM object a_mp = kim_api_get_model_buffer_f(pkim, ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_model_buffer_f", ier) goto 42 endif ! free storage deallocate(Icomp) call deallocate_meam_parameters call free(a_mp) ier = KIM_STATUS_OK 42 continue Destroy = ier return end function Destroy end module kcc_meam !------------------------------------------------------------------------------- ! ! Model initialization routine (REQUIRED) ! !------------------------------------------------------------------------------- integer function model_driver_init(pkim, byte_paramfile, nmstrlen, numparamfiles) use kcc_meam use meam_parameters use KIM_API implicit none !-- Transferred variables integer(kind=kim_intptr), intent(in) :: pkim integer, intent(in) :: nmstrlen integer, intent(in) :: numparamfiles byte, intent(in) :: byte_paramfile(nmstrlen*numparamfiles) !-- Local variables integer(kind=kim_intptr), parameter :: one=1 character(len=nmstrlen) paramfile_names(numparamfiles) integer ier, idum, return_error integer i !-- KIM variables double precision cutoff; pointer(pcutoff,cutoff) character(len=KIM_KEY_STRING_LENGTH) types(1); pointer(ptypes,types) integer num_types type(meam_params), pointer :: temp_mp ! assume all is well ! return_error = KIM_STATUS_OK ! generic code to process model parameter file names from byte string ! do i=0,numparamfiles-1 write(paramfile_names(i+1),'(1000a)') & char(byte_paramfile(i*nmstrlen+1: & i*nmstrlen+minloc(abs(byte_paramfile),dim=1)-1)) enddo ! store pointers in KIM object ! call kim_api_setm_method_f(pkim, ier, & "compute", one, loc(Compute_Energy_Forces), 1, & "reinit", one, loc(ReInit), 1, & "destroy", one, loc(Destroy), 1) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_setm_method_f", ier) return_error = ier goto 1000 endif ! get number and names of species from KIM API object ! ptypes = kim_api_get_test_partcl_typs_f(pkim,num_types,ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_partcl_typs_f",ier) return_error = ier goto 1000 endif Ncomp=num_types allocate(Icomp(Ncomp)) do i=1,Ncomp Icomp(i) = types(i)(1:index(types(i),char(0))-1) ! strip off null character !@ print *,i,Icomp(i) enddo ! which signals end of string ! set codes for species read in do i=1,Ncomp call kim_api_set_partcl_type_code_f(pkim,types(i),i,ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_set_partcl_type_code_f",ier) return_error = ier goto 1000 endif enddo ! in API ! allocate storage for MEAM parameter tables ! allocate(temp_mp) a_mp = loc(temp_mp) call allocate_meam_parameters(Ncomp,350001) ! Read in model parameters from parameter file ! call Define_EAM_Potential(paramfile_names(1),nmstrlen,Ncomp,Icomp,ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "Unable to read MEAM parameters, kimerror = ",KIM_STATUS_FAIL) return_error = ier goto 1000 endif ! store pointer to parameters in KIM API object ! call kim_api_set_model_buffer_f(pkim,loc(mp),ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_set_model_buffer_f",ier) return_error = ier goto 1000 endif ! store model cutoff in KIM object pcutoff = kim_api_get_data_f(pkim,"cutoff",ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_data", ier) return_error = ier goto 1000 endif cutoff = mp%Rcutoff 1000 continue model_driver_init = return_error return end function model_driver_init