onyx.f90

计算光子晶体能带的程序

! Defines the dielectric cylinder for the test calculation
! -----------------------------------------------------------------------------
subroutine defcell(eps,mu,sigma,sigma_m,u1,u2,u3)
use files
use parameters
use physconsts
implicit none

complex,pointer :: eps(:,:,:),mu(:,:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
real :: u1(3),u2(3),u3(3)

complex :: epsset,eps1,eps2
integer :: ix,iy,iz,i,j,l,m,isub,iatom,natom
real :: asub,test,r(3),radius,g(3,3),vol
real :: r1,r2,atom(3,17)

write(logfile,*) 'Define unit cell'

natom=14
!bclayer1=0
!bclayer2=0
bclayer1=50
bclayer2=50
epsref=(9.0,0.0)

eps=epsref
mu=(1.0,0.0)
sigma=(0.0,0.0)
sigma_m=(0.0,0.0)

g(1,1)=u1(1)*u1(1)+u1(2)*u1(2)+u1(3)*u1(3)
g(1,2)=u1(1)*u2(1)+u1(2)*u2(2)+u1(3)*u2(3)
g(1,3)=u1(1)*u3(1)+u1(2)*u3(2)+u1(3)*u3(3)
g(2,1)=u2(1)*u1(1)+u2(2)*u1(2)+u2(3)*u1(3)
g(2,2)=u2(1)*u2(1)+u2(2)*u2(2)+u2(3)*u2(3)
g(2,3)=u2(1)*u3(1)+u2(2)*u3(2)+u2(3)*u3(3)
g(3,1)=u3(1)*u1(1)+u3(2)*u1(2)+u3(3)*u1(3)
g(3,2)=u3(1)*u2(1)+u3(2)*u2(2)+u3(3)*u2(3)
g(3,3)=u3(1)*u3(1)+u3(2)*u3(2)+u3(3)*u3(3)

! Set up the default parameters for the unit cell

! radius = ratio of the radius of cylinder to the unit cell
! eps2 = dielectric constant for cylinder
! eps1 = dielectric constant for background medium
! isub = number of sub-divisions in average over each lattice point

eps1=(9.0,0.0)
eps2=(1.0,0.0)
radius=0.46
!radius=0.406698
isub=50

r1=radius**2

vol=0.0
! TEST PATCH
!!$atom(1,1)=0.5
!!$atom(2,1)=0.0
!!$atom(3,1)=0.5

atom(1,1)=0.5
atom(2,1)=0.0
atom(3,1)=0.5

atom(1,2)=0.5
atom(2,2)=0.0
atom(3,2)=1.5

atom(1,3)=0.5
atom(2,3)=0.0
atom(3,3)=2.5

atom(1,4)=0.5
atom(2,4)=0.0
atom(3,4)=3.5

atom(1,5)=0.0
atom(2,5)=0.0
atom(3,5)=0.0

atom(1,6)=0.0
atom(2,6)=0.0
atom(3,6)=1.0

atom(1,7)=1.0
atom(2,7)=0.0
atom(3,7)=0.0

atom(1,8)=1.0
atom(2,8)=0.0
atom(3,8)=1.0

atom(1,9)=0.0
atom(2,9)=0.0
atom(3,9)=2.0

atom(1,10)=1.0
atom(2,10)=0.0
atom(3,10)=2.0

atom(1,11)=0.0
atom(2,11)=0.0
atom(3,11)=3.0

atom(1,12)=1.0
atom(2,12)=0.0
atom(3,12)=3.0

atom(1,13)=0.0
atom(2,13)=0.0
atom(3,13)=4.0

atom(1,14)=1.0
atom(2,14)=0.0
atom(3,14)=4.0


do ix=1,ixmax
!do iz=1,izmax
do iz=izmax/2+1-40,izmax/2+4*ixmax-40
  asub=1.0/float(isub)
  test=0.0
  epsset=0.0
  do l=1,isub
  do m=1,isub
    do iatom=1,natom
      r(1)=Q(1)*((ix-1      +(l-1)*asub)/float(ixmax) -atom(1,iatom))
      r(2)=0.0
!      r(3)=Q(3)*((iz-1      +(m-1)*asub)/float(ixmax) -atom(3,iatom))
      r(3)=Q(3)*((iz-(izmax/2-40)+(m-1)*asub)/float(ixmax) -atom(3,iatom))
      r2=0.0
      do i=1,3
      do j=1,3
        r2=r2+(g(i,j)*r(i)*r(j))
      enddo
      enddo
      if (r2.lt.r1) epsset=epsset+1.0
      if (r2.lt.r1) test=test+1.0
    enddo
  enddo
  enddo
  epsset=epsset/real(isub**2)
  test=test/real(isub**2)
  do iy=1,iymax
    eps(ix,iy,iz)=epsset*eps2+(1.0-epsset)*eps1
  enddo
  vol=vol+test
enddo
enddo

write(logfile,*) 'test f=',vol/(4.0*ixmax*iymax*ixmax) &
&               ,'compare f=',2.0*pi*radius**2/(Q(1)*Q(2)*Q(3))

write(trans,'(A1,A12,F15.5)') '#','eps1=',real(eps1)
write(trans,'(A1,A12,F15.5)') '#','eps2=',real(eps2)
write(ref,'(A1,A12,F15.5)') '#','eps1=',real(eps1)
write(ref,'(A1,A12,F15.5)') '#','eps2=',real(eps2)
write(trans,'(A1,A12,F15.5)') '#','test f=',vol/(4.0*ixmax*iymax*ixmax)
write(ref,'(A1,A12,F15.5)') '#','test f=',vol/(4.0*ixmax*iymax*ixmax)
write(trans,'(A1,A12,F15.5)') '#','compare f=' &
&                             ,2.0*pi*radius**2/(Q(1)*Q(2)*Q(3))
write(ref,'(A1,A12,F15.5)') '#','compare f=' &
&                             ,2.0*pi*radius**2/(Q(1)*Q(2)*Q(3))

return
end subroutine defcell

! -----------------------------------------------------------------------------
! Subroutine init_store_pts_band
! Initialize the random points for storing the fields for the FFT
! For the band structure store fields at random set of points
! -----------------------------------------------------------------------------
subroutine init_store_pts_band(store_pts)
use files
use parameters
implicit none

integer,pointer :: store_pts(:,:)

integer :: i
real :: x

write(logfile,*) 'Initialize store_pts in init_store_pts_band'

do i=1,n_pts_store
  call random_number(x)
  store_pts(i,1)=floor(3.0*x+4.0)
  call random_number(x)
  store_pts(i,2)=floor(ixmax*x+1.0)
  call random_number(x)
  store_pts(i,3)=floor(iymax*x+1.0)
  call random_number(x)
  store_pts(i,4)=floor(izmax*x+1.0)
  if (store_pts(i,1)<1.or.store_pts(i,1)>6         &
&     .or.store_pts(i,2)<1.or.store_pts(i,2)>ixmax &
&     .or.store_pts(i,3)<1.or.store_pts(i,3)>iymax &
&     .or.store_pts(i,4)<1.or.store_pts(i,4)>izmax) call err(8)
enddo

return
end subroutine init_store_pts_band

! -----------------------------------------------------------------------------
! Subroutine init_store_pts_dos
! Initialize the random points for storing the fields for the FFT
! For the Green's function store the fields at r' (for G(w,r,r'))
! For the DOS r=r'=(ix_cur,iy_cur,iz_cur)
! -----------------------------------------------------------------------------
subroutine init_store_pts_dos(store_pts,ix_cur,iy_cur,iz_cur,i_pol)
use files
implicit none

integer,pointer :: store_pts(:,:)
integer,intent(in) :: ix_cur,iy_cur,iz_cur,i_pol

write(logfile,*) 'Initialize store_pts in init_store_pts_dos'

store_pts=1

store_pts(1,1)=i_pol
store_pts(1,2)=ix_cur
store_pts(1,3)=iy_cur
store_pts(1,4)=iz_cur

return
end subroutine init_store_pts_dos

! -----------------------------------------------------------------------------
! Subroutine init_store_pts_trans
! Initialize the random points for storing the fields for the FFT
! For the trans/ref store fields at the ends of the system
! -----------------------------------------------------------------------------
subroutine init_store_pts_trans(store_pts)
use files
use parameters
implicit none

integer,pointer :: store_pts(:,:)

integer :: i,ix,iy,ip

write(logfile,*) 'Initialize store_pts in init_store_pts_trans'

do ix=1,ixmax
do iy=1,iymax
do ip=1,4
  i=ip+4*(ix-1+ixmax*(iy-1))
if (ip<3) then
  store_pts(i,1)=ip
else
  store_pts(i,1)=ip+1
endif
  store_pts(i,2)=ix
  store_pts(i,3)=iy
  store_pts(i,4)=bclayer1+50
enddo
enddo
enddo

do ix=1,ixmax
do iy=1,iymax
do ip=1,4
  i=4*ixmax*iymax + ip+4*(ix-1+ixmax*(iy-1))
if (ip<3) then
  store_pts(i,1)=ip
else
  store_pts(i,1)=ip+1
endif
  store_pts(i,2)=ix
  store_pts(i,3)=iy
  store_pts(i,4)=izmax-bclayer2-25
enddo
enddo
enddo

return
end subroutine init_store_pts_trans

! -----------------------------------------------------------------------------
! Subroutine defmetric
! Define the metric and various other factors needed for
! generalised co-ordinate systems
! -----------------------------------------------------------------------------
subroutine defmetric(u1,u2,u3,g1,g2,g3,g,omega)
use parameters
use matrices
use files
implicit none

real,intent(in) :: u1(3),u2(3),u3(3)
real,intent(out) :: g1(3),g2(3),g3(3)
real,intent(out) :: g(3,3),omega

real :: tmp(3)

write(logfile,*) 'Define the metric'

! g is the metric. (g**-1)ij = ui.uj
! ui's are the real space lattice vectors while
! the gi's are the reciprocal lattice vectors 
! (apart from the 2*pi/Qi factor!)

! omega=|u1.u2xu3|

tmp=cross_prod(u2,u3)
omega=dot_product(u1,tmp)

! Define the reciprocal vectors

g1=tmp/omega
g2=cross_prod(u3,u1)/omega
g3=cross_prod(u1,u2)/omega

! Define the metric g

g(1,1)=dot_product(g1,g1)
g(1,2)=dot_product(g1,g2)
g(1,3)=dot_product(g1,g3)
g(2,1)=dot_product(g2,g1)
g(2,2)=dot_product(g2,g2)
g(2,3)=dot_product(g2,g3)
g(3,1)=dot_product(g3,g1)
g(3,2)=dot_product(g3,g2)
g(3,3)=dot_product(g3,g3)

! Q0 is scale factor equal to a typical length

Q0=Q(3)

return
end subroutine defmetric

! -----------------------------------------------------------------------------
! Subroutine renorm
! Sets up eps_hat, mu_hat, eps_inv, mu_inv, sigma and sigma_m
! based on the metric for our co-ordinate system
! -----------------------------------------------------------------------------
subroutine renorm(eps,mu,sigma,sigma_m,eps_inv,mu_inv,eps_hat,mu_hat,g, &
&                 omega)
use matrices
use parameters
use physconsts
use files
implicit none

complex,pointer :: eps(:,:,:),mu(:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:)
complex,pointer :: mu_inv(:,:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:)
complex,pointer :: mu_hat(:,:,:,:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
real,intent(in) :: g(3,3),omega

complex, dimension(3,3) :: eps_tmp,mu_tmp
integer :: i,j,ix,iy,iz
logical :: fail

complex :: EQ(3),HQ(3)

write(logfile,*) 'Define eps_inv, mu_inv for general co-ordinates'

do ix=1,ixmax
do iy=1,iymax
do iz=1,izmax

! Define effective eps & mu for the generalised co-ordineates

   EQ(:)=Q(:)
   HQ(:)=Q(:)

   do i=1,3
   do j=1,3
     eps_tmp(i,j)=eps(ix,iy,iz)*g(i,j)*omega*HQ(1)*HQ(2)*HQ(3)/(HQ(i)*EQ(j)*Q0)
     mu_tmp(i,j)=mu(ix,iy,iz)*g(i,j)*omega*EQ(1)*EQ(2)*EQ(3)/(EQ(i)*HQ(j)*Q0)
   enddo
   enddo

! Invert effective eps & mu

   eps_inv(1:3,1:3,ix,iy,iz)=matinv3(eps_tmp,emach,fail)
   if (fail) call err(3)
   mu_inv(1:3,1:3,ix,iy,iz)=matinv3(mu_tmp,emach,fail)
   if (fail) call err(3)

! Store the non-inverted eps & mu (used in calculating the div's)

   eps_hat(1:3,1:3,ix,iy,iz)=eps_tmp
   mu_hat(1:3,1:3,ix,iy,iz)=mu_tmp

! Define the effective condutivity

   sigma(ix,iy,iz)=sigma(ix,iy,iz)*dt/(eps0*eps(ix,iy,iz))
   sigma_m(ix,iy,iz)=sigma_m(ix,iy,iz)*dt/(mu0*mu(ix,iy,iz))

enddo
enddo
enddo

return
end subroutine renorm

! -----------------------------------------------------------------------------
! Subroutine postproc_band
! Post-process the stored time dependent fields to extract the bands at
! that k-point
! -----------------------------------------------------------------------------
subroutine postproc_band(fft_data,spectrum)
use files
use physconsts
use parameters
use interface4
implicit none

complex,pointer :: fft_data(:,:)
real,pointer :: spectrum(:)

integer :: i

write(logfile,*) 'Postprocessing in postproc_band'

! Call FFT routine to calculate the power spectrum
! For the Band Structure

spectrum=0
call power_spec(fft_data,spectrum)

if (ikmax==1) then

! Write out the spectrum

  do i=1,fft_size
    write(outfile,*) i,(spectrum(i))
  enddo

else

! Write out the Band structure

  do i=2,fft_size-1
    if (spectrum(i)-spectrum(i-1)>1.0e+9*emach.and.&
&       spectrum(i+1)-spectrum(i)<-1.0e+9*emach.and.&
&       (i/(dt*2.0*fft_size))*(Q(1)*ixmax)/c0<2.4) then

! Write out frequencies in dimensionless units

! Default
!      write(outfile,*) akx, &
!&                      (i/(dt*2.0*fft_size))*(Q(1)*ixmax)/c0

! Gamma - M for the triangular lattice
      write(outfile,*) akz, &
&                      (i/(dt*2.0*fft_size))*(Q(1)*ixmax)/c0

! Gamma - K for the triangular lattice
!      write(outfile,*) sqrt(3.0)*akx, &
!&                      (i/(dt*2.0*fft_size))*(Q(1)*ixmax)/c0

! M - K for the triangular lattice
!      write(outfile,*) pi-(3.0/4.0)*akz, &
!&                      (i/(dt*2.0*fft_size))*(Q(1)*ixmax)/c0

! Write out frequencies in electron volts
!      write(outfile,*) akx,(27.2*2.0*pi*i/(dt*2*fft_size))

    endif
  enddo

endif

return
end subroutine postproc_band

! -----------------------------------------------------------------------------
! Subroutine postproc_dos
! Post-process the stored time dependent fields to calculate contribution
! to the density of states
! -----------------------------------------------------------------------------