od_mkgz.f 4.14 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

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
      MODULE m_od_mkgz
      CONTAINS
      SUBROUTINE od_mkgz(
     >     z1,nmzxy,delz,
     >     nmzdf,jspins,rh1,rh2,rhdz1,rhdz2,rhdzz1,rhdzz2,
     <     agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,
     <     gzgr)

c     by use of cartesian cylindrical z,phi,r components of charge 
c     density gradients, makes the quantities
c     agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,
c     gzgr
c     used to calculate gradient contribution to xc potential and
c     energy.
c     for more comments on variables go to mkgxyz3_1.F and vvacxcg_1.F
c     Y.Mokrousov, 16-17 Oct 2002

      implicit none
      
      integer, intent (in) :: nmzdf,jspins,nmzxy
      real,    intent (in) :: z1,delz
      real,    intent (in) :: rh1(nmzdf),rh2(nmzdf) 
      real,    intent (in) :: rhdz1(nmzdf),rhdz2(nmzdf)
      real,    intent (in) ::  rhdzz1(nmzdf),rhdzz2(nmzdf)

      real,   intent (out) :: agrt(nmzdf),agru(nmzdf),agrd(nmzdf)
      real,   intent (out) :: g2rt(nmzdf),g2ru(nmzdf),g2rd(nmzdf)
      real,   intent (out) :: gggrt(nmzdf),gggru(nmzdf),gggrd(nmzdf)
      real,   intent (out) :: gzgr(nmzdf)

c     local

      integer :: i
      real    :: vlt,vlu,vld
      real    :: dvzt,dvzzt,dvzu
      real    :: dvzzu,dvzd,dvzzd
      real    :: dagrzt,dagrzu,dagrzd,dztadz
      real    :: sml,z

      sml = 1.e-14

      if (jspins.eq.1) then
         do 10 i = 1,nmzdf

            z = z1 + delz*(nmzxy+i)

            vlu = max(rh1(i)/2,sml)
            dvzu = rhdz1(i)/2
            dvzzu = rhdzz1(i)/2
            vld = vlu
            dvzd = dvzu
            dvzzd = dvzzu

            vlt = vlu+vld
            dvzt = dvzu+dvzd
            dvzzt = dvzzu+dvzzd

c     agr(up,down,total): abs(grad(rho))

            agrt(i) = max(abs(dvzt),sml)
            agru(i) = max(abs(dvzu),sml)
            agrd(i) = max(abs(dvzd),sml)

c     d(abs(grad(rho)))/dr

            dagrzt= dvzt*dvzzt/agrt(i)
            dagrzu= dvzu*dvzzu/agru(i)
            dagrzd= dvzd*dvzzd/agrd(i)

c     grad(rho)*grad(abs(grad(ro)))

            gggrt(i) = dvzt*dagrzt
            gggru(i) = dvzu*dagrzu
            gggrd(i) = dvzd*dagrzd

c     dztadz=d(zeta)/dz,..
            
            dztadz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
            
c     gzgr=grad(zeta)*grad(ro).
            
            gzgr(i) = dztadz*dvzt
            
c     g2rt: grad(grad(ro))
            
            g2rt(i)  = dvzzt + dvzt/z
            g2ru(i)  = dvzzu + dvzu/z
            g2rd(i)  = dvzzd + dvzd/z


 10    continue

      else
         
         do 20 i = 1,nmzdf

            z = z1 + delz*(nmzxy + i)
            
            vlu = max(rh1(i),sml)
            dvzu = rhdz1(i)
            dvzzu = rhdzz1(i)
            vld = max(rh2(i),sml)
            dvzd = rhdz2(i)
            dvzzd = rhdzz2(i)
            
            vlt = vlu+vld
            dvzt = dvzu+dvzd
            dvzzt = dvzzu+dvzzd
            
c     agrt: abs(grad(ro)), u,d for up and down.
            
            agrt(i) = max(abs(dvzt),sml)
            agru(i) = max(abs(dvzu),sml)
            agrd(i) = max(abs(dvzd),sml)
            
            dagrzt= dvzt*dvzzt/agrt(i)
            dagrzu= dvzu*dvzzu/agru(i)
            dagrzd= dvzd*dvzzd/agrd(i)
            
            gggrt(i) = dvzt*dagrzt
            gggru(i) = dvzu*dagrzu
            gggrd(i) = dvzd*dagrzd
            
c     dztadz=d(zeta)/dz,..
            
            dztadz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
            
c     gzgr=grad(zeta)*grad(ro).
            
            gzgr(i) = dztadz*dvzt
            
c     g2rt: grad(grad(ro))
            
            g2rt(i)  = dvzzt + dvzt/z
            g2ru(i)  = dvzzu + dvzu/z
            g2rd(i)  = dvzzd + dvzd/z
            
            
 20      continue

      ENDIF
      
      END SUBROUTINE od_mkgz
      END MODULE m_od_mkgz