Deleted some old versions of subroutines in eigen

eigen/hsint.F90

-!--------------------------------------------------------------------------------
-! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
-! This file is part of FLEUR and available as free software under the conditions
-! of the MIT license as expressed in the LICENSE file in more detail.
-!--------------------------------------------------------------------------------
-
-MODULE m_hsint
-CONTAINS
-  SUBROUTINE hsint(input,noco,jij,stars, vpw,lapw,jspin,&
-       n_size,n_rank,bkpt,cell,atoms,l_real,hamOvlp)
-    !*********************************************************************
-    !     initializes and sets up the hamiltonian and overlap matrices
-    !     for the interstitial. only the lower triangle of the hermitian
-    !     matrices are stored in compact real mode such that if h(i,j),
-    !     i.ge.j, is hermitian and a is real, then
-    !       a(i,j)=real( h(i,j) )  and  a(j,i)=aimag( h(i,j) )
-    !                    m. weinert  1986
-    !
-    ! For the eigenvector parallelization each pe calculates an equal share
-    ! of columns labeled nc. Then the starting element of a columns nc is
-    !
-    !     ii = (nc-1)*( n_rank - n_size + 1 ) + n_size*(nc-1)*nc/2
-    !
-    ! and, if a non-collinear matrix has to be set up, the starting column
-    ! for the second spin-direction is
-    !
-    !     nc = int( 1. + (nv - n_rank - 1)/n_size ) + 1 .
-    !
-    ! For this direction, the outer loop starts at
-    !
-    ! istart = n_rank + (nc - 1)*n_size - nv .                        gb99
-    !
-    ! for a lo-calculation nv has to be replaced by nv+nlotot         gb01
-    !
-    !*********************************************************************
-    USE m_types
-    IMPLICIT NONE
-    TYPE(t_input),INTENT(IN)      :: input
-    TYPE(t_noco),INTENT(IN)       :: noco
-    TYPE(t_jij),INTENT(IN)        :: jij
-    TYPE(t_stars),INTENT(IN)      :: stars
-    TYPE(t_cell),INTENT(IN)       :: cell
-    TYPE(t_atoms),INTENT(IN)      :: atoms
-    TYPE(t_lapw),INTENT(INOUT)    :: lapw
-    TYPE(t_hamOvlp),INTENT(INOUT) :: hamOvlp
-    !     ..
-    !     .. Scalar Arguments ..
-    INTEGER, INTENT (IN) :: n_size,n_rank,jspin
-    !     ..
-    !     .. Array Arguments ..
-    COMPLEX, INTENT (INOUT) :: vpw(stars%ng3)
-    REAL,    INTENT (IN)    :: bkpt(3) 
-    LOGICAL,INTENT(IN)      :: l_real
-    !     ..
-    !     .. Local Scalars ..
-    COMPLEX th,ts,phase
-    REAL b1(3),b2(3),r2
-    INTEGER i,i1,i2,i3,ii,in,j,ig3,ispin,l,iloc
-    INTEGER istart,nc
-
-    COMPLEX ust1,vp1
-    COMPLEX, ALLOCATABLE :: vpw1(:)  ! for J constants
-    !     ..
-    ! ..
-    !$OMP PARALLEL 
-    if (l_real) THEN
-       !$OMP DO 
-       do i = 1, size(hamOvlp%a_r)
-           hamOvlp%a_r(i)=0.0
-       end do
-       !OMP END DO
-       !$OMP DO 
-       do i = 1, size(hamOvlp%b_r)
-           hamOvlp%b_r(i)=0.0
-       end do
-       !OMP END DO
-    ELSE
-       !$OMP DO 
-       do i = 1, size(hamOvlp%a_c)
-           hamOvlp%a_c(i)=0.0
-       end do
-       !$OMP END DO
-       !$OMP DO 
-       do i = 1, size(hamOvlp%b_c)
-           hamOvlp%b_c(i)=0.0
-       end do
-       !$OMP END DO
-    ENDIF
-    !$OMP END PARALLEL 
-    ust1 = stars%ustep(1)
-    ispin = jspin
-    lapw%nmat = lapw%nv(ispin)
-
-    !---> pk non-collinear
-    IF (noco%l_noco) THEN
-       !---> determine spin-up spin-up part of Hamiltonian- and overlapp-matrix
-       !---> reload V_11
-       READ (25) (vpw(ig3),ig3=1,stars%ng3)
-
-       !--- J const
-       IF( jij%l_J) THEN
-          ALLOCATE ( vpw1(stars%ng3) )
-          READ (25) (vpw1(ig3),ig3=1,stars%ng3)
-       ENDIF
-       !--- J const
-
-       lapw%nmat = lapw%nv(1) + lapw%nv(2)
-       ispin = 1
-
-       !--- J const
-       IF (jij%l_J) THEN
-          DO i = 1,stars%ng3
-             vpw(i) = (vpw(i) + vpw1(i))/2.
-          END DO
-       ENDIF
-       !--- J const
-
-       vp1 = REAL(vpw(1))
-    ENDIF
-    !---> pk non-collinear
-
-    vp1 = vpw(1)
-    !---> loop over (k+g')
-    ii = 0
-    !$OMP PARALLEL DO SCHEDULE(dynamic) DEFAULT(none) &
-    !$OMP SHARED(n_rank,n_size,lapw,ispin,stars,input,bkpt,cell,vpw,ust1,vp1) &
-    !$OMP SHARED(l_real,hamOvlp)&
-    !$OMP PRIVATE(i,j,iloc,i1,i2,i3,in,phase,b1,b2,r2,th,ts)&
-    !$OMP FIRSTPRIVATE(ii)
-    DO  i = n_rank+1, lapw%nv(ispin), n_size
-       !--->    loop over (k+g)
-       DO  j = 1,i - 1
-          ii = 0
-          DO iloc = n_rank+1,i-n_size,n_size
-             ii = ii + iloc
-          ENDDO
-          ii = ii + j
-          !-->     determine index and phase factor
-          i1 = lapw%k1(i,ispin) - lapw%k1(j,ispin)
-          i2 = lapw%k2(i,ispin) - lapw%k2(j,ispin)
-          i3 = lapw%k3(i,ispin) - lapw%k3(j,ispin)
-          in = stars%ig(i1,i2,i3)
-          IF (in.EQ.0) CYCLE
-          phase = stars%rgphs(i1,i2,i3)
-          !+APW_LO
-          IF (input%l_useapw) THEN
-             b1(1) = bkpt(1)+lapw%k1(i,ispin) ; b2(1) = bkpt(1)+lapw%k1(j,ispin)
-             b1(2) = bkpt(2)+lapw%k2(i,ispin) ; b2(2) = bkpt(2)+lapw%k2(j,ispin)
-             b1(3) = bkpt(3)+lapw%k3(i,ispin) ; b2(3) = bkpt(3)+lapw%k3(j,ispin)
-             r2 = DOT_PRODUCT(MATMUL(b2,cell%bbmat),b1)   
-
-             th = phase*(0.5*r2*stars%ustep(in)+vpw(in))
-          ELSE
-             th = phase* (0.25* (lapw%rk(i,ispin)**2+lapw%rk(j,ispin)**2)*stars%ustep(in) + vpw(in))
-          ENDIF
-          !-APW_LO
-          !--->    determine matrix element and store
-          ts = phase*stars%ustep(in)
-          if (l_real) THEN
-          hamOvlp%a_r(ii) = REAL(th)
-          hamOvlp%b_r(ii) = REAL(ts)
-          else
-          hamOvlp%a_c(ii) = th
-          hamOvlp%b_c(ii) = ts
-          endif
-       ENDDO
-       !--->    diagonal term (g-g'=0 always first star)
-       ii = ii + 1
-       if (l_real) THEN
-       hamOvlp%a_r(ii) = 0.5*lapw%rk(i,ispin)*lapw%rk(i,ispin)*REAL(ust1) + REAL(vp1)
-       hamOvlp%b_r(ii) = REAL(ust1)
-       else
-       hamOvlp%a_c(ii) = 0.5*lapw%rk(i,ispin)*lapw%rk(i,ispin)*ust1 + vp1
-       hamOvlp%b_c(ii) = ust1
-       endif
-    ENDDO
-    !$OMP END PARALLEL DO
-
-    !---> pk non-collinear
-    IF (noco%l_noco) THEN
-       !+gb99
-       nc = INT( 1. + (lapw%nv(1)+atoms%nlotot - n_rank - 1)/n_size )
-       istart = n_rank + nc*n_size - (lapw%nv(1)+atoms%nlotot)
-       !      ii = (nv(1)+nlotot+1)*(nv(1)+nlotot+2)/2 - 1
-       ii = nc*(n_rank-n_size+1) + n_size*(nc+1)*nc/2 + lapw%nv(1)+atoms%nlotot
-       !-gb99
-       ispin = 2
-       !---> determine spin-down spin-down part of Hamiltonian- and ovlp-matrix
-       !---> reload V_22
-
-       !--- J constants 
-       IF(.NOT.jij%l_J) THEN
-          READ (25) (vpw(ig3),ig3=1,stars%ng3)
-          vp1 = REAL(vpw(1))
-       ENDIF
-       !--- J constants
-
-       !---> loop over (k+g')
-       DO i = istart+1, lapw%nv(ispin), n_size
-          nc = nc + 1
-          !--->    loop over (k+g)
-          DO j = 1,i - 1
-             !-gb99      ii = (nv(1)+i-1)*(nv(1)+i)/2 + nv(1) + j
-             ii = (nc-1)*( n_rank - n_size + 1 ) + n_size*(nc-1)*nc/2 + lapw%nv(1)+atoms%nlotot + j
-             !--->       determine index and phase factor
-             i1 = lapw%k1(i,ispin) - lapw%k1(j,ispin)
-             i2 = lapw%k2(i,ispin) - lapw%k2(j,ispin)
-             i3 = lapw%k3(i,ispin) - lapw%k3(j,ispin)
-             in = stars%ig(i1,i2,i3)
-             IF (in.EQ.0) THEN
-                WRITE (*,*) 'HSINT: G-G'' not in star i,j= ',i,j
-             ELSE
-                phase = stars%rgphs(i1,i2,i3)
-                !+APW_LO
-                IF (input%l_useapw) THEN
-                   b1(1) = bkpt(1)+lapw%k1(i,ispin) ; b2(1) = bkpt(1)+lapw%k1(j,ispin)
-                   b1(2) = bkpt(2)+lapw%k2(i,ispin) ; b2(2) = bkpt(2)+lapw%k2(j,ispin)
-                   b1(3) = bkpt(3)+lapw%k3(i,ispin) ; b2(3) = bkpt(3)+lapw%k3(j,ispin)
-                   r2 = DOT_PRODUCT(MATMUL(b2,cell%bbmat),b1)   
-
-                   th = phase*( 0.5*r2*stars%ustep(in) + vpw(in) )
-                ELSE
-                   th = phase* (0.25* (lapw%rk(i,ispin)**2+lapw%rk(j,ispin)**2)*stars%ustep(in) + vpw(in))
-                ENDIF
-                !-APW_LO
-                ts = phase*stars%ustep(in)
-                hamOvlp%a_c(ii) = th
-                hamOvlp%b_c(ii) = ts
-             ENDIF
-          ENDDO
-          !--->    diagonal term (g-g'=0 always first star)
-          !-gb99   ii = (nv(1)+i)*(nv(1)+i+1)/2
-          ii = ii + 1
-          hamOvlp%a_c(ii) = 0.5*lapw%rk(i,ispin)*lapw%rk(i,ispin)*ust1 + vp1
-          hamOvlp%b_c(ii) = ust1
-       ENDDO
-
-       !---> determine spin-down spin-up part of Hamiltonian- and ovlp-matrix
-       !---> reload real part of V_21
-       READ (25) (vpw(ig3),ig3=1,stars%ng3)
-       nc = INT( 1. + (lapw%nv(1)+atoms%nlotot - n_rank - 1)/n_size )
-       !
-       !---> loop over (k+g')
-       !
-       DO i = istart+1, lapw%nv(2), n_size
-          nc = nc + 1
-          !--->    loop over (k+g)
-          DO j = 1,lapw%nv(1)
-             !-gb99      ii = (nv(1)+i-1)*(nv(1)+i)/2 + j
-             ii = (nc-1)*( n_rank - n_size + 1 ) + n_size*(nc-1)*nc/2 + j
-             !--->       determine index and phase factor
-             i1 = lapw%k1(i,2) - lapw%k1(j,1)
-             i2 = lapw%k2(i,2) - lapw%k2(j,1)
-             i3 = lapw%k3(i,2) - lapw%k3(j,1)
-             in = stars%ig(i1,i2,i3)
-             IF (in.EQ.0) THEN
-                WRITE (*,*) 'HSINT: G-G'' not in star i,j= ',i,j
-             ELSE
-                hamOvlp%a_c(ii) = stars%rgphs(i1,i2,i3)*vpw(in) 
-                !--- J constants 
-                IF(jij%l_J) THEN
-                   hamOvlp%a_c(ii) = 0
-                ENDIF
-                !--- J constants
-
-             ENDIF
-          ENDDO
-       ENDDO
-       !---> pk non-collinear
-    ENDIF
-
-    IF (jij%l_J) DEALLOCATE (vpw1)
-
-    RETURN
-  END SUBROUTINE hsint
-END MODULE m_hsint

eigen/hsmt_blas.F90

-!--------------------------------------------------------------------------------
-! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
-! This file is part of FLEUR and available as free software under the conditions
-! of the MIT license as expressed in the LICENSE file in more detail.
-!--------------------------------------------------------------------------------
-MODULE m_hsmt_blas
-  use m_juDFT
-  implicit none
-CONTAINS
-  SUBROUTINE hsmt_blas(sym,atoms,isp,noco,cell,lapw,td,ud,gk,vk,fj,gj,smat,hmat)
-!Calculate overlap matrix
-    USE m_hsmt_ab
-    USE m_constants, ONLY : fpi_const,tpi_const
-    USE m_types
-    USE m_ylm
-    IMPLICIT NONE
-    TYPE(t_sym),INTENT(IN)      :: sym
-    TYPE(t_noco),INTENT(IN)     :: noco
-    TYPE(t_cell),INTENT(IN)     :: cell
-    TYPE(t_atoms),INTENT(IN)    :: atoms
-    TYPE(t_lapw),INTENT(IN)     :: lapw
-    TYPE(t_tlmplm),INTENT(IN)   :: td
-    TYPE(t_usdus),INTENT(IN)    :: ud
-    !     ..
-    !     .. Scalar Arguments ..
-    INTEGER, INTENT (IN) :: isp
-    !     ..
-    !     .. Array Arguments ..
-    REAL,INTENT(IN) :: gk(:,:,:),vk(:,:,:)
-    REAL,INTENT(IN) :: fj(:,0:,:,:),gj(:,0:,:,:)
-    TYPE(t_lapwmat),INTENT(INOUT)::smat,hmat
-
-    
-    INTEGER:: n,nn,na,aboffset,l,ll,m
-    COMPLEX,ALLOCATABLE:: ab(:,:),tmpdata(:,:),tmp_s(:,:),tmp_h(:,:),ab1(:,:)
-
-
-    ALLOCATE(ab(lapw%nv(isp),2*atoms%lmaxd*(atoms%lmaxd+2)+2),ab1(lapw%nv(isp),2*atoms%lmaxd*(atoms%lmaxd+2)+2))
-    
-    ALLOCATE(tmp_s(smat%matsize1,smat%matsize2),tmp_h(smat%matsize1,smat%matsize2))
-    tmp_s=0.0;tmp_h=0.0;ab=0.0;ab1=0.0
-
-    ntyploop: DO n=1,atoms%ntype
-       DO nn = 1,atoms%neq(n)
-          na = SUM(atoms%neq(:n-1))+nn
-          IF ((atoms%invsat(na)==0) .OR. (atoms%invsat(na)==1)) THEN
-             
-             !--->   Calculate Overlapp matrix
-             CALL timestart("ab-coefficients")
-             CALL hsmt_ab(sym,atoms,isp,n,na,cell,lapw,gk,vk,fj,gj,ab,aboffset)
-             CALL timestop("ab-coefficients")
-             CALL timestart("Overlapp")
-             CALL ZHERK("U","N",lapw%nv(isp),aboffset,1.,ab,SIZE(ab,1),1.0,tmp_s,SIZE(tmp_s,1))
-             DO l=0,atoms%lmax(n)
-                ll=l*(l+1)
-                DO m=-l,l
-                   ab1(:,1+ll+m)=SQRT(ud%ddn(l,n,isp))*ab(:,aboffset+1+ll+m)
-                ENDDO
-             ENDDO
-             CALL ZHERK("U","N",lapw%nv(isp),aboffset,1.,ab1,SIZE(ab,1),1.0,tmp_s,SIZE(tmp_s,1))
-             CALL timestop("Overlapp")
-             CALL timestart("Hamiltonian")
-         
-             !Calculate Hamiltonian
-             CALL zgemm("N","N",SIZE(ab,1),2*aboffset,2*aboffset,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(:,:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab1,SIZE(ab,1))
-!             CALL zgemm("N","C",lapw%nv(isp),lapw%nv(isp),2*aboffset,CMPLX(1.0,0.0),ab,SIZE(ab,1),ab1,SIZE(ab,1),CMPLX(1.0,0),tmp_h,SIZE(tmp_h,1))
-             CALL ZHERK("U","N",lapw%nv(isp),2*aboffset,1.,ab1,SIZE(ab,1),1.0,tmp_h,SIZE(tmp_h,1))
-             CALL timestop("Hamiltonian")
-          ENDIF
-       END DO
-    END DO ntyploop
-    !Copy tmp array back
-    IF (smat%l_real) THEN
-       smat%data_r=smat%data_r+tmp_s
-       hmat%data_r=hmat%data_r+tmp_h-td%e_shift*tmp_s
-    ELSE
-       smat%data_c=smat%data_c+tmp_s
-       hmat%data_c=hmat%data_c+tmp_h-td%e_shift*tmp_s
-    ENDIF
-  END SUBROUTINE hsmt_blas
-
-#if 1==2
-    !this version uses zherk for Hamiltonian
-    ntyploop: DO n=1,atoms%ntype
-       DO nn = 1,atoms%neq(n)
-          na = SUM(atoms%neq(:n-1))+nn
-          IF ((atoms%invsat(na)==0) .OR. (atoms%invsat(na)==1)) THEN
-             
-             !--->   Calculate Overlapp matrix
-             CALL timestart("ab-coefficients")
-             CALL hsmt_ab(sym,atoms,isp,n,na,cell,lapw,gk,vk,fj,gj,ab,aboffset)
-             CALL timestop("ab-coefficients")
-             CALL timestart("Overlapp")
-             CALL ZHERK("U","N",lapw%nv(isp),aboffset,1.,ab,SIZE(ab,1),1.0,tmp_s,SIZE(tmp_s,1))
-             DO l=0,atoms%lmax(n)
-                ll=l*(l+1)
-                DO m=-l,l
-                   ab1(:,1+ll+m)=SQRT(ud%ddn(l,n,isp))*ab(:,aboffset+1+ll+m)
-                ENDDO
-             ENDDO
-             CALL ZHERK("U","N",lapw%nv(isp),aboffset,1.,ab1,SIZE(ab,1),1.0,tmp_s,SIZE(tmp_s,1))
-             CALL timestop("Overlapp")
-             CALL timestart("Hamiltonian")
-         
-             !Calculate Hamiltonian
-             CALL zgemm("N","N",SIZE(ab,1),2*aboffset,2*aboffset,CMPLX(1.0,0.0),ab,SIZE(ab,1),td%h_loc(:,:,n,isp),SIZE(td%h_loc,1),CMPLX(0.,0.),ab1,SIZE(ab,1))
-             CALL ZHERK("U","N",lapw%nv(isp),2*aboffset,1.,ab1,SIZE(ab,1),1.0,tmp_h,SIZE(tmp_h,1))
-             CALL timestop("Hamiltonian")
-          ENDIF
-       END DO
-    END DO ntyploop
-    !Copy tmp array back
-    IF (smat%l_real) THEN
-       smat%data_r=smat%data_r+tmp_s
-       hmat%data_r=hmat%data_r+tmp_h-td%e_shift*tmp_s
-    ELSE
-       smat%data_c=smat%data_c+tmp_s
-       hmat%data_c=hmat%data_c+tmp_h-td%e_shift*tmp_s
-    ENDIF
-#endif
-  
-END MODULE m_hsmt_blas

eigen/hsmt_extra.F90

-!--------------------------------------------------------------------------------
-! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
-! This file is part of FLEUR and available as free software under the conditions
-! of the MIT license as expressed in the LICENSE file in more detail.
-!--------------------------------------------------------------------------------
-
-MODULE m_hsmt_extra
-  USE m_juDFT
-  IMPLICIT NONE
-CONTAINS
-  SUBROUTINE hsmt_extra(DIMENSION,atoms,sym,isp,n_size,n_rank,input,nintsp,sub_comm,&
-       hlpmsize,lmaxb,noco,l_socfirst, lapw,cell,el, fj,gj,gk,vk,tlmplm,usdus, vs_mmp,oneD,& !in
-       kveclo,l_real,aa_r,bb_r,aa_c,bb_c) !out/inout
-    USE m_constants, ONLY : tpi_const,fpi_const
-    USE m_uham
-    USE m_ylm
-    USE m_abccoflo
-    USE m_hlomat
-    USE m_slomat
-    USE m_vecforlo
-    USE m_setabc1lo
-    USE m_hsmt_spinor
-    USE m_hsmt_hlptomat
-    USE m_types
-    IMPLICIT NONE
-    TYPE(t_dimension),INTENT(IN):: DIMENSION
-    TYPE(t_oneD),INTENT(IN)     :: oneD
-    TYPE(t_input),INTENT(IN)    :: input
-    TYPE(t_noco),INTENT(IN)     :: noco
-    TYPE(t_sym),INTENT(IN)      :: sym
-    TYPE(t_cell),INTENT(IN)     :: cell
-    TYPE(t_atoms),INTENT(IN)    :: atoms
-    TYPE(t_lapw),INTENT(INOUT)  :: lapw !lapw%nmat is updated with lo-basisfcts
-    !     ..
-    !     .. Scalar Arguments ..
-    INTEGER, INTENT (IN) :: isp
-    INTEGER, INTENT (IN) :: n_size,n_rank ,nintsp,sub_comm
-    INTEGER, INTENT (IN) :: hlpmsize
-    INTEGER, INTENT (IN) :: lmaxb
-    LOGICAL, INTENT (IN) :: l_socfirst 
-
-
-    !     ..
-    !     .. Array Arguments ..
-
-    REAL,    INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd)
-    COMPLEX,INTENT(IN)   :: vs_mmp(-lmaxb:lmaxb,-lmaxb:lmaxb,atoms%n_u,input%jspins)
-    REAL,INTENT(IN),ALLOCATABLE :: fj(:,:,:,:),gj(:,:,:,:),gk(:,:,:),vk(:,:,:)
-
-
-    TYPE(t_usdus),INTENT(IN)   :: usdus
-
-    INTEGER, INTENT (OUT)    :: kveclo(atoms%nlotot)
-    TYPE(t_tlmplm),INTENT(INOUT)::tlmplm
-
-
-    REAL,    OPTIONAL,ALLOCATABLE,INTENT (INOUT) :: aa_r(:),bb_r(:)!(matsize)
-    COMPLEX, OPTIONAL,ALLOCATABLE,INTENT (INOUT) :: aa_c(:),bb_c(:)!(matsize)
-    LOGICAL,INTENT(IN)      :: l_real
-    !     ..
-    !     .. Local Scalars ..
-    REAL con1,termi,termr,th,invsfct
-
-
-    COMPLEX chi11,chi21,chi22
-    INTEGER k,i,spin2,l,ll1,lo,jd
-<<<<<<< HEAD
-    INTEGER m,n,na,nn,np,i_u,i_u_save
-=======
-    INTEGER m,n,na,nn,np,i_u
->>>>>>> hsmt_simple
-    INTEGER iiloh,iilos,nkvecprevath,nkvecprevats,iintsp,jintsp
-    INTEGER nc,locolh,locols,nkvecprevatu,iilou,locolu
-    INTEGER nkvecprevatuTemp,iilouTemp,locoluTemp
-    INTEGER ab_dim,nkvec_sv,fjstart
-    LOGICAL enough,l_lo1
-    !     ..
-    !     .. Local Arrays ..
-    INTEGER kvec(2*(2*atoms%llod+1),atoms%nlod )
-    REAL alo1(atoms%nlod),blo1(atoms%nlod),clo1(atoms%nlod)
-    REAL bmrot(3,3),gkrot(DIMENSION%nvd,3),vmult(3),v(3)
-    REAL qssbti(3)
-    REAL, ALLOCATABLE :: ar(:,:,:),ai(:,:,:),br(:,:,:),bi(:,:,:)
-    REAL, ALLOCATABLE :: rph(:,:),cph(:,:)
-    COMPLEX, ALLOCATABLE :: alo(:,:,:,:),blo(:,:,:,:),clo(:,:,:,:)
-    INTEGER, ALLOCATABLE :: nkvec(:,:)
-    COMPLEX,ALLOCATABLE :: aahlp(:),bbhlp(:)
-
-    COMPLEX ylm( (atoms%lmaxd+1)**2 )
-
-    COMPLEX chi(2,2)
-
-
-    REAL, PARAMETER :: eps = 1.0E-30
-
-    con1=fpi_const/SQRT(cell%omtil)
-    ab_dim=1
-    IF ( noco%l_noco .AND. (.NOT. noco%l_ss) )  ALLOCATE ( aahlp(hlpmsize),bbhlp(hlpmsize) )
-    IF (noco%l_ss) ab_dim=2
-    ALLOCATE(ar(DIMENSION%nvd,0:DIMENSION%lmd,ab_dim),ai(DIMENSION%nvd,0:DIMENSION%lmd,ab_dim))
-    ALLOCATE(br(DIMENSION%nvd,0:DIMENSION%lmd,ab_dim),bi(DIMENSION%nvd,0:DIMENSION%lmd,ab_dim))
-    ALLOCATE(alo(-atoms%llod:atoms%llod,2*(2*atoms%llod+1),atoms%nlod,ab_dim))
-    ALLOCATE(blo(-atoms%llod:atoms%llod,2*(2*atoms%llod+1),atoms%nlod,ab_dim))
-    ALLOCATE(clo(-atoms%llod:atoms%llod,2*(2*atoms%llod+1),atoms%nlod,ab_dim))
-    ALLOCATE(rph(DIMENSION%nvd,ab_dim),cph(DIMENSION%nvd,ab_dim))
-    ALLOCATE(nkvec(atoms%nlod,ab_dim))
-    na = 0
-    nkvecprevats = 0
-    nkvecprevath = 0
-    nkvecprevatu = 0
-    nkvec_sv = 0
-    !Determine index of first LO
-    locols = lapw%nv(1)
-    locolh = lapw%nv(1)
-#ifdef CPP_MPI
-    nc = 0
-    k = 0
-    DO  i=1+n_rank, lapw%nv(1), n_size
-       nc = nc + 1
-       k = k + n_size*(nc-1) + n_rank + 1
-    ENDDO
-    iilos = k
-    iiloh = k
-#else
-    iilos = lapw%nv(1)* (lapw%nv(1)+1)/2 
-    iiloh = lapw%nv(1)* (lapw%nv(1)+1)/2
-#endif
-
-    iilou = iilos
-    locolu = locols
-
-    i_u = 1
-    ntype_loop: DO n=1,atoms%ntype
-
-       IF (noco%l_noco) THEN
-          IF (.NOT.noco%l_ss) aahlp=CMPLX(0.,0.)
-          IF (.NOT.noco%l_ss) bbhlp=CMPLX(0.,0.)
-          CALL hsmt_spinor(isp,n,noco,input,chi,chi11,chi21,chi22)
-       ENDIF
-       DO nn = 1,atoms%neq(n)
-          na = na + 1
-          i_u_save = i_u
-          IF (atoms%lnonsph(n).LT.0) CYCLE ntype_loop
-          IF ((atoms%invsat(na).EQ.0) .OR. (atoms%invsat(na).EQ.1)) THEN
-             IF (atoms%invsat(na).EQ.0) invsfct = 1
-             IF (atoms%invsat(na).EQ.1) invsfct = 2
-             np = sym%invtab(atoms%ngopr(na))
-             IF (oneD%odi%d1) THEN
-                np = oneD%ods%ngopr(na)
-             END IF
-             CALL vec_for_lo(atoms,nintsp,sym,l_real,na, n,np,noco, lapw,cell, gk,vk, nkvec,kvec)
-             DO lo = 1,atoms%nlo(n)
-                kveclo(nkvec_sv+1:nkvec_sv+nkvec(lo,1)) = kvec(1:nkvec(lo,1),lo)
-                nkvec_sv = nkvec_sv+nkvec(lo,1)
-                nkvec(lo,:) = 0
-             ENDDO
-
-             !--->       loop over interstitial spins
-             DO iintsp = 1,nintsp
-                IF (noco%l_constr.OR.l_socfirst) THEN
-                   spin2=isp
-                ELSE
-                   spin2=iintsp
-                ENDIF
-                IF (iintsp==1) THEN
-                   qssbti(:) = - noco%qss(:)/2
-                ELSE
-                   qssbti(:) = + noco%qss(:)/2
-                ENDIF
-                !--->          set up phase factors
-                DO k = 1,lapw%nv(iintsp)
-                   th= DOT_PRODUCT((/lapw%k1(k,iintsp),lapw%k2(k,iintsp),lapw%k3(k,iintsp)/)+qssbti,atoms%taual(:,na))
-                   rph(k,iintsp) = COS(tpi_const*th)
-                   cph(k,iintsp) = -SIN(tpi_const*th)
-                END DO
-
-                !--->          set up the a,b and c  coefficients
-                !--->          for the local orbitals, if necessary.
-                IF (atoms%nlo(n).GE.1) THEN
-                   enough = .FALSE.
-                   CALL setabc1lo(atoms, n,usdus,isp, alo1,blo1,clo1) 
-                ELSE
-                   enough = .TRUE.
-                END IF
-
-                !--->          set up the a and b coefficients
-
-                IF (np==1) THEN
-                   gkrot( 1:lapw%nv(iintsp),:) = gk( 1:lapw%nv(iintsp),:,iintsp)
-                ELSE 
-                   IF (oneD%odi%d1) THEN
-                      bmrot=MATMUL(oneD%ods%mrot(:,:,np),cell%bmat)
-                   ELSE
-                      bmrot=MATMUL(1.*sym%mrot(:,:,np),cell%bmat)
-                   END IF
-                   DO k = 1,lapw%nv(iintsp)