diff --git a/lib/meam/meam_data.F b/lib/meam/meam_data.F index 389b2c64e..d94836392 100755 --- a/lib/meam/meam_data.F +++ b/lib/meam/meam_data.F @@ -1,83 +1,86 @@ 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 rho_ref_meam = background density for reference structure 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 zbl_meam = 1 if zbl potential for small r to be use, else 0 +c emb_lin_neg = 1 if linear embedding function for rhob to be used, else 0 +c bkgd_dyn = 1 if reference densities follows Dynamo, 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 ialloy = flag for newer alloy formulation (as in dynamo code) c mix_ref_t = flag to recover "old" way of computing t in reference config c erose_form = selection parameter for form of E_rose function 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) real*8 rho_ref_meam(maxelt) integer ibar_meam(maxelt),ielt_meam(maxelt) character*3 lattce_meam(maxelt,maxelt) integer nn2_meam(maxelt,maxelt) integer zbl_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, ialloy, mix_ref_t, erose_form + integer emb_lin_neg, bkgd_dyn 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 index 35378a83d..2d5bcef36 100755 --- a/lib/meam/meam_dens_final.F +++ b/lib/meam/meam_dens_final.F @@ -1,246 +1,282 @@ c Extern "C" declaration has the form: c c void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *, c int *, int *, int *, 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_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom, c &eng_vdwl,eatom,ntype,type,fmap, c &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0], c &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1, c dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag); c subroutine meam_dens_final(nlocal, nmax, $ eflag_either, eflag_global, eflag_atom, eng_vdwl, eatom, $ ntype, type, fmap, $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, $ Gamma, dGamma1, dGamma2, dGamma3, $ rho, rho0, rho1, rho2, rho3, fp, errorflag) use meam_data implicit none integer nlocal, nmax, eflag_either, eflag_global, eflag_atom integer ntype, type, fmap real*8 eng_vdwl, eatom, Arho1, Arho2 real*8 Arho2b, Arho3, Arho3b real*8 t_ave, tsq_ave real*8 Gamma, dGamma1, dGamma2, dGamma3 real*8 rho, rho0, rho1, rho2, rho3 real*8 fp integer errorflag dimension eatom(nmax) 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 tsq_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, rho_bkgd c Complete the calculation of density do i = 1,nlocal elti = fmap(type(i)) if (elti.gt.0) then 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 if (ialloy.eq.1) then t_ave(1,i) = t_ave(1,i)/tsq_ave(1,i) t_ave(2,i) = t_ave(2,i)/tsq_ave(2,i) t_ave(3,i) = t_ave(3,i)/tsq_ave(3,i) else if (ialloy.eq.2) then t_ave(1,i) = t1_meam(elti) t_ave(2,i) = t2_meam(elti) t_ave(3,i) = t3_meam(elti) else 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 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 Z = Z_meam(elti) 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)) if (mix_ref_t.eq.1) then gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) $ +t_ave(3,i)*shpi(3))/(Z*Z) else gam = (t1_meam(elti)*shp(1)+t2_meam(elti)*shp(2) $ +t3_meam(elti)*shp(3))/(Z*Z) endif call G_gam(gam,ibar_meam(elti),gsmooth_factor, $ Gbar,errorflag) endif rho(i) = rho0(i) * G if (mix_ref_t.eq.1) then 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 rho_bkgd = rho0_meam(elti)*Z*Gbar else - rho_bkgd = rho_ref_meam(elti) + if (bkgd_dyn.eq.1) then + rho_bkgd = rho0_meam(elti)*Z + else + rho_bkgd = rho_ref_meam(elti) + endif endif rhob = rho(i)/rho_bkgd denom = 1.d0/rho_bkgd call dG_gam(Gamma(i),ibar_meam(elti),gsmooth_factor,G,dG) 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 c dGamma3 is nonzero only if we are using the "mixed" rule for c computing t in the reference system (which is not correct, but c included for backward compatibility if (mix_ref_t.eq.1) then dGamma3(i) = rho0(i)*G*dGbar/(Gbar*Z*Z)*denom else dGamma3(i) = 0.0 endif B = A_meam(elti)*Ec_meam(elti,elti) if( rhob .ne. 0.d0 ) then - fp(i) = B*(log(rhob)+1.d0) + if (emb_lin_neg.eq.1 .and. rhob.le.0) then + fp(i) = -B + else + fp(i) = B*(log(rhob)+1.d0) + endif if (eflag_either.ne.0) then if (eflag_global.ne.0) then - eng_vdwl = eng_vdwl + B*rhob*log(rhob) + if (emb_lin_neg.eq.1 .and. rhob.le.0) then + eng_vdwl = eng_vdwl - B*rhob + else + eng_vdwl = eng_vdwl + B*rhob*log(rhob) + endif endif if (eflag_atom.ne.0) then - eatom(i) = eatom(i) + B*rhob*log(rhob) + if (emb_lin_neg.eq.1 .and. rhob.le.0) then + eatom(i) = eatom(i) - B*rhob + else + eatom(i) = eatom(i) + B*rhob*log(rhob) + endif endif endif else - fp(i) = B + if (emb_lin_neg.eq.1) then + fp(i) = -B + else + fp(i) = B + endif endif 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) +c -5 => G = +-sqrt(abs(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 if (ibar.eq.-5) then + if ((1.d0+Gamma).ge.0) then + G = sqrt(1.d0+Gamma) + else + G = -sqrt(-1.d0-Gamma) + endif 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) +c -5 => G = +-sqrt(abs(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 + else if (ibar.eq.-5) then + if ((1.d0+Gamma).ge.0) then + G = sqrt(1.d0+Gamma) + dG = 1.d0/(2.d0*G) + else + G = -sqrt(-1.d0-Gamma) + dG = -1.d0/(2.d0*G) + endif endif return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc diff --git a/lib/meam/meam_setup_done.F b/lib/meam/meam_setup_done.F index 597d9a6f2..b9f52bc06 100755 --- a/lib/meam/meam_setup_done.F +++ b/lib/meam/meam_setup_done.F @@ -1,1005 +1,1019 @@ 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 background densities for reference structure call compute_reference_density 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,rarat,scrn,scrn2 real*8 phiaa,phibb,phitmp real*8 C,s111,s112,s221,S11,S22 real*8, external :: phi_meam real*8, external :: zbl real*8, external :: compute_phi 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 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)) c The B1, B2, and L12 cases with NN2 have a trick to them; we need to c compute the contributions from second nearest neighbors, like a-a c pairs, but need to include NN2 contributions to those pairs as c well. if (lattce_meam(a,b).eq.'b1'.or. $ lattce_meam(a,b).eq.'b2'.or. $ lattce_meam(a,b).eq.'l12') then rarat = r*arat c phi_aa phiaa = phi_meam(rarat,a,a) call get_Zij(Z1,lattce_meam(a,a)) call get_Zij2(Z2,arat,scrn,lattce_meam(a,a), $ Cmin_meam(a,a,a),Cmax_meam(a,a,a)) nmax = 10 if (scrn.gt.0.0) then do n = 1,nmax phiaa = phiaa + $ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,a,a) enddo endif c phi_bb phibb = phi_meam(rarat,b,b) call get_Zij(Z1,lattce_meam(b,b)) call get_Zij2(Z2,arat,scrn,lattce_meam(b,b), $ Cmin_meam(b,b,b),Cmax_meam(b,b,b)) nmax = 10 if (scrn.gt.0.0) then do n = 1,nmax phibb = phibb + $ (-Z2*scrn/Z1)**n * phi_meam(rarat*arat**n,b,b) enddo endif if (lattce_meam(a,b).eq.'b1'. $ or.lattce_meam(a,b).eq.'b2') then c Add contributions to the B1 or B2 potential 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)) phir(j,nv2) = phir(j,nv2) - $ Z2*scrn/(2*Z1) * phiaa call get_Zij2(Z2,arat,scrn2,lattce_meam(a,b), $ Cmin_meam(b,b,a),Cmax_meam(b,b,a)) phir(j,nv2) = phir(j,nv2) - $ Z2*scrn2/(2*Z1) * phibb else if (lattce_meam(a,b).eq.'l12') then c The L12 case has one last trick; we have to be careful to compute c the correct screening between 2nd-neighbor pairs. 1-1 c second-neighbor pairs are screened by 2 type 1 atoms and two type c 2 atoms. 2-2 second-neighbor pairs are screened by 4 type 1 c atoms. C = 1.d0 call get_sijk(C,a,a,a,s111) call get_sijk(C,a,a,b,s112) call get_sijk(C,b,b,a,s221) S11 = s111 * s111 * s112 * s112 S22 = s221**4 phir(j,nv2) = phir(j,nv2) - $ 0.75*S11*phiaa - 0.25*S22*phibb endif else 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 endif c For Zbl potential: 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 if (zbl_meam(a,b).eq.1) then 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 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 real*8 arat,scrn,scrn2 integer Z12, errorflag integer n,nmax,Z1nn,Z2nn character*3 latta,lattb real*8 rho_bkgd1, rho_bkgd2 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 if densities are too small, numerical problems may result; just return zero if (rho01.le.1e-14.and.rho02.le.1e-14) then phi_m = 0.0 return endif c calculate average weighting factors for the reference structure if (lattce_meam(a,b).eq.'c11') then if (ialloy.eq.2) then t11av = t1_meam(a) t12av = t1_meam(b) t21av = t2_meam(a) t22av = t2_meam(b) t31av = t3_meam(a) t32av = t3_meam(b) else 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 endif 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), $ r,a,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 If using mixing rule for t, apply to reference structure; else c use precomputed values if (mix_ref_t.eq.1) then 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,errorflag) 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,errorflag) endif rho0_1 = rho0_meam(a)*Z1*G1 rho0_2 = rho0_meam(b)*Z2*G2 endif 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,errorflag) call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag) if (mix_ref_t.eq.1) then rho_bkgd1 = rho0_1 rho_bkgd2 = rho0_2 else - rho_bkgd1 = rho_ref_meam(a) - rho_bkgd2 = rho_ref_meam(b) + if (bkgd_dyn.eq.1) then + rho_bkgd1 = rho0_meam(a)*Z_meam(a) + rho_bkgd2 = rho0_meam(b)*Z_meam(b) + else + rho_bkgd1 = rho_ref_meam(a) + rho_bkgd2 = rho_ref_meam(b) + endif endif rhobar1 = rho01/rho_bkgd1*G1 rhobar2 = rho02/rho_bkgd2*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) + if (emb_lin_neg.eq.1 .and. rhobar1.le.0) then + F1 = -A_meam(a)*Ec_meam(a,a)*rhobar1 + else + F1 = A_meam(a)*Ec_meam(a,a)*rhobar1*log(rhobar1) + endif endif if (rhobar2.eq.0.d0) then F2 = 0.d0 else - F2 = A_meam(b)*Ec_meam(b,b)*rhobar2*log(rhobar2) + if (emb_lin_neg.eq.1 .and. rhobar2.le.0) then + F2 = -A_meam(b)*Ec_meam(b,b)*rhobar2 + else + F2 = A_meam(b)*Ec_meam(b,b)*rhobar2*log(rhobar2) + endif 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),erose_form) 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 else if (lattce_meam(a,b).eq.'l12') then phiaa = phi_meam(r,a,a) c account for second neighbor a-a potential here... call get_Zij(Z1nn,lattce_meam(a,a)) call get_Zij2(Z2nn,arat,scrn,lattce_meam(a,a), $ Cmin_meam(a,a,a),Cmax_meam(a,a,a)) nmax = 10 if (scrn.gt.0.0) then do n = 1,nmax phiaa = phiaa + $ (-Z2nn*scrn/Z1nn)**n * phi_meam(r*arat**n,a,a) enddo endif phi_m = Eu/3. - F1/4. - F2/12. - phiaa 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 Compute background density for reference structure of each element subroutine compute_reference_density use meam_data implicit none integer a,Z,Z2,errorflag real*8 gam,Gbar,shp(3) real*8 rho0,rho0_2nn,arat,scrn c loop over element types do a = 1,neltypes Z = Z_meam(a) if (ibar_meam(a).le.0) then Gbar = 1.d0 else call get_shpfcn(shp,lattce_meam(a,a)) gam = (t1_meam(a)*shp(1)+t2_meam(a)*shp(2) $ +t3_meam(a)*shp(3))/(Z*Z) call G_gam(gam,ibar_meam(a),gsmooth_factor, $ Gbar,errorflag) endif c The zeroth order density in the reference structure, with c equilibrium spacing, is just the number of first neighbors times c the rho0_meam coefficient... rho0 = rho0_meam(a)*Z c ...unless we have unscreened second neighbors, in which case we c add on the contribution from those (accounting for partial c screening) if (nn2_meam(a,a).eq.1) then call get_Zij2(Z2,arat,scrn,lattce_meam(a,a), $ Cmin_meam(a,a,a),Cmax_meam(a,a,a)) rho0_2nn = rho0_meam(a)*exp(-beta0_meam(a)*(arat-1)) rho0 = rho0 + Z2*rho0_2nn*scrn endif rho_ref_meam(a) = rho0*Gbar enddo 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'.or.latt.eq.'b2') 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, $ r,a,b,latt) use meam_data implicit none real*8 t11av,t21av,t31av,t12av,t22av,t32av real*8 t11,t21,t31,t12,t22,t32,r integer a,b character*3 latt real*8 rhoa01,rhoa02,a1,a2,rho01,rho02 c For ialloy = 2, no averaging is done if (ialloy.eq.2) then t11av = t11 t21av = t21 t31av = t31 t12av = t12 t22av = t22 t32av = t32 else if (latt.eq.'fcc'.or.latt.eq.'bcc'.or.latt.eq.'dia' $ .or.latt.eq.'hcp'.or.latt.eq.'b1' $ .or.latt.eq.'dim'.or.latt.eq.'b2') then c all neighbors are of the opposite type t11av = t12 t21av = t22 t31av = t32 t12av = t11 t22av = t21 t32av = t31 else 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) rhoa02 = rho0_meam(b)*exp(-beta0_meam(b)*a2) if (latt.eq.'l12') then rho01 = 8*rhoa01 + 4*rhoa02 t11av = (8*t11*rhoa01 + 4*t12*rhoa02)/rho01 t12av = t11 t21av = (8*t21*rhoa01 + 4*t22*rhoa02)/rho01 t22av = t21 t31av = (8*t31*rhoa01 + 4*t32*rhoa02)/rho01 t32av = t31 else c call error('Lattice not defined in get_tavref.') endif endif 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 if (latt.eq.'l12') then Zij = 12 else if (latt.eq.'b2') then Zij = 8 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 = 12 a = sqrt(2.d0) numscr = 2 else if (latt.eq.'l12') then Zij2 = 6 a = sqrt(2.d0) numscr = 4 else if (latt.eq.'b2') then Zij2 = 6 a = 2.d0/sqrt(3.d0) numscr = 4 else if (latt.eq.'dim') then c this really shouldn't be allowed; make sure screening is zero Zij2 = 0 a = 1 S = 0 return 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 subroutine get_sijk(C,i,j,k,sijk) use meam_data implicit none real*8 C,sijk integer i,j,k real*8 x x = (C-Cmin_meam(i,j,k))/(Cmax_meam(i,j,k)-Cmin_meam(i,j,k)) call fcut(x,sijk) 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 s(3) character*3 lat integer a,b integer Zij1nn,Zij2nn real*8 rhoa01nn,rhoa02nn real*8 arat,scrn,denom real*8 C,s111,s112,s221,S11,S22 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 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 if (lat.eq.'l12') then rho01 = 8*rhoa01 + 4*rhoa02 rho02 = 12*rhoa01 if (ialloy.eq.1) then rho21 = 8./3.*(rhoa21*t2_meam(a)-rhoa22*t2_meam(b))**2 denom = 8*rhoa01*t2_meam(a)**2 + 4*rhoa02*t2_meam(b)**2 if (denom.gt.0.) then rho21 = rho21/denom * rho01 endif else rho21 = 8./3.*(rhoa21-rhoa22)*(rhoa21-rhoa22) endif else if (lat.eq.'b2') then rho01 = 8.d0*rhoa02 rho02 = 8.d0*rhoa01 else c call error('Lattice not defined in get_densref.') endif if (nn2_meam(a,b).eq.1) then 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) if (lat.eq.'l12') then c As usual, L12 thinks it's special; we need to be careful computing c the screening functions C = 1.d0 call get_sijk(C,a,a,a,s111) call get_sijk(C,a,a,b,s112) call get_sijk(C,b,b,a,s221) S11 = s111 * s111 * s112 * s112 S22 = s221**4 rho01 = rho01 + 6*S11*rhoa01nn rho02 = rho02 + 6*S22*rhoa02nn else c For other cases, assume that second neighbor is of same type, c first neighbor may be of different type rho01 = rho01 + Zij2nn*scrn*rhoa01nn c Assume Zij2nn and arat don't depend on order, but scrn might call get_Zij2(Zij2nn,arat,scrn,lat, $ Cmin_meam(b,b,a),Cmax_meam(b,b,a)) rho02 = rho02 + Zij2nn*scrn*rhoa02nn endif 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,form) implicit none real*8 r,re,alpha,Ec,repuls,attrac,astar,a3 integer form 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 if (form.eq.1) then erose = -Ec*(1+astar+(-attrac+repuls/r)* $ (astar**3))*exp(-astar) else if (form.eq.2) then erose = -Ec * (1 +astar + a3*(astar**3))*exp(-astar) else erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar) endif 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 c--------------------------------------------------------------------- c Compute Rose energy function, I.16 c real*8 function compute_phi(rij, elti, eltj) use meam_data implicit none real*8 rij, pp integer elti, eltj, ind, kk ind = eltind(elti, eltj) pp = rij*rdrar + 1.0D0 kk = pp kk = min(kk,nrar-1) pp = pp - kk pp = min(pp,1.0D0) compute_phi = ((phirar3(kk,ind)*pp + phirar2(kk,ind))*pp $ + phirar1(kk,ind))*pp + phirar(kk,ind) return end diff --git a/lib/meam/meam_setup_global.F b/lib/meam/meam_setup_global.F index 548c97b25..97c4a8427 100755 --- a/lib/meam/meam_setup_global.F +++ b/lib/meam/meam_setup_global.F @@ -1,109 +1,111 @@ 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 zbl_meam(:,:) = 1 gsmooth_factor = 99.0 augt1 = 1 ialloy = 0 mix_ref_t = 0 + emb_lin_neg = 0 + bkgd_dyn = 0 erose_form = 0 return end diff --git a/lib/meam/meam_setup_param.F b/lib/meam/meam_setup_param.F index 70fccf8da..922eb2c57 100755 --- a/lib/meam/meam_setup_param.F +++ b/lib/meam/meam_setup_param.F @@ -1,150 +1,160 @@ 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 c 15 = ialloy c 16 = mixture_ref_t c 17 = erose_form c 18 = zbl_meam +c 19 = emb_lin_neg +c 20 = bkgd_dyn subroutine meam_setup_param(which, value, nindex, $ index, errorflag) use meam_data implicit none integer which, nindex, index(3), errorflag real*8 value integer i1, i2 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' else if (value.eq.7) then lattce_meam(index(1),index(2)) = 'l12' else if (value.eq.8) then lattce_meam(index(1),index(2)) = 'b2' 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 i1 = min(index(1),index(2)) i2 = max(index(1),index(2)) nn2_meam(i1,i2) = 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 c 15 = ialloy else if (which.eq.15) then ialloy = value c 16 = mixture_ref_t else if (which.eq.16) then mix_ref_t = value c 17 = erose_form else if (which.eq.17) then erose_form = value c 18 = zbl_meam else if (which.eq.18) then i1 = min(index(1),index(2)) i2 = max(index(1),index(2)) zbl_meam(i1,i2) = value +c 19 = emb_lin_neg + else if (which.eq.19) then + emb_lin_neg = value + +c 20 = bkgd_dyn + else if (which.eq.20) then + bkgd_dyn = value + else errorflag = 1 endif return end