LCOV - code coverage report
Current view: top level - xc-pot - xcvwn.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 110 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_xcvwn
       8             :    use m_juDFT
       9             : !-----------------------------------------------------------------------
      10             : !     Called in case of icorr=4 : spin-polarized exchange-correlation
      11             : !       potential from Ceperley-Alder monte carlo results with parametrization
      12             : !       of Vosko,Wilk,Nusair, Can. J. Phys. 58, 1200 (1980)
      13             : !         ( In case of non-spinpolarization the fit of the correlation
      14             : !           energy and potential is essentially equivallent to the
      15             : !           parametrization of Perdew,Zunger,PRB 23, 5048 (1981), but
      16             : !           in case of spin polarization the magnetization dependent
      17             : !           interpolation formula for correlation is better.
      18             : !           Part of the subroutine was supplied by M. Manninen. )
      19             : 
      20             : !     krla=1: Relativistic correction of exchange energy and potential
      21             : !             related to Dirac kinetic energy, according to:
      22             : !             A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
      23             : 
      24             : !     be careful: calculation in rydberg!
      25             : 
      26             : !     vxc   calculates the XC-potential and
      27             : !     exc   calculates the XC-energy
      28             : 
      29             : !     based on a subroutine by S. Bluegel;   r.pentcheva 22.01.96
      30             : !-----------------------------------------------------------------------
      31             : 
      32             :    USE m_constants, ONLY : pi_const
      33             :    USE m_relcor
      34             :    IMPLICIT NONE
      35             : 
      36             :    REAL, PARAMETER, PRIVATE :: cex = 0.91633058742  ! 3/2 * ( 3/(2*pi) )^(2/3)
      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 :: thrd = one/three , two = 2.0
      40             :    REAL, PARAMETER, PRIVATE :: fothrd = four * thrd
      41             :    REAL, PARAMETER, PRIVATE :: ap = 0.0621814 , xp0 = -0.10498
      42             :    REAL, PARAMETER, PRIVATE :: af = 0.0310907 , xf0 = -0.32500
      43             :    REAL, PARAMETER, PRIVATE :: al =-0.03377373, xl0 = -0.0047584
      44             :    REAL, PARAMETER, PRIVATE :: bp = 3.72744   , cp  = 12.9352
      45             :    REAL, PARAMETER, PRIVATE :: bf = 7.06042   , cf  = 18.0578
      46             :    REAL, PARAMETER, PRIVATE :: bl = 1.13107   , cl  = 13.0045
      47             :    REAL, PARAMETER, PRIVATE :: qp = 6.1519908 , fdd0 = 1.70992093
      48             :    REAL, PARAMETER, PRIVATE :: qf = 4.7309269 , ql = 7.123109
      49             : 
      50             : CONTAINS
      51             : !************************************************************************
      52           0 :    SUBROUTINE vxcvwn( &
      53             :       iofile,krla,jspins, &
      54           0 :       mgrid,ngrid,rh, &
      55           0 :       vx,vxc)
      56             : !************************************************************************
      57             : 
      58             : !     .. Scalar Arguments ..
      59             :       INTEGER, INTENT (IN) :: jspins
      60             :       INTEGER, INTENT (IN) :: krla        !  run mode parameters
      61             :       INTEGER, INTENT (IN) :: iofile      !  file number for read and write
      62             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
      63             : 
      64             : !     .. Array Arguments ..
      65             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      66             :       REAL, INTENT (OUT) :: vxc(mgrid,jspins)     ! xc potential
      67             :       REAL, INTENT (OUT) :: vx (mgrid,jspins)     ! x  potential
      68             : 
      69             : !     .. Local Scalars ..
      70             :       REAL :: s3, s4, c_2, cbrt1, cbrt2, dfds, decdrp, decdrf
      71             :       REAL :: dacdr, dbdr, dec1, dec2, cvx, mucp
      72             :       REAL :: rho, rh1, rh2 ! total, spin up & spin down  charge density
      73             :       REAL :: x, y1, y2, s, thfpi, c_1, rs, beta, bs41, alc
      74             :       REAL :: xpx, xfx, xlx, xpx0, xfx0, xlx0, ecp, ecf, fs, ec
      75             :       INTEGER :: ir
      76             : 
      77             : !     .. Local Arrays ..
      78           0 :       REAL, ALLOCATABLE :: psi(:)       ! relativistic exchange potential corr.
      79             : 
      80           0 :       thfpi  = three / ( four * pi_const )
      81             : 
      82             : !-----> evaluate relativistic corrections for exchange
      83             : 
      84           0 :       ALLOCATE ( psi(ngrid) )
      85             :       CALL relcor( &
      86             :          mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
      87           0 :          psi)
      88             : 
      89           0 :       cvx = fothrd * cex
      90           0 :       IF ( jspins == 2 ) THEN         ! spinpolarized calculation
      91             : 
      92             :          c_1 = one / ( two**fothrd - two )
      93           0 :          DO ir = 1,ngrid
      94           0 :             rh1 = max(d_15,rh(ir,1))
      95           0 :             rh2 = max(d_15,rh(ir,jspins))
      96           0 :             rho = rh1 + rh2
      97           0 :             y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1  ! s = (rh2 - rh1) / rho
      98           0 :             s3 = s**3 ; s4 = s*s3                      ! spin polarisation
      99           0 :             cbrt1 = (one-s) ** thrd
     100           0 :             cbrt2 = (one+s) ** thrd
     101           0 :             fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
     102           0 :             dfds = c_1 * fothrd * ( cbrt1-cbrt2 )
     103           0 :             rs = ( thfpi/rho )**thrd
     104           0 :             x = sqrt(rs)
     105           0 :             xpx = x*x + bp*x + cp
     106           0 :             xfx = x*x + bf*x + cf
     107           0 :             xlx = x*x + bl*x + cl
     108           0 :             xpx0 = xp0*xp0 + bp*xp0 + cp
     109           0 :             xfx0 = xf0*xf0 + bf*xf0 + cf
     110           0 :             xlx0 = xl0*xl0 + bl*xl0 + cl
     111           0 :             ecp  = fec(ap,x,xpx,xp0,xpx0,bp,qp)   ! paramagnetic correlation energy
     112           0 :             ecf  = fec(af,x,xfx,xf0,xfx0,bf,qf)   ! ferromagnetic correlation ener.
     113           0 :             alc  = fec(al,x,xlx,xl0,xlx0,bl,ql)   ! alpha_c
     114           0 :             beta = fdd0* (ecf-ecp)/alc - one      ! beta = (f"(0)*(ecf-ecp))/alc -1
     115           0 :             bs41 = one + beta * s**4
     116             : 
     117           0 :             ec = ecp + alc * fs / fdd0 * bs41        ! total correlation energy
     118           0 :             decdrp = fdedr(rho,ap,x,xpx,xp0,bp,cp)   ! d(ecp)/d(rho)
     119           0 :             decdrf = fdedr(rho,af,x,xfx,xf0,bf,cf)   ! d(ecf)/d(rho)
     120           0 :             dacdr  = fdedr(rho,al,x,xlx,xl0,bl,cl)   ! d(alc)/d(rho)
     121             : 
     122           0 :             dbdr =fdd0* ((decdrf-decdrp)*alc- (ecf-ecp)*dacdr)/alc**2     ! d(beta)/d(rho)
     123           0 :             dec1 =ec + rho* (decdrp+ (dacdr*fs*bs41+alc*fs*dbdr*s4)/fdd0) ! mucp= d(rho*ec)/d(rho)
     124           0 :             dec2 =two*alc/ (fdd0*rho)* (dfds*bs41+four*fs*beta*s3)        ! 2/rho*d(ec)/ds
     125             : 
     126           0 :             c_2 = cvx / rs * psi(ir)                   ! exchange potential muxp=-cvx/rs= 4/3*ex
     127           0 :             vxc(ir,1)     =dec1+dec2*rh2 - c_2*cbrt1   !                        muxp*(2x)**(1/3)
     128           0 :             vxc(ir,jspins)=dec1-dec2*rh1 - c_2*cbrt2   ! calculate exchange correlation potential
     129             : 
     130           0 :             vx (ir,1)     = - c_2*cbrt1
     131           0 :             vx (ir,jspins)= - c_2*cbrt2
     132             :          ENDDO
     133             : 
     134           0 :       ELSEIF ( jspins == 1 ) THEN     ! non-spinpolarized calculation
     135             : 
     136           0 :          DO ir = 1,ngrid
     137           0 :             rho = max(d_15,rh(ir,1))
     138           0 :             rs = ( thfpi/rho )**thrd
     139           0 :             x = sqrt(rs)
     140           0 :             xpx = x*x + bp*x + cp
     141           0 :             xpx0 = xp0*xp0 + bp*xp0 + cp
     142             : 
     143           0 :             ecp = fec(ap,x,xpx,xp0,xpx0,bp,qp)     ! correlation energy paramagn.
     144           0 :             decdrp = fdedr(rho,ap,x,xpx,xp0,bp,cp) ! d(ecp)/d(rho)
     145           0 :             mucp = ecp + rho*decdrp                ! d(rho*ec)/d(rho)
     146             : 
     147           0 :             vxc(ir,1)     = mucp - cvx/rs*psi(ir)  ! -1.221774/rs :exchange potential
     148             : 
     149           0 :             vx (ir,1)     = - cvx/rs*psi(ir)
     150             :          ENDDO
     151             : 
     152             :       ELSE
     153           0 :          WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
     154           0 :          CALL juDFT_error("vxcvwn",calledby="xcvwn")
     155             :       ENDIF
     156             : 
     157           0 :       DEALLOCATE (psi)
     158           0 :       RETURN
     159             : 
     160             :    END SUBROUTINE vxcvwn
     161             : !***********************************************************************
     162           0 :    SUBROUTINE excvwn( &
     163             :       iofile,krla,jspins, &
     164           0 :       mgrid,ngrid,rh, &
     165           0 :       exc)
     166             : !***********************************************************************
     167             : 
     168             : !     .. Scalar Arguments ..
     169             :       INTEGER, INTENT (IN) :: jspins
     170             :       INTEGER, INTENT (IN) :: krla        !  run mode parameters
     171             :       INTEGER, INTENT (IN) :: iofile      !  file number for read and write
     172             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
     173             : 
     174             : !     .. Array Arguments ..
     175             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
     176             :       REAL, INTENT (OUT) :: exc(mgrid)            ! xc energy
     177             : 
     178             : !     .. Local Scalars ..
     179             :       REAL ::  c_1,ex
     180             :       REAL :: rho, rh1, rh2     ! total, spin up & spin down  charge density
     181             :       REAL :: x, y1, y2, s, thfpi, rs, beta, bs41, alc
     182             :       REAL :: xpx, xfx, xlx, xpx0, xfx0, xlx0, ecp, ecf, fs, ec
     183             :       INTEGER :: ir
     184             : 
     185             : !     .. Local Arrays ..
     186           0 :       REAL, ALLOCATABLE :: phi(:)       ! relativistic exchange energy correct.
     187             : 
     188           0 :       thfpi  = three / ( four * pi_const )
     189             : 
     190           0 :       ALLOCATE ( phi(ngrid) )
     191             :       CALL relcor( &
     192             :          mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
     193           0 :          phi)
     194             : 
     195           0 :       IF ( jspins == 2 ) THEN         ! spinpolarized calculation
     196             : 
     197             :          c_1 = one / ( two**fothrd - two )
     198           0 :          DO ir = 1,ngrid
     199           0 :             rh1 = max(d_15,rh(ir,1))
     200           0 :             rh2 = max(d_15,rh(ir,jspins))
     201           0 :             rho = rh1 + rh2
     202           0 :             y1 = rh1/rho ; y2 = rh2/rho ; s = y2 - y1
     203           0 :             fs = c_1 * ( (one+s)**fothrd + (one-s)**fothrd - two )
     204           0 :             rs = ( thfpi/rho )**thrd
     205           0 :             x = sqrt(rs)
     206           0 :             xpx = x*x + bp*x + cp
     207           0 :             xfx = x*x + bf*x + cf
     208           0 :             xlx = x*x + bl*x + cl
     209           0 :             xpx0 = xp0*xp0 + bp*xp0 + cp
     210           0 :             xfx0 = xf0*xf0 + bf*xf0 + cf
     211           0 :             xlx0 = xl0*xl0 + bl*xl0 + cl
     212           0 :             ecp  = fec(ap,x,xpx,xp0,xpx0,bp,qp)   ! paramagnetic correlation energy
     213           0 :             ecf  = fec(af,x,xfx,xf0,xfx0,bf,qf)   ! ferromagnetic correlation ener.
     214           0 :             alc  = fec(al,x,xlx,xl0,xlx0,bl,ql)   ! alpha_c
     215           0 :             beta = fdd0* (ecf-ecp)/alc - one      ! beta = (f"(0)*(ecf-ecp))/alc -1
     216           0 :             bs41 = one + beta * s**4
     217             : 
     218           0 :             ec = ecp + alc * fs / fdd0 * bs41        ! total correlation energy
     219           0 :             ex = -cex/rs* (one + 0.2599210482 * fs)  ! ex = exp + (exf-exp)*f(s)
     220             :             ! exf - exp = (2**(1/3)-1) * exp
     221             :             !             = 0.25992105 * exp
     222           0 :             exc(ir) = ec + ex*phi(ir)
     223             :          ENDDO
     224             : 
     225           0 :       ELSEIF ( jspins == 1 ) THEN     ! non-spinpolarized calculation
     226             : 
     227           0 :          DO ir = 1,ngrid
     228           0 :             rho = max(d_15,rh(ir,1))
     229           0 :             rs = ( thfpi/rho )**thrd
     230           0 :             x = sqrt(rs)
     231           0 :             xpx = x*x + bp*x + cp
     232           0 :             xpx0 = xp0*xp0 + bp*xp0 + cp
     233           0 :             ecp = fec(ap,x,xpx,xp0,xpx0,bp,qp)  ! correlation energy paramagn.
     234           0 :             ex  = -cex / rs                     ! exp: paramag. exchange energy
     235             :             ! (like in wigner formula)
     236           0 :             exc(ir) = ecp + ex*phi(ir)
     237             :          ENDDO
     238             : 
     239             :       ELSE
     240           0 :          WRITE (iofile,'('' error in jspins, jspins ='',i2)') jspins
     241           0 :          CALL juDFT_error("excvwn",calledby="xcvwn")
     242             :       ENDIF
     243             : 
     244           0 :       DEALLOCATE (phi)
     245           0 :       RETURN
     246             : 
     247             :    END SUBROUTINE excvwn
     248             : 
     249             : !--------------------------------------------------------------------
     250           0 :    REAL FUNCTION fec(ai,z,ziz,zi0,ziz0,bi,qi)
     251             :       REAL :: ai,z,ziz,zi0,ziz0,bi,qi
     252             :       fec = ai* (alog(z*z/ziz)+ two*bi/qi * atan(qi/(two*z+bi)) - &
     253             :                  bi*zi0/ziz0* ( alog((z-zi0)**2/ziz) + &
     254           0 :                                two*(bi+ two*zi0)/qi*atan(qi/ (two*z+bi))) )
     255           0 :    END  FUNCTION fec
     256           0 :    REAL FUNCTION fdedr(rho,ai,z,ziz,zi0,bi,ci)
     257             :       REAL :: rho,ai,z,ziz,zi0,bi,ci
     258           0 :       fdedr = -ai*z/(three*rho*ziz)*(ci/z-bi*zi0/(z-zi0))
     259           0 :    END FUNCTION fdedr
     260             : !--------------------------------------------------------------------
     261             : 
     262             : END MODULE m_xcvwn

Generated by: LCOV version 1.13