TB-tbG/src/constants.f90 at df02ee0d2208fc19db3d97cfef7c79d33a19b26c · gawcista/TB-tbG · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
module constants
    implicit none
    logical, public :: timer_mpi = .true.
    logical, public :: timer_cpu = .true.
    logical, public :: timer_kcalc = .false.
    logical, public :: timer_debug = .false.
    logical, public :: write_hr = .false.
    integer, parameter, public :: prec = SELECTED_REAL_KIND(15)
    real(prec), parameter, public :: pi = 3.14159265358979323846
    real(prec), parameter, public :: k_B = 8.617333262145179e-05 ! Boltzmann constant in eV/K
    real(prec), parameter, public :: epsilon_0_SI = 8.8541878188e-12 !
    real(prec), public :: epsilon_r = 4.9 ! from PhysRevB.102.125403
    real(prec), parameter, public :: elec = 1.602176634e-19
    ! hoppings
    real(prec), public :: onsite = -0.78_prec ! on-site energy (eV)
    real(prec), public :: gamma0 = 2.7_prec ! 1-st NN hopping (eV)
    real(prec), public :: gamma1 = 0.48_prec ! interlayer coupling (eV)
    ! geometric parameters (tbG)
    real(prec), parameter :: a0 = 2.4684713048999796_prec ! lattice constant (AA)
    real(prec), parameter :: a_d = a0/3_prec**0.5_prec ! C-C distance (AA)
    real(prec), public :: d0 = 3.380180130823483_prec ! layer distance (AA)
    real(prec), parameter :: l_c = 0.265_prec ! cutoff ratio (AA)
    real(prec), parameter :: r_c = a0*2.5_prec ! cutoff distance
    ! geometric parameters (plasmon)
    logical :: calc_iqr
    real(prec) :: delta ! broadening
    real(prec) :: efermi = 0.0
    real(prec) :: temperature = 300
    ! decay constants
    real(prec), parameter :: f_decay = log(0.1_prec)/(a_d-a0)
    real(prec), parameter :: q_pi = a_d * f_decay

    ! global variants
    logical :: band_select
    logical :: kpt_select
    logical :: job_band
    logical :: write_band
    logical :: eigenvec_out
    logical :: write_matrix
    character(len=3) :: model_type
    character(len=64) :: f_input
    character(len=64) :: f_poscar
    character(len=64) :: f_eig
    character(len=64) :: f_eigvec
    character(len=64) :: f_kpoint
    character(len=64) :: f_wannier
    character(len=64) :: f_hr
    character(len=64) :: f_matrix
    character(len=64) :: f_eels
    character(len=64) :: f_dos
    character(len=1) :: q_tag,job_type
    real(prec) :: basis(3,3), basis_rec(3,3)
    real(prec), allocatable :: position_frac(:,:), position_cart(:,:),omegalist(:),qlist(:,:)
    integer :: nk_path,eels_mode
    integer :: nions,nbands,nomega,iband(2),nedos
    real(prec) :: sigma_dos, gaussian_cutoff, erange(2)
    logical :: dos_normalize

    integer ::  rank, nproc, ncache
contains
    subroutine init_constants()
        band_select = .false.
        kpt_select = .false.
        job_band = .false.
        write_band = .false.
        write_hr = .false.
        write_matrix = .false.
        f_input='input.in'
        f_poscar='cont.vasp'
        f_eig='tb_band.dat'
        f_eigvec='tb_wavef.dat'
        f_kpoint='KPOINTS'
        f_matrix='tb_eels.mat'
        f_hr='tb_hr.dat'
        f_eels='tb_eels.dat'
        f_dos='tb_dos.dat'
        nk_path = 10
        ncache = 0
        nomega = 21
        q_tag = "D"
        delta = 1e-3
        calc_iqr = .false.
        eigenvec_out = .false.
        model_type = "tbg"
        f_wannier = 'wannier90'
        eels_mode = 0
        ! DOS parameters
        sigma_dos = 0.01_prec
        gaussian_cutoff = 0.0_prec
        nedos = 1001
        erange = [-1.0_prec, 1.0_prec]
        dos_normalize = .true.
        call set_omegalist(0._prec,0.1_prec)

    end subroutine init_constants

    subroutine set_omegalist(omega_i,omega_f)
        real(prec), intent(in) :: omega_i,omega_f
        integer :: i
        if (allocated(omegalist)) deallocate(omegalist)
        allocate(omegalist(nomega))
        do i=1, nomega
            omegalist(i)=omega_i+(omega_f-omega_i)*(i-1)/(nomega-1)
        end do
    end subroutine set_omegalist

    function get_label(k_frac) result(label)
        real(prec), intent(in) :: k_frac(3)
        real(prec), parameter :: tol = 1.0e-5_prec
        character(len=3) :: label

        if (all(abs(k_frac - [0.0_prec, 0.0_prec, 0.0_prec]) < tol)) then
            label = ' Γ '
        else if (all(abs(k_frac - [1.0_prec/3.0_prec, 1.0_prec/3.0_prec, 0.0_prec]) < tol)) then
            label = ' K '
        else if (all(abs(k_frac - [0.5_prec, 0.0_prec, 0.0_prec]) < tol)) then
            label = ' M '
        else
            label = ' ? '
        end if

    end function get_label

    subroutine kpath_default(paths)
        real(prec), intent(out), allocatable :: paths(:,:,:)

        if (.not. allocated(paths)) then
            allocate(paths(3,2,3))
        end if

        paths(1,1,:) = [0.0_prec, 0.0_prec, 0.0_prec]
        paths(1,2,:) = [1.0_prec/3.0_prec, 1.0_prec/3.0_prec, 0.0_prec]
        paths(2,1,:) = [1.0_prec/3.0_prec, 1.0_prec/3.0_prec, 0.0_prec]
        paths(2,2,:) = [0.5_prec, 0.0_prec, 0.0_prec]
        paths(3,1,:) = [0.5_prec, 0.0_prec, 0.0_prec]
        paths(3,2,:) = [0.0_prec, 0.0_prec, 0.0_prec]

    end subroutine kpath_default

end module constants