metric.f90 7.82 KB
Newer Older
1 2 3 4 5 6 7
MODULE m_metric
  USE m_juDFT
  !*********************************************************
  !     multiplicates the vector s_in with the metric G
  !     output vector sout
  !********************************************************* 
CONTAINS
8 9
  SUBROUTINE metric(cell,atoms,vacuum,sphhar,input,noco,stars,sym,oneD,&
                    mmap,nmaph,mapmt,mapvac2,s_in,sout,lpot) 
10 11 12
    USE m_metrz0
    USE m_convol
    USE m_types
13

14
    IMPLICIT NONE
15

16
    TYPE(t_oneD),INTENT(IN)   :: oneD
17 18
    TYPE(t_input),INTENT(IN)  :: input
    TYPE(t_vacuum),INTENT(IN) :: vacuum
19
    TYPE(t_noco),INTENT(IN)   :: noco
20 21
    TYPE(t_sym),INTENT(IN)    :: sym
    TYPE(t_stars),INTENT(IN)  :: stars
22
    TYPE(t_cell),INTENT(IN)   :: cell
23 24 25 26 27 28 29
    TYPE(t_sphhar),INTENT(IN) :: sphhar
    TYPE(t_atoms),INTENT(IN)  :: atoms

    ! Scalar Arguments
    INTEGER, INTENT (IN)          :: mmap
    INTEGER, INTENT (IN)          :: mapmt,mapvac2
    INTEGER, INTENT (IN)          :: nmaph      
30
    LOGICAL, OPTIONAL,INTENT (IN) :: lpot !do we mix a potential??
31 32 33 34 35 36

    ! Array Arguments
    REAL,    INTENT (IN)          :: s_in(mmap)
    REAL,    INTENT (OUT)         :: sout(mmap)

    ! Local Scalars
37
    INTEGER              :: imap,ivac,iz,j,js,k2,l,n,iv2c,iv2,na,ioff,i
38 39 40 41
    REAL                 :: dvol,dxn,dxn2,dxn4,volnstr2
    LOGICAL              :: l_pot

    ! Local Arrays
42 43
    REAL,    ALLOCATABLE :: g(:),wght(:)
    COMPLEX, ALLOCATABLE :: ag3(:),fg3(:)
44 45 46 47 48 49 50

    ! calculate the coefficients of the g-matrix
    ! for the m.t. and vacuum region

    l_pot = .FALSE.
    IF (PRESENT(lpot)) l_pot = lpot ! for potential mixing

51
    ALLOCATE (g(mmap),wght(vacuum%nmzd),ag3(stars%ng3),fg3(stars%ng3))
52 53 54 55 56 57 58
    g = 0.0

    imap = 2 * stars%ng3 ! complex values without invs
    IF (sym%invs) imap = stars%ng3 ! only real values with invs
    iv2 = 2
    IF (sym%invs2) iv2 = 1

59 60 61 62 63
    ! metric for MT is r^2 dr/di di = r^3 dx ; divided by r^4 to 
    ! compensate use of n(r) * r^2 in array rho
    ! simpson integration used, weights for first and last point: 1
    ! weights forthe rest alternating: 2 or 4    
    na = 1
64 65 66 67 68 69
    DO n = 1, atoms%ntype
       dxn = atoms%neq(n) * atoms%dx(n) / 3.0e0
       dxn2 =2.0e0 * dxn
       dxn4 =4.0e0 * dxn

       DO l = 0, sphhar%nlh(atoms%ntypsy(na))
70
          imap = imap + 1
71
          g(imap) = dxn / atoms%rmsh(1,n)
72
          IF (.NOT.l_pot) THEN
73
             DO j = 2, atoms%jri(n) - 1, 2
74
                imap = imap + 2
75 76 77
                g(imap-1) = dxn4 / atoms%rmsh(j,n) 
                g(imap) = dxn2 / atoms%rmsh(j+1,n) 
             END DO
78 79
             ! CHANGE JR 96/12/01
             ! take care when jri(n) is even
80 81
             imap = imap + 1 - MOD(atoms%jri(n),2)
             g(imap) = dxn / atoms%rmsh(atoms%jri(n),n)
82
          ELSE
83 84
             ! for the potential multiply by r^4
             DO j = 2, atoms%jri(n) - 1, 2
85
                imap = imap + 2
86 87 88 89 90 91 92
                g(imap-1) = dxn4 * atoms%rmsh(j,n)**3 
                g(imap) = dxn2 * atoms%rmsh(j+1,n)**3
             END DO
             imap = imap + 1 - MOD(atoms%jri(n),2)
             g(imap) = dxn * atoms%rmsh(atoms%jri(n),n)**3
          END IF
       END DO
93
       na = na + atoms%neq(n)
94 95
    END DO

96 97 98
    ! vacuum contribution
    IF (input%film) THEN
       dvol = cell%area*vacuum%delz
99 100 101 102 103
       ! nvac=1 if (zrfs.or.invs)
       IF (vacuum%nvac.EQ.1) dvol = dvol + dvol
       IF (oneD%odi%d1) dvol = cell%area*vacuum%delz
       DO ivac = 1, vacuum%nvac
          ! G||=0 components
104
          !
105 106
          ! use 7-point simpson integration in accordance to intgz0.f
          ! calculate weights for integration
107
          CALL metr_z0(vacuum%nmz,wght)
108
          DO iz = 1, vacuum%nmz
109 110
             imap = imap + 1
             IF (oneD%odi%d1) THEN
111
                g(imap) = wght(iz) * dvol * (cell%z1+(iz-1)*vacuum%delz)
