← Back to branch summary

~siesta-pseudos-bases/siesta/trunk-psml

« back to all changes in this revision

Viewing changes to Src/m_ts_elec_se.F90

  • Committer: Alberto Garcia
  • Date: 2019-09-02 14:09:43 UTC
  • mfrom: (427.6.323 trunk)
  • Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/m_ts_elec_se.F90
71
71
    logical,  intent(in), optional :: non_Eq
72
72
    
73
73
    complex(dp) :: E
74
 
    integer :: nou, no, nq
 
74
    integer :: no_used, no_tot, nq
75
75
    logical :: lnon_Eq
76
76
    
77
77
    if ( cE%fake ) return
78
78
    
79
79
    call timer('ts_expand',1)
80
80
    
81
 
    nou = El%no_used
82
 
    no  = TotUsedOrbs(El)
83
 
    nq  = product(El%Bloch)
84
 
    if ( El%pre_expand > 0 .and. nq > 1 ) then
85
 
      nou = no
86
 
      nq = 1
 
81
    no_used = El%no_used
 
82
    no_tot = TotUsedOrbs(El)
 
83
    nq = El%Bloch%size()
 
84
    if ( nq > 1 ) then
 
85
      if ( El%pre_expand > 0 ) then
 
86
        no_used = no_tot
 
87
        nq = 1
 
88
      end if
87
89
    end if
88
90
    
89
91
    ! Save energy
105
107
#else
106
108
        E = cE%e
107
109
#endif
108
 
      endif
 
110
      end if
109
111
      call UC_expansion_Sigma_GammaT(E, &
110
 
          nou,no,El,nq, &
 
112
          no_used,no_tot,El,nq, &
111
113
          El%GA,El%Sigma,El%Gamma,nwork,work)
112
114
    else
113
115
      if ( El%Bulk ) then
