Skip to content

Commit

Permalink
Low spin ROKS with hybrid functionals (no ADMM) + regtests (#2582)
Browse files Browse the repository at this point in the history
  • Loading branch information
@juerghutter
juerghutter committed Feb 17, 2023
1 parent 8c73008 commit 04096ff
Show file tree
Hide file tree
Showing 8 changed files with 340 additions and 15 deletions.
14 changes: 13 additions & 1 deletion src/qs_environment.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1544,6 +1544,18 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
"Try UKS instead of ROKS")
END IF
END IF
IF (dft_control%low_spin_roks) THEN
SELECT CASE (dft_control%qs_control%method_id)
CASE DEFAULT
CASE (do_method_xtb, do_method_dftb)
CALL cp_abort(__LOCATION__, &
"xTB/DFTB methods are not compatible with low spin ROKS.")
CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pm3, &
do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
CALL cp_abort(__LOCATION__, &
"SE methods are not compatible with low spin ROKS.")
END SELECT
END IF

! in principle the restricted calculation could be performed
! using just one set of MOs and special casing most of the code
Expand Down Expand Up @@ -1576,7 +1588,7 @@ SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell
CALL set_qs_env(qs_env, mos=mos)
! allocate mos when switch_surf_dip is triggered [SGh]
! allocate mos when switch_surf_dip is triggered [SGh]
IF (dft_control%switch_surf_dip) THEN
ALLOCATE (mos_last_converged(dft_control%nspins))
DO ispin = 1, dft_control%nspins
Expand Down
2 changes: 1 addition & 1 deletion src/qs_ks_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -850,7 +850,7 @@ SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
IF (calculate_forces .AND. dft_control%do_admm) CALL calc_admm_ovlp_forces(qs_env)

! deal with low spin roks
CALL low_spin_roks(energy, qs_env, dft_control, just_energy, &
CALL low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
calculate_forces, auxbas_pw_pool)

! deal with sic on explicit orbitals
Expand Down
99 changes: 90 additions & 9 deletions src/qs_ks_utils.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,9 @@ MODULE qs_ks_utils
dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_dot, dbcsr_get_info, dbcsr_init_p, &
dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, dbcsr_scale_by_vector, &
dbcsr_set, dbcsr_type
USE hfx_admm_utils, ONLY: tddft_hfx_matrix
USE hfx_derivatives, ONLY: derivatives_four_center
USE hfx_types, ONLY: hfx_type
USE input_constants, ONLY: &
cdft_alpha_constraint, cdft_beta_constraint, cdft_charge_constraint, &
cdft_magnetization_constraint, do_admm_aux_exch_func_none, do_admm_exch_scaling_merlot, &
Expand Down Expand Up @@ -142,49 +145,74 @@ MODULE qs_ks_utils
!> \param energy ...
!> \param qs_env ...
!> \param dft_control ...
!> \param do_hfx ...
!> \param just_energy ...
!> \param calculate_forces ...
!> \param auxbas_pw_pool ...
! **************************************************************************************************
SUBROUTINE low_spin_roks(energy, qs_env, dft_control, just_energy, &
SUBROUTINE low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
calculate_forces, auxbas_pw_pool)

TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dft_control_type), POINTER :: dft_control
LOGICAL, INTENT(IN) :: just_energy, calculate_forces
LOGICAL, INTENT(IN) :: do_hfx, just_energy, calculate_forces
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool

CHARACTER(*), PARAMETER :: routineN = 'low_spin_roks'

