sorad.f90 11.8 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
MODULE m_sorad
  USE m_juDFT
  !*********************************************************************
  !     1. generates radial spin-orbit matrix elements
  !     based on m.weinert's radsra and radsrd subroutines
  !*********************************************************************
CONTAINS
  SUBROUTINE sorad(atoms,input,ntyp,vr,enpara,spav,&
       rsopp,rsopdpd,rsoppd,rsopdp,usdus,&
       rsoplop,rsoplopd,rsopdplo,rsopplo,rsoploplop)

    USE m_constants, ONLY : c_light
    USE m_intgr,     ONLY : intgr0
    USE m_sointg
    USE m_radsra
    USE m_radsrd
    USE m_radsrdn
    USE m_types
    IMPLICIT NONE
    TYPE(t_enpara),INTENT(IN)   :: enpara
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_atoms),INTENT(IN)   :: atoms
23
    TYPE(t_usdus),INTENT(INOUT)   :: usdus
24 25 26 27 28 29 30
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: ntyp
    LOGICAL, INTENT (IN) :: spav ! if T, spin-averaged pot is used
    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN) :: vr(:,:)!(atoms%jmtd,dimension%jspd),
31 32 33 34 35 36 37 38 39
    REAL,    INTENT (INOUT) :: rsopp  (atoms%ntypd,atoms%lmaxd,2,2)
    REAL,    INTENT (INOUT) :: rsoppd (atoms%ntypd,atoms%lmaxd,2,2)
    REAL,    INTENT (INOUT) :: rsopdp (atoms%ntypd,atoms%lmaxd,2,2)
    REAL,    INTENT (INOUT) :: rsopdpd(atoms%ntypd,atoms%lmaxd,2,2)
    REAL,    INTENT (INOUT) :: rsoplop (atoms%ntypd,atoms%nlod,2,2)
    REAL,    INTENT (INOUT) :: rsoplopd(atoms%ntypd,atoms%nlod,2,2)
    REAL,    INTENT (INOUT) :: rsopdplo(atoms%ntypd,atoms%nlod,2,2)
    REAL,    INTENT (INOUT) :: rsopplo (atoms%ntypd,atoms%nlod,2,2)
    REAL,    INTENT (INOUT) :: rsoploplop(atoms%ntypd,atoms%nlod,atoms%nlod,2,2)
40 41 42 43 44 45 46 47 48 49 50
    !     ..
    !     .. Local Scalars ..
    REAL ddn1,e ,ulops,dulops,duds1
    INTEGER i,j,ir,jspin,l,noded,nodeu,ilo,ilop
    !     ..
    !     .. Local Arrays ..
    REAL, ALLOCATABLE :: p(:,:),pd(:,:),q(:,:),qd(:,:),plo(:,:)
    REAL, ALLOCATABLE :: plop(:,:),glo(:,:),fint(:),pqlo(:,:)
    REAL, ALLOCATABLE :: filo(:,:)
    REAL, ALLOCATABLE :: v0(:),vso(:,:),qlo(:,:)
    !     ..
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
    IF (atoms%jri(ntyp)>atoms%jmtd)  CALL juDFT_error("atoms%jri(ntyp).GT.atoms%jmtd",calledby ="sorad")
    ALLOCATE ( p(atoms%jmtd,2),pd(atoms%jmtd,2),q(atoms%jmtd,2),plo(atoms%jmtd,2),fint(atoms%jmtd),&
         &   qlo(atoms%jmtd,2),plop(atoms%jmtd,2),qd(atoms%jmtd,2),v0(atoms%jmtd),vso(atoms%jmtd,2) )
    !
    DO l = 0,atoms%lmax(ntyp) 

       DO jspin = 1,input%jspins
          !
          !--->    calculate normalized function at e: p and q 
          !
          e = enpara%el0(l,ntyp,jspin)
          CALL radsra(&
               e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
               usdus%us(l,ntyp,jspin),usdus%dus(l,ntyp,jspin),&
               nodeu,p(:,jspin),q(:,jspin))
          !                     
          !--->    calculate orthogonal energy derivative at e : pd and qd
          !
          CALL radsrd(&
               e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
               usdus%uds(l,ntyp,jspin),usdus%duds(l,ntyp,jspin),&
               usdus%ddn(l,ntyp,jspin),noded,pd(:,jspin),qd(:,jspin),&
               p(:,jspin),q(:,jspin),usdus%dus(l,ntyp,jspin))

       END DO     ! end of spin loop
       !
       !---> in case of jspins=1
       !

       IF (input%jspins.EQ.1) THEN
          DO i = 1,atoms%jri(ntyp)
             p(i,2) =  p(i,1)
             pd(i,2) = pd(i,1)
          ENDDO
       ENDIF
       !
       !---> common spin-orbit integrant V   (average spin directions)
       !                                  SO
       v0(:) = 0.0
       IF (input%jspins.EQ.1) THEN
          v0(1:atoms%jri(ntyp)) = vr(1:atoms%jri(ntyp),1)
          e = enpara%el0(l,ntyp,1)
       ELSE
          DO i = 1,atoms%jri(ntyp)
             v0(i) = (vr(i,1)+vr(i,input%jspins))/2.
          END DO
          e = (enpara%el0(l,ntyp,1)+enpara%el0(l,ntyp,input%jspins))/2.
       END IF

       CALL sointg(ntyp,e,vr,v0,atoms,input,vso)
       IF (spav) THEN
          DO i= 1,atoms%jmtd
             vso(i,1)= (vso(i,1)+vso(i,2))/2.
             vso(i,2)= vso(i,1)
          ENDDO
       ENDIF

       !                        s       s'            .s       s'
       !-->  radial integrals <u  |V  |u  > = rsopp, <u  |V  |u  > = rsopdp etc.
       !                            SO                     SO

       IF (l.GT.0) THEN ! there is no spin-orbit for s-states
          DO i = 1, 2
             DO j = 1, 2
