! ! 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 ! ! !**************************************************************************** !** !** MODULE model__P_ !** !** 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)) module model_He_P_AbIn use KIM_API implicit none save private public Compute_Energy_Forces, & model_cutoff ! Below are the definitions and values of all Model parameters integer, parameter :: DIM=3 ! dimensionality of space integer, parameter :: speccode = 1 ! internal species code double precision, parameter :: model_cutoff = 8.0D0 ! cutoff radius ! in angstroms double precision, parameter :: empir_cutoff = 0.522562533D0 ! cutoff radius for ! empirical/ab initio switch in angstroms double precision, parameter :: model_cutsq = model_cutoff**2 !------------------------------------------------------------------------------- ! Below are the definitions and values of all additional model parameters ! ! Recall that the Fortran 90 format for declaring parameters is as follows: ! ! integer, parameter :: parname = value ! This defines an integer ! ! parameter called `parname' with a ! ! value equal to `value' (a number) ! ! double precision, parameter :: parname = value ! This defines a real double precision ! ! parameter called `parname' with a ! ! value equal to `value' (a number) !------------------------------------------------------------------------------- DOUBLE PRECISION,PARAMETER :: TOA = 0.529177249 ! Bohr to Angstrom DOUBLE PRECISION,PARAMETER :: A = 0.307092338615e7 ! Units K DOUBLE PRECISION,PARAMETER :: a1 = -0.201651289932e1/TOA ! Units 1/a0 DOUBLE PRECISION,PARAMETER :: a_1 = -0.431646276045e0*TOA ! Units a0 DOUBLE PRECISION,PARAMETER :: a2 = -0.459521265125e-1/(TOA*TOA) ! Units 1/a0/a0 DOUBLE PRECISION,PARAMETER :: a_2 = 0.138539045980e0*TOA*TOA ! Units a0*a0 DOUBLE PRECISION,PARAMETER :: d1 = 0.167127323768e-2 ! Dimensionless DOUBLE PRECISION,PARAMETER :: d2 = 0.178284243205e1/TOA ! Units 1/a0 DOUBLE PRECISION,PARAMETER :: d3 = 0.176635702255e1 ! Dimensionless DOUBLE PRECISION,PARAMETER :: b = 0.203625105759e1/TOA ! Units 1/a0 DOUBLE PRECISION,PARAMETER :: ep_kb = 10.997898 ! Units K DOUBLE PRECISION,PARAMETER :: R_ep = 5.608068*TOA ! Units a0 DOUBLE PRECISION,PARAMETER :: sig = 4.990672*TOA ! Units a0 DOUBLE PRECISION,PARAMETER :: TOK = 315776.177845143 ! AU to K DOUBLE PRECISION,PARAMETER :: TOEV = 27.21138386 ! AU to eV DOUBLE PRECISION,PARAMETER :: Z = 2.0D0 ! Number of protons in He DOUBLE PRECISION,PARAMETER :: aa = 0.88534/SQRT(2.0*Z**0.23) DOUBLE PRECISION,PARAMETER :: aa1 = 0.1818 DOUBLE PRECISION,PARAMETER :: aa2 = 0.50990 DOUBLE PRECISION,PARAMETER :: aa3 = 0.2802 DOUBLE PRECISION,PARAMETER :: aa4 = 0.02817 DOUBLE PRECISION,PARAMETER :: bb1 = 3.2 DOUBLE PRECISION,PARAMETER :: bb2 = 0.9423 DOUBLE PRECISION,PARAMETER :: bb3 = 0.4029 DOUBLE PRECISION,PARAMETER :: bb4 = 0.2016 REAL(KIND=8),PARAMETER,DIMENSION(17) :: C = (/0.d0,0.d0,0.d0,0.d0,0.d0,0.4616213781D6,0.d0, & & 0.4460565781D7,0.d0,0.5803352873D8,0.d0,0.1031677697D10,0.d0,0.2415716766D11,0.d0,0.7191492488D12,0.d0 /) !C(6) = 0.4616213781D6 !C(8) = 0.4460565781D7 !C(10) = 0.5803352873D8 !C(12) = 0.1031677697D10 !C(14) = 0.2415716766D11 !C(16) = 0.7191492488D12 contains !------------------------------------------------------------------------------- ! ! Calculate pair potential phi(r) ! !------------------------------------------------------------------------------- subroutine calc_phi(R,phi) implicit none !-- Transferred variables double precision, intent(in) :: R double precision, intent(out) :: phi !-- Local variables REAL(KIND=8) :: fact REAL(KIND=8) :: ret0 REAL(KIND=8) :: ret1 REAL(KIND=8) :: ret2 REAL(KIND=8) :: atom INTEGER :: n INTEGER :: k INTEGER :: i ret1 = 0.d0 IF (R .gt. model_cutoff) THEN ! Argument exceeds cutoff radius phi = 0.d0 ELSE IF ( (R .gt. empir_cutoff) .AND. (R .lt. model_cutoff) ) THEN ! Use ab initio potential ret0 = A*EXP(a1*R+a2*R**(2.0) + a_1*R**(-1.0) & & + a_2*R**(-2.0) + d1*SIN(d2*R+d3)) DO n=3,8 ret2 = 0.0D0 DO k=0,(2*n) fact = 0.0D0 IF ( (k == 0) .OR. (k == 1) ) THEN fact = 1.0D0 ELSE fact = 1.0D0 DO i=2,k fact = fact*i END DO ! i END IF ! if (k) ret2 = ret2 + (b*R)**k/fact END DO ! k ret1 = ret1 + (C(2*n)/(R**(2*n)))*(1.0D0-EXP(-b*R)*ret2) END DO ! n phi = (ret0 - ret1)/TOK ! [au] ELSE ! Use screened Coulomb potential atom = aa1*EXP(-bb1*(R/aa)) + aa2*EXP(-bb2*(R/aa)) & & + aa3*EXP(-bb3*(R/aa)) + aa4*EXP(-bb4*(R/aa)) phi = (Z*Z/R)*atom*10**1.37503/TOEV END IF ! Convert from au to ev phi = phi*TOEV end subroutine calc_phi !------------------------------------------------------------------------------- ! ! Calculate pair potential phi(r) and its derivative dphi(r) ! !------------------------------------------------------------------------------- subroutine calc_phi_dphi(R,phi,dphi) implicit none !-- Transferred variables double precision, intent(in) :: R double precision, intent(out) :: phi,dphi !-- Local variables REAL(KIND=8) :: fact REAL(KIND=8) :: ret0 REAL(KIND=8) :: ret1 REAL(KIND=8) :: ret2 REAL(KIND=8) :: dret0 REAL(KIND=8) :: dret1 REAL(KIND=8) :: dret2 REAL(KIND=8) :: atom REAL(KIND=8) :: datom REAL(KIND=8) :: ksum1 REAL(KIND=8) :: ksum2 INTEGER :: n INTEGER :: k INTEGER :: i ret1 = 0.d0 if (R .gt. model_cutoff) then ! Argument exceeds cutoff radius phi = 0.d0 dphi = 0.d0 else if ( (R .gt. empir_cutoff) .and. (R .lt. model_cutoff) ) then ! Use ab initio potential ret0 = A*EXP(a1*R+a2*R**(2.0) + a_1*R**(-1.0) & & + a_2*R**(-2.0) + d1*SIN(d2*R+d3)) DO n=3,8 ret2 = 0.0D0 DO k=0,(2*n) fact = 0.0D0 IF ( (k == 0) .OR. (k == 1) ) THEN fact = 1.0D0 ELSE fact = 1.0D0 DO i=2,k fact = fact*i END DO ! i END IF ! if (k) ret2 = ret2 + (b*R)**k/fact END DO ! k ret1 = ret1 + (C(2*n)/(R**(2*n)))*(1.0D0-EXP(-b*R)*ret2) END DO ! n phi = (ret0 - ret1)/TOK ! [au] ! Calculate derivative of potential dret0 = ret0*(a1+2.d0*a2*R-a_1/(R*R)-2.d0*a_2/(R*R*R)+ & & d1*d2*COS(d2*R+d3)) dret1 = 0.d0 DO n=3,8 ksum1 = 0.d0 ksum2 = 0.d0 DO k=0,(2*n) fact = 0.d0 IF ( (k.eq.0) .or. (k.eq.1) ) THEN fact = 1.d0 ELSE fact = 1.d0 DO i=2,k fact = fact*i END DO ! i END IF ksum1 = ksum1 + (b*R)**k/fact ksum2 = ksum2 + (b**k)*DBLE(k)*(R**(k-1))/fact END DO ! k dret1 = dret1 - 2.d0*DBLE(n)*C(2*n)*(R**(-2*n-1))* & & (1.d0-EXP(-b*R)*ksum1) + C(2*n)/(R**(2*n))* & & (b*EXP(-b*R)*ksum1 - EXP(-b*R)*ksum2) END DO ! n dphi = (dret0 - dret1)/TOK else ! Use screened Coulomb potential atom = aa1*EXP(-bb1*(R/aa)) + aa2*EXP(-bb2*(R/aa)) & & + aa3*EXP(-bb3*(R/aa)) + aa4*EXP(-bb4*(R/aa)) datom = -(aa1*bb1/aa)*EXP(-bb1*(R/aa)) & & -(aa2*bb2/aa)*EXP(-bb2*(R/aa)) & & -(aa3*bb3/aa)*EXP(-bb3*(R/aa)) & & -(aa4*bb4/aa)*EXP(-bb4*(R/aa)) phi = (Z*Z/R)*atom*10**1.37503/TOEV dphi = -(Z*Z/(R*R))*atom + (Z*Z/R)*datom dphi = dphi*10**1.37503/TOEV endif ! Convert to eV phi = phi*TOEV dphi = dphi*TOEV 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 double precision :: Rij(DIM) double precision :: r,Rsqij,phi,dphi,dEidr integer :: ier integer :: i,j,jj,numnei,atom_ret,comp_force,comp_enepot,comp_virial,comp_energy integer, allocatable, target :: nei1atom_substitute(:) character*80 :: error_message !-- KIM variables integer N; pointer(pN,N) 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 boxSideLengths(DIM); pointer(pboxSideLengths,boxSideLengths) double precision Rij_list(DIM,1); pointer(pRij_list,Rij_list) integer numContrib; pointer(pnumContrib,numContrib) integer nei1atom(1); pointer(pnei1atom,nei1atom) integer particleTypes(1); pointer(pparticleTypes,particleTypes) double precision virialdum(1); pointer(pvirial,virialdum) character*64 NBC_Method; pointer(pNBC_Method,NBC_Method) double precision, pointer :: coor(:,:),force(:,:),ene_pot(:),virial_global(:) integer IterOrLoca integer HalfOrFull integer NBC integer numberContrib integer idum ! Determine neighbor list boundary condition (NBC) ! and half versus full mode: ! ***************************** ! * HalfOrFull = 1 -- Half ! * = 2 -- Full ! ***************************** ! ! 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 if (index(NBC_Method,"CLUSTER").eq.1) then NBC = 0 HalfOrFull = 1 elseif (index(NBC_Method,"MI_OPBC_H").eq.1) then NBC = 1 HalfOrFull = 1 elseif (index(NBC_Method,"MI_OPBC_F").eq.1) then NBC = 1 HalfOrFull = 2 elseif (index(NBC_Method,"NEIGH_PURE_H").eq.1) then NBC = 2 HalfOrFull = 1 elseif (index(NBC_Method,"NEIGH_PURE_F").eq.1) then NBC = 2 HalfOrFull = 2 elseif (index(NBC_Method,"NEIGH_RVEC_F").eq.1) then NBC = 3 HalfOrFull = 2 else ier = KIM_STATUS_FAIL idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "Unknown NBC method", ier) goto 42 endif call free(pNBC_Method) ! don't forget to release the memory... ! Determine neighbor list handling mode ! if (NBC.ne.0) then !***************************** !* IterOrLoca = 1 -- Iterator !* = 2 -- Locator !***************************** IterOrLoca = kim_api_get_neigh_mode_f(pkim, ier) if (ier.lt.KIM_STATUS_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_neigh_mode_f", ier) goto 42 endif if (IterOrLoca.ne.1 .and. IterOrLoca.ne.2) then ier = KIM_STATUS_FAIL write(error_message,'(a,i1)') & 'Unsupported IterOrLoca mode = ',IterOrLoca idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & error_message, ier) goto 42 endif else IterOrLoca = 2 ! for CLUSTER NBC endif ! Check to see if we have been asked to compute the forces, energyperatom, ! energy 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 ! Unpack data from KIM object ! call kim_api_getm_data_f(pkim, ier, & "numberOfParticles", pN, 1, & "particleTypes", pparticleTypes, 1, & "coordinates", pcoor, 1, & "numberContributingParticles", pnumContrib, TRUEFALSE(HalfOrFull.eq.1), & "boxSideLengths", pboxSideLengths, TRUEFALSE(NBC.eq.1), & "energy", penergy, TRUEFALSE(comp_energy.eq.1), & "forces", pforce, TRUEFALSE(comp_force.eq.1), & "particleEnergy", penepot, TRUEFALSE(comp_enepot.eq.1), & "virial", pvirial, TRUEFALSE(comp_virial.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,coor,DIM,N) if (comp_force.eq.1) call KIM_to_F90_real_array_2d(forcedum,force,DIM,N) if (comp_enepot.eq.1) call KIM_to_F90_real_array_1d(enepotdum,ene_pot,N) if (comp_virial.eq.1) call KIM_to_F90_real_array_1d(virialdum,virial_global,6) if (HalfOrFull.eq.1) then if (NBC.ne.0) then ! non-CLUSTER cases numberContrib = numContrib else ! CLUSTER case numberContrib = N endif endif ! Check to be sure that the atom types are correct ! ier = KIM_STATUS_FAIL ! assume an error do i = 1,N if (particleTypes(i).ne.speccode) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "Unexpected species type detected", 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:N) = 0.d0 if (comp_energy.eq.1) energy = 0.d0 if (comp_force.eq.1) force(1:3,1:N) = 0.d0 if (comp_virial.eq.1) virial_global = 0.d0 ! Initialize neighbor handling for CLUSTER NBC ! if (NBC.eq.0) then allocate( nei1atom_substitute(N) ) pnei1atom = loc(nei1atom_substitute) endif ! Initialize neighbor handling for Iterator mode ! if (IterOrLoca.eq.1) then ier = kim_api_get_neigh_f(pkim,0,0,atom_ret,numnei,pnei1atom,pRij_list) ! check for successful initialization if (ier.ne.KIM_STATUS_NEIGH_ITER_INIT_OK) then idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_neigh_f", ier) ier = KIM_STATUS_FAIL goto 42 endif endif ! ! Compute energy and forces ! ! Loop over particles and compute energy and forces ! i = 0 do ! Set up neighbor list for next atom for all NBC methods ! if (IterOrLoca.eq.1) then ! ITERATOR mode ier = kim_api_get_neigh_f(pkim,0,1,atom_ret,numnei,pnei1atom,pRij_list) if (ier.eq.KIM_STATUS_NEIGH_ITER_PAST_END) exit ! incremented past the end of the list, ! terminate loop if (ier.lt.KIM_STATUS_OK) then ! some sort of problem, exit idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_neigh_f", ier) goto 42 endif i = atom_ret else ! LOCATOR mode i = i + 1 if (i.gt.N) exit ! incremented past end of list, ! terminate loop if (NBC.eq.0) then ! CLUSTER NBC method numnei = N - i ! number of neighbors in list i+1, ..., N nei1atom(1:numnei) = (/ (i+jj, jj = 1,numnei) /) ier = KIM_STATUS_OK else ier = kim_api_get_neigh_f(pkim,1,i,atom_ret,numnei,pnei1atom,pRij_list) if (ier.ne.KIM_STATUS_OK) then ! some sort of problem, exit idum = kim_api_report_error_f(__LINE__, THIS_FILE_NAME, & "kim_api_get_neigh_f", ier) ier = KIM_STATUS_FAIL goto 42 endif endif endif ! Loop over the neighbors of atom i ! do jj = 1, numnei j = nei1atom(jj) ! get neighbor ID ! compute relative position vector ! if (NBC.ne.3) then ! all methods except NEIGH_RVEC_F Rij(:) = coor(:,j) - coor(:,i) ! distance vector between i j else Rij(:) = Rij_list(:,jj) endif ! apply periodic boundary conditions if required ! if (NBC.eq.1) then where ( abs(Rij) .gt. 0.5d0*boxSideLengths ) ! periodic boundary conditions Rij = Rij - sign(boxSideLengths,Rij) ! applied where needed. end where ! endif ! compute energy and forces ! Rsqij = dot_product(Rij,Rij) ! compute square distance if ( Rsqij .lt. model_cutsq ) then ! particles are interacting? r = sqrt(Rsqij) ! compute distance if (comp_force.eq.1.or.comp_virial.eq.1) then call calc_phi_dphi(r,phi,dphi) ! compute pair potential ! and it derivative if ((HalfOrFull.eq.1) .and. & (j .le. numberContrib)) then ! HALF mode dEidr = dphi ! double contribution else ! FULL mode dEidr = 0.5d0*dphi ! regular contribution endif else call calc_phi(r,phi) ! compute just pair potential endif ! contribution to energy ! if (comp_enepot.eq.1) then ene_pot(i) = ene_pot(i) + 0.5d0*phi ! accumulate energy if ((HalfOrFull.eq.1) .and. & (j .le. numberContrib)) & ! HALF mode ene_pot(j) = ene_pot(j) + 0.5d0*phi! (i and j share it) endif if (comp_energy.eq.1) then if ((HalfOrFull.eq.1) .and. & (j .le. numberContrib)) then ! HALF mode energy = energy + phi ! add v to total energy else ! FULL mode energy = energy + 0.5d0*phi ! add half v to total energy endif endif ! contribution to virial tensor, virial(i,j)=r(i)*r(j)*(dV/dr)/r ! if (comp_virial.eq.1) then virial_global(1) = virial_global(1) + Rij(1)*Rij(1)*dEidr/r virial_global(2) = virial_global(2) + Rij(2)*Rij(2)*dEidr/r virial_global(3) = virial_global(3) + Rij(3)*Rij(3)*dEidr/r virial_global(4) = virial_global(4) + Rij(2)*Rij(3)*dEidr/r virial_global(5) = virial_global(5) + Rij(1)*Rij(3)*dEidr/r virial_global(6) = virial_global(6) + Rij(1)*Rij(2)*dEidr/r endif ! contribution to forces ! if (comp_force.eq.1) 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 enddo ! loop on jj enddo ! infinite do loop (terminated by exit statements above) ! Free temporary storage ! if (NBC.eq.0) deallocate( nei1atom_substitute ) ! Everything is great ! ier = KIM_STATUS_OK 42 continue Compute_Energy_Forces = ier return end function Compute_Energy_Forces end module model_He_P_AbIn !------------------------------------------------------------------------------- ! ! Model initialization routine (REQUIRED) ! !------------------------------------------------------------------------------- integer function model_init(pkim) use model_He_P_AbIn 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 cutoff; pointer(pcutoff,cutoff) ! 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 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) goto 42 endif cutoff = model_cutoff ier = KIM_STATUS_OK 42 continue model_init = ier return end function model_init