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

          Line data    Source code
       1             : MODULE m_excwb91
       2             : !.....-----------------------------------------------------------------
       3             : !.....pw91 exchange-correlation energy density in hartree.
       4             : !.....------------------------------------------------------------------
       5             : CONTAINS
       6           0 :    SUBROUTINE excwb91( &
       7           0 :       mirm,irmx,rh1,rh2,agr,agru,agrd, &
       8           0 :       g2r,g2ru,g2rd,gggr,gggru,gggrd,gzgr, &
       9           0 :       exc, &
      10             :       idsprs,isprsv,sprsv)
      11             : 
      12             :       USE m_corl91
      13             :       USE m_corg91
      14             :       USE m_xch91
      15             :       USE m_constants
      16             : 
      17             :       IMPLICIT NONE
      18             : 
      19             : ! .. Arguments ..
      20             :       INTEGER, INTENT (IN) :: irmx,mirm,idsprs,isprsv
      21             :       REAL,    INTENT (IN) :: sprsv
      22             : 
      23             :       REAL,    INTENT (IN) :: rh1(mirm),rh2(mirm)
      24             :       REAL,    INTENT (IN) :: agr(mirm),agru(mirm),agrd(mirm)
      25             :       REAL,    INTENT (IN) :: g2r(mirm),g2ru(mirm),g2rd(mirm)
      26             :       REAL,    INTENT (IN) :: gggr(mirm),gggru(mirm)
      27             :       REAL,    INTENT (IN) :: gggrd(mirm),gzgr(mirm)
      28             :       REAL,    INTENT (OUT) :: exc(mirm)
      29             : 
      30             : ! .. local variables
      31             :       REAL :: ro,zta,alf,alfc,c13,c23,c43,c53, &
      32             :               cedg,cedl,dbrod,dbrou,dsprs,ec,ecrs,eczta,fk,gz, &
      33             :               ro13,ro2,rod,rod3,rod43,rod53,rou,rou3,rou43,rou53, &
      34             :               rs,sd,sk,su,tc,td,tksg,tu,uc,ud,uu,vc, &
      35             :               vcgd,vcgu,vcld,vclu,vxgd,vxgu,vxld, &
      36             :               vxlu,wc,xced,xedg,xedgd,xedgu,xedl,xedld,xedlu
      37             : 
      38             :       INTEGER :: i
      39             :       REAL, PARAMETER :: sml = 1.e-14
      40             :       REAL, PARAMETER :: huges = 1.e+6
      41             : !.....------------------------------------------------------------------
      42             : !.....huges: if su(sdw).gt.huges, grad.expansion calc.non-excuted.
      43             : !.....-----------------------------------------------------------------
      44             : 
      45           0 :       do 300 i = 1,irmx
      46             : 
      47           0 :          ro=rh1(i)+rh2(i)
      48           0 :          zta=(rh1(i)-rh2(i))/ro
      49           0 :          if(zta > 1.0-sml) zta = 1.0 - sml
      50           0 :          if(zta < -1.0-sml) zta = -1.0 + sml
      51             : 
      52             :          !.....
      53             :          !       xedl,xedg: exchange energy density (local,grad.exp.) in ry.
      54             :          !       cedl,cedg: exchange energy density (local,grad.expnd.) in ry.
      55             :          !       all later in hartree.
      56             :          !.....
      57           0 :          xedl = 0.0e0
      58           0 :          cedl = 0.0e0
      59           0 :          xedg = 0.0e0
      60           0 :          cedg = 0.0e0
      61             : 
      62             :          !.....
      63           0 :          if(ro < sml) then
      64           0 :             ro = sml
      65           0 :             zta = 0.0
      66           0 :             go to 200
      67             :          endif
      68             :          !.....
      69           0 :          c13 = 1.e0/3.e0
      70           0 :          c23 = 2.e0/3.e0
      71           0 :          c43 = 4.e0/3.e0
      72           0 :          c53 = 5.e0/3.e0
      73             :          !.....  alf=-3*(3/4*pai)**(1/3).
      74           0 :          alf = -1.861051473e0
      75             :          !.....
      76           0 :          ro2 = ro*ro
      77           0 :          ro13 = ro**c13
      78             :          !.....
      79           0 :          rou = rh1(i)
      80           0 :          rou3 = rou**3
      81           0 :          rou43 = rou**c43
      82             :          !.....
      83           0 :          rod = rh2(i)
      84           0 :          rod3 = rod**3
      85           0 :          rod43 = rod**c43
      86             :          !.....
      87             :          !       gr2=drr*drr
      88             :          !       gr2u=drru**2
      89             :          !       drrd=drr-drru
      90             :          !       gr2d=drrd**2
      91             :          !       ddrrd=ddrr-ddrru
      92             :          !.....
      93             :          !.....  gz,gz2,gz3: for wang-perdew ssf.
      94           0 :          gz = ((1.0+zta)**c23+ (1.0-zta)**c23)/2.e0
      95             :          !.....
      96           0 :          rs = 0.620350491/ro13
      97             :          !.....
      98             :          !.....  xedl: exchange-energy-density in ry.
      99           0 :          xedl = alf* (rou43+rod43)
     100             :          !.....
     101             :          !....  .gradient correction.
     102             :          !.....
     103           0 :          if(abs(agr(i)) < sml) go to 200
     104             :          !.....
     105             : 
     106             :          !.....
     107           0 :          dsprs = 1.0
     108           0 :          if(idsprs == 1) dsprs = 1.e-19
     109           0 :          if(isprsv == 1) dsprs = dsprs*sprsv
     110             : 
     111             :          !.....
     112             :          !      agr,agru,agrd: abs(grad(rho)), for all, up, and down.
     113             :          !c     gr2,gr2u,gr2d: grad(rho_all)**2, grad(rho_up)**2, grad(rho_d)**2.
     114             :          !      g2r,g2ru,g2rd: laplacian rho_all, _up and _down.
     115             :          !      gggru,-d: grad(rho)*grad(abs(grad(rho))) for all,up and down.
     116             :          !      grgru,-d: grad(rho_all)*grad(rhor_up) and for down.
     117             : 
     118             :          !         g2r=ddrr+2*drr/rv
     119             :          !.....
     120           0 :          rou53 = rou**c53
     121             :          !.....
     122             :          !.....  edrru: d(abs(d(rou)/dr))/dr, edrrd for down.
     123             :          !         edrru=ddrru
     124             :          !         if(drru.lt.0.) edrru=-ddrru
     125             :          !.....
     126             :          !         agr,agbru,-d: abs(grad(rho)),for rou, rod.
     127             :          !         gggru,-d: grad(rho)*grad(abs(grad(rho))) for up and down.
     128             :          !.....  su:at ro=2*rou. 1/(2(3*pai**2)**(1/3))*|grad(rou)|/rou**(4/3).
     129           0 :          su = 0.128278244*agru(i)/rou43
     130           0 :          if(su > huges) go to 200
     131             :          !         g2ru=ddrru+2*drru/rv
     132           0 :          tu = .016455307e0*g2ru(i)/rou53
     133           0 :          uu = 0.002110857e0*gggru(i)/rou3
     134             : 
     135           0 :          dbrou = rou*2
     136             : 
     137           0 :          call xch91(dbrou,su,uu,tu,xedlu,xedgu,vxlu,vxgu)
     138             : 
     139           0 :          xedl = xedlu/2
     140           0 :          xedgu = dsprs*xedgu
     141             : 
     142             :          !.....
     143           0 :          rod53 = rod**c53
     144             :          !         edrrd=ddrrd
     145             :          !         if(drrd.lt.0.) edrrd=-ddrrd
     146             : 
     147           0 :          sd = 0.128278244*agrd(i)/rod43
     148           0 :          if(sd > huges) go to 200
     149             : 
     150             :          !         g2rd=ddrrd+2*drrd/rv
     151             : 
     152           0 :          td = .016455307e0*g2rd(i)/rod53
     153           0 :          ud = 0.002110857e0*gggrd(i)/rod3
     154             : 
     155           0 :          dbrod = rod*2
     156             : 
     157           0 :          call xch91(dbrod,sd,ud,td,xedld,xedgd,vxld,vxgd)
     158             : 
     159           0 :          xedl = xedl + xedld/2
     160           0 :          xedgd = dsprs*xedgd
     161             : 
     162           0 :          xedg = (xedgu+xedgd)/2
     163             : 
     164             :          !....   cro: c(n) of (6),phys.rev..b33,8822('86). in ry.
     165             :          !....   dcdr: d(cro)/d(ro).
     166             :          !.....  0.001625816=1.745*f(=0.11)*cro(rs=0).
     167             : 
     168           0 :          call corl91(rs,zta,ec,vclu,vcld,ecrs,eczta,alfc)
     169             : 
     170           0 :          cedl = ec*2.e0*ro
     171             : 
     172           0 :          fk = 1.91915829/rs
     173           0 :          sk = sqrt(4.e0*fk/pi_const)
     174           0 :          tksg = 2.e0*sk*gz
     175           0 :          tc = agr(i)/ (ro*tksg)
     176           0 :          uc = gggr(i)/ (ro2*tksg**3)
     177           0 :          vc = g2r(i)/ (ro*tksg**2)
     178           0 :          wc = gzgr(i)/ (ro*tksg**2)
     179             : 
     180             :          call corg91(fk,sk,gz,ec,ecrs,eczta,rs,zta,tc,uc,vc,wc,cedg, &
     181           0 :                      vcgu,vcgd)
     182             : 
     183           0 :          cedg = dsprs*cedg*ro*2.e0
     184             : 
     185             : 200      continue
     186             : 
     187           0 :          xced = (xedl+cedl+xedg+cedg)
     188           0 :          if(ro >= sml) xced = xced/ro
     189             : 
     190             :          !     change to hartree.
     191           0 :          exc(i) = xced/2.e0
     192             : 
     193           0 : 300   END DO
     194             : 
     195           0 :    END SUBROUTINE excwb91
     196             : END MODULE m_excwb91

Generated by: LCOV version 1.13