program lattice2.f
implicit none
c This program is to calculate the surface modes of a 2-d harmonic lattice
c with nearest neighbor interactions (springs with constant k1) and next
c nearest neighbor interactions (springs with spring constant k2) for
c any given crystal momentum q
c Inputs are k1,k2,q
c w the real part of the frequency
c The compilation command under gcc is
c
c gfortran lattice2.f
c
real q, k1 ,k2 ,eta ,pi, w ,s ,c
complex tr
c s sin(q)
c c cos(q)
c eta is the imaginary part of frequensy
c tr is the function giving the trace of a matrix
parameter(k1=1., k2=0.5, eta=0.01, pi=2.*asin(1.) )
integer i,j,k
c g is greens function
c Du is hoping up matrix
c Dd is hopoing down matrix
c Do is the dynamical matrix
c Del is the correction to Do for the surfacew term
c z is the complex frequency
c a,b,e,ee,f,g1,g2 an auxiliary matrices
c inv a function that inverts a 2x2 matrix
complex g(2,2), Du(2,2), Dd(2,2), D(2,2), z, a(2,2), b(2,2)
complex f(2,2), ee(2,2), e(2,2) , g1(2,2), g2(2,2), Del(2,2)
complex g3(2,2), g4(2,2)
print*, "give a value for q"
read*, q
c find sine and cosine of q
c=cos(q)
s=sin(q)
c initializing the matrices to zero
do i=1,2
do j=1,2
a(i,j)=(0.,0.)
b(i,j)=(0.,0.)
e(i,j)=(0.,0.)
ee(i,j)=(0.,0.)
f(i,j)=(0.,0.)
g(i,j)=(0.,0.)
g1(i,j)=(0.,0.)
g2(i,j)=(0.,0.)
g3(i,j)=(0.,0.)
g4(i,j)=(0.,0.)
enddo
enddo
c determining the dynamical matrices of the system
D(1,1)=cmplx(-2.*(k1+k2-k1*c),0.)
Du(1,1)=cmplx(k2*c,0.)
Dd(1,1)=Du(1,1)
Del(1,1)=cmplx(-k2,0.)
D(2,2)= cmplx(-2.*(k1+k2),0.)
Du(2,2)=cmplx(k2*c+k1,0.)
Dd(2,2)=Du(2,2)
Del(2,2)=cmplx(-k2-k1,0.)
D(2,1)= (0.,0.)
Du(2,1)=cmplx(0.,-k2*s)
Dd(1,2)=-Du(2,1)
Del(1,2)=(0.,0.)
D(1,2)= (0.,0.)
Du(1,2)=cmplx(0.,-k2*s)
Dd(2,1)=-Du(1,2)
Del(2,1)=(0.,0.)
c opening file to write the results
open(10,file='g.dat')
c$$$ call inv(g,a)
c$$$ call mult(g,a,b)
c$$$
c$$$ do i=1,2
c$$$ do j=1,2
c$$$
c$$$ print*, D(i,j), Dd(i,j), Du(i,j)
c$$$
c$$$ enddo
c$$$ enddo
c print*, tr(D)
c initialize fraquency
w=-0.0001
c frequency loop
do while (w .le. 3.)
w= w + 0.001
z=cmplx(w, eta)
c print*, z**2
c a= -(w+i eta)**2-D
do i=1,2
do j=1,2
if (i .eq. j) then
a(i,j)=-z**2-D(i,j)
else
a(i,j)=-D(i,j)
endif
enddo
enddo
c loop to find g for this w
do i=1,500
do j=1,2
do k=1,2
b(j,k)=(0.,0.)
ee(j,k)=(0.,0.)
e(j,k)=(0.,0.)
f(j,k)=(0.,0.)
g1(j,k)=(0.,0.)
g2(j,k)=(0.,0.)
enddo
enddo
c b=Du*g
call mult(Du,g,b)
c f=a-b=-(x=i eta)^2-D-Dug
do j=1,2
do k=1,2
f(j,k)=a(j,k)-b(j,k)
enddo
enddo
c e=f^-1=(-(x=i eta)^2-D-Dug)^-1
call inv(f,e)
c ee=e*Dd=(-(x=i eta)^2-D-Dug)^-1*Dd
call mult(e,Dd,ee)
c Storing the value of g
do j=1,2
do k=1,2
g(j,k)=ee(j,k)
enddo
enddo
enddo
call inv(Dd,g1)
call mult(g,g1,g2)
c Correcting g to obtain the surface term
call inv(g2,g3)
do j=1,2
do k=1,2
g3(j,k)=g3(j,k)+Del(j,k)
enddo
enddo
call inv(g3,g4)
write(10,6) w , g4(1,1), g4(2,2)
6 format(f10.4, 2e18.7, 2e18.7)
enddo
stop
end
subroutine inv(A, B)
implicit none
c This subroutine inverts a complex 2x2 matrix A and stores the inverse to
c B
c det is the determinant of A
complex A(2,2), B(2,2), det
det= A(1,1)*A(2,2)-A(1,2)*A(2,1)
B(1,1)= A(2,2)/det
B(2,2)= A(1,1)/det
B(1,2)=-A(1,2)/det
B(2,1)=-A(2,1)/det
end
subroutine mult(A, B, C)
implicit none
c This subroutine multiplies a complex 2x2 matrix A with another compex
c 2x2 martix and stores the product in C
complex A(2,2), B(2,2), C(2,2)
integer i,j,k
do i=1,2
do j=1,2
do k=1,2
C(i,j)= C(i,j)+A(i,k)*B(k,j)
enddo
enddo
enddo
end
complex function tr(a)
implicit none
c This is a function giving the trace of an 2x2 matrix
complex a(2,2)
tr= a(1,1) + a(2,2)
end