tran_lcr_2c2_build_ham Subroutine

private subroutine tran_lcr_2c2_build_ham(pl_warning)

Uses

  • proc~~tran_lcr_2c2_build_ham~~UsesGraph proc~tran_lcr_2c2_build_ham tran_lcr_2c2_build_ham module~w90_constants w90_constants proc~tran_lcr_2c2_build_ham->module~w90_constants module~w90_io w90_io proc~tran_lcr_2c2_build_ham->module~w90_io module~w90_parameters w90_parameters proc~tran_lcr_2c2_build_ham->module~w90_parameters module~w90_hamiltonian w90_hamiltonian proc~tran_lcr_2c2_build_ham->module~w90_hamiltonian module~w90_io->module~w90_constants module~w90_parameters->module~w90_constants module~w90_parameters->module~w90_io module~w90_hamiltonian->module~w90_constants module~w90_comms w90_comms module~w90_hamiltonian->module~w90_comms module~w90_comms->module~w90_constants module~w90_comms->module~w90_io

Arguments

Type IntentOptional AttributesName
logical, intent(in) :: pl_warning

Calls

proc~~tran_lcr_2c2_build_ham~~CallsGraph proc~tran_lcr_2c2_build_ham tran_lcr_2c2_build_ham proc~io_error io_error proc~tran_lcr_2c2_build_ham->proc~io_error proc~tran_reduce_hr tran_reduce_hr proc~tran_lcr_2c2_build_ham->proc~tran_reduce_hr proc~tran_cut_hr_one_dim tran_cut_hr_one_dim proc~tran_lcr_2c2_build_ham->proc~tran_cut_hr_one_dim proc~io_date io_date proc~tran_lcr_2c2_build_ham->proc~io_date proc~io_file_unit io_file_unit proc~tran_lcr_2c2_build_ham->proc~io_file_unit proc~tran_reduce_hr->proc~io_error

Called by

proc~~tran_lcr_2c2_build_ham~~CalledByGraph proc~tran_lcr_2c2_build_ham tran_lcr_2c2_build_ham proc~tran_main tran_main proc~tran_main->proc~tran_lcr_2c2_build_ham program~wannier wannier program~wannier->proc~tran_main proc~wannier_run wannier_run proc~wannier_run->proc~tran_main

Contents


