xcall.f 7.37 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
      MODULE m_xcall
      use m_juDFT
!***********************************************************************
!
!  vxcall ... driver subroutine for the XC potential and
!  excall ... (only in case of total = .true.)  the XC energy density
!
!  INPUT :  electron density rh()                        (hartree units)
!  OUTPUT:  effective exch-correlation potential vxc()   (hartree units)
!                               energy density   exc()
!
!   icorr selects various exchange correlation potentials:
!   icorr =           -3 means EXX
!                     -2 means HF method
!                      0 means X-alpha method
!                      1 means Wigner correlation
!                      2 means Moruzzi,Janak,Williams correlat
!                      3 means von Barth,Hedin correlation
!                      4 means Vosko,Wilk,Nusair correlation
!                      5 means Perdew,Zunger correlation
!
!   krla=1 : relativistic correction according to
!            A. H. MacDonald and S. H. Vosko, J. Phys. C12, 2977 (1979)
!
!   based on a subroutine from S.Bluegel,  R.Pentcheva, 22.01.96
!
!***********************************************************************
      IMPLICIT NONE
      REAL, PARAMETER, PRIVATE :: hrtr_half = 0.5

      CONTAINS
!***********************************************************************
      SUBROUTINE vxcall(
40
     >                  iofile,xcpot,jspins,
41 42 43 44 45 46 47 48 49
     >                  mgrid,ngrid,rh,
     <                  vx,vxc)
!***********************************************************************

      USE m_xcxal, ONLY : vxcxal
      USE m_xcwgn, ONLY : vxcwgn
      USE m_xcbh,  ONLY : vxcbh
      USE m_xcvwn, ONLY : vxcvwn
      USE m_xcpz,  ONLY : vxcpz
50
      USE m_types
51 52
!
!     .. Scalar Arguments ..
53 54 55
      INTEGER, INTENT (IN) :: iofile ! file number for read and write
      TYPE(t_xcpot),INTENT(IN)::xcpot
      INTEGER, INTENT (IN) :: jspins   ! run mode parameters
56 57 58 59 60 61 62 63 64
      INTEGER, INTENT (IN) :: ngrid,mgrid         ! mesh,number of mesh points
!
!     .. Array Arguments ..
      REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      REAL, INTENT (OUT) :: vx (mgrid,jspins)     ! x potential
      REAL, INTENT (OUT) :: vxc(mgrid,jspins)     ! xc energy  potential
!
!--> Determine exchange correlation energy density
!
65 66
      INTEGER :: ir

67
      IF (xcpot%is_name("x-a"))  THEN   ! X-alpha method
68
         CALL vxcxal(
69
     >               xcpot%krla,jspins,
70 71 72
     >               mgrid,ngrid,rh,
     <               vx,vxc)

73
      ELSEIF (xcpot%is_name("wign")) THEN    ! Wigner interpolation formula
74
         CALL vxcwgn(
75
     >               xcpot%krla,jspins,
76 77
     >               mgrid,ngrid,rh,
     <               vx,vxc)
78
      ELSEIF (xcpot%is_name("mjw").or.xcpot%is_name("bh")) THEN ! von Barth,Hedin correlation
79
        CALL vxcbh(
80
     >             iofile,xcpot,jspins,
81 82 83
     >             mgrid,ngrid,rh,
     <             vx,vxc)

84
      ELSEIF (xcpot%is_name("vwn")) THEN     ! Vosko,Wilk,Nusair correlation
85
        CALL vxcvwn(
86
     >              iofile,xcpot%krla,jspins,
87 88
     >              mgrid,ngrid,rh,
     <              vx,vxc)
89
      ELSEIF (xcpot%is_name("pz")) THEN     ! Perdew,Zunger correlation
90
        CALL vxcpz(                     
91
     >             iofile,xcpot%krla,jspins,
92 93
     >             mgrid,ngrid,rh,
     <             vx,vxc)
94
      ELSEIF (xcpot%is_name("hf")) THEN
95 96 97
      ! Hartree-Fock  calculation: X-alpha potential is added to generate a rational local potential,
      !                            later it is subtracted again
        CALL vxcxal(
98
     >               xcpot%krla,jspins,
99 100 101
     >               mgrid,ngrid,rh,
     <               vx,vxc)
