From 4334b5423c88ff688a746eee164a218c003f0708 Mon Sep 17 00:00:00 2001 From: QuanSheng Wu Date: Wed, 20 Sep 2017 14:01:35 +0200 Subject: [PATCH] add bulkband_plane_calc for Dirac cone plot --- soc/ek_bulk2D.f90 | 183 +++++++++++++++++++++++++++++++++++++++++-- soc/fermisurface.f90 | 138 +++++++++++++++++++++++--------- soc/main.f90 | 13 +++ soc/module.f90 | 4 +- soc/readinput.f90 | 4 +- 5 files changed, 294 insertions(+), 48 deletions(-) diff --git a/soc/ek_bulk2D.f90 b/soc/ek_bulk2D.f90 index 70ea1038..07fc3663 100644 --- a/soc/ek_bulk2D.f90 +++ b/soc/ek_bulk2D.f90 @@ -1,6 +1,175 @@ -! calculate bulk's energy band using wannier TB method -! Copyright (c) 2010 QuanSheng Wu. All rights reserved. - subroutine ek_bulk2D + subroutine ek_bulk_plane + ! Calculate bulk's energy band using wannier TB method for a given k plane + ! Copyright (c) 2017 QuanSheng Wu. All rights reserved. + ! Revised by QS.Wu on Sep 19 2017 @EPFL, Switzerland + ! This subroutine is usefull for Dirac cone plot + + use wmpi + use para + + implicit none + + integer :: ik, i, j, ib + integer :: knv3 + integer :: ierr + integer :: nwann + + integer :: nband_min + integer :: nband_max + integer :: nband_store + + real(dp) :: time_start, time_end + + real(Dp) :: k(3) + real(Dp) :: W(Num_wann) + real(dp) :: kxmin, kxmax, kymin, kymax + real(dp), allocatable :: kxy(:,:) + real(dp), allocatable :: kxy_shape(:,:) + real(dp), allocatable :: kxy_plane(:,:) + + ! Hamiltonian of bulk system + complex(Dp) :: Hamk_bulk(Num_wann,Num_wann) + + ! eigen value of H + real(dp), allocatable :: gap(:) + real(dp), allocatable :: gap_mpi(:) + real(dp), allocatable :: eigv(:,:) + real(dp), allocatable :: eigv_mpi(:,:) + + nband_min= Numoccupied- 1 + nband_max= Numoccupied+ 2 + if (nband_min< 1) nband_min= 1 + if (nband_max> Num_wann) nband_max= Num_wann + nband_store= nband_max- nband_min+ 1 + + + knv3= nk1*Nk2 + allocate( kxy(3, nk1*Nk2)) + allocate( kxy_shape(3, nk1*Nk2)) + allocate( kxy_plane(3, nk1*Nk2)) + + allocate( eigv (nband_store, knv3)) + allocate( eigv_mpi(nband_store, knv3)) + eigv = 0d0 + eigv_mpi= 0d0 + + kxy=0d0 + kxy_shape=0d0 + kxy_plane=0d0 + + if (Nk1<2 .or. Nk2<2) stop 'ERROR: I refuse to do this job because you give me so small Nk1 and Nk2' + if (Numoccupied> Num_wann) stop 'ERROR: please set correct Numoccupied, it should be small than num_wann' + + ik =0 + do i= 1, Nk1 + do j= 1, Nk2 + ik =ik +1 + kxy(:, ik)= K3D_start+ K3D_vec1*(i-1)/dble(Nk1-1)+ K3D_vec2*(j-1)/dble(Nk2-1) & + -(K3D_vec1+K3D_vec2)/2d0 + kxy_shape(:, ik)= kxy(1, ik)* Kua+ kxy(2, ik)* Kub+ kxy(3, ik)* Kuc + call rotate_k3_to_kplane(kxy_shape(:, ik), kxy_plane(:, ik)) + enddo + enddo + + time_start= 0d0 + time_end= 0d0 + do ik= 1+cpuid, knv3, num_cpu + if (cpuid==0.and. mod(ik/num_cpu, 500)==0) & + write(stdout, '(a, i12, a, i12, a, f10.2, a)') & + 'ek_bulk_plane, ik ', ik, ' knv3',knv3, ' time left', & + (knv3-ik)*(time_end- time_start)/num_cpu, ' s' + call now(time_start) + + k = kxy(:, ik) + + ! calculation bulk hamiltonian + Hamk_bulk= 0d0 + call ham_bulk(k, Hamk_bulk) + + !> diagonalization by call zheev in lapack + W= 0d0 + call eigensystem_c( 'N', 'U', Num_wann ,Hamk_bulk, W) + eigv(:, ik)= W(nband_min:nband_max) + + call now(time_end) + enddo ! ik + +#if defined (MPI) + call mpi_allreduce(eigv,eigv_mpi,size(eigv),& + mpi_dp,mpi_sum,mpi_cmw,ierr) +#else + eigv_mpi= eigv +#endif + + + !> write out the data in gnuplot format + outfileindex= outfileindex+ 1 + if (cpuid==0)then + open(unit=outfileindex, file='bulkek_plane.dat') + write(outfileindex, '(1000a19)')'# kx', 'ky', 'kz', 'k1', 'k2', & + 'E(Numoccupied-1)', 'E(Numoccupied)' , 'E(Numoccupied+1)', 'E(Numoccupied+2)' + do ik=1, knv3 + write(outfileindex, '(1000f19.9)')kxy_shape(:, ik), & + kxy_plane(:, ik), eigv_mpi(:, ik) + if (mod(ik, nk2)==0) write(outfileindex, *)' ' + enddo + close(outfileindex) + endif + + !> write out the data in gnuplot format + outfileindex= outfileindex+ 1 + if (cpuid==0)then + open(unit=outfileindex, file='bulkek_plane-matlab.dat') + write(outfileindex, '(1000a19)')'% kx', 'ky', 'kz', 'k1', 'k2', & + 'E(Numoccupied-1)', 'E(Numoccupied)' , 'E(Numoccupied+1)', 'E(Numoccupied+2)' + do ik=1, knv3 + write(outfileindex, '(1000f19.9)')kxy_shape(:, ik), & + kxy_plane(:, ik), eigv_mpi(:, ik) + enddo + close(outfileindex) + endif + + !> write out a script that can be used for gnuplot + outfileindex= outfileindex+ 1 + if (cpuid==0)then + open(unit=outfileindex, file='bulkek_plane.gnu') + write(outfileindex, '(a)')"set encoding iso_8859_1" + write(outfileindex, '(a)')'#set terminal postscript enhanced color' + write(outfileindex, '(a)')"#set output 'bulkek_plane.eps'" + write(outfileindex, '(3a)')'set terminal png truecolor enhanced', & + ' size 1920, 1680 font ",36"' + write(outfileindex, '(a)')"set output 'bulkek_plane.png'" + write(outfileindex, '(a)')'set palette rgbformulae 33,13,10' + write(outfileindex, '(a)')'unset key' + write(outfileindex, '(a)')'set pm3d' + write(outfileindex, '(a)')'set origin 0.2, 0' + write(outfileindex, '(a)')'set size 0.8, 1' + write(outfileindex, '(a)')'set border lw 3' + write(outfileindex, '(a)')'#set xtics font ",24"' + write(outfileindex, '(a)')'#set ytics font ",24"' + write(outfileindex, '(a)')'set size ratio -1' + write(outfileindex, '(a)')'set xtics' + write(outfileindex, '(a)')'set ytics' + write(outfileindex, '(a)')'set view 80,60' + write(outfileindex, '(a)')'set xlabel "k_1"' + write(outfileindex, '(a)')'set ylabel "k_2"' + write(outfileindex, '(a)')'set zlabel "Energy (eV)" rotate by 90' + write(outfileindex, '(a)')'unset colorbox' + write(outfileindex, '(a)')'set autoscale fix' + write(outfileindex, '(a)')'set pm3d interpolate 4,4' + write(outfileindex, '(2a)')"splot 'bulkek_plane.dat' u 4:5:8 w pm3d, \" + write(outfileindex, '(2a)')" 'bulkek_plane.dat' u 4:5:9 w pm3d" + + close(outfileindex) + + endif ! cpuid + + return + end subroutine ek_bulk_plane + + subroutine ek_bulk2D_spin + ! Calculate bulk's energy band using wannier TB method for a given k plane + ! Copyright (c) 2010 QuanSheng Wu. All rights reserved. use wmpi use para @@ -154,11 +323,11 @@ subroutine ek_bulk2D endif return - end subroutine ek_bulk2D + end subroutine ek_bulk2D_spin - !> using green's function - subroutine ek_bulk2D_spin + subroutine ek_bulk2D_spin_green + !> using green's function use wmpi use para @@ -347,4 +516,4 @@ subroutine ek_bulk2D_spin endif return - end subroutine ek_bulk2D_spin + end subroutine ek_bulk2D_spin_green diff --git a/soc/fermisurface.f90 b/soc/fermisurface.f90 index ba276218..deb5d5ac 100644 --- a/soc/fermisurface.f90 +++ b/soc/fermisurface.f90 @@ -21,6 +21,8 @@ subroutine fermisurface3D integer :: nband_min integer :: nband_max integer :: nband_store + + real(dp) :: time_start, time_end ! Hamiltonian of bulk system complex(Dp) :: Hamk_bulk(Num_wann,Num_wann) @@ -74,8 +76,13 @@ subroutine fermisurface3D allocate(eigval_mpi(nband_store, knv3)) eigval_mpi= 0d0 eigval= 0d0 + time_start= 0d0 + time_end= 0d0 do ik= 1+cpuid, knv3, num_cpu - if (cpuid==0.and. mod(ik,100)<3)write(stdout, *)'3DFS, ik, knv3', ik, knv3 + if (cpuid==0.and. mod(ik/num_cpu, 500)==0) & + write(stdout, *) '3DFS, ik ', ik, 'knv3',knv3, 'time left', & + (knv3-ik)*(time_end- time_start)/num_cpu, ' s' + call now(time_start) k(1) = kxyz(1, ik) k(2) = kxyz(2, ik) @@ -86,6 +93,7 @@ subroutine fermisurface3D call ham_bulk_old(k, Hamk_bulk) call eigensystem_c( 'N', 'U', Num_wann ,Hamk_bulk, W) eigval_mpi(:, ik)= W(nband_min:nband_max) + call now(time_end) enddo #if defined (MPI) @@ -154,24 +162,29 @@ subroutine fermisurface ! Hamiltonian of bulk system complex(Dp) :: Hamk_bulk(Num_wann,Num_wann) + real(dp) :: time_start, time_end real(dp) :: kxmin, kxmax, kymin, kymax real(dp) :: zmin, zmax real(dp) :: kxmin_shape, kxmax_shape, kymin_shape, kymax_shape real(dp), allocatable :: kxy(:,:) real(dp), allocatable :: kxy_shape(:,:) + real(dp), allocatable :: kxy_plane(:,:) real(dp), allocatable :: dos(:) real(dp), allocatable :: dos_mpi(:) complex(dp), allocatable :: ones(:,:) - nkx= Nk1 - nky= Nk2 + nkx= Nk + nky= Nk allocate( kxy(3, nkx*nky)) allocate( kxy_shape(3, nkx*nky)) + allocate( kxy_plane(3, nkx*nky)) kxy=0d0 kxy_shape=0d0 + kxy_plane=0d0 + ik =0 do i= 1, nkx @@ -179,6 +192,7 @@ subroutine fermisurface ik =ik +1 kxy(:, ik)= K3D_start+ K3D_vec1*(i-1)/dble(nkx-1)+ K3D_vec2*(j-1)/dble(nky-1) kxy_shape(:, ik)= kxy(1, ik)* Kua+ kxy(2, ik)* Kub+ kxy(3, ik)* Kuc + call rotate_k3_to_kplane(kxy_shape(:, ik), kxy_plane(:, ik)) enddo enddo @@ -188,6 +202,7 @@ subroutine fermisurface kxmax_shape=maxval(kxy_shape(i1,:)) kymin_shape=minval(kxy_shape(i2,:)) kymax_shape=maxval(kxy_shape(i2,:)) + knv3= nkx*nky @@ -202,20 +217,26 @@ subroutine fermisurface ones(i, i)= 1d0 enddo + time_start= 0d0 + time_end= 0d0 do ik= 1+cpuid, knv3, num_cpu - if (cpuid==0) write(stdout, *)'FS, ik, knv3' , ik, knv3 + if (cpuid==0.and. mod(ik/num_cpu, 500)==0) & + write(stdout, *) 'FS_Plane, ik ', ik, 'knv3',knv3, 'time left', & + (knv3-ik)*(time_end- time_start)/num_cpu, ' s' + call now(time_start) k = kxy(:, ik) ! calculation bulk hamiltonian Hamk_bulk= 0d0 - call ham_bulk_old(k, Hamk_bulk) + call ham_bulk(k, Hamk_bulk) Hamk_bulk= (E_arc -zi* eta_arc)* ones - Hamk_bulk call inv(Num_wann, Hamk_bulk) do i=1, Num_wann dos(ik)= dos(ik)+ aimag(Hamk_bulk(i, i))/pi enddo + call now(time_end) enddo @@ -229,9 +250,10 @@ subroutine fermisurface outfileindex= outfileindex+ 1 if (cpuid==0)then open(unit=outfileindex, file='fs.dat') + write(outfileindex, '(100a16)')'# kx', 'ky', 'kz', 'log(A(k,E))', 'kp1', 'kp2', 'kp3' do ik=1, knv3 - write(outfileindex, '(30f16.8)')kxy_shape(:, ik), log(dos_mpi(ik)) + write(outfileindex, '(30f16.8)')kxy_shape(:, ik), log(dos_mpi(ik)), kxy_plane(:, ik) if (mod(ik, nky)==0) write(outfileindex, *)' ' enddo close(outfileindex) @@ -241,39 +263,41 @@ subroutine fermisurface !> minimum and maximum value of energy bands + outfileindex= outfileindex+ 1 !> write script for gnuplot if (cpuid==0) then - open(unit=101, file='fs.gnu') - write(101, '(a)')"set encoding iso_8859_1" - write(101, '(a)')'#set terminal postscript enhanced color' - write(101, '(a)')"#set output 'fs.eps'" - write(101, '(3a)')'set terminal pngcairo truecolor enhanced', & + open(unit=outfileindex, file='fs.gnu') + write(outfileindex, '(a)')"set encoding iso_8859_1" + write(outfileindex, '(a)')'#set terminal postscript enhanced color' + write(outfileindex, '(a)')"#set output 'fs.eps'" + write(outfileindex, '(3a)')'set terminal pngcairo truecolor enhanced', & ' size 1920, 1680 font ",36"' - write(101, '(a)')"set output 'fs.png'" - write(101,'(a, f10.4, a, f10.4, a, f10.4, a)') & + write(outfileindex, '(a)')"set output 'fs.png'" + write(outfileindex,'(a, f10.4, 2a, f10.4, a)') & 'set palette defined ( ', zmin, ' "white", ', & - zmin+(zmin+zmax)/6d0, '"black", ', zmax,' "red" )' - write(101, '(a)')'#set palette rgbformulae 33,13,10' - write(101, '(a)')'unset ztics' - write(101, '(a)')'unset key' - write(101, '(a)')'set pm3d' - write(101, '(a)')'#set view equal xyz' - write(101, '(a)')'set view map' - write(101, '(a)')'set border lw 3' - write(101, '(a)')'#set xtics font ",24"' - write(101, '(a)')'#set ytics font ",24"' - write(101, '(a)')'set size ratio -1' - write(101, '(a)')'unset xtics' - write(101, '(a)')'unset ytics' - write(101, '(a)')'set colorbox' - !write(101, '(a, f10.5, a, f10.5, a)')'set xrange [', kxmin, ':', kxmax, ']' - !write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', kymin, ':', kymax, ']' - !write(101, '(a, f10.5, a, f10.5, a)')'set xrange [', kxmin_shape, ':', kxmax_shape, ']' - !write(101, '(a, f10.5, a, f10.5, a)')'set yrange [', kymin_shape, ':', kymax_shape, ']' - write(101, '(a)')'set pm3d interpolate 2,2' - write(101, '(2a)')"splot 'fs.dat' u 1:2:4 w pm3d" - - close(101) + '0 "black", ', zmax,' "red" )' + write(outfileindex, '(a)')'#set palette rgbformulae 33,13,10' + write(outfileindex, '(a)')'unset ztics' + write(outfileindex, '(a)')'unset key' + write(outfileindex, '(a)')'set pm3d' + write(outfileindex, '(a)')'#set view equal xyz' + write(outfileindex, '(a)')'set view map' + write(outfileindex, '(a)')'set border lw 3' + write(outfileindex, '(a)')'#set xtics font ",24"' + write(outfileindex, '(a)')'#set ytics font ",24"' + write(outfileindex, '(a)')'set size ratio -1' + write(outfileindex, '(a)')'unset xtics' + write(outfileindex, '(a)')'unset ytics' + write(outfileindex, '(a)')'set colorbox' + !write(outfileindex, '(a, f10.5, a, f10.5, a)')'set xrange [', kxmin, ':', kxmax, ']' + !write(outfileindex, '(a, f10.5, a, f10.5, a)')'set yrange [', kymin, ':', kymax, ']' + write(outfileindex, '(a)')'set autoscale fix' + !write(outfileindex, '(a, f10.5, a, f10.5, a)')'set xrange [', kxmin_shape, ':', kxmax_shape, ']' + !write(outfileindex, '(a, f10.5, a, f10.5, a)')'set yrange [', kymin_shape, ':', kymax_shape, ']' + write(outfileindex, '(a)')'set pm3d interpolate 2,2' + write(outfileindex, '(2a)')"splot 'fs.dat' u 1:2:4 w pm3d" + + close(outfileindex) endif @@ -371,7 +395,7 @@ subroutine gapshape3D ! calculation bulk hamiltonian Hamk_bulk= 0d0 - call ham_bulk_old(k, Hamk_bulk) + call ham_bulk(k, Hamk_bulk) call eigensystem_c( 'N', 'U', Num_wann ,Hamk_bulk, W) gap(1, ik)= W(Numoccupied+1)- W(Numoccupied) @@ -563,8 +587,8 @@ subroutine gapshape write(outfileindex, '(100a16)')'% kx', 'ky', 'kz', 'gap', 'Ev2', 'Ev1', 'Ec1', & 'Ec2', 'k1', 'k2', 'k3' do ik=1, knv3 - if (abs(gap_mpi(1, ik))< 0.10d0) then - write(outfileindex, '(8f16.8)')kxy_shape(:, ik), (gap_mpi(:, ik)), kxy(:, ik) + if (abs(gap_mpi(1, ik))< Gap_threshold) then + write(outfileindex, '(80f16.8)')kxy_shape(:, ik), (gap_mpi(:, ik)), kxy(:, ik) endif enddo close(outfileindex) @@ -600,7 +624,7 @@ subroutine gapshape write(outfileindex, '(a)')"set output 'GapPlane.png'" write(outfileindex,'(a, f10.4, a, f10.4, a, f10.4, a)') & 'set palette defined ( ', zmin, ' "black", ', & - zmin+ (zmin+zmax)/20d0,' "orange", ',zmax,' "white" )' + (zmin+zmax)/20d0,' "orange", ',zmax,' "white" )' write(outfileindex, '(a)')"set origin 0.10, 0.0" write(outfileindex, '(a)')"set size 0.85, 1.0" write(outfileindex, '(a)')'unset ztics' @@ -778,5 +802,41 @@ function fermi(omega, Beta) result(value) return end function fermi + subroutine rotate_k3_to_kplane(k3, kplane) + use para, only : dp, K3D_vec1, K3D_vec2, Kua, Kub, Kuc + implicit none + + real(dp), intent(in) :: k3(3) + real(dp), intent(out) :: kplane(3) + !> three new unit vectors + real(dp) :: Urot_t(3, 3) + + real(dp) :: kvec1(3), kvec2(3) + + real(dp), external :: norm + + kvec1= K3D_vec1(1)*Kua+ K3D_vec1(2)*Kub+ K3D_vec1(3)*Kuc + kvec2= K3D_vec2(1)*Kua+ K3D_vec2(2)*Kub+ K3D_vec2(3)*Kuc + + Urot_t(1, :)= kvec1/norm(kvec1) + + !> e_z' + Urot_t(3, 1)= (kvec1(2)*kvec2(3)- kvec1(3)*kvec2(2)) + Urot_t(3, 2)= (kvec1(3)*kvec2(1)- kvec1(1)*kvec2(3)) + Urot_t(3, 3)= (kvec1(1)*kvec2(2)- kvec1(2)*kvec2(1)) + Urot_t(3, :)= Urot_t(3, :)/norm(Urot_t(3, :)) + + !> e_y'= e_z'\cross e_x' + Urot_t(2, 1)= (Urot_t(3, 2)*Urot_t(1, 3)- Urot_t(3, 3)*Urot_t(1, 2)) + Urot_t(2, 2)= (Urot_t(3, 3)*Urot_t(1, 1)- Urot_t(3, 1)*Urot_t(1, 3)) + Urot_t(2, 3)= (Urot_t(3, 1)*Urot_t(1, 2)- Urot_t(3, 2)*Urot_t(1, 1)) + Urot_t(2, :)= Urot_t(2, :)/norm(Urot_t(2, :)) + + kplane(1)= Urot_t(1, 1)*k3(1)+ Urot_t(1, 2)*k3(2)+ Urot_t(1, 3)*k3(3) + kplane(2)= Urot_t(2, 1)*k3(1)+ Urot_t(2, 2)*k3(2)+ Urot_t(2, 3)*k3(3) + kplane(3)= Urot_t(3, 1)*k3(1)+ Urot_t(3, 2)*k3(2)+ Urot_t(3, 3)*k3(3) + + return + end subroutine diff --git a/soc/main.f90 b/soc/main.f90 index 83853657..64f814e6 100644 --- a/soc/main.f90 +++ b/soc/main.f90 @@ -140,6 +140,19 @@ program main if(cpuid.eq.0)write(stdout, *)'<< End of calculating bulk band' endif + + !> bulk band in a plane. For Dirac or Weyl cone + if (BulkBand_plane_calc) then + if(cpuid.eq.0)write(stdout, *)' ' + if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the bulk band in plane' + call now(time_start) + call ek_bulk_plane + call now(time_end) + call print_time_cost(time_start, time_end, 'BulkBand_plane') + if(cpuid.eq.0)write(stdout, *)'<< End of calculating the bulk band in plane' + endif + + !> Find nodes in BZ if (FindNodes_calc) then if(cpuid.eq.0)write(stdout, *)' ' diff --git a/soc/module.f90 b/soc/module.f90 index 500ec5bf..bc6e300a 100644 --- a/soc/module.f90 +++ b/soc/module.f90 @@ -26,6 +26,7 @@ module para !> control parameters logical :: BulkBand_calc ! Flag for bulk energy band calculation + logical :: BulkBand_plane_calc ! Flag for bulk energy band calculation logical :: BulkFS_calc ! Flag for bulk 3D fermi surface in 3D BZ calculation logical :: BulkFS_plane_calc ! Flag for bulk fermi surface for a fix k plane calculation logical :: BulkGap_cube_calc ! Flag for Gap_cube calculation @@ -53,7 +54,8 @@ module para SlabSS_calc, SlabArc_calc, SlabSpintexture_calc, & WannierCenter_calc,BerryPhase_calc,BerryCurvature_calc, & Z2_3D_calc, Chern_3D_calc, WeylChirality_calc, & - Dos_calc, JDos_calc, EffectiveMass_calc, SlabQPI_calc, FindNodes_calc + Dos_calc, JDos_calc, EffectiveMass_calc, SlabQPI_calc, FindNodes_calc, & + BulkBand_plane_calc ! double precision integer,parameter :: Dp=kind(1.0d0) diff --git a/soc/readinput.f90 b/soc/readinput.f90 index e11a5dda..4460df7e 100644 --- a/soc/readinput.f90 +++ b/soc/readinput.f90 @@ -54,6 +54,7 @@ subroutine readinput BulkBand_calc = .FALSE. + BulkBand_plane_calc = .FALSE. BulkFS_calc = .FALSE. BulkGap_cube_calc = .FALSE. BulkGap_plane_calc = .FALSE. @@ -85,7 +86,7 @@ subroutine readinput write(*, *)"SlabSpintexture,wanniercenter_calc" write(*, *)"BerryPhase_calc,BerryCurvature_calc, Z2_3D_calc" write(*, *)"Dos_calc, JDos_calc, FindNodes_calc " - write(*, *)"BulkFS_plane_calc" + write(*, *)"BulkFS_plane_calc, BulkBand_plane_calc" write(*, *)"Z2_3D_calc" write(*, *)"Chern_3D_calc" write(*, *)"WeylChirality_calc" @@ -99,6 +100,7 @@ subroutine readinput write(stdout, *) " " write(stdout, *) ">>>Control parameters: " write(stdout, *) "BulkBand_calc : ", BulkBand_calc + write(stdout, *) "BulkBand_plane_calc : ", BulkBand_plane_calc write(stdout, *) "BulkFS_calc : ", BulkFS_calc write(stdout, *) "BulkFS_Plane_calc : ", BulkFS_Plane_calc write(stdout, *) "BulkGap_cube_calc : ", BulkGap_cube_calc