LCOV - code coverage report
Current view: top level - types - types_lapw.F90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 251 320 78.4 %
Date: 2024-04-25 04:21:55 Functions: 11 13 84.6 %

          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_types_lapw
       8             :    USE m_judft
       9             :    use m_types_fleurinput
      10             :    use m_types_nococonv
      11             :    IMPLICIT NONE
      12             :    PRIVATE
      13             :    !These dimensions should be set once per call of FLEUR
      14             :    !They can be queried by the functions lapw%dim_nvd,...
      15             :    !You probably should avoid using the variables directly
      16             :    integer, save :: lapw_dim_nvd
      17             :    integer, save :: lapw_dim_nv2d
      18             :    integer, save :: lapw_dim_nbasfcn
      19             : 
      20             :    TYPE t_lapw
      21             :       INTEGER :: nv(2)
      22             :       INTEGER :: num_local_cols(2)
      23             :       INTEGER :: nv_tot
      24             :       INTEGER :: nmat
      25             :       INTEGER :: nlotot
      26             :       INTEGER, ALLOCATABLE:: k1(:, :)
      27             :       INTEGER, ALLOCATABLE:: k2(:, :)
      28             :       INTEGER, ALLOCATABLE:: k3(:, :)
      29             :       INTEGER, ALLOCATABLE:: gvec(:, :, :) !replaces k1,k2,k3
      30             :       INTEGER, ALLOCATABLE:: kp(:, :)
      31             :       REAL, ALLOCATABLE::rk(:, :)
      32             :       REAL, ALLOCATABLE::gk(:, :, :)
      33             :       REAL, ALLOCATABLE::vk(:, :, :)
      34             :       INTEGER, ALLOCATABLE::index_lo(:, :)
      35             :       INTEGER, ALLOCATABLE::kvec(:, :, :)
      36             :       INTEGER, ALLOCATABLE::nkvec(:, :)
      37             :       REAL   :: bkpt(3)
      38             :       REAL   :: qphon(3)
      39             :    CONTAINS
      40             :       procedure       :: lapw_init => t_lapw_init
      41             :       procedure       :: lapw_init_fi => t_lapw_init_fi
      42             :       GENERIC         :: init => lapw_init, lapw_init_fi
      43             :       PROCEDURE, PASS :: alloc => lapw_alloc
      44             :       PROCEDURE, PASS :: phase_factors => lapw_phase_factors
      45             :       procedure, pass :: hyb_num_bas_fun => hyb_num_bas_fun
      46             :       PROCEDURE, NOPASS:: dim_nvd
      47             :       PROCEDURE, NOPASS:: dim_nv2d
      48             :       PROCEDURE, NOPASS:: dim_nbasfcn
      49             :       PROCEDURE, NOPASS:: init_dim => lapw_init_dim
      50             :    END TYPE t_lapw
      51             :    PUBLIC :: t_lapw, lapw_dim_nbasfcn, lapw_dim_nvd, lapw_dim_nv2d
      52             : 
      53             : CONTAINS
      54         160 :    function hyb_num_bas_fun(lapw, fi) result(nbasfcn)
      55             :       implicit NONE
      56             :       class(t_lapw), intent(in)         :: lapw
      57             :       type(t_fleurinput), intent(in)    :: fi
      58             : 
      59             :       integer :: nbasfcn
      60         160 :       if (fi%noco%l_noco) then
      61           0 :          nbasfcn = lapw%nv(1) + lapw%nv(2) + 2*fi%atoms%nlotot
      62             :       else
      63         160 :          nbasfcn = lapw%nv(1) + fi%atoms%nlotot
      64             :       endif
      65         160 :    end function hyb_num_bas_fun
      66             : 
      67         160 :    subroutine lapw_init_dim(nvd_in, nv2d_in, nbasfcn_in)
      68             :       IMPLICIT NONE
      69             :       INTEGER, INTENT(IN)      :: nvd_in, nv2d_in, nbasfcn_in
      70         160 :       lapw_dim_nvd = nvd_in
      71         160 :       lapw_dim_nv2d = nv2d_in
      72         160 :       lapw_dim_nbasfcn = nbasfcn_in
      73         160 :    end subroutine
      74             : 
      75         795 :    PURE INTEGER function dim_nvd()
      76         795 :       dim_nvd = lapw_dim_nvd
      77         795 :    end function
      78       39526 :    PURE INTEGER function dim_nv2d()
      79       39526 :       dim_nv2d = lapw_dim_nv2d
      80       39526 :    end function
      81         694 :    PURE INTEGER function dim_nbasfcn()
      82         694 :       dim_nbasfcn = lapw_dim_nbasfcn
      83         694 :    end function
      84             : 
      85       16632 :    SUBROUTINE lapw_alloc(lapw, cell, input, noco, nococonv)
      86             :       !
      87             :       !*********************************************************************
      88             :       !     determines dimensions of the lapw basis set with |k+G|<rkmax.
      89             :       !     bkpt is the k-point given in internal units
      90             :       !*********************************************************************
      91             :       USE m_boxdim
      92             :       USE m_types_fleurinput
      93             :       USE m_types_nococonv
      94             :       IMPLICIT NONE
      95             :       TYPE(t_cell), INTENT(IN)      :: cell
      96             :       TYPE(t_input), INTENT(IN)     :: input
      97             :       TYPE(t_noco), INTENT(IN)      :: noco
      98             :       TYPE(t_nococonv), INTENT(IN)  :: nococonv
      99             :       CLASS(t_lapw), INTENT(INOUT)  :: lapw
     100             : 
     101             :       INTEGER j1, j2, j3, mk1, mk2, mk3, nv, addX, addY, addZ
     102             :       INTEGER ispin, nvh(2)
     103             : 
     104             :       REAL arltv1, arltv2, arltv3, rkm, rk2, r2, s(3), sonlyg(3)
     105             :       ! ..
     106             :       !
     107             :       !------->          ABBREVIATIONS
     108             :       !
     109             : 
     110             :       !   rkmax       : cut-off for |g+k|
     111             :       !   arltv(i)    : length of reciprical lattice vector along
     112             :       !                 direction (i)
     113             :       !
     114             :       !---> Determine rkmax box of size mk1, mk2, mk3,
     115             :       !     for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
     116             :       !
     117       16632 :       CALL boxdim(cell%bmat, arltv1, arltv2, arltv3)
     118             : 
     119             :       !     (add 1+1 due to integer rounding, strange k_vector in BZ)
     120       16632 :       mk1 = int(input%rkmax/arltv1) + 2
     121       16632 :       mk2 = int(input%rkmax/arltv2) + 2
     122       16632 :       mk3 = int(input%rkmax/arltv3) + 2
     123             : 
     124       16632 :       rkm = input%rkmax
     125       16632 :       rk2 = rkm*rkm
     126             :       !---> obtain vectors
     127             :       !---> in a spin-spiral calculation different basis sets are used for
     128             :       !---> the two spin directions, because the cutoff radius is defined
     129             :       !---> by |G + k +/- qss/2| < rkmax.
     130       16632 :       nvh(2) = 0
     131       49896 :       DO ispin = 1, MERGE(2, 1, noco%l_ss)
     132       16734 :          addX = abs(NINT((lapw%bkpt(1) + (2*ispin - 3)/2.0*nococonv%qss(1)+lapw%qphon(1))/arltv1))
     133       16734 :          addY = abs(NINT((lapw%bkpt(2) + (2*ispin - 3)/2.0*nococonv%qss(2)+lapw%qPhon(2))/arltv2))
     134       16734 :          addZ = abs(NINT((lapw%bkpt(3) + (2*ispin - 3)/2.0*nococonv%qss(3)+lapw%qphon(3))/arltv3))
     135       16734 :          nv = 0
     136      179032 :          DO j1 = -mk1 - addX, mk1 + addX
     137     1793598 :             DO j2 = -mk2 - addY, mk2 + addY
     138    19588354 :                DO j3 = -mk3 - addZ, mk3 + addZ
     139    71245960 :                   s = lapw%bkpt + (/j1, j2, j3/) + (2*ispin - 3)/2.0*nococonv%qss + lapw%qphon
     140    17811490 :                   sonlyg = (/j1, j2, j3/)
     141   284983840 :                   r2 = dot_PRODUCT(MATMUL(s, cell%bbmat), s)
     142             :                   !r2 = dot_PRODUCT(MATMUL(sonlyg, cell%bbmat), sonlyg)
     143    19426056 :                   IF (r2 .LE. rk2) nv = nv + 1
     144             :                END DO
     145             :             END DO
     146             :          END DO
     147       33366 :          nvh(ispin) = nv
     148             :       END DO
     149       16632 :       nv = MAX(nvh(1), nvh(2))
     150             : 
     151       16632 :       IF (ALLOCATED(lapw%rk)) THEN
     152       42813 :          IF (SIZE(lapw%rk) == nv) THEN
     153         129 :             RETURN !
     154             :          ELSE
     155       14142 :             DEALLOCATE (lapw%rk, lapw%gvec, lapw%vk, lapw%gk)
     156       14142 :             DEALLOCATE (lapw%k1, lapw%k2, lapw%k3)
     157             :          ENDIF
     158             :       ENDIF
     159       66012 :       ALLOCATE (lapw%k1(nv, input%jspins)) !should be removed
     160       49509 :       ALLOCATE (lapw%k2(nv, input%jspins)) !
     161       49509 :       ALLOCATE (lapw%k3(nv, input%jspins)) !
     162       66012 :       ALLOCATE (lapw%rk(nv, input%jspins))
     163       66012 :       ALLOCATE (lapw%gvec(3, nv, input%jspins))
     164       66012 :       ALLOCATE (lapw%vk(3, nv, input%jspins))
     165       49509 :       ALLOCATE (lapw%gk(3, nv, input%jspins))
     166             :      
     167    18051693 :       lapw%rk = 0; lapw%gvec = 0; lapw%nv = 0
     168             :    END SUBROUTINE lapw_alloc
     169             : 
     170         136 :    subroutine t_lapw_init_fi(lapw, fi, nococonv, nk, mpi, dfpt_q) 
     171             :       USE m_types_mpi
     172             :       use m_types_fleurinput
     173             :       implicit none 
     174             :       CLASS(t_lapw), INTENT(INOUT)    :: lapw
     175             :       type(t_fleurinput), intent(in)  :: fi
     176             :       TYPE(t_nococonv), INTENT(IN)    :: nococonv
     177             :       INTEGER, INTENT(IN) :: nk
     178             :       TYPE(t_mpi), INTENT(IN), OPTIONAL:: mpi
     179             :       REAL, INTENT(IN), OPTIONAL :: dfpt_q(3)
     180             : 
     181         136 :       IF (PRESENT(mpi)) THEN
     182           0 :          IF (PRESENT(dfpt_q)) THEN
     183           0 :             call lapw%lapw_init(fi%input, fi%noco, nococonv, fi%kpts, fi%atoms, fi%sym, nk, fi%cell, mpi, dfpt_q)
     184             :          ELSE
     185           0 :             call lapw%lapw_init(fi%input, fi%noco, nococonv, fi%kpts, fi%atoms, fi%sym, nk, fi%cell, mpi)
     186             :          END IF
     187             :       ELSE
     188         136 :          IF (PRESENT(dfpt_q)) THEN
     189           0 :             call lapw%lapw_init(fi%input, fi%noco, nococonv, fi%kpts, fi%atoms, fi%sym, nk, fi%cell, dfpt_q=dfpt_q)
     190             :          ELSE
     191         136 :             call lapw%lapw_init(fi%input, fi%noco, nococonv, fi%kpts, fi%atoms, fi%sym, nk, fi%cell)
     192             :          END IF
     193             :       END IF  
     194         136 :    end subroutine t_lapw_init_fi
     195             : 
     196       16632 :    SUBROUTINE t_lapw_init(lapw, input, noco, nococonv, kpts, atoms, sym, &
     197             :                         nk, cell,  mpi, dfpt_q)
     198             :       USE m_types_mpi
     199             :       USE m_sort
     200             :       USE m_boxdim
     201             :       USE m_types_fleurinput
     202             :       USE m_types_kpts
     203             :       USE m_types_nococonv
     204             :       IMPLICIT NONE
     205             : 
     206             : 
     207             :       TYPE(t_input), INTENT(IN)       :: input
     208             :       TYPE(t_noco), INTENT(IN)        :: noco
     209             :       TYPE(t_nococonv), INTENT(IN)    :: nococonv
     210             :       TYPE(t_cell), INTENT(IN)        :: cell
     211             :       TYPE(t_atoms), INTENT(IN)       :: atoms
     212             :       TYPE(t_sym), INTENT(IN)         :: sym
     213             :       TYPE(t_kpts), INTENT(IN)        :: kpts
     214             :       TYPE(t_mpi), INTENT(IN), OPTIONAL:: mpi
     215             :       CLASS(t_lapw), INTENT(INOUT)    :: lapw
     216             : 
     217             :       REAL, INTENT(IN), OPTIONAL :: dfpt_q(3)
     218             :       !     ..
     219             :       !     .. Scalar Arguments ..
     220             :       INTEGER, INTENT(IN) :: nk
     221             :       !LOGICAL, INTENT(IN)  :: l_zref
     222             :       !     ..
     223             :       !     .. Array Arguments ..
     224             :       !     ..
     225             :       !     .. Local Scalars ..
     226             :       REAL arltv1, arltv2, arltv3, r2, rk2, rkm, r2q, gla, eps, t, r2g, r2phon
     227             :       INTEGER i, j, j1, j2, j3, k, l, mk1, mk2, mk3, n, ispin, gmi, m, nred, n_inner, n_bound, itt(3), addX, addY, addZ
     228             :       !     ..
     229             :       !     .. Local Arrays ..
     230             :       REAL                :: s(3), sq(3), sg(3), qphon(3), sphon(3)
     231       16632 :       REAL, ALLOCATABLE    :: rk(:), rkq(:), rkqq(:), rg(:)
     232       16632 :       INTEGER, ALLOCATABLE :: gvec(:, :), index3(:)
     233             : 
     234       16632 :       call timestart("t_lapw_init")
     235             :       !     ..
     236             :       !---> in a spin-spiral calculation different basis sets are used for
     237             :       !---> the two spin directions, because the cutoff radius is defined
     238             :       !---> by |G + k +/- qss/2| < rkmax.
     239             : 
     240       66528 :       lapw%qphon = [0.0,0.0,0.0]
     241       16632 :       IF (PRESENT(dfpt_q)) lapw%qphon = dfpt_q
     242             :       
     243       16632 :       IF (nk > kpts%nkpt) THEN
     244        1616 :          lapw%bkpt(:) = kpts%bkf(:, nk)
     245             :       ELSE
     246       64912 :          lapw%bkpt(:) = kpts%bk(:, nk)
     247             :       ENDIF
     248             : 
     249       16632 :       CALL lapw%alloc(cell, input, noco, nococonv)
     250             : 
     251       49896 :       ALLOCATE (gvec(3, SIZE(lapw%gvec, 2)))
     252       83160 :       ALLOCATE (rk(SIZE(lapw%gvec, 2)), rkq(SIZE(lapw%gvec, 2)), rkqq(SIZE(lapw%gvec, 2)))
     253       33264 :       ALLOCATE (rg(SIZE(lapw%gvec, 2)))
     254       49896 :       ALLOCATE (index3(SIZE(lapw%gvec, 2)))
     255             : 
     256             :       !---> Determine rkmax box of size mk1, mk2, mk3,
     257             :       !     for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
     258             :       !     arltv(i) length of reciprical lattice vector along direction (i)
     259             :       !
     260       16632 :       CALL boxdim(cell%bmat, arltv1, arltv2, arltv3)
     261             : 
     262             :       !     (add 1+1 due to integer rounding, strange k_vector in BZ)
     263       16632 :       mk1 = int(input%rkmax/arltv1) + 4
     264       16632 :       mk2 = int(input%rkmax/arltv2) + 4
     265       16632 :       mk3 = int(input%rkmax/arltv3) + 4
     266             : 
     267       16632 :       rk2 = input%rkmax*input%rkmax
     268             :       !---> if too many basis functions, reduce rkmax
     269       18834 :       spinloop: DO ispin = 1, input%jspins
     270       16734 :          addX = abs(NINT((lapw%bkpt(1) + (2*ispin - 3)/2.0*nococonv%qss(1)+lapw%qphon(1))/arltv1))
     271       16734 :          addY = abs(NINT((lapw%bkpt(2) + (2*ispin - 3)/2.0*nococonv%qss(2)+lapw%qphon(2))/arltv2))
     272       16734 :          addZ = abs(NINT((lapw%bkpt(3) + (2*ispin - 3)/2.0*nococonv%qss(3)+lapw%qphon(3))/arltv3))
     273             :          !--->    obtain vectors
     274       16734 :          n = 0
     275      245968 :          DO j1 = -mk1 - addX, mk1 + addX
     276     3429862 :             DO j2 = -mk2 - addY, mk2 + addY
     277    50907506 :                DO j3 = -mk3 - addZ, mk3 + addZ
     278   189977512 :                   s = lapw%bkpt + (/j1, j2, j3/) + (2*ispin - 3)/2.0*nococonv%qss + lapw%qphon
     279   189977512 :                   sq = lapw%bkpt + (/j1, j2, j3/)
     280   189977512 :                   sg = (/j1, j2, j3/)
     281   759910048 :                   r2 = dot_PRODUCT(s, MATMUL(s, cell%bbmat))
     282   759910048 :                   r2q = dot_PRODUCT(sq, MATMUL(sq, cell%bbmat))
     283   759910048 :                   r2g = dot_PRODUCT(sg, MATMUL(sg, cell%bbmat))
     284    50678272 :                   IF (r2 .LE. rk2) THEN
     285             :                   !IF (r2g .LE. rk2) THEN
     286     2043850 :                      n = n + 1
     287     8175400 :                      gvec(:, n) = (/j1, j2, j3/)
     288     2043850 :                      rk(n) = SQRT(r2)
     289     2043850 :                      rkq(n) = SQRT(r2q)
     290     2043850 :                      rg(n) = SQRT(r2g)
     291             :                   END IF
     292             :                ENDDO
     293             :             ENDDO
     294             :          ENDDO
     295       16734 :          lapw%nv(ispin) = n
     296             : 
     297             :          !Sort according to k+g, first construct secondary sort key
     298     2060584 :          DO k = 1, lapw%nv(ispin)
     299             :             rkqq(k) = (mk1 + gvec(1, k)) + (mk2 + gvec(2, k))*(2*mk1 + 1) + &
     300     2060584 :                       (mk3 + gvec(3, k))*(2*mk1 + 1)*(2*mk2 + 1)
     301             :          ENDDO
     302       16734 :          CALL sort(index3(:lapw%nv(ispin)), rkq, rkqq)
     303             :          !CALL sort(index3(:lapw%nv(ispin)), rg, rkqq)
     304     2060584 :          DO n = 1, lapw%nv(ispin)
     305     8175400 :             lapw%gvec(:, n, ispin) = gvec(:, index3(n))
     306     2060584 :             lapw%rk(n, ispin) = rk(index3(n))
     307             :          ENDDO
     308             :          !--->    determine pairs of K-vectors, where K_z = K'_-z to use
     309             :          !--->    z-reflection
     310     2060584 :          DO k = 1, lapw%nv(ispin)
     311     8175400 :             lapw%vk(:, k, ispin) = lapw%bkpt + lapw%gvec(:, k, ispin) + (ispin - 1.5)*nococonv%qss + lapw%qphon
     312    14323684 :             lapw%gk(:, k, ispin) = MATMUL(TRANSPOSE(cell%bmat), lapw%vk(:, k, ispin))/MAX(lapw%rk(k, ispin), 1.0e-30)
     313             :          ENDDO
     314             : 
     315       18834 :          IF (.NOT. noco%l_ss .AND. input%jspins == 2) THEN
     316             :             !Second spin is the same
     317       14532 :             lapw%nv(2) = lapw%nv(1)
     318     6278380 :             lapw%gvec(:, :, 2) = lapw%gvec(:, :, 1)
     319     1580494 :             lapw%rk(:, 2) = lapw%rk(:, 1)
     320     6278380 :             lapw%vk(:, :, 2) = lapw%vk(:, :, 1)
     321     6278380 :             lapw%gk(:, :, 2) = lapw%gk(:, :, 1)
     322             :             EXIT spinloop
     323             :          END IF
     324             : 
     325             :       ENDDO spinloop
     326             :       !should be removed later...
     327     7315420 :       lapw%k1 = lapw%gvec(1, :, :)
     328     7315420 :       lapw%k2 = lapw%gvec(2, :, :)
     329     7315420 :       lapw%k3 = lapw%gvec(3, :, :)
     330             : 
     331             :       !Count No of lapw distributed to this PE
     332       49896 :       lapw%num_local_cols = 0
     333       47898 :       DO ispin = 1, input%jspins
     334       47898 :          IF (PRESENT(mpi)) THEN
     335       29698 :             DO k = mpi%n_rank + 1, lapw%nv(ispin), mpi%n_size
     336     2053065 :                lapw%num_local_cols(ispin) = lapw%num_local_cols(ispin) + 1
     337             :             END DO
     338             :          ELSE
     339        1568 :             lapw%num_local_cols(ispin) = lapw%nv(ispin)
     340             :          END IF
     341             :       END DO
     342             : 
     343       19324 :       IF (ANY(atoms%nlo > 0)) CALL priv_lo_basis_setup(lapw, atoms, input, sym, noco, nococonv, cell)
     344             : 
     345       16632 :       lapw%nv_tot = lapw%nv(1)
     346       16632 :       lapw%nmat = lapw%nv(1) + atoms%nlotot
     347       16632 :       IF (noco%l_noco) lapw%nv_tot = lapw%nv_tot + lapw%nv(2)
     348       16632 :       IF (noco%l_noco) lapw%nmat = lapw%nv_tot + 2*atoms%nlotot
     349             : 
     350             : 
     351       16632 :       call timestop("t_lapw_init")
     352             :    CONTAINS
     353             : 
     354       15192 :       SUBROUTINE priv_lo_basis_setup(lapw, atoms, input, sym, noco, nococonv, cell)
     355             :          USE m_types_fleurinput
     356             : 
     357             :          IMPLICIT NONE
     358             :          TYPE(t_lapw), INTENT(INOUT):: lapw
     359             :          TYPE(t_atoms), INTENT(IN)  :: atoms
     360             :          TYPE(t_input), INTENT(IN)  :: input
     361             :          TYPE(t_sym), INTENT(IN)    :: sym
     362             :          TYPE(t_cell), INTENT(IN)   :: cell
     363             :          TYPE(t_noco), INTENT(IN)   :: noco
     364             :          TYPE(t_nococonv), INTENT(IN)   :: nococonv
     365             : 
     366             :          INTEGER:: n, na, nn, np, lo, nkvec_sv, nkvec(atoms%nlod, 2), iindex
     367       15192 :          IF (.NOT. ALLOCATED(lapw%kvec)) THEN
     368        8730 :             ALLOCATE (lapw%kvec(2*(2*atoms%llod + 1), atoms%nlod, atoms%nat))
     369       15164 :             ALLOCATE (lapw%nkvec(atoms%nlod, atoms%nat));lapw%nkvec=0
     370        5238 :             ALLOCATE (lapw%index_lo(atoms%nlod, atoms%nat))
     371             :          ENDIF
     372       15192 :          iindex = 0
     373       15192 :          na = 0
     374       15192 :          nkvec_sv = 0
     375       48254 :          DO n = 1, atoms%ntype
     376       81724 :             DO nn = 1, atoms%neq(n)
     377       33470 :                na = na + 1
     378       33470 :                if (sym%invsat(na) > 1) cycle
     379       33288 :                np = sym%invtab(sym%ngopr(na))
     380       33288 :                CALL priv_vec_for_lo(atoms, input, sym, na, n, np, noco, nococonv, lapw, cell)
     381      103496 :                DO lo = 1, atoms%nlo(n)
     382       37146 :                   lapw%index_lo(lo, na) = iindex
     383       70616 :                   iindex = iindex + lapw%nkvec(lo, na)
     384             :                ENDDO
     385             :             ENDDO
     386             :          ENDDO
     387       15192 :       END SUBROUTINE priv_lo_basis_setup
     388             : 
     389             :    END SUBROUTINE t_lapw_init
     390             : 
     391       60895 :    SUBROUTINE lapw_phase_factors(lapw, iintsp, tau, qss, cph)
     392             :       USE m_constants
     393             :       USE m_types_fleurinput
     394             :       IMPLICIT NONE
     395             :       CLASS(t_lapw), INTENT(in):: lapw
     396             :       INTEGER, INTENT(IN)     :: iintsp
     397             :       REAL, INTENT(in)        :: tau(3), qss(3)
     398             :       COMPLEX, INTENT(out)    :: cph(:)
     399             : 
     400             :       INTEGER:: k
     401             :       REAL:: th
     402             : 
     403             :       !$OMP PARALLEL DO DEFAULT(none) &
     404             :       !$OMP& SHARED(lapw,iintsp,tau,qss,cph)&
     405       60895 :       !$OMP& PRIVATE(k,th)
     406             :       DO k = 1, lapw%nv(iintsp)
     407             :          th = DOT_PRODUCT(lapw%gvec(:, k, iintsp) + (iintsp - 1.5)*qss + lapw%bkpt + lapw%qphon, tau)
     408             :          cph(k) = CMPLX(COS(tpi_const*th), SIN(tpi_const*th))
     409             :       END DO
     410             :       !$OMP END PARALLEL DO
     411       60895 :    END SUBROUTINE lapw_phase_factors
     412             : 
     413             :    SUBROUTINE priv_vec_for_lo_old(atoms, input, sym, na, n, np, noco, nococonv, lapw, cell)
     414             : 
     415             :       USE m_constants
     416             :       USE m_orthoglo
     417             :       USE m_ylm
     418             :       USE m_types_fleurinput
     419             : 
     420             :       IMPLICIT NONE
     421             : 
     422             :       TYPE(t_noco), INTENT(IN)   :: noco
     423             :       TYPE(t_nococonv), INTENT(IN):: nococonv
     424             :       TYPE(t_sym), INTENT(IN)    :: sym
     425             :       TYPE(t_cell), INTENT(IN)   :: cell
     426             :       TYPE(t_atoms), INTENT(IN)  :: atoms
     427             :       TYPE(t_input), INTENT(IN)  :: input
     428             :       TYPE(t_lapw), INTENT(INOUT):: lapw
     429             :       !     ..
     430             :       !     .. Scalar Arguments ..
     431             :       INTEGER, INTENT(IN) :: na, n, np
     432             :       !     ..
     433             :       !     .. Array Arguments ..
     434             :       !     ..
     435             :       !     .. Local Scalars ..
     436             :       COMPLEX term1
     437             :       REAL th, con1
     438             :       INTEGER l, lo, mind, ll1, lm, iintsp, k, nkmin, ntyp, lmp, m, nintsp, k_start
     439             :       LOGICAL linind, enough, l_lo1
     440             :       !     ..
     441             :       !     .. Local Arrays ..
     442             :       INTEGER :: nkvec(atoms%nlod, 2)
     443             :       REAL qssbti(3), bmrot(3, 3), v(3), vmult(3)
     444             :       REAL :: gkrot(3, SIZE(lapw%gk, 2), 2)
     445             :       REAL :: rph(SIZE(lapw%gk, 2), 2)
     446             :       REAL :: cph(SIZE(lapw%gk, 2), 2)
     447             :       COMPLEX ylm((atoms%lmaxd + 1)**2)
     448             :       COMPLEX cwork(-2*atoms%llod:2*atoms%llod + 1, 2*(2*atoms%llod + 1), atoms%nlod, 2)
     449             :       !     ..
     450             :       !     .. Data statements ..
     451             :       REAL, PARAMETER :: eps = 1.0E-8
     452             :       REAL, PARAMETER :: linindq = 1.0e-6
     453             : 
     454             :       con1 = fpi_const/SQRT(cell%omtil)
     455             :       ntyp = n
     456             :       nintsp = MERGE(2, 1, noco%l_ss)
     457             :       DO iintsp = 1, nintsp
     458             :          IF (iintsp .EQ. 1) THEN
     459             :             qssbti = -nococonv%qss/2
     460             :          ELSE
     461             :             qssbti = +nococonv%qss/2
     462             :          ENDIF
     463             : 
     464             :          !--->    set up phase factors
     465             :          DO k = 1, lapw%nv(iintsp)
     466             :             th = tpi_const*DOT_PRODUCT((/lapw%k1(k, iintsp), lapw%k2(k, iintsp), lapw%k3(k, iintsp)/) + qssbti + lapw%qphon, atoms%taual(:, na))
     467             :             rph(k, iintsp) = COS(th)
     468             :             cph(k, iintsp) = -SIN(th)
     469             :          END DO
     470             : 
     471             :          IF (np .EQ. 1) THEN
     472             :             gkrot(:, :, iintsp) = lapw%gk(:, :, iintsp)
     473             :          ELSE
     474             :             bmrot = MATMUL(1.*sym%mrot(:, :, np), cell%bmat)
     475             :             DO k = 1, lapw%nv(iintsp)
     476             :                !-->           apply the rotation that brings this atom into the
     477             :                !-->           representative (this is the definition of ngopr(na))
     478             :                !-->           and transform to cartesian coordinates
     479             :                v(:) = lapw%vk(:, k, iintsp)
     480             :                gkrot(:, k, iintsp) = MATMUL(v, bmrot)
     481             :             END DO
     482             :          END IF
     483             :          !--->   end loop over interstitial spin
     484             :       ENDDO
     485             : 
     486             :       nkvec(:, :) = 0
     487             :       cwork(:, :, :, :) = CMPLX(0.0, 0.0)
     488             :       enough = .FALSE.
     489             : 
     490             :       IF (noco%l_ss) THEN
     491             :          k_start = 2  ! avoid k=1 !!! GB16
     492             :       ELSE
     493             :          k_start = 1
     494             :       ENDIF
     495             : 
     496             :       DO k = k_start, MIN(lapw%nv(1), lapw%nv(nintsp))
     497             : !       IF (ANY(lapw%rk(k,:nintsp).LT.eps)) CYCLE
     498             :          IF (.NOT. enough) THEN
     499             :             DO iintsp = 1, nintsp
     500             : 
     501             :                !-->        generate spherical harmonics
     502             :                vmult(:) = gkrot(:, k, iintsp)
     503             :                CALL ylm4(atoms%lnonsph(ntyp), vmult, ylm)
     504             :                l_lo1 = .false.
     505             :                IF ((lapw%rk(k, iintsp) .LT. eps) .AND. (.not. noco%l_ss)) THEN
     506             :                   l_lo1 = .true.
     507             :                ELSE
     508             :                   l_lo1 = .false.
     509             :                ENDIF
     510             :                ! --> here comes a part of abccoflo()
     511             :                IF (l_lo1) THEN
     512             :                   DO lo = 1, atoms%nlo(ntyp)
     513             :                      IF ((nkvec(lo, iintsp) .EQ. 0) .AND. (atoms%llo(lo, ntyp) .EQ. 0)) THEN
     514             :                         enough = .false.
     515             :                         nkvec(lo, iintsp) = 1
     516             :                         lapw%kvec(nkvec(lo, iintsp), lo, na) = k
     517             :                         term1 = con1*((atoms%rmt(ntyp)**2)/2)
     518             :                         cwork(0, 1, lo, iintsp) = term1/sqrt(2*tpi_const)
     519             :                         IF ((sym%invsat(na) .EQ. 1) .OR. (sym%invsat(na) .EQ. 2)) THEN
     520             :                            cwork(1, 1, lo, iintsp) = conjg(term1)/sqrt(2*tpi_const)
     521             :                         ENDIF
     522             :                      ENDIF
     523             :                   ENDDO
     524             :                ELSE
     525             :                   enough = .TRUE.
     526             :                   term1 = con1*((atoms%rmt(ntyp)**2)/2)*CMPLX(rph(k, iintsp), cph(k, iintsp))
     527             :                   DO lo = 1, atoms%nlo(ntyp)
     528             :                      IF (sym%invsat(na) .EQ. 0) THEN
     529             :                         IF ((nkvec(lo, iintsp)) .LT. (2*atoms%llo(lo, ntyp) + 1)) THEN
     530             :                            enough = .FALSE.
     531             :                            nkvec(lo, iintsp) = nkvec(lo, iintsp) + 1
     532             :                            l = atoms%llo(lo, ntyp)
     533             :                            ll1 = l*(l + 1) + 1
     534             :                            DO m = -l, l
     535             :                               lm = ll1 + m
     536             :                               cwork(m, nkvec(lo, iintsp), lo, iintsp) = term1*ylm(lm)
     537             :                            END DO
     538             :                            CALL orthoglo(input%l_real, atoms, nkvec(lo, iintsp), lo, l, linindq, .FALSE., cwork(-2*atoms%llod, 1, 1, iintsp), linind)
     539             :                            IF (linind) THEN
     540             :                               lapw%kvec(nkvec(lo, iintsp), lo, na) = k
     541             :                            ELSE
     542             :                               nkvec(lo, iintsp) = nkvec(lo, iintsp) - 1
     543             :                            ENDIF
     544             :                         ENDIF
     545             :                      ELSE
     546             :                         IF ((sym%invsat(na) .EQ. 1) .OR. (sym%invsat(na) .EQ. 2)) THEN
     547             :                            IF (nkvec(lo, iintsp) .LT. 2*(2*atoms%llo(lo, ntyp) + 1)) THEN
     548             :                               enough = .FALSE.
     549             :                               nkvec(lo, iintsp) = nkvec(lo, iintsp) + 1
     550             :                               l = atoms%llo(lo, ntyp)
     551             :                               ll1 = l*(l + 1) + 1
     552             :                               DO m = -l, l
     553             :                                  lm = ll1 + m
     554             :                                  mind = -l + m
     555             :                                  cwork(mind, nkvec(lo, iintsp), lo, iintsp) = term1*ylm(lm)
     556             :                                  mind = l + 1 + m
     557             :                                  lmp = ll1 - m
     558             :                                  cwork(mind, nkvec(lo, iintsp), lo, iintsp) = ((-1)**(l + m))*CONJG(term1*ylm(lmp))
     559             :                               END DO
     560             :                               CALL orthoglo(input%l_real, atoms, nkvec(lo, iintsp), lo, l, linindq, .TRUE., cwork(-2*atoms%llod, 1, 1, iintsp), linind)
     561             :                               IF (linind) THEN
     562             :                                  lapw%kvec(nkvec(lo, iintsp), lo, na) = k
     563             :                                  !                          write(*,*) nkvec(lo,iintsp),k,' <- '
     564             :                               ELSE
     565             :                                  nkvec(lo, iintsp) = nkvec(lo, iintsp) - 1
     566             :                               END IF
     567             :                            END IF
     568             :                         END IF
     569             :                      END IF
     570             :                   END DO
     571             :                   IF ((k .EQ. lapw%nv(iintsp)) .AND. (.NOT. enough)) THEN
     572             :                      WRITE (oUnit, FMT=*) 'vec_for_lo did not find enough linearly independent'
     573             :                      WRITE (oUnit, FMT=*) 'clo coefficient-vectors. the linear independence'
     574             :                      WRITE (oUnit, FMT=*) 'quality, linindq, is set: ', linindq
     575             :                      WRITE (oUnit, FMT=*) 'this value might be to large.'
     576             :                      WRITE (*, *) na, k, lapw%nv
     577             :                      CALL juDFT_error("not enough lin. indep. clo-vectors", calledby="vec_for_lo")
     578             :                   END IF
     579             :                   ! -- >        end of abccoflo-part
     580             :                ENDIF
     581             :             ENDDO
     582             :          ENDIF
     583             : 
     584             :          ! -->    check whether we have already enough k-vecs
     585             :          enough = .TRUE.
     586             :          DO lo = 1, atoms%nlo(ntyp)
     587             :             IF (nkvec(lo, 1) .EQ. nkvec(lo, nintsp)) THEN   ! k-vec accepted by both spin channels
     588             :                IF (sym%invsat(na) .EQ. 0) THEN
     589             :                   IF (nkvec(lo, 1) .LT. (2*atoms%llo(lo, ntyp) + 1)) THEN
     590             :                      enough = .FALSE.
     591             :                   ENDIF
     592             :                ELSE
     593             :                   IF (nkvec(lo, 1) .LT. (2*(2*atoms%llo(lo, ntyp) + 1))) THEN
     594             :                      enough = .FALSE.
     595             :                   ENDIF
     596             :                ENDIF
     597             :             ELSE
     598             :                nkmin = MIN(nkvec(lo, 1), nkvec(lo, nintsp)) ! try another k-vec
     599             :                nkvec(lo, 1) = nkmin; nkvec(lo, nintsp) = nkmin
     600             :                enough = .FALSE.
     601             :             ENDIF
     602             :          ENDDO
     603             :          IF (enough) THEN
     604             :             lapw%nkvec(:atoms%nlo(ntyp), na) = nkvec(:atoms%nlo(ntyp), 1)
     605             :             RETURN
     606             :          ENDIF
     607             :       ENDDO
     608             : 
     609             :    END SUBROUTINE priv_vec_for_lo_old
     610             : 
     611       33288 :    SUBROUTINE priv_vec_for_lo(atoms, input, sym, na, ntype, np, noco, nococonv, lapw, cell)
     612             : 
     613             :       USE m_constants
     614             :       USE m_orthoglo
     615             :       USE m_ylm
     616             :       USE m_types_fleurinput
     617             : 
     618             :       IMPLICIT NONE
     619             : 
     620             :       TYPE(t_noco), INTENT(IN)   :: noco
     621             :       TYPE(t_nococonv), INTENT(IN):: nococonv
     622             :       TYPE(t_sym), INTENT(IN)    :: sym
     623             :       TYPE(t_cell), INTENT(IN)   :: cell
     624             :       TYPE(t_atoms), INTENT(IN)  :: atoms
     625             :       TYPE(t_input), INTENT(IN)  :: input
     626             :       TYPE(t_lapw), INTENT(INOUT):: lapw
     627             :       !     ..
     628             :       !     .. Scalar Arguments ..
     629             :       INTEGER, INTENT(IN) :: na, ntype, np
     630             :       !     ..
     631             :       !     .. Array Arguments ..
     632             :       !     ..
     633             :       !     .. Local Scalars ..
     634             :       COMPLEX term1, norm
     635             :       REAL th, con1, linindq, linindqStart, linindqEnd, numK, stepSize
     636             :       INTEGER l, lo, mind, ll1, lm, iintsp, k, nkmin, lmp, m, nintsp, k_start, k_end, minIndex, maxIndex, increment
     637             :       INTEGER nApproach, nApproachEnd
     638             :       LOGICAL linind, enough, l_lo1, l_norm
     639             :       !     ..
     640             :       !     .. Local Arrays ..
     641       33288 :       INTEGER :: nkvec(atoms%nlod, 2)
     642             :       REAL qssbti(3), bmrot(3, 3), v(3), vmult(3)
     643       33288 :       REAL :: gkrot(3, SIZE(lapw%gk, 2), 2)
     644       33288 :       REAL :: rph(SIZE(lapw%gk, 2), 2)
     645       33288 :       REAL :: cph(SIZE(lapw%gk, 2), 2)
     646       33288 :       COMPLEX ylm((atoms%lmaxd + 1)**2)
     647       33288 :       COMPLEX cwork(-2*atoms%llod:2*atoms%llod + 1, 2*(2*atoms%llod + 1), atoms%nlod, 2)
     648             :       !     ..
     649             :       !     .. Data statements ..
     650             :       REAL, PARAMETER :: eps = 1.0E-8
     651             : 
     652       33288 :       con1 = fpi_const/SQRT(cell%omtil)
     653       33288 :       nintsp = MERGE(2, 1, noco%l_ss)
     654       66576 :       DO iintsp = 1, nintsp
     655       33288 :          IF (iintsp .EQ. 1) THEN
     656      133152 :             qssbti = -nococonv%qss/2
     657             :          ELSE
     658           0 :             qssbti = +nococonv%qss/2
     659             :          ENDIF
     660             : 
     661             :          !--->    set up phase factors
     662     4130628 :          DO k = 1, lapw%nv(iintsp)
     663    16389360 :             th = tpi_const*DOT_PRODUCT((/lapw%k1(k, iintsp), lapw%k2(k, iintsp), lapw%k3(k, iintsp)/) + qssbti + lapw%qphon, atoms%taual(:, na))
     664     4097340 :             rph(k, iintsp) = COS(th)
     665     4130628 :             cph(k, iintsp) = -SIN(th)
     666             :          END DO
     667             : 
     668       66576 :          IF (np .EQ. 1) THEN
     669    16073352 :             gkrot(:, :, iintsp) = lapw%gk(:, :, iintsp)
     670             :          ELSE
     671       11648 :             bmrot = MATMUL(1.*sym%mrot(:, :, np), cell%bmat)
     672       87492 :             DO k = 1, lapw%nv(iintsp)
     673             :                !-->           apply the rotation that brings this atom into the
     674             :                !-->           representative (this is the definition of ngopr(na))
     675             :                !-->           and transform to cartesian coordinates
     676      349072 :                v(:) = lapw%vk(:, k, iintsp)
     677     1396512 :                gkrot(:, k, iintsp) = MATMUL(v, bmrot)
     678             :             END DO
     679             :          END IF
     680             :          !--->   end loop over interstitial spin
     681             :       ENDDO
     682             : 
     683      194328 :       nkvec(:, :) = 0
     684     4402536 :       cwork(:, :, :, :) = CMPLX(0.0, 0.0)
     685       33288 :       enough = .FALSE.
     686             : 
     687             :       ! Typically the search for linearly independent K vectors starts from the
     688             :       ! top and then goes down. This seems to be more stable than the opposite
     689             :       ! direction. The exception is a calculation with a spin spiral. There it
     690             :       ! is important that the K vectors for both spins feature the same G
     691             :       ! vectors. When starting from the bottom this is typically simple to check
     692             :       ! by just comparing the indices of the two K vectors. From the top this
     693             :       ! would require a more elaborate comparison.
     694             : 
     695       33288 :       nApproachEnd = 8
     696             : 
     697       33288 :       IF (noco%l_ss) THEN
     698           0 :          nApproachEnd = 4
     699             :       END IF
     700             : 
     701       70434 :       DO lo = 1, atoms%nlo(ntype)
     702             :          enough = .FALSE.
     703             :          nApproach = 0
     704      107580 :          DO WHILE (.NOT.enough)
     705       37146 :             nApproach = nApproach + 1
     706      111438 :             nkvec(lo, :) = 0
     707     3468342 :             cwork(:,:,lo,:) = CMPLX(0.0, 0.0)
     708       37146 :             k_start = 2
     709       37146 :             k_end = 2
     710       37146 :             increment = 1
     711       37146 :             SELECT CASE (nApproach)
     712             :                CASE (1)
     713       37146 :                   k_start = 2  ! avoid k=1 !!! GB16
     714       37146 :                   k_end = MIN(lapw%nv(1), lapw%nv(nintsp))
     715       37146 :                   increment = 1
     716       37146 :                   l_norm = .FALSE.
     717       37146 :                   linindqStart = 1.0e-6
     718       37146 :                   linindqEnd = 1.0e-6
     719             :                CASE (2)
     720           0 :                   k_start = 2  ! avoid k=1 !!! GB16
     721           0 :                   k_end = MIN(lapw%nv(1), lapw%nv(nintsp))
     722           0 :                   increment = 1
     723           0 :                   l_norm = .TRUE.
     724           0 :                   linindqStart = 1.0e-6
     725           0 :                   linindqEnd = 1.0e-6
     726             :                CASE (3)
     727           0 :                   k_start = 2  ! avoid k=1 !!! GB16
     728           0 :                   k_end = MIN(lapw%nv(1), lapw%nv(nintsp))
     729           0 :                   increment = 1
     730           0 :                   l_norm = .FALSE.
     731           0 :                   linindqStart = 2.0e-5
     732           0 :                   linindqEnd = 1.0e-6
     733             :                CASE (4)
     734           0 :                   k_start = 2  ! avoid k=1 !!! GB16
     735           0 :                   k_end = MIN(lapw%nv(1), lapw%nv(nintsp))
     736           0 :                   increment = 1
     737           0 :                   l_norm = .TRUE.
     738           0 :                   linindqStart = 2.0e-5
     739           0 :                   linindqEnd = 1.0e-6
     740             :                CASE (5)
     741           0 :                   k_start = MIN(lapw%nv(1), lapw%nv(nintsp))
     742           0 :                   k_end = 1
     743           0 :                   increment = -1
     744           0 :                   l_norm = .FALSE.
     745           0 :                   linindqStart = 1.0e-6
     746           0 :                   linindqEnd = 1.0e-6
     747             :                CASE (6)
     748           0 :                   k_start = MIN(lapw%nv(1), lapw%nv(nintsp))
     749           0 :                   k_end = 1
     750           0 :                   increment = -1
     751           0 :                   l_norm = .TRUE.
     752           0 :                   linindqStart = 1.0e-6
     753           0 :                   linindqEnd = 1.0e-6
     754             :                CASE (7)
     755           0 :                   k_start = MIN(lapw%nv(1), lapw%nv(nintsp))
     756           0 :                   k_end = 1
     757           0 :                   increment = -1
     758           0 :                   l_norm = .FALSE.
     759           0 :                   linindqStart = 2.0e-5
     760           0 :                   linindqEnd = 1.0e-6
     761             :                CASE (8)
     762           0 :                   k_start = MIN(lapw%nv(1), lapw%nv(nintsp))
     763           0 :                   k_end = 1
     764           0 :                   increment = -1
     765           0 :                   l_norm = .TRUE.
     766           0 :                   linindqStart = 2.0e-5
     767       37146 :                   linindqEnd = 1.0e-6
     768             :             END SELECT
     769       37146 :             k = k_start - increment
     770       37146 :             DO WHILE (.NOT.enough)
     771       98442 :                enough = .FALSE.
     772       98442 :                k = k + increment
     773       98442 :                IF ((k.GT.MAX(k_start,k_end)).OR.(k.LT.MIN(k_start,k_end))) THEN
     774           0 :                   IF (nApproach.GT.nApproachEnd) THEN
     775           0 :                      WRITE (oUnit, FMT=*) 'vec_for_lo did not find enough linearly independent'
     776           0 :                      WRITE (oUnit, FMT=*) 'clo coefficient-vectors. the linear independence'
     777           0 :                      WRITE (oUnit, FMT=*) 'quality, linindq, is set: ', linindqEnd
     778           0 :                      WRITE (oUnit, FMT=*) 'this value might be to large.'
     779           0 :                      WRITE (*, *) 'Atom: ', na, 'type: ', ntype, 'nv: ', lapw%nv, 'lo: ', lo, 'l: ', atoms%llo(lo, ntype)
     780           0 :                      CALL juDFT_error("not enough lin. indep. clo-vectors", calledby="vec_for_lo")
     781             :                   END IF
     782             :                   EXIT
     783             :                END IF
     784             : 
     785      196884 :                DO iintsp = 1, nintsp
     786             : 
     787      393768 :                   vmult(:) = gkrot(:, k, iintsp)
     788       98442 :                   CALL ylm4(atoms%lnonsph(ntype), vmult, ylm)
     789       98442 :                   l_lo1 = .false.
     790       98442 :                   IF ((lapw%rk(k, iintsp) .LT. eps) .AND. (.not. noco%l_ss)) THEN
     791             :                      l_lo1 = .true.
     792             :                   ELSE
     793       98442 :                      l_lo1 = .false.
     794             :                   END IF
     795             : 
     796       98442 :                   numK = MERGE(2*atoms%llo(lo, ntype)+1,2*(2*atoms%llo(lo, ntype)+1),sym%invsat(na).EQ.0)
     797             : 
     798      196884 :                   IF (l_lo1) THEN
     799           0 :                      IF ((nkvec(lo, iintsp) .EQ. 0) .AND. (atoms%llo(lo, ntype) .EQ. 0)) THEN
     800           0 :                         nkvec(lo, iintsp) = 1
     801           0 :                         lapw%kvec(nkvec(lo, iintsp), lo, na) = k
     802           0 :                         term1 = con1 * (atoms%rmt(ntype)**2.0) / 2.0
     803           0 :                         cwork(0, 1, lo, iintsp) = term1 / sqrt(2.0*tpi_const)
     804           0 :                         norm = DOT_PRODUCT(cwork(0:0, 1, lo, iintsp),cwork(0:0, 1, lo, iintsp))
     805           0 :                         IF (l_norm) cwork(0:0, 1, lo, iintsp) = cwork(0:0, 1, lo, iintsp) / SQRT(norm)
     806           0 :                         IF ((sym%invsat(na) .EQ. 1) .OR. (sym%invsat(na) .EQ. 2)) THEN
     807           0 :                            cwork(1, 1, lo, iintsp) = conjg(term1) / sqrt(2.0*tpi_const)
     808             :                         END IF
     809             :                      
     810             :                      END IF
     811             :                   ELSE
     812       98442 :                      term1 = con1 * ((atoms%rmt(ntype)**2.0)/2.0) * CMPLX(rph(k, iintsp), cph(k, iintsp))
     813       98442 :                      IF (sym%invsat(na) .EQ. 0) THEN
     814       97240 :                         IF ((nkvec(lo, iintsp)) .LT. (2*atoms%llo(lo, ntype) + 1)) THEN
     815       97240 :                            nkvec(lo, iintsp) = nkvec(lo, iintsp) + 1
     816       97240 :                            l = atoms%llo(lo, ntype)
     817       97240 :                            ll1 = l*(l + 1) + 1
     818      389724 :                            DO m = -l, l
     819      292484 :                               lm = ll1 + m
     820      389724 :                               cwork(m, nkvec(lo, iintsp), lo, iintsp) = term1 * ylm(lm)
     821             :                            END DO
     822      389724 :                            norm = DOT_PRODUCT(cwork(-l:l, 1, lo, iintsp),cwork(-l:l, 1, lo, iintsp))
     823       97240 :                            IF (l_norm) cwork(-l:l, 1, lo, iintsp) = cwork(-l:l, 1, lo, iintsp) / SQRT(norm)
     824       97240 :                            stepSize = (REAL(linindqStart - linindqEnd)) / numK
     825       97240 :                            linindq = (numK - (REAL(nkvec(lo, iintsp) + 1))) * stepSize + linindqEnd
     826       97240 :                            CALL orthoglo(input%l_real, atoms, nkvec(lo, iintsp), lo, l, linindq, .FALSE., cwork(-2*atoms%llod, 1, 1, iintsp), linind)
     827       97240 :                            IF (linind) THEN
     828       87752 :                               lapw%kvec(nkvec(lo, iintsp), lo, na) = k
     829             :                            ELSE
     830        9488 :                               nkvec(lo, iintsp) = nkvec(lo, iintsp) - 1
     831             :                            END IF
     832             :                         END IF
     833        1202 :                      ELSE IF ((sym%invsat(na) .EQ. 1) .OR. (sym%invsat(na) .EQ. 2)) THEN
     834        1202 :                         IF (nkvec(lo, iintsp) .LT. 2*(2*atoms%llo(lo, ntype) + 1)) THEN
     835        1202 :                            nkvec(lo, iintsp) = nkvec(lo, iintsp) + 1
     836        1202 :                            l = atoms%llo(lo, ntype)
     837        1202 :                            ll1 = l*(l + 1) + 1
     838        4264 :                            DO m = -l, l
     839        3062 :                               lm = ll1 + m
     840        3062 :                               mind = -l + m
     841        3062 :                               cwork(mind, nkvec(lo, iintsp), lo, iintsp) = term1 * ylm(lm)
     842        3062 :                               mind = l + 1 + m
     843        3062 :                               lmp = ll1 - m
     844        4264 :                               cwork(mind, nkvec(lo, iintsp), lo, iintsp) = ((-1)**(l + m))*CONJG(term1*ylm(lmp))
     845             :                            END DO
     846        1202 :                            minIndex = -l - l
     847        1202 :                            maxIndex = l + l + 1
     848        7326 :                            norm = DOT_PRODUCT(cwork(-2*l:2*l+1, 1, lo, iintsp),cwork(-2*l:2*l+1, 1, lo, iintsp))
     849        1202 :                            IF (l_norm) cwork(-2*l:2*l+1, 1, lo, iintsp) = cwork(-2*l:2*l+1, 1, lo, iintsp) / SQRT(norm)
     850        1202 :                            stepSize = (REAL(linindqStart - linindqEnd)) / numK
     851        1202 :                            linindq = (numK - (REAL(nkvec(lo, iintsp) + 1))) * stepSize + linindqEnd
     852        1202 :                            CALL orthoglo(input%l_real, atoms, nkvec(lo, iintsp), lo, l, linindq, .TRUE., cwork(-2*atoms%llod, 1, 1, iintsp), linind)
     853        1202 :                            IF (linind) THEN
     854        1172 :                               lapw%kvec(nkvec(lo, iintsp), lo, na) = k
     855             :                            ELSE
     856          30 :                               nkvec(lo, iintsp) = nkvec(lo, iintsp) - 1
     857             :                            END IF
     858             :                         END IF
     859             : 
     860             :                      END IF
     861             :                   END IF
     862             :                END DO
     863             : 
     864       98442 :                enough = .TRUE.
     865       98442 :                IF (nkvec(lo, 1) .EQ. nkvec(lo, nintsp)) THEN   ! k-vec accepted by both spin channels
     866       98442 :                   IF (sym%invsat(na) .EQ. 0) THEN
     867       97240 :                      IF (nkvec(lo, 1) .LT. (2*atoms%llo(lo, ntype) + 1)) THEN
     868             :                         enough = .FALSE.
     869             :                      ENDIF
     870             :                   ELSE
     871        1202 :                      IF (nkvec(lo, 1) .LT. (2*(2*atoms%llo(lo, ntype) + 1))) THEN
     872             :                         enough = .FALSE.
     873             :                      ENDIF
     874             :                   ENDIF
     875             :                ELSE
     876           0 :                   nkmin = MIN(nkvec(lo, 1), nkvec(lo, nintsp)) ! try another k-vec
     877           0 :                   nkvec(lo, 1) = nkmin
     878           0 :                   nkvec(lo, nintsp) = nkmin
     879           0 :                   enough = .FALSE.
     880             :                END IF
     881             :             END DO
     882             :          END DO
     883             :       END DO
     884             : 
     885       70434 :       lapw%nkvec(:atoms%nlo(ntype), na) = nkvec(:atoms%nlo(ntype), 1)
     886             : 
     887       33288 :    END SUBROUTINE priv_vec_for_lo
     888             : 
     889           0 : END MODULE m_types_lapw

Generated by: LCOV version 1.14