diff --git a/sktwocnt/lib/twocnt.f90 b/sktwocnt/lib/twocnt.f90 index a3420957..9c63bf24 100644 --- a/sktwocnt/lib/twocnt.f90 +++ b/sktwocnt/lib/twocnt.f90 @@ -229,38 +229,39 @@ 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]) @@ -268,7 +269,7 @@ subroutine get_twocenter_integrals(inp, imap, skham, skover) 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 index 8f8a90d3..1c636b49 100644 --- 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,10 +102,10 @@ 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) @@ -113,7 +113,7 @@ subroutine readInput(inp, fname) 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 index 4568acc4..cda52119 100644 --- 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 index bb2b75fc..a5ddc130 100644 --- 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 index 8ef864b0..cb20cd34 100644 --- 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),&