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! * It has been heavily inspired by the original authors of the 12! Transiesta code (hence the references here are still remaining) * 13 14! This particular solution method relies on solving the GF 15! with the tri-diagonalization routine. 16! This will leverage memory usage and also the execution time. 17 18module m_ts_trik 19 20 use precision, only : dp 21 use m_region 22 23 use m_ts_sparse_helper 24 25 use m_ts_dm_update, only : init_DM 26 use m_ts_dm_update, only : update_DM, update_zDM 27 use m_ts_dm_update, only : add_k_DM 28 29 use m_ts_weight, only : weight_DM 30 use m_ts_weight, only : TS_W_K_METHOD 31 use m_ts_weight, only : TS_W_K_CORRELATED 32 use m_ts_weight, only : TS_W_K_UNCORRELATED 33 34 use m_ts_tri_init, only : c_Tri 35 36 use m_ts_method, only : orb_offset, no_Buf, r_pvt 37 use m_ts_method, only : ts_A_method, TS_BTD_A_COLUMN 38 39 implicit none 40 41 public :: ts_trik 42 43 private 44 45contains 46 47 subroutine ts_trik(N_Elec,Elecs, & 48 nq, uGF, ucell, nspin, na_u, lasto, & 49 sp_dist, sparse_pattern, & 50 no_u, n_nzs, & 51 Hs, Ss, DM, EDM, Ef, & 52 DE_NEGF) 53 54 use units, only : Pi 55 use parallel, only : Node, Nodes 56 57#ifdef MPI 58 use mpi_siesta 59#endif 60 61 use alloc, only : re_alloc, de_alloc 62 63 use class_OrbitalDistribution 64 use class_Sparsity 65 use class_zSpData1D 66 use class_dSpData2D 67 use class_zSpData2D 68 use class_zTriMat 69 70 use m_ts_electype 71 ! Self-energy read 72 use m_ts_gf 73 ! Self-energy expansion 74 use m_ts_elec_se 75 76 use m_ts_kpoints, only : ts_nkpnt, ts_kpoint, ts_kweight 77 78 use m_ts_options, only : Calc_Forces 79 use m_ts_options, only : N_mu, mus 80 81 use m_ts_options, only : IsVolt 82 83 use m_ts_sparse, only : ts_sp_uc, tsup_sp_uc 84 use m_ts_sparse, only : ltsup_sp_sc, ltsup_sc_pnt 85 use m_ts_sparse, only : sc_off 86 87 use m_ts_cctype 88 use m_ts_contour_eq, only : Eq_E, ID2idx, c2weight_eq, c2energy 89 use m_ts_contour_neq, only : nEq_E, has_cE_nEq 90 use m_ts_contour_neq, only : N_nEq_ID, c2weight_neq 91 use m_ts_contour_eq, only : N_Eq_E 92 use m_ts_contour_neq, only : N_nEq_E 93 94 use m_iterator 95 use m_mat_invert 96 97 use m_trimat_invert 98 99 ! Gf calculation 100 use m_ts_trimat_invert 101 102 use m_ts_tri_common, only : GFGGF_needed_worksize, nnzs_tri 103 104 ! Gf.Gamma.Gf 105 use m_ts_tri_scat 106 107 use ts_dq_m, only: ts_dq 108 109#ifdef TRANSIESTA_DEBUG 110 use m_ts_debug 111#endif 112 113! ******************** 114! * INPUT variables * 115! ******************** 116 integer, intent(in) :: N_Elec 117 type(Elec), intent(inout) :: Elecs(N_Elec) 118 integer, intent(in) :: nq(N_Elec), uGF(N_Elec) 119 real(dp), intent(in) :: ucell(3,3) 120 integer, intent(in) :: nspin, na_u, lasto(0:na_u) 121 type(OrbitalDistribution), intent(inout) :: sp_dist 122 type(Sparsity), intent(inout) :: sparse_pattern 123 integer, intent(in) :: no_u 124 integer, intent(in) :: n_nzs 125 real(dp), intent(in) :: Hs(n_nzs,nspin), Ss(n_nzs) 126 real(dp), intent(inout) :: DM(n_nzs,nspin), EDM(n_nzs,nspin) 127 real(dp), intent(in) :: Ef 128 real(dp), intent(inout) :: DE_NEGF 129 130! ******************* Computational arrays ******************* 131 integer :: nzwork, n_s 132 complex(dp), pointer :: zwork(:) 133 type(zTriMat) :: zwork_tri, GF_tri 134 ! A local orbital distribution class (this is "fake") 135 type(OrbitalDistribution) :: fdist 136 ! The Hamiltonian and overlap sparse matrices 137 type(zSpData1D) :: spH, spS 138 ! local sparsity pattern in local SC pattern 139 type(dSpData2D) :: spDM, spDMneq 140 type(dSpData2D) :: spEDM ! only used if calc_forces 141 ! The different sparse matrices that will surmount to the integral 142 ! These two lines are in global update sparsity pattern (UC) 143 type(zSpData2D) :: spuDM 144 type(zSpData2D) :: spuEDM ! only used if calc_forces 145 ! To figure out which parts of the tri-diagonal blocks we need 146 ! to calculate 147 logical, pointer :: calc_parts(:) => null() 148! ************************************************************ 149 150! ****************** Electrode variables ********************* 151 integer :: padding, GFGGF_size ! with IsVolt we need padding and work-array 152 complex(dp), pointer :: GFGGF_work(:) => null() 153! ************************************************************ 154 155! ******************* Computational variables **************** 156 type(ts_c_idx) :: cE 157 integer :: NE 158 integer :: index_dq !< Index for the current charge calculation @ E == mu 159 real(dp) :: kw, kpt(3), bkpt(3), dq_mu 160 complex(dp) :: W, ZW 161 type(tRgn) :: pvt 162! ************************************************************ 163 164! ******************** Loop variables ************************ 165 type(itt2) :: SpKp 166 integer, pointer :: ispin, ikpt 167 integer :: iEl, iID 168 integer :: iE, imu, io, idx 169 integer :: no 170! ************************************************************ 171 172#ifdef TRANSIESTA_DEBUG 173 call write_debug( 'PRE transiesta mem' ) 174#endif 175 176 ! Create the back-pivoting region 177 call rgn_init(pvt,nrows_g(sparse_pattern),val=0) 178 do io = 1 , r_pvt%n 179 pvt%r(r_pvt%r(io)) = io 180 end do 181 182 ! Number of supercells 183 n_s = size(sc_off,dim=2) 184 185 ! We do need the full GF AND a single work array to handle the 186 ! left-hand side of the inversion... 187 ! We will provide all work arrays as single dimension arrays. 188 ! This will make interfaces more stringent and allow for 189 ! re-use in several other places. 190 ! However, this comes at the cost of programmer book-keeping. 191 ! It should be expected that the work arrays return GARBAGE 192 ! on ALL routines, i.e. they are not used for anything other 193 ! than, yes, work. 194 195 ! The zwork is needed to construct the LHS for solving: G^{-1} G = I 196 ! Hence, we will minimum require this... 197 if ( ts_A_method == TS_BTD_A_COLUMN .and. IsVolt ) then 198 call GFGGF_needed_worksize(c_Tri%n,c_Tri%r, & 199 N_Elec, Elecs, padding, GFGGF_size) 200 else 201 padding = 0 202 GFGGF_size = 0 203 end if 204 ! Now figure out the required worksize for SE expansion 205 call UC_minimum_worksize(IsVolt, N_Elec, Elecs, idx) 206 io = nnzs_tri(c_Tri%n, c_Tri%r) 207 padding = max(padding, idx - io) 208 call newzTriMat(zwork_tri,c_Tri%n,c_Tri%r,'GFinv', & 209 padding=padding) 210 nzwork = elements(zwork_tri,all=.true.) 211 212 ! Initialize the tri-diagonal inversion routine 213 call init_TriMat_inversion(zwork_tri) 214 215 call newzTriMat(GF_tri,c_Tri%n,c_Tri%r,'GF') 216 217 ! initialize the matrix inversion tool 218 call init_mat_inversion(maxval(c_Tri%r)) 219 220 ! Allocate the logical array to handle calculated 221 ! entries in the block-matrix 222 call re_alloc(calc_parts,1,c_Tri%n) 223 ! initialize to calculate all blocks 224 calc_parts(:) = .true. 225 226 ! we use the GF as a placement for the self-energies 227 no = 0 228 zwork => val(GF_tri,all=.true.) 229 do iEl = 1 , N_Elec 230 231 ! This seems stupid, however, we never use the Sigma and 232 ! GF at the same time. Hence it will be safe 233 ! to have them point to the same array. 234 ! When the UC_expansion_Sigma_GammaT is called 235 ! first the Sigma is assigned and then 236 ! it is required that prepare_GF_inv is called 237 ! immediately (which it is) 238 ! Hence the GF must NOT be used in between these two calls! 239 io = TotUsedOrbs(Elecs(iEl)) ** 2 240 Elecs(iEl)%Sigma => zwork(no+1:no+io) 241 no = no + io 242 243 ! if we need the cross-terms we can not skip the blocks 244 ! that are fully inside the electrode 245 if ( Elecs(iEl)%DM_update /= 0 ) cycle 246 247 io = Elecs(iEl)%idx_o 248 io = io - orb_offset(io) 249 idx = io + TotUsedOrbs(Elecs(iEl)) - 1 250 251 end do 252 253 ! Save the work-space 254 ! Now the programmer should "keep a straight tongue" 255 ! The zwork points to the array in the zwork_tri 256 ! tri-diagonal array. This means that there are two 257 ! arrays that point to the same. 258 ! Generally the zwork need only to retain the value in 259 ! one call! 260 zwork => val(zwork_tri,all=.true.) 261 262 ! Create the Fake distribution 263 ! The Block-size is the number of orbitals, i.e. all on the first processor 264 ! Notice that we DO need it to be the SIESTA size. 265#ifdef MPI 266 call newDistribution(no_u,MPI_Comm_Self,fdist,name='TS-fake dist') 267#else 268 call newDistribution(no_u,-1 ,fdist,name='TS-fake dist') 269#endif 270 271 ! The Hamiltonian and overlap matrices (in Gamma calculations 272 ! we will not have any phases, hence, it makes no sense to 273 ! have the arrays in complex) 274 ! TODO move into a Data2D (could reduce overhead of COMM) 275 call newzSpData1D(ts_sp_uc,fdist,spH,name='TS spH') 276 call newzSpData1D(ts_sp_uc,fdist,spS,name='TS spS') 277 278 ! If we have a bias calculation we need additional arrays. 279 ! If not bias we don't need the update arrays (we already have 280 ! all information in tsup_sp_uc (spDMu)) 281 282 ! Allocate space for global sparsity arrays 283 no = max(N_mu,N_nEq_id) 284 call newzSpData2D(tsup_sp_uc,no,fdist, spuDM, name='TS spuDM') 285 if ( Calc_Forces ) then 286 call newzSpData2D(tsup_sp_uc,N_mu,fdist, spuEDM, name='TS spuEDM') 287 end if 288 289 if ( IsVolt ) then 290 ! Allocate space for update arrays, local sparsity arrays 291 call newdSpData2D(ltsup_sp_sc,N_mu, sp_dist,spDM ,name='TS spDM') 292 call newdSpData2D(ltsup_sp_sc,N_nEq_id,sp_dist,spDMneq,name='TS spDM-neq') 293 if ( Calc_Forces ) then 294 call newdSpData2D(ltsup_sp_sc,N_mu, sp_dist,spEDM ,name='TS spEDM') 295 end if 296 end if 297 298 if ( ts_A_method == TS_BTD_A_COLUMN .and. IsVolt ) then 299 ! we need only allocate one work-array for 300 ! Gf.G.Gf^\dagger 301 call re_alloc(GFGGF_work,1,GFGGF_size,routine='transiesta') 302 end if 303 304#ifdef TRANSIESTA_WEIGHT_DEBUG 305 do iel = 1 , n_mu 306 print '(a20,tr1,i3)','k '//trim(mus(iEl)%name),iel 307 end do 308#endif 309 310 ! Initialize the charge correction scheme (will return if not used) 311 call ts_dq%initialize_dq() 312 313 ! Total number of energy points 314 NE = N_Eq_E() + N_nEq_E() 315 316 ! start the itterators 317 call itt_init (SpKp,end1=nspin,end2=ts_nkpnt) 318 ! point to the index iterators 319 call itt_attach(SpKp,cur1=ispin,cur2=ikpt) 320 321 do while ( .not. itt_step(SpKp) ) 322 323 if ( itt_stepped(SpKp,1) ) then 324 call init_DM(sp_dist,sparse_pattern, & 325 n_nzs, DM(:,ispin), EDM(:,ispin), & 326 tsup_sp_uc, Calc_Forces) 327 do iEl = 1, N_Elec 328 call reread_Gamma_Green(Elecs(iEl), uGF(iEl), NE, ispin) 329 end do 330 end if 331 332 ! Include spin factor and 1/(2\pi) 333 kpt(:) = ts_kpoint(:,ikpt) 334 ! create the k-point in reciprocal space 335 call kpoint_convert(ucell,kpt,bkpt,1) 336 337 ! Include spin factor and 1/\pi 338 ! Since we are calculating G - G^\dagger in the equilibrium contour 339 ! we need *half* the weight 340 kw = ts_kweight(ikpt) / (Pi * nspin) 341 342#ifdef TRANSIESTA_TIMING 343 call timer('TS_HS',1) 344#endif 345 346 ! Work-arrays are for MPI distribution... 347 call create_HS(sp_dist,sparse_pattern, & 348 Ef, & 349 N_Elec, Elecs, no_u, n_s, & ! electrodes, SIESTA size 350 n_nzs, Hs(:,ispin), Ss, sc_off, & 351 spH, spS, kpt, & 352 nzwork, zwork) 353 354#ifdef TRANSIESTA_TIMING 355 call timer('TS_HS',2) 356#endif 357 358 359#ifdef TRANSIESTA_TIMING 360 call timer('TS_EQ',1) 361#endif 362 363 ! *************** 364 ! * EQUILIBRIUM * 365 ! *************** 366 call init_val(spuDM) 367 if ( Calc_Forces ) call init_val(spuEDM) 368 iE = Nodes - Node 369 cE = Eq_E(iE,step=Nodes) ! we read them backwards 370 do while ( cE%exist ) 371 372 ! ******************* 373 ! * prep Sigma * 374 ! ******************* 375 call read_next_GS(ispin, ikpt, bkpt, & 376 cE, N_Elec, uGF, Elecs, & 377 nzwork, zwork, .false., forward = .false. ) 378 do iEl = 1 , N_Elec 379 call UC_expansion(cE, Elecs(iEl), nzwork, zwork, & 380 non_Eq = .false. ) 381 end do 382 383 ! ******************* 384 ! * prep GF^-1 * 385 ! ******************* 386 call prepare_invGF(cE, zwork_tri, r_pvt, pvt, & 387 N_Elec, Elecs, & 388 spH=spH , spS=spS ) 389 390 ! ******************* 391 ! * calc GF * 392 ! ******************* 393 if ( .not. cE%fake ) then 394 call invert_TriMat(zwork_tri,GF_tri,calc_parts) 395 end if 396 397 ! ** At this point we have calculated the Green function 398 399 ! **************** 400 ! * save GF * 401 ! **************** 402 do imu = 1 , N_mu 403 if ( cE%fake ) cycle ! in case we implement an MPI communication solution... 404 call ID2idx(cE,mus(imu)%ID,idx) 405 if ( idx < 1 ) cycle 406 407 call c2weight_eq(cE,idx, kw, W ,ZW) 408 409 ! Figure out if this point is a charge-correction energy 410 index_dq = ts_dq%get_index(imu, iE) 411 if ( index_dq > 0 ) then 412 ! Accummulate charge at the electrodes chemical potentials 413 ! Note that this dq_mu does NOT have the prefactor 1/Pi 414 call add_DM( spuDM, W, spuEDM, ZW, & 415 GF_tri, r_pvt, pvt, & 416 N_Elec, Elecs, & 417 DMidx=mus(imu)%ID, & 418 spS=spS, q=dq_mu) 419 ts_dq%mus(imu)%dq(index_dq) = ts_dq%mus(imu)%dq(index_dq) + dq_mu * kw 420 else 421 call add_DM( spuDM, W, spuEDM, ZW, & 422 GF_tri, r_pvt, pvt, & 423 N_Elec, Elecs, & 424 DMidx=mus(imu)%ID) 425 end if 426 end do 427 428 ! step energy-point 429 iE = iE + Nodes 430 cE = Eq_E(iE,step=Nodes) ! we read them backwards 431 end do 432 433#ifdef TRANSIESTA_TIMING 434 call timer('TS_EQ',2) 435#endif 436 437#ifdef MPI 438 ! We need to reduce all the arrays 439 call MPI_Barrier(MPI_Comm_World,io) 440 call timer('TS_comm',1) 441 call AllReduce_SpData(spuDM,nzwork,zwork,N_mu) 442 if ( Calc_Forces ) then 443 call AllReduce_SpData(spuEDM,nzwork,zwork,N_mu) 444 end if 445 call timer('TS_comm',2) 446#endif 447 448 if ( .not. IsVolt ) then 449 call update_zDM(sp_dist,sparse_pattern, n_nzs, & 450 DM(:,ispin), spuDM, Ef, & 451 EDM(:,ispin), spuEDM, kpt, n_s, sc_off) 452 453 ! The remaining code segment only deals with 454 ! bias integration... So we skip instantly 455 do iEl = 1, N_Elec 456 call reread_Gamma_Green(Elecs(iEl), uGF(iEl), NE, ispin) 457 end do 458 459 cycle 460 461 end if 462 463 ! ***************** 464 ! * only things with non-Equilibrium contour... 465 ! ***************** 466 467 ! initialize to zero 468 ! local sparsity update patterns 469 if ( TS_W_K_METHOD == TS_W_K_UNCORRELATED ) then 470 call init_val(spDM) 471 call init_val(spDMneq) 472 if ( Calc_Forces ) call init_val(spEDM) 473 else if ( itt_stepped(SpKp,1) ) then 474 ! we only need to initialize once per spin 475 call init_val(spDM) 476 call init_val(spDMneq) 477 if ( Calc_Forces ) call init_val(spEDM) 478 end if 479 480 ! transfer data to local sparsity arrays 481 call add_k_DM(spDM, spuDM, N_mu, & 482 spEDM, spuEDM, N_mu, & 483 n_s, sc_off, kpt, non_Eq = .false. ) 484 485#ifdef TRANSIESTA_TIMING 486 call timer('TS_NEQ',1) 487#endif 488 489 ! ******************* 490 ! * NON-EQUILIBRIUM * 491 ! ******************* 492 493 call init_val(spuDM) 494 if ( Calc_Forces ) call init_val(spuEDM) 495 iE = Nodes - Node 496 cE = nEq_E(iE,step=Nodes) ! we read them backwards 497 do while ( cE%exist ) 498 499 ! ******************* 500 ! * prep Sigma * 501 ! ******************* 502 call read_next_GS(ispin, ikpt, bkpt, & 503 cE, N_Elec, uGF, Elecs, & 504 nzwork, zwork, .false., forward = .false. ) 505 do iEl = 1 , N_Elec 506 call UC_expansion(cE, Elecs(iEl), nzwork, zwork, & 507 non_Eq = .true. ) 508 end do 509 510 ! ******************* 511 ! * prep GF^-1 * 512 ! ******************* 513 call prepare_invGF(cE, zwork_tri, r_pvt, pvt, & 514 N_Elec, Elecs, & 515 spH =spH , spS =spS) 516 517 ! ******************* 518 ! * prep GF * 519 ! ******************* 520 if ( .not. cE%fake ) then 521 call invert_BiasTriMat_prep(zwork_tri,GF_tri, & 522 all_nn = .true. ) 523 end if 524 525 ! ** At this point we have calculated the needed 526 ! ** information to create the Green function column 527 ! ** for all the electrodes 528 529 ! **************** 530 ! * save GF * 531 ! **************** 532 do iEl = 1 , N_Elec 533 if ( cE%fake ) cycle ! in case we implement an MPI communication solution 534 535 ! ****************** 536 ! * calc GF-column * 537 ! ****************** 538 if ( ts_A_method == TS_BTD_A_COLUMN ) then 539 call invert_BiasTriMat_rgn(GF_tri,zwork_tri, & 540 r_pvt, pvt, Elecs(iEl)%o_inD) 541 542 call GF_Gamma_GF(zwork_tri, Elecs(iEl), Elecs(iEl)%o_inD%n, & 543 calc_parts, GFGGF_size, GFGGF_work) 544#ifdef TRANSIESTA_WEIGHT_DEBUG 545 print '(a7,tr1,i3,2(tr1,f10.5),tr5,2(tr1,f10.5))', & 546 trim(Elecs(iEl)%name),iE,zwork(index(zwork_tri,28,28)),cE%e 547#endif 548 549 else 550 call dir_GF_Gamma_GF(Gf_tri, zwork_tri, r_pvt, pvt, & 551 Elecs(iEl), calc_parts) 552 end if 553 554 do iID = 1 , N_nEq_ID 555 556 if ( .not. has_cE_nEq(cE,iEl,iID) ) cycle 557 558 call c2weight_nEq(cE,iID,kw,W,imu,ZW) 559 560#ifdef TRANSIESTA_WEIGHT_DEBUG 561 print '(a20,2(tr1,i3),2(tr1,e12.5))', & 562 trim(Elecs(iEl)%name),iID,imu,W 563#endif 564 565 call add_DM( spuDM, W, spuEDM, ZW, & 566 zwork_tri, r_pvt, pvt, & 567 N_Elec, Elecs, & 568 DMidx=iID, EDMidx=imu, is_eq = .false.) 569 end do 570 end do 571 572 ! step energy-point 573 iE = iE + Nodes 574 cE = nEq_E(iE,step=Nodes) ! we read them backwards 575 end do 576 577#ifdef TRANSIESTA_TIMING 578 call timer('TS_NEQ',2) 579#endif 580 581#ifdef MPI 582 ! We need to reduce all the arrays 583 call MPI_Barrier(MPI_Comm_World,io) 584 call timer('TS_comm',1) 585 call AllReduce_SpData(spuDM, nzwork, zwork, N_nEq_id) 586 if ( Calc_Forces ) then 587 call AllReduce_SpData(spuEDM, nzwork, zwork, N_mu) 588 end if 589 call timer('TS_comm',2) 590#endif 591 592#ifdef TRANSIESTA_TIMING 593 call timer('TS_weight',1) 594#endif 595 596 ! 1. move from global UC to local SC 597 ! 2. calculate the correct contribution by applying the weight 598 ! 3. add the density to the real arrays 599 call add_k_DM(spDMneq, spuDM, N_nEq_id, & 600 spEDM, spuEDM, N_mu, & 601 n_s, sc_off, kpt, non_Eq = .true. ) 602 603 if ( TS_W_K_METHOD == TS_W_K_UNCORRELATED ) then 604 call weight_DM( N_Elec, Elecs, N_mu, mus, na_u, lasto, & 605 sp_dist, sparse_pattern, Ss, & 606 spDM, spDMneq, spEDM, n_s, sc_off, DE_NEGF) 607 608#ifdef TRANSIESTA_WEIGHT_DEBUG 609 call die('') 610#endif 611 call update_DM(sp_dist,sparse_pattern, n_nzs, & 612 DM(:,ispin), spDM, Ef=Ef, & 613 EDM=EDM(:,ispin), spEDM=spEDM, ipnt=ltsup_sc_pnt) 614 else if ( itt_last(SpKp,2) ) then ! TS_W_K_METHOD == TS_W_K_CORRELATED 615 call weight_DM( N_Elec, Elecs, N_mu, mus, na_u, lasto, & 616 sp_dist, sparse_pattern, Ss, & 617 spDM, spDMneq, spEDM, n_s, sc_off, DE_NEGF) 618 619 call update_DM(sp_dist,sparse_pattern, n_nzs, & 620 DM(:,ispin), spDM, Ef=Ef, & 621 EDM=EDM(:,ispin), spEDM=spEDM, ipnt=ltsup_sc_pnt) 622 end if 623 624#ifdef TRANSIESTA_TIMING 625 call timer('TS_weight',2) 626#endif 627 628 do iEl = 1, N_Elec 629 call reread_Gamma_Green(Elecs(iEl), uGF(iEl), NE, ispin) 630 end do 631 632 end do ! spin, k-point 633 634 call itt_destroy(SpKp) 635 636#ifdef TRANSIESTA_DEBUG 637 write(*,*) 'Completed TRANSIESTA SCF' 638#endif 639 640!*********************** 641! CLEAN UP 642!*********************** 643 644 call de_alloc(calc_parts) 645 646 call delete(zwork_tri) 647 call delete(GF_tri) 648 649 call delete(spH) 650 call delete(spS) 651 652 call delete(spuDM) 653 call delete(spuEDM) 654 655 call delete(spDM) 656 call delete(spDMneq) 657 call delete(spEDM) 658 659 ! We can safely delete the orbital distribution, it is local 660 call delete(fdist) 661 662 call clear_TriMat_inversion() 663 if ( ts_A_method == TS_BTD_A_COLUMN .and. IsVolt ) then 664 call de_alloc(GFGGF_work, routine='transiesta') 665 end if 666 call clear_mat_inversion() 667 668 call rgn_delete(pvt) 669 670 ! Nullify external pointers 671 do iEl = 1, N_Elec 672 nullify(Elecs(iEl)%Sigma) 673 end do 674 675#ifdef TRANSIESTA_DEBUG 676 call write_debug( 'POS transiesta mem' ) 677#endif 678 679 end subroutine ts_trik 680 681 682 subroutine add_DM(DM, DMfact, EDM, EDMfact, & 683 GF_tri, r, pvt, & 684 N_Elec, Elecs, & 685 DMidx,EDMidx, & 686 spS, q, & 687 is_eq) 688 689 use intrinsic_missing, only: SFIND 690 691 use class_zSpData1D 692 use class_zSpData2D 693 use class_Sparsity 694 use class_zTriMat 695 696 use m_ts_electype 697 698 ! The DM and EDM equivalent matrices 699 type(zSpData2D), intent(inout) :: DM 700 complex(dp), intent(in) :: DMfact 701 type(zSpData2D), intent(inout) :: EDM 702 complex(dp), intent(in) :: EDMfact 703 704 ! The Green function 705 type(zTriMat), intent(inout) :: GF_tri 706 type(tRgn), intent(in) :: r, pvt 707 integer, intent(in) :: N_Elec 708 type(Elec), intent(in) :: Elecs(N_Elec) 709 ! the index of the partition 710 integer, intent(in) :: DMidx 711 integer, intent(in), optional :: EDMidx 712 !< Overlap matrix setup for a k-point is needed for calculating q 713 type(zSpData1D), intent(in), optional :: spS 714 !< Charge calculated at this energy-point 715 !! 716 !! This does not contain the additional factor 1/Pi 717 real(dp), intent(inout), optional :: q 718 logical, intent(in), optional :: is_eq 719 720 ! Arrays needed for looping the sparsity 721 type(Sparsity), pointer :: s 722 integer, pointer :: l_ncol(:), l_ptr(:), l_col(:) 723 integer, pointer :: s_ncol(:), s_ptr(:), s_col(:) 724 complex(dp), pointer :: D(:,:), E(:,:) 725 complex(dp), pointer :: Gf(:) 726 complex(dp), pointer :: Sk(:) 727 integer :: s_ptr_begin, s_ptr_end, sin 728 integer :: io, ind, iu, idx, i1, i2 729 logical :: calc_q 730 logical :: hasEDM, lis_eq 731 732 lis_eq = .true. 733 if ( present(is_eq) ) lis_eq = is_eq 734 735 calc_q = present(q) .and. present(spS) 736 737 s => spar(DM) 738 call attach(s, n_col=l_ncol,list_ptr=l_ptr,list_col=l_col) 739 D => val(DM) 740 hasEDM = initialized(EDM) 741 if ( hasEDM ) E => val(EDM) 742 743 i1 = DMidx 744 i2 = i1 745 if ( present(EDMidx) ) i2 = EDMidx 746 747 Gf => val(Gf_tri) 748 749 if ( lis_eq ) then 750 751 if ( calc_q ) then 752 q = 0._dp 753 s => spar(spS) 754 Sk => val(spS) 755 call attach(s, n_col=s_ncol, list_ptr=s_ptr, list_col=s_col) 756 end if 757 758 if ( calc_q .and. hasEDM ) then 759 760!$OMP parallel do default(shared), & 761!$OMP&private(io,iu,ind,idx,s_ptr_begin,s_ptr_end,sin) 762 do iu = 1 , r%n 763 io = r%r(iu) 764 if ( l_ncol(io) /= 0 ) then 765 766 s_ptr_begin = s_ptr(io) + 1 767 s_ptr_end = s_ptr(io) + s_ncol(io) 768 769 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 770 771 idx = index(Gf_tri,iu,pvt%r(l_col(ind))) 772 773 ! Search for overlap index 774 ! spS is transposed, so we have to conjugate the 775 ! S value, then we may take the imaginary part. 776 sin = s_ptr_begin - 1 + SFIND(s_col(s_ptr_begin:s_ptr_end), l_col(ind)) 777 778 ! Since we are calculating Gf - Gf^\dagger we need a factor two 779 ! The weights are taking care of this. 780 if ( sin >= s_ptr_begin ) q = q - aimag(GF(idx) * conjg(Sk(sin))) 781 D(ind,i1) = D(ind,i1) - GF(idx) * DMfact 782 E(ind,i2) = E(ind,i2) - GF(idx) * EDMfact 783 784 end do 785 end if 786 end do 787!$OMP end parallel do 788 789 else if ( hasEDM ) then 790 791!$OMP parallel do default(shared), & 792!$OMP&private(io,iu,ind,idx) 793 do iu = 1 , r%n 794 io = r%r(iu) 795 if ( l_ncol(io) /= 0 ) then 796 797 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 798 799 idx = index(Gf_tri,iu,pvt%r(l_col(ind))) 800 801 D(ind,i1) = D(ind,i1) - GF(idx) * DMfact 802 E(ind,i2) = E(ind,i2) - GF(idx) * EDMfact 803 804 end do 805 end if 806 end do 807!$OMP end parallel do 808 809 else if ( calc_q ) then 810 811!$OMP parallel do default(shared), & 812!$OMP&private(io,iu,ind,idx,s_ptr_begin,s_ptr_end,sin) 813 do iu = 1 , r%n 814 io = r%r(iu) 815 if ( l_ncol(io) /= 0 ) then 816 817 s_ptr_begin = s_ptr(io) + 1 818 s_ptr_end = s_ptr(io) + s_ncol(io) 819 820 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 821 822 idx = index(Gf_tri,iu,pvt%r(l_col(ind))) 823 sin = s_ptr_begin - 1 + SFIND(s_col(s_ptr_begin:s_ptr_end), l_col(ind)) 824 825 if ( sin >= s_ptr_begin ) q = q - aimag(GF(idx) * conjg(Sk(sin))) 826 D(ind,i1) = D(ind,i1) - GF(idx) * DMfact 827 828 end do 829 end if 830 end do 831!$OMP end parallel do 832 833 else 834 835!$OMP parallel do default(shared), & 836!$OMP&private(io,iu,ind,idx) 837 do iu = 1 , r%n 838 io = r%r(iu) 839 if ( l_ncol(io) /= 0 ) then 840 841 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 842 843 idx = index(Gf_tri,iu,pvt%r(l_col(ind))) 844 845 D(ind,i1) = D(ind,i1) - GF(idx) * DMfact 846 847 end do 848 end if 849 end do 850!$OMP end parallel do 851 852 end if 853 854 else 855 856 if ( hasEDM ) then 857 858!$OMP parallel do default(shared), & 859!$OMP&private(io,iu,ind,idx) 860 do iu = 1 , r%n 861 io = r%r(iu) 862 if ( l_ncol(io) /= 0 ) then 863 864 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 865 866 idx = index(Gf_tri,iu,pvt%r(l_col(ind))) 867 868 D(ind,i1) = D(ind,i1) + GF(idx) * DMfact 869 E(ind,i2) = E(ind,i2) + GF(idx) * EDMfact 870 871 end do 872 end if 873 end do 874!$OMP end parallel do 875 876 else 877 878!$OMP parallel do default(shared), & 879!$OMP&private(io,iu,ind,idx) 880 do iu = 1 , r%n 881 io = r%r(iu) 882 if ( l_ncol(io) /= 0 ) then 883 884 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 885 886 idx = index(Gf_tri,iu,pvt%r(l_col(ind))) 887 888 D(ind,i1) = D(ind,i1) + GF(idx) * DMfact 889 890 end do 891 end if 892 end do 893!$OMP end parallel do 894 895 end if 896 897 end if 898 899 if ( calc_q ) q = q * 2._dp 900 901 end subroutine add_DM 902 903 ! creation of the GF^{-1}. 904 ! this routine will insert the zS-H and \Sigma_{LR} terms in the GF 905 subroutine prepare_invGF(cE, GFinv_tri, r, pvt, & 906 N_Elec, Elecs, spH, spS) 907 908 use class_Sparsity 909 use class_zSpData1D 910 use class_zTriMat 911 use m_ts_electype 912 use m_ts_tri_scat, only : insert_Self_Energies 913 use m_ts_cctype, only : ts_c_idx 914 915 ! the current energy point 916 type(ts_c_idx), intent(in) :: cE 917 type(zTriMat), intent(inout) :: GFinv_tri 918 type(tRgn), intent(in) :: r, pvt 919 integer, intent(in) :: N_Elec 920 type(Elec), intent(in) :: Elecs(N_Elec) 921 ! The Hamiltonian and overlap sparse matrices 922 type(zSpData1D), intent(inout) :: spH, spS 923 924 ! Local variables 925 complex(dp) :: Z 926 type(Sparsity), pointer :: sp 927 integer, pointer :: l_ncol(:), l_ptr(:), l_col(:) 928 complex(dp), pointer :: H(:), S(:) 929 complex(dp), pointer :: Gfinv(:) 930 integer :: io, iu, ind, idx 931 932 if ( cE%fake ) return 933 934#ifdef TRANSIESTA_TIMING 935 call timer('TS-prep',1) 936#endif 937 938 Z = cE%e 939 940 sp => spar(spH) 941 H => val (spH) 942 S => val (spS) 943 944 call attach(sp, n_col=l_ncol, list_ptr=l_ptr, list_col=l_col) 945 946 Gfinv => val(Gfinv_tri) 947 ! Initialize 948 GFinv(:) = dcmplx(0._dp,0._dp) 949 950!$OMP parallel default(shared), private(io,iu,ind,idx) 951 952 ! We will only loop in the central region 953 ! We have constructed the sparse array to only contain 954 ! values in this part... 955!$OMP do 956 do iu = 1, r%n 957 io = r%r(iu) 958 if ( l_ncol(io) /= 0 ) then 959 960 do ind = l_ptr(io) + 1 , l_ptr(io) + l_ncol(io) 961 962 ! Notice that we transpose back here... 963 ! See symmetrize_HS_kpt 964 idx = index(Gfinv_tri,pvt%r(l_col(ind)),iu) 965 966 GFinv(idx) = Z * S(ind) - H(ind) 967 end do 968 end if 969 end do 970!$OMP end do 971 972 do io = 1 , N_Elec 973 call insert_Self_Energies(Gfinv_tri, Gfinv, pvt, Elecs(io)) 974 end do 975 976!$OMP end parallel 977 978#ifdef TRANSIESTA_TIMING 979 call timer('TS-prep',2) 980#endif 981 982 end subroutine prepare_invGF 983 984end module m_ts_trik 985