mkgz.f 3.39 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
      MODULE m_mkgz
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,
cc      gzgr
cc    used to calculate gradient contribution to xc potential and
cc    energy.
c.....------------------------------------------------------------------
      CONTAINS 
      SUBROUTINE mkgz(
12 13 14
     >     nmzdf,jspins,rh1,rh2,rhdz1,rhdz2,rhdzz1,rhdzz2,
     <     grad)
      USE m_types
15 16 17 18
      IMPLICIT NONE
      INTEGER, INTENT (IN) :: nmzdf,jspins
      REAL,    INTENT (IN) :: rh1(nmzdf),rhdz1(nmzdf),rhdzz1(nmzdf)
      REAL,    INTENT (IN) :: rh2(nmzdf),rhdz2(nmzdf),rhdzz2(nmzdf)
Daniel Wortmann's avatar
Daniel Wortmann committed
19
      TYPE(t_gradients),INTENT(INOUT)::grad
20 21 22 23 24 25 26

      INTEGER i
      REAL    vlt,dvzt,dvzzt,vlu,dvzu,dvzzu,vld,dvzd,dvzzd
      REAL    dagrzt,dagrzu,dagrzd,dztadz,sml

      sml = 1.e-14

27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
      IF (ALLOCATED(grad%sigma)) THEN
         IF(jspins==1) THEN
            DO i=1,nmzdf
               grad%sigma(1,i)=rhdz1(i)*rhdz1(i)
            ENDDO
         ELSE
             DO i=1,nmzdf
               grad%sigma(1,i)=rhdz1(i)*rhdz1(i)
               grad%sigma(2,i)=rhdz1(i)*rhdz2(i)
               grad%sigma(3,i)=rhdz2(i)*rhdz2(i)
            ENDDO
         ENDIF
         RETURN
      ENDIF
      
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
      IF (jspins == 1) THEN

        DO i = 1,nmzdf

          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         agrt: abs(grad(ro)), u,d for up and down.

59 60 61
          grad%agrt(i) = max(abs(dvzt),sml)
          grad%agru(i) = max(abs(dvzu),sml)
          grad%agrd(i) = max(abs(dvzd),sml)
62

63 64 65
          dagrzt= dvzt*dvzzt/grad%agrt(i)
          dagrzu= dvzu*dvzzu/grad%agru(i)
          dagrzd= dvzd*dvzzd/grad%agrd(i)
66

67 68 69
          grad%gggrt(i) = dvzt*dagrzt
          grad%gggru(i) = dvzu*dagrzu
          grad%gggrd(i) = dvzd*dagrzd
70 71 72 73 74 75 76

c         dztadz=d(zeta)/dz,..

          dztadz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2

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

77
          grad%gzgr(i) = dztadz*dvzt
78 79 80

c         g2rt: grad(grad(ro))

81 82 83
          grad%g2rt(i)  = dvzzt
          grad%g2ru(i)  = dvzzu
          grad%g2rd(i)  = dvzzd
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103

        ENDDO

      ELSE

        DO i = 1,nmzdf

          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.

104 105 106
          grad%agrt(i) = max(abs(dvzt),sml)
          grad%agru(i) = max(abs(dvzu),sml)
          grad%agrd(i) = max(abs(dvzd),sml)
107

108 109 110
          dagrzt= dvzt*dvzzt/grad%agrt(i)
          dagrzu= dvzu*dvzzu/grad%agru(i)
          dagrzd= dvzd*dvzzd/grad%agrd(i)
111

112 113 114
          grad%gggrt(i) = dvzt*dagrzt
          grad%gggru(i) = dvzu*dagrzu
          grad%gggrd(i) = dvzd*dagrzd
115 116 117 118 119 120 121

c         dztadz=d(zeta)/dz,..

          dztadz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2

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

122
          grad%gzgr(i) = dztadz*dvzt
123 124 125

c         g2rt: grad(grad(ro))

126 127 128
          grad%g2rt(i)  = dvzzt
          grad%g2ru(i)  = dvzzu
          grad%g2rd(i)  = dvzzd
129 130 131 132 133 134 135


        ENDDO
      ENDIF

      END SUBROUTINE mkgz
      END MODULE m_mkgz