114
 
        call UC_expansion_Sigma_Bulk(nou,no,El,nq, &
 
116
        call UC_expansion_Sigma_Bulk(no_used,no_tot,El,nq, &
115
117
            El%GA,El%Sigma,nwork,work)
116
118
      else
117
119
        if ( cE%idx(1) /= 1 ) then ! .not. CONTOUR_EQ
129
131
#endif
130
132
          endif
131
133
        end if
132
 
        call UC_expansion_Sigma(E,nou,no,El,nq, &
 
134
        call UC_expansion_Sigma(E,no_used,no_tot,El,nq, &
133
135
            El%GA,El%Sigma,nwork,work)
134
136
      end if
135
137
    end if
141
143
  subroutine UC_expansion_Sigma_Bulk(no_u,no_s,El,nq,&
142
144
      GS,Sigma,nwork,work)
143
145
    use intrinsic_missing, only : EYE
144
 
    use units, only : Pi
145
146
! ********************
146
147
! * INPUT variables  *
147
148
! ********************
237
238
        El%na_used,El%lasto_used,El%HA,El%SA,GS,nwork, &
238
239
        work(1,1,2), Sigma(1,1))
239
240
 
 
241
    ! We do not need to check for nq > 1 since
 
242
    ! the above call ensures correct handling
 
243
    
240
244
    if ( nq == 1 ) then
241
245
#ifndef TS_NOCHECKS
242
246
      if ( no_u /= no_s ) call die('no_E/=no_s')
263
267
    
264
268
    ! Do:
265
269
    ! \Sigma = Z*S - H - \Sigma_bulk
266
 
!$OMP parallel do default(shared), private(io,jo), collapse(2)
 
270
!$OMP parallel do default(shared), private(io,jo)
267
271
    do jo = 1 , no_s
268
272
      do io = 1 , no_s
269
273
        Sigma(io,jo) = work(io,jo,2) - Sigma(io,jo)
351
355
 
352
356
      ! Do:
353
357
      ! work = Z*S - H - (Z*S - H - \Sigma_bulk)
354
 
!$OMP do collapse(2)
 
358
!$OMP do
355
359
      do jo = 1 , no_s
356
360
        do io = 1 , no_s
357
361
          work(io,jo,2) = work(io,jo,2) - Sigma(io,jo)
392
396
       
393
397
      ! Do:
394
398
      ! \Sigma = Z*S - H - (Z*S - H - \Sigma_bulk)
395
 
!$OMP do collapse(2)
 
399
!$OMP do
396
400
      do jo = 1 , no_s
397
401
        do io = 1 , no_s
398
402
          Sigma(io,jo) = work(io,jo,2) - Sigma(io,jo)
457
461
! ********************
458
462
! * LOCAL variables  *
459
463
! ********************
460
 
    integer :: iuo, juo, iow
 
464
    integer :: iuo, juo, no
461
465
 
462
466
    if ( nq == 1 ) then
463
467
#ifndef TS_NOCHECKS
464
468
      if ( no_u /= no_s ) call die('no_E/=no_s')
465
469
#endif
466
470
 
467
 
      ! In case the pre-expansion is not done on H, S
468
 
      if ( El%pre_expand == 1 .and. product(El%Bloch) > 1 ) then
469
 
 
 
471
      ! nq == 1 for pre_expand > 0, hence we need to check whether HS
 
472
      ! needs to be expanded
 
473
      if ( El%pre_expand == 1 .and. El%bloch%size() > 1 ) then
 
474
        
470
475
        ! Note that this is because the interface for H and S
471
 
        iow = El%no_used
472
 
!$OMP parallel do default(shared), private(iuo,juo), collapse(2)
473
 
        do juo = 1 , iow
 
476
        no = El%no_used
 
477
!$OMP parallel do default(shared), private(iuo,juo)
 
478
        do juo = 1 , no
474
479
          do iuo = 1 , no_s
475
480
            GSE(iuo,juo) = ZEnergy * S(iuo,juo,1) - H(iuo,juo,1)
476
481
          end do
477
482
        end do
478
483
!$OMP end parallel do
479
484
        
480
 
        call update_UC_expansion_A(iow,no_s,El,nq,na_u,lasto,&
 
485
        call update_UC_expansion_A(no,no_s,El,El%bloch%size(),na_u,lasto,&
481
486
            GSE(1,1),HSE(1,1))
482
487
 
483
488
      else
484
489
 
485
 
        ! We do not need to copy over GS, as it is
486
 
        ! used correctly
487
 
        
488
 
!$OMP parallel do default(shared), private(iuo,juo), collapse(2)
 
490
!$OMP parallel do default(shared), private(iuo,juo)
489
491
        do juo = 1 , no_s
490
492
          do iuo = 1 , no_s
491
493
            !GSE(iuo,juo,1) = GS(iuo,juo,1)
495
497
!$OMP end parallel do
496
498
 
497
499
      end if
498
 
 
 
500
      
499
501
    else if ( El%repeat ) then
500
502
      call repeat(H, S, GS)
501
503
    else
502
 
      call tile(H, S, GS)
 
504
      call El%bloch%matrix_unfold_HS_G(El%bkpt_cur, no_u, H, S, GS, Zenergy, HSE, GSE)
503
505
    end if
504
506
 
505
507
  contains
506
508
 
507
509
    subroutine repeat(H, S, GS)
508
 
      complex(dp), dimension(no_u,no_u,El%Bloch(1),El%Bloch(2),El%Bloch(3)), intent(in) :: H, S, GS
 
510
      complex(dp), dimension(no_u,no_u,El%Bloch%B(1),El%Bloch%B(2),El%Bloch%B(3)), intent(in) :: H, S, GS
509
511
      integer :: B(3), i1, i2, i3
510
512
      integer :: iau,ia1,ia2,ia3
511
 
      integer :: jow,jau,ja1,ja2,ja3
 
513
      integer :: iow,jow,jau,ja1,ja2,ja3
512
514
      complex(dp) :: p(3), pZ, qPi
513
515
      real(dp) :: rPi(3), wq
514
516
      
518
520
!$OMP&  private(jow,jau,ja1,ja2,ja3,juo)
519
521
 
520
522
!$OMP workshare
 
523
      HSE(:,:) = 0._dp
521
524
      GSE(:,:) = 0._dp
522
 
      HSE(:,:) = 0._dp
523
525
!$OMP end workshare
524
526
 
525
527
      ! Save some multiplications
526
 
      B(:) = El%Bloch(:)
 
528
      B(:) = El%Bloch%B(:)
527
529
      wq = log(1._dp / real(nq,dp))
528
530
 
529
531
      do i3 = 1, B(3)
530
532
      do i2 = 1, B(2)
531
533
      do i1 = 1, B(1)
532
534
 
533
 
       rPi = 2._dp * Pi * (q_exp(El,i1,i2,i3) + El%bkpt_cur)
 
535
       rPi = 2._dp * Pi * (El%bloch%get_k(i1,i2,i3) + El%bkpt_cur)
534
536
       qPi = exp(cmplx(0._dp,rPi(1),dp))
535
537
 
536
 
!$OMP do collapse(4)
 
538
!$OMP do
537
539
       do iau = 1 , na_u
538
540
        do ia3 = 1 , B(3)
539
541
        do ia2 = 1 , B(2)
559
561
              do juo = 1 + lasto(jau-1) , lasto(jau)
560
562
               jow = jow + 1
561
563
                
 
564
               HSE(jow,iow) = HSE(jow,iow) + pZ * S(juo,iuo,i1,i2,i3) - p(1) * H(juo,iuo,i1,i2,i3)
562
565
               GSE(jow,iow) = GSE(jow,iow) + p(1) * GS(juo,iuo,i1,i2,i3)
563
 
               HSE(jow,iow) = HSE(jow,iow) + pZ * S(juo,iuo,i1,i2,i3) - p(1) * H(juo,iuo,i1,i2,i3)
564
 
                
 
566
       
565
567
              end do !juo
566
568
            end do !ja1
567
569
            end do !ja2
582
584
 
583
585
    end subroutine repeat
584
586
 
585
 
    subroutine tile(H, S, GS)
586
 
      complex(dp), dimension(no_u,no_u,El%Bloch(1),El%Bloch(2),El%Bloch(3)), intent(in) :: H, S, GS
587
 
      integer :: B(3), i1, i2, i3
588
 
      integer :: iau,ia1,ia2,ia3
589
 
      integer :: jow,jau,ja1,ja2,ja3
590
 
      complex(dp) :: p(3), pZ, qPi
591
 
      real(dp) :: rPi(3), wq
592
 
 
593
 
!$OMP parallel default(shared), private(wq,rPi,qPi,p,pZ), &
594
 
!$OMP&  private(B,i1,i2,i3), &
595
 
!$OMP&  private(iow,iau,ia1,ia2,ia3,iuo), &
596
 
!$OMP&  private(jow,jau,ja1,ja2,ja3,juo)
597
 
 
598
 
!$OMP workshare
599
 
      GSE(:,:) = 0._dp
600
 
      HSE(:,:) = 0._dp
601
 
!$OMP end workshare
602
 
 
603
 
      ! Save some multiplications
604
 
      B(:) = El%Bloch(:)
605
 
      wq = log(1._dp / real(nq,dp))
606
 
 
607
 
      do i3 = 1, B(3)
608
 
      do i2 = 1, B(2)
609
 
      do i1 = 1, B(1)
610
 
 
611
 
       rPi = 2._dp * Pi * (q_exp(El,i1,i2,i3) + El%bkpt_cur)
612
 
       qPi = exp(cmplx(0._dp,rPi(1),dp))
613
 
 
614
 
!$OMP do collapse(3)
615
 
       do ia3 = 1 , B(3)
616
 
       do ia2 = 1 , B(2)
617
 
       do ia1 = 1 , B(1)
618
 
 
619
 
        p(3) = exp(cmplx(wq,-ia1*rPi(1)-ia2*rPi(2)-ia3*rPi(3),kind=dp))
620
 
        iow = (( (ia3-1)*B(2) + (ia2-1) ) * B(1) + (ia1-1)) * no_u
621
 
        do iuo = 1, no_u
622
 
         iow = iow + 1
623
 
 
624
 
         jow = 0
625
 
         do ja3 = 1 , B(3)
626
 
         p(2) = p(3)*exp(cmplx(0._dp,ja3*rPi(3),kind=dp))
627
 
         do ja2 = 1 , B(2)
628
 
         p(1) = p(2)*exp(cmplx(0._dp,ja2*rPi(2),kind=dp))
629
 
         do ja1 = 1 , B(1)
630
 
         ! This takes one additional phase per iteration
631
 
         p(1) = p(1)*qPi
632
 
         pZ = p(1) * ZEnergy
633
 
         do juo = 1, no_u
634
 
          jow = jow + 1
635
 
                
636
 
          GSE(jow,iow) = GSE(jow,iow) + p(1) * GS(juo,iuo,i1,i2,i3)
637
 
          HSE(jow,iow) = HSE(jow,iow) + pZ * S(juo,iuo,i1,i2,i3) - p(1) * H(juo,iuo,i1,i2,i3)
638
 
                
639
 
         end do !juo
640
 
         end do !ja1
641
 
         end do !ja2
642
 
         end do !ja3
643
 
        end do !iuo
644
 
       end do !ia1
645
 
       end do !ia2
646
 
       end do !ia3
647
 
!$OMP end do
648
 
 
649
 
      end do !i1
650
 
      end do !i2
651
 
      end do !i3
652
 
 
653
 
!$OMP end parallel
654
 
 
655
 
    end subroutine tile
656
 
    
657
587
  end subroutine update_UC_expansion
658
588
 
659
589
  subroutine update_UC_expansion_A(no_u,no_s,El,nq, &
665
595
    integer,  intent(in) :: no_u, no_s
666
596
    type(Elec), intent(in) :: El
667
597
    integer,  intent(in) :: nq, na_u, lasto(0:na_u)
668
 
    complex(dp), intent(in) :: A(no_u,no_u,El%Bloch(1),El%Bloch(2),El%Bloch(3))
 
598
    complex(dp), intent(in) :: A(no_u,no_u,El%Bloch%B(1),El%Bloch%B(2),El%Bloch%B(3))
669
599
! ********************
670
600
! * OUTPUT variables *
671
601
! ********************
672
602
    complex(dp), intent(inout) :: AE(no_s,no_s)
673
603
 
674
 
    integer :: B(3), i1, i2, i3
675
 
    integer :: iau,iow,ia1,ia2,ia3,iuo
676
 
    integer :: jau,jow,ja1,ja2,ja3,juo
677
 
    complex(dp) :: p(3), qPi
678
 
    real(dp) :: rPi(3), wq
679
 
    
680
604
    if ( no_u == no_s ) then
681
605
      call die('update_UC_expansion_A: error!')
682
606
    else if ( El%repeat ) then
683
607
      call repeat()
684
608
    else
685
 
      call tile()
 
609
      call El%Bloch%matrix_unfold(El%bkpt_cur, no_u, A, AE)
686
610
    end if
687
611
    
688
612
  contains
694
618
    ! surface Green function
695
619
    
696
620
    subroutine repeat()
 
621
      integer :: B(3), i1, i2, i3
 
622
      integer :: iau,iow,ia1,ia2,ia3,iuo
 
623
      integer :: jau,jow,ja1,ja2,ja3,juo
 
624
      complex(dp) :: p(3), qPi
 
625
      real(dp) :: rPi(3), wq
697
626
      
698
627
!$OMP parallel default(shared), private(wq,rPi,qPi,p), &
699
628
!$OMP&  private(B,i1,i2,i3), &
705
634
!$OMP end workshare
706
635
 
707
636
      ! Save some multiplications
708
 
      B(:) = El%Bloch(:)
 
637
      B(:) = El%Bloch%B(:)
709
638
      wq = log(1._dp / real(nq,dp))
710
639
 
711
640
      do i3 = 1, B(3)
712
641
      do i2 = 1, B(2)
713
642
      do i1 = 1, B(1)
714
643
 
715
 
       rPi(:) = 2._dp * Pi * (q_exp(El,i1,i2,i3) + El%bkpt_cur)
716
 
       qPi = exp(cmplx(0._dp,rPi(1),kind=dp))
 
644
       rPi(:) = 2._dp * Pi * (El%bloch%get_k(i1,i2,i3) + El%bkpt_cur)
 
645
       qPi = exp(cmplx(0._dp,rPi(1),dp))
717
646
 
718
 
!$OMP do collapse(4)
 
647
!$OMP do
719
648
       do iau = 1 , na_u
720
649
        do ia3 = 1 , B(3)
721
650
        do ia2 = 1 , B(2)
760
689
!$OMP end parallel
761
690
 
762
691
    end subroutine repeat
763
 
   
764
 
    subroutine tile()
765
 
          
766
 
!$OMP parallel default(shared), private(wq,rPi,qPi,p), &
767
 
!$OMP&  private(B,i1,i2,i3), &
768
 
!$OMP&  private(iow,iau,ia1,ia2,ia3,iuo), &
769
 
!$OMP&  private(jow,jau,ja1,ja2,ja3,juo)
770
 
 
771
 
!$OMP workshare
772
 
      AE(:,:) = 0._dp
773
 
!$OMP end workshare
774
 
 
775
 
      ! Save some multiplications
776
 
      B(:) = El%Bloch(:)
777
 
      wq = log(1._dp / real(nq,dp))
778
 
 
779
 
      do i3 = 1, B(3)
780
 
      do i2 = 1, B(2)
781
 
      do i1 = 1, B(1)
782
 
 
783
 
       rPi(:) = 2._dp * Pi * (q_exp(El,i1,i2,i3) + El%bkpt_cur)
784
 
       qPi = exp(cmplx(0._dp,rPi(1),kind=dp))
785
 
 
786
 
!$OMP do collapse(3)
787
 
       do ia3 = 1 , B(3)
788
 
       do ia2 = 1 , B(2)
789
 
       do ia1 = 1 , B(1)
790
 
            
791
 
        p(3) = exp(cmplx(wq,-ia1*rPi(1)-ia2*rPi(2)-ia3*rPi(3),kind=dp))
792
 
        iow = (( (ia3-1)*B(2) + (ia2-1) ) * B(1) + (ia1-1)) * no_u
793
 
        do iuo = 1, no_u
794
 
         iow = iow + 1
795
 
           
796
 
         jow = 0
797
 
         do ja3 = 1 , B(3)
798
 
         p(2) = p(3)*exp(cmplx(0._dp,ja3*rPi(3),kind=dp))
799
 
         do ja2 = 1 , B(2)
800
 
         p(1) = p(2)*exp(cmplx(0._dp,ja2*rPi(2),kind=dp))
801
 
         do ja1 = 1 , B(1)
802
 
         p(1) = p(1)*qPi
803
 
          do juo = 1, no_u
804
 
           jow = jow + 1
805
 
               
806
 
           AE(jow,iow) = AE(jow,iow) + p(1) * A(juo,iuo,i1,i2,i3)
807
 
               
808
 
          end do !juo
809
 
         end do !ja1
810
 
         end do !ja2
811
 
         end do !ja3
812
 
        end do !iuo
813
 
       end do !ia1
814
 
       end do !ia2
815
 
       end do !ia3
816
 
!$OMP end do
817
 
 
818
 
      end do !i1
819
 
      end do !i2
820
 
      end do !i3
821
 
 
822
 
!$OMP end parallel
823
 
 
824
 
    end subroutine tile
825
 
   
 
692
 
826
693
  end subroutine update_UC_expansion_A
827
694
 
828
695
end module m_ts_elec_se