mkgxyz3.f90 8.14 KB
Newer Older
1 2 3 4 5
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
6
MODULE m_mkgxyz3
7
  USE m_judft
8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
  !.....------------------------------------------------------------------
  !c     by use of cartesian x,y,z components of charge density gradients,
  !c     make the quantities
  !cc      agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr
  !cc    used to calculate gradient contribution to xc potential and
  !cc    energy.
  !c.....------------------------------------------------------------------
CONTAINS
  SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy,grad)
    USE m_types
    IMPLICIT NONE
    REAL,    INTENT (IN) :: vl(:,:)
    REAL,    INTENT (IN) :: dvx(:,:),dvy(:,:),dvz(:,:)
    REAL, INTENT (IN) :: dvxx(:,:),dvyy(:,:),dvzz(:,:)
    REAL, INTENT (IN) :: dvyz(:,:),dvzx(:,:),dvxy(:,:)

    TYPE(t_gradients),INTENT(INOUT)::grad

    REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt,&
         vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu,&
         vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvzxd,dvxyd,&
         dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
30 31
         dagrzu,dzdx,dzdy,dzdz
    REAL, PARAMETER  :: sml = 1.e-14
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
    INTEGER i,js,jspins,nsp

    nsp=SIZE(dvx,1)
    jspins=SIZE(dvx,2)

    IF (ALLOCATED(grad%gr)) THEN
       !      Gradients for libxc
       DO js=1,jspins
          DO i=1,nsp
             grad%gr(:,i,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
          ENDDO
       END DO
       IF(ALLOCATED(grad%sigma)) THEN
          !Use only contracted gradients for libxc
          IF (jspins==1) THEN
             DO i=1,nsp
Matthias Redies's avatar
Matthias Redies committed
48
                grad%sigma(1,i) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
49 50 51
             ENDDO
          ELSE
             DO i=1,nsp
Matthias Redies's avatar
Matthias Redies committed
52 53 54
                grad%sigma(1,i) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
                grad%sigma(2,i) = dvx(i,1)*dvx(i,2) + dvy(i,1)*dvy(i,2) + dvz(i,1)*dvz(i,2)
                grad%sigma(3,i) = dvx(i,2)*dvx(i,2) + dvy(i,2)*dvy(i,2) + dvz(i,2)*dvz(i,2)
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
             ENDDO
          ENDIF
       END IF
       IF(ALLOCATED(grad%laplace)) THEN
          DO js=1,jspins
             DO i=1,nsp
                grad%laplace(i,js)= dvxx(i,js)+dvyy(i,js)+dvzz(i,js)
             ENDDO
          ENDDO
       ENDIF
       RETURN
    ENDIF

    IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
    
    DO i = 1,size(grad%agrt)
Matthias Redies's avatar
Matthias Redies committed
71 72 73
       grad%agrt(i)  = 0.0
       grad%agru(i)  = 0.0
       grad%agrd(i)  = 0.0
74 75 76
       grad%gggrt(i) = 0.0
       grad%gggru(i) = 0.0
       grad%gggrd(i) = 0.0
Matthias Redies's avatar
Matthias Redies committed
77 78 79 80
       grad%gzgr(i)  = 0.0
       grad%g2rt(i)  = 0.0
       grad%g2ru(i)  = 0.0
       grad%g2rd(i)  = 0.0
81 82 83 84
    ENDDO

    IF (jspins.eq.1) THEN

Matthias Redies's avatar
Matthias Redies committed
85
       DO i = 1,nsp
86

Matthias Redies's avatar
Matthias Redies committed
87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107
          vlu   = max(vl(i,1)/2,sml)
          dvxu  = dvx(i,1)/2
          dvyu  = dvy(i,1)/2
          dvzu  = dvz(i,1)/2
          dvxxu = dvxx(i,1)/2
          dvyyu = dvyy(i,1)/2
          dvzzu = dvzz(i,1)/2
          dvyzu = dvyz(i,1)/2
          dvzxu = dvzx(i,1)/2
          dvxyu = dvxy(i,1)/2

          vld   = vlu
          dvxd  = dvxu
          dvyd  = dvyu
          dvzd  = dvzu
          dvxxd = dvxxu
          dvyyd = dvyyu
          dvzzd = dvzzu
          dvyzd = dvyzu
          dvzxd = dvzxu
          dvxyd = dvxyu
108 109 110 111


          vlt = vlu + vld

Matthias Redies's avatar
Matthias Redies committed
112 113 114
          dvxt  = dvxu  + dvxd
          dvyt  = dvyu  + dvyd
          dvzt  = dvzu  + dvzd
115 116 117 118 119 120 121
          dvxxt = dvxxu + dvxxd
          dvyyt = dvyyu + dvyyd
          dvzzt = dvzzu + dvzzd
          dvyzt = dvyzu + dvyzd
          dvzxt = dvzxu + dvzxd
          dvxyt = dvxyu + dvxyd

122
          !         agr: abs(grad(ro)), t,u,d for total, up and down.
123

Matthias Redies's avatar
Matthias Redies committed
124 125 126
          grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
          grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
          grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
127

Matthias Redies's avatar
Matthias Redies committed
128 129 130
          dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
          dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
          dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
131

Matthias Redies's avatar
Matthias Redies committed
132 133 134
          dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
          dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
          dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
135

Matthias Redies's avatar
Matthias Redies committed
136 137 138
          dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
          dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
          dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
139

140 141 142
          grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
          grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
          grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
143

144
          !         dzdx=d(zeta)/dx,..
145 146 147 148 149

          dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
          dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
          dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2

150
          !         gzgr=grad(zeta)*grad(ro).
151

152
          grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
153

154
          !         g2r: grad(grad(ro))
155

156 157 158
          grad%g2rt(i) = dvxxt + dvyyt + dvzzt
          grad%g2ru(i) = dvxxu + dvyyu + dvzzu
          grad%g2rd(i) = dvxxd + dvyyd + dvzzd
159 160


Matthias Redies's avatar
Matthias Redies committed
161
       ENDDO
162

163
    ELSE
164

Matthias Redies's avatar
Matthias Redies committed
165
       DO i = 1,nsp
166

Matthias Redies's avatar
Matthias Redies committed
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187
          vlu   = max(vl(i,1),sml)
          dvxu  = dvx(i,1)
          dvyu  = dvy(i,1)
          dvzu  = dvz(i,1)
          dvxxu = dvxx(i,1)
          dvyyu = dvyy(i,1)
          dvzzu = dvzz(i,1)
          dvyzu = dvyz(i,1)
          dvzxu = dvzx(i,1)
          dvxyu = dvxy(i,1)

          vld   = max(vl(i,jspins),sml)
          dvxd  = dvx(i,jspins)
          dvyd  = dvy(i,jspins)
          dvzd  = dvz(i,jspins)
          dvxxd = dvxx(i,jspins)
          dvyyd = dvyy(i,jspins)
          dvzzd = dvzz(i,jspins)
          dvyzd = dvyz(i,jspins)
          dvzxd = dvzx(i,jspins)
          dvxyd = dvxy(i,jspins)
188 189 190

          vlt = vlu + vld

Matthias Redies's avatar
Matthias Redies committed
191 192 193
          dvxt  = dvxu  + dvxd
          dvyt  = dvyu  + dvyd
          dvzt  = dvzu  + dvzd
194 195 196 197 198 199 200
          dvxxt = dvxxu + dvxxd
          dvyyt = dvyyu + dvyyd
          dvzzt = dvzzu + dvzzd
          dvyzt = dvyzu + dvyzd
          dvzxt = dvzxu + dvzxd
          dvxyt = dvxyu + dvxyd

201
          !c         agr: abs(grad(ro)), t,u,d for total, up and down.
202

Matthias Redies's avatar
Matthias Redies committed
203 204 205
          grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
          grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
          grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
206

Matthias Redies's avatar
Matthias Redies committed
207 208 209
          dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
          dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
          dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
210

Matthias Redies's avatar
Matthias Redies committed
211 212 213
          dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
          dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
          dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
214

Matthias Redies's avatar
Matthias Redies committed
215 216 217
          dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
          dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
          dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
218

219 220 221
          grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
          grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
          grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
222

223
          !c         dzdx=d(zeta)/dx,..
224 225 226 227 228

          dzdx = (dvxu-dvxd)/vlt -  (vlu-vld)*dvxt/vlt**2
          dzdy = (dvyu-dvyd)/vlt -  (vlu-vld)*dvyt/vlt**2
          dzdz = (dvzu-dvzd)/vlt -  (vlu-vld)*dvzt/vlt**2

229
          !c         gzgr=grad(zeta)*grad(ro).
230

231
          grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
232

233
          !c         g2r: grad(grad(ro))
234

235 236 237
          grad%g2rt(i) = dvxxt + dvyyt + dvzzt
          grad%g2ru(i) = dvxxu + dvyyu + dvzzu
          grad%g2rd(i) = dvxxd + dvyyd + dvzzd
238

Matthias Redies's avatar
Matthias Redies committed
239
       ENDDO
240

241 242 243 244
    ENDIF

  END SUBROUTINE mkgxyz3
END MODULE m_mkgxyz3