Clean up functional selection

vanderhe · Apr 20, 2024 · 6bb04ed · 6bb04ed
1 parent 097a928
commit 6bb04ed
Show file tree

Hide file tree

Showing 5 changed files with 35 additions and 22 deletions.
diff --git a/sktwocnt/lib/twocnt.f90 b/sktwocnt/lib/twocnt.f90
@@ -229,46 +229,47 @@ subroutine get_twocenter_integrals(inp, imap, skham, skover)
     !! number of radial and angular integration abscissas
     integer :: nRad, nAng
 
-    if (inp%iXC == xcFunctional%LDA_PW91) then
+    select case (inp%iXC)
+    case(xcFunctional%LDA_PW91)
       call xc_f03_func_init(xcfunc_x, XC_LDA_X, XC_UNPOLARIZED)
       call xc_f03_func_init(xcfunc_c, XC_LDA_C_PW, XC_UNPOLARIZED)
-    elseif (inp%iXC == xcFunctional%GGA_PBE96) then
+    case(xcFunctional%GGA_PBE96)
       call xc_f03_func_init(xcfunc_x, XC_GGA_X_PBE, XC_UNPOLARIZED)
       call xc_f03_func_init(xcfunc_c, XC_GGA_C_PBE, XC_UNPOLARIZED)
-    elseif (inp%iXC == xcFunctional%GGA_BLYP) then
+    case(xcFunctional%GGA_BLYP)
       call xc_f03_func_init(xcfunc_x, XC_GGA_X_B88, XC_UNPOLARIZED)
       call xc_f03_func_init(xcfunc_c, XC_GGA_C_LYP, XC_UNPOLARIZED)
-    elseif (inp%iXC == xcFunctional%LCY_PBE96) then
+    case(xcFunctional%LCY_PBE96)
       call xc_f03_func_init(xcfunc_x, XC_GGA_X_SFAT_PBE, XC_UNPOLARIZED)
       call xc_f03_func_set_ext_params(xcfunc_x, [inp%omega])
       call xc_f03_func_init(xcfunc_c, XC_GGA_C_PBE, XC_UNPOLARIZED)
-    elseif (inp%iXC == xcFunctional%LCY_BNL) then
+    case(xcFunctional%LCY_BNL)
       call xc_f03_func_init(xcfunc_x, XC_LDA_X_YUKAWA, XC_UNPOLARIZED)
       call xc_f03_func_set_ext_params(xcfunc_x, [inp%omega])
       call xc_f03_func_init(xcfunc_c, XC_GGA_C_PBE, XC_UNPOLARIZED)
-    elseif (inp%iXC == 6) then
+    case(xcFunctional%HYB_PBE0)
       ! xpbe96
       call xc_f03_func_init(xcfunc_x, XC_GGA_X_PBE, XC_UNPOLARIZED)
       ! cpbe96
       call xc_f03_func_init(xcfunc_c, XC_GGA_C_PBE, XC_UNPOLARIZED)
-    elseif (inp%iXC == 7) then
+    case(xcFunctional%HYB_B3LYP)
       call xc_f03_func_init(xcfunc_xc, XC_HYB_GGA_XC_B3LYP, XC_UNPOLARIZED)
       ! 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
+    case(xcFunctional%CAMY_B3LYP)
       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,&
           & -inp%camBeta, inp%omega])
-    elseif (inp%iXC == 9) then
+    case(xcFunctional%CAMY_PBEh)
       ! short-range xpbe96
       call xc_f03_func_init(xcfunc_xsr, XC_GGA_X_SFAT_PBE, XC_UNPOLARIZED)
       call xc_f03_func_set_ext_params(xcfunc_xsr, [inp%omega])
       ! xpbe96
       call xc_f03_func_init(xcfunc_x, XC_GGA_X_PBE, XC_UNPOLARIZED)
       ! cpbe96
       call xc_f03_func_init(xcfunc_c, XC_GGA_C_PBE, XC_UNPOLARIZED)
-    end if
+    end select
 
     if (inp%tLC .or. inp%tCam) then
       beckeGridParams%nRadial = inp%nRadial
@@ -378,11 +379,11 @@ subroutine get_twocenter_integrals(inp, imap, skham, skover)
 
     ! finalize libxc objects
     if (inp%tGlobalHybrid .or. inp%tCam) then
-      if (inp%iXC == 9) then
+      if (inp%iXC == xcFunctional%CAMY_PBEh) then
         call xc_f03_func_end(xcfunc_xsr)
         call xc_f03_func_end(xcfunc_x)
         call xc_f03_func_end(xcfunc_c)
-      elseif (inp%iXC == 6) then
+      elseif (inp%iXC == xcFunctional%HYB_PBE0) then
         call xc_f03_func_end(xcfunc_x)
         call xc_f03_func_end(xcfunc_c)
       else
