gyrotropic_main Subroutine

public subroutine gyrotropic_main()

Uses

  • proc~~gyrotropic_main~~UsesGraph proc~gyrotropic_main gyrotropic_main module~w90_postw90_common w90_postw90_common proc~gyrotropic_main->module~w90_postw90_common module~w90_utility w90_utility proc~gyrotropic_main->module~w90_utility module~w90_constants w90_constants proc~gyrotropic_main->module~w90_constants module~w90_get_oper w90_get_oper proc~gyrotropic_main->module~w90_get_oper module~w90_io w90_io proc~gyrotropic_main->module~w90_io module~w90_comms w90_comms proc~gyrotropic_main->module~w90_comms module~w90_parameters w90_parameters proc~gyrotropic_main->module~w90_parameters module~w90_postw90_common->module~w90_constants module~w90_postw90_common->module~w90_comms module~w90_utility->module~w90_constants module~w90_get_oper->module~w90_constants module~w90_io->module~w90_constants module~w90_comms->module~w90_constants module~w90_comms->module~w90_io module~w90_parameters->module~w90_constants module~w90_parameters->module~w90_io

Computes the following quantities: (i) D tensor (ii) K tensor (iii) C tensor (iv) current-induced optical activity (v) natural optical activity

Arguments

None

Calls

proc~~gyrotropic_main~~CallsGraph proc~gyrotropic_main gyrotropic_main proc~get_hh_r get_HH_R proc~gyrotropic_main->proc~get_hh_r proc~get_bb_r get_BB_R proc~gyrotropic_main->proc~get_bb_r proc~get_aa_r get_AA_R proc~gyrotropic_main->proc~get_aa_r proc~utility_det3 utility_det3 proc~gyrotropic_main->proc~utility_det3 proc~get_ss_r get_SS_R proc~gyrotropic_main->proc~get_ss_r proc~gyrotropic_get_k_list gyrotropic_get_k_list proc~gyrotropic_main->proc~gyrotropic_get_k_list proc~get_cc_r get_CC_R proc~gyrotropic_main->proc~get_cc_r proc~io_stopwatch io_stopwatch proc~gyrotropic_main->proc~io_stopwatch proc~gyrotropic_outprint_tensor gyrotropic_outprint_tensor proc~gyrotropic_main->proc~gyrotropic_outprint_tensor proc~fourier_q_to_r fourier_q_to_R proc~get_hh_r->proc~fourier_q_to_r proc~io_error io_error proc~get_hh_r->proc~io_error proc~get_win_min get_win_min proc~get_hh_r->proc~get_win_min proc~io_file_unit io_file_unit proc~get_hh_r->proc~io_file_unit proc~get_bb_r->proc~io_error proc~get_bb_r->proc~get_win_min proc~get_bb_r->proc~io_file_unit proc~get_aa_r->proc~io_error proc~get_aa_r->proc~io_file_unit proc~get_ss_r->proc~io_file_unit proc~pw90common_fourier_r_to_k_vec pw90common_fourier_R_to_k_vec proc~gyrotropic_get_k_list->proc~pw90common_fourier_r_to_k_vec proc~spin_get_s spin_get_S proc~gyrotropic_get_k_list->proc~spin_get_s proc~gyrotropic_get_noa_k gyrotropic_get_NOA_k proc~gyrotropic_get_k_list->proc~gyrotropic_get_noa_k proc~wham_get_eig_deleig wham_get_eig_deleig proc~gyrotropic_get_k_list->proc~wham_get_eig_deleig proc~wham_get_d_h wham_get_D_h proc~gyrotropic_get_k_list->proc~wham_get_d_h proc~utility_w0gauss utility_w0gauss proc~gyrotropic_get_k_list->proc~utility_w0gauss proc~utility_rotate utility_rotate proc~gyrotropic_get_k_list->proc~utility_rotate proc~berry_get_imf_klist berry_get_imf_klist proc~gyrotropic_get_k_list->proc~berry_get_imf_klist proc~berry_get_imfgh_klist berry_get_imfgh_klist proc~gyrotropic_get_k_list->proc~berry_get_imfgh_klist proc~get_cc_r->proc~get_win_min proc~get_cc_r->proc~io_file_unit proc~io_stopwatch->proc~io_error proc~gyrotropic_outprint_tensor_w gyrotropic_outprint_tensor_w proc~gyrotropic_outprint_tensor->proc~gyrotropic_outprint_tensor_w proc~gyrotropic_outprint_tensor->proc~io_file_unit proc~pw90common_fourier_r_to_k pw90common_fourier_R_to_k proc~spin_get_s->proc~pw90common_fourier_r_to_k proc~utility_rotate_diag utility_rotate_diag proc~spin_get_s->proc~utility_rotate_diag proc~utility_diagonalize utility_diagonalize proc~spin_get_s->proc~utility_diagonalize proc~gyrotropic_get_noa_k->proc~utility_rotate proc~gyrotropic_get_noa_bnl_spin gyrotropic_get_NOA_Bnl_spin proc~gyrotropic_get_noa_k->proc~gyrotropic_get_noa_bnl_spin proc~gyrotropic_get_noa_bnl_orb gyrotropic_get_NOA_Bnl_orb proc~gyrotropic_get_noa_k->proc~gyrotropic_get_noa_bnl_orb proc~wham_get_eig_deleig->proc~get_hh_r proc~wham_get_eig_deleig->proc~pw90common_fourier_r_to_k proc~wham_get_eig_deleig->proc~utility_diagonalize proc~wham_get_d_h->proc~utility_rotate proc~berry_get_imf_klist->proc~berry_get_imfgh_klist proc~berry_get_imfgh_klist->proc~pw90common_fourier_r_to_k_vec proc~wham_get_eig_uu_hh_jjlist wham_get_eig_UU_HH_JJlist proc~berry_get_imfgh_klist->proc~wham_get_eig_uu_hh_jjlist proc~utility_im_tr_prod utility_im_tr_prod proc~berry_get_imfgh_klist->proc~utility_im_tr_prod proc~wham_get_occ_mat_list wham_get_occ_mat_list proc~berry_get_imfgh_klist->proc~wham_get_occ_mat_list proc~utility_re_tr_prod utility_re_tr_prod proc~berry_get_imfgh_klist->proc~utility_re_tr_prod proc~wham_get_eig_uu_hh_jjlist->proc~get_hh_r proc~wham_get_eig_uu_hh_jjlist->proc~utility_diagonalize proc~utility_zgemm_new utility_zgemm_new proc~utility_rotate_diag->proc~utility_zgemm_new proc~utility_matmul_diag utility_matmul_diag proc~utility_rotate_diag->proc~utility_matmul_diag proc~wham_get_occ_mat_list->proc~io_error proc~utility_diagonalize->proc~io_error

