EPW Forum

Code: Select all

  !-----------------------------------------------------------------------
  PROGRAM main
  !-----------------------------------------------------------------------
  !
  IMPLICIT NONE
  !
  DOUBLE PRECISION :: delta_min, delta_max, dbin, e, w, z, zn
  INTEGER :: ik, ibnd, nbin, ibin, iufil, j
  INTEGER :: nkftot, nkf1, nkf2, nkf3, nkfs, n, nbnd_, nbndfs
  DOUBLE PRECISION :: arg, ef, ef0, dosef, degaussw, fsthick, w0gauss, weight, w0gauss2, arg2
  DOUBLE PRECISION, ALLOCATABLE :: wkfs(:), xkfs(:,:), ekf_(:,:), ekfs(:,:), w0g(:,:), delta(:,:), delta_k_bin(:)
  DOUBLE PRECISION, ALLOCATABLE :: b(:,:), x(:), dist(:)
  DOUBLE PRECISION, PARAMETER :: ryd2ev = 13.6056981d0, zero = 0.d0
  CHARACTER (len=256) :: word
  !
  iufil = 10
  OPEN(unit=iufil, file="egnv", status='unknown', form='unformatted')
  READ(iufil) nkftot, nkf1, nkf2, nkf3, nkfs
  READ(iufil) nbnd_, ef, ef0, dosef, degaussw, fsthick
  !
  degaussw = degaussw * ryd2ev
  ef0 = ef0 * ryd2ev
  ef = ef * ryd2ev
  fsthick = fsthick * ryd2ev
  dosef = dosef / ryd2ev
  !
  WRITE(*,'(5i7)') nkftot, nkf1, nkf2, nkf3, nkfs
  WRITE(*,'(i7,5d24.15)') nbnd_, ef, ef0, dosef, degaussw, fsthick
  !
  IF (.not. ALLOCATED(wkfs)) ALLOCATE(wkfs(nkfs))
  IF (.not. ALLOCATED(xkfs)) ALLOCATE(xkfs(3, nkfs))
  wkfs(:) = 0.d0
  xkfs(:, :) = 0.d0
  IF ( .not. ALLOCATED(ekf_) ) ALLOCATE(ekf_(nbnd_, nkfs))
  ekf_(:, :) = 0.d0
  !
  nbndfs = 0
  DO ik = 1, nkfs ! loop over irreducible k-points
     READ(iufil) wkfs(ik), xkfs(1, ik), xkfs(2, ik), xkfs(3, ik)
     !WRITE(*,'(4f15.7)') xkfs(1, ik), xkfs(2, ik), xkfs(3, ik), wkfs(ik)
     DO ibnd = 1, nbnd_
        READ(iufil) ekf_(ibnd, ik)
     !WRITE(*,'(f15.7)') ekf_(ibnd, ik)
     ENDDO
     n = 0
     DO ibnd = 1, nbnd_
        ! go from Ryd to eV
        ekf_(ibnd, ik) = ekf_(ibnd, ik) * ryd2ev
        IF ( abs( ekf_(ibnd, ik) - ef0 ) .lt. fsthick ) THEN
           n = n + 1
           IF ( nbndfs .lt. n ) nbndfs = n
        ENDIF
     ENDDO
  ENDDO
  !WRITE(*,'(2i7)') nbnd_, nbndfs
  CLOSE(iufil)
  !
  IF (.not. ALLOCATED(ekfs)) ALLOCATE(ekfs(nbndfs, nkfs))
  IF (.not. ALLOCATED(w0g))  ALLOCATE(w0g(nbndfs, nkfs))
  !
  ekfs(:, :) = ef0 - 10.d0 * fsthick
  w0g(:, :)  = 0.d0
  !
  DO ik = 1, nkfs ! loop over k-points
     n = 0
     DO ibnd = 1, nbnd_
        IF ( abs( ekf_(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
           n = n + 1
           ekfs(n,ik) = ekf_(ibnd,ik)
           arg = ( ekfs(n,ik) - ef0 ) / degaussw
           w0gauss = exp( -arg*arg ) / sqrt(3.14159265358979323846d0)
           w0g(n,ik) =  w0gauss / degaussw
        ENDIF
     ENDDO
  ENDDO

  IF ( ALLOCATED(ekf_) ) DEALLOCATE(ekf_)
  IF ( ALLOCATED(xkfs) ) DEALLOCATE(xkfs)
  !
  IF ( .not. ALLOCATED(delta) ) ALLOCATE(delta(nbndfs,nkfs))
  !
  delta(:, :) = 0.d0
  iufil = 10
  OPEN(unit = iufil, file = "NbSe2.imag_aniso_004.00", form = 'formatted')
  READ(iufil,'(a)') word
  DO ik = 1, nkfs
     DO ibnd = 1, nbndfs
        IF ( abs( ekfs(ibnd,ik) - ef0 ) .lt. fsthick ) THEN
           READ(iufil,'(5ES20.10)') w, e, z, delta(ibnd,ik), zn
           !WRITE(*,'(5ES20.10)') w, e, z, delta(ibnd,ik), zn
        ENDIF
     ENDDO ! ibnd
  ENDDO ! ik
  CLOSE(iufil)
  !
  ! Adapted from io_eliashberg.f90/SUBROUTINE gap_distribution_FS(itemp, cname)
  delta_min = MINVAL(delta(:, :))
  delta_max = MAXVAL(delta(:, :))
  !
  WRITE(*, '(a, 2f12.6, a)') 'Min. / Max. values of superconducting gap = ', &
                delta_min * 1000.d0, delta_max * 1000.d0, ' meV'
  !
  IF (delta_min > zero) THEN
    delta_min = 0.9d0 * delta_min
  ELSE
    delta_min = 1.1d0 * delta_min
  ENDIF
  !
  IF (delta_max > zero) THEN
    delta_max = 1.1d0 * delta_max
  ELSE
    delta_max = 0.9d0 * delta_max
  ENDIF
  !
  !nbin = NINT((delta_max - delta_min) / eps4) + 1
  !dbin = (delta_max - delta_min) / DBLE(nbin)
  dbin = 3.0d-5 !eV
  nbin = NINT((delta_max - delta_min) / dbin) + 1
  WRITE(*,'(3d24.15)') maxval(delta(:, :)), dbin, dosef
  WRITE(*, *) "nbin, dbin", nbin, dbin
  IF ( .not. ALLOCATED(delta_k_bin) ) ALLOCATE ( delta_k_bin(nbin+1) )
  delta_k_bin(:) = 0.d0
  !
  !
  !TOTAL
  DO ik = 1, nkfs
     DO ibnd = 1, nbndfs
        IF ((ABS( ekfs(ibnd,ik) - ef0 ) .lt. fsthick) ) THEN
           ibin = NINT((delta(ibnd,ik) - delta_min) / dbin) + 1
           weight = w0g(ibnd, ik)
           delta_k_bin(ibin) = delta_k_bin(ibin) + weight
           !WRITE(*,'(3i7,4d24.15)') ik, ibnd, ibin, wkfs(ik), w0g(ibnd, ik), weight, delta_k_bin(ibin)
        ENDIF
     ENDDO
  ENDDO
  !
  OPEN( unit = iufil, file = "new_gap.dat", form = 'formatted')
  DO ibin = 1, nbin
     WRITE(iufil,'(2d24.15)') delta_k_bin(ibin) / maxval(delta_k_bin(:)), dbin * dble(ibin)
  ENDDO
  CLOSE(iufil)
  !
  DEALLOCATE(delta)
  DEALLOCATE(ekfs)
  DEALLOCATE(wkfs)
  DEALLOCATE(w0g)
  DEALLOCATE(delta_k_bin)
  !
  STOP
  !
  END PROGRAM main