116 117 118 119
                rsopp(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
                rsopdp(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i), p(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
                rsoppd(ntyp,l,i,j) = radso( p(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
                rsopdpd(ntyp,l,i,j) = radso(pd(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
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
             ENDDO
          ENDDO
       ENDIF ! l>0
       !
       !--->  Check for local orbitals with same l
       !
       DO ilo = 1, atoms%nlo(ntyp)
          IF (atoms%llo(ilo,ntyp).EQ.l) THEN

             DO jspin = 1,input%jspins
                e = enpara%ello0(ilo,ntyp,jspin)
                CALL radsra(&
                     e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
                     usdus%ulos(ilo,ntyp,jspin),usdus%dulos(ilo,ntyp,jspin),&
                     nodeu,plo(:,jspin),qlo(:,jspin))

                !+apw+lo
                IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN !  calculate energy derivative (of order atoms%ulo_der) at e
                   ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
                   pqlo(1:atoms%jri(ntyp),1)=plo(1:atoms%jri(ntyp),jspin)
                   pqlo(1:atoms%jri(ntyp),2)=qlo(1:atoms%jri(ntyp),jspin)
                   i = atoms%ulo_der(ilo,ntyp)
                   IF(atoms%l_dulo(ilo,ntyp)) i=1
                   CALL radsrdn(&
                        e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
                        usdus%ulos(ilo,ntyp,jspin),duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
                        pqlo,usdus%dulos(ilo,ntyp,jspin),i)
                   ddn1 = SQRT(ddn1)
                   IF(atoms%l_dulo(ilo,ntyp)) ddn1=1.0
                   plo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)/ddn1
                   qlo(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),2)/ddn1
                   usdus%dulos(ilo,ntyp,jspin) = duds1/ddn1
                   DEALLOCATE (glo,pqlo,filo)
                ENDIF
                !-apw+lo
             ENDDO

             IF (input%jspins.EQ.1) THEN
                plo(1:atoms%jri(ntyp),2) = plo(1:atoms%jri(ntyp),1)
                e = (enpara%ello0(ilo,ntyp,1) + enpara%el0(l,ntyp,1) )/2
             ELSE
                e = (enpara%ello0(ilo,ntyp,1) +  enpara%ello0(ilo,ntyp,input%jspins) +&
                     enpara%el0(l,ntyp,1) + enpara%el0(l,ntyp,input%jspins) )/4
             END IF
             CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
             IF (spav) THEN
                DO i= 1,atoms%jmtd
                   vso(i,1)= (vso(i,1)+vso(i,2))/2.
                   vso(i,2)= vso(i,1)
                ENDDO
             ENDIF

             DO i = 1, 2
                DO j = 1, 2
Daniel Wortmann's avatar
Daniel Wortmann committed
174 175 176 177
                   rsoplop (ntyp,ilo,i,j) = radso(plo(:atoms%jri(ntyp),i),p (:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
                   rsoplopd(ntyp,ilo,i,j) = radso(plo(:atoms%jri(ntyp),i),pd(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
                   rsopplo (ntyp,ilo,i,j) = radso(p (:atoms%jri(ntyp),i),plo(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
                   rsopdplo(ntyp,ilo,i,j) = radso(pd(:atoms%jri(ntyp),i),plo(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp), atoms%rmsh(1,ntyp))
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
                ENDDO
             ENDDO

             DO i = 1,input%jspins
                fint(:) = plo(:,i) *  p(:,i) + qlo(:,i) *  q(:,i)
                CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%uulon(ilo,ntyp,i))
                fint(:) = plo(:,i) * pd(:,i) + qlo(:,i) * qd(:,i)
                CALL intgr0(fint,atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),usdus%dulon(ilo,ntyp,i))
             ENDDO

             DO ilop = 1, atoms%nlo(ntyp)
                IF (atoms%llo(ilop,ntyp).EQ.l) THEN

                   DO jspin = 1,input%jspins
                      e = enpara%ello0(ilop,ntyp,jspin)
                      CALL radsra(&
                           e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
                           ulops,dulops,nodeu,plop(:,jspin),q(:,1))
                      !+apw+lo
                      IF (atoms%l_dulo(ilo,ntyp).OR.atoms%ulo_der(ilo,ntyp).GE.1) THEN ! calculate orthogonal energy derivative at e
                         ALLOCATE (glo(atoms%jmtd,2),pqlo(atoms%jmtd,2),filo(atoms%jmtd,2))
                         pqlo(1:atoms%jri(ntyp),1)=plop(1:atoms%jri(ntyp),jspin)
                         pqlo(1:atoms%jri(ntyp),2)=q(1:atoms%jri(ntyp),1)
                         i = atoms%ulo_der(ilo,ntyp)
                         IF(atoms%l_dulo(ilo,ntyp)) i=1
                         CALL radsrdn(&
                              e,l,vr(:,jspin),atoms%rmsh(1,ntyp),atoms%dx(ntyp),atoms%jri(ntyp),c_light(1.0),&
                              ulops,duds1,ddn1,noded,glo,filo,&!filo is a dummy array&
                              pqlo,dulops,i)
                         plop(1:atoms%jri(ntyp),jspin) = glo(1:atoms%jri(ntyp),1)
                         DEALLOCATE (glo,pqlo,filo)
                      ENDIF
                      !-apw+lo
                   ENDDO

                   IF (input%jspins.EQ.1) THEN
                      plop(1:atoms%jri(ntyp),2) = plop(1:atoms%jri(ntyp),1)
                      e = (enpara%ello0(ilo,ntyp,1) + enpara%ello0(ilop,ntyp,1) )/2
                   ELSE
                      e = (enpara%ello0(ilo,ntyp,1) +  enpara%ello0(ilo,ntyp,input%jspins) +  &
                           enpara%ello0(ilop,ntyp,1) + enpara%ello0(ilop,ntyp,input%jspins) )/4
                   END IF
                   CALL sointg(ntyp,e,vr,v0,atoms,input, vso)
                   IF (spav) THEN
                      DO i= 1,atoms%jmtd
                         vso(i,1)= (vso(i,1)+vso(i,2))/2.
                         vso(i,2)= vso(i,1)
                      ENDDO
                   ENDIF

                   DO i = 1, 2
                      DO j = 1, 2
                         rsoploplop(ntyp,ilo,ilop,i,j) =&
Daniel Wortmann's avatar
Daniel Wortmann committed
231
                              radso(plo(:atoms%jri(ntyp),i),plop(:atoms%jri(ntyp),j),vso(:atoms%jri(ntyp),i),atoms%dx(ntyp),atoms%rmsh(1,ntyp))
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 271 272 273 274
                      ENDDO
                   ENDDO

                ENDIF
             ENDDO

          ENDIF
       ENDDO ! end of lo-loop

    ENDDO ! end of l-loop

    DEALLOCATE ( p,pd,q,qd,plo,plop,qlo,fint,v0,vso )
    !      rsoplop (:,:,:,:) = 0.0
    !      rsoplopd(:,:,:,:) = 0.0
    !      rsopplo (:,:,:,:) = 0.0
    !      rsopdplo(:,:,:,:) = 0.0
    !      rsoploplop(:,:,:,:,:) = 0.0

  END SUBROUTINE sorad
  !--------------------------------------------------------------------
  REAL FUNCTION radso(a,b,vso,dx,r0)
    !
    !     compute radial spin-orbit integrals
    !
    USE m_intgr, ONLY : intgr0
    IMPLICIT NONE
    !
    !     .. Scalar Arguments ..
    REAL,    INTENT (IN) :: r0,dx
    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN) :: a(:),b(:),vso(:)
    !     ..
    !     .. Local Arrays ..
    REAL q(size(a))
    !     ..
    q = a*b*vso
    CALL intgr0(q,r0,dx,size(a),radso)

    RETURN
  END FUNCTION radso
  !--------------------------------------------------------------------
END MODULE m_sorad