Source Code

  subroutine tran_lcr_2c2_build_ham(pl_warning)
    !==============================================!
    ! Builds hamiltonians blocks required for the  !
    ! Greens function caclulations of the quantum  !
    ! conductance according to the 2c2 geometry.   !
    ! Leads are also symmetrised, in that unit cell!
    ! sub-blocks are copied to create truely ideal !
    ! leads.                                       !
    !==============================================!

    use w90_constants, only: dp, eps5
    use w90_io, only: io_error, stdout, seedname, io_file_unit, io_date, io_stopwatch
    use w90_parameters, only: tran_num_cell_ll, num_wann, tran_num_ll, kpt_cart, nfermi, fermi_energy_list, &
      tran_write_ht, tran_num_rr, tran_num_lc, tran_num_cr, tran_num_cc, &
      tran_num_bandc, timing_level, dist_cutoff_mode, dist_cutoff, &
      dist_cutoff_hc
    use w90_hamiltonian, only: wannier_centres_translated

    implicit none

    logical, intent(in)                     :: pl_warning

    integer                                :: i, j, k, num_wann_cell_ll, file_unit, ierr, band_size

    real(dp), allocatable, dimension(:, :)    :: sub_block
    real(dp)                               :: PL_length, dist, dist_vec(3)

    character(len=9)                       :: cdate, ctime

    if (timing_level > 1) call io_stopwatch('tran: lcr_2c2_build_ham', 1)

    if (nfermi > 1) call io_error("Error in tran_lcr_2c2_build_ham: nfermi>1. " &
                                  //"Set the fermi level using the input parameter 'fermi_evel'")

    allocate (hL0(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hL0 in tran_lcr_2c2_build_ham')
    allocate (hL1(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hL1 in tran_lcr_2c2_build_ham')
    allocate (hR0(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hR0 in tran_lcr_2c2_build_ham')
    allocate (hR1(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hR1 in tran_lcr_2c2_build_ham')
    allocate (hLC(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hLC in tran_lcr_2c2_build_ham')
    allocate (hCR(tran_num_ll, tran_num_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hCR in tran_lcr_2c2_build_ham')
    allocate (hC(num_wann - (2*tran_num_ll), num_wann - (2*tran_num_ll)), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating hC in tran_lcr_2c2_build_ham')
    !
    !This checks that only the gamma point is used in wannierisation
    !This is necessary since this calculation only makes sense if we
    !have periodicity over the supercell.
    !
    if ((size(kpt_cart, 2) .ne. 1) .and. (kpt_cart(1, 1) .eq. 0.0_dp) &
        .and. (kpt_cart(2, 1) .eq. 0.0_dp) &
        .and. (kpt_cart(3, 1) .eq. 0.0_dp)) then
      call io_error('Calculation must be performed at gamma only')
    endif

    num_wann_cell_ll = tran_num_ll/tran_num_cell_ll

    allocate (sub_block(num_wann_cell_ll, num_wann_cell_ll), stat=ierr)
    if (ierr /= 0) call io_error('Error in allocating sub_block in tran_lcr_2c2_build_ham')
    !
    !Build hL0 & hL1
    !
    hL0 = 0.0_dp
    hL1 = 0.0_dp
    !
    !Loop over the sub_blocks corresponding to distinct unit cells inside the principal layer
    !
    do i = 1, tran_num_cell_ll
      !
      !Each sub_block will be duplicated along the corresponding diagonal. This ensures the correct symmetry for the leads.
      !
      sub_block = 0.0_dp
      !
      !Extract matrix elements from hr_one_dim needed for hL0 (and lower triangular sub_blocks of hL1)
      !
      do j = 1, num_wann_cell_ll
        do k = 1, num_wann_cell_ll
          sub_block(j, k) = hr_one_dim(tran_sorted_idx(j), tran_sorted_idx((i - 1)*num_wann_cell_ll + k), 0)
        enddo
      enddo
      !
      !Filling up hL0 sub_block by sub_block
      !
      do j = 1, tran_num_cell_ll - i + 1
        !
        !Fill diagonal and upper diagonal sub_blocks
        !
        hL0((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
            (j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll) = sub_block
        !
        !Fill lower diagonal sub_blocks
        !
        if (i .gt. 1) then
          hL0((j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll, &
              (j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = transpose(sub_block)
        endif
      enddo
      !
      !Filling up non-diagonal hL1 sub_blocks (nothing need be done for i=1)
      !
      if (i .gt. 1) then
        do j = 1, i - 1
          hL1((tran_num_cell_ll - (i - j))*num_wann_cell_ll + 1:(tran_num_cell_ll - (i - 1 - j))*num_wann_cell_ll, &
              (j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = sub_block
        enddo
      endif
      !
      ! MS: Get diagonal and upper triangular sublocks for hL1 - use periodic image of PL4
      !
      sub_block = 0.0_dp
      !
      if (i == 1) then !Do diagonal only
        do j = 1, num_wann_cell_ll
          do k = 1, num_wann_cell_ll
            sub_block(j, k) = hr_one_dim( &
                              tran_sorted_idx(num_wann - tran_num_ll + j), &
                              tran_sorted_idx((i - 1)*num_wann_cell_ll + k), 0)
          enddo
        enddo
        !
        ! MS: Now fill subblocks of hL1
        !
        do j = 1, tran_num_cell_ll - i + 1
          hL1((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
              (j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)* &
              num_wann_cell_ll) = sub_block
        enddo
      endif
    enddo
    !
    !Special case tran_num_cell_ll=1, the diagonal sub-block of hL1 is hL1, so cannot be left as zero
    !
    if (tran_num_cell_ll .eq. 1) then
      do j = num_wann - num_wann_cell_ll + 1, num_wann
        do k = 1, num_wann_cell_ll
          hL1(j - num_wann + num_wann_cell_ll, k) = hr_one_dim(tran_sorted_idx(j), tran_sorted_idx(k), 0)
        enddo
      enddo
    endif

    !
    !Build hR0 & hR1
    !
    hR0 = 0.0_dp
    hR1 = 0.0_dp
    !
    !Loop over the sub_blocks corresponding to distinct unit cells inside the principal layer
    !
    do i = 1, tran_num_cell_ll
      !
      !Each sub_block will be duplicated along the corresponding diagonal. This ensures the correct symmetry for the leads.
      !
      sub_block = 0.0_dp
      !
      !Extract matrix elements from hr_one_dim needed for hR0 (and lower triangular sub_blocks of hR1)
      !
      do j = 1, num_wann_cell_ll
        do k = 1, num_wann_cell_ll
          sub_block(j, k) = hr_one_dim(tran_sorted_idx(num_wann - i*(num_wann_cell_ll) + j), &
                                       tran_sorted_idx(num_wann - num_wann_cell_ll + k), 0)
        enddo
      enddo
      !
      !Filling up hR0 sub_block by sub_block
      !
      do j = 1, tran_num_cell_ll - i + 1
        !
        !Fill diagonal and upper diagonal sub_blocks
        !
        hR0((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
            (j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll) = sub_block
        !
        !Fill lower diagonal sub_blocks
        !
        if (i .gt. 1) then
          hR0((j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll, &
              (j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = transpose(sub_block)
        endif
      enddo
      !
      !Filling up non-diagonal hR1 sub_blocks (nothing need be done for i=1)
      !
      if (i .gt. 1) then
        do j = 1, i - 1
          hR1((tran_num_cell_ll - (i - j))*num_wann_cell_ll + 1:(tran_num_cell_ll - (i - 1 - j))*num_wann_cell_ll, &
              (j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll) = sub_block
        enddo
      endif
      !
      ! MS: Get diagonal and upper triangular sublocks for hR1 - use periodic image of PL1
      !
      sub_block = 0.0_dp
      !
      if (i == 1) then  !Do diagonal only
        do j = 1, num_wann_cell_ll
          do k = 1, num_wann_cell_ll
            sub_block(j, k) = hr_one_dim(tran_sorted_idx((i - 1)*num_wann_cell_ll + k), &
                                         tran_sorted_idx(num_wann - tran_num_ll + j), 0)
          enddo
        enddo
        !
        ! MS: Now fill subblocks of hR1
        !
        do j = 1, tran_num_cell_ll - i + 1
          hR1((j - 1)*num_wann_cell_ll + 1:j*num_wann_cell_ll, &
              (j - 1)*num_wann_cell_ll + 1 + (i - 1)*num_wann_cell_ll:j*num_wann_cell_ll + (i - 1)*num_wann_cell_ll) = sub_block
        enddo
      endif
    enddo
    !
    !Special case tran_num_cell_ll=1, the diagonal sub-block of hR1 is hR1, so cannot be left as zero
    !
    if (tran_num_cell_ll .eq. 1) then
      do j = 1, num_wann_cell_ll
        do k = num_wann - num_wann_cell_ll + 1, num_wann
          hR1(k - num_wann + num_wann_cell_ll, j) = hr_one_dim(tran_sorted_idx(j), tran_sorted_idx(k), 0)
        enddo
      enddo
    endif

    !
    !Building hLC
    !
    hLC = 0.0_dp
    do i = 1, tran_num_ll
      do j = tran_num_ll + 1, 2*tran_num_ll
        hLC(i, j - tran_num_ll) = hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
      enddo
    enddo
!----!
! MS ! Rely on dist_cutoff doing the work here, as it cuts element-wise, not block wise (incorrect)
!----!
!    if (tran_num_cell_ll .gt. 1) then
!        do j=1,tran_num_cell_ll
!            do k=1,tran_num_cell_ll
!                if (k .ge. j) then
!                    hLC((j-1)*num_wann_cell_ll+1:j*num_wann_cell_ll,(k-1)*num_wann_cell_ll+1:k*num_wann_cell_ll)=0.0_dp
!                endif
!            enddo
!        enddo
!    endif
!---!
!end!
!---!

    !
    !Building hC
    !
    hC = 0.0_dp
    !
    band_size = 0
    if (dist_cutoff_hc .ne. dist_cutoff) then
      dist_cutoff = dist_cutoff_hc
      write (stdout, *) 'Applying dist_cutoff_hc to Hamiltonian for construction of hC'
      deallocate (hr_one_dim, stat=ierr)
      if (ierr /= 0) call io_error('Error deallocating hr_one_dim in tran_lcr_2c2_sort')
      call tran_reduce_hr()
      call tran_cut_hr_one_dim()
    endif

    do i = tran_num_ll + 1, num_wann - tran_num_ll
      do j = tran_num_ll + 1, num_wann - tran_num_ll
        hC(i - tran_num_ll, j - tran_num_ll) = hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
        !
        ! Impose a ham_cutoff of 1e-4 eV to reduce tran_num_bandc (and in turn hCband, and speed up transport)
        !
        if (abs(hC(i - tran_num_ll, j - tran_num_ll)) .lt. 10.0_dp*eps5) then
          hC(i - tran_num_ll, j - tran_num_ll) = 0.0_dp
          band_size = max(band_size, abs(i - j))
        endif
      enddo
    enddo
    !
    !Building hCR
    !
    hCR = 0.0_dp
    do i = num_wann - 2*tran_num_ll + 1, num_wann - tran_num_ll
      do j = num_wann - tran_num_ll + 1, num_wann
        hCR(i - (num_wann - 2*tran_num_ll), j - (num_wann - tran_num_ll)) = hr_one_dim(tran_sorted_idx(i), tran_sorted_idx(j), 0)
      enddo
    enddo
!----!
! MS ! Rely on dist_cutoff doing the work here, as it cuts element-wise, not block wise (incorrect)
!----!
!    if (tran_num_cell_ll .gt. 1) then
!        do j=1,tran_num_cell_ll
!            do k=1,tran_num_cell_ll
!                if (k .ge. j) then
!                    hCR((j-1)*num_wann_cell_ll+1:j*num_wann_cell_ll,(k-1)*num_wann_cell_ll+1:k*num_wann_cell_ll)=0.0_dp
!                endif
!            enddo
!        enddo
!    endif
!---!
!end!
!---!

    !
    !Subtract the Fermi energy from the diagonal elements of hC,hL0,hR0
    !
    do i = 1, tran_num_ll
      hL0(i, i) = hL0(i, i) - fermi_energy_list(1)
      hR0(i, i) = hR0(i, i) - fermi_energy_list(1)
    enddo
    do i = 1, num_wann - (2*tran_num_ll)
      hC(i, i) = hC(i, i) - fermi_energy_list(1)
    enddo
    !
    !Define tran_num_** parameters that are used later in tran_lcr
    !
    tran_num_rr = tran_num_ll
    tran_num_lc = tran_num_ll
    tran_num_cr = tran_num_ll
    tran_num_cc = num_wann - (2*tran_num_ll)

    !
    ! Set appropriate tran_num_bandc if has not been set (0.0_dp is default value)
    !
    if (tran_num_bandc .eq. 0.0_dp) then
      tran_num_bandc = min(band_size + 1, (tran_num_cc + 1)/2 + 1)
    endif

    !
    ! MS: Find and print effective PL length
    !
    if (.not. pl_warning) then
      PL_length = 0.0_dp
      do i = 1, tran_num_ll
        do j = 1, tran_num_ll
          if (abs(hL1(i, j)) .gt. 0.0_dp) then
            if (index(dist_cutoff_mode, 'one_dim') .gt. 0) then
              dist = abs(wannier_centres_translated(coord(1), tran_sorted_idx(i)) &
                         - wannier_centres_translated(coord(1), tran_sorted_idx(j + tran_num_ll)))
            else
              dist_vec(:) = wannier_centres_translated(:, tran_sorted_idx(i)) &
                            - wannier_centres_translated(:, tran_sorted_idx(j + tran_num_ll))
              dist = sqrt(dot_product(dist_vec, dist_vec))
            endif
            PL_length = max(PL_length, dist)
          endif
          if (abs(hR1(i, j)) .gt. 0.0_dp) then
            if (index(dist_cutoff_mode, 'one_dim') .gt. 0) then
              dist = abs(wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - 2*tran_num_ll + i)) &
                         - wannier_centres_translated(coord(1), tran_sorted_idx(num_wann - tran_num_ll + j)))
            else
              dist_vec(:) = wannier_centres_translated(:, tran_sorted_idx(num_wann - 2*tran_num_ll + i)) &
                            - wannier_centres_translated(:, tran_sorted_idx(num_wann - tran_num_ll + j))
              dist = sqrt(dot_product(dist_vec, dist_vec))
            endif
            PL_length = max(PL_length, dist)
          endif
        enddo
      enddo
      write (stdout, '(1x,a,f12.6,a)') 'Approximate effective principal layer length is: ', PL_length, ' Ang.'
    endif

    !
    !Writing to file:
    !
    if (tran_write_ht) then
      write (stdout, *) '------------------------------- Writing ht files  ----------------------------'
      !
      file_unit = io_file_unit()
      open (file_unit, file=trim(seedname)//'_htL.dat', status='unknown', form='formatted', action='write')

      call io_date(cdate, ctime)
      write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
      write (file_unit, '(I6)') tran_num_ll
      write (file_unit, '(6F12.6)') ((hL0(j, i), j=1, tran_num_ll), i=1, tran_num_ll)
      write (file_unit, '(I6)') tran_num_ll
      write (file_unit, '(6F12.6)') ((hL1(j, i), j=1, tran_num_ll), i=1, tran_num_ll)

      close (file_unit)
      write (stdout, *) ' '//trim(seedname)//'_htL.dat  written'
      !
      !hR
      !
      file_unit = io_file_unit()
      open (file_unit, file=trim(seedname)//'_htR.dat', status='unknown', form='formatted', action='write')

      call io_date(cdate, ctime)
      write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
      write (file_unit, '(I6)') tran_num_rr
      write (file_unit, '(6F12.6)') ((hR0(j, i), j=1, tran_num_rr), i=1, tran_num_rr)
      write (file_unit, '(I6)') tran_num_rr
      write (file_unit, '(6F12.6)') ((hR1(j, i), j=1, tran_num_rr), i=1, tran_num_rr)

      close (file_unit)
      write (stdout, *) ' '//trim(seedname)//'_htR.dat  written'
      !
      !hLC
      !
      file_unit = io_file_unit()
      open (file_unit, file=trim(seedname)//'_htLC.dat', status='unknown', form='formatted', action='write')

      call io_date(cdate, ctime)
      write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
      write (file_unit, '(2I6)') tran_num_ll, tran_num_lc
      write (file_unit, '(6F12.6)') ((hLC(j, i), j=1, tran_num_lc), i=1, tran_num_lc)

      close (file_unit)
      write (stdout, *) ' '//trim(seedname)//'_htLC.dat written'
      !
      !hCR
      !
      file_unit = io_file_unit()
      open (file_unit, file=trim(seedname)//'_htCR.dat', status='unknown', form='formatted', action='write')

      call io_date(cdate, ctime)
      write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
      write (file_unit, '(2I6)') tran_num_cr, tran_num_rr
      write (file_unit, '(6F12.6)') ((hCR(j, i), j=1, tran_num_cr), i=1, tran_num_cr)

      close (file_unit)
      write (stdout, *) ' '//trim(seedname)//'_htCR.dat written'
      !
      !hC
      !
      file_unit = io_file_unit()
      open (file_unit, file=trim(seedname)//'_htC.dat', status='unknown', form='formatted', action='write')

      call io_date(cdate, ctime)
      write (file_unit, *) 'written on '//cdate//' at '//ctime ! Date and time
      write (file_unit, '(I6)') tran_num_cc
      write (file_unit, '(6F12.6)') ((hC(j, i), j=1, tran_num_cc), i=1, tran_num_cc)

      close (file_unit)
      write (stdout, *) ' '//trim(seedname)//'_htC.dat  written'

      write (stdout, *) '------------------------------------------------------------------------------'
    end if

    deallocate (sub_block, stat=ierr)
    if (ierr /= 0) call io_error('Error deallocating sub_block in tran_lcr_2c2_build_ham')

    if (timing_level > 1) call io_stopwatch('tran: lcr_2c2_build_ham', 2)

    return

  end subroutine tran_lcr_2c2_build_ham