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, 2013, nickpapior@gmail.com
10! Please conctact the author, prior to re-using this code.
11
12module m_ts_full_scat
13
14  use precision, only : dp
15
16  use m_ts_electype
17  use m_ts_cctype
18
19  implicit none
20
21  private
22
23  public :: calc_GF
24  public :: calc_GF_Bias
25  public :: calc_GF_Part
26  public :: GF_Gamma_GF
27  public :: insert_Self_Energies
28
29contains
30
31
32! Full converted GF.G.GF^\dagger routine for speed.
33! This routine is extremely fast compared to any previous implementation.
34! It relies on the fact that Gf only contains the electrode columns.
35  subroutine GF_Gamma_GF(El, no_u_TS, no, GF, &
36       GGG,nwork,work)
37
38!  This routine returns GGG=GF.Gamma.GF^\dagger, where GF is a (no_u)x(no)
39!  matrix and the states
40!  corresponds to the (no) electrode states
41!  Gamma is a (no)x(no) matrix.
42
43! *********************
44! * INPUT variables   *
45! *********************
46    ! electrode self-energy
47    type(Elec), intent(in) :: El
48    integer, intent(in) :: no_u_TS ! no. states in contact region
49    integer, intent(in) :: no      ! no. states for this electrode
50    ! The Green function (it has to be the column that corresponds to the electrode)
51    complex(dp), intent(inout) :: GF(no_u_TS,no)
52    ! A work array for doing the calculation... (nwork has to be larger than no_u_TS)
53    integer,     intent(in)    :: nwork
54    complex(dp), intent(inout) :: work(nwork)
55
56! *********************
57! * OUTPUT variables  *
58! *********************
59    complex(dp), intent(out) :: GGG(no_u_TS*no_u_TS)    !GF.Gamma.GF^\dagger
60
61! *********************
62! * LOCAL variables   *
63! *********************
64    complex(dp), parameter :: z0 = dcmplx(0._dp, 0._dp)
65    complex(dp), parameter :: z1 = dcmplx(1._dp, 0._dp)
66    complex(dp), parameter :: zi = dcmplx(0._dp, 1._dp)
67
68    integer :: i, NB, ind, iB
69
70#ifdef TRANSIESTA_DEBUG
71    call write_debug( 'PRE GFGammaGF' )
72#endif
73
74    call timer("GFGGF",1)
75
76    ! Number of times we can divide the large matrix
77    NB = no_u_TS / no
78
79    ! Loop over bottom row matrix
80    do iB = 0 , NB - 1
81
82       ! Collect the top row of complex conjugated Gf
83       ind = no_u_TS * no * iB + 1
84       do i = 1 , no
85          GGG(ind:ind-1+no) = dconjg(Gf(iB*no+1:(iB+1)*no,i))
86          ind = ind + no
87       end do
88       ind = no_u_TS * no * iB + 1
89
90       ! Do Gamma.Gf^\dagger
91#ifdef USE_GEMM3M
92       call zgemm3m( &
93#else
94       call zgemm( &
95#endif
96            'T','T',no,no,no,z1, &
97            El%Gamma, no, &
98            GGG(ind), no, &
99            z0, work,no)
100
101       ! Calculate the Gf.Gamma.Gf^\dagger product for the entire column
102#ifdef USE_GEMM3M
103       call zgemm3m( &
104#else
105       call zgemm( &
106#endif
107            'N','N',no_u_TS,no,no,z1, &
108            Gf(1,1), no_u_TS, &
109            work   ,      no, &
110            z0, GGG(ind),no_u_TS)
111
112    end do
113
114    ! in case the block size does not match the matrix order
115    if ( NB * no /= no_u_TS ) then
116
117       ! The size of the remaining block
118       iB = no_u_TS - NB * no
119
120       ! Copy over the block
121       ind = no_u_TS * no * NB + 1
122       do i = 1 , no
123          ! So this is the complex conjugated of the iB'th block
124          GGG(ind:ind-1+iB) = dconjg(Gf(NB*no+1:NB*no+iB,i))
125          ind = ind + iB
126       end do
127       ind = no_u_TS * no * NB + 1
128
129       ! Do Gamma.Gf^\dagger
130#ifdef USE_GEMM3M
131       call zgemm3m( &
132#else
133       call zgemm( &
134#endif
135            'T','T',no,iB,no,z1, &
136            El%Gamma, no, &
137            GGG(ind), iB, &
138            z0, work,no)
139
140       ! Calculate the Gf.Gamma.Gf^\dagger product for the entire column
141#ifdef USE_GEMM3M
142       call zgemm3m( &
143#else
144       call zgemm( &
145#endif
146            'N','N',no_u_TS,iB,no,z1, &
147            Gf(1,1), no_u_TS, &
148            work   ,      no, &
149            z0, GGG(ind),no_u_TS)
150
151    end if
152
153
154#ifdef TRANSIESTA_31
155    ! Lets try and impose symmetry...
156    call my_symmetrize(no_u_TS,GGG)
157#endif
158
159    call timer("GFGGF",2)
160
161#ifdef TRANSIESTA_DEBUG
162    call write_debug( 'POS GFGammaGF' )
163#endif
164
165#ifdef TRANSIESTA_31
166  contains
167    subroutine my_symmetrize(N,M)
168      integer    , intent(in) :: N
169      complex(dp), intent(inout) :: M(N,N)
170      integer :: i,j
171      do j = 1 , N
172        do i = 1 , j
173          M(j,i) = aimag(M(i,j))
174          M(i,j) = M(j,i)
175        end do
176      end do
177    end subroutine my_symmetrize
178#endif
179
180  end subroutine GF_Gamma_GF
181
182
183! ##################################################################
184! ## Calculating Full Green functions of                          ##
185! ##                                                              ##
186! ##                            By                                ##
187! ##              Mads Brandbyge, mbr@mic.dtu.dk                  ##
188! ##                                                              ##
189! ##                                                              ##
190! ## Completely restructured to be able to handle sparse matrices ##
191! ##                                                              ##
192! ##                                                              ##
193! ##  Modified by Nick Papior Andersen                            ##
194! ##################################################################
195  subroutine calc_GF(cE,no_u_TS,GFinv,GF)
196
197    use intrinsic_missing, only: EYE
198    use precision, only: dp
199
200    implicit none
201
202! *********************
203! * INPUT variables   *
204! *********************
205    type(ts_c_idx), intent(in) :: cE
206    ! Sizes of the different regions...
207    integer, intent(in) :: no_u_TS
208    ! Work should already contain Z*S - H
209    ! This may seem strange, however, it will clean up this routine extensively
210    ! as we dont need to make two different routines for real and complex
211    ! Hamiltonian values.
212    complex(dp), intent(in out) :: GFinv(no_u_TS,no_u_TS) ! the inverted GF
213    complex(dp), intent(out) :: GF(no_u_TS,no_u_TS)
214
215! Local variables
216    integer :: ipvt(no_u_TS), ierr
217
218    if ( cE%fake ) return
219
220#ifdef TRANSIESTA_DEBUG
221    call write_debug( 'PRE getGF' )
222#endif
223
224    call timer('GFT',1)
225
226    call EYE(no_u_TS,GF)
227
228    ! Invert directly
229    call zgesv(no_u_TS,no_u_TS,GFinv,no_u_TS,ipvt,GF,no_u_TS,ierr)
230    if ( ierr /= 0 ) call die('GF: Could not invert the Green function')
231
232    call timer('GFT',2)
233
234#ifdef TRANSIESTA_DEBUG
235    call write_debug( 'POS getGF' )
236#endif
237
238  end subroutine calc_GF
239
240
241! ##################################################################
242! ## Calculating Green functions in the reigon of the electrodes  ##
243! ##                                                              ##
244! ##  Fully created by Nick Papior Andersen, nickpapior@gmail.com ##
245! ##################################################################
246  subroutine calc_GF_Bias(cE,no_u_TS,no_Els,N_Elec,Elecs,GFinv,GF)
247
248    use precision, only: dp
249
250    use m_ts_method, only : orb_offset
251
252    implicit none
253
254! *********************
255! * INPUT variables   *
256! *********************
257    type(ts_c_idx), intent(in) :: cE
258    ! Sizes of the different regions...
259    integer, intent(in) :: no_u_TS, no_Els
260    ! Electrodes
261    integer, intent(in) :: N_Elec
262    type(Elec), intent(in) :: Elecs(N_Elec)
263    ! Work should already contain Z*S - H
264    ! This may seem strange, however, it will clean up this routine extensively
265    ! as we dont need to make two different routines for real and complex
266    ! Hamiltonian values.
267    complex(dp), intent(in out) :: GFinv(no_u_TS,no_u_TS) ! the inverted GF
268    ! We only need Gf in the left and right blocks...
269    complex(dp), intent(out) :: GF(no_u_TS,no_Els)
270
271! Local variables
272    integer :: ipvt(no_u_TS)
273    integer :: i, o, iEl, off_row
274
275    if ( cE%fake ) return
276
277#ifdef TRANSIESTA_DEBUG
278    call write_debug( 'PRE getGF' )
279#endif
280
281    call timer('GFTB',1)
282
283    ! Create the RHS for inversion...
284    GF(:,:) = dcmplx(0._dp,0._dp)
285
286    o = 0
287    do iEl = 1 , N_Elec
288      i = Elecs(iEl)%idx_o
289      off_row = i - orb_offset(i) - 1
290      do i = 1 , TotUsedOrbs(Elecs(iEl))
291        o = o + 1
292        GF(off_row+i,o) = dcmplx(1._dp,0._dp)
293      end do
294    end do
295
296    if ( o /= no_Els ) call die('GFB: Error in sizes of electrodes')
297
298    ! Invert directly
299    call zgesv(no_u_TS,no_Els,GFinv,no_u_TS,ipvt,GF,no_u_TS,i)
300    if ( i /= 0 ) call die('GFB: Could not invert the Green function')
301
302    call timer('GFTB',2)
303
304#ifdef TRANSIESTA_DEBUG
305    call write_debug( 'POS getGF' )
306#endif
307
308  end subroutine calc_GF_Bias
309
310
311  ! Calculate a sub-part of the full Green function.
312  ! This routine calculates the dense part that does not overlap with
313  ! electrodes.
314  ! I.e. this can *only* be used for equilibrium contour points
315  ! since it removes the columns for electrodes where DM_update == 0.
316  subroutine calc_GF_Part(cE,no_u, no_u_TS, no_col, N_Elec, Elecs, & ! Size of the problem
317       GFinv,GF)
318
319    use intrinsic_missing, only: EYE
320    use precision, only: dp
321    use m_ts_method, only : orb_offset, orb_type, TYP_BUFFER
322
323    implicit none
324
325! *********************
326! * INPUT variables   *
327! *********************
328    type(ts_c_idx), intent(in) :: cE
329    ! Sizes of the different regions...
330    integer, intent(in) :: no_u, no_u_TS, no_col
331    integer, intent(in) :: N_Elec
332    type(Elec), intent(in) :: Elecs(N_Elec)
333    ! Work should already contain Z*S - H
334    ! This may seem strange, however, it will clean up this routine extensively
335    ! as we dont need to make two different routines for real and complex
336    ! Hamiltonian values.
337    complex(dp), intent(in out) :: GFinv(no_u_TS,no_u_TS) ! the inverted GF
338    complex(dp), intent(out) :: GF(no_u_TS,no_col)
339
340! Local variables
341    integer :: ipvt(no_u_TS)
342    integer :: jo, i, j
343    logical :: dm_update_0(N_elec)
344
345    if ( cE%fake ) return
346
347#ifdef TRANSIESTA_DEBUG
348    call write_debug( 'PRE getGF' )
349#endif
350
351    call timer('GFT_P',1)
352
353    dm_update_0(:) = Elecs(:)%DM_update == 0
354
355    ! initialize
356    GF(:,:) = dcmplx(0._dp,0._dp)
357
358    j = 0
359    i = 0
360    do jo = 1, no_u
361
362      ! if buffer, skip!
363      if ( orb_type(jo) == TYP_BUFFER ) cycle
364      if ( any(OrbInElec(Elecs, jo) .and. dm_update_0(:)) ) then
365        i = i + 1
366        cycle
367      end if
368      i = i + 1
369      j = j + 1
370
371      GF(i,j) = dcmplx(1._dp,0._dp)
372    end do
373
374    if ( i /= no_u_TS .or. j /= no_col ) call die('GFP: error in constructing part GF')
375
376    ! Invert directly
377    call zgesv(no_u_TS,no_col,GFinv,no_u_TS,ipvt,GF,no_u_TS,i)
378    if ( i /= 0 ) call die('GFP: Could not invert the Green function')
379
380    call timer('GFT_P',2)
381
382#ifdef TRANSIESTA_DEBUG
383    call write_debug( 'POS getGF' )
384#endif
385
386  end subroutine calc_GF_Part
387
388  subroutine insert_Self_Energies(no_u, Gfinv, El)
389    use m_ts_method, only : orb_offset
390    integer, intent(in) :: no_u
391    complex(dp), intent(in out) :: GFinv(no_u,no_u)
392    type(Elec), intent(in) :: El
393
394    integer :: i, j, ii, jj, iii, off, no
395
396    no = TotUsedOrbs(El)
397    off = El%idx_o - orb_offset(El%idx_o) - 1
398
399    if ( El%Bulk ) then
400!$OMP do private(j,jj,i,ii)
401       do j = 0 , no - 1
402          jj = off + j + 1
403          ii = j * no
404          do i = 1 , no
405             Gfinv(off+i,jj) = El%Sigma(ii+i)
406          end do
407       end do
408!$OMP end do nowait
409    else
410!$OMP do private(j,jj,i,ii,iii)
411       do j = 0 , no - 1
412          jj = off + j + 1
413          ii = j * no
414          do i = 1 , no
415             iii = off + i
416             Gfinv(iii,jj) = Gfinv(iii,jj) - El%Sigma(ii+i)
417          end do
418       end do
419!$OMP end do nowait
420    end if
421
422  end subroutine insert_Self_Energies
423
424end module m_ts_full_scat
425