INTEGER :: handle, ispin, iterm, k, k_alpha, &
INTEGER :: handle, irep, ispin, iterm, k, k_alpha, &
k_beta, n_rep, Nelectron, Nspin, Nterms
INTEGER, DIMENSION(:), POINTER :: ivec
INTEGER, DIMENSION(:, :, :), POINTER :: occupations
LOGICAL :: compute_virial, in_range, &
uniform_occupation
REAL(KIND=dp) :: exc
REAL(KIND=dp) :: ehfx, exc
REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp
REAL(KIND=dp), DIMENSION(:), POINTER :: energy_scaling, rvec, scaling
TYPE(cell_type), POINTER :: cell
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_p, mo_derivs, rho_ao
TYPE(cp_para_env_type), POINTER :: para_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_hfx, matrix_p, mdummy, &
mo_derivs, rho_ao
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p2
TYPE(dbcsr_type), POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, &
mo_coeff
TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_type), POINTER :: xc_pw_pool
TYPE(pw_type) :: work_v_rspace
TYPE(pw_type), DIMENSION(:), POINTER :: rho_g, rho_r, tau, vxc, vxc_tau
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_vals_type), POINTER :: input, low_spin_roks_section, xc_section
TYPE(section_vals_type), POINTER :: hfx_section, input, &
low_spin_roks_section, xc_section
TYPE(virial_type), POINTER :: virial

IF (.NOT. dft_control%low_spin_roks) RETURN
NULLIFY (ks_env, rho_ao)

CALL timeset(routineN, handle)

NULLIFY (ks_env, rho_ao)

! Test for not compatible options
IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
CALL cp_abort(__LOCATION__, "GAPW/GAPW_XC are not compatible with low spin ROKS method.")
END IF
IF (dft_control%do_admm) THEN
CALL cp_abort(__LOCATION__, "ADMM not compatible with low spin ROKS method.")
END IF
IF (dft_control%do_admm) THEN
IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
CALL cp_abort(__LOCATION__, "ADMM with XC correction functional "// &
"not compatible with low spin ROKS method.")
END IF
END IF
IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
dft_control%qs_control%xtb) THEN
CALL cp_abort(__LOCATION__, "SE/xTB/DFTB are not compatible with low spin ROKS method.")
END IF

CALL get_qs_env(qs_env, &
ks_env=ks_env, &
mo_derivs=mo_derivs, &
Expand All @@ -199,6 +227,7 @@ SUBROUTINE low_spin_roks(energy, qs_env, dft_control, just_energy, &

compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")

! some assumptions need to be checked
! we have two spins
Expand All @@ -209,6 +238,9 @@ SUBROUTINE low_spin_roks(energy, qs_env, dft_control, just_energy, &
CPASSERT(uniform_occupation)
CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
CPASSERT(uniform_occupation)
IF (do_hfx .AND. calculate_forces .AND. compute_virial) THEN
CALL cp_abort(__LOCATION__, "ROKS virial with HFX not available.")
END IF

NULLIFY (dbcsr_deriv)
CALL dbcsr_init_p(dbcsr_deriv)
Expand Down Expand Up @@ -267,6 +299,16 @@ SUBROUTINE low_spin_roks(energy, qs_env, dft_control, just_energy, &
CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
END DO

IF (do_hfx) THEN
NULLIFY (matrix_hfx)
CALL dbcsr_allocate_matrix_set(matrix_hfx, Nspin)
DO ispin = 1, Nspin
ALLOCATE (matrix_hfx(ispin)%matrix)
CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
name="HFX matrix low spin roks")
END DO
END IF

! grids in real and g space for rho and vxc
! tau functionals are not supported
NULLIFY (tau, vxc_tau, vxc)
Expand Down Expand Up @@ -328,21 +370,57 @@ SUBROUTINE low_spin_roks(energy, qs_env, dft_control, just_energy, &

energy%exc = energy%exc + energy_scaling(iterm)*exc

IF (do_hfx) THEN
! Add Hartree-Fock contribution
DO ispin = 1, Nspin
CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
END DO
ehfx = energy%ex
CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, &
recalc_integrals=.FALSE., update_energy=.TRUE.)
energy%ex = ehfx + energy_scaling(iterm)*energy%ex
END IF

! add the corresponding derivatives to the MO derivatives
IF (.NOT. just_energy) THEN
! get the potential in matrix form
DO ispin = 1, Nspin
CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
! use a work_v_rspace
work_v_rspace%cr3d = (energy_scaling(iterm)*vxc(ispin)%pw_grid%dvol)* &
vxc(ispin)%cr3d
! zero first ?!
CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
qs_env=qs_env, calculate_forces=calculate_forces)
CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc(ispin))
END DO
DEALLOCATE (vxc)

IF (do_hfx) THEN
! add HFX contribution
DO ispin = 1, Nspin
CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
1.0_dp, energy_scaling(iterm))
END DO
IF (calculate_forces) THEN
CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
IF (x_data(1, 1)%n_rep_hf /= 1) THEN
CALL cp_abort(__LOCATION__, "Multiple HFX section forces not compatible "// &
"with low spin ROKS method.")
END IF
IF (x_data(1, 1)%do_hfx_ri) THEN
CALL cp_abort(__LOCATION__, "HFX_RI forces not compatible with low spin ROKS method.")
ELSE
irep = 1
NULLIFY (mdummy)
matrix_p2(1:Nspin, 1:1) => matrix_p(1:Nspin)
CALL derivatives_four_center(qs_env, matrix_p2, mdummy, hfx_section, para_env, &
irep, compute_virial, &
adiabatic_rescale_factor=energy_scaling(iterm))
END IF
END IF

