Codes for constructing matrix of Hamiltonion of 2D square Tight-bonding model

Hello,
Does anyone know how to write codes for constructing the matrix of Hamiltonion of 2D square Tight-bonding model. Then diagonalize it.
Any language is OK. (Fortran is best)

 PhysOrg.com physics news on PhysOrg.com >> Promising doped zirconia>> Nanocrystals grow from liquid interface>> New insights into how materials transfer heat could lead to improved electronics
 Hello, This program shows the principle of diagonalizing a matrix in fortran, using a 2D geometry. To compile, write gfortran -llapack 2D.f90 then write ./a.out to run the program Or ifort -llapack 2D.f90 Code: module shared !every function and routine which implicit none !includes the module shared can use its real*8, allocatable :: Hdiag(:,:),eigenvalues(:) !variables end module shared program diag use shared !now the program diag knows about Hdiag and eigenvalues implicit none !this makes sure one has do define everything explicitly integer :: nsize !size of the matrix to be diagonalized integer :: i,j,icount,jcount,norb integer :: i1,i2,j1,j2 real*8 :: v norb = 5 !size of system nsize = norb*norb !size of hamiltonian allocate(Hdiag(nsize,nsize)) !This is the matrix to be diagonalized allocate(eigenvalues(nsize)) !after diagonalization, this vector !contains the eigenvalues v = -1.0 !this is the hopping parameter !Now we construct the hamiltonian, and we do that by first looping over !j, and then i. This are non-periodic boundary conditions. icount = 0 jcount = 0 do i1=1,norb do j1=1,norb jcount = 0 icount = icount + 1 do i2=1,norb do j2=1,norb jcount = jcount + 1 if (abs(i1-i2) + abs(j1-j2) == 1)Hdiag(icount,jcount) = v !if distance == 1, put hopping parameter = v write(11,*)'hop(',icount,jcount,')=',Hdiag(icount,jcount) !this prints the hamiltonian to disc end do end do end do end do call diagonalize(nsize) !now Hdiag will be diagonalized. !After diagonalization, Hdiag(:,:) will contain the eigenvectors, !Where Hdiag(:,1) is the first eigenvector, Hdiag(:,2) the second and so on !and eigenvalues(1) will contain the eigenvalue corresponding to !eigenvector 1, and so on. do i=1,nsize write(*,*)'eigenvalue:',i,eigenvalues(i) !writing eigenvalues to screen end do end program diag !program is finished subroutine diagonalize(ndiag) !Diagonalizes the matrix Hdiag use shared implicit none real*8, allocatable :: work(:) integer, allocatable :: iwork(:) character job, uplo external dsyevd integer nmax, lda, lwork, liwork, info,ndiag nmax=ndiag lda=nmax lwork=4*nmax*nmax+100 liwork=10*nmax allocate(work(lwork),iwork(liwork)) job='v' uplo='l' call dsyevd(job,uplo,ndiag,Hdiag,lda,eigenvalues,work,lwork,iwork,liwork,info) if (info > 0) then write (*,*) 'failure to converge.' stop endif deallocate(work,iwork) return end