LCOV - code coverage report
Current view: top level - xc-pot - xcwgn.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 44 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 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_xcwgn
       8             : !-----------------------------------------------------------------------
       9             : !     Called in case of icorr=1 : non-spinpolarized exchange-correlation
      10             : !                                     from wigners interpolation formula.
      11             : 
      12             : !     krla=1: Relativistic correction of exchange energy and potential
      13             : !             related to Dirac kinetic energy, according to:
      14             : !             A.H. MacDonald and S.H. Vosko, J. Phys. C12, 2977 (1979)
      15             : 
      16             : !     be careful: calculation in rydberg!
      17             : 
      18             : !     vxc   calculates the XC-potential and
      19             : !     exc   calculates the XC-energy
      20             : 
      21             : !     based on a subroutine by S. Bluegel;   r.pentcheva 22.01.96
      22             : !-----------------------------------------------------------------------
      23             : 
      24             :    USE m_constants, ONLY : pi_const
      25             :    USE m_relcor
      26             :    IMPLICIT NONE
      27             : 
      28             :    REAL, PARAMETER, PRIVATE :: one = 1.0 , three = 3.0 , four = 4.0
      29             :    REAL, PARAMETER, PRIVATE :: thrd = one/three , d_15 = 1.e-15
      30             :    REAL, PARAMETER, PRIVATE :: cex = 0.91633058742  ! 3/2 * ( 3/(2*pi) )^(2/3)
      31             : 
      32             : CONTAINS
      33             : !************************************************************************
      34           0 :    SUBROUTINE vxcwgn( &
      35             :       krla,jspins, &
      36           0 :       mgrid,ngrid,rh, &
      37           0 :       vx,vxc)
      38             : !************************************************************************
      39             : 
      40             : !     .. Scalar Arguments ..
      41             :       INTEGER, INTENT (IN) :: jspins
      42             :       INTEGER, INTENT (IN) :: krla        !  run mode parameters
      43             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
      44             : 
      45             : !     .. Array Arguments ..
      46             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
      47             :       REAL, INTENT (OUT) :: vx(mgrid,jspins)      ! x  potential
      48             :       REAL, INTENT (OUT) :: vxc(mgrid,jspins)     ! xc potential
      49             : 
      50             : !     .. Local Scalars ..
      51             :       REAL :: fothrd,vxp,vcp
      52             :       REAL ::  rs, rho, thfpi, exp, ecp
      53             :       INTEGER :: ir
      54             : 
      55             : !     .. Local Arrays ..
      56           0 :       REAL, ALLOCATABLE :: psi(:)       ! relativistic exchange potential corr.
      57             : 
      58             : !-----s Intrinsic Functions
      59             :       INTRINSIC alog,max
      60             : 
      61           0 :       thfpi  = three / ( four * pi_const )
      62           0 :       fothrd = four/three
      63             : 
      64             : !-----> evaluate relativistic corrections for exchange
      65             : 
      66           0 :       ALLOCATE ( psi(ngrid) )
      67             :       CALL relcor( &
      68             :          mgrid,ngrid,jspins,krla, .TRUE. ,rh, &
      69           0 :          psi)
      70             : 
      71           0 :       DO ir = 1,ngrid
      72           0 :          IF (jspins == 1) THEN
      73           0 :             rho = max(d_15,rh(ir,1))
      74             :          ELSE
      75           0 :             rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
      76             :          ENDIF
      77           0 :          rs  = (thfpi/rho)**thrd
      78           0 :          exp = - cex / rs              ! exchange energy = -0.9163306/rs
      79           0 :          ecp = - 0.88 / (rs + 7.8)     ! correlation energy = -0.88/(rs + 7.8)
      80             : 
      81           0 :          vxp = fothrd * exp            ! exchange potential = 4/3 exp
      82           0 :          vcp = fothrd * ecp + 2.288 / ( rs + 7.8 ) ** 2
      83             : 
      84           0 :          vxc(ir,1) = vcp + vxp * psi(ir)
      85             : 
      86           0 :          vx(ir,1)  = vxp * psi(ir)
      87             :       ENDDO
      88           0 :       IF ( jspins == 2 ) THEN        ! spinpolarized calculation
      89           0 :          DO ir = 1,ngrid
      90           0 :             vxc(ir,jspins) = vxc(ir,1)
      91             : 
      92           0 :             vx(ir,jspins)  = vx(ir,1)
      93             :          ENDDO
      94             :       ENDIF
      95             : 
      96           0 :       DEALLOCATE (psi)
      97           0 :       RETURN
      98             : 
      99             :    END SUBROUTINE vxcwgn
     100             : !***********************************************************************
     101           0 :    SUBROUTINE excwgn( &
     102             :       iofile,krla,jspins, &
     103           0 :       mgrid,ngrid,rh, &
     104           0 :       exc)
     105             : !***********************************************************************
     106             : 
     107             : !     .. Scalar Arguments ..
     108             :       INTEGER, INTENT (IN) :: jspins
     109             :       INTEGER, INTENT (IN) :: krla        !  run mode parameters
     110             :       INTEGER, INTENT (IN) :: iofile      !  file number for read and write
     111             :       INTEGER, INTENT (IN) :: mgrid,ngrid !  mesh points
     112             : 
     113             : !     .. Array Arguments ..
     114             :       REAL, INTENT (IN)  :: rh(mgrid,jspins)      ! charge density
     115             :       REAL, INTENT (OUT) :: exc(mgrid)            ! xc energy
     116             : 
     117             : !     .. Local Arrays ..
     118           0 :       REAL, ALLOCATABLE :: phi(:)       ! relativistic exchange energy correct.
     119             : 
     120             :       REAL ::  rs, rho, thfpi, exp, ecp
     121             :       INTEGER :: ir
     122             : 
     123           0 :       thfpi  = three / ( four * pi_const )
     124             : 
     125           0 :       ALLOCATE ( phi(ngrid) )
     126             :       CALL relcor( &
     127             :          mgrid,ngrid,jspins,krla, .FALSE. ,rh, &
     128           0 :          phi)
     129             : 
     130           0 :       DO ir = 1,ngrid
     131           0 :          IF (jspins == 1) THEN
     132           0 :             rho = max(d_15,rh(ir,1))
     133             :          ELSE
     134           0 :             rho = max(d_15,rh(ir,1)) + max(d_15,rh(ir,jspins))
     135             :          ENDIF
     136           0 :          rs  = (thfpi/rho)**thrd
     137           0 :          exp = - cex / rs              ! exchange energy = -0.9163306/rs
     138           0 :          ecp = - 0.88 / (rs + 7.8)     ! correlation energy = -0.88/(rs + 7.8)
     139             : 
     140           0 :          exc(ir) = ecp + exp * phi(ir)
     141             :       ENDDO
     142           0 :       IF (jspins == 2) THEN
     143             :          WRITE (iofile,'('' WARNING: Wigner correlation !'', &
     144           0 :                 ''only applicable for non-spinpolarized calculations'')')
     145             :       ENDIF
     146             : 
     147           0 :       DEALLOCATE (phi)
     148           0 :       RETURN
     149             : 
     150             :    END SUBROUTINE excwgn
     151             : END MODULE m_xcwgn

Generated by: LCOV version 1.13