!**************************************************************************** ! ! 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. ! !**************************************************************************** ! ! MODULE Pair_Exp6_Hogervorst ! ! ! Compute energy and forces on an isolated cluster of atoms using ! the exp-6 potential shifted to have zero energy at cutoff. ! Parameters due to Hogervorst [1] and mixing rule between species ! due to Kong and Chakrabarty [2]. ! ! Author : E. B. Tadmor (10-MAR-12) ! ! Support atom types: ! ! 1 = "Ar" ! 2 = "Ne" ! ! References: ! ! [1] W. Hogervorst, Physica, vol 51, 77-89 (1971). (Like interactions.) ! ! [2] C. L. Kong and M. R. Chakrabarty, J. Phys. Chem., Vol. 77, 2668-2670 ! (1973). (Mixing rule.) ! !**************************************************************************** #include "KIM_API_status.h" #define THIS_FILE_NAME __FILE__ #define TRUEFALSE(TRUTH) merge(1,0,(TRUTH)) module Pair_Exp6_Hogervorst use KIM_API implicit none save private public Compute_Energy_Forces,cutoff ! ! Define potential parameters ! double precision, parameter :: kb = 8.6173e-5 ! Boltzmann's constant [eV/K] double precision, parameter :: cutoff = 8.15d0 ! cutoff radius [A] (arbitrary value) double precision, parameter :: eps11 = 138.d0*kb ! Ar-Ar epsilon parameter [eV] double precision, parameter :: rm11 = 3.77d0 ! Ar-Ar rm parameter [A] double precision, parameter :: alf11 = 14.8d0 ! Ar-Ar alpha parameter [-] double precision, parameter :: eps22 = 43.d0*kb ! Ne-Ne epsilon parameter [eV] double precision, parameter :: rm22 = 3.03d0 ! Ne-Ne rm parameter [A] double precision, parameter :: alf22 = 16.d0 ! Ne-Ne alpha parameter [-] double precision, parameter :: eps12 = 68.89*kb ! Ar-Ne epsilon parameter [eV] double precision, parameter :: rm12 = 3.447d0 ! Ar-Ne rm parameter [A] double precision, parameter :: alf12 = 15.52d0 ! Ar-Ne alpha parameter [-] contains !------------------------------------------------------------------------------- ! ! Calculate exp-6 pair potential phi(r) and its derivative dphi(r) ! Note: potential is shifted to have zero energy at r=cutoff ! !------------------------------------------------------------------------------- subroutine calc_phi_dphi(r,phi,dphi,eps,rm,alf,cutoff) implicit none !-- Transferred variables double precision, intent(in) :: r double precision, intent(in) :: eps,rm,alf,cutoff double precision, intent(out) :: phi,dphi !-- Local variables double precision soa,amp,ror,cor,phicut if (r .gt. cutoff) then ! Argument exceeds cutoff radius phi = 0.d0 dphi = 0.d0 else soa = 6.d0/alf amp = eps/(1.d0-soa) cor = cutoff/rm ror = r/rm phicut = amp*(soa*exp(alf*(1.d0-cor))-cor**(-6)) phi = amp*(soa*exp(alf*(1.d0-ror))-ror**(-6)) - phicut dphi = -6.d0*amp/rm*(exp(alf*(1.d0-ror))-ror**(-7)) endif return end subroutine calc_phi_dphi !------------------------------------------------------------------------------- ! ! Compute energy and forces on atoms from the positions. ! !------------------------------------------------------------------------------- integer function Compute_Energy_Forces(pkim) implicit none !--Transferred Variables integer(kind=kim_intptr), intent(in) :: pkim !--Local Variables integer ier integer i,j double precision Rij(3),Rsqij,cutoffsq,r,phi,dphi !-- KIM variables integer N; pointer(pN,N) double precision energy; pointer(penergy,energy) double precision coordum(3,1); pointer(pcoor,coordum) double precision forcedum(3,1); pointer(pforce,forcedum) integer particleTypes(1); pointer(pparticleTypes,particleTypes) double precision, pointer :: coors(:,:),forces(:,:) integer, pointer :: types(:) integer idum integer comp_energy,comp_force ! ! Check to see if we have been asked to compute the energy and forces ! call kim_api_getm_compute_f(pkim, ier, & "energy", comp_energy, 1, & "forces", comp_force, 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 ! ! Unpack data from KIM object ! call kim_api_getm_data_f(pkim, ier, & "numberOfParticles", pN, 1, & "particleTypes", pparticleTypes, 1, & "coordinates", pcoor, 1, & "energy", penergy, TRUEFALSE(comp_energy.eq.1), & "forces", pforce, TRUEFALSE(comp_force.eq.1)) 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_real_array_2d(coordum,coors,3,N) if (comp_force.eq.1) call KIM_to_F90_real_array_2d(forcedum,forces,3,N) call KIM_to_F90_int_array_1d(particleTypes,types,N) ! ! Verify input ! ier = KIM_STATUS_OK ! assume all is OK do i=1,N if (types(i).ne.1 .and. types(i).ne.2) ier = KIM_STATUS_FAIL enddo if (ier.ne.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "Unexpected species type detected", ier) goto 42 endif ! ! Compute energy and forces (use half-list approach) ! if (comp_energy.eq.1) energy = 0.d0 if (comp_force.eq.1) forces = 0.d0 cutoffsq = cutoff**2 do i=1,N-1 do j=i+1,N Rij(:) = coors(:,j) - coors(:,i) Rsqij = Rij(1)*Rij(1) + Rij(2)*Rij(2) + Rij(3)*Rij(3) if ( Rsqij.lt.cutoffsq ) then r = sqrt(Rsqij) if (types(i).eq.1 .and. types(j).eq.1) then call calc_phi_dphi(r,phi,dphi,eps11,rm11,alf11,cutoff) elseif (types(i).eq.2 .and. types(j).eq.2) then call calc_phi_dphi(r,phi,dphi,eps22,rm22,alf22,cutoff) else call calc_phi_dphi(r,phi,dphi,eps12,rm12,alf12,cutoff) endif if (comp_energy.eq.1) energy = energy + phi if (comp_force.eq.1) then forces(:,i) = forces(:,i) + dphi*Rij/r forces(:,j) = forces(:,j) - dphi*Rij/r endif endif enddo enddo ier = KIM_STATUS_OK 42 continue Compute_Energy_Forces = ier return end function Compute_Energy_Forces end module Pair_Exp6_Hogervorst !------------------------------------------------------------------------------- ! ! Model initialization routine (REQUIRED) ! !------------------------------------------------------------------------------- integer function model_init(pkim) use Pair_Exp6_Hogervorst use KIM_API implicit none !-- Transferred variables integer(kind=kim_intptr), intent(in) :: pkim !-- Local variables integer(kind=kim_intptr), parameter :: one=1 integer ier, idum !-- KIM variables double precision rcut; pointer(prcut,rcut) ! store pointer to compute function in KIM object ier = kim_api_set_method_f(pkim,"compute",one,loc(Compute_Energy_Forces)) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_set_method", ier) goto 42 endif ! store model cutoff in KIM object prcut = 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) goto 42 endif rcut = cutoff ier = KIM_STATUS_OK 42 continue model_init = ier return end function model_init