mkgxyz3.f 7.51 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 7 8 9 10 11 12 13 14 15
      MODULE m_mkgxyz3
c.....------------------------------------------------------------------
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(
16
     >                   ndm,jsdm,ng3,jspins,vl,
17
     >                   dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy,
18 19
     <                   grad)
      USE m_types
20
      IMPLICIT NONE
21
      INTEGER, INTENT (IN) :: ndm,ng3,jsdm,jspins
22 23 24 25 26
      REAL,    INTENT (IN) :: vl(ndm,jsdm)
      REAL,    INTENT (IN) :: dvx(ndm,jsdm),dvy(ndm,jsdm),dvz(ndm,jsdm)
      REAL, INTENT (IN) :: dvxx(ndm,jsdm),dvyy(ndm,jsdm),dvzz(ndm,jsdm)
      REAL, INTENT (IN) :: dvyz(ndm,jsdm),dvzx(ndm,jsdm),dvxy(ndm,jsdm)

Daniel Wortmann's avatar
Daniel Wortmann committed
27
      TYPE(t_gradients),INTENT(INOUT)::grad
28 29 30 31 32 33 34 35 36 37 38

      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,
     +     dagrzu,dzdx,dzdy,dzdz,
     +     sml
      INTEGER i

      sml = 1.e-14

39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
      IF(ALLOCATED(grad%sigma)) THEN
!Use only contracted gradients for libxc
         if (jspins==1) THEN
            DO i=1,ng3
               grad%sigma(1,i)=
     +             dvx(i,1)*dvx(i,1)+dvy(i,1)*dvy(i,1)+dvz(i,1)*dvz(i,1)
            ENDDO
         ELSE
            DO i=1,ng3
               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)
            ENDDO
         ENDIF            
         RETURN
      ENDIF

59
      DO i = 1,ndm
60 61 62 63 64 65 66 67 68 69
          grad%agrt(i) = 0.0
          grad%agru(i) = 0.0
          grad%agrd(i) = 0.0
          grad%gggrt(i) = 0.0
          grad%gggru(i) = 0.0
          grad%gggrd(i) = 0.0
          grad%gzgr(i) = 0.0
          grad%g2rt(i) = 0.0
          grad%g2ru(i) = 0.0
          grad%g2rd(i) = 0.0
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
      ENDDO

      IF (jspins.eq.1) THEN

        DO 10 i = 1,ng3

          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


          vlt = vlu + vld

          dvxt = dvxu + dvxd
          dvyt = dvyu + dvyd
          dvzt = dvzu + dvzd
          dvxxt = dvxxu + dvxxd
          dvyyt = dvyyu + dvyyd
          dvzzt = dvzzu + dvzzd
          dvyzt = dvyzu + dvyzd
          dvzxt = dvzxu + dvzxd
          dvxyt = dvxyu + dvxyd

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

113 114 115
          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)
116

117 118 119
          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)
120

121 122 123
          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)
124

125 126 127
          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)
128

129 130 131
          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
132 133 134 135 136 137 138 139 140

c         dzdx=d(zeta)/dx,..

          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

c         gzgr=grad(zeta)*grad(ro).

141
          grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
142 143 144

c         g2r: grad(grad(ro))

145 146 147
          grad%g2rt(i) = dvxxt + dvyyt + dvzzt
          grad%g2ru(i) = dvxxu + dvyyu + dvzzu
          grad%g2rd(i) = dvxxd + dvyyd + dvzzd
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


  10    ENDDO

      ELSE

        DO 20 i = 1,ng3

          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)

          vlt = vlu + vld

          dvxt = dvxu + dvxd
          dvyt = dvyu + dvyd
          dvzt = dvzu + dvzd
          dvxxt = dvxxu + dvxxd
          dvyyt = dvyyu + dvyyd
          dvzzt = dvzzu + dvzzd
          dvyzt = dvyzu + dvyzd
          dvzxt = dvzxu + dvzxd
          dvxyt = dvxyu + dvxyd

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

192 193 194
          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)
195

196 197 198
          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)
199

200 201 202
          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)
203

204 205 206
          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)
207

208 209 210
          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
211 212 213 214 215 216 217 218 219

c         dzdx=d(zeta)/dx,..

          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

c         gzgr=grad(zeta)*grad(ro).

220
          grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
221 222 223

c         g2r: grad(grad(ro))

224 225 226
          grad%g2rt(i) = dvxxt + dvyyt + dvzzt
          grad%g2ru(i) = dvxxu + dvyyu + dvzzu
          grad%g2rd(i) = dvxxd + dvyyd + dvzzd
227 228 229 230

  20    ENDDO

      ENDIF
231
    
232 233
      END SUBROUTINE mkgxyz3
      END MODULE m_mkgxyz3