onyx.f90

计算光子晶体能带的程序

! For the transmission / reflection calculation over-ride n_pts_store
n_pts_store=4*2*ixmax*iymax

! Allocate arrays for fields etc

allocate(e(3,0:ixmax+1,0:iymax+1,0:izmax+1))
allocate(h(3,0:ixmax+1,0:iymax+1,0:izmax+1))
allocate(fft_data(n_pts_store,itmax))
allocate(spectrum(fft_size))
allocate(store_pts(n_pts_store,4))
spectrum=0.0
fft_data=0.0

! Define the unit vectors for the co-ordinate system

! Default
u1=(/1.0,0.0,0.0/)
u2=(/0.0,1.0,0.0/)
u3=(/0.0,0.0,1.0/)

! Co-ordinate system for test calculation
!u1=(/1.0,0.0,0.0/)
!u2=(/0.0,1.0,0.0/)
!u3=(/0.5,0.0,sqrt(3.0)/2.0/)

call defmetric(u1,u2,u3,g1,g2,g3,g,omega)

! Initialise position dependent eps & mu

allocate (eps(ixmax,iymax,izmax))
allocate (mu(ixmax,iymax,izmax))
allocate (sigma(ixmax,iymax,izmax))
allocate (sigma_m(ixmax,iymax,izmax))

call defcell(eps,mu,sigma,sigma_m,u1,u2,u3)

! Define the effectve eps**-1, mu**-1 for the generalised co-ordinates

allocate (eps_inv(3,3,ixmax,iymax,izmax))
allocate (mu_inv(3,3,ixmax,iymax,izmax))
allocate (eps_hat(3,3,ixmax,iymax,izmax))
allocate (mu_hat(3,3,ixmax,iymax,izmax))

call renorm(eps,mu,sigma,sigma_m,eps_inv,mu_inv,eps_hat,mu_hat,g,omega)

! Default loop definitions
!k_loop: do ik=0,ikmax-1,1
!k_loop2: do ik2=0,ikmax-1,1
!k_loop3: do ik3=0,ikmax-1,1

! Loop required for test band structure calculation
!k_loop: do ik=0,ikmax-1,1
! Loops required for test DOS calcualtion
!k_loop: do ik=0,ikmax-1,1
!k_loop3: do ik3=0,ikmax-1,1

! Set the value for akx, aky, akz. akx=kx*(Q1*ixmax) etc.

! Default definitions of akx, aky, akz
!if (ikmax>1) akx=pi*ik/real(ikmax-1)
!if (ikmax>1) aky=pi*ik2/real(ikmax-1) 
!if (ikmax>1) akz=pi*ik/real(ikmax-1)

! Definitions of akx, akz for test band structure calculation: Gamma - M
!if (ikmax>1) akx=pi*ik/real(ikmax-1) 
!if (ikmax>1) akz=0.0

! Definitions of akx, akz for test band structure calculation: Gamma - K
!if (ikmax>1) akx=-(2.0/3.0)*pi*ik/real(ikmax-1) 
!if (ikmax>1) akz=(2.0/3.0)*pi*ik/real(ikmax-1)

! Definitions of akx, akz for test band structure calculation: M - K
!if (ikmax>1) akx=pi-(1.0/3.0)*pi*ik/real(ikmax-1) 
!if (ikmax>1) akz=-(2.0/3.0)*pi*ik/real(ikmax-1)

! Definitions of akx, akz required for test DOS calculation
!if (ikmax>1) akx=pi*ik/real(ikmax-1)
!if (ikmax>1) akz=pi*ik3/real(ikmax-1)

! ix,iy,iz define the point r - For the DOS we calc G(r,r)
! ipol determines which component of the fields we are looking at

iz_cur=bclayer1+20
!iz_cur=1
iy_cur=1
ix_cur=1
i_pol=1

! Loop over i_pol is required for the DOS calculation
!pol_loop: do i_pol=1,6

! Initialise store_pts which hold the locations of the points
! at which we store the fields

! For the band structure
!call init_store_pts_band(store_pts)

! For trans/ref
call init_store_pts_trans(store_pts)

! For the DOS
!call init_store_pts_dos(store_pts,ix_cur,iy_cur,iz_cur,i_pol)

! Initialise fields

! For the band structure
!call initfields(g1,g2,g3,e,h,eps_hat,mu_hat,ix_cur,iy_cur,iz_cur,i_pol)

! For the trans/ref
call initfields_trans(g1,g2,g3,e,h,eps_hat,mu_hat,ix_cur,iy_cur, &
&                     iz_cur,i_pol)

! For the DOS
!call initfields_dos(g1,g2,g3,e,h,eps_hat,mu_hat,ix_cur,iy_cur,iz_cur,i_pol)

! Transfer control to central driver

