diff --git a/lib/meam/Makefile b/lib/meam/Makefile new file mode 100755 index 0000000000..7997be14dd --- /dev/null +++ b/lib/meam/Makefile @@ -0,0 +1,51 @@ +# * +# *_________________________________________________________________________* +# * MEAM: MODEFIED EMBEDDED ATOM METHOD * +# * DESCRIPTION: SEE READ-ME * +# * FILE NAME: Makefile * +# * AUTHORS: Greg Wagner, Sandia National Laboratories * +# * CONTACT: gjwagne@sandia.gov * +# *_________________________________________________________________________*/ + +SHELL = /bin/sh + +# ------ FILES ------ + +SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F + +FILES = $(SRC) Makefile + +# ------ DEFINITIONS ------ + +LIB = libmeam.a +OBJ = $(SRC:.F=.o) + +# ------ SETTINGS ------ + +F90 = gfortran +F90FLAGS = -g -fno-second-underscore +#F90FLAGS = -O +ARCHIVE = ar +ARCHFLAG = -rc +LINK = g++ +LINKFLAGS = -O +USRLIB = +SYSLIB = + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) + +# ------ COMPILE RULES ------ + +%.o:%.F + $(F90) $(F90FLAGS) -c $< + +# ------ CLEAN ------ + +clean: + -rm *.o *.mod *~ $(LIB) + +tar: + -tar -cvf ../MEAM.tar $(FILES) \ No newline at end of file diff --git a/lib/meam/Makefile.g95 b/lib/meam/Makefile.g95 new file mode 100644 index 0000000000..0125e4183d --- /dev/null +++ b/lib/meam/Makefile.g95 @@ -0,0 +1,48 @@ +# * +# *_________________________________________________________________________* +# * MEAM: MODEFIED EMBEDDED ATOM METHOD * +# * DESCRIPTION: SEE READ-ME * +# * FILE NAME: Makefile * +# * AUTHORS: Greg Wagner, Sandia National Laboratories * +# * CONTACT: gjwagne@sandia.gov * +# *_________________________________________________________________________*/ + +SHELL = /bin/sh + +# ------ FILES ------ + +SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F + +FILES = $(SRC) Makefile + +# ------ DEFINITIONS ------ + +LIB = libmeam.a +OBJ = $(SRC:.F=.o) + +# ------ SETTINGS ------ + +F90 = g95 +F90FLAGS = -O +ARCHIVE = ar +ARCHFLAG = -rc +USRLIB = +SYSLIB = + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) + +# ------ COMPILE RULES ------ + +%.o:%.F + $(F90) $(F90FLAGS) -c $< + +# ------ CLEAN ------ + +clean: + -rm *.o *.mod *~ $(LIB) + +tar: + -tar -cvf ../MEAM.tar $(FILES) \ No newline at end of file diff --git a/lib/meam/Makefile.gfortran b/lib/meam/Makefile.gfortran new file mode 100644 index 0000000000..7997be14dd --- /dev/null +++ b/lib/meam/Makefile.gfortran @@ -0,0 +1,51 @@ +# * +# *_________________________________________________________________________* +# * MEAM: MODEFIED EMBEDDED ATOM METHOD * +# * DESCRIPTION: SEE READ-ME * +# * FILE NAME: Makefile * +# * AUTHORS: Greg Wagner, Sandia National Laboratories * +# * CONTACT: gjwagne@sandia.gov * +# *_________________________________________________________________________*/ + +SHELL = /bin/sh + +# ------ FILES ------ + +SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F + +FILES = $(SRC) Makefile + +# ------ DEFINITIONS ------ + +LIB = libmeam.a +OBJ = $(SRC:.F=.o) + +# ------ SETTINGS ------ + +F90 = gfortran +F90FLAGS = -g -fno-second-underscore +#F90FLAGS = -O +ARCHIVE = ar +ARCHFLAG = -rc +LINK = g++ +LINKFLAGS = -O +USRLIB = +SYSLIB = + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) + +# ------ COMPILE RULES ------ + +%.o:%.F + $(F90) $(F90FLAGS) -c $< + +# ------ CLEAN ------ + +clean: + -rm *.o *.mod *~ $(LIB) + +tar: + -tar -cvf ../MEAM.tar $(FILES) \ No newline at end of file diff --git a/lib/meam/Makefile.ifort b/lib/meam/Makefile.ifort new file mode 100644 index 0000000000..66807b817d --- /dev/null +++ b/lib/meam/Makefile.ifort @@ -0,0 +1,48 @@ +# * +# *_________________________________________________________________________* +# * MEAM: MODEFIED EMBEDDED ATOM METHOD * +# * DESCRIPTION: SEE READ-ME * +# * FILE NAME: Makefile * +# * AUTHORS: Greg Wagner, Sandia National Laboratories * +# * CONTACT: gjwagne@sandia.gov * +# *_________________________________________________________________________*/ + +SHELL = /bin/sh + +# ------ FILES ------ + +SRC = meam_data.F meam_setup_done.F meam_setup_global.F meam_setup_param.F meam_dens_init.F meam_dens_final.F meam_force.F + +FILES = $(SRC) Makefile + +# ------ DEFINITIONS ------ + +LIB = libmeam.a +OBJ = $(SRC:.F=.o) + +# ------ SETTINGS ------ + +F90 = ifort +F90FLAGS = -O +ARCHIVE = ar +ARCHFLAG = -rc +USRLIB = +SYSLIB = + +# ------ MAKE PROCEDURE ------ + +lib: $(OBJ) + $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) + +# ------ COMPILE RULES ------ + +%.o:%.F + $(F90) $(F90FLAGS) -c $< + +# ------ CLEAN ------ + +clean: + -rm *.o *.mod *~ $(LIB) + +tar: + -tar -cvf ../MEAM.tar $(FILES) \ No newline at end of file diff --git a/lib/meam/README b/lib/meam/README new file mode 100644 index 0000000000..f6d051419c --- /dev/null +++ b/lib/meam/README @@ -0,0 +1,20 @@ +MEAM (modified embedded atom method) library + +Greg Wagner, Sandia National Labs +gjwagne at sandia.gov +Jan 2007 + +-------------- + +This library is in implementation of the MEAM potential, specifically +designed to work with LAMMPS. + +This library must be built with a F90 compiler, before LAMMPS is +built, so LAMMPS can link against it. + +Build the library using one of the provided Makefiles or creating your +own, specific to your compiler and system. For example: + +make -f Makefile.g95 + +If the build is successful, you should end up with a libmeam.a file. diff --git a/lib/meam/meam_data.F b/lib/meam/meam_data.F new file mode 100755 index 0000000000..9269c7ccb7 --- /dev/null +++ b/lib/meam/meam_data.F @@ -0,0 +1,76 @@ + + module meam_data + + integer, parameter :: maxelt = 5 + +c cutforce = force cutoff +c cutforcesq = force cutoff squared + + real*8 cutforce,cutforcesq + +c Ec_meam = cohesive energy +c re_meam = nearest-neighbor distance +c Omega_meam = atomic volume +c B_meam = bulk modulus +c Z_meam = number of first neighbors for reference structure +c ielt_meam = atomic number of element +c A_meam = adjustable parameter +c alpha_meam = sqrt(9*Omega*B/Ec) +c rho0_meam = density scaling parameter +c delta_meam = heat of formation for alloys +c beta[0-3]_meam = electron density constants +c t[0-3]_meam = coefficients on densities in Gamma computation +c ibar_meam(i) = selection parameter for Gamma function for elt i, +c lattce_meam(i,j) = lattce configuration for elt i or alloy (i,j) +c neltypes = maximum number of element type defined +c eltind = index number of pair (similar to Voigt notation; ij = ji) +c phir = pair potential function array +c phirar[1-6] = spline coeffs +c attrac_meam = attraction parameter in Rose energy +c repuls_meam = repulsion parameter in Rose energy +c nn2_meam = 1 if second nearest neighbors are to be computed, else 0 +c Cmin_meam, Cmax_meam = min and max values in screening cutoff +c rc_meam = cutoff distance for meam +c delr_meam = cutoff region for meam +c ebound_meam = factor giving maximum boundary of sceen fcn ellipse +c augt1 = flag for whether t1 coefficient should be augmented +c gsmooth_factor = factor determining length of G smoothing region +c vind[23]D = Voight notation index maps for 2 and 3D +c v2D,v3D = array of factors to apply for Voight notation + +c nr,dr = pair function discretization parameters +c nrar,rdrar = spline coeff array parameters + + real*8 Ec_meam(maxelt,maxelt),re_meam(maxelt,maxelt) + real*8 Omega_meam(maxelt),Z_meam(maxelt) + real*8 A_meam(maxelt),alpha_meam(maxelt,maxelt),rho0_meam(maxelt) + real*8 delta_meam(maxelt,maxelt) + real*8 beta0_meam(maxelt),beta1_meam(maxelt) + real*8 beta2_meam(maxelt),beta3_meam(maxelt) + real*8 t0_meam(maxelt),t1_meam(maxelt) + real*8 t2_meam(maxelt),t3_meam(maxelt) + integer ibar_meam(maxelt),ielt_meam(maxelt) + character*3 lattce_meam(maxelt,maxelt) + integer nn2_meam(maxelt,maxelt) + integer eltind(maxelt,maxelt) + integer neltypes + + real*8, allocatable :: phir(:,:) + + real*8, allocatable :: phirar(:,:),phirar1(:,:),phirar2(:,:), + $ phirar3(:,:),phirar4(:,:),phirar5(:,:),phirar6(:,:) + + real*8 attrac_meam(maxelt,maxelt),repuls_meam(maxelt,maxelt) + + real*8 Cmin_meam(maxelt,maxelt,maxelt) + real*8 Cmax_meam(maxelt,maxelt,maxelt) + real*8 rc_meam,delr_meam,ebound_meam(maxelt,maxelt) + integer augt1 + real*8 gsmooth_factor + integer vind2D(3,3),vind3D(3,3,3) + integer v2D(6),v3D(10) + + integer nr,nrar + real*8 dr,rdrar + + end module diff --git a/lib/meam/meam_dens_final.F b/lib/meam/meam_dens_final.F new file mode 100755 index 0000000000..cf0088bcf4 --- /dev/null +++ b/lib/meam/meam_dens_final.F @@ -0,0 +1,201 @@ +c Extern "C" declaration has the form: +c +c void meam_dens_final_(int *, int *, int *, double *, int *, int *, int *, +c double *, double *, double *, double *, double *, double *, +c double *, double *, double *, double *, double *, double *, +c double *, double *, double *, double *, int *); +c +c Call from pair_meam.cpp has the form: +c +c meam_dens_final_(&nlocal,&nmax,&eflag,&eng_vdwl,ntype,type,fmap, +c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], +c &arho3b[0][0],&t_ave[0][0],gamma,dgamma1, +c dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); +c + + subroutine meam_dens_final(nlocal, nmax, eflag, eng_vdwl, + $ ntype, type, fmap, + $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, + $ Gamma, dGamma1, dGamma2, dGamma3, + $ rho, rho0, rho1, rho2, rho3, fp, errorflag) + + use meam_data + implicit none + + integer nlocal, nmax, eflag, ntype, type, fmap + real*8 eng_vdwl, Arho1, Arho2 + real*8 Arho2b, Arho3, Arho3b + real*8 t_ave + real*8 Gamma, dGamma1, dGamma2, dGamma3 + real*8 rho, rho0, rho1, rho2, rho3 + real*8 fp + integer errorflag + + dimension type(nmax), fmap(ntype) + dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax) + dimension Arho3(10,nmax), Arho3b(3,nmax), t_ave(3,nmax) + dimension Gamma(nmax), dGamma1(nmax), dGamma2(nmax) + dimension dGamma3(nmax), rho(nmax), rho0(nmax) + dimension rho1(nmax), rho2(nmax), rho3(nmax) + dimension fp(nmax) + + integer i, elti + integer m + real*8 rhob, G, dG, Gbar, dGbar, gam, shp(3), shpi(3), Z + real*8 B, denom + +c Complete the calculation of density + + do i = 1,nlocal + + elti = fmap(type(i)) + rho1(i) = 0.d0 + rho2(i) = -1.d0/3.d0*Arho2b(i)*Arho2b(i) + rho3(i) = 0.d0 + do m = 1,3 + rho1(i) = rho1(i) + Arho1(m,i)*Arho1(m,i) + rho3(i) = rho3(i) - 3.d0/5.d0*Arho3b(m,i)*Arho3b(m,i) + enddo + do m = 1,6 + rho2(i) = rho2(i) + v2D(m)*Arho2(m,i)*Arho2(m,i) + enddo + do m = 1,10 + rho3(i) = rho3(i) + v3D(m)*Arho3(m,i)*Arho3(m,i) + enddo + + if( rho0(i) .gt. 0.0 ) then + t_ave(1,i) = t_ave(1,i)/rho0(i) + t_ave(2,i) = t_ave(2,i)/rho0(i) + t_ave(3,i) = t_ave(3,i)/rho0(i) + endif + + Gamma(i) = t_ave(1,i)*rho1(i) + $ + t_ave(2,i)*rho2(i) + t_ave(3,i)*rho3(i) + + if( rho0(i) .gt. 0.0 ) then + Gamma(i) = Gamma(i)/(rho0(i)*rho0(i)) + end if + + call G_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,errorflag) + if (errorflag.ne.0) return + if (ibar_meam(elti).le.0) then + Gbar = 1.d0 + else + call get_shpfcn(shp,lattce_meam(elti,elti)) + gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2) + $ +t_ave(3,i)*shp(3))/(Z_meam(elti)**2) + call G_gam(gam,ibar_meam(elti),gsmooth_factor, + $ Gbar,errorflag) + endif + rho(i) = rho0(i) * G + rhob = rho(i)/(rho0_meam(elti)*Z_meam(elti)*Gbar) + + Z = Z_meam(elti) + call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG) + if (ibar_meam(elti).le.0) then + Gbar = 1.d0 + dGbar = 0.d0 + else + call get_shpfcn(shpi,lattce_meam(elti,elti)) + gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) + $ +t_ave(3,i)*shpi(3))/(Z*Z) + call dG_gam(gam,ibar_meam(elti),gsmooth_factor, + $ Gbar,dGbar) + endif + denom = 1.d0/(rho0_meam(elti)*Z*Gbar) + dGamma1(i) = (G - 2*dG*Gamma(i))*denom + + if( rho0(i) .ne. 0.d0 ) then + dGamma2(i) = (dG/rho0(i))*denom + else + dGamma2(i) = 0.d0 + end if + + dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom + + B = A_meam(elti)*Ec_meam(elti,elti) + + if( rhob .ne. 0.d0 ) then + fp(i) = B*(log(rhob)+1.d0) + if (eflag.eq.1) eng_vdwl = eng_vdwl + B*rhob*log(rhob) + else + fp(i) = B + endif + + enddo + + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine G_gam(Gamma,ibar,gsmooth_factor,G,errorflag) +c Compute G(Gamma) based on selection flag ibar: +c 0 => G = sqrt(1+Gamma) +c 1 => G = exp(Gamma/2) +c 2 => not implemented +c 3 => G = 2/(1+exp(-Gamma)) +c 4 => G = sqrt(1+Gamma) + implicit none + real*8 Gamma,G + real*8 gsmooth_factor, gsmooth_switchpoint + integer ibar, errorflag + if (ibar.eq.0.or.ibar.eq.4) then + gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1) + if (Gamma.lt.gsmooth_switchpoint) then +c e.g. gsmooth_factor is 99, then: +c gsmooth_switchpoint = -0.99 +c G = 0.01*(-0.99/Gamma)**99 + G = 1/(gsmooth_factor+1) + $ *(gsmooth_switchpoint/Gamma)**gsmooth_factor + G = sqrt(G) + else + G = sqrt(1.d0+Gamma) + endif + else if (ibar.eq.1) then + G = exp(Gamma/2.d0) + else if (ibar.eq.3) then + G = 2.d0/(1.d0+exp(-Gamma)) + else + errorflag = 1 + endif + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine dG_gam(Gamma,ibar,gsmooth_factor,G,dG) +c Compute G(Gamma) and dG(gamma) based on selection flag ibar: +c 0 => G = sqrt(1+Gamma) +c 1 => G = exp(Gamma/2) +c 2 => not implemented +c 3 => G = 2/(1+exp(-Gamma)) +c 4 => G = sqrt(1+Gamma) + real*8 Gamma,G,dG + real*8 gsmooth_factor, gsmooth_switchpoint + integer ibar + if (ibar.eq.0.or.ibar.eq.4) then + gsmooth_switchpoint = -gsmooth_factor / (gsmooth_factor+1) + if (Gamma.lt.gsmooth_switchpoint) then +c e.g. gsmooth_factor is 99, then: +c gsmooth_switchpoint = -0.99 +c G = 0.01*(-0.99/Gamma)**99 + G = 1/(gsmooth_factor+1) + $ *(gsmooth_switchpoint/Gamma)**gsmooth_factor + G = sqrt(G) + dG = -gsmooth_factor*G/(2.0*Gamma) + else + G = sqrt(1.d0+Gamma) + dG = 1.d0/(2.d0*G) + endif + else if (ibar.eq.1) then + G = exp(Gamma/2.d0) + dG = G/2.d0 + else if (ibar.eq.3) then + G = 2.d0/(1.d0+exp(-Gamma)) + dG = G*(2.d0-G)/2 + endif + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/lib/meam/meam_dens_init.F b/lib/meam/meam_dens_init.F new file mode 100755 index 0000000000..bc992ea95d --- /dev/null +++ b/lib/meam/meam_dens_init.F @@ -0,0 +1,521 @@ +c Extern "C" declaration has the form: +c +c void meam_dens_init_(int *, int *, int *, double *, int *, int *, int *, double *, +c int *, int *, int *, int *, +c double *, double *, double *, double *, double *, double *, +c double *, double *, double *, double *, int *); +c +c +c Call from pair_meam.cpp has the form: +c +c meam_dens_init_(&i,&nmax,&eflag,&eng_vdwl,ntype,type,fmap,&x[0][0], +c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], +c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], +c rho0,&arho1[0][0],&arho2[0][0],arho2b, +c &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&errorflag); +c + + subroutine meam_dens_init(i, nmax, eflag, eng_vdwl, + $ ntype, type, fmap, x, + $ numneigh, firstneigh, + $ numneigh_full, firstneigh_full, + $ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b, + $ arho3, arho3b, t_ave, errorflag) + + use meam_data + implicit none + + integer i, nmax, eflag, ntype, type, fmap + real*8 eng_vdwl + real*8 x + integer numneigh, firstneigh, numneigh_full, firstneigh_full + real*8 scrfcn, dscrfcn, fcpair + real*8 rho0, arho1, arho2 + real*8 arho2b, arho3, arho3b, t_ave + integer errorflag + integer j,jn + + dimension x(3,nmax) + dimension type(nmax), fmap(ntype) + dimension firstneigh(numneigh), firstneigh_full(numneigh_full) + dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh) + dimension rho0(nmax), arho1(3,nmax), arho2(6,nmax) + dimension arho2b(nmax), arho3(10,nmax), arho3b(3,nmax) + dimension t_ave(3,nmax) + + errorflag = 0 + +c Compute screening function and derivatives + call getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x, + $ numneigh, firstneigh, + $ numneigh_full, firstneigh_full, + $ ntype, type, fmap) + +c Calculate intermediate density terms to be communicated + call calc_rho1(i, nmax, ntype, type, fmap, x, + $ numneigh, firstneigh, + $ scrfcn, fcpair, rho0, arho1, arho2, arho2b, + $ arho3, arho3b, t_ave) + + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine getscreen(i, nmax, scrfcn, dscrfcn, fcpair, x, + $ numneigh, firstneigh, + $ numneigh_full, firstneigh_full, + $ ntype, type, fmap) + + use meam_data + implicit none + + integer i, nmax + real*8 scrfcn, dscrfcn, fcpair, x + integer numneigh, firstneigh, numneigh_full, firstneigh_full + integer ntype, type, fmap + + dimension scrfcn(numneigh), dscrfcn(numneigh) + dimension fcpair(numneigh), x(3,nmax) + dimension firstneigh(numneigh), firstneigh_full(numneigh_full) + dimension type(nmax), fmap(ntype) + + integer jn,j,kn,k + integer elti,eltj,eltk + real*8 xitmp,yitmp,zitmp,delxij,delyij,delzij,rij2,rij + real*8 xjtmp,yjtmp,zjtmp,delxik,delyik,delzik,rik2,rik + real*8 xktmp,yktmp,zktmp,delxjk,delyjk,delzjk,rjk2,rjk + real*8 xik,xjk,sij,fcij,sfcij,dfcij,sikj,dfikj,cikj + real*8 Cmin,Cmax,delc,ebound,rbound,a,coef1,coef2 + real*8 coef1a,coef1b,coef2a,coef2b + real*8 dcikj + real*8 dC1a,dC1b,dC2a,dC2b + real*8 rnorm,fc,dfc,drinv + + drinv = 1.d0/delr_meam + elti = fmap(type(i)) + + xitmp = x(1,i) + yitmp = x(2,i) + zitmp = x(3,i) + + do jn = 1,numneigh + j = firstneigh(jn) + +c First compute screening function itself, sij + xjtmp = x(1,j) + yjtmp = x(2,j) + zjtmp = x(3,j) + delxij = xjtmp - xitmp + delyij = yjtmp - yitmp + delzij = zjtmp - zitmp + rij2 = delxij*delxij + delyij*delyij + delzij*delzij + rij = sqrt(rij2) + if (rij.gt.rc_meam) then + fcij = 0.0 + dfcij = 0.d0 + sij = 0.d0 + else + rnorm = (rc_meam-rij)*drinv + call screen(i, j, nmax, x, rij2, sij, + $ numneigh_full, firstneigh_full, ntype, type, fmap) + call dfcut(rnorm,fc,dfc) + fcij = fc + dfcij = dfc*drinv + endif + +c Now compute derivatives + dscrfcn(jn) = 0.d0 + sfcij = sij*fcij + if (sfcij.eq.0.d0.or.sfcij.eq.1.d0) goto 100 + eltj = fmap(type(j)) + rbound = ebound_meam(elti,eltj) * rij2 + do kn = 1,numneigh_full + k = firstneigh_full(kn) + if (k.eq.j) goto 10 + eltk = fmap(type(k)) + xktmp = x(1,k) + yktmp = x(2,k) + zktmp = x(3,k) + delxjk = xktmp - xjtmp + delyjk = yktmp - yjtmp + delzjk = zktmp - zjtmp + rjk2 = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk + if (rjk2.gt.rbound) goto 10 + delxik = xktmp - xitmp + delyik = yktmp - yitmp + delzik = zktmp - zitmp + rik2 = delxik*delxik + delyik*delyik + delzik*delzik + if (rik2.gt.rbound) goto 10 + xik = rik2/rij2 + xjk = rjk2/rij2 + a = 1 - (xik-xjk)*(xik-xjk) + if (a.eq.0.d0) goto 10 + cikj = (2.d0*(xik+xjk) + a - 2.d0)/a + Cmax = Cmax_meam(elti,eltj,eltk) + Cmin = Cmin_meam(elti,eltj,eltk) + if (cikj.ge.Cmax.or.cikj.lt.0.d0) then + goto 10 +c Note that we never have 0 ebound*rijsq, atom k is definitely outside the ellipse + rbound = ebound_meam(elti,eltj)*rijsq + + do nk = 1,numneigh_full + k = firstneigh_full(nk) + eltk = fmap(type(k)) + if (k.eq.j) goto 10 + delxjk = x(1,k) - x(1,j) + delyjk = x(2,k) - x(2,j) + delzjk = x(3,k) - x(3,j) + rjksq = delxjk*delxjk + delyjk*delyjk + delzjk*delzjk + if (rjksq.gt.rbound) goto 10 + delxik = x(1,k) - x(1,i) + delyik = x(2,k) - x(2,i) + delzik = x(3,k) - x(3,i) + riksq = delxik*delxik + delyik*delyik + delzik*delzik + if (riksq.gt.rbound) goto 10 + xik = riksq/rijsq + xjk = rjksq/rijsq + a = 1 - (xik-xjk)*(xik-xjk) + if (a.eq.0.d0) goto 10 + cikj = (2.d0*(xik+xjk) + a - 2.d0)/a + Cmax = Cmax_meam(elti,eltj,eltk) + Cmin = Cmin_meam(elti,eltj,eltk) + if (cikj.ge.Cmax.or.cikj.lt.0.d0) then + goto 10 + else if (cikj.le.Cmin) then + sij = 0.d0 + goto 20 + else + delc = Cmax - Cmin + cikj = (cikj-Cmin)/delc + call fcut(cikj,sikj) + endif + sij = sij * sikj + 10 continue + enddo + + 20 continue + + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2, + $ ntype,type,fmap,x,scrfcn,fcpair) +c Inputs: i,j,k = id's of 3 atom triplet +c jn = id of i-j pair +c rij2 = squared distance between i and j +c Outputs: dsij1 = deriv. of sij w.r.t. rik +c dsij2 = deriv. of sij w.r.t. rjk + use meam_data + implicit none + integer i,j,k,jn,nmax,numneigh + integer elti,eltj,eltk + real*8 rij2,rik2,rjk2,dsij1,dsij2 + integer ntype, type, fmap + real*8 x, scrfcn, fcpair + + dimension type(nmax), fmap(ntype) + dimension x(3,nmax), scrfcn(numneigh), fcpair(numneigh) + + real*8 dxik,dyik,dzik + real*8 dxjk,dyjk,dzjk + real*8 rbound,delc,sij,xik,xjk,cikj,sikj,dfc,a + real*8 Cmax,Cmin,dCikj1,dCikj2 + + sij = scrfcn(jn)*fcpair(jn) + elti = fmap(type(i)) + eltj = fmap(type(j)) + eltk = fmap(type(k)) + Cmax = Cmax_meam(elti,eltj,eltk) + Cmin = Cmin_meam(elti,eltj,eltk) + + dsij1 = 0.d0 + dsij2 = 0.d0 + if ((sij.ne.0.d0).and.(sij.ne.1.d0)) then + rbound = rij2*ebound_meam(elti,eltj) + delc = Cmax-Cmin + dxjk = x(1,k) - x(1,j) + dyjk = x(2,k) - x(2,j) + dzjk = x(3,k) - x(3,j) + rjk2 = dxjk*dxjk + dyjk*dyjk + dzjk*dzjk + if (rjk2.le.rbound) then + dxik = x(1,k) - x(1,i) + dyik = x(2,k) - x(2,i) + dzik = x(3,k) - x(3,i) + rik2 = dxik*dxik + dyik*dyik + dzik*dzik + if (rik2.le.rbound) then + xik = rik2/rij2 + xjk = rjk2/rij2 + a = 1 - (xik-xjk)*(xik-xjk) + if (a.ne.0.d0) then + cikj = (2.d0*(xik+xjk) + a - 2.d0)/a + if (cikj.ge.Cmin.and.cikj.le.Cmax) then + cikj = (cikj-Cmin)/delc + call dfcut(cikj,sikj,dfc) + call dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2) + a = sij/delc*dfc/sikj + dsij1 = a*dCikj1 + dsij2 = a*dCikj2 + endif + endif + endif + endif + endif + + return + end + + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine fcut(xi,fc) +c cutoff function + implicit none + real*8 xi,fc + real*8 a + if (xi.ge.1.d0) then + fc = 1.d0 + else if (xi.le.0.d0) then + fc = 0.d0 + else + a = 1.d0-xi + a = a*a + a = a*a + a = 1.d0-a + fc = a*a +c fc = xi + endif + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine dfcut(xi,fc,dfc) +c cutoff function and its derivative + implicit none + real*8 xi,fc,dfc,a,a3,a4 + if (xi.ge.1.d0) then + fc = 1.d0 + dfc = 0.d0 + else if (xi.le.0.d0) then + fc = 0.d0 + dfc = 0.d0 + else + a = 1.d0-xi + a3 = a*a*a + a4 = a*a3 + fc = (1.d0-a4)**2 + dfc = 8*(1.d0-a4)*a3 +c fc = xi +c dfc = 1.d0 + endif + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine dCfunc(rij2,rik2,rjk2,dCikj) +c Inputs: rij,rij2,rik2,rjk2 +c Outputs: dCikj = derivative of Cikj w.r.t. rij + implicit none + real*8 rij2,rik2,rjk2,dCikj + real*8 rij4,a,b,denom + + rij4 = rij2*rij2 + a = rik2-rjk2 + b = rik2+rjk2 + denom = rij4 - a*a + denom = denom*denom + dCikj = -4*(-2*rij2*a*a + rij4*b + a*a*b)/denom + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + subroutine dCfunc2(rij2,rik2,rjk2,dCikj1,dCikj2) +c Inputs: rij,rij2,rik2,rjk2 +c Outputs: dCikj1 = derivative of Cikj w.r.t. rik +c dCikj2 = derivative of Cikj w.r.t. rjk + implicit none + real*8 rij2,rik2,rjk2,dCikj1,dCikj2 + real*8 rij4,rik4,rjk4,a,b,denom + + rij4 = rij2*rij2 + rik4 = rik2*rik2 + rjk4 = rjk2*rjk2 + a = rik2-rjk2 + b = rik2+rjk2 + denom = rij4 - a*a + denom = denom*denom + dCikj1 = 4*rij2*(rij4 + rik4 + 2*rik2*rjk2 - 3*rjk4 - 2*rij2*a)/ + $ denom + dCikj2 = 4*rij2*(rij4 - 3*rik4 + 2*rik2*rjk2 + rjk4 + 2*rij2*a)/ + $ denom + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + + + + diff --git a/lib/meam/meam_force.F b/lib/meam/meam_force.F new file mode 100755 index 0000000000..d2a910c7e7 --- /dev/null +++ b/lib/meam/meam_force.F @@ -0,0 +1,438 @@ +c Extern "C" declaration has the form: +c +c void meam_force_(int *, int *, int *, double *, int *, int *, int *, double *, +c int *, int *, int *, int *, double *, double *, +c double *, double *, double *, double *, double *, double *, +c double *, double *, double *, double *, double *, double *, +c double *, double *, double *, double *, double *, int *); +c +c Call from pair_meam.cpp has the form: +c +c meam_force_(&i,&nmax,&eflag,&eng_vdwl,&ntype,type,fmap,&x[0][0], +c &numneigh[i],firstneigh[i],&numneigh_full[i],firstneigh_full[i], +c &scrfcn[offset],&dscrfcn[offset],&fcpair[offset], +c dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop, +c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0], +c &t_ave[0][0],&f[0][0],&strssa[0][0][0],&errorflag); +c + + subroutine meam_force(i, nmax, eflag, eng_vdwl, + $ ntype, type, fmap, x, + $ numneigh, firstneigh, numneigh_full, firstneigh_full, + $ scrfcn, dscrfcn, fcpair, + $ dGamma1, dGamma2, dGamma3, rho0, rho1, rho2, rho3, fp, + $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, f, strssa, + $ errorflag) + + use meam_data + implicit none + + integer eflag, nmax, ntype, type, fmap + real*8 eng_vdwl, x + integer numneigh, firstneigh, numneigh_full, firstneigh_full + real*8 scrfcn, dscrfcn, fcpair + real*8 dGamma1, dGamma2, dGamma3 + real*8 rho0, rho1, rho2, rho3, fp + real*8 Arho1, Arho2, Arho2b + real*8 Arho3, Arho3b + real*8 t_ave, f, strssa + integer errorflag + + dimension type(nmax), fmap(ntype) + dimension x(3,nmax) + dimension firstneigh(numneigh), firstneigh_full(numneigh_full) + dimension scrfcn(numneigh), dscrfcn(numneigh), fcpair(numneigh) + dimension dGamma1(nmax), dGamma2(nmax), dGamma3(nmax) + dimension rho0(nmax), rho1(nmax), rho2(nmax), rho3(nmax), fp(nmax) + dimension Arho1(3,nmax), Arho2(6,nmax), Arho2b(nmax) + dimension Arho3(10,nmax), Arho3b(3,nmax) + dimension t_ave(3,nmax), f(3,nmax), strssa(3,3,nmax) + + integer i,j,jn,k,kn,kk,m,n,p,q,strscomp + integer nv2,nv3,elti,eltj,ind + real*8 xitmp,yitmp,zitmp,delij(3),delref(3),rij2,rij,rij3 + real*8 delik(3),deljk(3) + real*8 Eu,astar,astarp + real*8 pp,phiforce,dUdrij,dUdsij,dUdrijm(3),force,forcem + real*8 B,r,recip,phi,phip,rhop,a + real*8 sij,fcij,dfcij,ds(3) + real*8 a0,a1,a1i,a1j,a2,a2i,a2j + real*8 a3i,a3j,a3i1,a3i2,a3j1,a3j2 + real*8 G,dG,Gbar,dGbar,gam,shpi(3),shpj(3),Z,denom + real*8 ai,aj,ro0i,ro0j,invrei,invrej + real*8 b0,rhoa0j,drhoa0j,rhoa0i,drhoa0i + real*8 b1,rhoa1j,drhoa1j,rhoa1i,drhoa1i + real*8 b2,rhoa2j,drhoa2j,rhoa2i,drhoa2i + real*8 a3,a3a,b3,rhoa3j,drhoa3j,rhoa3i,drhoa3i + real*8 drho0dr1,drho0dr2,drho0ds1,drho0ds2 + real*8 drho1dr1,drho1dr2,drho1ds1,drho1ds2 + real*8 drho1drm1(3),drho1drm2(3) + real*8 drho2dr1,drho2dr2,drho2ds1,drho2ds2 + real*8 drho2drm1(3),drho2drm2(3) + real*8 drho3dr1,drho3dr2,drho3ds1,drho3ds2 + real*8 drho3drm1(3),drho3drm2(3) + real*8 dt1dr1,dt1dr2,dt1ds1,dt1ds2 + real*8 dt2dr1,dt2dr2,dt2ds1,dt2ds2 + real*8 dt3dr1,dt3dr2,dt3ds1,dt3ds2 + real*8 drhodr1,drhodr2,drhods1,drhods2,drhodrm1(3),drhodrm2(3) + real*8 arg,arg1,arg2 + real*8 arg1i1,arg1j1,arg1i2,arg1j2,arg2i2,arg2j2 + real*8 arg1i3,arg1j3,arg2i3,arg2j3,arg3i3,arg3j3 + real*8 dsij1,dsij2,force1,force2 + real*8 t1i,t2i,t3i,t1j,t2j,t3j + + errorflag = 0 + + strscomp = 0 + +c Compute forces atom i + + elti = fmap(type(i)) + xitmp = x(1,i) + yitmp = x(2,i) + zitmp = x(3,i) + +c Treat each pair + do jn = 1,numneigh + + if (scrfcn(jn).ne.0.d0) then + j = firstneigh(jn) + + sij = scrfcn(jn)*fcpair(jn) + delij(1) = x(1,j) - xitmp + delij(2) = x(2,j) - yitmp + delij(3) = x(3,j) - zitmp + rij2 = delij(1)*delij(1) + delij(2)*delij(2) + $ + delij(3)*delij(3) + if (rij2.lt.cutforcesq) then + rij = sqrt(rij2) + r = rij + eltj = fmap(type(j)) + +c Compute phi and phip + ind = eltind(elti,eltj) + pp = rij*rdrar + 1.0D0 + kk = pp + kk = min(kk,nrar-1) + pp = pp - kk + pp = min(pp,1.0D0) + phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp + $ + phirar1(kk,ind))*pp + phirar(kk,ind) + phip = (phirar6(kk,ind)*pp + phirar5(kk,ind))*pp + $ + phirar4(kk,ind) + recip = 1.0d0/r + + if (eflag.eq.1) eng_vdwl = eng_vdwl + phi*sij + +c write(1,*) "force_meamf: phi: ",phi +c write(1,*) "force_meamf: phip: ",phip + +c Compute pair densities and derivatives + invrei = 1.d0/re_meam(elti,elti) + ai = rij*invrei - 1.d0 + ro0i = rho0_meam(elti) + rhoa0i = ro0i*exp(-beta0_meam(elti)*ai) + drhoa0i = -beta0_meam(elti)*invrei*rhoa0i + rhoa1i = ro0i*exp(-beta1_meam(elti)*ai) + drhoa1i = -beta1_meam(elti)*invrei*rhoa1i + rhoa2i = ro0i*exp(-beta2_meam(elti)*ai) + drhoa2i = -beta2_meam(elti)*invrei*rhoa2i + rhoa3i = ro0i*exp(-beta3_meam(elti)*ai) + drhoa3i = -beta3_meam(elti)*invrei*rhoa3i + + if (elti.ne.eltj) then + invrej = 1.d0/re_meam(eltj,eltj) + aj = rij*invrej - 1.d0 + ro0j = rho0_meam(eltj) + rhoa0j = ro0j*exp(-beta0_meam(eltj)*aj) + drhoa0j = -beta0_meam(eltj)*invrej*rhoa0j + rhoa1j = ro0j*exp(-beta1_meam(eltj)*aj) + drhoa1j = -beta1_meam(eltj)*invrej*rhoa1j + rhoa2j = ro0j*exp(-beta2_meam(eltj)*aj) + drhoa2j = -beta2_meam(eltj)*invrej*rhoa2j + rhoa3j = ro0j*exp(-beta3_meam(eltj)*aj) + drhoa3j = -beta3_meam(eltj)*invrej*rhoa3j + else + rhoa0j = rhoa0i + drhoa0j = drhoa0i + rhoa1j = rhoa1i + drhoa1j = drhoa1i + rhoa2j = rhoa2i + drhoa2j = drhoa2i + rhoa3j = rhoa3i + drhoa3j = drhoa3i + endif + + nv2 = 1 + nv3 = 1 + arg1i1 = 0.d0 + arg1j1 = 0.d0 + arg1i2 = 0.d0 + arg1j2 = 0.d0 + arg1i3 = 0.d0 + arg1j3 = 0.d0 + arg3i3 = 0.d0 + arg3j3 = 0.d0 + do n = 1,3 + do p = n,3 + do q = p,3 + arg = delij(n)*delij(p)*delij(q)*v3D(nv3) + arg1i3 = arg1i3 + Arho3(nv3,i)*arg + arg1j3 = arg1j3 - Arho3(nv3,j)*arg + nv3 = nv3+1 + enddo + arg = delij(n)*delij(p)*v2D(nv2) + arg1i2 = arg1i2 + Arho2(nv2,i)*arg + arg1j2 = arg1j2 + Arho2(nv2,j)*arg + nv2 = nv2+1 + enddo + arg1i1 = arg1i1 + Arho1(n,i)*delij(n) + arg1j1 = arg1j1 - Arho1(n,j)*delij(n) + arg3i3 = arg3i3 + Arho3b(n,i)*delij(n) + arg3j3 = arg3j3 - Arho3b(n,j)*delij(n) + enddo + +c rho0 terms + drho0dr1 = drhoa0j * sij + drho0dr2 = drhoa0i * sij + +c rho1 terms + a1 = 2*sij/rij + drho1dr1 = a1*(drhoa1j-rhoa1j/rij)*arg1i1 + drho1dr2 = a1*(drhoa1i-rhoa1i/rij)*arg1j1 + a1 = 2.d0*sij/rij + do m = 1,3 + drho1drm1(m) = a1*rhoa1j*Arho1(m,i) + drho1drm2(m) = -a1*rhoa1i*Arho1(m,j) + enddo + +c rho2 terms + a2 = 2*sij/rij2 + drho2dr1 = a2*(drhoa2j - 2*rhoa2j/rij)*arg1i2 + $ - 2.d0/3.d0*Arho2b(i)*drhoa2j*sij + drho2dr2 = a2*(drhoa2i - 2*rhoa2i/rij)*arg1j2 + $ - 2.d0/3.d0*Arho2b(j)*drhoa2i*sij + a2 = 4*sij/rij2 + do m = 1,3 + drho2drm1(m) = 0.d0 + drho2drm2(m) = 0.d0 + do n = 1,3 + drho2drm1(m) = drho2drm1(m) + $ + Arho2(vind2D(m,n),i)*delij(n) + drho2drm2(m) = drho2drm2(m) + $ - Arho2(vind2D(m,n),j)*delij(n) + enddo + drho2drm1(m) = a2*rhoa2j*drho2drm1(m) + drho2drm2(m) = -a2*rhoa2i*drho2drm2(m) + enddo + +c rho3 terms + rij3 = rij*rij2 + a3 = 2*sij/rij3 + a3a = 6.d0/5.d0*sij/rij + drho3dr1 = a3*(drhoa3j - 3*rhoa3j/rij)*arg1i3 + $ - a3a*(drhoa3j - rhoa3j/rij)*arg3i3 + drho3dr2 = a3*(drhoa3i - 3*rhoa3i/rij)*arg1j3 + $ - a3a*(drhoa3i - rhoa3i/rij)*arg3j3 + a3 = 6*sij/rij3 + a3a = 6*sij/(5*rij) + do m = 1,3 + drho3drm1(m) = 0.d0 + drho3drm2(m) = 0.d0 + nv2 = 1 + do n = 1,3 + do p = n,3 + arg = delij(n)*delij(p)*v2D(nv2) + drho3drm1(m) = drho3drm1(m) + $ + Arho3(vind3D(m,n,p),i)*arg + drho3drm2(m) = drho3drm2(m) + $ + Arho3(vind3D(m,n,p),j)*arg + nv2 = nv2 + 1 + enddo + enddo + drho3drm1(m) = (a3*drho3drm1(m) - a3a*Arho3b(m,i)) + $ *rhoa3j + drho3drm2(m) = (-a3*drho3drm2(m) + a3a*Arho3b(m,j)) + $ *rhoa3i + enddo + +c Compute derivatives of weighting functions t wrt rij + t1i = t_ave(1,i) + t2i = t_ave(2,i) + t3i = t_ave(3,i) + t1j = t_ave(1,j) + t2j = t_ave(2,j) + t3j = t_ave(3,j) + + ai = 0.d0 + if( rho0(i) .ne. 0.d0 ) then + ai = drhoa0j*sij/rho0(i) + end if + aj = 0.d0 + if( rho0(j) .ne. 0.d0 ) then + aj = drhoa0i*sij/rho0(j) + end if + + dt1dr1 = ai*(t1_meam(eltj)-t1i) + dt1dr2 = aj*(t1_meam(elti)-t1j) + dt2dr1 = ai*(t2_meam(eltj)-t2i) + dt2dr2 = aj*(t2_meam(elti)-t2j) + dt3dr1 = ai*(t3_meam(eltj)-t3i) + dt3dr2 = aj*(t3_meam(elti)-t3j) + +c Compute derivatives of total density wrt rij, sij and rij(3) + call get_shpfcn(shpi,lattce_meam(elti,elti)) + call get_shpfcn(shpj,lattce_meam(eltj,eltj)) + drhodr1 = dGamma1(i)*drho0dr1 + $ + dGamma2(i)* + $ (dt1dr1*rho1(i)+t1i*drho1dr1 + $ + dt2dr1*rho2(i)+t2i*drho2dr1 + $ + dt3dr1*rho3(i)+t3i*drho3dr1) + $ - dGamma3(i)* + $ (shpi(1)*dt1dr1+shpi(2)*dt2dr1+shpi(3)*dt3dr1) + drhodr2 = dGamma1(j)*drho0dr2 + $ + dGamma2(j)* + $ (dt1dr2*rho1(j)+t1j*drho1dr2 + $ + dt2dr2*rho2(j)+t2j*drho2dr2 + $ + dt3dr2*rho3(j)+t3j*drho3dr2) + $ - dGamma3(j)* + $ (shpj(1)*dt1dr2+shpj(2)*dt2dr2+shpj(3)*dt3dr2) + do m = 1,3 + drhodrm1(m) = 0.d0 + drhodrm2(m) = 0.d0 + drhodrm1(m) = dGamma2(i)* + $ (t1i*drho1drm1(m) + $ + t2i*drho2drm1(m) + $ + t3i*drho3drm1(m)) + drhodrm2(m) = dGamma2(j)* + $ (t1j*drho1drm2(m) + $ + t2j*drho2drm2(m) + $ + t3j*drho3drm2(m)) + enddo + +c Compute derivatives wrt sij, but only if necessary + if (dscrfcn(jn).ne.0.d0) then + drho0ds1 = rhoa0j + drho0ds2 = rhoa0i + a1 = 2.d0/rij + drho1ds1 = a1*rhoa1j*arg1i1 + drho1ds2 = a1*rhoa1i*arg1j1 + a2 = 2.d0/rij2 + drho2ds1 = a2*rhoa2j*arg1i2 + $ - 2.d0/3.d0*Arho2b(i)*rhoa2j + drho2ds2 = a2*rhoa2i*arg1j2 + $ - 2.d0/3.d0*Arho2b(j)*rhoa2i + a3 = 2.d0/rij3 + a3a = 6.d0/(5.d0*rij) + drho3ds1 = a3*rhoa3j*arg1i3 - a3a*rhoa3j*arg3i3 + drho3ds2 = a3*rhoa3i*arg1j3 - a3a*rhoa3i*arg3j3 + + ai = 0.d0 + if( rho0(i) .ne. 0.d0 ) then + ai = rhoa0j/rho0(i) + end if + aj = 0.d0 + if( rho0(j) .ne. 0.d0 ) then + aj = rhoa0i/rho0(j) + end if + + dt1ds1 = ai*(t1_meam(eltj)-t1i) + dt1ds2 = aj*(t1_meam(elti)-t1j) + dt2ds1 = ai*(t2_meam(eltj)-t2i) + dt2ds2 = aj*(t2_meam(elti)-t2j) + dt3ds1 = ai*(t3_meam(eltj)-t3i) + dt3ds2 = aj*(t3_meam(elti)-t3j) + drhods1 = dGamma1(i)*drho0ds1 + $ + dGamma2(i)* + $ (dt1ds1*rho1(i)+t1i*drho1ds1 + $ + dt2ds1*rho2(i)+t2i*drho2ds1 + $ + dt3ds1*rho3(i)+t3i*drho3ds1) + $ - dGamma3(i)* + $ (shpi(1)*dt1ds1+shpi(2)*dt2ds1+shpi(3)*dt3ds1) + drhods2 = dGamma1(j)*drho0ds2 + $ + dGamma2(j)* + $ (dt1ds2*rho1(j)+t1j*drho1ds2 + $ + dt2ds2*rho2(j)+t2j*drho2ds2 + $ + dt3ds2*rho3(j)+t3j*drho3ds2) + $ - dGamma3(j)* + $ (shpj(1)*dt1ds2+shpj(2)*dt2ds2+shpj(3)*dt3ds2) + endif + +c Compute derivatives of energy wrt rij, sij and rij(3) + dUdrij = phip*sij + $ + fp(i)*drhodr1 + fp(j)*drhodr2 + dUdsij = 0.d0 + if (dscrfcn(jn).ne.0.d0) then + dUdsij = phi + $ + fp(i)*drhods1 + fp(j)*drhods2 + endif + do m = 1,3 + dUdrijm(m) = fp(i)*drhodrm1(m) + fp(j)*drhodrm2(m) + enddo + +c Add the part of the force due to dUdrij and dUdsij + force = dUdrij*recip + dUdsij*dscrfcn(jn) + do m = 1,3 + forcem = delij(m)*force + dUdrijm(m) + f(m,i) = f(m,i) + forcem + f(m,j) = f(m,j) - forcem + if (strscomp>0) then + do n = 1,3 + arg = delij(n)*forcem + strssa(n,m,i) = strssa(n,m,i) + arg + strssa(n,m,j) = strssa(n,m,j) + arg + enddo + endif + enddo +c write(1,*)"forcei1: ",f(1,i)," ",f(2,i)," ",f(3,i) +c write(1,*)"forcej1: ",f(1,j)," ",f(2,j)," ",f(3,j) + +c Now compute forces on other atoms k due to change in sij + if (sij.eq.0.d0.or.sij.eq.1.d0) goto 100 + do kn = 1,numneigh_full + k = firstneigh_full(kn) + if (k.ne.j) then + call dsij(i,j,k,jn,nmax,numneigh,rij2,dsij1,dsij2, + $ ntype,type,fmap,x,scrfcn,fcpair) + if (dsij1.ne.0.d0.or.dsij2.ne.0.d0) then + force1 = dUdsij*dsij1 + force2 = dUdsij*dsij2 + do m = 1,3 + delik(m) = x(m,k) - x(m,i) + deljk(m) = x(m,k) - x(m,j) + enddo + do m = 1,3 + f(m,i) = f(m,i) + force1*delik(m) + f(m,j) = f(m,j) + force2*deljk(m) + f(m,k) = f(m,k) - force1*delik(m) + $ - force2*deljk(m) + if (strscomp>0) then + do n = m,3 + arg1 = force1*delik(m)*delik(n) + arg2 = force2*deljk(m)*deljk(n) + strssa(m,n,i) = strssa(m,n,i) + arg1 + strssa(m,n,j) = strssa(m,n,j) + arg2 + strssa(m,n,k) = strssa(m,n,k) + + $ arg1 + arg2 + if (n.ne.m) then + strssa(n,m,i) = strssa(n,m,i) + arg1 + strssa(n,m,j) = strssa(n,m,j) + arg2 + strssa(n,m,k) = strssa(n,m,k) + $ + arg1 + arg2 + endif + enddo + endif + enddo + + endif + endif +c end of k loop + enddo + endif + 100 continue + endif +c end of j loop + enddo + + return + end diff --git a/lib/meam/meam_setup_done.F b/lib/meam/meam_setup_done.F new file mode 100755 index 0000000000..fe40bb45d6 --- /dev/null +++ b/lib/meam/meam_setup_done.F @@ -0,0 +1,733 @@ +c Declaration in pair_meam.h: +c +c void meam_setup_done(double *) +c +c Call from pair_meam.cpp: +c +c meam_setup_done(&cutmax) +c + + subroutine meam_setup_done(cutmax) + use meam_data + implicit none + + real*8 cutmax + + integer nv2, nv3, m, n, p + +c Force cutoff + cutforce = rc_meam + cutforcesq = cutforce*cutforce + +c Pass cutoff back to calling program + cutmax = cutforce + +c Augment t1 term + t1_meam(:) = t1_meam(:) + augt1 * 3.d0/5.d0 * t3_meam(:) + +c Compute off-diagonal alloy parameters + call alloyparams + +c indices and factors for Voight notation + nv2 = 1 + nv3 = 1 + do m = 1,3 + do n = m,3 + vind2D(m,n) = nv2 + vind2D(n,m) = nv2 + nv2 = nv2+1 + do p = n,3 + vind3D(m,n,p) = nv3 + vind3D(m,p,n) = nv3 + vind3D(n,m,p) = nv3 + vind3D(n,p,m) = nv3 + vind3D(p,m,n) = nv3 + vind3D(p,n,m) = nv3 + nv3 = nv3+1 + enddo + enddo + enddo + + v2D(1) = 1 + v2D(2) = 2 + v2D(3) = 2 + v2D(4) = 1 + v2D(5) = 2 + v2D(6) = 1 + + v3D(1) = 1 + v3D(2) = 3 + v3D(3) = 3 + v3D(4) = 3 + v3D(5) = 6 + v3D(6) = 3 + v3D(7) = 1 + v3D(8) = 3 + v3D(9) = 3 + v3D(10) = 1 + + nv2 = 1 + do m = 1,neltypes + do n = m,neltypes + eltind(m,n) = nv2 + eltind(n,m) = nv2 + nv2 = nv2+1 + enddo + enddo + +c Compute pair potentials and setup arrays for interpolation + nr = 1000 + dr = 1.1*rc_meam/nr + call compute_pair_meam + + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc +c Fill off-diagonal alloy parameters + subroutine alloyparams + use meam_data + implicit none + integer i,j,k + real*8 eb + +c Loop over pairs + do i = 1,neltypes + do j = 1,neltypes +c Treat off-diagonal pairs +c If i>j, set all equal to ij, set equal to the iebound, +c atom k definitely lies outside the screening function ellipse (so +c there is no need to calculate its effects). Here, compute it for all +c triplets (i,j,k) so that ebound(i,j) is the maximized over k + do i = 1,neltypes + do j = 1,neltypes + do k = 1,neltypes + eb = (Cmax_meam(i,j,k)*Cmax_meam(i,j,k)) + $ /(4.d0*(Cmax_meam(i,j,k)-1.d0)) + ebound_meam(i,j) = max(ebound_meam(i,j),eb) + enddo + enddo + enddo + + return + end + +c----------------------------------------------------------------------- +c compute MEAM pair potential for each pair of element types +c + + subroutine compute_pair_meam + use meam_data + implicit none + + real*8 r, temp + integer j,a,b,nv2 + real*8 astar,frac,phizbl + integer n,nmax,Z1,Z2 + real*8 arat,scrn + real*8 phitmp + + real*8, external :: phi_meam + real*8, external :: zbl + + +c allocate memory for array that defines the potential + allocate(phir(nr,(neltypes*(neltypes+1))/2)) + +c allocate coeff memory + + allocate(phirar(nr,(neltypes*(neltypes+1))/2)) + allocate(phirar1(nr,(neltypes*(neltypes+1))/2)) + allocate(phirar2(nr,(neltypes*(neltypes+1))/2)) + allocate(phirar3(nr,(neltypes*(neltypes+1))/2)) + allocate(phirar4(nr,(neltypes*(neltypes+1))/2)) + allocate(phirar5(nr,(neltypes*(neltypes+1))/2)) + allocate(phirar6(nr,(neltypes*(neltypes+1))/2)) + +c loop over pairs of element types + nv2 = 0 + do a = 1,neltypes + do b = a,neltypes + nv2 = nv2 + 1 + +c loop over r values and compute + do j = 1,nr + + r = (j-1)*dr + + phir(j,nv2) = phi_meam(r,a,b) + +c if using second-nearest neighbor, solve recursive problem +c (see Lee and Baskes, PRB 62(13):8564 eqn.(21)) + if (nn2_meam(a,b).eq.1) then + if (a.ne.b) then + write(6,*) 'Second-neighbor currently only valid ' + write(6,*) 'if element types are the same.' +c call error(' ') + endif + call get_Zij(Z1,lattce_meam(a,b)) + call get_Zij2(Z2,arat,scrn,lattce_meam(a,b), + $ Cmin_meam(a,a,b),Cmax_meam(a,a,b)) + + nmax = 10 + do n = 1,nmax + phir(j,nv2) = phir(j,nv2) + + $ (-Z2*scrn/Z1)**n * phi_meam(r*arat**n,a,b) + enddo + + endif + +c if astar <= -3 +c potential is zbl potential +c else if -3 < astar < -1 +c potential is linear combination with zbl potential +c endif + astar = alpha_meam(a,b) * (r/re_meam(a,b) - 1.d0) + if (astar.le.-3.d0) then + phir(j,nv2) = zbl(r,ielt_meam(a),ielt_meam(b)) + else if (astar.gt.-3.d0.and.astar.lt.-1.d0) then + call fcut(1-(astar+1.d0)/(-3.d0+1.d0),frac) + phizbl = zbl(r,ielt_meam(a),ielt_meam(b)) + phir(j,nv2) = frac*phir(j,nv2) + (1-frac)*phizbl + endif + + enddo + +c call interpolation + call interpolate_meam(nv2) + + enddo + enddo + + return + end + + +c----------------------------------------------------------------------c +c Compute MEAM pair potential for distance r, element types a and b +c + real*8 recursive function phi_meam(r,a,b)result(phi_m) + use meam_data + implicit none + + + integer a,b + real*8 r + real*8 a1,a2,a12 + real*8 t11av,t21av,t31av,t12av,t22av,t32av + real*8 G1,G2,s1(3),s2(3),s12(3),rho0_1,rho0_2 + real*8 Gam1,Gam2,Z1,Z2 + real*8 rhobar1,rhobar2,F1,F2 + real*8 rhoa01,rhoa11,rhoa21,rhoa31 + real*8 rhoa02,rhoa12,rhoa22,rhoa32 + real*8 rho01,rho11,rho21,rho31 + real*8 rho02,rho12,rho22,rho32 + real*8 scalfac,phiaa,phibb + real*8 Eu + integer Z12 + character*3 latta,lattb + + real*8, external :: erose + +c Equation numbers below refer to: +c I. Huang et.al., Modelling simul. Mater. Sci. Eng. 3:615 + +c get number of neighbors in the reference structure +c Nref(i,j) = # of i's neighbors of type j + call get_Zij(Z12,lattce_meam(a,b)) + + call get_densref(r,a,b,rho01,rho11,rho21,rho31, + $ rho02,rho12,rho22,rho32) + +c calculate average weighting factors for the reference structure + if (lattce_meam(a,b).eq.'c11') then + scalfac = 1.0/(rho01+rho02) + t11av = scalfac*(t1_meam(a)*rho01 + t1_meam(b)*rho02) + t12av = t11av + t21av = scalfac*(t2_meam(a)*rho01 + t2_meam(b)*rho02) + t22av = t21av + t31av = scalfac*(t3_meam(a)*rho01 + t3_meam(b)*rho02) + t32av = t31av + else +c average weighting factors for the reference structure, eqn. I.8 + call get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, + $ t1_meam(a),t2_meam(a),t3_meam(a), + $ t1_meam(b),t2_meam(b),t3_meam(b),lattce_meam(a,b)) + endif + +c for c11b structure, calculate background electron densities + if (lattce_meam(a,b).eq.'c11') then + latta = lattce_meam(a,a) + if (latta.eq.'dia') then + rhobar1 = ((Z12/2)*(rho02+rho01))**2 + + $ t11av*(rho12-rho11)**2 + + $ t21av/6.0*(rho22+rho21)**2 + + $ 121.0/40.*t31av*(rho32-rho31)**2 + rhobar1 = sqrt(rhobar1) + rhobar2 = (Z12*rho01)**2 + 2.0/3.0*t21av*rho21**2 + rhobar2 = sqrt(rhobar2) + else + rhobar2 = ((Z12/2)*(rho01+rho02))**2 + + $ t12av*(rho11-rho12)**2 + + $ t22av/6.0*(rho21+rho22)**2 + + $ 121.0/40.*t32av*(rho31-rho32)**2 + rhobar2 = sqrt(rhobar2) + rhobar1 = (Z12*rho02)**2 + 2.0/3.0*t22av*rho22**2 + rhobar1 = sqrt(rhobar1) + endif + else +c for other structures, use formalism developed in Huang's paper +c +c composition-dependent scaling, equation I.7 +c -- note: The shape factors for the individual elements are used, since +c this function should cancel the F(rho) terms when the atoms are +c in the reference structure. + Z1 = Z_meam(a) + Z2 = Z_meam(b) + if (ibar_meam(a).le.0) then + G1 = 1.d0 + else + call get_shpfcn(s1,lattce_meam(a,a)) + Gam1 = (s1(1)*t11av+s1(2)*t21av+s1(3)*t31av)/(Z1*Z1) + call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1) + endif + if (ibar_meam(b).le.0) then + G2 = 1.d0 + else + call get_shpfcn(s2,lattce_meam(b,b)) + Gam2 = (s2(1)*t12av+s2(2)*t22av+s2(3)*t32av)/(Z2*Z2) + call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2) + endif + rho0_1 = rho0_meam(a)*Z1*G1 + rho0_2 = rho0_meam(b)*Z2*G2 + +c get number of neighbors in the reference structure +c Nref(i,j) = # of i's neighbors of type j +c call get_Zij(Z12,lattce_meam(a,b)) +c +c call get_densref(r,a,b,rho01,rho11,rho21,rho31, +c $ rho02,rho12,rho22,rho32) + +c compute total background density, eqn I.6 + Gam1 = (t11av*rho11+t21av*rho21+t31av*rho31) + Gam1 = Gam1/(rho01*rho01) + Gam2 = (t12av*rho12+t22av*rho22+t32av*rho32) + Gam2 = Gam2/(rho02*rho02) + call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1) + call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2) + rhobar1 = rho01/rho0_1*G1 + rhobar2 = rho02/rho0_2*G2 + endif +c compute embedding functions, eqn I.5 + if (rhobar1.eq.0.d0) then + F1 = 0.d0 + else + F1 = A_meam(a)*Ec_meam(a,a)*rhobar1*log(rhobar1) + endif + if (rhobar2.eq.0.d0) then + F2 = 0.d0 + else + F2 = A_meam(b)*Ec_meam(b,b)*rhobar2*log(rhobar2) + endif + +c compute Rose function, I.16 + Eu = erose(r,re_meam(a,b),alpha_meam(a,b), + $ Ec_meam(a,b),repuls_meam(a,b),attrac_meam(a,b)) + +c calculate the pair energy + if (lattce_meam(a,b).eq.'c11') then + latta = lattce_meam(a,a) + if (latta.eq.'dia') then + phiaa = phi_meam(r,a,a) + phi_m = (3*Eu - F2 - 2*F1 - 5*phiaa)/Z12 + else + phibb = phi_meam(r,b,b) + phi_m = (3*Eu - F1 - 2*F2 - 5*phibb)/Z12 + endif +c phi_m = 0.d0 + else +c +c potential is computed from Rose function and embedding energy + phi_m = (2*Eu - F1 - F2)/Z12 +c + endif + +c if r = 0, just return 0 + if (r.eq.0.d0) then + phi_m = 0.d0 + endif + + return + end + +c----------------------------------------------------------------------c +c Shape factors for various configurations + subroutine get_shpfcn(s,latt) + implicit none + real*8 s(3) + character*3 latt + if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'b1') then + s(1) = 0.d0 + s(2) = 0.d0 + s(3) = 0.d0 + else if (latt.eq.'hcp') then + s(1) = 0.d0 + s(2) = 0.d0 + s(3) = 1.d0/3.d0 + else if (latt.eq.'dia') then + s(1) = 0.d0 + s(2) = 0.d0 + s(3) = 32.d0/9.d0 + else if (latt.eq.'dim') then + s(1) = 1.d0 + s(2) = 2.d0/3.d0 +c s(3) = 1.d0 + s(3) = 0.4d0 + else + s(1) = 0.0 +c call error('Lattice not defined in get_shpfcn.') + endif + return + end +c------------------------------------------------------------------------------c +c Average weighting factors for the reference structure + subroutine get_tavref(t11av,t21av,t31av,t12av,t22av,t32av, + $ t11,t21,t31,t12,t22,t32,latt) + implicit none + real*8 t11av,t21av,t31av,t12av,t22av,t32av + real*8 t11,t21,t31,t12,t22,t32 + character*3 latt + if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'dia' + $ .or.latt.eq.'hcp'.or.latt.eq.'b1' + $ .or.latt.eq.'dim') then +c all neighbors are of the opposite type + t11av = t12 + t21av = t22 + t31av = t32 + t12av = t11 + t22av = t21 + t32av = t31 + else +c call error('Lattice not defined in get_tavref.') + endif + return + end +c------------------------------------------------------------------------------c +c Number of neighbors for the reference structure + subroutine get_Zij(Zij,latt) + implicit none + integer Zij + character*3 latt + if (latt.eq.'fcc') then + Zij = 12 + else if (latt.eq.'bcc') then + Zij = 8 + else if (latt.eq.'hcp') then + Zij = 12 + else if (latt.eq.'b1') then + Zij = 6 + else if (latt.eq.'dia') then + Zij = 4 + else if (latt.eq.'dim') then + Zij = 1 + else if (latt.eq.'c11') then + Zij = 10 + else +c call error('Lattice not defined in get_Zij.') + endif + return + end + +c------------------------------------------------------------------------------c +c Zij2 = number of second neighbors, a = distance ratio R1/R2, and S = second +c neighbor screening function for lattice type "latt" + + subroutine get_Zij2(Zij2,a,S,latt,cmin,cmax) + implicit none + integer Zij2 + real*8 a,S,cmin,cmax + character*3 latt + real*8 rratio,C,x,sijk + integer numscr + + if (latt.eq.'bcc') then + Zij2 = 6 + a = 2.d0/sqrt(3.d0) + numscr = 4 + else if (latt.eq.'fcc') then + Zij2 = 6 + a = sqrt(2.d0) + numscr = 4 + else if (latt.eq.'dia') then + Zij2 = 0 + a = sqrt(8.d0/3.d0) + numscr = 4 + if (cmin.lt.0.500001) then +c call error('can not do 2NN MEAM for dia') + endif + else if (latt.eq.'hcp') then + Zij2 = 6 + a = sqrt(2.d0) + numscr = 4 + else if (latt.eq.'b1') then + Zij2 = 6 + a = sqrt(2.d0) + numscr = 4 + else +c call error('Lattice not defined in get_Zij2.') + endif + +c Compute screening for each first neighbor + C = 4.d0/(a*a) - 1.d0 + x = (C-cmin)/(cmax-cmin) + call fcut(x,sijk) +c There are numscr first neighbors screening the second neighbors + S = sijk**numscr + + return + end + +c------------------------------------------------------------------------------c +c Calculate density functions, assuming reference configuration + subroutine get_densref(r,a,b,rho01,rho11,rho21,rho31, + $ rho02,rho12,rho22,rho32) + use meam_data + implicit none + real*8 r,rho01,rho11,rho21,rho31,rho02,rho12,rho22,rho32 + real*8 a1,a2 + real*8 rhoa01,rhoa11,rhoa21,rhoa31,rhoa02,rhoa12,rhoa22,rhoa32 + real*8 C,Cmin,Cmax,fc,s(3) + character*3 lat + integer a,b + integer Zij1nn,Zij2nn + real*8 rhoa01nn,rhoa02nn + real*8 arat,scrn + + a1 = r/re_meam(a,a) - 1.d0 + a2 = r/re_meam(b,b) - 1.d0 + + rhoa01 = rho0_meam(a)*exp(-beta0_meam(a)*a1) + rhoa11 = rho0_meam(a)*exp(-beta1_meam(a)*a1) + rhoa21 = rho0_meam(a)*exp(-beta2_meam(a)*a1) + rhoa31 = rho0_meam(a)*exp(-beta3_meam(a)*a1) + rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2) + rhoa12 = rho0_meam(b)*exp(-beta1_meam(b)*a2) + rhoa22 = rho0_meam(b)*exp(-beta2_meam(b)*a2) + rhoa32 = rho0_meam(b)*exp(-beta3_meam(b)*a2) + + lat = lattce_meam(a,b) + + rho11 = 0.d0 + rho21 = 0.d0 + rho31 = 0.d0 + rho12 = 0.d0 + rho22 = 0.d0 + rho32 = 0.d0 + + call get_Zij(Zij1nn,lat) + + if (lat.eq.'fcc') then + rho01 = 12.d0*rhoa02 + rho02 = 12.d0*rhoa01 + else if (lat.eq.'bcc') then + rho01 = 8.d0*rhoa02 + rho02 = 8.d0*rhoa01 + else if (lat.eq.'b1') then +c Warning: this assumes that atoms of element 1 are not completely +c screened from each other, but that atoms of element 2 are. So +c Cmin_meam(1,1,2) can be less than 1.0, but Cmin_meam(2,2,1) cannot. +c This should be made more general in the future. +c Cmin = Cmin_meam(1,1,2) +c Cmax = Cmax_meam(1,1,2) +c C = (1.0-Cmin)/(Cmax-Cmin) +c call fcut(C,fc) +c a1 = sqrt(2.0)*r/re_meam(1,1) - 1.d0 +c rho01 = 6*rhoa02 + 12*fc*fc*rho0_meam(1)*exp(-beta0_meam(1)*a1) +c rho02 = 6*rhoa01 + rho01 = 6*rhoa02 + rho02 = 6*rhoa01 + else if (lat.eq.'dia') then + rho01 = 4*rhoa02 + rho02 = 4*rhoa01 + rho31 = 32.d0/9.d0*rhoa32*rhoa32 + rho32 = 32.d0/9.d0*rhoa31*rhoa31 + else if (lat.eq.'hcp') then + rho01 = 12*rhoa02 + rho02 = 12*rhoa01 + rho31 = 1.d0/3.d0*rhoa32*rhoa32 + rho32 = 1.d0/3.d0*rhoa31*rhoa31 + else if (lat.eq.'dim') then + call get_shpfcn(s,'dim') + rho01 = rhoa02 + rho02 = rhoa01 + rho11 = s(1)*rhoa12*rhoa12 + rho12 = s(1)*rhoa11*rhoa11 + rho21 = s(2)*rhoa22*rhoa22 + rho22 = s(2)*rhoa21*rhoa21 + rho31 = s(3)*rhoa32*rhoa32 + rho32 = s(3)*rhoa31*rhoa31 + else if (lat.eq.'c11') then + rho01 = rhoa01 + rho02 = rhoa02 + rho11 = rhoa11 + rho12 = rhoa12 + rho21 = rhoa21 + rho22 = rhoa22 + rho31 = rhoa31 + rho32 = rhoa32 + else +c call error('Lattice not defined in get_densref.') + endif + + if (nn2_meam(a,b).eq.1) then +c Assume that second neighbor is of same type, first neighbor +c may be of different type + + if (lat.ne.'bcc') then +c call error('Second-neighbor not defined for this lattice') + endif + + call get_Zij2(Zij2nn,arat,scrn,lat, + $ Cmin_meam(a,a,b),Cmax_meam(a,a,b)) + + a1 = arat*r/re_meam(a,a) - 1.d0 + a2 = arat*r/re_meam(b,b) - 1.d0 + + rhoa01nn = rho0_meam(a)*exp(-beta0_meam(a)*a1) + rhoa02nn = rho0_meam(b)*exp(-beta0_meam(b)*a2) + + rho01 = rho01 + Zij2nn*scrn*rhoa02nn + rho02 = rho02 + Zij2nn*scrn*rhoa01nn + + endif + + return + end + +c--------------------------------------------------------------------- +c Compute ZBL potential +c + real*8 function zbl(r,z1,z2) + implicit none + integer i,z1,z2 + real*8 r,c,d,a,azero,cc,x + dimension c(4),d(4) + data c /0.028171,0.28022,0.50986,0.18175/ + data d /0.20162,0.40290,0.94229,3.1998/ + data azero /0.4685/ + data cc /14.3997/ +c azero = (9pi^2/128)^1/3 (0.529) Angstroms + a = azero/(z1**0.23+z2**0.23) + zbl = 0.0 + x = r/a + do i=1,4 + zbl = zbl + c(i)*exp(-d(i)*x) + enddo + if (r.gt.0.d0) zbl = zbl*z1*z2/r*cc + return + end + +c--------------------------------------------------------------------- +c Compute Rose energy function, I.16 +c + real*8 function erose(r,re,alpha,Ec,repuls,attrac) + implicit none + real*8 r,re,alpha,Ec,repuls,attrac,astar,a3 + + erose = 0.d0 + + if (r.gt.0.d0) then + astar = alpha * (r/re - 1.d0) + a3 = 0.d0 + if (astar.ge.0) then + a3 = attrac + else if (astar.lt.0) then + a3 = repuls + endif +c erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar) + erose = -Ec*(1+astar+(-attrac+repuls/r)*(astar**3))*exp(-astar) + endif + + return + end + +c ----------------------------------------------------------------------- + + subroutine interpolate_meam(ind) + use meam_data + implicit none + + integer j,ind + real*8 drar + +c map to coefficient space + + nrar = nr + drar = dr + rdrar = 1.0D0/drar + +c phir interp + do j = 1,nrar + phirar(j,ind) = phir(j,ind) + enddo + + phirar1(1,ind) = phirar(2,ind)-phirar(1,ind) + phirar1(2,ind) = 0.5D0*(phirar(3,ind)-phirar(1,ind)) + phirar1(nrar-1,ind) = 0.5D0*(phirar(nrar,ind) + $ -phirar(nrar-2,ind)) + phirar1(nrar,ind) = 0.0D0 + do j = 3,nrar-2 + phirar1(j,ind) = ((phirar(j-2,ind)-phirar(j+2,ind)) + + $ 8.0D0*(phirar(j+1,ind)-phirar(j-1,ind)))/12. + enddo + + do j = 1,nrar-1 + phirar2(j,ind) = 3.0D0*(phirar(j+1,ind)-phirar(j,ind)) - + $ 2.0D0*phirar1(j,ind) - phirar1(j+1,ind) + phirar3(j,ind) = phirar1(j,ind) + phirar1(j+1,ind) - + $ 2.0D0*(phirar(j+1,ind)-phirar(j,ind)) + enddo + phirar2(nrar,ind) = 0.0D0 + phirar3(nrar,ind) = 0.0D0 + + do j = 1,nrar + phirar4(j,ind) = phirar1(j,ind)/drar + phirar5(j,ind) = 2.0D0*phirar2(j,ind)/drar + phirar6(j,ind) = 3.0D0*phirar3(j,ind)/drar + enddo + + end diff --git a/lib/meam/meam_setup_global.F b/lib/meam/meam_setup_global.F new file mode 100755 index 0000000000..1223e34b22 --- /dev/null +++ b/lib/meam/meam_setup_global.F @@ -0,0 +1,105 @@ +c +c declaration in pair_meam.h: +c +c void meam_setup_global(int *, int *, double *, int *, double *, double *, +c double *, double *, double *, double *, double *, +c double *, double *, double *, double *, double *, +c double *, double *, int *); +c +c call in pair_meam.cpp: +c +c meam_setup_global(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3, +c alat,esub,asub,t0,t1,t2,t3,rozero,ibar); +c +c + + subroutine meam_setup_global(nelt, lat, z, ielement, atwt, alpha, + $ b0, b1, b2, b3, alat, esub, asub, + $ t0, t1, t2, t3, rozero, ibar) + + use meam_data + implicit none + + integer nelt, lat, ielement, ibar + real*8 z, atwt, alpha, b0, b1, b2, b3 + real*8 alat, esub, asub, t0, t1, t2, t3 + real*8 rozero + + dimension lat(nelt), ielement(nelt), ibar(nelt) + dimension z(nelt), atwt(nelt), alpha(nelt) + dimension b0(nelt), b1(nelt), b2(nelt), b3(nelt) + dimension alat(nelt), esub(nelt), asub(nelt) + dimension t0(nelt), t1(nelt), t2(nelt), t3(nelt), rozero(nelt) + + integer i + real*8 tmplat(maxelt) + + neltypes = nelt + + do i = 1,nelt + + if (lat(i).eq.0) then + lattce_meam(i,i) = 'fcc' + else if (lat(i).eq.1) then + lattce_meam(i,i) = 'bcc' + else if (lat(i).eq.2) then + lattce_meam(i,i) = 'hcp' + else if (lat(i).eq.3) then + lattce_meam(i,i) = 'dim' + else if (lat(i).eq.4) then + lattce_meam(i,i) = 'dia' + else +c unknown + endif + + Z_meam(i) = z(i) + ielt_meam(i) = ielement(i) + alpha_meam(i,i) = alpha(i) + beta0_meam(i) = b0(i) + beta1_meam(i) = b1(i) + beta2_meam(i) = b2(i) + beta3_meam(i) = b3(i) + tmplat(i) = alat(i) + Ec_meam(i,i) = esub(i) + A_meam(i) = asub(i) + t0_meam(i) = t0(i) + t1_meam(i) = t1(i) + t2_meam(i) = t2(i) + t3_meam(i) = t3(i) + rho0_meam(i) = rozero(i) + ibar_meam(i) = ibar(i) + + if (lattce_meam(i,i).eq.'fcc') then + re_meam(i,i) = tmplat(i)/sqrt(2.d0) + elseif (lattce_meam(i,i).eq.'bcc') then + re_meam(i,i) = tmplat(i)*sqrt(3.d0)/2.d0 + elseif (lattce_meam(i,i).eq.'hcp') then + re_meam(i,i) = tmplat(i) + elseif (lattce_meam(i,i).eq.'dim') then + re_meam(i,i) = tmplat(i) + elseif (lattce_meam(i,i).eq.'dia') then + re_meam(i,i) = tmplat(i)*sqrt(3.d0)/4.d0 + else +c error + endif + + enddo + + +c Set some defaults + rc_meam = 4.0 + delr_meam = 0.1 + attrac_meam(:,:) = 0.0 + repuls_meam(:,:) = 0.0 + Cmax_meam(:,:,:) = 2.8 + Cmin_meam(:,:,:) = 2.0 + ebound_meam(:,:) = (2.8d0**2)/(4.d0*(2.8d0-1.d0)) + delta_meam(:,:) = 0.0 + nn2_meam(:,:) = 0 + gsmooth_factor = 99.0 + augt1 = 1 + + return + end + + diff --git a/lib/meam/meam_setup_param.F b/lib/meam/meam_setup_param.F new file mode 100755 index 0000000000..7a26ba71dd --- /dev/null +++ b/lib/meam/meam_setup_param.F @@ -0,0 +1,121 @@ +c +c Declaration in pair_meam.h: +c +c void meam_setup_param(int *, double *, int *, int *, int *); +c +c Call in pair_meam.cpp +c +c meam_setup_param(&which,&value,&nindex,index,&errorflag); +c +c +c +c The "which" argument corresponds to the index of the "keyword" array +c in pair_meam.cpp: +c +c 0 = Ec_meam +c 1 = alpha_meam +c 2 = rho0_meam +c 3 = delta_meam +c 4 = lattce_meam +c 5 = attrac_meam +c 6 = repuls_meam +c 7 = nn2_meam +c 8 = Cmin_meam +c 9 = Cmax_meam +c 10 = rc_meam +c 11 = delr_meam +c 12 = augt1 +c 13 = gsmooth_factor +c 14 = re_meam + + subroutine meam_setup_param(which, value, nindex, + $ index, errorflag) + + use meam_data + implicit none + + integer which, nindex, index(3), errorflag + real*8 value + + errorflag = 0 + +c 0 = Ec_meam + if (which.eq.0) then + Ec_meam(index(1),index(2)) = value + +c 1 = alpha_meam + else if (which.eq.1) then + alpha_meam(index(1),index(2)) = value + +c 2 = rho0_meam + else if (which.eq.2) then + rho0_meam(index(1)) = value + +c 3 = delta_meam + else if (which.eq.3) then + delta_meam(index(1),index(2)) = value + +c 4 = lattce_meam + else if (which.eq.4) then + if (value.eq.0) then + lattce_meam(index(1),index(2)) = "fcc" + else if (value.eq.1) then + lattce_meam(index(1),index(2)) = "bcc" + else if (value.eq.2) then + lattce_meam(index(1),index(2)) = "hcp" + else if (value.eq.3) then + lattce_meam(index(1),index(2)) = "dim" + else if (value.eq.4) then + lattce_meam(index(1),index(2)) = "dia" + else if (value.eq.5) then + lattce_meam(index(1),index(2)) = 'b1' + else if (value.eq.6) then + lattce_meam(index(1),index(2)) = 'c11' + endif + +c 5 = attrac_meam + else if (which.eq.5) then + attrac_meam(index(1),index(2)) = value + +c 6 = repuls_meam + else if (which.eq.6) then + repuls_meam(index(1),index(2)) = value + +c 7 = nn2_meam + else if (which.eq.7) then + nn2_meam(index(1),index(2)) = value + +c 8 = Cmin_meam + else if (which.eq.8) then + Cmin_meam(index(1),index(2),index(3)) = value + +c 9 = Cmax_meam + else if (which.eq.9) then + Cmax_meam(index(1),index(2),index(3)) = value + +c 10 = rc_meam + else if (which.eq.10) then + rc_meam = value + +c 11 = delr_meam + else if (which.eq.11) then + delr_meam = value + +c 12 = augt1 + else if (which.eq.12) then + augt1 = value + +c 13 = gsmooth + else if (which.eq.13) then + gsmooth_factor = value + +c 14 = re_meam + else if (which.eq.14) then + re_meam(index(1),index(2)) = value + + else + errorflag = 1 + endif + + return + end