Back
to ABE
!-------------------------------------------------
! Thermal Diffuse Scattering from Mettalic Sodium
! "Diffuse X-ray Reflections from Crystals" ! W. A. Wooster, Dorver, 1999
!------------------------------------------------
PROGRAM Huang
integer, parameter :: dp = selected_real_kind(10)
Real (kind=dp), dimension ( : , : ), allocatable :: TDS, HDS, a, b
Real, dimension ( : , : ), allocatable :: x, y, z
Real (kind=dp), dimension ( : ), allocatable :: c, fp, fpp, Fr
real, dimension (3) :: h
REAL (kind=dp) :: pi, dg, phi, DW, Fh2, Freal, Fimag, Vc, AW, Btotal
Real (kind=dp) :: S, DT, Temp, comp, a0, b0, c0, c11, c12, c44, hm
REAL :: dX, dY, dZ
integer :: nX, nY, nZ, nAtom, nA, nWL
Character (Len=50) :: DR1, DR2, File_output, File_f0, File_f12
pi = 3.1415926535
!------------------------------------------------------
! nWL=1 : Cr(Ka) ! nWL=2 : Fe(Ka) ! nWL=3 : Cu(Ka) ! nWL=4 : Mo(Ka) ! nWL=5 : Ag(Ka) nWL=4 ! WL = 1.542 ! Wave Length [A]
WL = 0.71069 ! Wave Length [A] !-------------------------------------------------
DR1 = 'B:\QB\X-ray\Asf\'
DR2 = 'B:\ABE\Na\'
! file_output='Na310hk'; h(1)=3.0; h(2)=1.0; h(3)=0.0 ! (h1, h2, h3) reciprocal lattice vector !
file_output = 'Na220hk'; h(1)=2.0; h(2)=2.0; h(3)=0.0 ! (h1, h2, h3) reciprocal lattice vector
file_output = 'Na420hk'; h(1)=4.0; h(2)=2.0; h(3)=0.0 ! (h1, h2, h3) reciprocal lattice vector
file_output = trim(DR2) // trim(file_output) // '.TDS'
!--------------------------------------------------
! (1-c)*A + c*B
comp = 0.00 ! Composition
a0=4.29; b0=4.29; c0=4.29 ! Lattice constant [A]
Vc = a0*b0*c0 ! Volume of the unit cell [A^3]
DT = 158.0 ! Debye Temperature [K]
AW = 22.9898 ! Atomic Weight of Sodium ! AW = (1-comp)*AW1 + comp*AW2
c11=7.39 ! [10^10 dyne/cm^2]
c12=6.22
c44=4.19 ! elastic constant of Metallic sodium
Temp = 290.0 ! Setting temperature [K]
!---------------------------------
! Phi(x) determination
dg = 0.001
phi = phi_X( DT, Temp, dg )
Btotal = 11492.0 * Temp / AW / DT**2 * phi + 2873 / AW / DT !B-factor [A^2]
write(*,'(A, F10.4)') ' Debye Temperature =',DT
write(*,'(A, F10.4)') ' T =',Temp
write(*,'(A, F10.4)') ' Atomic Weight =',AW
write(*,'(A, F10.4)') ' composition =',comp
write(*,'(A, F10.4)') ' phi =',phi
write(*,'(A, F10.4)') ' B(total) =',Btotal
!-----------------------------------------------------
! Read the Atomic Scattering Factor
nAtom = 1
ALLOCATE( a(nAtom, 4), b(nAtom, 4), c(nAtom), fp(nAtom), fpp(nAtom), Fr(nAtom) )
File_f0 = trim(DR1) // 'Na.f0'
File_f12= trim(DR1) // 'Na.f12'
nA=1
CALL f_open (File_f0, File_f12, a, b, c, fp, fpp, nAtom, nA, nWL)
write(*,*) '------------------------------'
write(*,*) ' Atomic Scattering Factor'
write(*,*) '------------------------------'
do i=1,nA
do j=1,4
write(*,'(2F12.6)') a(i,j), b(i,j)
end do
write(*,'(F12.6)') c(i)
write(*,'(2F12.6)') fp(i), fpp(i)
end do
!------------------------------------
! Calculate Form Factor
! s = sqrt( (h(1)/a0)**2 + (h(2)/b0)**2 + (h(3)/c0)**2 )
! s = sin(the)/WL for Orthorombic
hm = h(1)**2 + h(2)**2 + h(3)**2
s = sqrt( hm )/2.0/a0
! s = sin(the)/WL for Cubic
! Atomic Scattering Factor
do k=1,nAtom
Fr(k)=0.0
do l=1,4
Fr(k) = Fr(k) + a(k,l)*exp( -b(k,l)*s**2 )
end do
Fr(k) = Fr(k) + c(k) + fp(k)
end do
Freal = Fr(1)*( 1.0 + (-1)**( h(1)+h(2)+h(3) ) )
! Structure Factor for bcc
Fimag = (1-comp)*fpp(1) + comp*fpp(2)
!---------------------------------------------------------
! Debye-Waller Factor exp(-M) = exp(-B*(sin(the)/WL)**2)
DW = exp(-Btotal*s**2)
Fh2 = ( Freal**2 + Fimag**2 ) * DW**2
!-----------------------------------------------------
! (h, k, 0) plane
!-----------------------------------------------------
write(*,*) '----------------------------------------'
write(*,'(A, 3F10.3)') ' ( h, k, l )',h(1),h(2),h(3)
write(*,*) ' File name = ', file_output
write(*,*)
nX=31; nY=31; nZ=1 n1=nX; n2=nY
dX=0.01; dY=0.01; dZ=0.000
ALLOCATE( TDS(n1, n2), HDS(n1, n2), x(n1, n2), y(n1, n2), z(n1, n2) )
CALL Index_hk (h, X, Y, Z, dX, dY, dZ, nX, nY, nZ, n1, n2)
!----------------------------------------------------
! TDS calculation
!----------------------------------------------------
write(*,*) 'start TDS calculation'
READ(*,*) aaaaa
CALL TDS_2D(TDS, Temp, n1, n2, a0, b0, c0, h, X, Y, Z, c11, c12, c44, Fh2)
!----------------------------------------
! Output DATA file
!----------------------------------------
open(8, FILE = file_output)
do i=1,n1
do j=1,n2
TDS(i,j) = TDS(i,j) / 1000.0
write(*,'(3F10.5, F15.6)') x(i,j),y(i,j),z(i,j),TDS(i,j)
write(8,'(3F10.5, F15.6)') x(i,j),y(i,j),z(i,j),TDS(i,j)
end do
end do
close(8)
END PROGRAM Huang
!************************************************************************
subroutine f_open (file_f0, file_f12, a, b, c, fp, fpp, nAtom, nA, nWL)
!************************************************************************
integer, parameter :: dp = selected_real_kind(10)
integer,intent(in) :: nAtom, nA, nWL
character (Len=50) , intent(in) :: File_f0, File_f12
real (kind=dp), dimension ( : , : ), intent(OUT) :: a(nAtom, 4), b(nAtom, 4)
real (kind=dp), dimension( : ), intent(OUT) :: c(nAtom), fp(nAtom), fpp(nAtom)
character (Len=50) :: dummy
open( 2, FILE=file_f0 )
do i=1,4
read(2,*) a(nA, i)
read(2,*) b(nA, i)
end do
read(2,*) c(nA)
close(2)
open( 2, FILE=file_f12 )
do i = 1, 5
read(2,*) f1, f2, dummy
if (nWL== i) then
fp(nA) =f1
fpp(nA)=f2
end if
end do
close(2)
end subroutine f_open
!****************************************
FUNCTION phi_X(DT,Temp,dg) result(phi)
!****************************************
integer, parameter :: dp = selected_real_kind(10)
real (kind=dp), intent(in) :: DT, Temp, dg
real(kind=dp) :: X
g = 0.0
phi = 0.0
if (Temp > 0.0) then
X = DT / Temp
do
g=anint(1000.0*(g+dg))/1000.0
phi = phi + g/(exp(g)-1)*dg
if (g>=X) exit
end do
phi = phi / X
end if
end function phi_X
!******************************************************
SUBROUTINE Index_hk (h,X,Y,Z,dX,dY,dZ,nX,nY,nZ,n1,n2)
!******************************************************
real, dimension ( : ), intent(in) :: h(3)
Integer, intent (in) :: nX, nY, nZ, n1, n2
real, intent (in) :: dX, dY, dZ
real, dimension ( : , : ), intent(OUT) :: X(n1, n2), Y(n1, n2), Z(n1, n2)
Xs = h(1) - (nX-1)/2.0*dX
Ys = h(2) - (nY-1)/2.0*dY
Zs = h(3) - (nZ-1)/2.0*dZ
x1=Xs; x2=Ys ;x3=Zs+dZ
do i=1, n1
do j=1, n2
x1=anint(1000*x1)/1000
x2=anint(1000*x2)/1000
x3=anint(1000*x3)/1000
X(i,j)=x1; Y(i,j)=x2; Z(i,j)=x3
! write(*,*) X(i,j),Y(i,j),Z(i,j)
x1=x1+dX
end do
x1=Xs
x2=x2+dY
! read(*,*) aaaaa
end do
end subroutine Index_hk
!*******************************************************************
SUBROUTINE TDS_2D(TDS,Temp,n1,n2,a0,b0,c0,h,X,Y,Z,c11,c12,c44,Fh2)
!*******************************************************************
integer, parameter :: dp = selected_real_kind(10)
real (kind=dp), dimension( : , : ), intent(OUT) :: TDS(n1, n2)
real, dimension ( : , : ), intent(IN) :: X(n1, n2),Y(n1, n2),Z(n1, n2)
real, dimension( : ), intent(IN) :: h(3)
integer, intent ( in ) :: n1, n2
real (kind=dp), intent(in) :: Temp, a0, b0, c0, c11, c12, c44, Fh2
real (kind=dp), dimension( : ) :: q(3), uq(3), uh(3)
real (kind=dp) :: hm, qm
DO i=1, n1
DO j=1, n2
q(1)=anint(1000*(X(i, j)-h(1)))/1000
q(2)=anint(1000*(Y(i, j)-h(2)))/1000
q(3)=anint(1000*(Z(i, j)-h(3)))/1000
qm=sqrt( q(1)**2 + q(2)**2 + q(3)**2 )
if (qm/=0.0) then
hm = h(1)**2 + h(2)**2 + h(3)**2
hm = sqrt(hm)
do k=1, 3
uq(k)=q(k)/qm ! unit vector of q
uh(k)=h(k)/hm ! unit vector of h
end do
TDS(i, j)=Temp/a0/b0/c0*Fh2 * TDS_calc(c11,c12,c44,uh,uq)*(hm/qm)**2
else
TDS(i, j)=0.0
end if
! write(*,*) X(i, j), Y(i, j), Z(i, j), TDS(i, j)
end do
! read(*,*) aaa
end do
END SUBROUTINE TDS_2D
!************************************************
FUNCTION TDS_calc(c11,c12,c44,uh,uq) result(Sq)
!************************************************
integer, parameter :: dp=selected_real_kind(10)
real (kind=dp), dimension ( : ), intent(in) :: uq(3), uh(3)
real (kind=dp), intent(in) :: c11, c12, c44
real (kind=dp) :: dA
real (kind=dp) :: A11, A22, A33, A12, A13, A23
real (kind=dp) :: Ai11, Ai22, Ai33, Ai12, Ai13, Ai23
!--------------------
! A-matrix
! | A11 A12 A13 |
! | A12 A22 A23 |
! | A13 A23 A33 |
!------------------------------------
! A(i,j) calculation for Tetragonal
! A11 = c11*uq(1)**2 + c66*uq(2)**2 + c44*uq(3)**2
! A22 = c66*uq(1)**2 + c11*uq(2)**2 + c44*uq(3)**2
! A33 = c44*uq(1)**2 + c44*uq(2)**2 + c33*uq(3)**2
! A12 = (c12+c66)*uq(1)*uq(2)
! A13 = (c13+c44)*uq(1)*uq(3)
! A23 = (c13+c44)*uq(2)*uq(3)
!---------------------------------
! A(i,j) calculation for Cubic
A11 = c11*uq(1)**2 + c44*uq(2)**2 + c44*uq(3)**2
A22 = c44*uq(1)**2 + c11*uq(2)**2 + c44*uq(3)**2
A33 = c44*uq(1)**2 + c44*uq(2)**2 + c11*uq(3)**2
A12 = (c12+c44)*uq(1)*uq(2)
A13 = (c12+c44)*uq(1)*uq(3)
A23 = (c12+c44)*uq(2)*uq(3)
!----------------------------------
! The determinant of the A-matrix
dA = A11*A22*A33 + 2*A12*A23*A13 - A11*A23**2 - A22*A13**2 - A33*A12**2
!--------------------------------------------------------------------------
! A(i,j)^(-1)
Ai11 = ( A22*A33 - A23**2 )
Ai22 = ( A11*A33 - A13**2 )
Ai33 = ( A11*A22 - A12**2 )
Ai12 = -( A33*A12 - A13*A23 )
Ai13 = -( A11*A23 - A12*A13 )
Ai23 = ( A12*A23 - A22*A13 )
Sq = Ai11*uh(1)**2 + Ai22*uh(2)**2 + Ai33*uh(3)**2 &
& +2*(Ai12*uh(1)*uh(2) + Ai23*uh(2)*uh(3) + Ai13*uh(3)*uh(1) )
Sq = Sq / dA
end function TDS_calc
References
Back
to Department of Materials Science and Engineering
Last Modified: April 1, 2009