Line data    Source code

       1             :     MODULE m_corl91
       2             : !.....-----------------------------------------------------------------
       3             : !     uniform-gas correlation of perdew and wang 1991
       4             : !.....-----------------------------------------------------------------
       5             :     CONTAINS
       6           0 :     SUBROUTINE corl91( &
       7             :     rs,zta, &
       8             :     ec,vcup,vcdn,ecrs,eczta,alfc)
       9             : !.....-----------------------------------------------------------------
      10             : !     input: seitz radius (rs), relative spin polarization (zta)
      11             : !     output: correlation energy per electron (ec),
      12             : !             up- and down-spin potentials (vcup,vcdn),
      13             : !             derivatives of ec wrt rs (ecrs) &zta (eczta).
      14             : !     output: correlation contribution (alfc) to the spin stiffness
      15             : !.....-----------------------------------------------------------------
      16             :     IMPLICIT NONE
      17             : 
      18             :     REAL, INTENT (IN)  :: rs,zta
      19             :     REAL, INTENT (OUT) :: alfc,ec,ecrs,eczta,vcdn,vcup
      20             : 
      21             :     REAL :: alfm,alfrsm,comm,ep,eprs,eu,eurs,f,fz,z4
      22             :     REAL :: fzz,gam,thrd,thrd4
      23             : !.....-----------------------------------------------------------------
      24             : !     ..
      25             :     DATA gam,fzz/0.5198421,1.709921/
      26             :     DATA thrd,thrd4/0.333333333333e0,1.333333333333e0/
      27             : !.....-----------------------------------------------------------------
      28           0 :     f = ((1.0+zta)**thrd4+ (1.0-zta)**thrd4-2.e0)/gam
      29             : 
      30             :     CALL gcor91(0.0310907,0.21370,7.5957,3.5876,1.6382, &
      31           0 :     &           0.49294,1.00,rs,eu,eurs)
      32             : 
      33             :     CALL gcor91(0.01554535,0.20548,14.1189,6.1977,3.3662, &
      34           0 :     &           0.62517,1.00,rs,ep,eprs)
      35             : 
      36             :     CALL gcor91(0.0168869,0.11125,10.357,3.6231,0.88026, &
      37           0 :     &           0.49671,1.00,rs,alfm,alfrsm)
      38             : 
      39             : !  alfm is minus the spin stiffness alfc
      40           0 :     alfc = -alfm
      41           0 :     z4 = zta**4
      42           0 :     ec = eu* (1.0-f*z4) + ep*f*z4 - alfm*f* (1.0-z4)/fzz
      43             : !  energy done. now the potential:
      44           0 :     ecrs = eurs* (1.0-f*z4) + eprs*f*z4 - alfrsm*f* (1.0-z4)/fzz
      45           0 :     fz = thrd4* ((1.0+zta)**thrd- (1.0-zta)**thrd)/gam
      46             :     eczta = 4.e0* (zta**3)*f* (ep-eu+alfm/fzz) + &
      47           0 :     fz* (z4*ep-z4*eu- (1.0-z4)*alfm/fzz)
      48           0 :     comm = ec - rs*ecrs/3.e0 - zta*eczta
      49           0 :     vcup = comm + eczta
      50           0 :     vcdn = comm - eczta
      51             : 
      52           0 :     END SUBROUTINE corl91
      53             : !.....-----------------------------------------------------------------
      54             : !     called by corl91
      55             : !.....-----------------------------------------------------------------
      56           0 :     SUBROUTINE gcor91( &
      57             :     a,a1,b1,b2,b3,b4,p,rs, &
      58             :     gg,ggrs)
      59             : 
      60             :     IMPLICIT NONE
      61             :     REAL, INTENT (IN)  :: a,a1,b1,b2,b3,b4,p,rs
      62             :     REAL, INTENT (OUT) :: gg,ggrs
      63             : 
      64             :     REAL :: p1,q0,q1,q2,q3,rs12,rs32,rsp
      65             : !.....-----------------------------------------------------------------
      66             : !     ..
      67           0 :     p1 = p + 1.0
      68           0 :     q0 = -2.e0*a* (1.0+a1*rs)
      69           0 :     rs12 = sqrt(rs)
      70           0 :     rs32 = rs12**3
      71           0 :     rsp = rs**p
      72           0 :     q1 = 2.e0*a* (b1*rs12+b2*rs+b3*rs32+b4*rs*rsp)
      73           0 :     q2 = log(1.0+1.0/q1)
      74           0 :     gg = q0*q2
      75           0 :     q3 = a* (b1/rs12+2.e0*b2+3.e0*b3*rs12+2.e0*b4*p1*rsp)
      76           0 :     ggrs = -2.e0*a*a1*q2 - q0*q3/ (q1**2+q1)
      77             : 
      78           0 :     END SUBROUTINE gcor91
      79             : 
      80             :     END MODULE m_corl91