END IF

! add this to the mo_derivs, again based on the alpha mo_coeff
DO ispin = 1, Nspin
CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
Expand All @@ -366,6 +444,9 @@ SUBROUTINE low_spin_roks(energy, qs_env, dft_control, just_energy, &
DEALLOCATE (rho_r, rho_g)
CALL dbcsr_deallocate_matrix_set(matrix_p)
CALL dbcsr_deallocate_matrix_set(matrix_h)
IF (do_hfx) THEN
CALL dbcsr_deallocate_matrix_set(matrix_hfx)
END IF

CALL pw_pool_give_back_pw(auxbas_pw_pool, work_v_rspace)

Expand Down
7 changes: 4 additions & 3 deletions tests/QS/regtest-lsroks/O2.inp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,10 @@
MULTIP 3
ROKS
&LOW_SPIN_ROKS
ENERGY_SCALING 1.0 -1.0
SPIN_CONFIGURATION 1 1
SPIN_CONFIGURATION 1 2
! Singlet: E(s) = E(t) - 2*E(t) + 2*E(m) = 2*E(m) - E(t)
ENERGY_SCALING -2.0 2.0
SPIN_CONFIGURATION 1 1 ! (t)
SPIN_CONFIGURATION 1 2 ! (m)
&END
&MGRID
CUTOFF 280
Expand Down
54 changes: 54 additions & 0 deletions tests/QS/regtest-lsroks/O2_pbe0.inp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,54 @@
&GLOBAL
PROJECT O2
RUN_TYPE ENERGY
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
POTENTIAL_FILE_NAME GTH_POTENTIALS
MULTIP 3
ROKS
&LOW_SPIN_ROKS
! Singlet: E(s) = E(t) - 2*E(t) + 2*E(m) = 2*E(m) - E(t)
ENERGY_SCALING -2.0 2.0
SPIN_CONFIGURATION 1 1 ! (t)
SPIN_CONFIGURATION 1 2 ! (m)
&END
&MGRID
CUTOFF 280
&END MGRID
&QS
&END QS
&SCF
SCF_GUESS ATOMIC
EPS_SCF 1.0E-6
MAX_SCF 10
&OT
ROTATION
&END
&OUTER_SCF
EPS_SCF 1.0E-6
MAX_SCF 1
&END
&END SCF
&XC
&XC_FUNCTIONAL PBE0
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 4.0 4.0 6.0
&END CELL
&KIND O
BASIS_SET DZVP-GTH
POTENTIAL GTH-PADE-q6
&END KIND
&COORD
O 0.000000 0.000000 0.608000
O 0.000000 0.000000 -0.608000
&END COORD
&END SUBSYS
&END FORCE_EVAL
5 changes: 4 additions & 1 deletion tests/QS/regtest-lsroks/TEST_FILES
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,5 @@
O2.inp 1 2e-13 -31.43642997434278
O2.inp 1 1e-12 -31.38777641185785
O2_pbe0.inp 1 1e-12 -31.93896892726076
ch2o_rs.inp 1 2e-08 -22.5843713873
ch2o_hf.inp 1 2e-08 -22.0854638839
#EOF

0 comments on commit 04096ff

Please sign in to comment.