vacfun.f90 9.31 KB
Newer Older
1 2 3 4
MODULE m_vacfun
  use m_juDFT
CONTAINS
  SUBROUTINE vacfun(&
5
       vacuum,stars,input,nococonv,jspin1,jspin2,&
6
       sym, cell,ivac,evac,bkpt, vxy,vz,kvac,nv2,&
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
       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_input),INTENT(IN)       :: input
    TYPE(t_vacuum),INTENT(IN)      :: vacuum
24
    TYPE(t_nococonv),INTENT(IN)   :: nococonv
25 26 27 28 29
    TYPE(t_sym),INTENT(IN)         :: sym
    TYPE(t_stars),INTENT(IN)       :: stars
    TYPE(t_cell),INTENT(IN)        :: cell
    !     ..
    !     .. Scalar Arguments ..
30
    INTEGER, INTENT (IN) :: ivac,jspin1,jspin2
31 32 33
    REAL,    INTENT (OUT) :: wronk
    !     ..
    !     .. Array Arguments ..
34
    INTEGER, INTENT (IN) :: nv2(:)!(input%jspins)
35
    INTEGER, INTENT (IN) :: kvac(:,:,:)!(2,lapw%dim_nv2d(),input%jspins)
36
    COMPLEX, INTENT (IN) :: vxy(:,:,:,:) !(vacuum%nmzxyd,stars%ng2-1,nvac,:)
37 38
    COMPLEX, INTENT (OUT):: tddv(:,:),tduv(:,:)!(lapw%dim_nv2d(),lapw%dim_nv2d())
    COMPLEX, INTENT (OUT):: tudv(:,:),tuuv(:,:)!(lapw%dim_nv2d(),lapw%dim_nv2d())
39
    REAL,    INTENT (IN) :: vz(:,:,:) !(vacuum%nmzd,2,4) ,
40
    REAL,    INTENT (IN) :: evac(:,:)!(2,input%jspins)
41
    REAL,    INTENT (IN) :: bkpt(3)
42 43 44
    REAL,    INTENT (OUT):: udz(:,:),uz(:,:)!(lapw%dim_nv2d(),input%jspins)
    REAL,    INTENT (OUT):: dudz(:,:),duz(:,:)!(lapw%dim_nv2d(),input%jspins)
    REAL,    INTENT (OUT):: ddnv(:,:)!(lapw%dim_nv2d(),input%jspins)
45 46
    !     ..
    !     .. Local Scalars ..
47
    REAL ev,scale,xv,yv,vzero,fac
48
    COMPLEX phase
49
    INTEGER i,i1,i2,i3,ik,ind2,ind3,jk,np1,jspin,ipot
50 51 52
    LOGICAL tail
    !     ..
    !     .. Local Arrays ..
53
    REAL u(vacuum%nmzd,size(duz,1),input%jspins),ud(vacuum%nmzd,size(duz,1),input%jspins)
54 55
    REAL v(3),x(vacuum%nmzd), qssbti(2,2)
    !     ..
56
    fac=MERGE(1.0,-1.0,jspin1>=jspin2)
57 58
    ipot=MERGE(jspin1,3,jspin1==jspin2)

59 60
    tuuv=0.0;tudv=0.0;tddv=0.0;tduv=0.0
    udz=0.0;duz=0.0;ddnv=0.0;udz=0.;uz=0.
61 62 63 64 65
    tail = .true.
    np1 = vacuum%nmzxy + 1
    !--->    wronksian for the schrodinger equation given by an identity
    wronk = 2.0
    !---> setup the spin-spiral q-vector
66 67
    qssbti(1:2,1) = - nococonv%qss(1:2)/2
    qssbti(1:2,2) = + nococonv%qss(1:2)/2
68
    !--->    generate basis functions for each 2-d k+g
69
    DO jspin = MIN(jspin1,jspin2),MAX(jspin1,jspin2)
70
       DO  ik = 1,nv2(jspin)
71
          v(1:2) = bkpt(1:2) + kvac(:,ik,jspin) + qssbti(1:2,jspin)
72 73 74
          v(3) = 0.0
          ev = evac(ivac,jspin) - 0.5*dot_product(v,matmul(v,cell%bbmat))
          vzero = vz(vacuum%nmzd,ivac,jspin)
