LCOV - code coverage report
Current view: top level - hybrid - symmetrizeh.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 330 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 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_symmetrizeh
       8             : 
       9             : ! symmetrize the Hamiltonian according to the symmetry operations of the little group of k
      10             : 
      11             : CONTAINS
      12             : 
      13           0 :    SUBROUTINE symmetrizeh(atoms, bk, DIMENSION, jsp, lapw, sym, kveclo, cell, nsymop, psym, hmat)
      14             : 
      15             :       USE m_constants
      16             :       USE m_types
      17             : 
      18             :       IMPLICIT NONE
      19             : 
      20             :       TYPE(t_dimension), INTENT(IN)    :: DIMENSION
      21             :       TYPE(t_sym), INTENT(IN)    :: sym
      22             :       TYPE(t_cell), INTENT(IN)    :: cell
      23             :       TYPE(t_atoms), INTENT(IN)    :: atoms
      24             :       TYPE(t_lapw), INTENT(IN)    :: lapw
      25             :       TYPE(t_mat), INTENT(INOUT) :: hmat
      26             : 
      27             :       ! scalars
      28             :       INTEGER, INTENT(IN)    :: nsymop, jsp
      29             : 
      30             :       ! arrays
      31             :       INTEGER, INTENT(IN)    :: kveclo(atoms%nlotot)
      32             :       INTEGER, INTENT(IN)    :: psym(nsymop)
      33             :       REAL, INTENT(IN)    :: bk(3)
      34             : 
      35             :       ! local scalars
      36             :       INTEGER               ::  ilotot, itype, itype1, ilo, ilo1
      37             :       INTEGER               ::  iatom, iatom1, iiatom, iiatom1
      38             :       INTEGER               ::  i, ieq, ieq1, m
      39             :       INTEGER               ::  igpt_lo, igpt_lo1, igpt_lo2, igpt1_lo1
      40             :       INTEGER               ::  igpt1_lo2, isym, iop, ic, ic1, ic2
      41             :       INTEGER               ::  igpt, igpt1, igpt2, igpt3
      42             :       INTEGER               ::  invsfct, invsfct1, idum
      43             :       INTEGER               ::  l, l1, lm, j, ok, ratom, ratom1, nrgpt
      44             :       COMPLEX, PARAMETER     ::  img = (0.0, 1.0)
      45             :       COMPLEX               ::  cdum, cdum2
      46             : 
      47             :       ! local arrays
      48           0 :       INTEGER               ::  l_lo(atoms%nlotot)
      49           0 :       INTEGER               ::  itype_lo(atoms%nlotot)
      50           0 :       INTEGER               ::  gpt_lo(3, atoms%nlotot), gpthlp(3), g(3)
      51             :       INTEGER               ::  indx(DIMENSION%nbasfcn, DIMENSION%nbasfcn)
      52           0 :       INTEGER               ::  lo_indx(atoms%nlod, atoms%nat)
      53           0 :       INTEGER               ::  rot(3, 3, nsymop), rrot(3, 3, nsymop)
      54             : 
      55           0 :       INTEGER, ALLOCATABLE   ::  pointer_apw(:, :)
      56           0 :       INTEGER, ALLOCATABLE   ::  ipiv(:)
      57           0 :       INTEGER, ALLOCATABLE   ::  map(:, :)
      58             : 
      59             :       REAL                  ::  rtaual(3), kghlp(3)
      60             :       REAL                  ::  rotkpthlp(3), rotkpt(3)
      61           0 :       REAL                  ::  trans(3, nsymop)
      62           0 :       COMPLEX, ALLOCATABLE   ::  c_lo(:, :, :, :), c_rot(:, :, :, :, :), y(:)
      63           0 :       COMPLEX, ALLOCATABLE   ::  cfac(:, :), chelp(:, :)
      64             : 
      65           0 :       LOGICAL               ::  ldum(lapw%nv(jsp) + atoms%nlotot, lapw%nv(jsp) + atoms%nlotot)
      66             : 
      67           0 :       TYPE(t_mat) :: hmatTemp
      68             : 
      69           0 :       CALL hmatTemp%init(hmat%l_real, hmat%matsize1, hmat%matsize2)
      70           0 :       IF (hmat%l_real) THEN
      71           0 :          hmatTemp%data_r = hmat%data_r
      72             :       ELSE
      73           0 :          hmatTemp%data_c = CONJG(hmat%data_c)
      74             :       END IF
      75             : 
      76             :       ! calculate rotations in reciprocal space
      77           0 :       DO isym = 1, nsymop
      78           0 :          iop = psym(isym)
      79           0 :          IF (iop <= sym%nop) THEN
      80           0 :             rrot(:, :, isym) = TRANSPOSE(sym%mrot(:, :, sym%invtab(iop)))
      81             :          ELSE
      82           0 :             rrot(:, :, isym) = -TRANSPOSE(sym%mrot(:, :, sym%invtab(iop - sym%nop)))
      83             :          END IF
      84             :       END DO
      85             : 
      86             :       ! calculate rotations in real space (internal coordinates)
      87           0 :       DO isym = 1, nsymop
      88           0 :          iop = psym(isym)
      89           0 :          IF (iop <= sym%nop) THEN
      90           0 :             rot(:, :, isym) = sym%mrot(:, :, iop)
      91           0 :             trans(:, isym) = sym%tau(:, iop)
      92             :          ELSE
      93           0 :             rot(:, :, isym) = sym%mrot(:, :, iop - sym%nop)
      94           0 :             trans(:, isym) = sym%tau(:, iop - sym%nop)
      95             :          END IF
      96             :       END DO
      97             : 
      98             :       ! caclulate mapping of atoms
      99           0 :       ALLOCATE (map(nsymop, atoms%nat))
     100           0 :       map = 0
     101           0 :       iatom = 0
     102           0 :       iiatom = 0
     103           0 :       DO itype = 1, atoms%ntype
     104           0 :          DO ieq = 1, atoms%neq(itype)
     105           0 :             iatom = iatom + 1
     106           0 :             DO isym = 1, nsymop
     107           0 :                rtaual = MATMUL(rot(:, :, isym), atoms%taual(:, iatom)) + trans(:, isym)
     108           0 :                iatom1 = 0
     109           0 :                DO ieq1 = 1, atoms%neq(itype)
     110           0 :                   IF (ALL(ABS(MODULO(rtaual - atoms%taual(:, iiatom + ieq1) + 1e-12, 1.0)) < 1e-10)) THEN  !The 1e-12 is a dirty fix.
     111           0 :                      iatom1 = iiatom + ieq1
     112             :                   END IF
     113             :                END DO
     114           0 :                IF (iatom1 == 0) STOP 'symmetrizeh_new: error finding rotated atomic position'
     115           0 :                map(isym, iatom) = iatom1
     116             :             END DO
     117             :          END DO
     118           0 :          iiatom = iiatom + atoms%neq(itype)
     119             :       END DO
     120             : 
     121             :       ! initialze pointer_apw and the apw part of cfac
     122           0 :       ALLOCATE (pointer_apw(lapw%nv(jsp), nsymop), cfac(lapw%nv(jsp) + atoms%nlotot, nsymop), stat=ok)
     123           0 :       IF (ok /= 0) STOP 'symmetrizeh_new: failure allocation pointer_apw,cfac'
     124             : 
     125           0 :       pointer_apw = 0
     126           0 :       cfac = 0
     127             : 
     128           0 :       DO isym = 1, nsymop
     129             : 
     130             :          ! determine vector g, which map Rk back into BZ
     131             :          ! Rk - G = k => G = Rk-k
     132           0 :          rotkpt = MATMUL(rrot(:, :, isym), bk(:))
     133           0 :          g = NINT(rotkpt - bk)
     134             : 
     135           0 :          DO igpt = 1, lapw%nv(jsp)
     136             :             !rotate G vector corresponding to isym
     137           0 :             gpthlp = MATMUL(rrot(:, :, isym), lapw%gvec(:, igpt, jsp)) + g
     138             :             ! determine number of gpthlp
     139           0 :             nrgpt = 0
     140           0 :             DO i = 1, lapw%nv(jsp)
     141           0 :                IF (MAXVAL(ABS(gpthlp - lapw%gvec(:, i, jsp))) <= 1E-06) THEN
     142             :                   nrgpt = i
     143             :                   EXIT
     144             :                END IF
     145             :             END DO
     146           0 :             IF (nrgpt == 0) THEN
     147           0 :                PRINT *, igpt
     148           0 :                PRINT *, lapw%gvec(:, igpt, jsp)
     149           0 :                PRINT *, gpthlp
     150           0 :                PRINT *, g
     151           0 :                PRINT *, bk
     152           0 :                DO i = 1, lapw%nv(jsp)
     153           0 :                   WRITE (6, *) i, lapw%gvec(:, i, jsp)
     154             :                ENDDO
     155           0 :                STOP 'symmetrizeh_new: rotated G point not found'
     156             :             END IF
     157           0 :             pointer_apw(igpt, isym) = nrgpt
     158           0 :             cfac(igpt, isym) = EXP(-2*pi_const*img*(dot_PRODUCT(bk(:) + gpthlp(:), trans(:, isym))))
     159             :          END DO
     160             :       END DO
     161             : 
     162             :       ! average apw-part of symmetry-equivalent matrix elements
     163             : 
     164           0 :       ldum = .TRUE.
     165           0 :       DO i = 1, lapw%nv(jsp)
     166           0 :          DO j = 1, i
     167           0 :             cdum = 0
     168           0 :             ic = 0
     169           0 :             DO isym = 1, nsymop
     170           0 :                iop = psym(isym)
     171           0 :                igpt = pointer_apw(i, isym)
     172           0 :                igpt1 = pointer_apw(j, isym)
     173             : 
     174           0 :                IF (iop <= sym%nop) THEN
     175           0 :                   IF ((igpt /= 0) .AND. (igpt1 /= 0)) THEN
     176           0 :                      ic = ic + 1
     177           0 :                      IF (hmatTemp%l_real) THEN
     178           0 :                         cdum = cdum + CONJG(cfac(i, isym))*hmatTemp%data_r(igpt1, igpt)*cfac(j, isym)
     179             :                      ELSE
     180           0 :                         cdum = cdum + CONJG(cfac(i, isym))*hmatTemp%data_c(igpt1, igpt)*cfac(j, isym)
     181             :                      END IF
     182             :                   END IF
     183             :                ELSE
     184           0 :                   IF ((igpt /= 0) .AND. (igpt1 /= 0)) THEN
     185           0 :                      ic = ic + 1
     186           0 :                      IF (hmatTemp%l_real) THEN
     187           0 :                         cdum = cdum + CONJG(CONJG(cfac(i, isym))*hmatTemp%data_r(igpt1, igpt)*cfac(j, isym))
     188             :                      ELSE
     189           0 :                         cdum = cdum + CONJG(CONJG(cfac(i, isym))*hmatTemp%data_c(igpt1, igpt)*cfac(j, isym))
     190             :                      END IF
     191             :                   END IF
     192             :                END IF
     193             :             END DO
     194             : 
     195             :             cdum = cdum!/ic
     196           0 :             DO isym = 1, nsymop
     197           0 :                iop = psym(isym)
     198           0 :                igpt = pointer_apw(i, isym)
     199           0 :                igpt1 = pointer_apw(j, isym)
     200           0 :                IF ((igpt == 0) .OR. (igpt1 == 0)) CYCLE
     201           0 :                IF (igpt1 > igpt) CYCLE
     202           0 :                IF (ldum(igpt, igpt1)) THEN
     203           0 :                   IF (hmat%l_real) THEN
     204           0 :                      IF (iop <= sym%nop) THEN
     205           0 :                         hmat%data_r(igpt1, igpt) = cdum/(CONJG(cfac(i, isym))*cfac(j, isym))
     206           0 :                         ldum(igpt, igpt1) = .FALSE.
     207             :                      ELSE
     208           0 :                         hmat%data_r(igpt1, igpt) = CONJG(cdum/(CONJG(cfac(i, isym))*cfac(j, isym)))
     209           0 :                         ldum(igpt, igpt1) = .FALSE.
     210             :                      END IF
     211           0 :                      hmat%data_r(igpt, igpt1) = hmat%data_r(igpt1, igpt)
     212             :                   ELSE
     213           0 :                      IF (iop <= sym%nop) THEN
     214           0 :                         hmat%data_c(igpt1, igpt) = cdum/(CONJG(cfac(i, isym))*cfac(j, isym))
     215           0 :                         ldum(igpt, igpt1) = .FALSE.
     216             :                      ELSE
     217           0 :                         hmat%data_c(igpt1, igpt) = CONJG(cdum/(CONJG(cfac(i, isym))*cfac(j, isym)))
     218           0 :                         ldum(igpt, igpt1) = .FALSE.
     219             :                      END IF
     220           0 :                      hmat%data_c(igpt, igpt1) = CONJG(hmat%data_c(igpt1, igpt))
     221             :                   END IF
     222             :                END IF
     223             :             END DO
     224             :          END DO
     225             :       END DO
     226             : 
     227             :       ! average lo-part of matrix elements
     228           0 :       IF (ANY(atoms%nlo /= 0)) THEN
     229             : 
     230             :          ! preparations
     231           0 :          ilotot = 0
     232           0 :          iatom = 0
     233           0 :          DO itype = 1, atoms%ntype
     234           0 :             DO ieq = 1, atoms%neq(itype)
     235           0 :                iatom = iatom + 1
     236           0 :                IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN
     237           0 :                   IF (atoms%invsat(iatom) == 0) invsfct = 1
     238           0 :                   IF (atoms%invsat(iatom) == 1) invsfct = 2
     239           0 :                   DO ilo = 1, atoms%nlo(itype)
     240           0 :                      l = atoms%llo(ilo, itype)
     241           0 :                      DO m = 1, invsfct*(2*l + 1)
     242           0 :                         ilotot = ilotot + 1
     243           0 :                         l_lo(ilotot) = l
     244           0 :                         itype_lo(ilotot) = itype
     245           0 :                         gpt_lo(:, ilotot) = lapw%gvec(:, kveclo(ilotot), jsp)
     246             :                      END DO
     247             :                   END DO
     248             :                END IF
     249             :             END DO
     250             :          END DO
     251             : 
     252             :          ! calculate expansion coefficients for local orbitals
     253           0 :          IF (hmat%l_real) THEN
     254           0 :             ALLOCATE (c_lo(4*MAXVAL(l_lo) + 2, 4*MAXVAL(l_lo) + 2, atoms%nlod, atoms%nat), stat=ok)
     255             :          ELSE
     256           0 :             ALLOCATE (c_lo(2*MAXVAL(l_lo) + 1, 2*MAXVAL(l_lo) + 1, atoms%nlod, atoms%nat), stat=ok)
     257             :          END IF
     258             : 
     259           0 :          IF (ok /= 0) STOP 'symmetrizeh_new: failure allocation c_lo'
     260             : 
     261             :          iatom = 0
     262             :          ilotot = 0
     263           0 :          lo_indx = 0
     264           0 :          DO itype = 1, atoms%ntype
     265           0 :             DO ieq = 1, atoms%neq(itype)
     266           0 :                iatom = iatom + 1
     267           0 :                IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN
     268           0 :                   IF (atoms%invsat(iatom) == 0) invsfct = 1
     269           0 :                   IF (atoms%invsat(iatom) == 1) invsfct = 2
     270             : 
     271           0 :                   DO ilo = 1, atoms%nlo(itype)
     272           0 :                      l = atoms%llo(ilo, itype)
     273           0 :                      ALLOCATE (y((l + 1)**2))
     274           0 :                      lo_indx(ilo, iatom) = ilotot + 1
     275             : 
     276           0 :                      DO igpt_lo = 1, invsfct*(2*l + 1)
     277           0 :                         ilotot = ilotot + 1
     278           0 :                         kghlp = bk(:) + gpt_lo(:, ilotot)
     279             : 
     280             :                         !generate spherical harmonics
     281           0 :                         CALL harmonicsr(y, MATMUL(kghlp, cell%bmat), l)
     282             : 
     283           0 :                         lm = l**2
     284           0 :                         cdum = EXP(2*pi_const*img*dot_PRODUCT(kghlp, atoms%taual(:, iatom)))
     285             : 
     286           0 :                         IF (invsfct == 1) THEN
     287           0 :                            DO m = 1, 2*l + 1
     288           0 :                               lm = lm + 1
     289           0 :                               c_lo(m, igpt_lo, ilo, iatom) = cdum*CONJG(y(lm))
     290             :                            END DO
     291             :                         ELSE
     292           0 :                            DO m = 1, 2*l + 1
     293           0 :                               lm = lm + 1
     294           0 :                               c_lo(m, igpt_lo, ilo, iatom) = cdum*CONJG(y(lm))
     295           0 :                               c_lo(4*l + 3 - m, igpt_lo, ilo, iatom) = (-1)**(m - 1)*CONJG(cdum)*y(lm)
     296             :                            END DO
     297             :                         END IF
     298             :                      END DO
     299           0 :                      DEALLOCATE (y)
     300             :                   END DO
     301             :                END IF
     302             :             END DO
     303             :          END DO
     304             : 
     305           0 :          IF (ilotot /= atoms%nlotot) STOP 'symmetrizeh_new: failure counting local orbitals(ilotot)'
     306             : 
     307           0 :          IF (hmat%l_real) THEN
     308           0 :             ALLOCATE (c_rot(4*MAXVAL(l_lo) + 2, 4*MAXVAL(l_lo) + 2, atoms%nlod, atoms%nat, nsymop))
     309           0 :             ALLOCATE (chelp(4*MAXVAL(l_lo) + 2, 4*MAXVAL(l_lo) + 2))
     310             :          ELSE
     311           0 :             ALLOCATE (c_rot(2*MAXVAL(l_lo) + 1, 2*MAXVAL(l_lo) + 1, atoms%nlod, atoms%nat, nsymop))
     312           0 :             ALLOCATE (chelp(2*MAXVAL(l_lo) + 1, 2*MAXVAL(l_lo) + 1))
     313             :          END IF
     314             : 
     315           0 :          DO isym = 1, nsymop
     316             :             ! determine vector g, which map Rk back into BZ
     317             :             ! Rk - G = k => G = Rk-k
     318           0 :             rotkpt = MATMUL(rrot(:, :, isym), bk(:))
     319           0 :             g = NINT(rotkpt - bk)
     320             : 
     321           0 :             ilotot = 0
     322           0 :             iatom = 0
     323           0 :             DO itype = 1, atoms%ntype
     324           0 :                DO ieq = 1, atoms%neq(itype)
     325           0 :                   iatom = iatom + 1
     326           0 :                   ratom = map(isym, iatom)
     327             : 
     328           0 :                   IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN
     329           0 :                      IF (atoms%invsat(iatom) == 0) invsfct = 1
     330           0 :                      IF (atoms%invsat(iatom) == 1) THEN
     331           0 :                         IF (atoms%invsat(ratom) == 2) THEN
     332           0 :                            ratom = sym%invsatnr(ratom)
     333             :                         END IF
     334             :                         invsfct = 2
     335             :                      END IF
     336             : 
     337           0 :                      DO ilo = 1, atoms%nlo(itype)
     338           0 :                         l = atoms%llo(ilo, itype)
     339           0 :                         ALLOCATE (y((l + 1)**2))
     340             : 
     341           0 :                         DO igpt_lo = 1, invsfct*(2*l + 1)
     342           0 :                            ilotot = ilotot + 1
     343             : 
     344             :                            !rotate G_lo corresponding to iop
     345           0 :                            gpthlp = MATMUL(rrot(:, :, isym), gpt_lo(:, ilotot)) + g
     346           0 :                            kghlp = bk(:) + MATMUL(rrot(:, :, isym), gpt_lo(:, ilotot)) + g
     347             : 
     348             :                            !generate spherical harmonics
     349           0 :                            CALL harmonicsr(y, MATMUL(kghlp, cell%bmat), l)
     350             : 
     351           0 :                            lm = l**2
     352           0 :                            cdum = EXP(2*pi_const*img*dot_PRODUCT(kghlp, atoms%taual(:, ratom)))
     353           0 :                            IF (invsfct == 1) THEN
     354           0 :                               DO m = 1, 2*l + 1
     355           0 :                                  lm = lm + 1
     356           0 :                                  c_rot(m, igpt_lo, ilo, ratom, isym) = cdum*CONJG(y(lm))
     357             :                               END DO
     358             :                            ELSE
     359           0 :                               DO m = 1, 2*l + 1
     360           0 :                                  lm = lm + 1
     361           0 :                                  c_rot(m, igpt_lo, ilo, ratom, isym) = cdum*CONJG(y(lm))
     362           0 :                                  c_rot(4*l + 3 - m, igpt_lo, ilo, ratom, isym) = (-1)**(m - 1)*(CONJG(cdum))*y(lm)
     363             :                               END DO
     364             :                            END IF
     365           0 :                            cfac(lapw%nv(jsp) + ilotot, isym) = EXP(-2*pi_const*img*(dot_PRODUCT(kghlp, trans(:, isym))))
     366             :                         END DO
     367             : 
     368           0 :                         idum = invsfct*(2*l + 1)
     369             : 
     370           0 :                         ALLOCATE (ipiv(idum))
     371             : 
     372           0 :                         chelp(:idum, :idum) = c_lo(:idum, :idum, ilo, ratom)
     373             : 
     374           0 :                         CALL ZGESV(idum, idum, chelp(:idum, :idum), idum, ipiv, c_rot(:idum, :idum, ilo, ratom, isym), idum, ok)
     375             : 
     376           0 :                         IF (ok /= 0) STOP 'symmetrizeh_new: failure zgesv'
     377           0 :                         DEALLOCATE (ipiv, y)
     378             :                      END DO
     379             :                   END IF
     380             :                END DO
     381             :             END DO
     382             :          END DO
     383             : 
     384             :          ! symmetrize local-orbital-apw-part of the matrix
     385           0 :          i = 0
     386           0 :          iatom = 0
     387           0 :          DO itype = 1, atoms%ntype
     388           0 :             DO ieq = 1, atoms%neq(itype)
     389           0 :                iatom = iatom + 1
     390           0 :                IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN
     391           0 :                   IF (atoms%invsat(iatom) == 0) invsfct = 1
     392           0 :                   IF (atoms%invsat(iatom) == 1) invsfct = 2
     393             : 
     394           0 :                   DO ilo = 1, atoms%nlo(itype)
     395           0 :                      l = atoms%llo(ilo, itype)
     396           0 :                      DO igpt = 1, invsfct*(2*l + 1)
     397           0 :                         i = i + 1
     398           0 :                         DO j = 1, lapw%nv(jsp)
     399           0 :                            cdum = 0
     400           0 :                            ic = 0
     401           0 :                            DO isym = 1, nsymop
     402           0 :                               ic = ic + 1
     403           0 :                               iop = psym(isym)
     404           0 :                               ratom = map(isym, iatom)
     405           0 :                               IF (invsfct == 2) THEN
     406           0 :                                  IF (atoms%invsat(ratom) == 2) THEN
     407           0 :                                     ratom = sym%invsatnr(ratom)
     408             :                                  END IF
     409             :                               END IF
     410             : 
     411           0 :                               igpt_lo1 = lo_indx(ilo, ratom)
     412           0 :                               igpt_lo2 = igpt_lo1 + invsfct*2*l
     413           0 :                               IF (invsfct == 2) igpt_lo2 = igpt_lo2 + 1
     414           0 :                               igpt1 = pointer_apw(j, isym)
     415             : 
     416           0 :                               cdum2 = 0
     417           0 :                               ic1 = 0
     418           0 :                               DO igpt2 = igpt_lo1, igpt_lo2
     419           0 :                                  ic1 = ic1 + 1
     420           0 :                                  IF (hmatTemp%l_real) THEN
     421           0 :                                     cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))*hmatTemp%data_r(igpt1, lapw%nv(jsp) + igpt2)
     422             :                                  ELSE
     423           0 :                                     cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))*hmatTemp%data_c(igpt1, lapw%nv(jsp) + igpt2)
     424             :                                  END IF
     425             :                               END DO
     426             : 
     427           0 :                               IF (iop <= sym%nop) THEN
     428           0 :                                  cdum = cdum + cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(j, isym)
     429             :                               ELSE
     430           0 :                                  cdum = cdum + CONJG(cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(j, isym))
     431             :                               END IF
     432             :                            END DO
     433           0 :                            IF (hmat%l_real) THEN
     434           0 :                               hmat%data_r(j, lapw%nv(jsp) + i) = cdum!/ic
     435             :                            ELSE
     436           0 :                               hmat%data_c(j, lapw%nv(jsp) + i) = cdum!/ic
     437             :                            END IF
     438             :                         END DO
     439             :                      END DO
     440             :                   END DO
     441             :                END IF
     442             :             END DO
     443             :          END DO
     444             : 
     445             :          !lo's - lo's
     446           0 :          i = 0
     447           0 :          iatom = 0
     448           0 :          DO itype = 1, atoms%ntype
     449           0 :             DO ieq = 1, atoms%neq(itype)
     450           0 :                iatom = iatom + 1
     451           0 :                IF ((atoms%invsat(iatom) == 0) .OR. (atoms%invsat(iatom) == 1)) THEN
     452           0 :                   IF (atoms%invsat(iatom) == 0) invsfct = 1
     453           0 :                   IF (atoms%invsat(iatom) == 1) invsfct = 2
     454             : 
     455           0 :                   DO ilo = 1, atoms%nlo(itype)
     456           0 :                      l = atoms%llo(ilo, itype)
     457           0 :                      DO igpt = 1, invsfct*(2*l + 1)
     458           0 :                         i = i + 1
     459           0 :                         j = 0
     460           0 :                         iatom1 = 0
     461           0 :                         DO itype1 = 1, atoms%ntype
     462           0 :                            DO ieq1 = 1, atoms%neq(itype1)
     463           0 :                               iatom1 = iatom1 + 1
     464           0 :                               IF ((atoms%invsat(iatom1) == 0) .OR. (atoms%invsat(iatom1) == 1)) THEN
     465           0 :                                  IF (atoms%invsat(iatom1) == 0) invsfct1 = 1
     466           0 :                                  IF (atoms%invsat(iatom1) == 1) invsfct1 = 2
     467             : 
     468           0 :                                  DO ilo1 = 1, atoms%nlo(itype1)
     469           0 :                                     l1 = atoms%llo(ilo1, itype1)
     470           0 :                                     DO igpt1 = 1, invsfct1*(2*l1 + 1)
     471           0 :                                        j = j + 1
     472           0 :                                        IF (j > i) CYCLE
     473           0 :                                        cdum = 0
     474           0 :                                        ic = 0
     475           0 :                                        DO isym = 1, nsymop
     476           0 :                                           iop = psym(isym)
     477           0 :                                           ratom = map(isym, iatom)
     478           0 :                                           ratom1 = map(isym, iatom1)
     479             : 
     480           0 :                                           IF (invsfct == 2) THEN
     481           0 :                                              IF (atoms%invsat(ratom) == 2) THEN
     482           0 :                                                 ratom = sym%invsatnr(ratom)
     483             :                                              END IF
     484             :                                           END IF
     485           0 :                                           IF (invsfct1 == 2) THEN
     486           0 :                                              IF (atoms%invsat(ratom1) == 2) THEN
     487           0 :                                                 ratom1 = sym%invsatnr(ratom1)
     488             :                                              END IF
     489             :                                           END IF
     490             : 
     491           0 :                                           igpt_lo1 = lo_indx(ilo, ratom)
     492           0 :                                           igpt_lo2 = igpt_lo1 + invsfct*2*l
     493           0 :                                           IF (invsfct == 2) igpt_lo2 = igpt_lo2 + 1
     494             : 
     495           0 :                                           igpt1_lo1 = lo_indx(ilo1, ratom1)
     496           0 :                                           igpt1_lo2 = igpt1_lo1 + invsfct1*2*l1
     497           0 :                                           IF (invsfct1 == 2) igpt1_lo2 = igpt1_lo2 + 1
     498             : 
     499           0 :                                           cdum2 = 0
     500           0 :                                           ic1 = 0
     501           0 :                                           DO igpt2 = igpt_lo1, igpt_lo2
     502           0 :                                              ic1 = ic1 + 1
     503           0 :                                              ic2 = 0
     504           0 :                                              DO igpt3 = igpt1_lo1, igpt1_lo2
     505           0 :                                                 ic2 = ic2 + 1
     506           0 :                                                 IF (hmatTemp%l_real) THEN
     507             :                                                    cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))* &
     508             :                                                            hmatTemp%data_r(lapw%nv(jsp) + igpt3, lapw%nv(jsp) + igpt2)* &
     509           0 :                                                            c_rot(ic2, igpt1, ilo1, ratom1, isym)
     510             :                                                 ELSE
     511             :                                                    cdum2 = cdum2 + CONJG(c_rot(ic1, igpt, ilo, ratom, isym))* &
     512             :                                                            hmatTemp%data_c(lapw%nv(jsp) + igpt3, lapw%nv(jsp) + igpt2)* &
     513           0 :                                                            c_rot(ic2, igpt1, ilo1, ratom1, isym)
     514             :                                                 END IF
     515             :                                              END DO
     516             :                                           END DO
     517           0 :                                           ic = ic + 1
     518           0 :                                           IF (iop <= sym%nop) THEN
     519           0 :                                              cdum = cdum + cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(lapw%nv(jsp) + j, isym)
     520             :                                           ELSE
     521           0 :                                              cdum = cdum + CONJG(cdum2*CONJG(cfac(lapw%nv(jsp) + i, isym))*cfac(lapw%nv(jsp) + j, isym))
     522             :                                           END IF
     523             :                                        END DO
     524           0 :                                        IF (hmat%l_real) THEN
     525           0 :                                           hmat%data_r(lapw%nv(jsp) + j, lapw%nv(jsp) + i) = cdum!/ic
     526             :                                        ELSE
     527           0 :                                           hmat%data_c(lapw%nv(jsp) + j, lapw%nv(jsp) + i) = cdum!/ic
     528             :                                        END IF
     529             :                                     END DO  ! igpt_lo1
     530             :                                  END DO  ! ilo1
     531             :                               END IF
     532             :                            END DO  !ieq1
     533             :                         END DO  !itype1
     534             :                      END DO  ! igpt_lo
     535             :                   END DO  ! ilo
     536             :                END IF
     537             :             END DO  !ieq
     538             :          END DO  !itype
     539             : 
     540             :       END IF ! ANY(atoms%nlo.NE.0)
     541             : 
     542             :    CONTAINS
     543             : 
     544             :       ! Returns the spherical harmonics Y_lm(^rvec) for l = 0,...,ll in Y(1,...,(ll+1)**2).
     545           0 :       SUBROUTINE harmonicsr(Y, rvec, ll)
     546             :          use m_judft
     547             :          IMPLICIT NONE
     548             :          INTEGER, INTENT(IN)    :: ll
     549             :          REAL, INTENT(IN)       :: rvec(3)
     550             :          COMPLEX, INTENT(OUT)   :: Y((ll + 1)**2)
     551             :          REAL                  :: stheta, ctheta, sphi, cphi, r, rvec1(3)
     552             :          INTEGER               :: l, lm
     553             :          COMPLEX               :: c
     554             :          COMPLEX, PARAMETER     :: img = (0.0, 1.0)
     555             : 
     556           0 :          call timestart("symm harmonics")
     557           0 :          Y(1) = 0.282094791773878
     558           0 :          IF (ll == 0) RETURN
     559             : 
     560             :          stheta = 0
     561             :          ctheta = 0
     562             :          sphi = 0
     563             :          cphi = 0
     564           0 :          r = SQRT(SUM(rvec**2))
     565           0 :          IF (r > 1e-16) THEN
     566           0 :             rvec1 = rvec/r
     567           0 :             ctheta = rvec1(3)
     568           0 :             stheta = SQRT(rvec1(1)**2 + rvec1(2)**2)
     569           0 :             IF (stheta > 1e-16) THEN
     570           0 :                cphi = rvec1(1)/stheta
     571           0 :                sphi = rvec1(2)/stheta
     572             :             END IF
     573             :          ELSE
     574           0 :             Y(2:) = 0.0
     575             :             RETURN
     576             :          END IF
     577             : 
     578             :          ! define Y,l,-l and Y,l,l
     579           0 :          r = Y(1)
     580           0 :          c = 1
     581           0 :          DO l = 1, ll
     582           0 :             r = r*stheta*SQRT(1.0 + 1.0/(2*l))
     583           0 :             c = c*(cphi + img*sphi)
     584           0 :             Y(l**2 + 1) = r*CONJG(c)  ! l,-l
     585           0 :             Y((l + 1)**2) = r*c*(-1)**l ! l,l
     586             :          END DO
     587             : 
     588             :          ! define Y,l,-l+1 and Y,l,l-1
     589           0 :          Y(3) = 0.48860251190292*ctheta
     590           0 :          DO l = 2, ll
     591           0 :             r = SQRT(2.0*l + 1)*ctheta
     592           0 :             Y(l**2 + 2) = r*Y((l - 1)**2 + 1) ! l,-l+1
     593           0 :             Y(l*(l + 2)) = r*Y(l**2)       ! l,l-1
     594             :          END DO
     595             : 
     596             :          ! define Y,l,m, |m|<l-1
     597           0 :          DO l = 2, ll
     598           0 :             lm = l**2 + 2
     599           0 :             DO m = -l + 2, l - 2
     600           0 :                lm = lm + 1
     601           0 :                Y(lm) = SQRT((2.0*l + 1)/(l + m)/(l - m))*(SQRT(2.0*l - 1)*ctheta*Y(lm - 2*l) - SQRT((l + m - 1.0)*(l - m - 1)/(2*l - 3))*Y(lm - 4*l + 2))
     602             :             END DO
     603             :          END DO
     604           0 :          call timestop("symm harmonics")
     605             :       END SUBROUTINE harmonicsr
     606             : 
     607             :    END SUBROUTINE symmetrizeh
     608             : 
     609             : END MODULE m_symmetrizeh

Generated by: LCOV version 1.13