write(logfile,*) 'Begin main calculation'
call driver(e,h,eps_inv,mu_inv,eps_hat,mu_hat,sigma,sigma_m, &
&           fft_data,store_pts)
write(logfile,*) 'End main calculation'

! Postprocess if necessary

! For the band structure
!call postproc_band(fft_data,spectrum)

! For the DOS

!call postproc_dos(fft_data,spectrum)

! End of loop over i_pol
!enddo pol_loop

! End of loops over k
!enddo k_loop3
!enddo k_loop2
!enddo k_loop

! Tidy up!
deallocate(e,h,eps,mu,sigma,sigma_m,eps_hat,mu_hat)

! Final postprocessing if necessary

! For the DOS
!call finalproc_dos(spectrum)

! For trans/ref
call postproc_trans(fft_data,eps_inv,mu_inv)

write(logfile,*) 'Execution terminated normally'

! Tidy up!
close(unit=logfile)
close(unit=outfile)
close(unit=fields)
close(unit=ref)
close(unit=trans)
close(unit=infile)

deallocate(spectrum)
deallocate(fft_data,eps_inv,mu_inv)
deallocate(store_pts)

stop
end program onyx

! -----------------------------------------------------------------------------
! Subroutine driver
! Heart of the calculation. This subroutine actually does the time integration
! loop. Should be made to be as flexible as possible
! -----------------------------------------------------------------------------
subroutine driver(e_cur,h_cur,eps_inv,mu_inv,eps_hat,mu_hat,sigma,sigma_m, &
&                 fft_data,store_pts)
use interface2
use interface3
use interface5
use interface6
use parameters
use files
use physconsts
implicit none

