vacfun.f90 9.2 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
MODULE m_vacfun
  use m_juDFT
CONTAINS
  SUBROUTINE vacfun(&
       vacuum,dimension,stars, jsp,input,noco,ipot,&
       sym, cell,ivac,evac,bkpt, vxy,vz,kvac1,kvac2,nv2,&
       tuuv,tddv,tudv,tduv,uz,duz,udz,dudz,ddnv,wronk)
    !*********************************************************************
    !     determines the necessary values and derivatives on the vacuum
    !     boundary (ivac=1 upper vacuum; ivac=2, lower) for energy
    !     parameter evac.  also sets up the 2d hamiltonian matrices
    !     necessary to update the full hamiltonian matrix.
    !               m. weinert
    !*********************************************************************

    USE m_intgr, ONLY : intgz0
    USE m_vacuz
    USE m_vacudz
    USE m_types
    IMPLICIT NONE

    TYPE(t_dimension),INTENT(IN)   :: dimension
    TYPE(t_input),INTENT(IN)       :: input
    TYPE(t_vacuum),INTENT(IN)      :: vacuum
    TYPE(t_noco),INTENT(IN)        :: noco
    TYPE(t_sym),INTENT(IN)         :: sym
    TYPE(t_stars),INTENT(IN)       :: stars
    TYPE(t_cell),INTENT(IN)        :: cell
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: jsp ,ivac,ipot
    REAL,    INTENT (OUT) :: wronk
    !     ..
    !     .. Array Arguments ..
    INTEGER, INTENT (IN) :: nv2(dimension%jspd)
    INTEGER, INTENT (IN) :: kvac1(dimension%nv2d,dimension%jspd),kvac2(dimension%nv2d,dimension%jspd)
    COMPLEX, INTENT (IN) :: vxy(vacuum%nmzxyd,stars%n2d-1)
    COMPLEX, INTENT (OUT):: tddv(dimension%nv2d,dimension%nv2d),tduv(dimension%nv2d,dimension%nv2d)
    COMPLEX, INTENT (OUT):: tudv(dimension%nv2d,dimension%nv2d),tuuv(dimension%nv2d,dimension%nv2d)
    REAL,    INTENT (IN) :: vz(vacuum%nmzd,2,4) ,evac(2,dimension%jspd)
    REAL,    INTENT (IN) :: bkpt(3) 
    REAL,    INTENT (OUT):: udz(dimension%nv2d,dimension%jspd),uz(dimension%nv2d,dimension%jspd)
    REAL,    INTENT (OUT):: dudz(dimension%nv2d,dimension%jspd),duz(dimension%nv2d,dimension%jspd)
    REAL,    INTENT (OUT):: ddnv(dimension%nv2d,dimension%jspd)
    !     ..
    !     .. Local Scalars ..
    REAL ev,phase,scale,xv,yv,vzero
    INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin,jsp1,jsp2
    LOGICAL tail
    !     ..
    !     .. Local Arrays ..
    REAL u(vacuum%nmzd,dimension%nv2d,dimension%jspd),ud(vacuum%nmzd,dimension%nv2d,dimension%jspd)
    REAL v(3),x(vacuum%nmzd), qssbti(2,2)
    !     ..
    tail = .true.
    np1 = vacuum%nmzxy + 1
    !--->    wronksian for the schrodinger equation given by an identity
    wronk = 2.0
    !---> setup the spin-spiral q-vector
    qssbti(1,1) = - noco%qss(1)/2
    qssbti(2,1) = - noco%qss(2)/2
    qssbti(1,2) = + noco%qss(1)/2
    qssbti(2,2) = + noco%qss(2)/2
    !--->    generate basis functions for each 2-d k+g
    DO jspin = 1,input%jspins
       DO  ik = 1,nv2(jspin)
          v(1) = bkpt(1) + kvac1(ik,jspin) + qssbti(1,jspin)
          v(2) = bkpt(2) + kvac2(ik,jspin) + qssbti(2,jspin)
          v(3) = 0.0
          ev = evac(ivac,jspin) - 0.5*dot_product(v,matmul(v,cell%bbmat))
          vzero = vz(vacuum%nmzd,ivac,jspin)
          CALL vacuz(ev,vz(1,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
               uz(ik,jspin),duz(ik,jspin),u(1,ik,jspin))
          CALL vacudz(ev,vz(1,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
               udz(ik,jspin),dudz(ik,jspin),ddnv(ik,jspin),&
               ud(1,ik,jspin),duz(ik,jspin),u(1,ik,jspin))
          !--->       make sure the solutions satisfy the wronksian
          scale = wronk/ (udz(ik,jspin)*duz(ik,jspin)-&
               &                         dudz(ik,jspin)*uz(ik,jspin))
          udz(ik,jspin) = scale*udz(ik,jspin)
          dudz(ik,jspin) = scale*dudz(ik,jspin)
          ddnv(ik,jspin) = scale*ddnv(ik,jspin)
          ud(:,ik,jspin) = scale*ud(:,ik,jspin)
       enddo
    ENDDO
    !--->    set up the tuuv, etc. matrices
    IF (noco%l_noco) THEN
       IF (ipot.EQ.1) THEN
          jsp1 = 1
          jsp2 = 1
       ELSEIF (ipot.EQ.2) THEN
          jsp1 = 2
          jsp2 = 2
       ELSEIF (ipot.EQ.3) THEN
          jsp1 = 2
          jsp2 = 1
       ENDIF
    ELSE
       jsp1 = jsp
       jsp2 = jsp
    ENDIF
    DO  ik = 1,nv2(jsp1)
       DO  jk = 1,nv2(jsp2)

          !--->     determine the warping component of the potential
          i1 = kvac1(ik,jsp1) - kvac1(jk,jsp2)
          i2 = kvac2(ik,jsp1) - kvac2(jk,jsp2)
          i3 = 0
          ind3 = stars%ig(i1,i2,i3)
          IF (ind3.EQ.0) CYCLE
          phase = stars%rgphs(i1,i2,i3)
          ind2 = stars%ig2(ind3)
          IF (ind2.EQ.0) THEN
             WRITE (16,FMT=8000) ik,jk
             WRITE (6,FMT=8000) ik,jk
8000         FORMAT (' **** error in map2 for 2-d stars',2i5)
             CALL juDFT_error("error in map2 for 2-d stars",calledby ="vacfun")
          END IF
          !--->     get the proper warping index (vxy starts with the 2nd star)
          ind2 = ind2 - 1
          IF (ind2.NE.0) THEN
             !--->       only the warping part, 1st star (G=0) is done later

             !--->       obtain the warping matrix elements
             !--->       note that the tail correction (tail=.true.) is included for
             !--->       the integrals, i.e. the integrand is from infinity inward

             !--->       tuuv
             DO  i = 1,vacuum%nmzxy
                x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
                x(np1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tuuv(ik,jk) = phase*cmplx(xv,yv)

             !--->       tddv
             DO  i = 1,vacuum%nmzxy
                x(np1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
                x(np1-i) =ud(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tddv(ik,jk) = phase*cmplx(xv,yv)

             !--->       tudv
             DO  i = 1,vacuum%nmzxy
                x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*real(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
                x(np1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*aimag(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tudv(ik,jk) = phase*cmplx(xv,yv)

             !--->       tduv
             DO  i = 1,vacuum%nmzxy
                x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*real(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
                x(np1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*aimag(vxy(i,ind2))
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tduv(ik,jk) = phase*cmplx(xv,yv)

          ELSE

             !--->       diagonal (film muffin-tin) terms
             IF ((ipot.EQ.1) .OR. (ipot.EQ.2)) THEN
                tuuv(ik,ik) = cmplx(evac(ivac,jsp1),0.0)
                tddv(ik,ik) = cmplx(evac(ivac,jsp1)*ddnv(ik,jsp1),0.0)
                tudv(ik,ik) = cmplx(0.5,0.0)
                tduv(ik,ik) = cmplx(0.5,0.0)
             ELSE

                !--->          tuuv
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = u(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tuuv(ik,jk) = cmplx(xv,yv)

                !--->          tddv
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tddv(ik,jk) = cmplx(xv,yv)

                !--->          tudv
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,3)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = u(i,ik,jsp1)*ud(i,jk,jsp2)*vz(i,ivac,4)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tudv(ik,jk) = cmplx(xv,yv)

                !--->          tduv
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,3)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
                   x(vacuum%nmz+1-i) = ud(i,ik,jsp1)*u(i,jk,jsp2)*vz(i,ivac,4)
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tduv(ik,jk) = cmplx(xv,yv)
             ENDIF

          ENDIF
       enddo
    enddo

  END SUBROUTINE vacfun
END MODULE m_vacfun