LCOV - code coverage report
Current view: top level - xc-pot - vxcwb91.f90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 0.0 % 83 0
Test Date: 2025-07-19 04:34:44 Functions: 0.0 % 1 0

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

Generated by: LCOV version 2.0-1