Fix bug in sktwocnt for the global hybrid B3LYP

sktwocnt/lib/twocnt.f90

@@ -253,7 +253,9 @@ subroutine get_twocenter_integrals(inp, imap, skham, skover)
       call xc_f03_func_init(xcfunc_c, XC_GGA_C_PBE, XC_UNPOLARIZED)
     elseif (inp%iXC == 7) then
       call xc_f03_func_init(xcfunc_xc, XC_HYB_GGA_XC_B3LYP, XC_UNPOLARIZED)
-      call xc_f03_func_set_ext_params(xcfunc_xc, [0.20_dp, 0.72_dp, 0.81_dp])
+      ! Adjustable fraction of Fock-type exchange, otherwise standard parametrization taken from
+      ! J. Phys. Chem. 1994, 98, 45, 11623-11627; DOI: 10.1021/j100096a001
+      call xc_f03_func_set_ext_params(xcfunc_xc, [inp%camAlpha, 0.72_dp, 0.81_dp])
     elseif (inp%iXC == 8) then
       call xc_f03_func_init(xcfunc_xc, XC_HYB_GGA_XC_CAMY_B3LYP, XC_UNPOLARIZED)
       call xc_f03_func_set_ext_params(xcfunc_xc, [0.81_dp, inp%camAlpha + inp%camBeta,&
@@ -657,17 +659,11 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
       ! total density: \int (|\phi_1|^2 + |\phi_2|^2)
       dens = getDensity(radval1(:, i1), radval2(:, i2), spherval1, spherval2, weights)
 
-      if (iXC == xcFunctional%HYB_B3LYP) then
-        ! full-range Hartree-Fock exchange contribution
-        frx = 0.5_dp * getFullRangeHFContribution(radialHFQuadrature%xx, rr3, ll_max, atom1, atom2,&
-            & imap, ii, r1, theta1, r2, theta2, weights)
-        ! add up full-range exchange to the Hamiltonian
-        integ1 = integ1 - frx
-      elseif (iXC == xcFunctional%HYB_PBE0) then
+      if (iXC == xcFunctional%HYB_PBE0 .or. iXC == xcFunctional%HYB_B3LYP) then
         ! full-range Hartree-Fock exchange contribution
         frx = 0.5_dp * camAlpha * getFullRangeHFContribution(radialHFQuadrature%xx, rr3, ll_max,&
             & atom1, atom2, imap, ii, r1, theta1, r2, theta2, weights)
-        ! add up full-/long-range exchange to the Hamiltonian
+        ! add up full-range exchange to the Hamiltonian
         integ1 = integ1 - frx
       elseif (tLC) then
         ! long-range Hartree-Fock exchange contribution

sktwocnt/prog/input.f90

@@ -113,10 +113,14 @@ subroutine readInput(inp, fname)
     inp%iXC = iXC
 
     if (inp%iXC == xcFunctional%HYB_B3LYP) then
+      ! 20% HFX hard-coded at the moment
+      inp%camAlpha = 0.2_dp
+      inp%camBeta = 0.0_dp
       call nextline_(fp, iLine, line)
       read(line, *, iostat=iErr) inp%nRadial, inp%nAngular, inp%ll_max, inp%rm
       call checkerror_(fname, line, iLine, iErr)
     elseif (inp%iXC == xcFunctional%HYB_PBE0) then
+      inp%camBeta = 0.0_dp
       call nextline_(fp, iLine, line)
       ! currently only HYB-PBE0 does support arbitrary HFX portions (HYB-B3LYP does not)
       read(line, *, iostat=iErr) inp%camAlpha

slateratom/lib/xcfunctionals.f90

@@ -772,6 +772,8 @@ subroutine getExcVxc_HYB_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, exc, vxc)
     vxcsigma(:,:) = 0.0_dp
 
     call xc_f03_func_init(xcfunc_xc, XC_HYB_GGA_XC_B3LYP, XC_POLARIZED)
+    ! Standard parametrization of B3LYP taken from
+    ! J. Phys. Chem. 1994, 98, 45, 11623-11627; DOI: 10.1021/j100096a001
     call xc_f03_func_set_ext_params(xcfunc_xc, [0.20_dp, 0.72_dp, 0.81_dp])
 
     ! exchange + correlation