LCOV - code coverage report
Current view: top level - xc-pot - xcbh.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 89 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 4 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_xcbh
       8             :    use m_juDFT
       9             : !-----------------------------------------------------------------------
      10             : !     Called in case of icorr=2,3 : spin-polarized exchange-correlation
      11             : !       potential of U. von Barth and L. Hedin, J.Phys.C5,1629 (1972)
      12             : !       icorr = 2: parametrization of Moruzzi,Janak,Williams
      13             : !       icorr = 3: parametrization of von Barth and Hedin
      14             : 
      15             : !     krla=1: Relativistic correction of exchange energy and potential
      16             : !             related to Dirac kinetic energy, according to:
      17             : !             A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
      18             : 
      19             : !     be careful: calculation in rydberg!
      20             : 
      21             : !     vxcbh calculates the XC-potential and
      22             : !     excbh calculates the XC-energy
      23             : 
      24             : !     based on a subroutine by S. Bluegel;   r.pentcheva 22.01.96
      25             : !-----------------------------------------------------------------------
      26             : 
      27             :    USE m_constants, ONLY : pi_const
      28             :    USE m_relcor
      29             :    IMPLICIT NONE
      30             : 
      31             :    REAL, PARAMETER, PRIVATE :: ff  = 3.847322101863  ! 1 / ( 2^(1/3) - 1 )
      32             :    REAL, PARAMETER, PRIVATE :: cvx = 1.221774115422  ! 2 * ( 3/(2*pi) )^(2/3)
      33             :    REAL, PARAMETER, PRIVATE :: cpmjw = 0.045  , cfmjw = 0.0225
      34             :    REAL, PARAMETER, PRIVATE :: cpvbh = 0.0504 , cfvbh = 0.0254
      35             :    REAL, PARAMETER, PRIVATE :: rpmjw = 21.0 , rfmjw = 52.916684096
      36             :    REAL, PARAMETER, PRIVATE :: rpvbh = 30.0 , rfvbh = 75.0
      37             :    REAL, PARAMETER, PRIVATE :: d_15 = 1.e-15
      38             :    REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0
      39             :    REAL, PARAMETER, PRIVATE :: half = 0.5 , thrd = one/three
      40             :    REAL, PARAMETER, PRIVATE :: hfthrd = 0.79370052705 ! 2^(-1/3)
      41             :    REAL, PARAMETER, PRIVATE :: thrhalf = three * half
      42             :    REAL, PARAMETER, PRIVATE :: fothrd = four * thrd , two = 2.0
      43             : 
      44             : CONTAINS
      45             : !************************************************************************
      46           0 :    SUBROUTINE vxcbh &
      47             :       (iofile,xcpot,jspins, &
      48           0 :        mgrid,ngrid,rh, &
      49           0 :        vx,vxc)
      50             : !************************************************************************
      51             :       USE m_types_xcpot_data
      52             : 
      53             : !     .. Scalar Arguments ..
      54             :       INTEGER, INTENT (IN) :: jspins
      55             :       TYPE(t_xcpot_data), INTENT (IN) :: xcpot  !  run mode parameters
      56             :       INTEGER, INTENT (IN) :: iofile      !  file number for read and write
      57             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
      58             : 
      59             : !     .. Array Arguments ..
      60             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      61             :       REAL, INTENT (OUT) :: vxc(mgrid,jspins)     ! xc potential
      62             :       REAL, INTENT (OUT) :: vx (mgrid,jspins)     ! x  potential
      63             : 
      64             : !     .. Local Scalars ..
      65             :       REAL :: txthrd,tythrd,muxp,mucp,mucf,ecfmp,tauc,mucnm
      66             :       REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
      67             :       REAL :: x, y, cp, cf, rp, rf, rs, ecprs, ecfrs
      68             :       INTEGER :: ir
      69             : 
      70             : !     .. Local Arrays ..
      71           0 :       REAL, ALLOCATABLE :: psi(:)       ! relativistic exchange potential corr.
      72             : 
      73             : !---- Intrinsic Functions
      74             :       INTRINSIC alog,max
      75             : 
      76             : !-----> evaluate relativistic corrections for exchange
      77             : 
      78           0 :       ALLOCATE ( psi(ngrid) )
      79             :       CALL relcor( &
      80             :          mgrid,ngrid,jspins,xcpot%krla, .TRUE. ,rh, &
      81           0 :          psi)
      82             : 
      83             : !-----> select exchange correlation potential
      84             : 
      85           0 :       IF (xcpot%is_mjw) THEN
      86             :          cp = cpmjw ; cf = cfmjw
      87             :          rp = rpmjw ; rf = rfmjw
      88           0 :       ELSEIF (xcpot%is_bh) THEN
      89             :          cp = cpvbh ; cf = cfvbh
      90             :          rp = rpvbh ; rf = rfvbh
      91             :       ELSE
      92           0 :          WRITE (iofile,FMT=2000)
      93           0 :          CALL juDFT_error("BUG:vxcbh",calledby="xcbh")
      94             :       END IF
      95             : 2000  FORMAT (13x,'set key for exchange-correlation potential')
      96             : 
      97             : !-----> calculate exchange correlation potential
      98             : 
      99           0 :       IF ( jspins == 2) THEN               ! spinpolarized calculation
     100             : 
     101           0 :          DO ir = 1,ngrid                        ! loop over realspace gridpoints
     102           0 :             rh1 = max(d_15,rh(ir,1))
     103           0 :             rh2 = max(d_15,rh(ir,jspins))
     104           0 :             rho = rh1 + rh2
     105           0 :             x = rh1/rho
     106           0 :             y = rh2/rho
     107           0 :             txthrd = (2*x)**thrd
     108           0 :             tythrd = (2*y)**thrd
     109           0 :             rs= (four*pi_const*rho/three)**thrd
     110           0 :             rs = 1/rs
     111             : 
     112           0 :             ecprs = -cp*fc(rs/rp)            ! calculate correlation energy
     113           0 :             ecfrs = -cf*fc(rs/rf)            ! p : paramagnetic, f : ferromagnetic
     114             :             ! x : exchange,     c : correlation
     115           0 :             muxp = -psi(ir)* (cvx/rs)        ! paramagnetic exchange potential
     116             :             !       (psi contains rel. corr.)
     117           0 :             mucp = -cp*alog(one+rp/rs)       ! calculate correlation potential
     118           0 :             mucf = -cf*alog(one+rf/rs)
     119           0 :             ecfmp = fothrd * (ecfrs-ecprs)
     120           0 :             tauc = mucf - mucp - ecfmp
     121           0 :             mucnm = mucp + tauc*fex(x) - ff*ecfmp
     122             : 
     123           0 :             vxc(ir,1)      = mucnm + (muxp+ff*ecfmp)*txthrd  ! collect correlation
     124           0 :             vxc(ir,jspins) = mucnm + (muxp+ff*ecfmp)*tythrd  ! and exchange parts
     125             : 
     126           0 :             vx (ir,1)      = muxp*txthrd
     127           0 :             vx (ir,jspins) = muxp*tythrd
     128             :          ENDDO
     129             : 
     130           0 :       ELSEIF ( jspins == 1 ) THEN        ! non - spinpolarized calculation
     131             : 
     132           0 :          DO ir = 1, ngrid                   ! loop over realspace gridpoints
     133           0 :             rh1 = max(d_15,rh(ir,1))
     134           0 :             rs = (four*pi_const*rh1/three)**thrd
     135           0 :             rs = 1/rs
     136           0 :             muxp = -psi(ir) * (cvx/rs)       ! paramagnetic exchange potential
     137             :             !       (psi contains rel. corr.)
     138           0 :             mucp = -cp* alog(one+rp/rs)      ! calculate correlation potential
     139           0 :             vxc(ir,1)     = mucp + muxp      ! collect correlation & exchange part
     140             : 
     141           0 :             vx (ir,1)     = muxp
     142             :          ENDDO
     143             : 
     144             :       ELSE
     145           0 :          WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
     146           0 :          CALL juDFT_error("vxcbh",calledby="xcbh")
     147             :       ENDIF
     148             : 
     149           0 :       DEALLOCATE (psi)
     150           0 :       RETURN
     151             : 
     152             :    END SUBROUTINE vxcbh
     153             : !***********************************************************************
     154           0 :    SUBROUTINE excbh &
     155             :       (iofile,xcpot,jspins, &
     156           0 :        mgrid,ngrid,rh, &
     157           0 :        exc)
     158             : !***********************************************************************
     159             :       USE m_types_xcpot_data
     160             : 
     161             : !     .. Scalar Arguments ..
     162             :       INTEGER, INTENT (IN) :: jspins
     163             :       TYPE(t_xcpot_data), INTENT (IN) :: xcpot !  run mode parameters
     164             :       INTEGER, INTENT (IN) :: iofile      !  file number for read and write
     165             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
     166             : 
     167             : !     .. Array Arguments ..
     168             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
     169             :       REAL, INTENT (OUT) :: exc(mgrid)            ! xc energy
     170             : 
     171             : !     .. Local Scalars ..
     172             :       REAL :: thfpi,thrquart,exprs,exfrs,excprs,excfrs
     173             :       REAL :: rho, rh1, rh2 ! total, spin up & spin down charge density
     174             :       REAL :: x, y, cp, cf, rp, rf, rs, ecprs, ecfrs
     175             :       INTEGER :: ir
     176             : 
     177             : !     .. Local Arrays ..
     178           0 :       REAL, ALLOCATABLE :: phi(:)       ! relativistic exchange energy correct.
     179             : 
     180             : !-----> Intrinsic Functions
     181             :       INTRINSIC alog,max
     182             : 
     183           0 :       thrquart = 0.75
     184           0 :       thfpi = thrquart/pi_const
     185             : 
     186           0 :       ALLOCATE ( phi(ngrid) )
     187             :       CALL relcor( &
     188             :          mgrid,ngrid,jspins,xcpot%krla, .FALSE. ,rh, &
     189           0 :          phi)
     190             : 
     191             : !-----> select exchange correlation potential
     192             : 
     193           0 :       IF (xcpot%is_mjw) THEN
     194             :          cp = cpmjw ; cf = cfmjw
     195             :          rp = rpmjw ; rf = rfmjw
     196           0 :       ELSEIF (xcpot%is_bh) THEN
     197             :          cp = cpvbh ; cf = cfvbh
     198             :          rp = rpvbh ; rf = rfvbh
     199             :       ELSE
     200           0 :          WRITE (iofile,FMT=2000)
     201           0 :          CALL juDFT_error("excbh",calledby="xcbh")
     202             :       END IF
     203             : 2000  FORMAT (13x,'set key for exchange-correlation potential')
     204             : 
     205           0 :       IF ( jspins == 2) THEN       ! spinpolarized calculation
     206             : 
     207           0 :          DO  ir = 1,ngrid                  ! loop over realspace gridpoints
     208           0 :             rh1 = max(d_15,rh(ir,1))
     209           0 :             rh2 = max(d_15,rh(ir,jspins))
     210           0 :             rho = rh1 + rh2
     211           0 :             x = rh1/rho
     212           0 :             rs= (thfpi/rho)**thrd
     213             : 
     214           0 :             exprs = -phi(ir)*thrquart*cvx/rs    ! first exchange energy
     215           0 :             exfrs = (2.0**thrd)*exprs           ! phi contains rel. corr.
     216             : 
     217           0 :             ecprs = -cp*fc(rs/rp)               ! calculate correlation energy
     218           0 :             ecfrs = -cf*fc(rs/rf)               ! p: paramagnetic, f: ferromagn.
     219             : 
     220           0 :             excprs = exprs + ecprs              ! now add correlation energy
     221           0 :             excfrs = exfrs + ecfrs
     222             : 
     223           0 :             exc(ir) = excprs + (excfrs-excprs)*fex(x) ! collect all terms
     224             :          ENDDO
     225             : 
     226           0 :       ELSEIF ( jspins == 1 ) THEN  ! non - spinpolarized calculation
     227             : 
     228           0 :          DO ir = 1,ngrid                    ! loop over realspace gridpoints
     229           0 :             rh1 = max(d_15,rh(ir,1))
     230           0 :             rs = (thfpi/rh1)**thrd
     231           0 :             exprs = -phi(ir)*thrquart*cvx/rs ! exchange energy ; phi contains
     232             :             ! relativistic correctionS
     233           0 :             ecprs = -cp*fc(rs/rp)            ! calculate correlation energy
     234           0 :             exc(ir) = exprs + ecprs          ! add correlation energy
     235             :          ENDDO
     236             : 
     237             :       ELSE
     238           0 :          WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
     239           0 :          CALL juDFT_error("excbh",calledby="xcbh")
     240             :       ENDIF
     241             : 
     242           0 :       DEALLOCATE (phi)
     243           0 :       RETURN
     244             : 
     245             :    END SUBROUTINE excbh
     246             : !--------------------------------------------------------------------
     247           0 :    REAL FUNCTION fc(x)
     248             :       REAL :: x
     249             :       fc = (one+(x)*(x)*(x))*alog(one+one/(x)) &
     250           0 :            + half*(x) - (x)*(x) - thrd
     251           0 :    END  FUNCTION fc
     252           0 :    REAL FUNCTION fex(x)
     253             :       REAL :: x
     254           0 :       fex = ff/hfthrd*((x)**fothrd +(1-(x))**fothrd - hfthrd)
     255           0 :    END FUNCTION fex
     256             : !--------------------------------------------------------------------
     257             : 
     258             : END MODULE m_xcbh

Generated by: LCOV version 1.13