LCOV - code coverage report
Current view: top level - types - types_lapw.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 206 241 85.5 %
Date: 2019-09-08 04:53:50 Functions: 5 7 71.4 %

          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             :   PRIVATE
      10             :   TYPE t_lapw
      11             :      INTEGER :: nv(2)
      12             :      INTEGER :: num_local_cols(2)
      13             :      INTEGER :: nv_tot
      14             :      INTEGER :: nmat
      15             :      INTEGER :: nlotot
      16             :      INTEGER,ALLOCATABLE:: k1(:,:)
      17             :      INTEGER,ALLOCATABLE:: k2(:,:)
      18             :      INTEGER,ALLOCATABLE:: k3(:,:)
      19             :      INTEGER,ALLOCATABLE CPP_MANAGED:: gvec(:,:,:) !replaces k1,k2,k3
      20             :      INTEGER,ALLOCATABLE:: kp(:,:)
      21             :      REAL,ALLOCATABLE::rk(:,:)
      22             :      REAL,ALLOCATABLE CPP_MANAGED::gk(:,:,:)
      23             :      REAL,ALLOCATABLE::vk(:,:,:)
      24             :      INTEGER,ALLOCATABLE::matind(:,:)
      25             :      INTEGER,ALLOCATABLE::index_lo(:,:)
      26             :      INTEGER,ALLOCATABLE::kvec(:,:,:)
      27             :      INTEGER,ALLOCATABLE::nkvec(:,:)
      28             :      REAL   :: bkpt(3)
      29             :    CONTAINS
      30             :      PROCEDURE,PASS :: init =>lapw_init
      31             :      PROCEDURE,PASS :: alloc =>lapw_alloc
      32             :      PROCEDURE,PASS :: phase_factors =>lapw_phase_factors
      33             :   END TYPE t_lapw
      34             :   PUBLIC :: t_lapw
      35             : 
      36             : 
      37             : CONTAINS
      38        3080 :   SUBROUTINE lapw_alloc(lapw,cell,input,noco)
      39             :     !
      40             :     !*********************************************************************
      41             :     !     determines dimensions of the lapw basis set with |k+G|<rkmax.
      42             :     !     bkpt is the k-point given in internal units
      43             :     !*********************************************************************
      44             :     USE m_boxdim
      45             :     USE m_types_setup
      46             : 
      47             :     IMPLICIT NONE
      48             :     TYPE(t_cell),INTENT(IN)      :: cell
      49             :     TYPE(t_input),INTENT(IN)     :: input
      50             :     TYPE(t_noco),INTENT(IN)      :: noco
      51             :     CLASS(t_lapw),INTENT(INOUT)  :: lapw
      52             : 
      53             : 
      54             :     INTEGER j1,j2,j3,mk1,mk2,mk3,nv,addX,addY,addZ
      55             :     INTEGER ispin,nvh(2)
      56             : 
      57             :     REAL arltv1,arltv2,arltv3,rkm,rk2,r2,s(3)
      58             :     ! ..
      59             :     !
      60             :     !------->          ABBREVIATIONS
      61             :     !
      62             : 
      63             :     !   rkmax       : cut-off for |g+k|
      64             :     !   arltv(i)    : length of reciprical lattice vector along
      65             :     !                 direction (i)
      66             :     !
      67             :     !---> Determine rkmax box of size mk1, mk2, mk3,
      68             :     !     for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
      69             :     !
      70        3080 :     CALL boxdim(cell%bmat,arltv1,arltv2,arltv3)
      71             : 
      72             :     !     (add 1+1 due to integer rounding, strange k_vector in BZ)
      73        3080 :     mk1 = int(input%rkmax/arltv1)+2
      74        3080 :     mk2 = int(input%rkmax/arltv2)+2
      75        3080 :     mk3 = int(input%rkmax/arltv3)+2
      76             : 
      77        3080 :     rkm = input%rkmax
      78        3080 :     rk2 = rkm*rkm
      79             :     !---> obtain vectors
      80             :     !---> in a spin-spiral calculation different basis sets are used for
      81             :     !---> the two spin directions, because the cutoff radius is defined
      82             :     !---> by |G + k +/- qss/2| < rkmax.
      83        3080 :     nvh(2)=0
      84        6280 :     DO ispin = 1,MERGE(2,1,noco%l_ss)
      85        3200 :        addX = abs(NINT((lapw%bkpt(1)+ (2*ispin - 3)/2.0*noco%qss(1))/arltv1))
      86        3200 :        addY = abs(NINT((lapw%bkpt(2)+ (2*ispin - 3)/2.0*noco%qss(2))/arltv2))
      87        3200 :        addZ = abs(NINT((lapw%bkpt(3)+ (2*ispin - 3)/2.0*noco%qss(3))/arltv3))
      88        3200 :        nv = 0
      89       33960 :        DO  j1 = -mk1-addX,mk1+addX
      90      335096 :           DO  j2 = -mk2-addY,mk2+addY
      91     4977036 :              DO  j3 = -mk3-addZ,mk3+addZ
      92     4645140 :                 s = lapw%bkpt + (/j1,j2,j3/) + (2*ispin - 3)/2.0*noco%qss
      93     4645140 :                 r2 = dot_PRODUCT(MATMUL(s,cell%bbmat),s)
      94     4946276 :                 IF (r2.LE.rk2)  nv = nv + 1
      95             :              END DO
      96             :           END DO
      97             :        END DO
      98        6280 :        nvh(ispin)  = nv
      99             :     END DO
     100        3080 :     nv  = MAX(nvh(1),nvh(2))
     101             : 
     102        3080 :     IF (ALLOCATED(lapw%rk)) THEN
     103        2118 :        IF (SIZE(lapw%rk)==nv) THEN
     104          60 :           RETURN !
     105             :        ELSE
     106        2058 :           DEALLOCATE(lapw%rk,lapw%gvec,lapw%vk,lapw%gk,lapw%matind)
     107        2058 :           DEALLOCATE(lapw%k1,lapw%k2,lapw%k3)
     108             :        ENDIF
     109             :     ENDIF
     110        3020 :     ALLOCATE(lapw%rk(nv,input%jspins) )
     111        3020 :     ALLOCATE(lapw%gvec(3,nv,input%jspins))
     112        3020 :     ALLOCATE(lapw%vk(3,nv,input%jspins))
     113        3020 :     ALLOCATE(lapw%gk(3,nv,input%jspins))
     114        3020 :     ALLOCATE(lapw%k1(nv,input%jspins)) !shpuld be removed
     115        3020 :     ALLOCATE(lapw%k2(nv,input%jspins)) !
     116        3020 :     ALLOCATE(lapw%k3(nv,input%jspins)) !
     117        3020 :     ALLOCATE(lapw%matind(nv,2))
     118             : 
     119        8596 :     lapw%rk = 0 ; lapw%gvec = 0 ;lapw%nv=0
     120             :   END SUBROUTINE lapw_alloc
     121             : 
     122             : 
     123        3080 :   SUBROUTINE lapw_init(lapw,input,noco,kpts,atoms,sym,&
     124             :        nk,cell,l_zref,mpi)
     125             :     USE m_types_mpi
     126             :     USE m_sort
     127             :     USE m_boxdim
     128             :     USE m_types_setup
     129             :     USE m_types_kpts
     130             :     IMPLICIT NONE
     131             : 
     132             :     TYPE(t_input),INTENT(IN)       :: input
     133             :     TYPE(t_noco),INTENT(IN)        :: noco
     134             :     TYPE(t_cell),INTENT(IN)        :: cell
     135             :     TYPE(t_atoms),INTENT(IN)       :: atoms
     136             :     TYPE(t_sym),INTENT(IN)         :: sym
     137             :     TYPE(t_kpts),INTENT(IN)        :: kpts
     138             :     TYPE(t_mpi),INTENT(IN),OPTIONAL:: mpi
     139             :     CLASS(t_lapw),INTENT(INOUT)    :: lapw
     140             :     !     .. 
     141             :     !     .. Scalar Arguments ..
     142             :     INTEGER, INTENT  (IN) :: nk
     143             :     LOGICAL, INTENT (IN)  :: l_zref
     144             :     !     ..
     145             :     !     .. Array Arguments ..
     146             :     !     ..
     147             :     !     .. Local Scalars ..
     148             :     REAL arltv1,arltv2,arltv3,r2,rk2,rkm,r2q,gla,eps,t
     149             :     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
     150             :     !     ..
     151             :     !     .. Local Arrays ..
     152             :     REAL                :: s(3),sq(3)
     153        3080 :     REAL,ALLOCATABLE    :: rk(:),rkq(:),rkqq(:)
     154        3080 :     INTEGER,ALLOCATABLE :: gvec(:,:),index3(:)
     155             : 
     156             : 
     157             :     !     ..
     158             :     !---> in a spin-spiral calculation different basis sets are used for
     159             :     !---> the two spin directions, because the cutoff radius is defined
     160             :     !---> by |G + k +/- qss/2| < rkmax.
     161             : 
     162        3080 :     IF (nk>kpts%nkpt) THEN
     163           0 :        lapw%bkpt(:)=kpts%bkf(:,nk)
     164             :     ELSE
     165        3080 :        lapw%bkpt(:) = kpts%bk(:,nk)
     166             :     ENDIF
     167             : 
     168             : 
     169        3080 :     CALL lapw%alloc(cell,input,noco)
     170             : 
     171        3080 :     ALLOCATE(gvec(3,SIZE(lapw%gvec,2)))
     172        3080 :     ALLOCATE(rk(SIZE(lapw%gvec,2)),rkq(SIZE(lapw%gvec,2)),rkqq(SIZE(lapw%gvec,2)))
     173        3080 :     ALLOCATE(index3(SIZE(lapw%gvec,2)))
     174             : 
     175             : 
     176             :     !---> Determine rkmax box of size mk1, mk2, mk3,
     177             :     !     for which |G(mk1,mk2,mk3) + (k1,k2,k3)| < rkmax
     178             :     !     arltv(i) length of reciprical lattice vector along direction (i)
     179             :     !
     180        3080 :     CALL boxdim(cell%bmat,arltv1,arltv2,arltv3)
     181             : 
     182             :     !     (add 1+1 due to integer rounding, strange k_vector in BZ)
     183        3080 :     mk1 = int( input%rkmax/arltv1 )+4
     184        3080 :     mk2 = int( input%rkmax/arltv2 )+4
     185        3080 :     mk3 = int( input%rkmax/arltv3 )+4
     186             : 
     187        3080 :     rk2 = input%rkmax*input%rkmax
     188             :     !---> if too many basis functions, reduce rkmax
     189        3844 :     spinloop:DO ispin = 1,input%jspins
     190        3200 :        addX = abs(NINT((lapw%bkpt(1)+ (2*ispin - 3)/2.0*noco%qss(1))/arltv1))
     191        3200 :        addY = abs(NINT((lapw%bkpt(2)+ (2*ispin - 3)/2.0*noco%qss(2))/arltv2))
     192        3200 :        addZ = abs(NINT((lapw%bkpt(3)+ (2*ispin - 3)/2.0*noco%qss(3))/arltv3))
     193             :        !--->    obtain vectors
     194        3200 :        n = 0
     195       46760 :        DO  j1 = -mk1-addX,mk1+addX
     196      645272 :           DO  j2 = -mk2-addY,mk2+addY
     197    12064236 :              DO  j3 = -mk3-addZ,mk3+addZ
     198    11422164 :                 s=lapw%bkpt+(/j1,j2,j3/)+(2*ispin - 3)/2.0*noco%qss
     199    11422164 :                 sq = lapw%bkpt+ (/j1,j2,j3/)
     200    11422164 :                 r2 = dot_PRODUCT(s,MATMUL(s,cell%bbmat))
     201    11422164 :                 r2q = dot_PRODUCT(sq,MATMUL(sq,cell%bbmat))
     202    12020676 :                 IF (r2.LE.rk2) THEN
     203      559601 :                    n = n + 1
     204      559601 :                    gvec(:,n) = (/j1,j2,j3/)
     205      559601 :                    rk(n) = SQRT(r2)
     206      559601 :                    rkq(n) = SQRT(r2q)
     207             :                 END IF
     208             :              ENDDO
     209             :           ENDDO
     210             :        ENDDO
     211        3200 :        lapw%nv(ispin) = n
     212             : 
     213             :        !Sort according to k+g, first construct secondary sort key
     214      562801 :        DO  k = 1,lapw%nv(ispin)
     215             :           rkqq(k) = (mk1+gvec(1,k)) + (mk2+gvec(2,k))*(2*mk1+1) +&
     216      562801 :                (mk3+gvec(3,k))*(2*mk1+1)*(2*mk2+1)
     217             :        ENDDO
     218        3200 :        CALL sort(index3(:lapw%nv(ispin)),rkq,rkqq)
     219      562801 :        DO n=1,lapw%nv(ispin)
     220      559601 :           lapw%gvec(:,n,ispin) = gvec(:,index3(n))
     221      562801 :           lapw%rk(n,ispin)     = rk(index3(n))
     222             :        ENDDO
     223             :        !--->    determine pairs of K-vectors, where K_z = K'_-z to use 
     224             :        !--->    z-reflection
     225        3200 :        IF (l_zref) THEN
     226          68 :           n=0
     227        8638 :           DO i=1,lapw%nv(ispin)
     228     1257956 :              DO j=1,i
     229      624659 :                 IF (ALL(lapw%gvec(1:2,i,ispin).EQ.lapw%gvec(1:2,j,ispin)).AND.&
     230        8570 :                      (lapw%gvec(3,i,ispin).EQ.-lapw%gvec(3,j,ispin))) THEN
     231        5070 :                    n=n+1 
     232        5070 :                    lapw%matind(n,1)=i
     233        5070 :                    lapw%matind(n,2)=j
     234             :                 ENDIF
     235             :              ENDDO
     236             :           ENDDO
     237          68 :           nred=n
     238          68 :           IF (PRESENT(mpi)) THEN
     239          68 :              IF (mpi%n_size.GT.1) THEN
     240             :                 !
     241             :                 !--->     order K's in sequence K_1,...K_n | K_0,... | K_-1,....K_-n
     242             :                 !
     243           0 :                 n_inner = lapw%nv(ispin) - nred
     244           0 :                 IF (MOD(nred,mpi%n_size).EQ.0) THEN
     245             :                    n_bound = nred
     246             :                 ELSE
     247           0 :                    n_bound = (1+INT( nred/mpi%n_size ))*mpi%n_size
     248             :                 ENDIF
     249           0 :                 IF (lapw%nv(ispin) - nred + n_bound.GT.SIZE(lapw%gvec,2)) THEN
     250           0 :                    CALL juDFT_error("BUG:z-ref & ev || : dimension too small!" ,calledby ="types_lapw")
     251             :                 ENDIF
     252             : 
     253             :                 i = 1
     254             :                 j = 1
     255           0 :                 DO n = 1, nred 
     256           0 :                    IF (lapw%matind(n,1).EQ.lapw%matind(n,2)) THEN
     257           0 :                       index3(lapw%matind(n,1)) = n_inner + i
     258           0 :                       i = i + 1
     259             :                    ELSE
     260           0 :                       index3(lapw%matind(n,1)) = j
     261           0 :                       index3(lapw%matind(n,2)) = j + n_bound
     262           0 :                       j = j + 1
     263             :                    ENDIF
     264             :                 ENDDO
     265             :                 !--->          resort the rk,k1,k2,k3 and lapw%matind arrays:
     266           0 :                 DO n = 1, lapw%nv(ispin)
     267           0 :                    rk(n)  = lapw%rk(n,ispin)
     268           0 :                    gvec(:,n) = lapw%gvec(:,n,ispin)
     269             :                 ENDDO
     270           0 :                 DO n = lapw%nv(ispin), 1, -1
     271           0 :                    lapw%rk(index3(n),ispin) = rk(n)
     272           0 :                    lapw%gvec(:,index3(n),ispin) = gvec(:,n)
     273             :                 ENDDO
     274           0 :                 DO n = nred + 1, n_bound
     275           0 :                    lapw%rk(n,ispin) = lapw%rk(lapw%nv(ispin),ispin)
     276           0 :                    lapw%gvec(:,n,ispin) = lapw%gvec(:,lapw%nv(ispin),ispin)
     277             :                 ENDDO
     278           0 :                 lapw%nv(ispin) = lapw%nv(ispin) - nred + n_bound
     279             :              ENDIF
     280             :           ENDIF
     281             :        ENDIF
     282      562801 :        DO k=1,lapw%nv(ispin)
     283      559601 :           lapw%vk(:,k,ispin)=lapw%bkpt+lapw%gvec(:,k,ispin)+(ispin-1.5)*noco%qss
     284      562801 :           lapw%gk(:,k,ispin)=MATMUL(TRANSPOSE(cell%bmat),lapw%vk(:,k,ispin))/MAX (lapw%rk(k,ispin),1.0e-30)
     285             :        ENDDO
     286             : 
     287             : 
     288        3844 :        IF (.NOT.noco%l_ss.AND.input%jspins==2) THEN
     289             :           !Second spin is the same
     290        2436 :           lapw%nv(2)=lapw%nv(1)
     291        2436 :           lapw%gvec(:,:,2)=lapw%gvec(:,:,1)
     292        2436 :           lapw%rk(:,2)=lapw%rk(:,1)
     293        2436 :           lapw%vk(:,:,2)=lapw%vk(:,:,1)
     294        2436 :           lapw%gk(:,:,2)=lapw%gk(:,:,1)
     295             :           EXIT spinloop
     296             :        END IF
     297             : 
     298             :     ENDDO spinloop
     299             :     !should be removed later...
     300        3080 :     lapw%k1=lapw%gvec(1,:,:)
     301        3080 :     lapw%k2=lapw%gvec(2,:,:)
     302        3080 :     lapw%k3=lapw%gvec(3,:,:)
     303             : 
     304             :     !Count No of lapw distributed to this PE
     305        3080 :     lapw%num_local_cols=0
     306        8716 :     DO ispin=1,input%jspins
     307        8716 :        IF (PRESENT(mpi)) THEN
     308      575266 :           DO k=mpi%n_rank+1,lapw%nv(ispin),mpi%n_size
     309        5636 :              lapw%num_local_cols(ispin)=lapw%num_local_cols(ispin)+1
     310             :           END DO
     311             :        ELSE
     312           0 :           lapw%num_local_cols(ispin) = lapw%nv(ispin)
     313             :        END IF
     314             :     END DO
     315             : 
     316        3080 :     IF (ANY(atoms%nlo>0)) CALL priv_lo_basis_setup(lapw,atoms,sym,noco,cell)
     317             : 
     318        3080 :     lapw%nv_tot=lapw%nv(1)
     319        3080 :     lapw%nmat=lapw%nv(1)+atoms%nlotot
     320        3080 :     IF (noco%l_noco) lapw%nv_tot=lapw%nv_tot+lapw%nv(2)
     321        3080 :     IF (noco%l_noco) lapw%nmat=lapw%nv_tot+2*atoms%nlotot
     322             : 
     323             :   CONTAINS
     324             : 
     325         486 :     SUBROUTINE priv_lo_basis_setup(lapw,atoms,sym,noco,cell)
     326             :       USE m_types_setup
     327             : 
     328             :       IMPLICIT NONE
     329             :       TYPE(t_lapw),INTENT(INOUT):: lapw
     330             :       TYPE(t_atoms),INTENT(IN)  :: atoms
     331             :       TYPE(t_sym),INTENT(IN)    :: sym
     332             :       TYPE(t_cell),INTENT(IN)   :: cell
     333             :       TYPE(t_noco),INTENT(IN)   :: noco
     334             : 
     335             : 
     336             :       INTEGER:: n,na,nn,np,lo,nkvec_sv,nkvec(atoms%nlod,2),iindex
     337         486 :       IF (.NOT.ALLOCATED(lapw%kvec)) THEN
     338          94 :          ALLOCATE(lapw%kvec(2*(2*atoms%llod+1),atoms%nlod,atoms%nat))
     339          94 :          ALLOCATE(lapw%nkvec(atoms%nlod,atoms%nat))
     340          94 :          ALLOCATE(lapw%index_lo(atoms%nlod,atoms%nat))
     341             :       ENDIF
     342         486 :       iindex=0
     343         486 :       na=0
     344         486 :       nkvec_sv=0
     345        2544 :       DO n=1,atoms%ntype
     346        6086 :          DO nn=1,atoms%neq(n)
     347        3542 :             na=na+1
     348        3542 :             if (atoms%invsat(na)>1) cycle
     349             :             !np = MERGE(oneD%ods%ngopr(na),sym%invtab(atoms%ngopr(na)),oneD%odi%d1)
     350        2078 :             np=sym%invtab(atoms%ngopr(na))
     351        2078 :             CALL priv_vec_for_lo(atoms,sym,na,n,np,noco,lapw,cell)
     352        7978 :             DO lo = 1,atoms%nlo(n)
     353        3842 :                lapw%index_lo(lo,na)=iindex
     354        7384 :                iindex=iindex+lapw%nkvec(lo,na)
     355             :             ENDDO
     356             :          ENDDO
     357             :       ENDDO
     358         486 :     END SUBROUTINE priv_lo_basis_setup
     359             : 
     360             :   END SUBROUTINE lapw_init
     361             : 
     362             : 
     363        6620 :   SUBROUTINE lapw_phase_factors(lapw,iintsp,tau,qss,cph)
     364             :     USE m_constants
     365             :     USE m_types_setup
     366             :     IMPLICIT NONE
     367             :     CLASS(t_lapw),INTENT(in):: lapw
     368             :     INTEGER,INTENT(IN)     :: iintsp
     369             :     REAL,INTENT(in)        :: tau(3),qss(3)
     370             :     COMPLEX,INTENT(out)    :: cph(:)
     371             : 
     372             :     INTEGER:: k
     373             :     REAL:: th
     374     2048884 :     DO k = 1,lapw%nv(iintsp)
     375     2042264 :        th= DOT_PRODUCT(lapw%gvec(:,k,iintsp)+(iintsp-1.5)*qss,tau)
     376     2048884 :        cph(k) = CMPLX(COS(tpi_const*th),-SIN(tpi_const*th))
     377             :     END DO
     378        6620 :   END SUBROUTINE lapw_phase_factors
     379             : 
     380             :   
     381        2078 :   SUBROUTINE priv_vec_for_lo(atoms,sym,na,&
     382             :        n,np,noco, lapw,cell)
     383             :     USE m_constants,ONLY: tpi_const,fpi_const
     384             :     USE m_orthoglo
     385             :     USE m_ylm
     386             :     USE m_types_setup
     387             :     IMPLICIT NONE
     388             :     TYPE(t_noco),INTENT(IN)   :: noco
     389             :     TYPE(t_sym),INTENT(IN)    :: sym
     390             :     TYPE(t_cell),INTENT(IN)   :: cell
     391             :     TYPE(t_atoms),INTENT(IN)  :: atoms
     392             :     TYPE(t_lapw),INTENT(INOUT):: lapw
     393             :     !     ..
     394             :     !     .. Scalar Arguments ..
     395             :     INTEGER, INTENT (IN) :: na,n,np 
     396             :     !     ..
     397             :     !     .. Array Arguments ..
     398             :     !     ..
     399             :     !     .. Local Scalars ..
     400             :     COMPLEX term1 
     401             :     REAL th,con1
     402             :     INTEGER l,lo ,mind,ll1,lm,iintsp,k,nkmin,ntyp,lmp,m,nintsp
     403             :     LOGICAL linind,enough,l_lo1
     404             :     !     ..
     405             :     !     .. Local Arrays ..
     406        4156 :     INTEGER :: nkvec(atoms%nlod,2)
     407             :     REAL qssbti(3),bmrot(3,3),v(3),vmult(3)
     408        4156 :     REAL :: gkrot(3,SIZE(lapw%gk,2),2)
     409        4156 :     REAL :: rph(SIZE(lapw%gk,2),2)
     410        4156 :     REAL :: cph(SIZE(lapw%gk,2),2)
     411        4156 :     COMPLEX ylm( (atoms%lmaxd+1)**2 )
     412        4156 :     COMPLEX cwork(-2*atoms%llod:2*atoms%llod+1,2*(2*atoms%llod+1),atoms%nlod ,2)
     413             :     !     ..
     414             :     !     .. Data statements ..
     415             :     REAL, PARAMETER :: eps = 1.0E-8
     416             :     REAL, PARAMETER :: linindq = 1.0e-4
     417             : 
     418        2078 :     con1=fpi_const/SQRT(cell%omtil)
     419        2078 :     ntyp = n
     420        2078 :     nintsp=MERGE(2,1,noco%l_ss)
     421        4156 :     DO iintsp = 1,nintsp
     422        2078 :        IF (iintsp.EQ.1) THEN
     423        8312 :           qssbti = - noco%qss/2
     424             :        ELSE
     425           0 :           qssbti = + noco%qss/2
     426             :        ENDIF
     427             : 
     428             :        !--->    set up phase factors
     429     1133880 :        DO k = 1,lapw%nv(iintsp)
     430     1131802 :           th= tpi_const*DOT_PRODUCT((/lapw%k1(k,iintsp),lapw%k2(k,iintsp),lapw%k3(k,iintsp)/)+qssbti,atoms%taual(:,na))
     431     1131802 :           rph(k,iintsp) = COS(th)
     432     1133880 :           cph(k,iintsp) = -SIN(th)
     433             :        END DO
     434             : 
     435        4156 :        IF (np.EQ.1) THEN
     436        2058 :           gkrot(:,:,iintsp)=lapw%gk(:,:,iintsp)
     437             :        ELSE
     438          20 :           bmrot=MATMUL(1.*sym%mrot(:,:,np),cell%bmat)
     439       13820 :           DO k = 1,lapw%nv(iintsp)
     440             :              !-->           apply the rotation that brings this atom into the
     441             :              !-->           representative (this is the definition of ngopr(na))
     442             :              !-->           and transform to cartesian coordinates
     443       13800 :              v(:) = lapw%vk(:,k,iintsp)
     444       13820 :              gkrot(:,k,iintsp) = MATMUL(v,bmrot)
     445             :           END DO
     446             :        END IF
     447             :        !--->   end loop over interstitial spin
     448             :     ENDDO
     449             : 
     450        6234 :     nkvec(:,:) = 0
     451        6234 :     cwork(:,:,:,:) = CMPLX(0.0,0.0)
     452        2078 :     enough=.FALSE.
     453       15164 :     DO k = 1,MIN(lapw%nv(1),lapw%nv(nintsp))
     454       15164 :        IF (ANY(lapw%rk(k,:nintsp).LT.eps)) CYCLE
     455       15154 :        IF (.NOT.enough) THEN
     456       45462 :           DO iintsp = 1,nintsp
     457             : 
     458             :              !-->        generate spherical harmonics
     459       15154 :              vmult(:) =  gkrot(:,k,iintsp)
     460       15154 :              CALL ylm4(atoms%lnonsph(ntyp),vmult, ylm)
     461       15154 :                 enough = .TRUE.
     462       15154 :                 term1 = con1* ((atoms%rmt(ntyp)**2)/2)* CMPLX(rph(k,iintsp),cph(k,iintsp))
     463       44720 :                 DO lo = 1,atoms%nlo(ntyp)
     464       44720 :                    IF (atoms%invsat(na).EQ.0) THEN
     465        5230 :                       IF ((nkvec(lo,iintsp)).LT. (2*atoms%llo(lo,ntyp)+1)) THEN
     466        3286 :                          enough = .FALSE.
     467        3286 :                          nkvec(lo,iintsp) = nkvec(lo,iintsp) + 1
     468        3286 :                          l = atoms%llo(lo,ntyp)
     469        3286 :                          ll1 = l*(l+1) + 1
     470       12624 :                          DO m = -l,l
     471        9338 :                             lm = ll1 + m
     472       12624 :                             cwork(m,nkvec(lo,iintsp),lo,iintsp) = term1*ylm(lm)
     473             :                          END DO
     474             :                          CALL orthoglo(&
     475        3286 :                               sym%invs.and..not.noco%l_noco,atoms,nkvec(lo,iintsp),lo,l,linindq,.FALSE., cwork(-2*atoms%llod,1,1,iintsp),linind)
     476        3286 :                          IF (linind) THEN
     477        2062 :                             lapw%kvec(nkvec(lo,iintsp),lo,na) = k
     478             :                          ELSE
     479        1224 :                             nkvec(lo,iintsp) = nkvec(lo,iintsp) - 1
     480             :                          ENDIF
     481             :                       ENDIF
     482             :                    ELSE
     483       24336 :                       IF ((atoms%invsat(na).EQ.1) .OR. (atoms%invsat(na).EQ.2)) THEN
     484       24336 :                          IF (nkvec(lo,iintsp).LT.2*(2*atoms%llo(lo,ntyp)+1)) THEN
     485       15096 :                             enough = .FALSE.
     486       15096 :                             nkvec(lo,iintsp) = nkvec(lo,iintsp) + 1
     487       15096 :                             l = atoms%llo(lo,ntyp)
     488       15096 :                             ll1 = l*(l+1) + 1
     489       54528 :                             DO m = -l,l
     490       39432 :                                lm = ll1 + m
     491       39432 :                                mind = -l + m
     492       39432 :                                cwork(mind,nkvec(lo,iintsp),lo,iintsp) = term1*ylm(lm)
     493       39432 :                                mind = l + 1 + m
     494       39432 :                                lmp = ll1 - m
     495       54528 :                                cwork(mind,nkvec(lo,iintsp),lo,iintsp) = ((-1)** (l+m))*CONJG(term1*ylm(lmp))
     496             :                             END DO
     497             :                             CALL orthoglo(&
     498       15096 :                                  sym%invs.and..not.noco%l_noco,atoms,nkvec(lo,iintsp),lo,l,linindq,.TRUE., cwork(-2*atoms%llod,1,1,iintsp),linind)
     499       15096 :                             IF (linind) THEN
     500       11712 :                                lapw%kvec(nkvec(lo,iintsp),lo,na) = k
     501             :                                !                          write(*,*) nkvec(lo,iintsp),k,' <- '
     502             :                             ELSE
     503        3384 :                                nkvec(lo,iintsp) = nkvec(lo,iintsp) - 1
     504             :                             END IF
     505             :                          END IF
     506             :                       END IF
     507             :                    END IF
     508             :                 END DO
     509       30308 :                 IF ((k.EQ.lapw%nv(iintsp)) .AND. (.NOT.enough)) THEN
     510           0 :                    WRITE (6,FMT=*) 'vec_for_lo did not find enough linearly independent'
     511           0 :                    WRITE (6,FMT=*) 'clo coefficient-vectors. the linear independence'
     512           0 :                    WRITE (6,FMT=*) 'quality, linindq, is set: ',linindq
     513           0 :                    WRITE (6,FMT=*) 'this value might be to large.'
     514           0 :                    WRITE(*,*) na,k,lapw%nv 
     515           0 :                    CALL juDFT_error("not enough lin. indep. clo-vectors" ,calledby ="vec_for_lo")
     516             :                 END IF
     517             :              ! -- >        end of abccoflo-part           
     518             :           ENDDO
     519             :        ENDIF
     520             : 
     521             :        ! -->    check whether we have already enough k-vecs
     522       15154 :        enough=.TRUE.
     523       44720 :        DO lo = 1,atoms%nlo(ntyp)
     524       44720 :           IF (nkvec(lo,1).EQ.nkvec(lo,nintsp)) THEN   ! k-vec accepted by both spin channels
     525       29566 :              IF (atoms%invsat(na).EQ.0) THEN
     526        5230 :                 IF ( nkvec(lo,1).LT.(2*atoms%llo(lo,ntyp)+1) ) THEN 
     527        2372 :                    enough=.FALSE.
     528             :                 ENDIF
     529             :              ELSE
     530       24336 :                 IF ( nkvec(lo,1).LT.(2*(2*atoms%llo(lo,ntyp)+1))  ) THEN
     531       12168 :                    enough=.FALSE.
     532             :                 ENDIF
     533             :              ENDIF
     534             :           ELSE
     535           0 :              nkmin = MIN(nkvec(lo,1),nkvec(lo,nintsp)) ! try another k-vec
     536           0 :              nkvec(lo,1) = nkmin ; nkvec(lo,nintsp) = nkmin
     537           0 :              enough=.FALSE.
     538             :           ENDIF
     539             :        ENDDO
     540       15154 :        IF ( enough ) THEN
     541        2078 :           lapw%nkvec(:atoms%nlo(ntyp),na)=nkvec(:atoms%nlo(ntyp),1)
     542             :           RETURN
     543             :        ENDIF
     544             :     ENDDO
     545             : 
     546             :   END SUBROUTINE priv_vec_for_lo
     547           0 : END MODULE m_types_lapw

Generated by: LCOV version 1.13