diff --git a/lib/meam/meam_setup_done.F b/lib/meam/meam_setup_done.F
index 6225a1627..7907cdf96 100755
--- a/lib/meam/meam_setup_done.F
+++ b/lib/meam/meam_setup_done.F
@@ -1,813 +1,854 @@
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 i<j case (which has aready been set,
c here or in the input file)
if (i.gt.j) then
re_meam(i,j) = re_meam(j,i)
Ec_meam(i,j) = Ec_meam(j,i)
alpha_meam(i,j) = alpha_meam(j,i)
lattce_meam(i,j) = lattce_meam(j,i)
c If i<j and term is unset, use default values (e.g. mean of i-i and j-j)
else if (j.gt.i) then
if (Ec_meam(i,j).eq.0.d0) then
if (lattce_meam(i,j).eq.'l12') then
Ec_meam(i,j) = (3*Ec_meam(i,i)+Ec_meam(j,j))/4.d0
$ - delta_meam(i,j)
else if (lattce_meam(i,j).eq.'c11') then
if (lattce_meam(i,i).eq.'dia') then
Ec_meam(i,j) = (2*Ec_meam(i,i)+Ec_meam(j,j))/3.d0
$ - delta_meam(i,j)
else
Ec_meam(i,j) = (Ec_meam(i,i)+2*Ec_meam(j,j))/3.d0
$ - delta_meam(i,j)
endif
else
Ec_meam(i,j) = (Ec_meam(i,i)+Ec_meam(j,j))/2.d0
$ - delta_meam(i,j)
endif
endif
if (alpha_meam(i,j).eq.0.d0) then
alpha_meam(i,j) = (alpha_meam(i,i)+alpha_meam(j,j))/2.d0
endif
if (re_meam(i,j).eq.0.d0) then
re_meam(i,j) = (re_meam(i,i)+re_meam(j,j))/2.d0
endif
endif
enddo
enddo
c Cmin(i,k,j) is symmetric in i-j, but not k. For all triplets
c where i>j, set equal to the i<j element. Likewise for Cmax.
do i = 2,neltypes
do j = 1,i-1
do k = 1,neltypes
Cmin_meam(i,j,k) = Cmin_meam(j,i,k)
Cmax_meam(i,j,k) = Cmax_meam(j,i,k)
enddo
enddo
enddo
c ebound gives the squared distance such that, for rik2 or rjk2>ebound,
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,scrn2
- real*8 phitmp
+ real*8 arat,rarat,scrn,scrn2
+ real*8 phiaa,phibb,phitmp
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
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(' ')
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))
+c The B1 case with NN2 has a trick to it; we need to compute the
+c contributions from second nearest neighbors, like a-a pairs, but
+c need to include NN2 contributions to those pairs as well.
if (lattce_meam(a,b).eq.'b1') 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 phibb
+ 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
+
+c Now add contributions to the B1 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) * phi_meam(r*arat,a,a)
+ $ 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) * phi_meam(r*arat,b,b)
+ $ Z2*scrn2/(2*Z1) * phibb
+
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 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, errorflag
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 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
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),
$ 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 -- 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
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
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,errorflag)
call G_gam(Gam2,ibar_meam(b),gsmooth_factor,G2,errorflag)
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 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 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,
$ 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
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
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
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
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
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,denom
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.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
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
c if (lat.ne.'bcc') then
c call error('Second-neighbor not defined for this lattice')
c 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*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
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
erose = -Ec * (1 + astar + a3*(astar**3)/(r/re)) * exp(-astar)
c 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
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/src/USER-ACKLAND/compute_ackland_atom.cpp b/src/USER-ACKLAND/compute_ackland_atom.cpp
index a49034e40..20e3cb1ff 100644
--- a/src/USER-ACKLAND/compute_ackland_atom.cpp
+++ b/src/USER-ACKLAND/compute_ackland_atom.cpp
@@ -1,393 +1,393 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: G. Ziegenhain, gerolf@ziegenhain.com
Copyright (C) 2007
------------------------------------------------------------------------- */
#include "string.h"
#include "compute_ackland_atom.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "force.h"
#include "pair.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include <math.h>
using namespace LAMMPS_NS;
enum{UNKNOWN,BCC,FCC,HCP,ICO};
/* ---------------------------------------------------------------------- */
ComputeAcklandAtom::ComputeAcklandAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all("Illegal compute ackland/atom command");
peratom_flag = 1;
size_peratom_cols = 0;
nmax = 0;
structure = NULL;
maxneigh = 0;
distsq = NULL;
nearest = NULL;
nearest_n0 = NULL;
nearest_n1 = NULL;
}
/* ---------------------------------------------------------------------- */
ComputeAcklandAtom::~ComputeAcklandAtom()
{
memory->sfree(structure);
memory->sfree(distsq);
memory->sfree(nearest);
memory->sfree(nearest_n0);
memory->sfree(nearest_n1);
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::init()
{
// need an occasional full neighbor list
int irequest = neighbor->request((void *) this);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->requests[irequest]->occasional = 1;
int count = 0;
for (int i = 0; i < modify->ncompute; i++)
if (strcmp(modify->compute[i]->style,"ackland/atom") == 0) count++;
if (count > 1 && comm->me == 0)
error->warning("More than one compute ackland/atom");
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::init_list(int id, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::compute_peratom()
{
int i,j,ii,jj,k,n,inum,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rsq;
int *ilist,*jlist,*numneigh,**firstneigh;
int chi[8];
invoked_peratom = update->ntimestep;
// grow structure array if necessary
if (atom->nlocal > nmax) {
memory->sfree(structure);
nmax = atom->nmax;
structure = (double *)
memory->smalloc(nmax*sizeof(double),"compute/ackland/atom:ackland");
vector_atom = structure;
}
- // invoke half neighbor list (will copy or build if necessary)
+ // invoke full neighbor list (will copy or build if necessary)
neighbor->build_one(list->index);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// compute structure parameter for each atom in group
// use full neighbor list
double **x = atom->x;
int *mask = atom->mask;
int nall = atom->nlocal + atom->nghost;
double cutsq = force->pair->cutforce * force->pair->cutforce;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = firstneigh[i];
jnum = numneigh[i];
// ensure distsq and nearest arrays are long enough
if (jnum > maxneigh) {
memory->sfree(distsq);
memory->sfree(nearest);
memory->sfree(nearest_n0);
memory->sfree(nearest_n1);
maxneigh = jnum;
distsq = (double *) memory->smalloc(maxneigh*sizeof(double),
"compute/ackland/atom:distsq");
nearest = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest");
nearest_n0 = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest_n0");
nearest_n1 = (int *) memory->smalloc(maxneigh*sizeof(int),
"compute/ackland/atom:nearest_n1");
}
// loop over list of all neighbors within force cutoff
// distsq[] = distance sq to each
// nearest[] = atom indices of neighbors
n = 0;
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
if (j >= nall) j %= nall;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cutsq) {
distsq[n] = rsq;
nearest[n++] = j;
}
}
// Select 6 nearest neighbors
select2(6,n,distsq,nearest);
// Mean squared separation
double r0_sq = 0.;
for (j = 0; j < 6; j++)
r0_sq += distsq[j];
r0_sq /= 6.;
// n0 near neighbors with: distsq<1.45*r0_sq
// n1 near neighbors with: distsq<1.55*r0_sq
double n0_dist_sq = 1.45*r0_sq,
n1_dist_sq = 1.55*r0_sq;
int n0 = 0, n1 = 0;
for (j = 0; j < n; j++) {
if (distsq[j] < n1_dist_sq) {
nearest_n1[n1++] = nearest[j];
if (distsq[j] < n0_dist_sq) {
nearest_n0[n0++] = nearest[j];
}
}
}
// Evaluate all angles <(r_ij,rik) forall n0 particles with: distsq<1.45*r0_sq
double bond_angle;
double norm_j, norm_k;
chi[0] = chi[1] = chi[2] = chi[3] = chi[4] = chi[5] = chi[6] = chi[7] = 0;
double x_ij, y_ij, z_ij, x_ik, y_ik, z_ik;
for (j = 0; j < n0; j++) {
x_ij = x[i][0]-x[nearest_n0[j]][0];
y_ij = x[i][1]-x[nearest_n0[j]][1];
z_ij = x[i][2]-x[nearest_n0[j]][2];
norm_j = sqrt (x_ij*x_ij + y_ij*y_ij + z_ij*z_ij);
if (norm_j <= 0.) continue;
for (k = j+1; k < n0; k++) {
x_ik = x[i][0]-x[nearest_n0[k]][0];
y_ik = x[i][1]-x[nearest_n0[k]][1];
z_ik = x[i][2]-x[nearest_n0[k]][2];
norm_k = sqrt (x_ik*x_ik + y_ik*y_ik + z_ik*z_ik);
if (norm_k <= 0.)
continue;
bond_angle = (x_ij*x_ik + y_ij*y_ik + z_ij*z_ik) / (norm_j*norm_k);
// Histogram for identifying the relevant peaks
if (-1. <= bond_angle && bond_angle < -0.945) { chi[0]++; }
else if (-0.945 <= bond_angle && bond_angle < -0.915) { chi[1]++; }
else if (-0.915 <= bond_angle && bond_angle < -0.755) { chi[2]++; }
else if (-0.755 <= bond_angle && bond_angle < -0.195) { chi[3]++; }
else if (-0.195 <= bond_angle && bond_angle < 0.195) { chi[4]++; }
else if (0.195 <= bond_angle && bond_angle < 0.245) { chi[5]++; }
else if (0.245 <= bond_angle && bond_angle < 0.795) { chi[6]++; }
else if (0.795 <= bond_angle && bond_angle < 1.) { chi[7]++; }
}
}
// Deviations from the different lattice structures
double delta_bcc = 0.35*chi[4]/(double)(chi[5]+chi[6]-chi[4]),
delta_cp = fabs(1.-(double)chi[6]/24.),
delta_fcc = 0.61*(fabs((double)(chi[0]+chi[1]-6.))+(double)chi[2])/6.,
delta_hcp = (fabs((double)chi[0]-3.)+fabs((double)chi[0]+(double)chi[1]+(double)chi[2]+(double)chi[3]-9.))/12.;
// Identification of the local structure according to the reference
if (chi[0] == 7) { delta_bcc = 0.; }
else if (chi[0] == 6) { delta_fcc = 0.; }
else if (chi[0] <= 3) { delta_hcp = 0.; }
if (chi[7] > 0.)
structure[i] = UNKNOWN;
else
if (chi[4] < 3.)
{
if (n1 > 13 || n1 < 11)
structure[i] = UNKNOWN;
else
structure[i] = ICO;
} else
if (delta_bcc <= delta_cp)
{
if (n1 < 11)
structure[i] = UNKNOWN;
else
structure[i] = BCC;
} else
if (n1 > 12 || n1 < 11)
structure[i] = UNKNOWN;
else
if (delta_fcc < delta_hcp)
structure[i] = FCC;
else
structure[i] = HCP;
} else structure[i] = 0.0;
}
}
/* ----------------------------------------------------------------------
2 select routines from Numerical Recipes (slightly modified)
find k smallest values in array of length n
2nd routine sorts auxiliary array at same time
------------------------------------------------------------------------- */
#define SWAP(a,b) tmp = a; a = b; b = tmp;
#define ISWAP(a,b) itmp = a; a = b; b = itmp;
void ComputeAcklandAtom::select(int k, int n, double *arr)
{
int i,ir,j,l,mid;
double a,tmp;
arr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
}
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
}
i = l+1;
j = ir;
a = arr[l+1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
}
arr[l+1] = arr[j];
arr[j] = a;
if (j >= k) ir = j-1;
if (j <= k) l = i;
}
}
}
/* ---------------------------------------------------------------------- */
void ComputeAcklandAtom::select2(int k, int n, double *arr, int *iarr)
{
int i,ir,j,l,mid,ia,itmp;
double a,tmp;
arr--;
iarr--;
l = 1;
ir = n;
for (;;) {
if (ir <= l+1) {
if (ir == l+1 && arr[ir] < arr[l]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
}
return;
} else {
mid=(l+ir) >> 1;
SWAP(arr[mid],arr[l+1])
ISWAP(iarr[mid],iarr[l+1])
if (arr[l] > arr[ir]) {
SWAP(arr[l],arr[ir])
ISWAP(iarr[l],iarr[ir])
}
if (arr[l+1] > arr[ir]) {
SWAP(arr[l+1],arr[ir])
ISWAP(iarr[l+1],iarr[ir])
}
if (arr[l] > arr[l+1]) {
SWAP(arr[l],arr[l+1])
ISWAP(iarr[l],iarr[l+1])
}
i = l+1;
j = ir;
a = arr[l+1];
ia = iarr[l+1];
for (;;) {
do i++; while (arr[i] < a);
do j--; while (arr[j] > a);
if (j < i) break;
SWAP(arr[i],arr[j])
ISWAP(iarr[i],iarr[j])
}
arr[l+1] = arr[j];
arr[j] = a;
iarr[l+1] = iarr[j];
iarr[j] = ia;
if (j >= k) ir = j-1;
if (j <= k) l = i;
}
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeAcklandAtom::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}
diff --git a/src/version.h b/src/version.h
index c66d90571..c061af61c 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "5 Jun 2010"
+#define LAMMPS_VERSION "13 Jun 2010"