Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Clean up functional selection #76

Merged
merged 1 commit into from
May 3, 2024
Merged

Clean up functional selection #76

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 14 additions & 13 deletions sktwocnt/lib/twocnt.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand Down
12 changes: 7 additions & 5 deletions sktwocnt/prog/input.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion slateratom/lib/dft.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
bhourahine marked this conversation as resolved.
Show resolved Hide resolved
case(xcFunctional%CAMY_B3LYP)
call getExcVxc_CAMY_B3LYP(abcissa, dz, dzdr, rho, drho, sigma, omega, camAlpha, camBeta, exc,&
& vxc)
Expand Down
7 changes: 7 additions & 0 deletions slateratom/lib/input.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand Down
9 changes: 6 additions & 3 deletions slateratom/lib/xcfunctionals.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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(:)
Expand All @@ -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(:)

Expand Down Expand Up @@ -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),&
Expand Down
Loading