From 8a6024b3ea937beb95a4db9d54aa7db237d08511 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 18 Aug 2009 17:02:09 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3084 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/meam/meam_data.F | 3 +- lib/meam/meam_dens_final.F | 21 +++-- lib/meam/meam_dens_init.F | 40 +++++++-- lib/meam/meam_force.F | 160 ++++++++++++++++++++++++++------- lib/meam/meam_setup_done.F | 167 +++++++++++++++++++++++++---------- lib/meam/meam_setup_global.F | 1 + lib/meam/meam_setup_param.F | 7 ++ 7 files changed, 300 insertions(+), 99 deletions(-) diff --git a/lib/meam/meam_data.F b/lib/meam/meam_data.F index 9269c7ccb7..8562e9f03b 100755 --- a/lib/meam/meam_data.F +++ b/lib/meam/meam_data.F @@ -34,6 +34,7 @@ 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 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 @@ -65,7 +66,7 @@ c nrar,rdrar = spline coeff array parameters 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 + integer augt1, ialloy real*8 gsmooth_factor integer vind2D(3,3),vind3D(3,3,3) integer v2D(6),v3D(10) diff --git a/lib/meam/meam_dens_final.F b/lib/meam/meam_dens_final.F index 04f530ce9c..d39de8ba66 100755 --- a/lib/meam/meam_dens_final.F +++ b/lib/meam/meam_dens_final.F @@ -4,21 +4,21 @@ 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 *, int *); +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],gamma,dgamma1, +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, + $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, $ Gamma, dGamma1, dGamma2, dGamma3, $ rho, rho0, rho1, rho2, rho3, fp, errorflag) @@ -29,7 +29,7 @@ c integer ntype, type, fmap real*8 eng_vdwl, eatom, Arho1, Arho2 real*8 Arho2b, Arho3, Arho3b - real*8 t_ave + real*8 t_ave, tsq_ave real*8 Gamma, dGamma1, dGamma2, dGamma3 real*8 rho, rho0, rho1, rho2, rho3 real*8 fp @@ -39,6 +39,7 @@ c 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) @@ -70,9 +71,15 @@ c Complete the calculation of density 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) + 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 + 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) diff --git a/lib/meam/meam_dens_init.F b/lib/meam/meam_dens_init.F index df62ecabd4..2a4b2be68a 100755 --- a/lib/meam/meam_dens_init.F +++ b/lib/meam/meam_dens_init.F @@ -3,7 +3,7 @@ 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 double *, double *, double *, double *, double *, int *); c c c Call from pair_meam.cpp has the form: @@ -12,7 +12,7 @@ c meam_dens_init_(&i,&nmax,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 &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],&errorflag); c subroutine meam_dens_init(i, nmax, @@ -20,7 +20,7 @@ c $ numneigh, firstneigh, $ numneigh_full, firstneigh_full, $ scrfcn, dscrfcn, fcpair, rho0, arho1, arho2, arho2b, - $ arho3, arho3b, t_ave, errorflag) + $ arho3, arho3b, t_ave, tsq_ave, errorflag) use meam_data implicit none @@ -30,7 +30,7 @@ c integer numneigh, firstneigh, numneigh_full, firstneigh_full real*8 scrfcn, dscrfcn, fcpair real*8 rho0, arho1, arho2 - real*8 arho2b, arho3, arho3b, t_ave + real*8 arho2b, arho3, arho3b, t_ave, tsq_ave integer errorflag integer j,jn @@ -40,7 +40,7 @@ c 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) + dimension t_ave(3,nmax), tsq_ave(3,nmax) errorflag = 0 @@ -54,7 +54,7 @@ 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) + $ arho3, arho3b, t_ave, tsq_ave) return end @@ -194,7 +194,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine calc_rho1(i, nmax, ntype, type, fmap, x, $ numneigh, firstneigh, $ scrfcn, fcpair, rho0, arho1, arho2, arho2b, - $ arho3, arho3b, t_ave) + $ arho3, arho3b, t_ave, tsq_ave) use meam_data implicit none @@ -203,14 +203,14 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc real*8 x integer numneigh, firstneigh real*8 scrfcn, fcpair, rho0, arho1, arho2 - real*8 arho2b, arho3, arho3b, t_ave + real*8 arho2b, arho3, arho3b, t_ave, tsq_ave dimension type(nmax), fmap(ntype), x(3,nmax) dimension firstneigh(numneigh) dimension scrfcn(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) + dimension t_ave(3,nmax), tsq_ave(3,nmax) integer jn,j,m,n,p,elti,eltj integer nv2,nv3 @@ -248,6 +248,14 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc rhoa1i = ro0i*exp(-beta1_meam(elti)*ai)*sij rhoa2i = ro0i*exp(-beta2_meam(elti)*ai)*sij rhoa3i = ro0i*exp(-beta3_meam(elti)*ai)*sij + if (ialloy.eq.1) then + rhoa1j = rhoa1j * t1_meam(eltj) + rhoa2j = rhoa2j * t2_meam(eltj) + rhoa3j = rhoa3j * t3_meam(eltj) + rhoa1i = rhoa1i * t1_meam(elti) + rhoa2i = rhoa2i * t2_meam(elti) + rhoa3i = rhoa3i * t3_meam(elti) + endif rho0(i) = rho0(i) + rhoa0j rho0(j) = rho0(j) + rhoa0i t_ave(1,i) = t_ave(1,i) + t1_meam(eltj)*rhoa0j @@ -256,6 +264,20 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc t_ave(1,j) = t_ave(1,j) + t1_meam(elti)*rhoa0i t_ave(2,j) = t_ave(2,j) + t2_meam(elti)*rhoa0i t_ave(3,j) = t_ave(3,j) + t3_meam(elti)*rhoa0i + if (ialloy.eq.1) then + tsq_ave(1,i) = tsq_ave(1,i) + + $ t1_meam(eltj)*t1_meam(eltj)*rhoa0j + tsq_ave(2,i) = tsq_ave(2,i) + + $ t2_meam(eltj)*t2_meam(eltj)*rhoa0j + tsq_ave(3,i) = tsq_ave(3,i) + + $ t3_meam(eltj)*t3_meam(eltj)*rhoa0j + tsq_ave(1,j) = tsq_ave(1,j) + + $ t1_meam(elti)*t1_meam(elti)*rhoa0i + tsq_ave(2,j) = tsq_ave(2,j) + + $ t2_meam(elti)*t2_meam(elti)*rhoa0i + tsq_ave(3,j) = tsq_ave(3,j) + + $ t3_meam(elti)*t3_meam(elti)*rhoa0i + endif Arho2b(i) = Arho2b(i) + rhoa2j Arho2b(j) = Arho2b(j) + rhoa2i diff --git a/lib/meam/meam_force.F b/lib/meam/meam_force.F index 6b22d81c3d..05329a4b41 100755 --- a/lib/meam/meam_force.F +++ b/lib/meam/meam_force.F @@ -4,7 +4,7 @@ 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 double *, double *, double *, double *, double *, double *, int *); c c Call from pair_meam.cpp has the form: c @@ -14,7 +14,7 @@ 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],&vatom[0][0],&errorflag); +c &t_ave[0][0],&tsq_ave[0][0],&f[0][0],&vatom[0][0],&errorflag); c subroutine meam_force(i, nmax, @@ -23,8 +23,8 @@ c $ 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, vatom, - $ errorflag) + $ Arho1, Arho2, Arho2b, Arho3, Arho3b, t_ave, tsq_ave, f, + $ vatom, errorflag) use meam_data implicit none @@ -38,7 +38,7 @@ c real*8 rho0, rho1, rho2, rho3, fp real*8 Arho1, Arho2, Arho2b real*8 Arho3, Arho3b - real*8 t_ave, f, vatom + real*8 t_ave, tsq_ave, f, vatom integer errorflag dimension eatom(nmax) @@ -50,7 +50,7 @@ c 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), vatom(6,nmax) + dimension t_ave(3,nmax), tsq_ave(3,nmax), f(3,nmax), vatom(6,nmax) integer i,j,jn,k,kn,kk,m,n,p,q integer nv2,nv3,elti,eltj,eltk,ind @@ -176,6 +176,21 @@ c Compute pair densities and derivatives drhoa3j = drhoa3i endif + if (ialloy.eq.1) then + rhoa1j = rhoa1j * t1_meam(eltj) + rhoa2j = rhoa2j * t2_meam(eltj) + rhoa3j = rhoa3j * t3_meam(eltj) + rhoa1i = rhoa1i * t1_meam(elti) + rhoa2i = rhoa2i * t2_meam(elti) + rhoa3i = rhoa3i * t3_meam(elti) + drhoa1j = drhoa1j * t1_meam(eltj) + drhoa2j = drhoa2j * t2_meam(eltj) + drhoa3j = drhoa3j * t3_meam(eltj) + drhoa1i = drhoa1i * t1_meam(elti) + drhoa2i = drhoa2i * t2_meam(elti) + drhoa3i = drhoa3i * t3_meam(elti) + endif + nv2 = 1 nv3 = 1 arg1i1 = 0.d0 @@ -277,21 +292,59 @@ c Compute derivatives of weighting functions t wrt rij 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 + if (ialloy.eq.1) then - 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) + a1i = 0.d0 + a1j = 0.d0 + a2i = 0.d0 + a2j = 0.d0 + a3i = 0.d0 + a3j = 0.d0 + if ( tsq_ave(1,i) .ne. 0.d0 ) then + a1i = drhoa0j*sij/tsq_ave(1,i) + endif + if ( tsq_ave(1,j) .ne. 0.d0 ) then + a1j = drhoa0i*sij/tsq_ave(1,j) + endif + if ( tsq_ave(2,i) .ne. 0.d0 ) then + a2i = drhoa0j*sij/tsq_ave(2,i) + endif + if ( tsq_ave(2,j) .ne. 0.d0 ) then + a2j = drhoa0i*sij/tsq_ave(2,j) + endif + if ( tsq_ave(3,i) .ne. 0.d0 ) then + a3i = drhoa0j*sij/tsq_ave(3,i) + endif + if ( tsq_ave(3,j) .ne. 0.d0 ) then + a3j = drhoa0i*sij/tsq_ave(3,j) + endif + + dt1dr1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2) + dt1dr2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2) + dt2dr1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2) + dt2dr2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2) + dt3dr1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2) + dt3dr2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2) + + else + + 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) + + endif c Compute derivatives of total density wrt rij, sij and rij(3) call get_shpfcn(shpi,lattce_meam(elti,elti)) @@ -340,21 +393,60 @@ c Compute derivatives wrt sij, but only if necessary 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 + if (ialloy.eq.1) then + + a1i = 0.d0 + a1j = 0.d0 + a2i = 0.d0 + a2j = 0.d0 + a3i = 0.d0 + a3j = 0.d0 + if ( tsq_ave(1,i) .ne. 0.d0 ) then + a1i = rhoa0j/tsq_ave(1,i) + endif + if ( tsq_ave(1,j) .ne. 0.d0 ) then + a1j = rhoa0i/tsq_ave(1,j) + endif + if ( tsq_ave(2,i) .ne. 0.d0 ) then + a2i = rhoa0j/tsq_ave(2,i) + endif + if ( tsq_ave(2,j) .ne. 0.d0 ) then + a2j = rhoa0i/tsq_ave(2,j) + endif + if ( tsq_ave(3,i) .ne. 0.d0 ) then + a3i = rhoa0j/tsq_ave(3,i) + endif + if ( tsq_ave(3,j) .ne. 0.d0 ) then + a3j = rhoa0i/tsq_ave(3,j) + endif + + dt1ds1 = a1i*(t1_meam(eltj)-t1i*t1_meam(eltj)**2) + dt1ds2 = a1j*(t1_meam(elti)-t1j*t1_meam(elti)**2) + dt2ds1 = a2i*(t2_meam(eltj)-t2i*t2_meam(eltj)**2) + dt2ds2 = a2j*(t2_meam(elti)-t2j*t2_meam(elti)**2) + dt3ds1 = a3i*(t3_meam(eltj)-t3i*t3_meam(eltj)**2) + dt3ds2 = a3j*(t3_meam(elti)-t3j*t3_meam(elti)**2) + + else + + 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) + + endif - 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 diff --git a/lib/meam/meam_setup_done.F b/lib/meam/meam_setup_done.F index 95160cd24f..3bf5962945 100755 --- a/lib/meam/meam_setup_done.F +++ b/lib/meam/meam_setup_done.F @@ -158,11 +158,12 @@ c integer j,a,b,nv2 real*8 astar,frac,phizbl integer n,nmax,Z1,Z2 - real*8 arat,scrn + real*8 arat,scrn,scrn2 real*8 phitmp real*8, external :: phi_meam real*8, external :: zbl + real*8, external :: compute_phi c allocate memory for array that defines the potential @@ -178,6 +179,10 @@ c allocate coeff memory allocate(phirar5(nr,(neltypes*(neltypes+1))/2)) allocate(phirar6(nr,(neltypes*(neltypes+1))/2)) +c HACK for debug: compute phi_meam for Ni3Si equilibrium spacing + r = 2.47770216 + phitmp = phi_meam(r,1,2) + c loop over pairs of element types nv2 = 0 do a = 1,neltypes @@ -194,20 +199,29 @@ c loop over r values and compute 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 if (a.ne.b) then +c write(6,*) 'Second-neighbor currently only valid ' +c write(6,*) 'if element types are the same.' c call error(' ') - endif +c 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 + if (lattce_meam(a,b).eq.'b1') then + phir(j,nv2) = phir(j,nv2) - + $ Z2*scrn/(2*Z1) * phi_meam(r*arat,a,a) + 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) * phi_meam(r*arat,b,b) + 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 @@ -258,7 +272,7 @@ c real*8 rho02,rho12,rho22,rho32 real*8 scalfac,phiaa,phibb real*8 Eu - integer Z12 + integer Z12, errorflag character*3 latta,lattb real*8, external :: erose @@ -285,8 +299,9 @@ c calculate average weighting factors for the reference structure 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)) + $ 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 @@ -316,6 +331,12 @@ 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. +c -- GJW (6/12/09): I suspect there may be a problem here... since we should be +c using the pure element reference structure, we should also +c be using the element "t" values (not the averages compute +c above). It doesn't matter for reference structures with +c s(1)=s(2)=s(3)=0, but might cause diffs otherwise. For now +c what's here is consistent with what's done in meam_dens_final. Z1 = Z_meam(a) Z2 = Z_meam(b) if (ibar_meam(a).le.0) then @@ -323,14 +344,14 @@ c in the reference structure. 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) + 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) + 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 @@ -347,8 +368,8 @@ c compute total background density, eqn I.6 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) + call G_gam(Gam1,ibar_meam(a),gsmooth_factor,G1,errorflag) + call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag) rhobar1 = rho01/rho0_1*G1 rhobar2 = rho02/rho0_2*G2 endif @@ -370,17 +391,20 @@ c compute Rose function, I.16 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 + 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 if (lattce_meam(a,b).eq.'l12') then + phiaa = phi_meam(r,a,a) + phi_m = Eu/3. - F1/4. - F2/12. - phiaa else -c +c c potential is computed from Rose function and embedding energy phi_m = (2*Eu - F1 - F2)/Z12 c @@ -426,11 +450,15 @@ c call error('Lattice not defined in get_shpfcn.') 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) + $ 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 + real*8 t11,t21,t31,t12,t22,t32,r + integer a,b character*3 latt + real*8 rhoa01,rhoa02,a1,a2,rho01,rho02 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 @@ -442,7 +470,21 @@ c all neighbors are of the opposite type t22av = t21 t32av = t31 else -c call error('Lattice not defined in get_tavref.') + 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 return end @@ -466,6 +508,8 @@ c Number of neighbors for the reference structure Zij = 1 else if (latt.eq.'c11') then Zij = 10 + else if (latt.eq.'l12') then + Zij = 12 else c call error('Lattice not defined in get_Zij.') endif @@ -504,9 +548,9 @@ c call error('can not do 2NN MEAM for dia') a = sqrt(2.d0) numscr = 4 else if (latt.eq.'b1') then - Zij2 = 6 + Zij2 = 12 a = sqrt(2.d0) - numscr = 4 + numscr = 2 else c call error('Lattice not defined in get_Zij2.') endif @@ -535,7 +579,7 @@ c Calculate density functions, assuming reference configuration integer a,b integer Zij1nn,Zij2nn real*8 rhoa01nn,rhoa02nn - real*8 arat,scrn + real*8 arat,scrn,denom a1 = r/re_meam(a,a) - 1.d0 a2 = r/re_meam(b,b) - 1.d0 @@ -567,17 +611,6 @@ c Calculate density functions, assuming reference configuration 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 @@ -609,6 +642,18 @@ c rho02 = 6*rhoa01 rho22 = rhoa22 rho31 = rhoa31 rho32 = rhoa32 + else if (lat.eq.'l12') then + rho01 = 8*rhoa01 + 4*rhoa02 + rho02 = 12*rhoa01 + if (ialloy.eq.0) then + rho21 = 8./3.*(rhoa21-rhoa22)*(rhoa21-rhoa22) + else + 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 + endif else c call error('Lattice not defined in get_densref.') endif @@ -617,9 +662,9 @@ c call error('Lattice not defined in get_densref.') c Assume that second neighbor is of same type, first neighbor c may be of different type - if (lat.ne.'bcc') then +c if (lat.ne.'bcc') then c call error('Second-neighbor not defined for this lattice') - endif +c endif call get_Zij2(Zij2nn,arat,scrn,lat, $ Cmin_meam(a,a,b),Cmax_meam(a,a,b)) @@ -630,8 +675,12 @@ c call error('Second-neighbor not defined for this lattice') 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 + 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 @@ -731,3 +780,25 @@ c phir interp 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 1223e34b22..9603c4ecfc 100755 --- a/lib/meam/meam_setup_global.F +++ b/lib/meam/meam_setup_global.F @@ -98,6 +98,7 @@ c Set some defaults nn2_meam(:,:) = 0 gsmooth_factor = 99.0 augt1 = 1 + ialloy = 0 return end diff --git a/lib/meam/meam_setup_param.F b/lib/meam/meam_setup_param.F index 7a26ba71dd..d70afa74c8 100755 --- a/lib/meam/meam_setup_param.F +++ b/lib/meam/meam_setup_param.F @@ -27,6 +27,7 @@ c 11 = delr_meam c 12 = augt1 c 13 = gsmooth_factor c 14 = re_meam +c 15 = ialloy subroutine meam_setup_param(which, value, nindex, $ index, errorflag) @@ -71,6 +72,8 @@ c 4 = lattce_meam 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' endif c 5 = attrac_meam @@ -113,6 +116,10 @@ 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 + else errorflag = 1 endif