hsmt_extra.F90 11.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 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 47 48 49 50 51 52 53 54 55 56 57 58 59 60 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 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
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
       kveclo,aa,bb) !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
#ifndef CPP_INVERSION
    USE m_hsmt_hlptomat
#endif
    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 ..

    REAL,    INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntypd,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

#ifdef CPP_INVERSION
    REAL,    INTENT (INOUT) :: aa(:),bb(:)!(matsize)
#else
    COMPLEX, INTENT (INOUT) :: aa(:),bb(:)!(matsize)

#endif
    !     ..
    !     .. 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
#ifndef CPP_INVERSION
                   CALL slomat(&
                        atoms, n,na,lapw,con1,n_size,n_rank, gk,rph,cph,&
                        fj(:,0:,:,fjstart:), gj(:,0:,:,fjstart:),&
                        kvec,isp,usdus,alo1,blo1,clo1,noco, ab_dim,1,1,chi11,chi22,chi21,&
                        iilos,locols,nkvecprevats,bbhlp)
                   CALL hlomat(atoms,isp,isp,n_size,n_rank,&
                        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,&
                        iiloh,locolh,nkvecprevath,tlmplm,aahlp)
#endif
                ELSE
                   jd = 1 ; IF (noco%l_noco) jd = isp
                   DO iintsp = 1,nintsp
                      DO jintsp = 1,nintsp
                         CALL slomat(atoms,n,na,lapw,con1,n_size,n_rank,&
                              gk,rph,cph,fj,gj, kvec,isp,usdus,alo1,blo1,clo1,noco,&
                              ab_dim,iintsp,jintsp,chi11,chi22,chi21,&
                              iilos,locols,nkvecprevats,bb)
                         CALL hlomat(atoms,isp,jd,n_size,n_rank,&
                              n,na,lapw,ar(:,0:,1),br(:,0:,1),ai(:,0:,1),bi(:,0:,1),&
                              el(:,n,isp),alo,blo,clo,usdus, noco,iintsp,jintsp,chi11,chi22,chi21,&
                              iiloh,locolh,nkvecprevath,tlmplm,aa)
                      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 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309
                IF ( noco%l_noco .AND. (.NOT.noco%l_ss) ) THEN
#ifndef CPP_INVERSION
                   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,&
                        nkvecprevatu,iilou,locolu,aahlp)
#endif
                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,&
                              nkvecprevatu,iilou,locolu,aa)
                      ENDDO
                   ENDDO
                ENDIF
             ENDIF

          ENDIF                                ! atoms%invsat(na) = 0 or 1
          !--->    end loop over equivalent atoms
       END DO
#ifndef CPP_INVERSION
       IF ( noco%l_noco .AND. (.NOT. noco%l_ss) ) CALL hsmt_hlptomat(atoms%nlotot,lapw%nv,sub_comm,chi11,chi21,chi22,aahlp,aa,bbhlp,bb)
#endif
       !---> end loop over atom types (ntype)
    ENDDO ntype_loop
    lapw%nmat = lapw%nmat + nkvecprevats

    RETURN
  END SUBROUTINE hsmt_extra
END MODULE m_hsmt_extra