```
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"

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

```