@@ -572,7 +573,7 @@ subroutine getskintegrals(beckeInt, radialHFQuadrature, nRad, nAng, atom1, atom2
       end if
 
       ! CAMY-PBEh is assembled manually
-      if (iXC == 9) then
+      if (iXC == xcFunctional%CAMY_PBEh) then
         allocate(vxsigma_sr(nGrid))
         vxsigma_sr(:) = 0.0_dp
         allocate(vx_sr(nGrid))

diff --git a/sktwocnt/prog/input.f90 b/sktwocnt/prog/input.f90
@@ -90,10 +90,10 @@ subroutine readInput(inp, fname)
     case(xcFunctional%GGA_BLYP)
       ! GGA-BLYP
     case(xcFunctional%LCY_PBE96)
-      ! LCY-PBE96 (long-range corrected)
+      ! LCY-PBE96 (purely long-range corrected)
       inp%tLC = .true.
     case(xcFunctional%LCY_BNL)
-      ! LCY-BNL (long-range corrected)
+      ! LCY-BNL (purely long-range corrected)
       inp%tLC = .true.
     case(xcFunctional%HYB_PBE0)
       ! PBE0 (global hybrid)
@@ -102,18 +102,18 @@ subroutine readInput(inp, fname)
       ! B3LYP (global hybrid)
       inp%tGlobalHybrid = .true.
     case(xcFunctional%CAMY_B3LYP)
-      ! CAMY-B3LYP (CAM-functional)
+      ! CAMY-B3LYP (general CAM form)
       inp%tCam = .true.
     case(xcFunctional%CAMY_PBEh)
-      ! CAMY-PBEh (CAM-functional)
+      ! CAMY-PBEh (general CAM form)
       inp%tCam = .true.
     case default
       call error_("Unknown exchange-correlation functional!", fname, line, iline)
     end select
     inp%iXC = iXC
 
     if (inp%iXC == xcFunctional%HYB_B3LYP) then
-      ! 20% HFX hard-coded at the moment
+      ! 20% fraction of HFX hard-coded at the moment
       inp%camAlpha = 0.2_dp
       inp%camBeta = 0.0_dp
       call nextline_(fp, iLine, line)
@@ -129,6 +129,8 @@ subroutine readInput(inp, fname)
       read(line, *, iostat=iErr) inp%nRadial, inp%nAngular, inp%ll_max, inp%rm
       call checkerror_(fname, line, iLine, iErr)
     elseif (inp%tLC) then
+      inp%camAlpha = 0.0_dp
+      inp%camBeta = 1.0_dp
       call nextline_(fp, iLine, line)
       read(line, *, iostat=iErr) inp%omega
       if (inp%omega < 1.0e-08_dp) then

diff --git a/slateratom/lib/dft.f90 b/slateratom/lib/dft.f90
@@ -210,7 +210,7 @@ subroutine density_grid(pp, max_l, num_alpha, poly_order, alpha, num_mesh_points
     case(xcFunctional%HYB_PBE0)
       call getExcVxc_HYB_PBE0(abcissa, dz, dzdr, rho, drho, sigma, camAlpha, exc, vxc)
     case(xcFunctional%HYB_B3LYP)
-      call getExcVxc_HYB_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, exc, vxc)
+      call getExcVxc_HYB_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, camAlpha, exc, vxc)
     case(xcFunctional%CAMY_B3LYP)
       call getExcVxc_CAMY_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, omega, camAlpha, camBeta, exc,&
           & vxc)

diff --git a/slateratom/lib/input.f90 b/slateratom/lib/input.f90
@@ -101,6 +101,10 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     !! auxiliary variables
     integer :: ii, jj
 
+    omega = 0.0_dp
+    camAlpha = 0.0_dp
+    camBeta = 0.0_dp
+
     write(*, '(A)') 'Enter nuclear charge, maximal angular momentum (s=0), max. SCF, SCF tol., ZORA'
     read(*,*) nuc, max_l, maxiter, scftol, tZora
 
@@ -115,12 +119,15 @@ subroutine read_input_1(nuc, max_l, occ_shells, maxiter, scftol, poly_order, min
     end if
 
     if (xcFunctional%isLongRangeCorrected(xcnr)) then
+      camBeta = 1.0_dp
       write(*, '(A)') 'Enter range-separation parameter:'
       read(*,*) omega
     elseif (xcnr == xcFunctional%HYB_PBE0) then
       ! currently only HYB-PBE0 does support arbitrary HFX portions (HYB-B3LYP does not)
       write(*, '(A)') 'Enter portion of HFX (CAM alpha):'
       read(*,*) camAlpha
+    elseif (xcnr == xcFunctional%HYB_B3LYP) then
+      camAlpha = 0.2_dp
     elseif (xcFunctional%isCAMY(xcnr)) then
       write(*, '(A)') 'Enter range-separation parameter, CAM alpha, CAM beta:'
       read(*,*) omega, camAlpha, camBeta

diff --git a/slateratom/lib/xcfunctionals.f90 b/slateratom/lib/xcfunctionals.f90
@@ -713,7 +713,7 @@ end subroutine getExcVxc_HYB_PBE0
 
 
   !> Calculates exc and vxc for the HYB-B3LYP xc-functional.
-  subroutine getExcVxc_HYB_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, exc, vxc)
+  subroutine getExcVxc_HYB_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, camAlpha, exc, vxc)
 
     !> numerical integration abcissas
     real(dp), intent(in) :: abcissa(:)
@@ -733,6 +733,9 @@ subroutine getExcVxc_HYB_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, exc, vxc)
     !> contracted gradients of the density
     real(dp), intent(in), allocatable :: sigma(:,:)
 
+    !> CAM alpha parameter
+    real(dp), intent(in) :: camAlpha
+
     !> exc energy density on grid
     real(dp), intent(out) :: exc(:)
 
@@ -772,9 +775,9 @@ 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
+    ! 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, [0.20_dp, 0.72_dp, 0.81_dp])
+    call xc_f03_func_set_ext_params(xcfunc_xc, [camAlpha, 0.72_dp, 0.81_dp])
 
     ! exchange + correlation
     call xc_f03_gga_exc_vxc(xcfunc_xc, nn, rhor(1, 1), sigma(1, 1), exc_tmp(1), vxc_tmp(1, 1),&