Contents

Source Code


Source Code

  subroutine gyrotropic_main
    !============================================================!
    !                                                            !
    !! Computes the following quantities:
    !!   (i) D tensor
    !!  (ii) K tensor
    !! (iii) C tensor
    !!  (iv) current-induced optical activity
    !!   (v) natural optical activity
    !                                                            !
    !============================================================!

    use w90_constants, only: dp, cmplx_0, elem_charge_SI, hbar_SI, &
      eV_au, bohr, elec_mass_SI, twopi, eps0_SI
    use w90_comms, only: on_root, num_nodes, my_node_id, comms_reduce
    use w90_utility, only: utility_det3
    use w90_io, only: io_error, stdout, io_file_unit, seedname, &
      io_stopwatch
    use w90_postw90_common, only: nrpts, irvec, num_int_kpts_on_node, int_kpts, &
      weight
    use w90_parameters, only: timing_level, iprint, num_wann, gyrotropic_kmesh, &
      cell_volume, transl_inv, gyrotropic_task, &
      gyrotropic_nfreq, gyrotropic_freq_list, nfermi, &
      fermi_energy_list, gyrotropic_box, gyrotropic_box_corner, spinors
    use w90_get_oper, only: get_HH_R, get_AA_R, get_BB_R, get_CC_R, &
      get_SS_R

    real(kind=dp), allocatable    :: gyro_K_spn(:, :, :)
    real(kind=dp), allocatable    :: gyro_DOS(:)
    real(kind=dp), allocatable    :: gyro_K_orb(:, :, :)
    real(kind=dp), allocatable    :: gyro_C(:, :, :)
    real(kind=dp), allocatable    :: gyro_D(:, :, :)
    real(kind=dp), allocatable    :: gyro_Dw(:, :, :, :)
    real(kind=dp), allocatable    :: gyro_NOA_spn(:, :, :, :)
    real(kind=dp), allocatable    :: gyro_NOA_orb(:, :, :, :)

    character(len=30) :: f_out_name_tmp
    character(len=30) :: units_tmp
    character(len=120) :: comment_tmp

    real(kind=dp)     :: kweight, kpt(3), &
                         db1, db2, db3, fac, freq
    integer           :: n, i, j, k, ikpt, if, ierr, loop_x, loop_y, loop_z, &
                         loop_xyz, ifreq, &
                         file_unit
    logical           :: eval_K, eval_C, eval_D, eval_Dw, eval_NOA, eval_spn, eval_DOS

    if (nfermi == 0) call io_error( &
      'Must specify one or more Fermi levels when gyrotropic=true')

    if (timing_level > 1 .and. on_root) call io_stopwatch('gyrotropic: prelims', 1)

    ! Mesh spacing in reduced coordinates
    !
    db1 = 1.0_dp/real(gyrotropic_kmesh(1), dp)
    db2 = 1.0_dp/real(gyrotropic_kmesh(2), dp)
    db3 = 1.0_dp/real(gyrotropic_kmesh(3), dp)

    eval_K = .false.
    eval_C = .false.
    eval_D = .false.
    eval_Dw = .false.
    eval_spn = .false.
    eval_NOA = .false.
    eval_DOS = .false.

    if (index(gyrotropic_task, '-k') > 0) eval_K = .true.
    if (index(gyrotropic_task, '-c') > 0) eval_C = .true.
    if (index(gyrotropic_task, '-d0') > 0) eval_D = .true.
    if (index(gyrotropic_task, '-dw') > 0) eval_Dw = .true.
    if (index(gyrotropic_task, '-spin') > 0) eval_spn = .true.
    if (index(gyrotropic_task, '-noa') > 0) eval_NOA = .true.
    if (index(gyrotropic_task, '-dos') > 0) eval_DOS = .true.
    if (index(gyrotropic_task, 'all') > 0) then
      eval_K = .true.
      eval_C = .true.
      eval_D = .true.
      eval_Dw = .true.
      if (spinors) eval_spn = .true.
      eval_NOA = .true.
      eval_DOS = .true.
    endif

    if (.not. (eval_K .or. eval_noa)) eval_spn = .false.

    if ((.not. spinors) .and. eval_spn) call io_error( &
      "spin contribution requested for gyrotropic, but the wavefunctions are not spinors")

    ! Wannier matrix elements, allocations and initializations

    call get_HH_R
    if (eval_D .or. eval_Dw .or. eval_K .or. eval_NOA) then
      call get_AA_R
    endif

    if (eval_spn) then
      call get_SS_R
    endif

    if (eval_K) then
      call get_BB_R
      call get_CC_R
      allocate (gyro_K_orb(3, 3, nfermi))
      gyro_K_orb = 0.0_dp
      if (eval_spn) then
        allocate (gyro_K_spn(3, 3, nfermi))
        gyro_K_spn = 0.0_dp
      endif
    endif

    if (eval_D) then
      allocate (gyro_D(3, 3, nfermi))
      gyro_D = 0.0_dp
    endif

    if (eval_DOS) then
      allocate (gyro_DOS(nfermi))
      gyro_DOS = 0.0_dp
    endif

    if (eval_C) then
      allocate (gyro_C(3, 3, nfermi))
      gyro_C = 0.0_dp
    endif

    if (eval_Dw) then
      allocate (gyro_Dw(3, 3, nfermi, gyrotropic_nfreq))
      gyro_Dw = 0.0_dp
    endif

    if (eval_NOA) then
      allocate (gyro_NOA_orb(3, 3, nfermi, gyrotropic_nfreq))
      gyro_NOA_orb = 0.0_dp
      if (eval_spn) then
        allocate (gyro_NOA_spn(3, 3, nfermi, gyrotropic_nfreq))
        gyro_NOA_spn = 0.0_dp
      endif
    endif

    if (on_root) then
      flush(stdout)
      write (stdout, '(/,/,1x,a)') 'Properties calculated in module  g y r o t r o p i c'
      write (stdout, '(1x,a)') '------------------------------------------'

      if (eval_D) write (stdout, '(/,3x,a)') '* D-tensor  --- Eq.2 of TAS17 '

      if (eval_dos) write (stdout, '(/,3x,a)') '* density of states '

      if (eval_K) then
        write (stdout, '(/,3x,a)') '* K-tensor  --- Eq.3 of TAS17 '
        if (eval_spn) then
          write (stdout, '(3x,a)') '    * including spin component '
        else
          write (stdout, '(3x,a)') '    * excluding spin component '
        endif
      endif

      if (eval_Dw) write (stdout, '(/,3x,a)') '* Dw-tensor  --- Eq.12 of TAS17 '

      if (eval_C) write (stdout, '(/,3x,a)') '* C-tensor  --- Eq.B6 of TAS17 '

      if (eval_NOA) then
        write (stdout, '(/,3x,a)') '* gamma-tensor of NOA --- Eq.C12 of TAS17 '
        if (eval_spn) then
          write (stdout, '(3x,a)') '    * including spin component '
        else
          write (stdout, '(3x,a)') '    * excluding spin component '
        endif
      endif

      if (transl_inv) then
        if (eval_K) &
          call io_error('transl_inv=T disabled for K-tensor')
        write (stdout, '(/,1x,a)') &
          'Using a translationally-invariant discretization for the'
        write (stdout, '(1x,a)') &
          'band-diagonal Wannier matrix elements of r, etc.'
      endif

      if (timing_level > 1) then
        call io_stopwatch('gyrotropic: prelims', 2)
        call io_stopwatch('gyrotropic: k-interpolation', 1)
      endif

      write (stdout, '(1x,a20,3(i0,1x))') 'Interpolation grid: ', gyrotropic_kmesh(1:3)

      flush(stdout)

    end if !on_root

    ! Do not read 'kpoint.dat'. Loop over a regular grid in the full BZ

    kweight = db1*db2*db3*utility_det3(gyrotropic_box)

    do loop_xyz = my_node_id, PRODUCT(gyrotropic_kmesh) - 1, num_nodes
      loop_x = loop_xyz/(gyrotropic_kmesh(2)*gyrotropic_kmesh(3))
      loop_y = (loop_xyz - loop_x*(gyrotropic_kmesh(2) &
                                   *gyrotropic_kmesh(3)))/gyrotropic_kmesh(3)
      loop_z = loop_xyz - loop_x*(gyrotropic_kmesh(2)*gyrotropic_kmesh(3)) &
               - loop_y*gyrotropic_kmesh(3)
      kpt(1) = loop_x*db1
      kpt(2) = loop_y*db2
      kpt(3) = loop_z*db3
      kpt(:) = gyrotropic_box_corner(:) + matmul(kpt, gyrotropic_box)

      call gyrotropic_get_k_list(kpt, kweight, &
                                 gyro_K_spn, gyro_K_orb, gyro_D, gyro_Dw, gyro_C, &
                                 gyro_DOS, gyro_NOA_orb, gyro_NOA_spn, &
                                 eval_K, eval_D, eval_Dw, eval_NOA, eval_spn, eval_C, eval_dos)

    end do !loop_xyz

    ! Collect contributions from all nodes
    !
    if (eval_K) then
      call comms_reduce(gyro_K_orb(1, 1, 1), 3*3*nfermi, 'SUM')
      if (eval_spn) call comms_reduce(gyro_K_spn(1, 1, 1), 3*3*nfermi, 'SUM')
    endif

    if (eval_D) &
      call comms_reduce(gyro_D(1, 1, 1), 3*3*nfermi, 'SUM')

    if (eval_C) &
      call comms_reduce(gyro_C(1, 1, 1), 3*3*nfermi, 'SUM')

    if (eval_Dw) &
      call comms_reduce(gyro_Dw(1, 1, 1, 1), 3*3*nfermi*gyrotropic_nfreq, 'SUM')

    if (eval_dos) &
      call comms_reduce(gyro_DOS(1), nfermi, 'SUM')

    if (eval_NOA) then
      call comms_reduce(gyro_NOA_orb(1, 1, 1, 1), 3*3*nfermi*gyrotropic_nfreq, 'SUM')
      if (eval_spn) call comms_reduce(gyro_NOA_spn(1, 1, 1, 1), 3*3*nfermi*gyrotropic_nfreq, 'SUM')
    endif

    if (on_root) then

      if (timing_level > 1) call io_stopwatch('gyrotropic: k-interpolation', 2)
      write (stdout, '(1x,a)') ' '
      write (stdout, *) 'Calculation finished, writing results'
      flush(stdout)

      if (eval_K) then
        if (eval_spn) then
          ! At this point gme_spn_list contains
          ! (1/N) sum_k delta(E_kn-E_f).(d E_{kn}/d k_i).sigma_{kn,j}
          ! (units of length) in Angstroms.
          !        ====================================
          ! To get K in units of Ampere do the following:
          !        ====================================
          !   * Divide by V_c in Ang^3 to get a quantity with units of [L]^{-2}
          !   * Multiply by 10^20 to convert to SI
          !   * Multiply by -g_s.e.hbar/(4m_e) \simeq e.hbar/(2.m_e) in SI units
          ! ==============================
          ! fac = 10^20*e*hbar/(2.m_e.V_c)
          ! ==============================
          fac = -1.0e20_dp*elem_charge_SI*hbar_SI/(2.*elec_mass_SI*cell_volume)
          gyro_K_spn(:, :, :) = gyro_K_spn(:, :, :)*fac
          f_out_name_tmp = 'K_spin'
          units_tmp = "Ampere"
          comment_tmp = "spin part of the K tensor -- Eq. 3 of TAS17"
          call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_K_spn, units=units_tmp, &
                                          comment=comment_tmp)
        endif  ! eval_K && eval_spin

        ! At this point gme_orb_list contains
        ! (1/N)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i)
        !                           .Im[<del_k u_kn| x (H_k-E_kn)|del_k u_kn>]
        ! (units of energy times length^3) in eV.Ang^3.
        !        ====================================
        ! To get K  in units of Ampere do the following:
        !        ====================================
        !   * Divide by V_c in Ang^3 to get a quantity with units of eV
        !   * Multiply by 'e' in SI to convert to SI (Joules)
        !   * Multiply by e/(2.hbar) to get K in Ampere
        ! ====================================
        ! fac = e^2/(2.hbar.V_c)
        ! ====================================
        fac = elem_charge_SI**2/(2.*hbar_SI*cell_volume)
        gyro_K_orb(:, :, :) = gyro_K_orb(:, :, :)*fac

        f_out_name_tmp = 'K_orb'
        units_tmp = "Ampere"
        comment_tmp = "orbital part of the K tensor -- Eq. 3 of TAS17"
        call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_K_orb, units=units_tmp, &
                                        comment=comment_tmp)
      endif ! eval_K

      if (eval_D) then
        fac = 1./cell_volume
        gyro_D(:, :, :) = gyro_D(:, :, :)*fac

        f_out_name_tmp = 'D'
        units_tmp = "dimensionless"
        comment_tmp = "the D tensor -- Eq. 2 of TAS17"
        call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_D, units=units_tmp, &
                                        comment=comment_tmp)
      endif

      if (eval_Dw) then
        fac = 1./cell_volume
        gyro_Dw(:, :, :, :) = gyro_Dw(:, :, :, :)*fac

        f_out_name_tmp = 'tildeD'
        units_tmp = "dimensionless"
        comment_tmp = "the tildeD tensor -- Eq. 12 of TAS17"
        call gyrotropic_outprint_tensor(f_out_name_tmp, arrEfW=gyro_Dw, units=units_tmp, &
                                        comment=comment_tmp)
      endif

      if (eval_C) then
        ! At this point gyro_C contains
        ! (1/N)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i).(d E_{kn}/d k_j)
        ! (units of energy*length^2) in eV*Ang^2
        !
        ! To get it in Cab = e/h * (1/N*V_cell)sum_{k,n} delta(E_kn-E_f).(d E_{kn}/d k_i).(d E_{kn}/d k_j)
        ! in units Ampere/cm
        !
        ! divide by V_c in Ang^3 to get  eV/Ang
        ! multiply by 10^8*e in SI to get J/cm
        ! multiply by e/h in SI
        !
        fac = 1.0e+8_dp*elem_charge_SI**2/(twopi*hbar_SI*cell_volume)
        gyro_C(:, :, :) = gyro_C(:, :, :)*fac

        f_out_name_tmp = 'C'
        units_tmp = "Ampere/cm"
        comment_tmp = "the C tensor -- Eq. B6 of TAS17"
        call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf=gyro_C, units=units_tmp, &
                                        comment=comment_tmp)
      endif

      if (eval_noa) then
        ! at this point gyro_NOA_orb  is in eV^-1.Ang^3   !
        !  We want the result in angstrems  !
        !   * Divide by V_c in Ang^3 to make it eV^{-1}
        !   * Divide by e in SI to get J^{-1}
        !   * multiply by e^2/eps_0 to get meters
        !   *multiply dy 1e10 to get Ang
        fac = 1e+10_dp*elem_charge_SI/(cell_volume*eps0_SI)
        gyro_NOA_orb = gyro_NOA_orb*fac
        f_out_name_tmp = 'NOA_orb'
        units_tmp = "Ang"
        comment_tmp = "the tensor $gamma_{abc}^{orb}$ (Eq. C12,C14 of TAS17)"
        call gyrotropic_outprint_tensor(f_out_name_tmp, arrEfW=gyro_NOA_orb, units=units_tmp, &
                                        comment=comment_tmp, symmetrize=.false.)

        if (eval_spn) then
          ! at this point gyro_NOA_spn  is in eV^-2.Ang   !
          !  We want the result in angstrems  !
          !   * Divide by V_c in Ang^3 to make it (eV.Ang)^{-2}
          !   * multiply by 1e20/e^2 in SI to get (J.m)^{-2}
          !   * multiply by e^2/eps_0 to get (J.m)^{-1}
          !   *multiply dy hbar^2/m_e to get m
          !   *multiply by 1e10 to get Ang
          fac = 1e+30_dp*hbar_SI**2/(cell_volume*eps0_SI*elec_mass_SI)
          gyro_NOA_spn = gyro_NOA_spn*fac
          f_out_name_tmp = 'NOA_spin'
          units_tmp = "Ang"
          comment_tmp = "the tensor $gamma_{abc}^{spin}$ (Eq. C12,C15 of TAS17)"
          call gyrotropic_outprint_tensor(f_out_name_tmp, arrEfW=gyro_NOA_spn, units=units_tmp, &
                                          comment=comment_tmp, symmetrize=.false.)
        endif
      endif  !eval_NOA

      if (eval_DOS) then
        ! At this point gyro_C contains
        ! (1/N)sum_{k,n} delta(E_kn-E_f)
        ! in units of eV^{-1}
        ! divide by V_c in Ang^3 to get units 1./(eV*Ang^3)
        gyro_DOS(:) = gyro_DOS(:)/cell_volume
        f_out_name_tmp = 'DOS'
        units_tmp = "eV^{-1}.Ang^{-3}"
        comment_tmp = "density of states"
        call gyrotropic_outprint_tensor(f_out_name_tmp, arrEf1d=gyro_DOS, units=units_tmp, &
                                        comment=comment_tmp)
      endif

    end if !on_root

  end subroutine gyrotropic_main