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