1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8! This code segment has been fully created by:
9! Nick Papior Andersen, 2012, nickpapior@gmail.com
10!
11
12module m_ts_contour
13!
14! Routines that are used to setup the contour for integration of the GFs
15!
16! Use the type associated with the contour
17! Maybe they should be collected to this module.
18! However, I like this partition.
19  use m_ts_electype
20  use m_ts_chem_pot
21
22  use m_ts_cctype
23  use m_ts_io_contour
24
25  use m_ts_contour_eq,  only : CONTOUR_EQ
26  use m_ts_contour_neq, only : CONTOUR_NEQ
27
28  use precision, only : dp
29
30  implicit none
31
32  private
33
34  public :: read_contour_options
35  public :: print_contour_options
36  public :: print_contour_block
37  public :: io_contour
38  public :: sort_contour
39
40contains
41
42  subroutine io_contour(IsVolt, mus, slabel,suffix)
43    use m_ts_contour_eq, only : io_contour_Eq
44    use m_ts_contour_neq, only : io_contour_nEq
45    logical, intent(in) :: IsVolt
46    type(ts_mu), intent(in) :: mus(:)
47    character(len=*), intent(in) :: slabel
48    character(len=*), intent(in), optional :: suffix
49
50    call io_contour_Eq(mus,slabel,suffix)
51
52    if ( IsVolt ) then
53       call io_contour_nEq(slabel,suffix)
54    end if
55
56  end subroutine io_contour
57
58  subroutine read_contour_options(N_Elec, Elecs, N_mu, mus, kT, IsVolt, Volt)
59    use m_ts_contour_eq, only : read_contour_eq_options
60    use m_ts_contour_neq, only : read_contour_neq_options
61
62    integer, intent(in) :: N_Elec
63    type(Elec), intent(inout) :: Elecs(N_Elec)
64    integer, intent(in) :: N_mu
65    type(ts_mu), intent(inout) :: mus(N_mu)
66    ! SIESTA electronic temperature
67    real(dp), intent(in) :: kT
68    logical, intent(in) :: IsVolt
69    real(dp), intent(in) :: Volt
70
71    call read_contour_eq_options(N_mu,mus,Volt)
72
73    if ( IsVolt ) then
74       call read_contour_neq_options(N_Elec,Elecs,N_mu,mus,kT,Volt)
75    end if
76
77  end subroutine read_contour_options
78
79  subroutine print_contour_options(prefix,IsVolt)
80    use m_ts_contour_eq, only : print_contour_eq_options
81    use m_ts_contour_neq, only : print_contour_neq_options
82
83    character(len=*), intent(in) :: prefix
84    logical, intent(in) :: IsVolt
85
86    call print_contour_eq_options(prefix)
87
88    if ( IsVolt ) then
89       call print_contour_neq_options(prefix)
90    end if
91
92  end subroutine print_contour_options
93
94  subroutine print_contour_block(prefix, IsVolt)
95    use m_ts_contour_eq, only : print_contour_eq_block
96    use m_ts_contour_neq, only : print_contour_neq_block
97
98    character(len=*), intent(in) :: prefix
99    logical, intent(in) :: IsVolt
100
101    call print_contour_eq_block(prefix)
102
103    if ( IsVolt ) then
104       call print_contour_neq_block(prefix)
105    end if
106
107  end subroutine print_contour_block
108
109  ! In order to retain the best numerical accuracy we introduce
110  ! an ordering of the energy contour by weights.
111  ! It is the only reasonable thing to access as the functional is
112  ! non-deterministic.
113  ! An example of the importance of this sorting:
114  ! Take +200 voltage energy points.
115  ! The weights will be something like this:
116  ! 1e-5,...,1e-4,...,1e-3,...,1e-4,...,1e-5
117  ! Which means that the summation up till the half works great.
118  ! But when we reach the last weight we could be in the situation where:
119  !  1._dp + 1e-12_dp which will limit the accuracy obtained for that energy point.
120
121  ! In order to circumvent this we simply sort by weights.
122  subroutine sort_contour(NC,ce,cw)
123    integer, intent(in) :: NC
124    complex(dp), intent(inout) :: ce(NC), cw(NC)
125
126    ! Local variables
127    real(dp) :: tmp
128    complex(dp) :: ctmp
129    integer :: i,j
130    ! As we will only do this once we dont need a fancy sorting
131    ! algorithm...
132    do i = 1 , NC - 1
133       tmp = real(cw(i)*conjg(cw(i)),dp)
134       do j = i+1, NC
135          if ( tmp > real(cw(j)*conjg(cw(j))) ) then
136             ctmp  = ce(i)
137             ce(i) = ce(j)
138             ce(j) = ctmp
139             ctmp  = cw(i)
140             cw(i) = cw(j)
141             cw(j) = ctmp
142          end if
143       end do
144    end do
145
146  end subroutine sort_contour
147
148end module m_ts_contour
149