75
          CALL vacuz(ev,vz(1:,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
76
               uz(ik,jspin),duz(ik,jspin),u(1,ik,jspin))
77
          CALL vacudz(ev,vz(1:,ivac,jspin),vzero,vacuum%nmz,vacuum%delz,&
78 79 80 81 82
               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))
83
          udz(ik,jspin)  = scale*udz(ik,jspin)
84 85 86 87 88 89
          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
90 91 92 93 94
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP& COLLAPSE(2) &
    !$OMP& SHARED(tuuv,tddv,tudv,tduv,ddnv,vz) &
    !$OMP& SHARED(stars,jspin1,jspin2,evac,nv2,kvac,vacuum,u,vxy,tail,fac,np1,ivac,ipot,ud) &
    !$OMP& PRIVATE(i1,i2,i3,ind3,phase,ind2,x,xv,yv)
95 96
    DO  ik = 1,nv2(jspin1)
       DO  jk = 1,nv2(jspin2)
97 98

          !--->     determine the warping component of the potential
99 100
          i1 = kvac(1,ik,jspin1) - kvac(1,jk,jspin2)
          i2 = kvac(2,ik,jspin1) - kvac(2,jk,jspin2)
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121
          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 (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
122
                x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
123 124 125
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
126
                x(np1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
127 128 129 130 131 132
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tuuv(ik,jk) = phase*cmplx(xv,yv)

             !--->       tddv
             DO  i = 1,vacuum%nmzxy
133
                x(np1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*REAL(vxy(i,ind2,ivac,ipot))
134 135 136
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
137
                x(np1-i) =ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
138 139 140 141 142 143
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tddv(ik,jk) = phase*cmplx(xv,yv)

             !--->       tudv
             DO  i = 1,vacuum%nmzxy
144
                x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
145 146 147
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
148
                x(np1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
149 150 151 152 153 154
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tudv(ik,jk) = phase*cmplx(xv,yv)

             !--->       tduv
             DO  i = 1,vacuum%nmzxy
155
                x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*real(vxy(i,ind2,ivac,ipot))
156 157 158
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,xv,tail)
             DO  i = 1,vacuum%nmzxy
159
                x(np1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*AIMAG(vxy(i,ind2,ivac,ipot))
160 161 162 163 164 165
             enddo
             CALL intgz0(x,vacuum%delz,vacuum%nmzxy,yv,tail)
             tduv(ik,jk) = phase*cmplx(xv,yv)

          ELSE
             !--->       diagonal (film muffin-tin) terms
166 167 168
             IF (jspin1==jspin2) THEN
                tuuv(ik,ik) = cmplx(evac(ivac,jspin1),0.0)
                tddv(ik,ik) = cmplx(evac(ivac,jspin1)*ddnv(ik,jspin1),0.0)
169 170 171 172 173
                tudv(ik,ik) = cmplx(0.5,0.0)
                tduv(ik,ik) = cmplx(0.5,0.0)
             ELSE
                !--->          tuuv
                DO i = 1,vacuum%nmz
174
                   x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
175 176 177
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
178
                   x(vacuum%nmz+1-i) = u(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
179 180 181 182 183 184
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tuuv(ik,jk) = cmplx(xv,yv)

                !--->          tddv
                DO i = 1,vacuum%nmz
185
                   x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
186 187 188
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
189
                   x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
190 191 192 193 194 195
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tddv(ik,jk) = cmplx(xv,yv)

                !--->          tudv
                DO i = 1,vacuum%nmz
196
                   x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*vz(i,ivac,3)
197 198 199
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
200
                   x(vacuum%nmz+1-i) = u(i,ik,jspin1)*ud(i,jk,jspin2)*fac*vz(i,ivac,4)
201 202 203 204 205 206
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tudv(ik,jk) = cmplx(xv,yv)

                !--->          tduv
                DO i = 1,vacuum%nmz
207
                   x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*vz(i,ivac,3)
208 209 210
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,xv,tail)
                DO i = 1,vacuum%nmz
211
                   x(vacuum%nmz+1-i) = ud(i,ik,jspin1)*u(i,jk,jspin2)*fac*vz(i,ivac,4)
212 213 214 215 216 217
                ENDDO
                CALL intgz0(x,vacuum%delz,vacuum%nmz,yv,tail)
                tduv(ik,jk) = cmplx(xv,yv)
             ENDIF

          ENDIF
218 219 220
       ENDDO
    ENDDO
    !$OMP END PARALLEL DO
221 222 223

  END SUBROUTINE vacfun
END MODULE m_vacfun