complex,pointer :: e_cur(:,:,:,:)
complex,pointer :: h_cur(:,:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:)
complex,pointer :: mu_inv(:,:,:,:,:)
complex,pointer :: eps_hat(:,:,:,:,:)
complex,pointer :: mu_hat(:,:,:,:,:)
complex,pointer :: fft_data(:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
integer,pointer :: store_pts(:,:)

complex,pointer :: e_prev(:,:,:,:)
complex,pointer :: h_prev(:,:,:,:)
complex,pointer :: div(:,:,:)
real,pointer :: rho(:,:,:)
integer :: iblock
integer :: i

complex,pointer :: intcurle_1(:,:,:),intcurlh_1(:,:,:)
complex,pointer :: intcurle_2(:,:,:),intcurlh_2(:,:,:)
real,pointer :: wz_1(:),wz_2(:)

allocate(e_prev(3,0:ixmax+1,0:iymax+1,0:izmax+1))
allocate(h_prev(3,0:ixmax+1,0:iymax+1,0:izmax+1))
allocate(div(ixmax,iymax,bclayer1:izmax-bclayer2))
allocate(rho(ixmax,iymax,bclayer1:izmax-bclayer2))

! Need these for the absorbing layers
allocate(wz_1(bclayer1))
allocate(intcurle_1(ixmax,iymax,bclayer1))
allocate(intcurlh_1(ixmax,iymax,bclayer1))
allocate(wz_2(bclayer2))
allocate(intcurle_2(ixmax,iymax,bclayer2))
allocate(intcurlh_2(ixmax,iymax,bclayer2))

! Initialize 2nd set of fields

e_prev=(0.0,0.0)
h_prev=(0.0,0.0)

intcurle_1=(0.0,0.0)
intcurlh_1=(0.0,0.0)
intcurle_2=(0.0,0.0)
intcurlh_2=(0.0,0.0)

! Initialize the absorbing boundary layer PML

call set_wz(wz_1,wz_2)

!write(fields,*)
!write(fields,*) (i,real(h_cur(1,1,1,i)),i=1,izmax)

maintime: do iblock=1,n_block

!write(logfile,'(A10,I10)') ' Time step',1+(iblock-1)*block_size

call int(e_cur,h_cur,e_prev,h_prev,eps_inv,mu_inv,sigma,sigma_m,     &
&        block_size,iblock,fft_data,store_pts,intcurle_1,intcurle_2, &
&        intcurlh_1,intcurlh_2,wz_1,wz_2)

! Various tests on the fields can be performed here...

! Update the remain boundaries prior to calculating divs

!call bc_xmax_metal(e_cur,h_cur)
!call bc_ymax_metal(e_cur,h_cur)
!call bc_zmax_metal(e_cur,h_cur)
call bc_xmax_bloch(e_cur,h_cur)
call bc_ymax_bloch(e_cur,h_cur)
call bc_zmax_bloch(e_cur,h_cur)

! Test for energy conservation

call calc_energy_density(e_cur,e_prev,h_prev,rho,eps_hat,mu_hat)
write(logfile,*) 'Total energy ',sum(rho(:,:,:))

! Test for conservation of div D, div B

!call calc_div_B(h_cur,mu_hat,div)
!write(logfile,*) 'Sum div B ', sum(abs(div(:,:,:)))
!call calc_div_D(e_cur,eps_hat,div)
!write(logfile,*) 'Sum div D ', sum(abs(div(:,:,:)))

! Write out the fields to a file

!write(fields,*)
!write(fields,*) (i,real(h_cur(1,1,1,i)),i=1,izmax)

!write(fields,*)
!write(fields,*) (i,-real(h_cur(2,1,1,i)),i=1,izmax)

enddo maintime

deallocate(e_prev,h_prev)
deallocate(div,rho)
deallocate(wz_1,wz_2)
deallocate(intcurle_1,intcurle_2)
deallocate(intcurlh_1,intcurlh_2)

return
end subroutine driver

! -----------------------------------------------------------------------------
! Subroutine int
! Do the time integration for nt time steps
! -----------------------------------------------------------------------------
subroutine int(e_cur,h_cur,e_prev,h_prev,eps_inv,mu_inv,sigma,sigma_m, &
&              nt,iblock,fft_data,store_pts,intcurle_1,intcurle_2,     &
&              intcurlh_1,intcurlh_2,wz_1,wz_2)
use interface3
use interface5
use interface6
use physconsts
use parameters
implicit none

integer,intent(in) :: nt,iblock
complex,pointer :: e_cur(:,:,:,:)
complex,pointer :: e_prev(:,:,:,:)
complex,pointer :: h_cur(:,:,:,:)
complex,pointer :: h_prev(:,:,:,:)
complex,pointer :: eps_inv(:,:,:,:,:)
complex,pointer :: mu_inv(:,:,:,:,:)
complex,pointer :: fft_data(:,:)
real,pointer :: sigma(:,:,:),sigma_m(:,:,:)
integer,intent(in) :: store_pts(:,:)

complex,pointer :: intcurle_1(:,:,:),intcurlh_1(:,:,:)
complex,pointer :: intcurle_2(:,:,:),intcurlh_2(:,:,:)
real,intent(in) :: wz_1(:),wz_2(:)

integer :: ix,iy,iz,it,i_pts
complex :: curl(3)
real :: dtcq2
complex,pointer :: temp(:,:,:,:)

dtcq2=(dt*c0/Q0)**2

time: do it=1,nt

! Current fields become previous fields

temp=>e_cur
e_cur=>e_prev
e_prev=>temp

temp=>h_cur
h_cur=>h_prev
h_prev=>temp

! First the E-fields

do iz=1+bclayer1,izmax-bclayer2
   do iy=1,iymax
      do ix=1,ixmax

! Define Curl H

         curl(1)=h_prev(3,ix,iy,iz)-h_prev(3,ix,iy-1,iz) &
&               -h_prev(2,ix,iy,iz)+h_prev(2,ix,iy,iz-1)
         curl(2)=h_prev(1,ix,iy,iz)-h_prev(1,ix,iy,iz-1) &
&               -h_prev(3,ix,iy,iz)+h_prev(3,ix-1,iy,iz)
         curl(3)=h_prev(2,ix,iy,iz)-h_prev(2,ix-1,iy,iz) &
&               -h_prev(1,ix,iy,iz)+h_prev(1,ix,iy-1,iz)

! Integrate fields in time

         e_cur(1,ix,iy,iz)=(1.0-sigma(ix,iy,iz))*e_prev(1,ix,iy,iz)+  &
&        (eps_inv(1,1,ix,iy,iz)*curl(1)  &
&        +eps_inv(1,2,ix,iy,iz)*curl(2)  &
&        +eps_inv(1,3,ix,iy,iz)*curl(3))

         e_cur(2,ix,iy,iz)=(1.0-sigma(ix,iy,iz))*e_prev(2,ix,iy,iz)+  &
&        (eps_inv(2,1,ix,iy,iz)*curl(1)  &
&        +eps_inv(2,2,ix,iy,iz)*curl(2)  &
&        +eps_inv(2,3,ix,iy,iz)*curl(3))

         e_cur(3,ix,iy,iz)=(1.0-sigma(ix,iy,iz))*e_prev(3,ix,iy,iz)+  &
&        (eps_inv(3,1,ix,iy,iz)*curl(1)  &
&        +eps_inv(3,2,ix,iy,iz)*curl(2)  &
&        +eps_inv(3,3,ix,iy,iz)*curl(3))

      enddo
   enddo
enddo

! Update Berenger type PML equations for the E-field (Absorbing boundary layer)

call PML_e(e_cur,e_prev,h_prev,intcurlh_1,intcurlh_2,wz_1,wz_2,eps_inv)

! Update the boundary conditions needed for the E-fields

!call bc_xmax_metal(e_cur,h_cur)
!call bc_ymax_metal(e_cur,h_cur)
!call bc_zmax_metal(e_cur,h_cur)
call bc_xmax_bloch(e_cur,h_cur)
call bc_ymax_bloch(e_cur,h_cur)
call bc_zmax_bloch(e_cur,h_cur)

! Then the H-fields

do iz=1+bclayer1,izmax-bclayer2
   do iy=1,iymax
      do ix=1,ixmax

! Define Curl E

         curl(1)=e_cur(3,ix,iy+1,iz)-e_cur(3,ix,iy,iz) &
&               -e_cur(2,ix,iy,iz+1)+e_cur(2,ix,iy,iz)
         curl(2)=e_cur(1,ix,iy,iz+1)-e_cur(1,ix,iy,iz) &
&               -e_cur(3,ix+1,iy,iz)+e_cur(3,ix,iy,iz)
         curl(3)=e_cur(2,ix+1,iy,iz)-e_cur(2,ix,iy,iz) &
&               -e_cur(1,ix,iy+1,iz)+e_cur(1,ix,iy,iz)

! Integrate the fields in time

         h_cur(1,ix,iy,iz)=1.0/(1.0+sigma_m(ix,iy,iz))* &
&        (h_prev(1,ix,iy,iz)-dtcq2*                     &
&        (mu_inv(1,1,ix,iy,iz)*curl(1)                  &
&        +mu_inv(1,2,ix,iy,iz)*curl(2)                  &
&        +mu_inv(1,3,ix,iy,iz)*curl(3)))

         h_cur(2,ix,iy,iz)=1.0/(1.0+sigma_m(ix,iy,iz))* &
&        (h_prev(2,ix,iy,iz)-dtcq2*                     &
&        (mu_inv(2,1,ix,iy,iz)*curl(1)                  &
&        +mu_inv(2,2,ix,iy,iz)*curl(2)                  &
&        +mu_inv(2,3,ix,iy,iz)*curl(3)))

         h_cur(3,ix,iy,iz)=1.0/(1.0+sigma_m(ix,iy,iz))* &
&        (h_prev(3,ix,iy,iz)-dtcq2*                     &
&        (mu_inv(3,1,ix,iy,iz)*curl(1)                  &
&        +mu_inv(3,2,ix,iy,iz)*curl(2)                  &
&        +mu_inv(3,3,ix,iy,iz)*curl(3)))

      enddo
   enddo
enddo

! Update Berenger type PML equations for the H-field (Absorbing boundary layer)

call PML_h(e_cur,h_cur,h_prev,intcurle_1,intcurle_2,wz_1,wz_2,mu_inv)

! Update the boundary conditions needed for the H-fields

!call bc_xmin_metal(e_cur,h_cur)
!call bc_ymin_metal(e_cur,h_cur)
!call bc_zmin_metal(e_cur,h_cur)
call bc_xmin_bloch(e_cur,h_cur)
call bc_ymin_bloch(e_cur,h_cur)
call bc_zmin_bloch(e_cur,h_cur)

! Store  points for Fourier transform later

do i_pts=1,n_pts_store
  if (store_pts(i_pts,1)<4) then
    fft_data(i_pts,it+(iblock-1)*nt)=e_cur(store_pts(i_pts,1),     &
&            store_pts(i_pts,2),store_pts(i_pts,3),store_pts(i_pts,4))
  else
    fft_data(i_pts,it+(iblock-1)*nt)=h_cur(store_pts(i_pts,1)-3,  &
&            store_pts(i_pts,2),store_pts(i_pts,3),store_pts(i_pts,4))
  endif
enddo

enddo time

return
end subroutine int

! -----------------------------------------------------------------------------
! Subroutine setparam
! Set the key parameters from a config file
! -----------------------------------------------------------------------------
subroutine setparam()
use parameters
use physconsts
use files
implicit none
integer :: nbits,set_bits,i

write(logfile,*) 'Read in parameters'

read(infile,*) ixmax
if (ixmax<1) call err(2)
write(logfile,'(A1,A12,I10)') '#','ixmax=',ixmax
write(outfile,'(A1,A12,I10)') '#','ixmax=',ixmax
write(fields,'(A1,A12,I10)') '#','ixmax=',ixmax
write(trans,'(A1,A12,I10)') '#','ixmax=',ixmax
write(ref,'(A1,A12,I10)') '#','ixmax=',ixmax

read(infile,*) iymax
if (iymax<1) call err(2)
write(logfile,'(A1,A12,I10)') '#','iymax=',iymax
write(outfile,'(A1,A12,I10)') '#','iymax=',iymax
write(fields,'(A1,A12,I10)') '#','iymax=',iymax
write(trans,'(A1,A12,I10)') '#','iymax=',iymax
write(ref,'(A1,A12,I10)') '#','iymax=',iymax