LCOV - code coverage report
Current view: top level - hybrid - trafo.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 769 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 11 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_trafo
       8             : 
       9             : CONTAINS
      10             : 
      11           0 :    SUBROUTINE waveftrafo_symm(cmt_out, z_out, cmt, l_real, z_r, z_c, bandi, ndb, &
      12             :                               nk, iop, atoms, hybrid, kpts, sym, &
      13             :                               jsp, dimension, cell, lapw)
      14             : 
      15             :       USE m_constants
      16             :       USE m_util, ONLY: modulo1
      17             :       USE m_wrapper
      18             :       USE m_types
      19             :       IMPLICIT NONE
      20             : 
      21             :       TYPE(t_dimension), INTENT(IN)   :: dimension
      22             :       TYPE(t_hybrid), INTENT(IN)      :: hybrid
      23             :       TYPE(t_sym), INTENT(IN)         :: sym
      24             :       TYPE(t_cell), INTENT(IN)        :: cell
      25             :       TYPE(t_kpts), INTENT(IN)        :: kpts
      26             :       TYPE(t_atoms), INTENT(IN)       :: atoms
      27             :       TYPE(t_lapw), INTENT(IN)        :: lapw
      28             : 
      29             : !     - scalars -
      30             :       INTEGER, INTENT(IN)      :: nk, jsp, ndb
      31             :       INTEGER, INTENT(IN)      ::  bandi, iop
      32             : 
      33             : !     - arrays -
      34             :       COMPLEX, INTENT(IN)      ::  cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
      35             :       LOGICAL, INTENT(IN)      ::  l_real
      36             :       REAL, INTENT(IN)         ::  z_r(dimension%nbasfcn, dimension%neigd)
      37             :       COMPLEX, INTENT(IN)      ::  z_c(dimension%nbasfcn, dimension%neigd)
      38             :       COMPLEX, INTENT(OUT)     ::  cmt_out(hybrid%maxlmindx, atoms%nat, ndb)
      39             :       COMPLEX, INTENT(OUT)     ::  z_out(lapw%nv(jsp), ndb)
      40             : 
      41             : !     - local -
      42             : 
      43             : !     - scalars -
      44             :       INTEGER                 ::  iatom, iatom1, iiatom, itype, igpt, igpt1, ieq, ieq1, iiop
      45             :       INTEGER                 ::  i, l, n, nn, lm0, lm1, lm2, m1, m2
      46             :       COMPLEX                 ::  cdum, tpiimg
      47             :       COMPLEX, PARAMETER       ::  img = (0.0, 1.0)
      48             : 
      49             : !     - arrays -
      50             :       INTEGER                 ::  rrot(3, 3), invrrot(3, 3)
      51             :       INTEGER                 ::  g(3), g1(3)
      52             :       REAL                    ::  tau1(3), rtaual(3), rkpt(3), rkpthlp(3), trans(3)
      53           0 :       COMPLEX                 ::  cmthlp(2*atoms%lmaxd + 1)
      54             :       LOGICAL                 ::  trs
      55             : 
      56           0 :       tpiimg = -tpi_const*img
      57             : 
      58           0 :       if (l_real) THEN
      59           0 :          rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
      60           0 :          invrrot = transpose(sym%mrot(:, :, iop))
      61           0 :          trans = sym%tau(:, iop)
      62             :       else
      63           0 :          IF (iop <= sym%nop) THEN
      64           0 :             trs = .false.
      65           0 :             rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
      66           0 :             invrrot = transpose(sym%mrot(:, :, iop))
      67           0 :             trans = sym%tau(:, iop)
      68             :          ELSE
      69           0 :             trs = .true.
      70           0 :             iiop = iop - sym%nop
      71           0 :             rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
      72           0 :             invrrot = -transpose(sym%mrot(:, :, iiop))
      73           0 :             trans = sym%tau(:, iiop)
      74             :          END IF
      75             :       endif
      76             : 
      77           0 :       rkpt = matmul(rrot, kpts%bk(:, nk))
      78           0 :       rkpthlp = rkpt
      79           0 :       rkpt = modulo1(rkpt, kpts%nkpt3)
      80           0 :       g1 = nint(rkpt - rkpthlp)
      81             : 
      82             : ! MT coefficients
      83           0 :       cmt_out = 0
      84           0 :       iatom = 0
      85           0 :       iiatom = 0
      86             : 
      87           0 :       DO itype = 1, atoms%ntype
      88           0 :          DO ieq = 1, atoms%neq(itype)
      89           0 :             iatom = iatom + 1
      90             : 
      91           0 :             iatom1 = hybrid%map(iatom, iop)
      92           0 :             tau1 = hybrid%tvec(:, iatom, iop)
      93             : 
      94           0 :             cdum = exp(tpiimg*dotprod(rkpt, tau1))
      95             : 
      96           0 :             lm0 = 0
      97           0 :             DO l = 0, atoms%lmax(itype)
      98           0 :                nn = hybrid%nindx(l, itype)
      99           0 :                DO n = 1, nn
     100           0 :                   lm1 = lm0 + n
     101           0 :                   lm2 = lm0 + n + 2*l*nn
     102           0 :                   DO i = 1, ndb
     103           0 :                      if (l_real) THEN
     104             :                         cmt_out(lm1:lm2:nn, iatom1, i) = cdum* &
     105             :      &                       matmul(cmt(bandi + i - 1, lm1:lm2:nn, iatom),&
     106           0 :      &                       sym%d_wgn(-l:l, -l:l, l, iop))
     107             :                      else
     108           0 :                         IF (trs) THEN
     109           0 :                            cmthlp(:2*l + 1) = CONJG(cmt(bandi + i - 1, lm1:lm2:nn, iatom))
     110             :                         ELSE
     111           0 :                            cmthlp(:2*l + 1) = cmt(bandi + i - 1, lm1:lm2:nn, iatom)
     112             :                         ENDIF
     113           0 :                         cmt_out(lm1:lm2:nn, iatom1, i) = cdum*matmul(cmthlp(:2*l + 1), sym%d_wgn(-l:l, -l:l, l, iop))
     114             :                      endif
     115             :                   END DO
     116             :                END DO
     117           0 :                lm0 = lm2
     118             :             END DO
     119             :          END DO
     120           0 :          iiatom = iiatom + atoms%neq(itype)
     121             :       END DO
     122             : 
     123             : ! PW coefficients
     124           0 :       z_out = 0
     125             : 
     126           0 :       DO igpt = 1, lapw%nv(jsp)
     127           0 :          g = matmul(invrrot, lapw%gvec(:, igpt, jsp) + g1)
     128             : !determine number of g
     129           0 :          igpt1 = 0
     130           0 :          DO i = 1, lapw%nv(jsp)
     131           0 :             IF (maxval(abs(g - lapw%gvec(:, i, jsp))) <= 1E-06) THEN
     132             :                igpt1 = i
     133             :                EXIT
     134             :             END IF
     135             :          END DO
     136           0 :          IF (igpt1 == 0) THEN
     137           0 :             STOP 'wavetrafo_symm: rotated G vector not found'
     138             :          END IF
     139           0 :          cdum = exp(tpiimg*dotprod(rkpt + lapw%gvec(:, igpt, jsp), trans(:)))
     140           0 :          if (l_real) THEN
     141           0 :             z_out(igpt, 1:ndb) = cdum*z_r(igpt1, bandi:bandi + ndb - 1)
     142             :          else
     143           0 :             IF (trs) THEN
     144           0 :                z_out(igpt, 1:ndb) = cdum*CONJG(z_c(igpt1, bandi:bandi + ndb - 1))
     145             :             ELSE
     146           0 :                z_out(igpt, 1:ndb) = cdum*z_c(igpt1, bandi:bandi + ndb - 1)
     147             :             END IF
     148             :          endif
     149             :       END DO
     150             : 
     151           0 :    END SUBROUTINE waveftrafo_symm
     152             : 
     153           0 :    SUBROUTINE waveftrafo_genwavf( &
     154           0 :       cmt_out, z_rout, z_cout, cmt, l_real, z_r, z_c, nk, iop, atoms, &
     155             :       hybrid, kpts, sym, jsp, dimension, nbands, &
     156             :       cell, lapw_nk, lapw_rkpt, phase)
     157             : 
     158             :       use m_juDFT
     159             :       USE m_constants
     160             :       USE m_util, ONLY: modulo1
     161             :       USE m_wrapper
     162             :       USE m_types
     163             :       IMPLICIT NONE
     164             : 
     165             :       TYPE(t_dimension), INTENT(IN)   :: dimension
     166             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
     167             :       TYPE(t_sym), INTENT(IN)   :: sym
     168             :       TYPE(t_cell), INTENT(IN)   :: cell
     169             :       TYPE(t_kpts), INTENT(IN)   :: kpts
     170             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     171             :       TYPE(t_lapw), INTENT(IN)    :: lapw_nk, lapw_rkpt
     172             : !     - scalars -
     173             :       INTEGER, INTENT(IN)      :: nk, jsp, nbands
     174             :       INTEGER, INTENT(IN)      ::  iop
     175             :       LOGICAL                 ::  phase
     176             : !     - arrays -
     177             :       COMPLEX, INTENT(IN)      ::  cmt(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     178             :       LOGICAL, INTENT(IN)      :: l_real
     179             :       REAL, INTENT(IN)         ::  z_r(dimension%nbasfcn, dimension%neigd)
     180             :       REAL, INTENT(OUT)        ::  z_rout(dimension%nbasfcn, dimension%neigd)
     181             :       COMPLEX, INTENT(IN)      ::  z_c(dimension%nbasfcn, dimension%neigd)
     182             :       COMPLEX, INTENT(OUT)     ::  z_cout(dimension%nbasfcn, dimension%neigd)
     183             : 
     184             :       COMPLEX, INTENT(OUT)    ::  cmt_out(dimension%neigd, hybrid%maxlmindx, atoms%nat)
     185             : !        - local -
     186             : 
     187             : !     - scalars -
     188             :       INTEGER                 ::  itype, iatom, iatom1, iiatom, igpt, igpt1, ieq, ieq1, iiop
     189             :       INTEGER                 ::  i, l, n, nn, lm0, lm1, lm2, m1, m2
     190             :       COMPLEX                 ::  cdum, tpiimg
     191             :       COMPLEX, PARAMETER       ::  img = (0.0, 1.0)
     192             :       LOGICAL                 ::  trs
     193             : 
     194             : !     - arrays -
     195             :       INTEGER                 ::  rrot(3, 3), invrrot(3, 3)
     196             :       INTEGER                 ::  g(3), g1(3)
     197             :       REAL                    ::  tau1(3), rkpt(3), rkpthlp(3), trans(3)
     198           0 :       COMPLEX                 ::  zhlp(dimension%nbasfcn, dimension%neigd)
     199           0 :       COMPLEX                 ::  cmthlp(2*atoms%lmaxd + 1)
     200             : 
     201           0 :       tpiimg = -tpi_const*img
     202           0 :       if (l_real) THEN
     203           0 :          rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
     204           0 :          invrrot = transpose(sym%mrot(:, :, iop))
     205           0 :          trans = sym%tau(:, iop)
     206             :       else
     207           0 :          IF (iop <= sym%nop) THEN
     208           0 :             trs = .false.
     209           0 :             rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
     210           0 :             invrrot = transpose(sym%mrot(:, :, iop))
     211           0 :             trans = sym%tau(:, iop)
     212             :          ELSE
     213             : ! in the case of SOC (l_soc=.true.)
     214             : ! time reversal symmetry is not valid anymore;
     215             : ! nsym should thus equal nop
     216           0 :             trs = .true.
     217           0 :             iiop = iop - sym%nop
     218           0 :             rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
     219           0 :             invrrot = -transpose(sym%mrot(:, :, iiop))
     220           0 :             trans = sym%tau(:, iiop)
     221             :          END IF
     222             :       endif
     223             : 
     224           0 :       rkpt = matmul(rrot, kpts%bk(:, nk))
     225           0 :       rkpthlp = rkpt
     226           0 :       rkpt = modulo1(rkpt, kpts%nkpt3)
     227           0 :       g1 = nint(rkpt - rkpthlp)
     228             : 
     229             :       ! MT coefficients
     230           0 :       cmt_out = 0
     231           0 :       iatom = 0
     232           0 :       iiatom = 0
     233             : 
     234           0 :       DO itype = 1, atoms%ntype
     235           0 :          DO ieq = 1, atoms%neq(itype)
     236           0 :             iatom = iatom + 1
     237             : 
     238           0 :             iatom1 = hybrid%map(iatom, iop)
     239           0 :             tau1 = hybrid%tvec(:, iatom, iop)
     240             : 
     241           0 :             cdum = exp(tpiimg*dotprod(rkpt, tau1))
     242             : 
     243           0 :             lm0 = 0
     244           0 :             DO l = 0, atoms%lmax(itype)
     245           0 :                nn = hybrid%nindx(l, itype)
     246           0 :                DO n = 1, nn
     247           0 :                   lm1 = lm0 + n
     248           0 :                   lm2 = lm0 + n + 2*l*nn
     249             : 
     250           0 :                   DO i = 1, nbands
     251           0 :                      if (l_real) THEN
     252             :                         cmt_out(i, lm1:lm2:nn, iatom1) = cdum*matmul(cmt(i, lm1:lm2:nn, iatom),&
     253           0 :              &                             hybrid%d_wgn2(-l:l, -l:l, l, iop))
     254             :                      else
     255           0 :                         IF (trs) THEN
     256           0 :                            cmthlp(:2*l + 1) = conjg(cmt(i, lm1:lm2:nn, iatom))
     257             :                         ELSE
     258           0 :                            cmthlp(:2*l + 1) = cmt(i, lm1:lm2:nn, iatom)
     259             :                         END IF
     260           0 :                         cmt_out(i, lm1:lm2:nn, iatom1) = cdum*matmul(cmthlp(:2*l + 1), hybrid%d_wgn2(-l:l, -l:l, l, iop))
     261             :                      endif
     262             : 
     263             :                   END DO
     264             :                END DO
     265           0 :                lm0 = lm2
     266             :             END DO
     267             :          END DO
     268           0 :          iiatom = iiatom + atoms%neq(itype)
     269             :       END DO
     270             : 
     271             :       ! PW coefficients
     272             : 
     273           0 :       zhlp = 0
     274           0 :       DO igpt = 1, lapw_rkpt%nv(jsp)
     275           0 :          g = matmul(invrrot, (/lapw_rkpt%k1(igpt, jsp), lapw_rkpt%k2(igpt, jsp), lapw_rkpt%k3(igpt, jsp)/) + g1)
     276             :          !determine number of g
     277           0 :          igpt1 = 0
     278           0 :          DO i = 1, lapw_nk%nv(jsp)
     279           0 :             IF (maxval(abs(g - (/lapw_nk%k1(i, jsp), lapw_nk%k2(i, jsp), lapw_nk%k3(i, jsp)/))) <= 1E-06) THEN
     280           0 :                igpt1 = i
     281           0 :                EXIT
     282             :             END IF
     283             :          END DO
     284           0 :          IF (igpt1 == 0) CYCLE
     285           0 :          cdum = exp(tpiimg*dotprod(rkpt + (/lapw_rkpt%k1(igpt, jsp), lapw_rkpt%k2(igpt, jsp), lapw_rkpt%k3(igpt, jsp)/), trans))
     286           0 :          if (l_real) THEN
     287           0 :             zhlp(igpt, :nbands) = cdum*z_r(igpt1, :nbands)
     288             :          else
     289           0 :             IF (trs) THEN
     290           0 :                zhlp(igpt, :nbands) = cdum*conjg(z_c(igpt1, :nbands))
     291             :             ELSE
     292           0 :                zhlp(igpt, :nbands) = cdum*z_c(igpt1, :nbands)
     293             :             END IF
     294             :          endif
     295             :       END DO
     296             : 
     297             :       ! If phase and inversion-sym. is true,
     298             :       ! define the phase such that z_out is real.
     299             : 
     300           0 :       IF (phase) THEN
     301           0 :          DO i = 1, nbands
     302           0 :             if (l_real) THEN
     303           0 :                CALL commonphase(cdum, zhlp(:, i), dimension%nbasfcn)
     304             : 
     305           0 :                IF (any(abs(aimag(zhlp(:, i)/cdum)) > 1e-8)) THEN
     306           0 :                   WRITE (*, *) maxval(abs(aimag(zhlp(:, i)/cdum)))
     307           0 :                   WRITE (*, *) zhlp
     308           0 :                   STOP 'waveftrafo1: Residual imaginary part.'
     309             :                END IF
     310           0 :                z_rout(:, i) = zhlp(:, i)/cdum
     311           0 :                cmt_out(i, :, :) = cmt_out(i, :, :)/cdum
     312             :             else
     313           0 :                z_cout(:, i) = zhlp(:, i)
     314             :             endif
     315             :          END DO
     316             :       ELSE
     317           0 :          if (l_real) THEN
     318           0 :             z_rout = zhlp
     319             :          else
     320           0 :             z_cout = zhlp
     321             :          endif
     322             :       END IF
     323             : 
     324           0 :    END SUBROUTINE waveftrafo_genwavf
     325             : 
     326             : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     327             :    ! Symmetrizes MT part of input matrix according to inversion symmetry.
     328             :    ! This is achieved by a transformation to
     329             :    !        1/sqrt(2) * ( exp(ikR) Y_lm(r-R) + (-1)**(l+m) exp(-ikR) Y_l,-m(r+R) )
     330             :    ! and                                                                                 if R /=0 or m<0
     331             :    !        i/sqrt(2) * ( exp(ikR) Y_lm(r-R) - (-1)**(l+m) exp(-ikR) Y_l,-m(r+R) ) .
     332             :    !
     333             :    !  or
     334             :    !        i*Y_l,0(r)                                                                   if R=0,m=0 and l odd
     335             :    ! These functions have the property f(-r)=f(r)* which makes the output matrix real symmetric.
     336             :    ! (Array mat is overwritten! )
     337             : 
     338           0 :    SUBROUTINE symmetrize(mat, dim1, dim2, imode, lreal, &
     339           0 :                          atoms, lcutm, maxlcutm, nindxm, sym)
     340             :       USE m_types
     341             :       IMPLICIT NONE
     342             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     343             :       TYPE(t_sym), INTENT(IN)     :: sym
     344             : 
     345             : !     - scalars -
     346             :       INTEGER, INTENT(IN)    ::  imode, dim1, dim2
     347             :       INTEGER, INTENT(IN)    :: maxlcutm
     348             :       LOGICAL, INTENT(IN)    ::  lreal
     349             : 
     350             : !     - arrays -
     351             :       INTEGER, INTENT(IN)    :: lcutm(atoms%ntype)
     352             :       INTEGER, INTENT(IN)    ::  nindxm(0:maxlcutm, atoms%ntype)
     353             :       COMPLEX, INTENT(INOUT) ::  mat(dim1, dim2)
     354             : 
     355             : !     -local scalars -
     356             :       INTEGER               ::  i, j, itype, ieq, ic, ic1, i1, i2, l, m, n, nn, ifac, ishift
     357             :       REAL                  ::  rfac, rdum, rmax
     358             :       COMPLEX               ::  img = (0.0, 1.0)
     359             : 
     360             : !     - local arrays -
     361           0 :       COMPLEX               ::  carr(max(dim1, dim2)), cfac
     362             : 
     363           0 :       rfac = sqrt(0.5)
     364           0 :       cfac = sqrt(0.5)*img
     365           0 :       ic = 0
     366           0 :       i = 0
     367             : 
     368           0 :       DO itype = 1, atoms%ntype
     369           0 :          nn = sum((/((2*l + 1)*nindxm(l, itype), l=0, lcutm(itype))/))
     370           0 :          DO ieq = 1, atoms%neq(itype)
     371           0 :             ic = ic + 1
     372           0 :             IF (atoms%invsat(ic) == 0) THEN
     373             : ! if the structure is inversion-symmetric, but the equivalent atom belongs to a different unit cell
     374             : ! invsat(atom) = 0, invsatnr(atom) = 0
     375             : ! but we need invsatnr(atom) = natom
     376             :                ic1 = ic
     377             :             ELSE
     378           0 :                ic1 = sym%invsatnr(ic)
     379             :             END IF
     380             : !ic1 = invsatnr(ic)
     381           0 :             IF (ic1 < ic) THEN
     382           0 :                i = i + nn
     383           0 :                CYCLE
     384             :             END IF
     385             : !     IF( ic1 .lt. ic ) cycle
     386           0 :             DO l = 0, lcutm(itype)
     387           0 :                ifac = -1
     388           0 :                DO m = -l, l
     389           0 :                   ifac = -ifac
     390           0 :                   ishift = (ic1 - ic)*nn - 2*m*nindxm(l, itype)
     391           0 :                   DO n = 1, nindxm(l, itype)
     392           0 :                      i = i + 1
     393           0 :                      j = i + ishift
     394           0 :                      IF (ic1 /= ic .or. m < 0) THEN
     395           0 :                         IF (iand(imode, 1) /= 0) THEN
     396           0 :                            carr(:dim2) = mat(i, :)
     397           0 :                            mat(i, :) = (carr(:dim2) + ifac*mat(j, :))*rfac
     398           0 :                            mat(j, :) = (carr(:dim2) - ifac*mat(j, :))*(-cfac)
     399             :                         END IF
     400           0 :                         IF (iand(imode, 2) /= 0) THEN
     401           0 :                            carr(:dim1) = mat(:, i)
     402           0 :                            mat(:, i) = (carr(:dim1) + ifac*mat(:, j))*rfac
     403           0 :                            mat(:, j) = (carr(:dim1) - ifac*mat(:, j))*cfac
     404             :                         END IF
     405           0 :                      ELSE IF (m == 0 .and. ifac == -1) THEN
     406           0 :                         IF (iand(imode, 1) /= 0) THEN
     407           0 :                            mat(i, :) = -img*mat(i, :)
     408             :                         END IF
     409           0 :                         IF (iand(imode, 2) /= 0) THEN
     410           0 :                            mat(:, i) = img*mat(:, i)
     411             :                         END IF
     412             :                      END IF
     413             :                   END DO
     414             :                END DO
     415             :             END DO
     416             :          END DO
     417             :       END DO
     418             : 
     419           0 :       IF (lreal) THEN
     420             : !     ! Determine common phase factor and devide by it to make the output matrix real.
     421             : !     rmax = 0
     422             : !     DO i = 1,dim1
     423             : !     DO j = 1,dim2
     424             : !     rdum = abs(real(mat(i,j)))+abs(aimag(mat(i,j)))
     425             : !     IF(rdum.gt.1e-6) THEN
     426             : !     cfac = mat(i,j)/abs(mat(i,j))
     427             : !     GOTO 1
     428             : !     ELSE IF(rdum.gt.rmax) THEN
     429             : !     cfac = mat(i,j)/abs(mat(i,j))
     430             : !     rmax = rdum
     431             : !     END IF
     432             : !     END DO
     433             : !     END DO
     434             : !     IF(1-abs(cfac)   .gt.1e-8) THEN ; mat = 0 ; RETURN ; END IF
     435             : !     1      IF(abs(1-cfac**2).gt.1e-8) mat = mat/cfac
     436             : !
     437             : !     IF(any(abs(aimag(mat)).gt.1e-8)) THEN
     438             : !     WRITE(*,*) maxval(aimag(mat))
     439             : !     STOP 'symmetrize: Residual imaginary part. Symmetrization failed.'
     440             : 
     441             : ! Determine common phase factor and divide by it to make the output matrix real.
     442           0 :          CALL commonphase(cfac, mat, dim1*dim2)
     443           0 :          mat = mat/cfac
     444           0 :          IF (any(abs(aimag(mat)) > 1e-8)) &
     445           0 :      &STOP 'symmetrize: Residual imaginary part. Symmetrization failed.'
     446             :       END IF
     447             : 
     448           0 :    END SUBROUTINE symmetrize
     449             : 
     450             : ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     451             : 
     452             :    ! Undoes symmetrization with routine symmetrize.
     453           0 :    SUBROUTINE desymmetrize(mat, dim1, dim2, imode, &
     454           0 :                            atoms, lcutm, maxlcutm, nindxm, sym)
     455             : 
     456             :       USE m_types
     457             :       IMPLICIT NONE
     458             :       TYPE(t_sym), INTENT(IN)   :: sym
     459             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     460             : 
     461             : !     - scalars -
     462             :       INTEGER, INTENT(IN)      ::  imode, dim1, dim2
     463             :       INTEGER, INTENT(IN)      :: maxlcutm
     464             : 
     465             : !     - arrays -
     466             :       INTEGER, INTENT(IN)      :: lcutm(atoms%ntype)
     467             :       INTEGER, INTENT(IN)      ::  nindxm(0:maxlcutm, atoms%ntype)
     468             :       COMPLEX, INTENT(INOUT)   ::  mat(dim1, dim2)
     469             : 
     470             : !     - local scalars -
     471             :       INTEGER                 ::  ifac, i, j, itype, ieq, ic, ic1, i1, i2, l, m, n, nn, ishift
     472             :       REAL                    ::  rfac1, rfac2
     473             :       COMPLEX                 ::  img = (0.0, 1.0)
     474             : !     - local arrays -
     475           0 :       COMPLEX                 ::  carr(max(dim1, dim2))
     476             : 
     477           0 :       rfac1 = sqrt(0.5)
     478           0 :       ic = 0
     479           0 :       i = 0
     480           0 :       DO itype = 1, atoms%ntype
     481           0 :          nn = sum((/((2*l + 1)*nindxm(l, itype), l=0, lcutm(itype))/))
     482           0 :          DO ieq = 1, atoms%neq(itype)
     483           0 :             ic = ic + 1
     484           0 :             IF (atoms%invsat(ic) == 0) THEN
     485             :                ! if the structure is inversion-symmetric, but the equivalent atom belongs to a different unit cell
     486             :                ! invsat(atom) = 0, invsatnr(atom) =0
     487             :                ! but we need invsatnr(atom) = natom
     488             :                ic1 = ic
     489             :             ELSE
     490           0 :                ic1 = sym%invsatnr(ic)
     491             :             END IF
     492             :             !ic1 = invsatnr(ic)
     493             :             !IF( ic1 .lt. ic ) cycle
     494           0 :             IF (ic1 < ic) THEN
     495           0 :                i = i + nn
     496           0 :                CYCLE
     497             :             END IF
     498           0 :             DO l = 0, lcutm(itype)
     499           0 :                ifac = -1
     500           0 :                DO m = -l, l
     501           0 :                   ifac = -ifac
     502           0 :                   rfac2 = rfac1*ifac
     503           0 :                   ishift = (ic1 - ic)*nn - 2*m*nindxm(l, itype)
     504           0 :                   DO n = 1, nindxm(l, itype)
     505           0 :                      i = i + 1
     506           0 :                      j = i + ishift
     507           0 :                      IF (ic1 /= ic .or. m < 0) THEN
     508           0 :                         IF (iand(imode, 1) /= 0) THEN
     509           0 :                            carr(:dim2) = mat(i, :)
     510           0 :                            mat(i, :) = (carr(:dim2) + img*mat(j, :))*rfac1
     511           0 :                            mat(j, :) = (carr(:dim2) - img*mat(j, :))*rfac2
     512             :                         END IF
     513           0 :                         IF (iand(imode, 2) /= 0) THEN
     514           0 :                            carr(:dim1) = mat(:, i)
     515           0 :                            mat(:, i) = (carr(:dim1) - img*mat(:, j))*rfac1
     516           0 :                            mat(:, j) = (carr(:dim1) + img*mat(:, j))*rfac2
     517             :                         END IF
     518           0 :                      ELSE IF (m == 0 .and. ifac == -1) THEN
     519           0 :                         IF (iand(imode, 1) /= 0) THEN
     520           0 :                            mat(i, :) = img*mat(i, :)
     521             :                         END IF
     522           0 :                         IF (iand(imode, 2) /= 0) THEN
     523           0 :                            mat(:, i) = -img*mat(:, i)
     524             :                         END IF
     525             :                      END IF
     526             :                   END DO
     527             :                END DO
     528             :             END DO
     529             :          END DO
     530             :       END DO
     531             : 
     532           0 :    END SUBROUTINE desymmetrize
     533             : 
     534             :    ! bra_trafo1 rotates cprod at ikpt0(<=> not irreducible k-point) to cprod at ikpt1 (bkp(ikpt0)), which is the
     535             :    ! symmetrie equivalent one
     536             :    ! isym maps ikpt0 on ikpt1
     537             : 
     538           0 :    SUBROUTINE bra_trafo2( &
     539           0 :       l_real, vecout_r, vecin_r, vecout_c, vecin_c, &
     540             :       dim, nobd, nbands, ikpt0, ikpt1, iop, sym, &
     541             :       hybrid, kpts, cell, atoms, &
     542           0 :       phase)
     543             : 
     544             :       !  ikpt0  ::  parent of ikpt1
     545             :       !  iop maps ikpt0 on ikpt1
     546             : 
     547             :       USE m_constants
     548             :       USE m_dwigner
     549             :       USE m_util
     550             :       USE m_types
     551             :       IMPLICIT NONE
     552             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
     553             :       TYPE(t_sym), INTENT(IN)   :: sym
     554             :       TYPE(t_cell), INTENT(IN)   :: cell
     555             :       TYPE(t_kpts), INTENT(IN)   :: kpts
     556             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     557             : 
     558             : !     - scalars -
     559             :       INTEGER, INTENT(IN)      ::  ikpt0, ikpt1, iop, dim, nobd, nbands
     560             : 
     561             : !     - arrays -
     562             : 
     563             :       LOGICAL, INTENT(IN)      :: l_real
     564             : 
     565             :       REAL, INTENT(IN)         ::  vecin_r(dim, nobd, nbands)
     566             :       REAL, INTENT(OUT)        ::  vecout_r(dim, nobd, nbands)
     567             :       COMPLEX, INTENT(IN)      ::  vecin_c(dim, nobd, nbands)
     568             :       COMPLEX, INTENT(OUT)     ::  vecout_c(dim, nobd, nbands)
     569             :       COMPLEX, INTENT(OUT)     ::  phase(nobd, nbands)
     570             : 
     571             : !          - local -
     572             : 
     573             : !     - scalars -
     574             :       INTEGER                 ::  nrkpt, rcent, itype, ieq, ic, l, n, i, j, nn, i1, i2, j1, j2, m1, m2, ok
     575             :       INTEGER                 ::  igptm, igptm2, igptp, icent1, iiatom, iiop, inviop
     576             :       COMPLEX                 ::  cexp, cdum
     577             :       COMPLEX, PARAMETER       ::  img = (0.0, 1.0)
     578             : !     - arrays -
     579             : 
     580             :       INTEGER                 ::  rrot(3, 3), invrot(3, 3)
     581           0 :       INTEGER                 ::  pnt(hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%nat)
     582             :       INTEGER                 ::  g(3), g1(3)
     583             :       REAL                    ::  rkpt(3), rkpthlp(3), rtaual(3), trans(3)
     584             :       REAL                    ::  arg
     585             :       COMPLEX                 ::  dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1,&
     586           0 :      &                                 -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1)
     587             : !       COMPLEX                 ::  vecin1(dim,nobd,nbands),vecout1(dim,nobd,nbands)
     588           0 :       COMPLEX, ALLOCATABLE    ::  vecin1(:, :, :), vecout1(:, :, :)
     589             : 
     590           0 :       call timestart("bra trafo")
     591             : 
     592             :       ALLOCATE (vecin1(dim, nobd, nbands), &
     593           0 :      &           vecout1(dim, nobd, nbands), stat=ok)
     594           0 :       IF (ok /= 0) &
     595           0 :      &             STOP 'bra_trafo2: error allocating vecin1 or vecout1'
     596           0 :       vecin1 = 0; vecout1 = 0
     597             : 
     598           0 :       IF (hybrid%maxlcutm1 > atoms%lmaxd) STOP 'bra_trafo2: maxlcutm > atoms%lmaxd'   ! very improbable case
     599             : 
     600             : !     transform back to unsymmetrized product basis in case of inversion symmetry
     601           0 :       if (l_real) THEN
     602           0 :          vecin1 = vecin_r
     603           0 :          DO i = 1, nbands
     604           0 :             DO j = 1, nobd
     605             :                CALL desymmetrize(vecin1(:hybrid%nbasp, j, i), hybrid%nbasp, 1, 1, &
     606           0 :                                  atoms, hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, sym)
     607             :             END DO
     608             :          END DO
     609             :       else
     610           0 :          vecin1 = vecin_c
     611             :       endif
     612             : 
     613           0 :       IF (iop <= sym%nop) THEN
     614           0 :          inviop = sym%invtab(iop)
     615           0 :          rrot = transpose(sym%mrot(:, :, sym%invtab(iop)))
     616             :          invrot = sym%mrot(:, :, sym%invtab(iop))
     617           0 :          trans = sym%tau(:, iop)
     618             : 
     619             :          dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1) &
     620           0 :             = hybrid%d_wgn2(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1, inviop)
     621             : 
     622             :       ELSE
     623           0 :          iiop = iop - sym%nop
     624           0 :          inviop = sym%invtab(iiop) + sym%nop
     625           0 :          rrot = -transpose(sym%mrot(:, :, sym%invtab(iiop)))
     626             :          invrot = sym%mrot(:, :, sym%invtab(iiop))
     627           0 :          trans = sym%tau(:, iiop)
     628             : 
     629             :          dwgn(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1) &
     630           0 :             = conjg(hybrid%d_wgn2(-hybrid%maxlcutm1:hybrid%maxlcutm1, -hybrid%maxlcutm1:hybrid%maxlcutm1, 0:hybrid%maxlcutm1, inviop))
     631             : 
     632             :       END IF
     633             : 
     634           0 :       rkpt = matmul(rrot, kpts%bkf(:, ikpt0))
     635           0 :       rkpthlp = rkpt
     636           0 :       rkpt = modulo1(rkpt, kpts%nkpt3)
     637           0 :       g = nint(rkpthlp - rkpt)
     638             : 
     639             : #ifdef CPP_DEBUG
     640             :       !test
     641             :       nrkpt = 0
     642             :       DO i = 1, kpts%nkptf
     643             :          IF (maxval(abs(rkpt - kpts%bkf(:, i))) <= 1E-06) THEN
     644             :             nrkpt = i
     645             :             EXIT
     646             :          END IF
     647             :       END DO
     648             :       IF (nrkpt /= ikpt1) THEN
     649             :          PRINT *, ikpt0, ikpt1
     650             :          PRINT *, kpts%bkf(:, ikpt1)
     651             :          PRINT *, kpts%bkf(:, ikpt0)
     652             :          PRINT *, rkpt
     653             : 
     654             :          STOP 'bra_trafo2: rotation failed'
     655             :       ENDIF
     656             : #endif
     657             : 
     658             : !     Define pointer to first mixed-basis functions (with m = -l)
     659           0 :       i = 0
     660           0 :       ic = 0
     661           0 :       DO itype = 1, atoms%ntype
     662           0 :          DO ieq = 1, atoms%neq(itype)
     663           0 :             ic = ic + 1
     664           0 :             DO l = 0, hybrid%lcutm1(itype)
     665           0 :                DO n = 1, hybrid%nindxm1(l, itype)
     666           0 :                   i = i + 1
     667           0 :                   pnt(n, l, ic) = i
     668             :                END DO
     669           0 :                i = i + hybrid%nindxm1(l, itype)*2*l
     670             :             END DO
     671             :          END DO
     672             :       END DO
     673             : 
     674             : !     Multiplication
     675             :       ! MT
     676           0 :       cexp = exp(img*tpi_const*dot_product(kpts%bkf(:, ikpt1) + g, trans(:)))
     677           0 :       ic = 0
     678           0 :       iiatom = 0
     679           0 :       DO itype = 1, atoms%ntype
     680           0 :          DO ieq = 1, atoms%neq(itype)
     681           0 :             ic = ic + 1
     682             : 
     683           0 :             rcent = hybrid%map(ic, iop)
     684             : 
     685           0 :             cdum = cexp*exp(-img*tpi_const*dot_product(g, atoms%taual(:, rcent)))
     686             : 
     687           0 :             DO l = 0, hybrid%lcutm1(itype)
     688           0 :                nn = hybrid%nindxm1(l, itype)
     689           0 :                DO n = 1, nn
     690             : 
     691           0 :                   i1 = pnt(n, l, ic)
     692           0 :                   i2 = i1 + nn*2*l
     693           0 :                   j1 = pnt(n, l, rcent)
     694           0 :                   j2 = j1 + nn*2*l
     695             : 
     696           0 :                   DO i = 1, nbands
     697           0 :                      DO j = 1, nobd
     698           0 :                         vecout1(i1:i2:nn, j, i) = cdum*matmul(vecin1(j1:j2:nn, j, i), dwgn(-l:l, -l:l, l))
     699             :                      END DO
     700             :                   END DO
     701             : 
     702             :                END DO
     703             :             END DO
     704             :          END DO
     705           0 :          iiatom = iiatom + atoms%neq(itype)
     706             :       END DO
     707             : 
     708             :       ! PW
     709           0 :       DO igptm = 1, hybrid%ngptm(ikpt0)
     710           0 :          igptp = hybrid%pgptm(igptm, ikpt0)
     711           0 :          g1 = matmul(rrot, hybrid%gptm(:, igptp)) + g
     712           0 :          igptm2 = 0
     713           0 :          DO i = 1, hybrid%ngptm(ikpt1)
     714           0 :             IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt1)))) <= 1E-06) THEN
     715             :                igptm2 = i
     716             :                EXIT
     717             :             END IF
     718             :          END DO
     719           0 :          IF (igptm2 == 0) THEN
     720           0 :             WRITE (*, *) ikpt0, ikpt1, g1
     721           0 :             WRITE (*, *) hybrid%ngptm(ikpt0), hybrid%ngptm(ikpt1)
     722           0 :             WRITE (*, *)
     723           0 :             WRITE (*, *) igptp, hybrid%gptm(:, igptp)
     724           0 :             WRITE (*, *) g
     725           0 :             WRITE (*, *) rrot
     726           0 :             WRITE (*, *) "Failed tests:", g1
     727           0 :             DO i = 1, hybrid%ngptm(ikpt1)
     728           0 :                WRITE (*, *) hybrid%gptm(:, hybrid%pgptm(i, ikpt1))
     729             :             ENDDO
     730           0 :             STOP 'bra_trafo2: G-point not found in G-point set.'
     731             :          END IF
     732           0 :          cdum = exp(img*tpi_const*dot_product(kpts%bkf(:, ikpt1) + g1, trans(:)))
     733             : 
     734           0 :          vecout1(hybrid%nbasp + igptm, :, :) = cdum*vecin1(hybrid%nbasp + igptm2, :, :)
     735             :       END DO
     736             : 
     737           0 :       DEALLOCATE (vecin1)
     738             : 
     739           0 :       if (l_real) THEN
     740           0 :          DO i = 1, nbands
     741           0 :             DO j = 1, nobd
     742             : 
     743             :                CALL symmetrize(vecout1(:, j, i), dim, 1, 1, .false., &
     744           0 :                                atoms, hybrid%lcutm1, hybrid%maxlcutm1, hybrid%nindxm1, sym)
     745             : 
     746           0 :                CALL commonphase(phase(j, i), vecout1(:, j, i), dim)
     747           0 :                vecout1(:, j, i) = vecout1(:, j, i)/phase(j, i)
     748           0 :                IF (any(abs(aimag(vecout1(:, j, i))) > 1e-8)) THEN
     749           0 :                   WRITE (*, *) vecout1(:, j, i)
     750           0 :                   STOP 'bra_trafo2: Residual imaginary part.'
     751             :                END IF
     752             : 
     753             :             END DO
     754             :          END DO
     755             :       else
     756           0 :          phase = (1.0, 0.0)
     757             :       endif
     758             : 
     759           0 :       if (l_real) THEN
     760           0 :          vecout_r = vecout1
     761             :       else
     762           0 :          vecout_c = vecout1
     763             :       endif
     764           0 :       DEALLOCATE (vecout1)
     765           0 :       call timestop("bra trafo")
     766           0 :    END SUBROUTINE bra_trafo2
     767             : 
     768             : !     This routine is not very fast at the moment.
     769           0 :    SUBROUTINE matrixtrafo(matout, matin, ikpt0, isym, lsymmetrize, atoms, &
     770           0 :                           kpts, sym, hybrid, cell, maxindxm, nindxm, nbasm, ngptmall, nbasp, &
     771           0 :                           lcutm, maxlcutm)
     772             : 
     773             :       USE m_wrapper
     774             :       USE m_dwigner
     775             :       USE m_util, ONLY: modulo1
     776             :       USE m_constants
     777             :       USE m_types
     778             :       IMPLICIT NONE
     779             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
     780             :       TYPE(t_sym), INTENT(IN)   :: sym
     781             :       TYPE(t_cell), INTENT(IN)   :: cell
     782             :       TYPE(t_kpts), INTENT(IN)   :: kpts
     783             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     784             : 
     785             :       ! - scalars -
     786             :       INTEGER, INTENT(IN)  ::  ikpt0, isym
     787             :       INTEGER, INTENT(IN)  ::  nbasp
     788             :       INTEGER, INTENT(IN)  ::  maxlcutm, maxindxm
     789             :       LOGICAL, INTENT(IN)  ::  lsymmetrize
     790             : 
     791             :       ! - arrays -
     792             :       INTEGER, INTENT(IN)  ::  nbasm(kpts%nkpt)
     793             :       INTEGER, INTENT(IN)  ::  nindxm(0:maxlcutm, atoms%ntype)
     794             :       INTEGER, INTENT(IN)  ::  ngptmall
     795             :       INTEGER, INTENT(IN)  ::  lcutm(atoms%ntype)
     796             : 
     797             :       COMPLEX, INTENT(IN)  ::  matin(nbasm(ikpt0), nbasm(ikpt0))
     798             :       COMPLEX, INTENT(OUT) ::  matout(nbasm(ikpt0), nbasm(ikpt0))
     799             :       ! - local scalars -
     800             :       INTEGER               ::  iatom, iatom1, iiatom, itype, ieq, ieq1, ic,&
     801             :      &                          l, n, i, nn, i1, i2, j1, j2
     802             :       INTEGER               ::  igptm, igptm1, igptm2, igptp, igptp1,&
     803             :      &                          ikpt1, isymi, iisym
     804             :       INTEGER               ::  m1, m2
     805             : 
     806             :       COMPLEX               ::  cexp, cdum
     807             :       COMPLEX, PARAMETER   ::  img = (0.0, 1.0)
     808             : 
     809             :       ! - local arrays -
     810           0 :       INTEGER               ::  pnt(maxindxm, 0:maxlcutm, atoms%nat),&
     811           0 :      &                          g(3), g1(3), iarr(hybrid%ngptm(ikpt0))
     812             :       INTEGER               ::  rot(3, 3), invrot(3, 3), rrot(3, 3), invrrot(3, 3)
     813             : 
     814             :       REAL                  ::  rkpt(3), rkpthlp(3), rtaual(3)
     815             :       REAL                  ::  trans(3)
     816             : 
     817           0 :       COMPLEX               ::  matin1(nbasm(ikpt0), nbasm(ikpt0))
     818           0 :       COMPLEX               ::  matout1(nbasm(ikpt0), nbasm(ikpt0))
     819             :       COMPLEX               ::  dwgn(-maxlcutm:maxlcutm,&
     820             :      &                                  -maxlcutm:maxlcutm,&
     821           0 :      &                                          0:maxlcutm)
     822             :       COMPLEX               ::  dwgninv(-maxlcutm:maxlcutm,&
     823             :      &                                  -maxlcutm:maxlcutm,&
     824           0 :      &                                          0:maxlcutm)
     825           0 :       COMPLEX               ::  carr(hybrid%ngptm(ikpt0))
     826             : 
     827             : !     Transform back to unsymmetrized product basis in case of inversion symmetry.
     828           0 :       matin1 = matin
     829           0 :       IF (lsymmetrize) THEN
     830             :          CALL desymmetrize(matin1, nbasm(ikpt0), nbasm(ikpt0), 3, &
     831           0 :                            atoms, lcutm, maxlcutm, nindxm, sym)
     832             :       END IF
     833             : 
     834           0 :       IF (isym <= sym%nop) THEN
     835           0 :          iisym = isym
     836           0 :          rot = sym%mrot(:, :, iisym)
     837           0 :          invrot = sym%mrot(:, :, sym%invtab(iisym))
     838           0 :          rrot = transpose(sym%mrot(:, :, sym%invtab(iisym)))
     839           0 :          invrrot = transpose(sym%mrot(:, :, iisym))
     840           0 :          rkpt = matmul(rrot, kpts%bk(:, ikpt0))
     841           0 :          rkpthlp = modulo1(rkpt, kpts%nkpt3)
     842           0 :          g = nint(rkpt - rkpthlp)
     843             : 
     844           0 :          CALL d_wigner(invrot, cell%bmat, maxlcutm, dwgn(:, :, 1:maxlcutm))
     845           0 :          dwgn(0, 0, 0) = 1
     846             : 
     847           0 :          DO l = 0, maxlcutm
     848           0 :             dwgn(:, :, l) = transpose(dwgn(:, :, l))
     849             :          END DO
     850             :       ELSE
     851           0 :          iisym = isym - sym%nop
     852           0 :          rot = sym%mrot(:, :, iisym)
     853           0 :          invrot = sym%mrot(:, :, sym%invtab(iisym))
     854           0 :          rrot = -transpose(sym%mrot(:, :, sym%invtab(iisym)))
     855           0 :          invrrot = -transpose(sym%mrot(:, :, iisym))
     856           0 :          rkpt = matmul(rrot, kpts%bk(:, ikpt0))
     857           0 :          rkpthlp = modulo1(rkpt, kpts%nkpt3)
     858           0 :          g = nint(rkpt - rkpthlp)
     859           0 :          matin1 = conjg(matin1)
     860             : 
     861           0 :          CALL d_wigner(invrot, cell%bmat, maxlcutm, dwgn(:, :, 1:maxlcutm))
     862           0 :          dwgn(0, 0, 0) = 1
     863             : 
     864           0 :          DO l = 0, maxlcutm
     865           0 :             dwgn(:, :, l) = transpose(dwgn(:, :, l))
     866             :          END DO
     867             : 
     868           0 :          DO l = 0, maxlcutm
     869           0 :             DO m1 = -l, l
     870           0 :                DO m2 = -l, -1
     871           0 :                   cdum = dwgn(m1, m2, l)
     872           0 :                   dwgn(m1, m2, l) = dwgn(m1, -m2, l)*(-1)**m2
     873           0 :                   dwgn(m1, -m2, l) = cdum*(-1)**m2
     874             :                END DO
     875             :             END DO
     876             :          END DO
     877             : 
     878             :       END IF
     879             :       ! determine number of rotated k-point bk(:,ikpt) -> ikpt1
     880             :       !
     881           0 :       DO i = 1, kpts%nkpt
     882           0 :          IF (maxval(abs(rkpthlp - kpts%bk(:, i))) <= 1E-06) THEN
     883             :             ikpt1 = i
     884             :             EXIT
     885             :          END IF
     886             :       END DO
     887             : 
     888           0 :       DO l = 0, maxlcutm
     889           0 :          dwgninv(-l:l, -l:l, l) = conjg(transpose(dwgn(-l:l, -l:l, l)))
     890             :       END DO
     891             : 
     892             : !     Define pointer to first mixed-basis functions (with m = -l)
     893           0 :       i = 0
     894           0 :       ic = 0
     895           0 :       DO itype = 1, atoms%ntype
     896           0 :          DO ieq = 1, atoms%neq(itype)
     897           0 :             ic = ic + 1
     898           0 :             DO l = 0, lcutm(itype)
     899           0 :                DO n = 1, nindxm(l, itype)
     900           0 :                   i = i + 1
     901           0 :                   pnt(n, l, ic) = i
     902             :                END DO
     903           0 :                i = i + nindxm(l, itype)*2*l
     904             :             END DO
     905             :          END DO
     906             :       END DO
     907             : 
     908             : !     Right-multiplication
     909             :       ! MT
     910           0 :       cexp = exp(img*tpi_const*dot_product(kpts%bk(:, ikpt1) + g, sym%tau(:, iisym)))
     911           0 :       iatom = 0
     912           0 :       iiatom = 0
     913           0 :       DO itype = 1, atoms%ntype
     914           0 :          DO ieq = 1, atoms%neq(itype)
     915           0 :             iatom = iatom + 1
     916           0 :             cdum = cexp*exp(-img*tpi_const*dot_product(g, atoms%taual(:, iatom)))
     917             : 
     918             : !           rtaual = matmul(invrot,taual(:,iatom)) + tau(:,invtab(iisym))
     919             : !           iatom1 = 0
     920             : !           DO ieq1 = 1,neq(itype)
     921             : !             IF( all(abs(modulo(rtaual-taual(:,iiatom+ieq1)+1e-12,1.0))
     922             : !      &                                               .lt. 1e-10) ) THEN ! The 1e-12 is a dirty fix.
     923             : !               iatom1 = iiatom + ieq1
     924             : !             END IF
     925             : !           END DO
     926             : !           IF( iatom1 .eq. 0 ) STOP 'matrixtrafo atom not found'
     927             : 
     928           0 :             iatom1 = hybrid%map(iatom, sym%invtab(iisym))
     929           0 :             DO l = 0, lcutm(itype)
     930           0 :                nn = nindxm(l, itype)
     931           0 :                DO n = 1, nn
     932             : 
     933           0 :                   i1 = pnt(n, l, iatom)
     934           0 :                   i2 = i1 + nn*2*l
     935           0 :                   j1 = pnt(n, l, iatom1)
     936           0 :                   j2 = j1 + nn*2*l
     937             : 
     938             :                   matout1(:, i1:i2:nn) = cdum*matmat(matin1(:, j1:j2:nn), &
     939           0 :          &                                             dwgn(-l:l, -l:l, l))
     940             : 
     941             :                END DO
     942             :             END DO
     943             :          END DO
     944           0 :          iiatom = iiatom + atoms%neq(itype)
     945             :       END DO
     946             : 
     947             :       ! PW
     948           0 :       DO igptm = 1, hybrid%ngptm(ikpt1)
     949           0 :          igptp = hybrid%pgptm(igptm, ikpt1)
     950           0 :          g1 = matmul(invrrot, hybrid%gptm(:, igptp) - g)
     951           0 :          igptm2 = 0
     952           0 :          DO i = 1, hybrid%ngptm(ikpt0)
     953           0 :             IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt0)))) <= 1E-06) THEN
     954             :                igptm2 = i
     955             :                EXIT
     956             :             END IF
     957             :          END DO
     958             : !         igptm2 = pntgptm(g1(1),g1(2),g1(3),ikpt0)
     959           0 :          IF (igptm2 == 0) STOP 'matrixtrafo: G point not found in G-point set.'
     960             : 
     961           0 :          cdum = exp(img*tpi_const*dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp), sym%tau(:, iisym)))
     962             : 
     963           0 :          matout1(:, nbasp + igptm) = cdum*matin1(:, nbasp + igptm2)
     964             : 
     965             :       END DO
     966             : 
     967             : !     Left-multiplication
     968             :       ! MT
     969           0 :       matin1 = matout1
     970           0 :       cexp = conjg(cexp)
     971           0 :       iatom = 0
     972           0 :       iiatom = 0
     973           0 :       DO itype = 1, atoms%ntype
     974           0 :          DO ieq = 1, atoms%neq(itype)
     975           0 :             iatom = iatom + 1
     976           0 :             cdum = cexp*exp(img*tpi_const*dot_product(g, atoms%taual(:, iatom)))
     977             : 
     978           0 :             iatom1 = hybrid%map(iatom, sym%invtab(iisym))
     979             : 
     980           0 :             DO l = 0, lcutm(itype)
     981           0 :                nn = nindxm(l, itype)
     982           0 :                DO n = 1, nn
     983             : 
     984           0 :                   i1 = pnt(n, l, iatom)
     985           0 :                   i2 = i1 + nn*2*l
     986           0 :                   j1 = pnt(n, l, iatom1)
     987           0 :                   j2 = j1 + nn*2*l
     988             : 
     989           0 :                   matout1(i1:i2:nn, :) = cdum*matmat(dwgninv(-l:l, -l:l, l), matin1(j1:j2:nn, :))
     990             : 
     991             :                END DO
     992             :             END DO
     993             :          END DO
     994           0 :          iiatom = iiatom + atoms%neq(itype)
     995             :       END DO
     996             : 
     997             :       ! PW
     998           0 :       DO igptm = 1, hybrid%ngptm(ikpt1)
     999           0 :          igptp = hybrid%pgptm(igptm, ikpt1)
    1000           0 :          g1 = matmul(invrrot, hybrid%gptm(:, igptp) - g)
    1001           0 :          igptm2 = 0
    1002           0 :          DO i = 1, hybrid%ngptm(ikpt0)
    1003           0 :             IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt0)))) <= 1E-06) THEN
    1004             :                igptm2 = i
    1005             :                EXIT
    1006             :             END IF
    1007             :          END DO
    1008             : !         igptm2 = pntgptm(g1(1),g1(2),g1(3),ikpt0)
    1009           0 :          IF (igptm2 == 0) STOP 'matrixtrafo: G point not found in G-point set.'
    1010           0 :          iarr(igptm) = igptm2
    1011             :          carr(igptm) = exp(-img*tpi_const* &
    1012           0 :       &              dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp), sym%tau(:, iisym)))
    1013             : !        cdum  = exp(-img * 2*pi * dot_product(bk(:,ikpt1)+gptm(:,igptp),tau(:,isym)))
    1014             : !        matout1(nbasp+igptm,:) = cdum * matin1(nbasp+igptm2,:)
    1015             :       END DO
    1016           0 :       DO i2 = 1, nbasm(ikpt1)
    1017           0 :          DO i1 = 1, hybrid%ngptm(ikpt1)
    1018           0 :             matout1(nbasp + i1, i2) = carr(i1)*matin1(nbasp + iarr(i1), i2)
    1019             :          END DO
    1020             :       END DO
    1021             : 
    1022             :       ! If inversion symmetry is applicable, symmetrize to make the values real.
    1023             : # ifdef CPP_INVERSION
    1024             :       IF (lsymmetrize) THEN
    1025             :          CALL symmetrize(matout1, nbasm(ikpt0), nbasm(ikpt0), 3, .false., &
    1026             :                          atoms, lcutm, maxlcutm, nindxm, sym)
    1027             :       END IF
    1028             : # endif
    1029           0 :       matout = matout1
    1030           0 :    END SUBROUTINE matrixtrafo
    1031             : 
    1032           0 :    SUBROUTINE matrixtrafo1( &
    1033           0 :       matout, matin, ikpt0, isym, lsymmetrize, atoms, kpts, sym, &
    1034           0 :       hybrid, cell, maxindxm, nindxm, nbasm, ngptmall, nbasp, lcutm, maxlcutm)
    1035             : 
    1036             :       USE m_wrapper
    1037             :       USE m_dwigner
    1038             :       USE m_util, ONLY: modulo1
    1039             :       USE m_constants
    1040             :       USE m_types
    1041             :       IMPLICIT NONE
    1042             : 
    1043             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
    1044             :       TYPE(t_sym), INTENT(IN)   :: sym
    1045             :       TYPE(t_cell), INTENT(IN)   :: cell
    1046             :       TYPE(t_kpts), INTENT(IN)   :: kpts
    1047             :       TYPE(t_atoms), INTENT(IN)   :: atoms
    1048             : 
    1049             :       ! - scalars -
    1050             :       INTEGER, INTENT(IN)  ::  ikpt0, isym
    1051             :       INTEGER, INTENT(IN)  ::  nbasp
    1052             :       INTEGER, INTENT(IN)  ::  maxlcutm, maxindxm
    1053             :       LOGICAL, INTENT(IN)  ::  lsymmetrize
    1054             : 
    1055             :       ! - arrays -
    1056             :       INTEGER, INTENT(IN)  ::  nbasm(kpts%nkpt)
    1057             :       INTEGER, INTENT(IN)  ::  nindxm(0:maxlcutm, atoms%ntype)
    1058             :       INTEGER, INTENT(IN)  ::  ngptmall
    1059             :       INTEGER, INTENT(IN)  ::  lcutm(atoms%ntype)
    1060             : 
    1061             :       COMPLEX, INTENT(IN)  ::  matin(nbasm(ikpt0), nbasm(ikpt0))
    1062             :       COMPLEX, INTENT(OUT) ::  matout(nbasm(ikpt0), nbasm(ikpt0))
    1063             :       ! - local scalars -
    1064             :       INTEGER               ::  iatom, iatom1, iiatom, itype, ieq, ieq1, ic, l,&
    1065             :      &                          n, i, nn, i1, i2, j1, j2
    1066             :       INTEGER               ::  igptm, igptm1, igptm2, igptp, igptp1, ikpt1,&
    1067             :      &                          isymi, iisym
    1068             :       INTEGER               ::  m1, m2
    1069             : 
    1070             :       COMPLEX               ::  cexp, cdum
    1071             :       COMPLEX, PARAMETER   ::  img = (0.0, 1.0)
    1072             : 
    1073             :       ! - local arrays -
    1074           0 :       INTEGER               ::  pnt(maxindxm, 0:maxlcutm, atoms%nat), g(3),&
    1075           0 :      &                          g1(3), iarr(hybrid%ngptm(ikpt0))
    1076             :       INTEGER               ::  rrot(3, 3)
    1077             : 
    1078             :       REAL                  ::  rkpt(3), rkpthlp(3), rtaual(3)
    1079             :       REAL                  ::  trans(3)
    1080             : 
    1081           0 :       COMPLEX               ::  matin1(nbasm(ikpt0), nbasm(ikpt0))
    1082           0 :       COMPLEX               ::  matout1(nbasm(ikpt0), nbasm(ikpt0))
    1083             :       COMPLEX               ::  dwgn(-maxlcutm:maxlcutm,&
    1084             :      &                                  -maxlcutm:maxlcutm,&
    1085           0 :      &                                          0:maxlcutm)
    1086             :       COMPLEX               ::  dwgninv(-maxlcutm:maxlcutm,&
    1087             :      &                                  -maxlcutm:maxlcutm,&
    1088           0 :      &                                          0:maxlcutm)
    1089           0 :       COMPLEX               ::  carr(hybrid%ngptm(ikpt0))
    1090             : 
    1091           0 :       IF (maxlcutm > atoms%lmaxd) STOP 'matrixtrafo1: maxlcutm .gt. atoms%lmaxd'
    1092             : 
    1093             : !     Transform back to unsymmetrized product basis in case of inversion symmetry.
    1094           0 :       matin1 = matin
    1095           0 :       IF (lsymmetrize) THEN
    1096             :          CALL desymmetrize(matin1, nbasm(ikpt0), nbasm(ikpt0), 3, &
    1097           0 :                            atoms, lcutm, maxlcutm, nindxm, sym)
    1098             :       END IF
    1099             : 
    1100           0 :       IF (isym <= sym%nop) THEN
    1101           0 :          iisym = isym
    1102           0 :          rrot = transpose(sym%mrot(:, :, sym%invtab(iisym)))
    1103             :       ELSE
    1104           0 :          iisym = isym - sym%nop
    1105           0 :          rrot = -transpose(sym%mrot(:, :, sym%invtab(iisym)))
    1106             :       END IF
    1107             : 
    1108           0 :       DO l = 0, maxlcutm
    1109             :          dwgn(-maxlcutm:maxlcutm, -maxlcutm:maxlcutm, l) = &
    1110           0 :             transpose(hybrid%d_wgn2(-maxlcutm:maxlcutm, -maxlcutm:maxlcutm, l, isym))
    1111             :       END DO
    1112             : 
    1113           0 :       rkpt = matmul(rrot, kpts%bk(:, ikpt0))
    1114           0 :       rkpthlp = modulo1(rkpt, kpts%nkpt3)
    1115           0 :       g = nint(rkpt - rkpthlp)
    1116             : 
    1117             :       ! determine number of rotated k-point bk(:,ikpt) -> ikpt1
    1118           0 :       DO i = 1, kpts%nkpt
    1119           0 :          IF (maxval(abs(rkpthlp - kpts%bk(:, i))) <= 1E-06) THEN
    1120             :             ikpt1 = i
    1121             :             EXIT
    1122             :          END IF
    1123             :       END DO
    1124             : 
    1125             : !     Define pointer to first mixed-basis functions (with m = -l)
    1126           0 :       i = 0
    1127           0 :       ic = 0
    1128           0 :       DO itype = 1, atoms%ntype
    1129           0 :          DO ieq = 1, atoms%neq(itype)
    1130           0 :             ic = ic + 1
    1131           0 :             DO l = 0, lcutm(itype)
    1132           0 :                DO n = 1, nindxm(l, itype)
    1133           0 :                   i = i + 1
    1134           0 :                   pnt(n, l, ic) = i
    1135             :                END DO
    1136           0 :                i = i + nindxm(l, itype)*2*l
    1137             :             END DO
    1138             :          END DO
    1139             :       END DO
    1140             : 
    1141             : !     Right-multiplication
    1142             :       ! MT
    1143           0 :       cexp = exp(-img*tpi_const*dot_product(kpts%bk(:, ikpt1) + g, sym%tau(:, iisym)))
    1144           0 :       iatom = 0
    1145           0 :       iiatom = 0
    1146           0 :       DO itype = 1, atoms%ntype
    1147           0 :          DO ieq = 1, atoms%neq(itype)
    1148           0 :             iatom = iatom + 1
    1149             : 
    1150           0 :             iatom1 = hybrid%map(iatom, iisym)
    1151           0 :             cdum = cexp*exp(img*tpi_const*dot_product(g, atoms%taual(:, iatom1)))
    1152             : 
    1153           0 :             DO l = 0, lcutm(itype)
    1154           0 :                nn = nindxm(l, itype)
    1155           0 :                DO n = 1, nn
    1156             : 
    1157           0 :                   i1 = pnt(n, l, iatom)
    1158           0 :                   i2 = i1 + nn*2*l
    1159           0 :                   j1 = pnt(n, l, iatom1)
    1160           0 :                   j2 = j1 + nn*2*l
    1161             : 
    1162           0 :                   matout1(:, i1:i2:nn) = cdum*matmat(matin1(:, j1:j2:nn), dwgn(-l:l, -l:l, l))
    1163             : 
    1164             :                END DO
    1165             :             END DO
    1166             :          END DO
    1167           0 :          iiatom = iiatom + atoms%neq(itype)
    1168             :       END DO
    1169             : 
    1170             :       ! PW
    1171           0 :       DO igptm = 1, hybrid%ngptm(ikpt0)
    1172           0 :          igptp = hybrid%pgptm(igptm, ikpt0)
    1173           0 :          g1 = matmul(rrot, hybrid%gptm(:, igptp)) + g
    1174           0 :          igptm1 = 0
    1175           0 :          DO i = 1, hybrid%ngptm(ikpt1)
    1176           0 :             IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt1)))) <= 1E-06) THEN
    1177             :                igptm1 = i
    1178             :                igptp1 = hybrid%pgptm(i, ikpt1)
    1179             :                EXIT
    1180             :             END IF
    1181             :          END DO
    1182           0 :          IF (igptm1 == 0) STOP 'matrixtrafo1: G point not found in G-point set.'
    1183             : 
    1184           0 :          cdum = exp(-img*tpi_const*dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp1), sym%tau(:, iisym)))
    1185             : 
    1186           0 :          matout1(:, nbasp + igptm) = cdum*matin1(:, nbasp + igptm1)
    1187             : 
    1188             :       END DO
    1189             : 
    1190             : !     Left-multiplication
    1191             :       ! MT
    1192             : 
    1193           0 :       DO l = 0, maxlcutm
    1194           0 :          dwgninv(-l:l, -l:l, l) = conjg(transpose(dwgn(-l:l, -l:l, l)))
    1195             :       END DO
    1196             : 
    1197           0 :       matin1 = matout1
    1198           0 :       cexp = conjg(cexp)
    1199             : 
    1200           0 :       iatom = 0
    1201           0 :       iiatom = 0
    1202           0 :       DO itype = 1, atoms%ntype
    1203           0 :          DO ieq = 1, atoms%neq(itype)
    1204           0 :             iatom = iatom + 1
    1205             : 
    1206           0 :             iatom1 = hybrid%map(iatom, iisym)
    1207           0 :             cdum = cexp*exp(-img*tpi_const*dot_product(g, atoms%taual(:, iatom1)))
    1208             : 
    1209           0 :             DO l = 0, lcutm(itype)
    1210           0 :                nn = nindxm(l, itype)
    1211           0 :                DO n = 1, nn
    1212             : 
    1213           0 :                   i1 = pnt(n, l, iatom)
    1214           0 :                   i2 = i1 + nn*2*l
    1215           0 :                   j1 = pnt(n, l, iatom1)
    1216           0 :                   j2 = j1 + nn*2*l
    1217             : 
    1218           0 :                   matout1(i1:i2:nn, :) = cdum*matmat(dwgninv(-l:l, -l:l, l), matin1(j1:j2:nn, :))
    1219             : 
    1220             :                END DO
    1221             :             END DO
    1222             :          END DO
    1223           0 :          iiatom = iiatom + atoms%neq(itype)
    1224             :       END DO
    1225             : 
    1226             :       ! PW
    1227           0 :       DO igptm = 1, hybrid%ngptm(ikpt0)
    1228           0 :          igptp = hybrid%pgptm(igptm, ikpt0)
    1229           0 :          g1 = matmul(rrot, hybrid%gptm(:, igptp)) + g
    1230           0 :          igptm1 = 0
    1231           0 :          DO i = 1, hybrid%ngptm(ikpt1)
    1232           0 :             IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt1)))) <= 1E-06) THEN
    1233             :                igptm1 = i
    1234             :                igptp1 = hybrid%pgptm(i, ikpt1)
    1235             :                EXIT
    1236             :             END IF
    1237             :          END DO
    1238           0 :          IF (igptm1 == 0) STOP 'matrixtrafo1: G point not found in G-point set.'
    1239           0 :          iarr(igptm) = igptm1
    1240           0 :          carr(igptm) = exp(img*tpi_const*dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp1), sym%tau(:, iisym)))
    1241             :       END DO
    1242           0 :       DO i2 = 1, nbasm(ikpt0)
    1243           0 :          DO i1 = 1, hybrid%ngptm(ikpt0)
    1244           0 :             matout1(nbasp + i1, i2) = carr(i1)*matin1(nbasp + iarr(i1), i2)
    1245             :          END DO
    1246             :       END DO
    1247             : 
    1248             :       ! If inversion symmetry is applicable, symmetrize to make the values real.
    1249           0 :       IF (lsymmetrize) THEN
    1250             :          CALL symmetrize(matout1, nbasm(ikpt0), nbasm(ikpt0), 3, .false., &
    1251           0 :                          atoms, lcutm, maxlcutm, nindxm, sym)
    1252             :       END IF
    1253             : 
    1254           0 :       IF (isym <= sym%nop) THEN
    1255           0 :          matout = matout1
    1256             :       ELSE
    1257           0 :          matout = conjg(matout1)
    1258             :       END IF
    1259             : 
    1260           0 :    END SUBROUTINE matrixtrafo1
    1261             : 
    1262           0 :    SUBROUTINE ket_trafo(&
    1263           0 :   &        vecout, vecin, ikpt0, isym, lreal, lsymmetrize,&
    1264             :   &        atoms, kpts, sym,&
    1265           0 :   &        hybrid, cell, maxindxm, nindxm, nbasm,&
    1266           0 :   &        ngptmall, nbasp, lcutm, maxlcutm)
    1267             : 
    1268             :       USE m_constants
    1269             :       USE m_util, ONLY: modulo1
    1270             :       USE m_dwigner
    1271             :       USE m_types
    1272             :       IMPLICIT NONE
    1273             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
    1274             :       TYPE(t_sym), INTENT(IN)   :: sym
    1275             :       TYPE(t_cell), INTENT(IN)   :: cell
    1276             :       TYPE(t_kpts), INTENT(IN)   :: kpts
    1277             :       TYPE(t_atoms), INTENT(IN)   :: atoms
    1278             : 
    1279             :       ! -scalars -
    1280             :       INTEGER, INTENT(IN)  ::  ikpt0, isym
    1281             :       INTEGER, INTENT(IN)  ::  nbasp
    1282             :       INTEGER, INTENT(IN)  ::  maxlcutm, maxindxm
    1283             :       LOGICAL, INTENT(IN)  ::  lreal, lsymmetrize
    1284             : 
    1285             :       ! - arrays -
    1286             : 
    1287             :       INTEGER, INTENT(IN)  ::  nbasm(kpts%nkpt)
    1288             :       INTEGER, INTENT(IN)  ::  nindxm(0:maxlcutm, atoms%ntype)
    1289             :       INTEGER, INTENT(IN)  ::  ngptmall
    1290             :       INTEGER, INTENT(IN)  ::  lcutm(atoms%ntype)
    1291             : 
    1292             :       COMPLEX, INTENT(IN)  ::  vecin(nbasm(ikpt0))
    1293             :       COMPLEX, INTENT(OUT) ::  vecout(nbasm(ikpt0))
    1294             : 
    1295             :       ! -local scalars -
    1296             :       INTEGER               ::  iatom, iatom1, iiatom, itype, ieq, ieq1, ic, l,&
    1297             :      &                          n, i, nn, i1, i2, j1, j2
    1298             :       INTEGER               ::  igptm, igptm1, igptm2, igptp, igptp1, ikpt1,&
    1299             :      &                          isymi, iisym
    1300             :       INTEGER               ::  m1, m2
    1301             : 
    1302             :       COMPLEX               ::  cexp, cdum
    1303             :       COMPLEX, PARAMETER   ::  img = (0.0, 1.0)
    1304             : 
    1305             :       ! - local arrays -
    1306           0 :       INTEGER               ::  pnt(maxindxm, 0:maxlcutm, atoms%nat), g(3), g1(3)
    1307             :       INTEGER               ::  rot(3, 3), invrot(3, 3),&
    1308             :      &                          rrot(3, 3), invrrot(3, 3)
    1309             : 
    1310             :       REAL                  ::  rkpt(3), rkpthlp(3), rtaual(3)
    1311             : 
    1312           0 :       COMPLEX               ::  vecin1(nbasm(ikpt0))
    1313           0 :       COMPLEX               ::  vecout1(nbasm(ikpt0))
    1314             :       COMPLEX               ::  dwgn(-maxlcutm:maxlcutm,&
    1315             :      &                                  -maxlcutm:maxlcutm,&
    1316           0 :      &                                          0:maxlcutm)
    1317             :       COMPLEX               ::  dwgninv(-maxlcutm:maxlcutm,&
    1318             :      &                                  -maxlcutm:maxlcutm,&
    1319             :      &                                          0:maxlcutm)
    1320             : 
    1321             : !     Transform back to unsymmetrized product basis in case of inversion symmetry.
    1322           0 :       vecin1 = vecin!(:nbasp)
    1323           0 :       IF (lsymmetrize) THEN
    1324             :          CALL desymmetrize(vecin1(:nbasp), 1, nbasp, 2, &
    1325           0 :                            atoms, lcutm, maxlcutm, nindxm, sym)
    1326             :       END IF
    1327             : 
    1328           0 :       IF (isym <= sym%nop) THEN
    1329           0 :          iisym = isym
    1330           0 :          rot = sym%mrot(:, :, iisym)
    1331           0 :          invrot = sym%mrot(:, :, sym%invtab(iisym))
    1332           0 :          rrot = transpose(sym%mrot(:, :, sym%invtab(iisym)))
    1333           0 :          invrrot = transpose(sym%mrot(:, :, iisym))
    1334           0 :          rkpt = matmul(rrot, kpts%bk(:, ikpt0))
    1335           0 :          rkpthlp = modulo1(rkpt, kpts%nkpt3)
    1336           0 :          g = nint(rkpt - rkpthlp)
    1337             : 
    1338           0 :          CALL d_wigner(invrot, cell%bmat, maxlcutm, dwgn(:, :, 1:maxlcutm))
    1339           0 :          dwgn(0, 0, 0) = 1
    1340             : 
    1341           0 :          DO l = 0, maxlcutm
    1342           0 :             dwgn(:, :, l) = transpose(dwgn(:, :, l))
    1343             :          END DO
    1344             :       ELSE
    1345           0 :          iisym = isym - sym%nop
    1346           0 :          rot = sym%mrot(:, :, iisym)
    1347           0 :          rrot = -transpose(sym%mrot(:, :, sym%invtab(iisym)))
    1348           0 :          invrot = sym%mrot(:, :, sym%invtab(iisym))
    1349           0 :          invrrot = -transpose(sym%mrot(:, :, iisym))
    1350           0 :          rkpt = matmul(rrot, kpts%bk(:, ikpt0))
    1351           0 :          rkpthlp = modulo1(rkpt, kpts%nkpt3)
    1352           0 :          g = nint(rkpt - rkpthlp)
    1353           0 :          vecin1 = conjg(vecin1)
    1354             : 
    1355           0 :          CALL d_wigner(invrot, cell%bmat, maxlcutm, dwgn(:, :, 1:maxlcutm))
    1356           0 :          dwgn(0, 0, 0) = 1
    1357             : 
    1358           0 :          DO l = 0, maxlcutm
    1359           0 :             dwgn(:, :, l) = transpose(dwgn(:, :, l))
    1360             :          END DO
    1361             : 
    1362           0 :          DO l = 0, maxlcutm
    1363           0 :             DO m1 = -l, l
    1364           0 :                DO m2 = -l, -1
    1365           0 :                   cdum = dwgn(m1, m2, l)
    1366           0 :                   dwgn(m1, m2, l) = dwgn(m1, -m2, l)*(-1)**m2
    1367           0 :                   dwgn(m1, -m2, l) = cdum*(-1)**m2
    1368             :                END DO
    1369             :             END DO
    1370             :          END DO
    1371             :       END IF
    1372             : 
    1373             :       ! determine number of rotated k-point bk(:,ikpt) -> ikpt1
    1374             :       !
    1375           0 :       DO i = 1, kpts%nkpt
    1376           0 :          IF (maxval(abs(rkpthlp - kpts%bk(:, i))) <= 1E-06) THEN
    1377             :             ikpt1 = i
    1378             :             EXIT
    1379             :          END IF
    1380             :       END DO
    1381             : 
    1382             : !       DO l = 0,maxlcutm
    1383             : !         dwgninv(-l:l,-l:l,l) = conjg(transpose(dwgn(-l:l,-l:l,l)))
    1384             : !       END DO
    1385             : 
    1386             : !     Define pointer to first mixed-basis functions (with m = -l)
    1387           0 :       i = 0
    1388           0 :       ic = 0
    1389           0 :       DO itype = 1, atoms%ntype
    1390           0 :          DO ieq = 1, atoms%neq(itype)
    1391           0 :             ic = ic + 1
    1392           0 :             DO l = 0, lcutm(itype)
    1393           0 :                DO n = 1, nindxm(l, itype)
    1394           0 :                   i = i + 1
    1395           0 :                   pnt(n, l, ic) = i
    1396             :                END DO
    1397           0 :                i = i + nindxm(l, itype)*2*l
    1398             :             END DO
    1399             :          END DO
    1400             :       END DO
    1401             : 
    1402             : !     Multiplication
    1403             :       ! MT
    1404           0 :       cexp = exp(img*tpi_const*dot_product(kpts%bk(:, ikpt1) + g, sym%tau(:, iisym)))
    1405           0 :       iatom = 0
    1406           0 :       iiatom = 0
    1407           0 :       DO itype = 1, atoms%ntype
    1408           0 :          DO ieq = 1, atoms%neq(itype)
    1409           0 :             iatom = iatom + 1
    1410           0 :             cdum = cexp*exp(-img*tpi_const*dot_product(g, atoms%taual(:, iatom)))
    1411             : 
    1412           0 :             rtaual = matmul(rot, atoms%taual(:, iatom)) + sym%tau(:, iisym)
    1413           0 :             iatom1 = 0
    1414           0 :             DO ieq1 = 1, atoms%neq(itype)
    1415           0 :                IF (all(abs(modulo(rtaual - atoms%taual(:, iiatom + ieq1) + 1e-12, 1.0))&
    1416           0 :         &                                               < 1e-10)) THEN ! The 1e-12 is a dirty fix.
    1417           0 :                   iatom1 = iiatom + ieq1
    1418             :                END IF
    1419             :             END DO
    1420           0 :             IF (iatom1 == 0) STOP 'ket_trafo: rotated atom not found'
    1421             : 
    1422           0 :             DO l = 0, lcutm(itype)
    1423           0 :                nn = nindxm(l, itype)
    1424           0 :                DO n = 1, nn
    1425             : 
    1426           0 :                   i1 = pnt(n, l, iatom)
    1427           0 :                   i2 = i1 + nn*2*l
    1428           0 :                   j1 = pnt(n, l, iatom1)
    1429           0 :                   j2 = j1 + nn*2*l
    1430             : 
    1431             :                   vecout1(i1:i2:nn) = cdum*matmul(vecin1(j1:j2:nn),&
    1432           0 :          &                                           dwgn(-l:l, -l:l, l))
    1433             : 
    1434             :                END DO
    1435             :             END DO
    1436             :          END DO
    1437           0 :          iiatom = iiatom + atoms%neq(itype)
    1438             :       ENDDO
    1439             : 
    1440             :       ! PW
    1441           0 :       DO igptm = 1, hybrid%ngptm(ikpt1)
    1442           0 :          igptp = hybrid%pgptm(igptm, ikpt1)
    1443           0 :          g1 = matmul(invrrot, hybrid%gptm(:, igptp) - g)
    1444           0 :          igptm2 = 0
    1445           0 :          DO i = 1, hybrid%ngptm(ikpt0)
    1446           0 :             IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt0)))) <= 1E-06) THEN
    1447             :                igptm2 = i
    1448             :                EXIT
    1449             :             END IF
    1450             :          END DO
    1451           0 :          IF (igptm2 == 0) &
    1452           0 :       &               STOP 'ket_trafo: G point not found in G-point set.'
    1453             : 
    1454             :          cdum = exp(img*tpi_const* &
    1455           0 :       &              dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp), sym%tau(:, iisym)))
    1456             : 
    1457           0 :          vecout1(nbasp + igptm) = cdum*vecin1(nbasp + igptm2)
    1458             : 
    1459             :       END DO
    1460             : 
    1461             :       ! If inversion symmetry is applicable, define the phase of vecout and symmetrize to make the values real.
    1462           0 :       IF (lsymmetrize) CALL symmetrize(vecout1, 1, nbasm(ikpt1), 2, lreal, &
    1463           0 :                                        atoms, lcutm, maxlcutm, nindxm, sym)
    1464           0 :       vecout = vecout1
    1465             : 
    1466           0 :    END SUBROUTINE ket_trafo
    1467             : 
    1468           0 :    SUBROUTINE ket_trafo1(vecout, vecin, ikpt0, isym, lreal, lsymmetrize, &
    1469             :                          atoms, kpts, sym, hybrid, &
    1470           0 :                          cell, maxindxm, nindxm, nbasm, ngptmall, nbasp, lcutm, maxlcutm)
    1471             : 
    1472             :       USE m_constants
    1473             :       USE m_util, ONLY: modulo1
    1474             :       USE m_dwigner
    1475             :       USE m_types
    1476             :       IMPLICIT NONE
    1477             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
    1478             :       TYPE(t_sym), INTENT(IN)   :: sym
    1479             :       TYPE(t_cell), INTENT(IN)   :: cell
    1480             :       TYPE(t_kpts), INTENT(IN)   :: kpts
    1481             :       TYPE(t_atoms), INTENT(IN)   :: atoms
    1482             : 
    1483             :       ! -scalars -
    1484             :       INTEGER, INTENT(IN)  ::  ikpt0, isym
    1485             :       INTEGER, INTENT(IN)  ::  nbasp
    1486             :       INTEGER, INTENT(IN)  ::  maxlcutm, maxindxm
    1487             :       LOGICAL, INTENT(IN)  ::  lreal, lsymmetrize
    1488             : 
    1489             :       ! - arrays -
    1490             : 
    1491             :       INTEGER, INTENT(IN)  ::  nbasm(kpts%nkpt)
    1492             :       INTEGER, INTENT(IN)  ::  nindxm(0:maxlcutm, atoms%ntype)
    1493             :       INTEGER, INTENT(IN)  ::  ngptmall
    1494             :       INTEGER, INTENT(IN)  ::  lcutm(atoms%ntype)
    1495             : 
    1496             :       COMPLEX, INTENT(IN)  ::  vecin(nbasm(ikpt0))
    1497             :       COMPLEX, INTENT(OUT) ::  vecout(nbasm(ikpt0))
    1498             : 
    1499             :       ! -local scalars -
    1500             :       INTEGER               ::  iatom, iatom1, iiatom, itype, ieq, ieq1, ic, l,&
    1501             :      &                          n, i, nn, i1, i2, j1, j2
    1502             :       INTEGER               ::  igptm, igptm1, igptm2, igptp, igptp1, ikpt1,&
    1503             :      &                          isymi, iisym
    1504             :       INTEGER               ::  m1, m2
    1505             : 
    1506             :       COMPLEX               ::  cexp, cdum
    1507             :       COMPLEX, PARAMETER   ::  img = (0.0, 1.0)
    1508             : 
    1509             :       ! - local arrays -
    1510           0 :       INTEGER               ::  pnt(maxindxm, 0:maxlcutm, atoms%nat), g(3), g1(3)
    1511             :       INTEGER               ::  rrot(3, 3)
    1512             : 
    1513             :       REAL                  ::  rkpt(3), rkpthlp(3), rtaual(3)
    1514             : 
    1515           0 :       COMPLEX               ::  vecin1(nbasm(ikpt0))
    1516           0 :       COMPLEX               ::  vecout1(nbasm(ikpt0))
    1517             :       COMPLEX               ::  dwgn(-maxlcutm:maxlcutm,&
    1518             :      &                               -maxlcutm:maxlcutm,&
    1519           0 :      &                                       0:maxlcutm)
    1520             : 
    1521           0 :       IF (maxlcutm > atoms%lmaxd) STOP 'kettrafo1: maxlcutm > atoms%lmaxd'
    1522             : 
    1523             : !     Transform back to unsymmetrized product basis in case of inversion symmetry.
    1524           0 :       vecin1 = vecin!(:nbasp)
    1525           0 :       IF (lsymmetrize) CALL desymmetrize(vecin1(:nbasp), 1, nbasp, 2, &
    1526           0 :                                          atoms, lcutm, maxlcutm, nindxm, sym)
    1527             : 
    1528           0 :       IF (isym <= sym%nop) THEN
    1529           0 :          iisym = isym
    1530           0 :          rrot = transpose(sym%mrot(:, :, sym%invtab(iisym)))
    1531             :       ELSE
    1532           0 :          iisym = isym - sym%nop
    1533           0 :          rrot = -transpose(sym%mrot(:, :, sym%invtab(iisym)))
    1534             :       END IF
    1535             : 
    1536           0 :       rkpt = matmul(rrot, kpts%bk(:, ikpt0))
    1537           0 :       rkpthlp = modulo1(rkpt, kpts%nkpt3)
    1538           0 :       g = nint(rkpt - rkpthlp)
    1539             : 
    1540           0 :       DO l = 0, maxlcutm
    1541             :          dwgn(-maxlcutm:maxlcutm, -maxlcutm:maxlcutm, l) = &
    1542           0 :             transpose(hybrid%d_wgn2(-maxlcutm:maxlcutm, -maxlcutm:maxlcutm, l, isym))
    1543             :       END DO
    1544             : 
    1545             :       !
    1546             :       ! determine number of rotated k-point bk(:,ikpt) -> ikpt1
    1547             :       !
    1548           0 :       DO i = 1, kpts%nkpt
    1549           0 :          IF (maxval(abs(rkpthlp - kpts%bk(:, i))) <= 1E-06) THEN
    1550             :             ikpt1 = i
    1551             :             EXIT
    1552             :          END IF
    1553             :       END DO
    1554             : 
    1555             : !     Define pointer to first mixed-basis functions (with m = -l)
    1556           0 :       i = 0
    1557           0 :       ic = 0
    1558           0 :       DO itype = 1, atoms%ntype
    1559           0 :          DO ieq = 1, atoms%neq(itype)
    1560           0 :             ic = ic + 1
    1561           0 :             DO l = 0, lcutm(itype)
    1562           0 :                DO n = 1, nindxm(l, itype)
    1563           0 :                   i = i + 1
    1564           0 :                   pnt(n, l, ic) = i
    1565             :                END DO
    1566           0 :                i = i + nindxm(l, itype)*2*l
    1567             :             END DO
    1568             :          END DO
    1569             :       END DO
    1570             : 
    1571             : !     Multiplication
    1572             :       ! MT
    1573           0 :       cexp = exp(-img*tpi_const*dot_product(kpts%bk(:, ikpt1) + g, sym%tau(:, iisym)))
    1574           0 :       iatom = 0
    1575           0 :       iiatom = 0
    1576           0 :       DO itype = 1, atoms%ntype
    1577           0 :          DO ieq = 1, atoms%neq(itype)
    1578           0 :             iatom = iatom + 1
    1579             : 
    1580           0 :             iatom1 = hybrid%map(iatom, iisym)
    1581           0 :             cdum = cexp*exp(img*tpi_const*dot_product(g, atoms%taual(:, iatom1)))
    1582             : 
    1583           0 :             DO l = 0, lcutm(itype)
    1584           0 :                nn = nindxm(l, itype)
    1585           0 :                DO n = 1, nn
    1586             : 
    1587           0 :                   i1 = pnt(n, l, iatom)
    1588           0 :                   i2 = i1 + nn*2*l
    1589           0 :                   j1 = pnt(n, l, iatom1)
    1590           0 :                   j2 = j1 + nn*2*l
    1591             : 
    1592           0 :                   vecout1(i1:i2:nn) = cdum*matmul(vecin1(j1:j2:nn), dwgn(-l:l, -l:l, l))
    1593             : 
    1594             :                END DO
    1595             :             END DO
    1596             :          END DO
    1597           0 :          iiatom = iiatom + atoms%neq(itype)
    1598             :       ENDDO
    1599             : 
    1600             :       ! PW
    1601           0 :       DO igptm = 1, hybrid%ngptm(ikpt0)
    1602           0 :          igptp = hybrid%pgptm(igptm, ikpt0)
    1603           0 :          g1 = matmul(rrot, hybrid%gptm(:, igptp)) + g
    1604           0 :          igptm1 = 0
    1605           0 :          DO i = 1, hybrid%ngptm(ikpt1)
    1606           0 :             IF (maxval(abs(g1 - hybrid%gptm(:, hybrid%pgptm(i, ikpt1)))) <= 1E-06) THEN
    1607             :                igptm1 = i
    1608             :                igptp1 = hybrid%pgptm(i, ikpt1)
    1609             :                EXIT
    1610             :             END IF
    1611             :          END DO
    1612           0 :          IF (igptm1 == 0) STOP 'ket_trafo: G point not found in G-point set.'
    1613             : 
    1614           0 :          cdum = exp(-img*tpi_const*dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp1), sym%tau(:, iisym)))
    1615             : 
    1616           0 :          vecout1(nbasp + igptm) = cdum*vecin1(nbasp + igptm1)
    1617             : 
    1618             :       END DO
    1619             : 
    1620             :       ! If inversion symmetry is applicable, define the phase of vecout and symmetrize to make the values real.
    1621           0 :       IF (lsymmetrize) CALL symmetrize(vecout1, 1, nbasm(ikpt1), 2, lreal, &
    1622           0 :                                        atoms, lcutm, maxlcutm, nindxm, sym)
    1623           0 :       IF (isym <= sym%nop) THEN
    1624           0 :          vecout = vecout1
    1625             :       ELSE
    1626           0 :          vecout = conjg(vecout1)
    1627             :       END IF
    1628             : 
    1629           0 :    END SUBROUTINE ket_trafo1
    1630             : 
    1631             :    ! Determines common phase factor (with unit norm)
    1632           0 :    SUBROUTINE commonphase(cfac, carr, n)
    1633             :       IMPLICIT NONE
    1634             :       INTEGER, INTENT(IN)      :: n
    1635             :       COMPLEX, INTENT(IN)      :: carr(n)
    1636             :       COMPLEX, INTENT(OUT)     :: cfac
    1637             :       REAL                    :: rdum, rmax
    1638             :       INTEGER                 :: i
    1639             : 
    1640             : !       IF( all( abs(carr) .lt. 1E-12 ) ) THEN
    1641             : !         cfac = 1
    1642             : !         RETURN
    1643             : !       END IF
    1644             : 
    1645           0 :       cfac = 0
    1646           0 :       rmax = 0
    1647           0 :       DO i = 1, n
    1648           0 :          rdum = abs(carr(i))
    1649           0 :          IF (rdum > 1e-6) THEN; cfac = carr(i)/rdum; EXIT
    1650           0 :          ELSE IF (rdum > rmax) THEN; cfac = carr(i)/rdum; rmax = rdum
    1651             :          END IF
    1652             :       END DO
    1653           0 :       IF (cfac == 0 .and. all(carr /= 0)) THEN
    1654           0 :          WRITE (999, *) carr
    1655           0 :          STOP 'commonphase: Could not determine common phase factor. (Wrote carr to fort.999)'
    1656             :       END IF
    1657           0 :    END SUBROUTINE commonphase
    1658             : 
    1659           0 :    SUBROUTINE bramat_trafo( &
    1660           0 :       vecout, igptm_out, vecin, igptm_in, ikpt0, iop, writevec, pointer, sym, &
    1661           0 :       rrot, invrrot, hybrid, kpts, maxlcutm, atoms, lcutm, nindxm, maxindxm, dwgn, nbasp, nbasm)
    1662             : 
    1663             :       USE m_constants
    1664             :       USE m_util
    1665             :       USE m_types
    1666             :       IMPLICIT NONE
    1667             :       TYPE(t_hybrid), INTENT(IN)   :: hybrid
    1668             :       TYPE(t_sym), INTENT(IN)   :: sym
    1669             :       TYPE(t_kpts), INTENT(IN)   :: kpts
    1670             :       TYPE(t_atoms), INTENT(IN)   :: atoms
    1671             : 
    1672             : !     - scalars
    1673             :       INTEGER, INTENT(IN)      ::  ikpt0, igptm_in, iop, maxindxm
    1674             :       INTEGER, INTENT(IN)      ::  maxlcutm
    1675             :       INTEGER, INTENT(IN)      ::  nbasp
    1676             :       LOGICAL, INTENT(IN)      ::  writevec
    1677             :       INTEGER, INTENT(OUT)     ::  igptm_out
    1678             : !     - arrays -
    1679             :       INTEGER, INTENT(IN)      ::  rrot(3, 3), invrrot(3, 3)
    1680             :       INTEGER, INTENT(IN)      :: lcutm(atoms%ntype),&
    1681             :      &                            nindxm(0:maxlcutm, atoms%ntype)
    1682             :       INTEGER, INTENT(IN)      :: nbasm(kpts%nkptf)
    1683             :       INTEGER, INTENT(IN)      ::  pointer(&
    1684             :      &                          minval(hybrid%gptm(1, :)) - 1:maxval(hybrid%gptm(1, :)) + 1,&
    1685             :      &                          minval(hybrid%gptm(2, :)) - 1:maxval(hybrid%gptm(2, :)) + 1,&
    1686             :      &                          minval(hybrid%gptm(3, :)) - 1:maxval(hybrid%gptm(3, :)) + 1)
    1687             : 
    1688             :       COMPLEX, INTENT(IN)      ::  vecin(nbasm(ikpt0))
    1689             :       COMPLEX, INTENT(IN)      ::  dwgn(-maxlcutm:maxlcutm,&
    1690             :      &                                 -maxlcutm:maxlcutm,&
    1691             :      &                                         0:maxlcutm)
    1692             :       COMPLEX, INTENT(OUT)     ::  vecout(nbasm(ikpt0))
    1693             : 
    1694             : !     - private scalars -
    1695             :       INTEGER                 ::  itype, ieq, ic, l, n, i, nn, i1, i2, j1, j2, m1
    1696             :       INTEGER                 ::  m2, igptm, igptm2, igptp, iiop, isym
    1697             :       INTEGER                 ::  ikpt1, isymi, rcent
    1698             :       LOGICAL                 ::  trs
    1699             :       COMPLEX, PARAMETER       ::  img = (0.0, 1.0)
    1700             :       COMPLEX                 ::  cexp, cdum
    1701             : !     - private arrays -
    1702           0 :       INTEGER                 ::  pnt(maxindxm, 0:maxlcutm, atoms%nat), g(3),&
    1703           0 :      &                            g1(3), iarr(hybrid%ngptm(ikpt0))
    1704             :       REAL                    ::  rkpt(3), rkpthlp(3), trans(3)
    1705           0 :       COMPLEX                 ::  vecin1(nbasm(ikpt0))
    1706           0 :       COMPLEX                 ::  carr(hybrid%ngptm(ikpt0))
    1707             : 
    1708           0 :       IF (iop <= sym%nop) THEN
    1709           0 :          isym = iop
    1710           0 :          trs = .false.
    1711           0 :          trans = sym%tau(:, isym)
    1712             :       ELSE
    1713           0 :          isym = iop - sym%nop
    1714           0 :          trs = .true.
    1715           0 :          trans = sym%tau(:, isym)
    1716             :       END IF
    1717             : 
    1718           0 :       rkpthlp = matmul(rrot, kpts%bk(:, ikpt0))
    1719           0 :       rkpt = modulo1(rkpthlp, kpts%nkpt3)
    1720           0 :       g = nint(rkpthlp - rkpt)
    1721             :       !
    1722             :       ! determine number of rotated k-point bk(:,ikpt) -> ikpt1
    1723             :       !
    1724           0 :       DO i = 1, kpts%nkpt
    1725           0 :          IF (maxval(abs(rkpt - kpts%bk(:, i))) <= 1E-06) THEN
    1726             :             ikpt1 = i
    1727             :             EXIT
    1728             :          END IF
    1729             :       END DO
    1730             : 
    1731           0 :       DO igptm = 1, hybrid%ngptm(ikpt1)
    1732           0 :          igptp = hybrid%pgptm(igptm, ikpt1)
    1733           0 :          g1 = matmul(invrrot, hybrid%gptm(:, igptp) - g)
    1734           0 :          igptm2 = pointer(g1(1), g1(2), g1(3))
    1735           0 :          IF (igptm2 == igptm_in) THEN
    1736           0 :             igptm_out = igptm
    1737           0 :             IF (writevec) THEN
    1738           0 :                cdum = exp(img*tpi_const*dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp), trans))
    1739           0 :                EXIT
    1740             :             ELSE
    1741             :                RETURN
    1742             :             END IF
    1743             :          END IF
    1744             :       END DO
    1745             : 
    1746             : !     Transform back to unsymmetrized product basis in case of inversion symmetry.
    1747           0 :       vecout = vecin
    1748           0 :       if (sym%invs) CALL desymmetrize(vecout, nbasp, 1, 1, &
    1749           0 :                                       atoms, lcutm, maxlcutm, nindxm, sym)
    1750             : 
    1751             : !     Right-multiplication
    1752             :       ! PW
    1753           0 :       IF (trs) THEN; vecin1 = cdum*conjg(vecout)
    1754           0 :       ELSE; vecin1 = cdum*vecout
    1755             :       END IF
    1756             : 
    1757             : !     Define pointer to first mixed-basis functions (with m = -l)
    1758           0 :       i = 0
    1759           0 :       ic = 0
    1760           0 :       DO itype = 1, atoms%ntype
    1761           0 :          DO ieq = 1, atoms%neq(itype)
    1762           0 :             ic = ic + 1
    1763           0 :             DO l = 0, lcutm(itype)
    1764           0 :                DO n = 1, nindxm(l, itype)
    1765           0 :                   i = i + 1
    1766           0 :                   pnt(n, l, ic) = i
    1767             :                END DO
    1768           0 :                i = i + nindxm(l, itype)*2*l
    1769             :             END DO
    1770             :          END DO
    1771             :       END DO
    1772             : 
    1773             : !     Left-multiplication
    1774             :       ! MT
    1775           0 :       cexp = exp(-img*tpi_const*dot_product(kpts%bk(:, ikpt1) + g, trans))
    1776           0 :       ic = 0
    1777           0 :       DO itype = 1, atoms%ntype
    1778           0 :          DO ieq = 1, atoms%neq(itype)
    1779           0 :             ic = ic + 1
    1780           0 :             rcent = hybrid%map(ic, sym%invtab(isym))
    1781           0 :             cdum = cexp*exp(img*tpi_const*dot_product(g, atoms%taual(:, ic))) !rcent)))
    1782           0 :             cdum = conjg(cdum)
    1783           0 :             DO l = 0, lcutm(itype)
    1784           0 :                nn = nindxm(l, itype)
    1785           0 :                DO n = 1, nn
    1786             : 
    1787           0 :                   i1 = pnt(n, l, ic)
    1788           0 :                   i2 = i1 + nn*2*l
    1789           0 :                   j1 = pnt(n, l, rcent)
    1790           0 :                   j2 = j1 + nn*2*l
    1791             : 
    1792           0 :                   vecout(i1:i2:nn) = cdum*matmul(dwgn(-l:l, -l:l, l), vecin1(j1:j2:nn))
    1793             : 
    1794             :                END DO
    1795             :             END DO
    1796             :          END DO
    1797             :       END DO
    1798             : 
    1799             :       ! PW
    1800           0 :       DO igptm = 1, hybrid%ngptm(ikpt1)
    1801           0 :          igptp = hybrid%pgptm(igptm, ikpt1)
    1802           0 :          g1 = matmul(invrrot, hybrid%gptm(:, igptp) - g)
    1803           0 :          iarr(igptm) = pointer(g1(1), g1(2), g1(3))
    1804           0 :          carr(igptm) = exp(-img*tpi_const*dot_product(kpts%bk(:, ikpt1) + hybrid%gptm(:, igptp), trans))
    1805             :       END DO
    1806           0 :       DO i1 = 1, hybrid%ngptm(ikpt1)
    1807           0 :          vecout(nbasp + i1) = carr(i1)*vecin1(nbasp + iarr(i1))
    1808             :       END DO
    1809             : 
    1810             :       ! If inversion symmetry is applicable, symmetrize to make the values real.
    1811           0 :       if (sym%invs) CALL symmetrize(vecout, nbasp, 1, 1, .false., &
    1812           0 :                                     atoms, lcutm, maxlcutm, nindxm, sym)
    1813             : 
    1814             :    END SUBROUTINE bramat_trafo
    1815             : 
    1816             : END MODULE m_trafo

Generated by: LCOV version 1.13