LCOV - code coverage report
Current view: top level - xc-pot - xcpz.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 51 101 50.5 %
Date: 2019-09-08 04:53:50 Functions: 4 6 66.7 %

          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_xcpz
       8             :    use m_juDFT
       9             : !-----------------------------------------------------------------------
      10             : !     Called in case of icorr=5 : spin-polarized exchange-correlation
      11             : !        from Ceperley-Alder monte carlo results with parametrization
      12             : !        of Perdew,Zunger,PRB 23, 5048 (1981)
      13             : 
      14             : !     krla=1: Relativistic correction of exchange energy and potential
      15             : !             related to Dirac kinetic energy, according to:
      16             : !             A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
      17             : 
      18             : !     be careful: calculation in rydberg!
      19             : 
      20             : !     vxcpz   calculates the XC-potential and
      21             : !     excpz   calculates the XC-energy
      22             : 
      23             : !     based on a subroutine by S. Bluegel;   r.pentcheva 22.01.96
      24             : !-----------------------------------------------------------------------
      25             : 
      26             :    USE m_constants, ONLY : pi_const
      27             :    USE m_relcor
      28             :    IMPLICIT NONE
      29             : 
      30             :    REAL, PARAMETER, PRIVATE :: cvx = 1.221774115422  ! 2 * ( 3/(2*pi) )^(2/3)
      31             :    REAL, PARAMETER, PRIVATE :: d_15 = 1.e-15 , c76 = 7.0 / 6.0
      32             :    REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0
      33             :    REAL, PARAMETER, PRIVATE :: half = 0.5 , thrd = one/three
      34             :    REAL, PARAMETER, PRIVATE :: thrhalf = three * half , two = 2.0
      35             :    REAL, PARAMETER, PRIVATE :: c23 = two * thrd , c43 = four * thrd
      36             :    REAL, PARAMETER, PRIVATE :: ap  =  0.0622 , af  =  0.0311
      37             :    REAL, PARAMETER, PRIVATE :: bp  = -0.0960 , bf  = -0.0538
      38             :    REAL, PARAMETER, PRIVATE :: cp  =  0.0040 , cf  =  0.0014
      39             :    REAL, PARAMETER, PRIVATE :: dp  = -0.0232 , df  = -0.0096
      40             :    REAL, PARAMETER, PRIVATE :: gp  = -0.2846 , gf  = -0.1686
      41             :    REAL, PARAMETER, PRIVATE :: b1p =  1.0529 , b1f =  1.3981
      42             :    REAL, PARAMETER, PRIVATE :: b2p =  0.3334 , b2f =  0.2611
      43             : 
      44             : CONTAINS
      45             : !************************************************************************
      46         222 :    SUBROUTINE vxcpz( &
      47             :       iofile,krla,jspins, &
      48         222 :       mgrid,ngrid,rh, &
      49         222 :       vx,vxc)
      50             : !************************************************************************
      51             : 
      52             : !     .. Scalar Arguments ..
      53             :       INTEGER, INTENT (IN) :: jspins
      54             :       INTEGER, INTENT (IN) :: krla        !  run mode parameters
      55             :       INTEGER, INTENT (IN) :: iofile      !  file number for read and write
      56             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
      57             : 
      58             : !     .. Array Arguments ..
      59             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)                   ! charge density
      60             :       REAL, INTENT (OUT) :: vx(mgrid,jspins),vxc(mgrid,jspins) ! x/xc potential
      61             : 
      62             : !     .. Local Scalars ..
      63             :       REAL :: c_2, cbrt1, cbrt2, dfds, decds, dec1, dec2, vcf, vcp
      64             : 
      65             : !     .. Local Arrays ..
      66         222 :       REAL, ALLOCATABLE :: psi(:)       ! relativistic exchange potential corr.
      67             : 
      68             : !-----s Intrinsic Functions
      69             :       INTRINSIC max
      70             :       REAL :: rho, rh1, rh2 ! total, spin up & spin down  charge density
      71             :       REAL :: fothrd, thfpi, c_1, y1, y2, s, fs, rs
      72             :       REAL :: ecp, ecf
      73             :       INTEGER :: ir
      74             : 
      75         222 :       fothrd = c43
      76         222 :       thfpi  = three / ( four * pi_const )
      77             : 
      78             : !-----> evaluate relativistic corrections for exchange
      79             : 
      80         222 :       ALLOCATE ( psi(ngrid) )
      81             :       CALL relcor( &
      82             :          mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
      83         222 :          psi)
      84             : 
      85         222 :       IF ( jspins == 2 ) THEN         ! spinpolarized calculation
      86             : 
      87             :          c_1 = one / ( two**fothrd - two )
      88           0 :          DO ir = 1,ngrid
      89           0 :             rh1 = max(d_15,rh(ir,1))
      90           0 :             rh2 = max(d_15,rh(ir,jspins))
      91           0 :             rho = rh1 + rh2
      92           0 :             y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1  ! s = (rh2 - rh1) / rho
      93           0 :             cbrt1 = (one-s) ** thrd
      94           0 :             cbrt2 = (one+s) ** thrd
      95           0 :             fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
      96           0 :             dfds = c_1 *  fothrd * (cbrt1 - cbrt2)
      97           0 :             rs = ( thfpi/rho )**thrd
      98             : 
      99           0 :             IF (rs >= one) THEN
     100           0 :                ecp = fecl(rs,gp,b1p,b2p)  ! correlation energy paramagnetic
     101           0 :                ecf = fecl(rs,gf,b1f,b2f)  ! correlation energy ferromagnetic
     102           0 :                vcp = fvcl(ecp,rs,b1p,b2p) ! d(rho*ecp)/d(rho) = ecp - rs/3*d(ecp)/d(rs)
     103           0 :                vcf = fvcl(ecf,rs,b1f,b2f) ! d(rho*ecf)/d(rho)
     104             :             ELSE
     105           0 :                ecp = fecs(rs,ap,bp,cp,dp)
     106           0 :                ecf = fecs(rs,af,bf,cf,df)
     107           0 :                vcp = fvcs(rs,ap,bp,cp,dp)
     108           0 :                vcf = fvcs(rs,af,bf,cf,df)
     109             :             ENDIF
     110           0 :             decds = (ecf-ecp)*dfds     ! = d(ec)/d(rho)
     111           0 :             dec1  = vcp + (vcf-vcp)*fs ! = d(rho*ec)/d(rho)
     112           0 :             dec2  = two/rho*decds      ! = 2/rho*d(ec)/ds
     113             : 
     114           0 :             c_2 = cvx / rs * psi(ir)                   ! exchange potential muxp=-cvx/rs= 4/3*ex
     115           0 :             vxc(ir,1)     =dec1+dec2*rh2 - c_2*cbrt1   !                        muxp*(2x)**(1/3)
     116           0 :             vxc(ir,jspins)=dec1-dec2*rh1 - c_2*cbrt2   ! calculate exchange correlation potential
     117             :             ! vc = ec +vcp + (vcf-vcp)f(s) - (ecf-ecp)df/ds*(s+/-1)
     118             : 
     119           0 :             vx (ir,1)     = - c_2*cbrt1
     120           0 :             vx (ir,jspins)= - c_2*cbrt2
     121             :          ENDDO
     122             : 
     123         222 :       ELSEIF ( jspins == 1 ) THEN     ! non-spinpolarized calculation
     124             : 
     125    11749322 :          DO ir = 1,ngrid
     126     5874550 :             rho = max(d_15,rh(ir,1))
     127     5874550 :             rs = ( thfpi/rho )**thrd
     128     5874550 :             IF (rs >= one) THEN
     129      674212 :                ecp = fecl(rs,gp,b1p,b2p)
     130      674212 :                vcp = fvcl(ecp,rs,b1p,b2p)
     131             :             ELSE
     132     5200338 :                ecp = fecs(rs,ap,bp,cp,dp)
     133     5200338 :                vcp = fvcs(rs,ap,bp,cp,dp)
     134             :             ENDIF
     135     5874550 :             vxc(ir,1) = vcp - cvx / rs * psi(ir)
     136             : 
     137     5874772 :             vx (ir,1) = cvx / rs * psi(ir)
     138             :          ENDDO
     139             : 
     140             :       ELSE
     141           0 :          WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
     142           0 :          CALL juDFT_error("vxcpz",calledby="xcpz")
     143             :       ENDIF
     144             : 
     145         222 :       DEALLOCATE (psi)
     146         222 :       RETURN
     147             : 
     148             :    END SUBROUTINE vxcpz
     149             : !***********************************************************************
     150          24 :    SUBROUTINE excpz( &
     151             :       iofile,krla,jspins, &
     152          24 :       mgrid,ngrid,rh, &
     153          24 :       exc)
     154             : !***********************************************************************
     155             : 
     156             : !     .. Scalar Arguments ..
     157             :       INTEGER, INTENT (IN) :: jspins
     158             :       INTEGER, INTENT (IN) :: krla        !  run mode parameters
     159             :       INTEGER, INTENT (IN) :: iofile      !  file number for read and write
     160             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
     161             : 
     162             : !     .. Array Arguments ..
     163             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
     164             :       REAL, INTENT (OUT) :: exc(mgrid)            ! xc energy
     165             : 
     166             : !     .. Local Scalars ..
     167             :       REAL :: ec, ex, cex
     168             : 
     169             : !     .. Local Arrays ..
     170          24 :       REAL, ALLOCATABLE :: phi(:)       ! relativistic exchange energy correct.
     171             : 
     172             : !-----> Intrinsic Functions
     173             :       INTRINSIC max
     174             : 
     175             :       REAL :: rho, rh1, rh2 ! total, spin up & spin down  charge density
     176             :       REAL :: fothrd, thfpi, c_1, y1, y2, s, fs, rs
     177             :       REAL :: ecp, ecf
     178             :       INTEGER :: ir
     179             : 
     180          24 :       fothrd = c43
     181          24 :       thfpi  = three / ( four * pi_const )
     182          24 :       cex = cvx / c43
     183             : 
     184          24 :       ALLOCATE ( phi(ngrid) )
     185             :       CALL relcor( &
     186             :          mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
     187          24 :          phi)
     188             : 
     189          24 :       IF ( jspins == 2 ) THEN         ! spinpolarized calculation
     190             : 
     191             :          c_1 = one / ( two**fothrd - two )
     192           0 :          DO ir = 1,ngrid
     193           0 :             rh1 = max(d_15,rh(ir,1))
     194           0 :             rh2 = max(d_15,rh(ir,jspins))
     195           0 :             rho = rh1 + rh2
     196           0 :             rs = ( thfpi/rho )**thrd
     197           0 :             y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1  ! s = (rh2 - rh1) / rho
     198           0 :             fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
     199           0 :             IF (rs >= one) THEN
     200           0 :                ecp = fecl(rs,gp,b1p,b2p)  ! correlation energy paramagnetic
     201           0 :                ecf = fecl(rs,gf,b1f,b2f)  ! correlation energy ferromagnetic
     202             :             ELSE
     203           0 :                ecp = fecs(rs,ap,bp,cp,dp)
     204           0 :                ecf = fecs(rs,af,bf,cf,df)
     205             :             ENDIF
     206             : 
     207           0 :             ec = ecp + (ecf-ecp)*fs                ! total correlation energy
     208           0 :             ex = -cex/rs* (one + 0.2599210482*fs)  ! ex = exp + (exf-exp)*f(s)
     209             :             ! exf-exp = (2**(1/3)-1) * exp
     210           0 :             exc(ir) = ec + ex*phi(ir)
     211             :          ENDDO
     212             : 
     213          24 :       ELSEIF ( jspins == 1 ) THEN     ! non-spinpolarized calculation
     214             : 
     215    11308376 :          DO ir = 1,ngrid
     216     5654176 :             rho = max(d_15,rh(ir,1))
     217     5654176 :             rs = ( thfpi/rho )**thrd
     218     5654176 :             IF (rs >= one) THEN
     219      630976 :                ecp = fecl(rs,gp,b1p,b2p)
     220             :             ELSE
     221     5023200 :                ecp = fecs(rs,ap,bp,cp,dp)
     222             :             ENDIF
     223     5654176 :             ex = -cex/rs
     224     5654200 :             exc(ir) = ecp + ex*phi(ir)
     225             :          ENDDO
     226             : 
     227             :       ELSE
     228           0 :          WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
     229           0 :          CALL juDFT_error("vxcpz",calledby="xcpz")
     230             :       ENDIF
     231             : 
     232          24 :       DEALLOCATE (phi)
     233          24 :       RETURN
     234             : 
     235             :    END SUBROUTINE excpz
     236             : 
     237             : !--------------------------------------------------------------------
     238           0 :    REAL FUNCTION fecl(r,g,b1,b2)
     239             :       REAL :: r,g,b1,b2
     240     1305188 :       fecl = g / ( one + b1*sqrt(r) + b2*r )
     241           0 :    END  FUNCTION fecl
     242      674212 :    REAL FUNCTION fvcl(ce,r,b1,b2)
     243             :       REAL :: ce,r,b1,b2
     244      674212 :       fvcl = ce* (one+c76*b1*sqrt(r)+c43*b2*r)/(one+b1*sqrt(r)+b2*r)
     245      674212 :    END FUNCTION fvcl
     246           0 :    REAL FUNCTION fecs(r,a,b,c,d)
     247             :       REAL :: r,a,b,c,d
     248             :       INTRINSIC alog
     249    10223538 :       fecs = a*alog(r) + b + c*r*alog(r) + d*r
     250           0 :    END FUNCTION fecs
     251     5200338 :    REAL FUNCTION fvcs(r,a,b,c,d)
     252             :       REAL :: r,a,b,c,d
     253             :       INTRINSIC alog
     254             :       fvcs = a*alog(r) + (b-a/three) + c23*c*r*alog(r) + &
     255     5200338 :              (two*d-c)*r/three
     256     5200338 :    END FUNCTION fvcs
     257             : !--------------------------------------------------------------------
     258             : 
     259             : END MODULE m_xcpz

Generated by: LCOV version 1.13