hsmt_extra.F90 12.4 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

7 8 9 10 11 12
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,gwc,noco,l_socfirst, lapw,cell,el, fj,gj,gk,vk,tlmplm,usdus, vs_mmp,oneD,& !in
13
       kveclo,l_real,aa_r,bb_r,aa_c,bb_c) !out/inout
14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
    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
    INTEGER, INTENT (IN) :: gwc
    LOGICAL, INTENT (IN) :: l_socfirst 


    !     ..
    !     .. Array Arguments ..

Daniel Wortmann's avatar
Daniel Wortmann committed
47
    REAL,    INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd)
48 49 50 51 52 53 54 55 56 57
    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


58 59 60
    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
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
    !     ..
    !     .. Local Scalars ..
    REAL con1,termi,termr,th,invsfct


    COMPLEX chi11,chi21,chi22
    INTEGER k,i,spin2,&
         l,ll1,lo,jd,&
         m,n,na,nn,np,&
         iiloh,iilos,nkvecprevath,nkvecprevats,&
         iintsp,jintsp
    INTEGER nc,locolh,locols,nkvecprevatu,iilou,locolu
    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
    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


    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
          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,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)
                      !-->                 apply the rotation that brings this atom into the
                      !-->                 representative (this is the definition of ngopr(na)
                      !-->                 and transform to cartesian coordinates
                      v(:) = vk(k,:,iintsp)
                      gkrot(k,:) = MATMUL(TRANSPOSE(bmrot),v)
                   END DO
                END IF
                DO k = 1,lapw%nv(iintsp)
                   !-->    generate spherical harmonics
                   vmult(:) =  gkrot(k,:)
                   CALL ylm4(atoms%lnonsph(n),vmult,ylm)
                   IF (.NOT.enough) THEN
                      l_lo1=.FALSE. 
                      IF ((lapw%rk(k,iintsp).LT.eps).AND.(.NOT.noco%l_ss)) THEN
                         l_lo1=.TRUE.
                      ELSE
                         l_lo1=.FALSE. 
                      ENDIF
                      CALL abccoflo(&
                           atoms,con1, rph(k,iintsp),cph(k,iintsp),ylm,n,&
                           na,k,lapw%nv(iintsp),l_lo1,alo1,blo1,clo1,&
                           nkvec(1,iintsp),enough, alo(-atoms%llod:,:,:,iintsp),blo(-atoms%llod:,:,:,iintsp),&
                           clo(-atoms%llod:,:,:,iintsp),kvec)
                      !-lo
                   ENDIF
                   !-->              synthesize the complex conjugates of a and b
                   DO l = 0,atoms%lnonsph(n)
                      ll1 = l* (l+1)
                      DO m = -l,l
                         termr = rph(k,iintsp)* REAL(ylm(ll1+m+1)) -&
                              cph(k,iintsp)*AIMAG(ylm(ll1+m+1))
                         termi = rph(k,iintsp)*AIMAG(ylm(ll1+m+1)) +&
                              cph(k,iintsp)* REAL(ylm(ll1+m+1))
                         ar(k,ll1+m,iintsp) = fj(k,l,n,spin2)*termr
                         ai(k,ll1+m,iintsp) = fj(k,l,n,spin2)*termi
                         br(k,ll1+m,iintsp) = gj(k,l,n,spin2)*termr
                         bi(k,ll1+m,iintsp) = gj(k,l,n,spin2)*termi
                      END DO
                   END DO

                END DO
                !--->       end loop over interstitial spin
             ENDDO

             !--->       add the local orbital contribution to the overlap and
             !--->       hamiltonian matrix, if they are used for this atom.

             IF (atoms%nlo(n).GE.1) THEN
                IF (atoms%n_u.GT.0) THEN
                   nkvecprevatu = nkvecprevats
                   iilou = iilos ; locolu = locols
                ENDIF
                IF ( noco%l_noco .AND. (.NOT.noco%l_ss) ) THEN
                   fjstart = 1
                   IF(l_socfirst) fjstart = isp
                   CALL slomat(&
246
                        input,atoms, n,na,lapw,con1,n_size,n_rank, gk,rph,cph,&
247 248
                        fj(:,0:,:,fjstart:), gj(:,0:,:,fjstart:),&
                        kvec,isp,usdus,alo1,blo1,clo1,noco, ab_dim,1,1,chi11,chi22,chi21,&
249
                        iilos,locols,nkvecprevats,.false.,bb_c=bbhlp)
250
                   CALL hlomat(input,atoms,isp,isp,n_size,n_rank,&
251 252
                        n,na,lapw,ar(:,0:,1),br(:,0:,1),ai(:,0:,1),bi(:,0:,1),&
                        el(:,n,isp),alo,blo,clo,usdus, noco,1,1,chi11,chi22,chi21,&
253
                        iiloh,locolh,nkvecprevath,tlmplm,.false.,aa_c=aahlp)
254 255 256 257
                ELSE
                   jd = 1 ; IF (noco%l_noco) jd = isp
                   DO iintsp = 1,nintsp
                      DO jintsp = 1,nintsp
258
                         CALL slomat(input,atoms,n,na,lapw,con1,n_size,n_rank,&
259 260
                              gk,rph,cph,fj,gj, kvec,isp,usdus,alo1,blo1,clo1,noco,&
                              ab_dim,iintsp,jintsp,chi11,chi22,chi21,&
261
                              iilos,locols,nkvecprevats,l_real,bb_r,bb_c)
262
                         CALL hlomat(input,atoms,isp,jd,n_size,n_rank,&
263
                              n,na,lapw,ar(:,0:,jintsp),br(:,0:,jintsp),ai(:,0:,jintsp),bi(:,0:,jintsp),&
264
                              el(:,n,isp),alo,blo,clo,usdus, noco,iintsp,jintsp,chi11,chi22,chi21,&
265
                              iiloh,locolh,nkvecprevath,tlmplm,l_real,aa_r,aa_c)
266 267 268 269 270
                      ENDDO
                   ENDDO
                ENDIF
             END IF

271
             IF (atoms%n_u>0.and.atoms%lda_u(n)%l.GE.0.AND.gwc.EQ.1) THEN
272 273 274 275 276 277 278
                IF ( noco%l_noco .AND. (.NOT.noco%l_ss) ) THEN
                   CALL u_ham(&
                        atoms,input,lapw,isp,n,invsfct,&
                        ar,ai,br,bi,vs_mmp,lmaxb,&
                        alo,blo,clo,&
                        n_size,n_rank,isp,usdus,noco,&
                        1,1,chi11,chi22,chi21,&
279
                        nkvecprevatu,iilou,locolu,.false.,aa_c=aahlp)
280 281 282 283 284 285 286 287 288
                ELSE
                   DO iintsp = 1,nintsp
                      DO jintsp = 1,iintsp
                         CALL u_ham(&
                              atoms,input,lapw,isp,n,invsfct,&
                              ar,ai,br,bi,vs_mmp,lmaxb,&
                              alo,blo,clo,&
                              n_size,n_rank,isp,usdus,noco,&
                              iintsp,jintsp,chi11,chi22,chi21,&
289
                              nkvecprevatu,iilou,locolu,l_real,aa_r,aa_c)
290 291 292 293 294 295 296 297
                      ENDDO
                   ENDDO
                ENDIF
             ENDIF

          ENDIF                                ! atoms%invsat(na) = 0 or 1
          !--->    end loop over equivalent atoms
       END DO
298
       IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) CALL hsmt_hlptomat(atoms%nlotot,lapw%nv,sub_comm,chi11,chi21,chi22,aahlp,aa_c,bbhlp,bb_c)
299 300 301 302 303 304 305
       !---> end loop over atom types (ntype)
    ENDDO ntype_loop
    lapw%nmat = lapw%nmat + nkvecprevats

    RETURN
  END SUBROUTINE hsmt_extra
END MODULE m_hsmt_extra