112
             ELSE
113 114 115 116 117 118
                g(imap) = wght(iz) * dvol
             END IF
          END DO
          ! G||.ne.0 components
          !
          ! calculate weights for integration
119
          CALL metr_z0(vacuum%nmzxy,wght)
120 121
          DO iv2c = 1, iv2
             DO k2 = 1, oneD%odi%nq2 - 1
122 123 124
                IF (oneD%odi%d1) THEN
                   DO iz = 1,vacuum%nmzxy
                      imap = imap + 1
125 126
                      g(imap) = wght(iz) * oneD%odi%nst2(k2) * dvol * (cell%z1+(iz-1)*vacuum%delz)
                   END DO
127
                ELSE
128 129
                   volnstr2 = dvol * stars%nstr2(k2)
                   DO iz = 1, vacuum%nmzxy
130
                      imap = imap + 1
131 132
                      g(imap) = wght(iz) * volnstr2
                   END DO
133
                END IF
134 135 136
             END DO
          END DO
       END DO
137
    END IF
138

139 140
    ! Apply metric to interstitial region (metric here = step function)
    ! + multiply metric g with s_in for MT and vacuum contributions (store in sout)
141 142
    DO js = 1, input%jspins
       ! map s_in on a complex help array ag3
143
       IF (sym%invs) THEN
144
          DO imap = 1, stars%ng3
145
             ag3(imap) = CMPLX(s_in(imap+nmaph*(js-1)),0.0)
146
          END DO
147
       ELSE
148 149 150
          DO imap = 1, stars%ng3
             ag3(imap) = CMPLX(s_in(imap+nmaph*(js-1)),s_in(imap+stars%ng3+nmaph*(js-1)))
          END DO
151
       ENDIF
152 153 154

       CALL convol(stars,fg3,ag3,stars%ufft)

155
       IF (sym%invs) THEN
156
          ! interstitial
157 158 159
          DO imap = 1, stars%ng3
             sout(imap+nmaph*(js-1)) = cell%omtil * REAL(fg3(imap))
          END DO
160
          ! MT + vacuum
161 162 163
          DO imap = stars%ng3+1, nmaph
             sout(imap+nmaph*(js-1)) = g(imap) * s_in(imap+nmaph*(js-1))
          END DO
164
       ELSE
165
          ! interstitial
166 167 168 169
          DO imap = 1, stars%ng3
             sout(imap+nmaph*(js-1))           = cell%omtil * REAL(fg3(imap))
             sout(imap+stars%ng3+nmaph*(js-1)) = cell%omtil * AIMAG(fg3(imap))
          END DO
170
          ! MT + vacuum
171 172 173 174 175
          DO imap = 2 * stars%ng3 + 1, nmaph
             sout(imap+nmaph*(js-1)) = g(imap) * s_in(imap+nmaph*(js-1))
          END DO
       END IF
    END DO
176 177

    IF (noco%l_noco) THEN
178 179 180 181 182 183 184 185 186 187
       DO imap = 1, stars%ng3
          ag3(imap) = CMPLX(s_in(2*nmaph + imap),s_in(2*nmaph + stars%ng3 + imap))
       END DO

       CALL convol(stars,fg3,ag3,stars%ufft)

       DO imap = 1, stars%ng3
          sout(2*nmaph + imap)             = cell%omtil * REAL(fg3(imap))
          sout(2*nmaph + stars%ng3 + imap) = cell%omtil * AIMAG(fg3(imap))
       END DO
188
       IF (input%film) THEN
189 190 191 192 193 194 195
          ! js runs over the real and imaginary part of the vacuum density
          ! coefficients (not the spin).
          ioff = 2*stars%ng3 + mapmt ! (for complex values)
          IF (sym%invs) ioff = stars%ng3 + mapmt ! (real values if sym%invs)

          DO js = 1, 2
             DO imap = 1, mapvac2 / 2
196
                sout(2*nmaph + 2*stars%ng3 + mapvac2/2*(js-1) + imap) = &
197 198 199 200
                   g(ioff + imap) * s_in(2*nmaph + 2*stars%ng3 + mapvac2/2*(js-1) + imap)
             END DO
          END DO
       END IF
201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
       !mt offdiagonal part
       imap=2*nmaph + 2*stars%ng3 + mapvac2/2*(js-1) + imap-1
       ioff=MERGE(stars%ng3,2*stars%ng2,sym%invs)+1
       j=0
       IF (noco%l_mtnocopot) THEN
          na = 1
          DO n = 1,atoms%ntype
             DO l = 0,sphhar%nlh(atoms%ntypsy(na))
                DO i = 1,atoms%jri(n)
                   j = j + 1
                   sout(imap+j) = g(ioff+(j-1)/2)*s_in(imap+j)
                   j = j + 1
                   sout(imap+j) = g(ioff+(j-1)/2)*s_in(imap+j)
                END DO
             END DO
             na = na + atoms%neq(n)
          END DO
       END IF
219
    END IF
220

221 222 223
    ! density matrix
    IF (atoms%n_u > 0) THEN
       j = input%jspins * nmaph
224
       IF (noco%l_noco) THEN
225
          j = j + 2 * stars%ng3
226 227
          IF (input%film) THEN
             j = j + mapvac2
228 229
          END IF
       END IF
230 231
       DO imap = j+1, j+49*2*input%jspins*atoms%n_u
          sout(imap) = s_in(imap)
232 233 234
       END DO
    END IF

235 236 237 238
    DEALLOCATE (g,wght,ag3,fg3)

  END SUBROUTINE metric
END MODULE m_metric