!         vxc=0
102
      ELSEIF (xcpot%is_name("exx")) THEN
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
        ! if exact exchange calculation do nothing
        vxc = 0
      ELSE
          CALL juDFT_error("vxcall",calledby="xcall")
 9000    FORMAT (13x,'set key for exchange-correlation potential',i2)
      ENDIF
!
!---> convert to hartree units
!
      IF (jspins.EQ.2) THEN
         DO ir = 1,ngrid
            vxc(ir,1)      = hrtr_half * vxc(ir,1)
            vxc(ir,jspins) = hrtr_half * vxc(ir,jspins)
         
            vx(ir,1)       = hrtr_half * vx(ir,1)
            vx(ir,jspins)  = hrtr_half * vx(ir,jspins)
         ENDDO
      ELSE IF (jspins.EQ.1) THEN
         DO ir = 1,ngrid
            vxc(ir,1) = hrtr_half * vxc(ir,1)
         
            vx(ir,1)  = hrtr_half * vx(ir,1)
         ENDDO
      ELSE
         WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
          CALL juDFT_error("vxcall",calledby="xcall")
      END IF

      RETURN
      END SUBROUTINE vxcall

!***********************************************************************
      SUBROUTINE  excall(
136
     >                   iofile,xcpot,jspins,
137 138 139 140 141 142 143 144 145
     >                   mgrid,ngrid,rh,
     <                   exc)
!***********************************************************************

      USE m_xcxal, ONLY : excxal
      USE m_xcwgn, ONLY : excwgn
      USE m_xcbh,  ONLY : excbh
      USE m_xcvwn, ONLY : excvwn
      USE m_xcpz,  ONLY : excpz
146
      USE m_types
147 148
!
!     .. Scalar Arguments ..
149 150 151
      type(t_xcpot), INTENT (IN) :: xcpot
      INTEGER, INTENT (IN) :: iofile                    ! file number for read and write
      INTEGER, INTENT (IN) :: jspins   ! run mode parameters
152 153 154 155 156
      INTEGER, INTENT (IN) :: ngrid,mgrid         ! mesh,number of mesh points
!
!     .. Array Arguments ..
      REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      REAL, INTENT (OUT) :: exc(mgrid)            ! xc energy  density
157 158

      INTEGER:: ir
159 160 161
!
!--> Determine exchange correlation energy density
!
162
      IF (xcpot%is_name("x-a"))  THEN   ! X-alpha method
163
         CALL excxal(
164
     >               iofile,xcpot%krla,jspins,
165 166 167
     >               mgrid,ngrid,rh,
     <               exc)
 
168
      ELSEIF (xcpot%is_name("wign")) THEN    ! Wigner interpolation formula
169
         CALL excwgn(
170
     >               iofile,xcpot%krla,jspins,
171 172
     >               mgrid,ngrid,rh,
     <               exc)
173
      ELSEIF (xcpot%is_name("mjw").or.xcpot%is_name("bh")) THEN ! von Barth,Hedin correlation
174
        CALL excbh(
175
     >             iofile,xcpot,jspins,
176 177 178
     >             mgrid,ngrid,rh,
     <             exc)

179
      ELSEIF (xcpot%is_name("vwn")) THEN     ! Vosko,Wilk,Nusair correlation
180
        CALL excvwn(
181
     >              iofile,xcpot%krla,jspins,
182 183
     >              mgrid,ngrid,rh,
     <              exc)
184
      ELSEIF (xcpot%is_name("pz")) THEN     ! Perdew,Zunger correlation
185
         CALL excpz(                      
186
     >              iofile,xcpot%krla,jspins,
187 188
     >              mgrid,ngrid,rh,
     <              exc)
189
      ELSEIF (xcpot%is_name("hf") .or. xcpot%is_name("exx")) THEN
190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
        exc=0
      ELSE
         WRITE (iofile,FMT=9001)
          CALL juDFT_error("excall",calledby="xcall")
 9001    FORMAT (13x,'set key for exchange-correlation potential')
      ENDIF
!
!---> convert to hartree units
!
      DO ir = 1,ngrid
         exc(ir) = hrtr_half * exc(ir)
      ENDDO

      RETURN
      END SUBROUTINE excall
      END MODULE m_xcall