From 04096ff2157fb22a84c717b3717398fe2432229b Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Fri, 17 Feb 2023 12:33:58 +0100 Subject: [PATCH] Low spin ROKS with hybrid functionals (no ADMM) + regtests (#2582) --- src/qs_environment.F | 14 +++- src/qs_ks_methods.F | 2 +- src/qs_ks_utils.F | 99 ++++++++++++++++++++++++++--- tests/QS/regtest-lsroks/O2.inp | 7 +- tests/QS/regtest-lsroks/O2_pbe0.inp | 54 ++++++++++++++++ tests/QS/regtest-lsroks/TEST_FILES | 5 +- tests/QS/regtest-lsroks/ch2o_hf.inp | 92 +++++++++++++++++++++++++++ tests/QS/regtest-lsroks/ch2o_rs.inp | 82 ++++++++++++++++++++++++ 8 files changed, 340 insertions(+), 15 deletions(-) create mode 100644 tests/QS/regtest-lsroks/O2_pbe0.inp create mode 100644 tests/QS/regtest-lsroks/ch2o_hf.inp create mode 100644 tests/QS/regtest-lsroks/ch2o_rs.inp diff --git a/src/qs_environment.F b/src/qs_environment.F index ceda228819..4e8f503298 100644 --- a/src/qs_environment.F +++ b/src/qs_environment.F @@ -1544,6 +1544,18 @@ CONTAINS "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 @@ -1576,7 +1588,7 @@ CONTAINS 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 diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index d83d89faad..43b2a191ad 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -850,7 +850,7 @@ CONTAINS 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 diff --git a/src/qs_ks_utils.F b/src/qs_ks_utils.F index 9c07b59aae..0495f4ef28 100644 --- a/src/qs_ks_utils.F +++ b/src/qs_ks_utils.F @@ -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, & @@ -142,34 +145,39 @@ CONTAINS !> \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 @@ -177,14 +185,34 @@ CONTAINS 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, & @@ -199,6 +227,7 @@ CONTAINS 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 @@ -209,6 +238,9 @@ CONTAINS 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) @@ -267,6 +299,16 @@ CONTAINS 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) @@ -328,21 +370,57 @@ CONTAINS 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, & @@ -366,6 +444,9 @@ CONTAINS 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) diff --git a/tests/QS/regtest-lsroks/O2.inp b/tests/QS/regtest-lsroks/O2.inp index 8c7635b56e..4b34962df7 100644 --- a/tests/QS/regtest-lsroks/O2.inp +++ b/tests/QS/regtest-lsroks/O2.inp @@ -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 diff --git a/tests/QS/regtest-lsroks/O2_pbe0.inp b/tests/QS/regtest-lsroks/O2_pbe0.inp new file mode 100644 index 0000000000..67e9196360 --- /dev/null +++ b/tests/QS/regtest-lsroks/O2_pbe0.inp @@ -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 diff --git a/tests/QS/regtest-lsroks/TEST_FILES b/tests/QS/regtest-lsroks/TEST_FILES index 8ec173767d..fda84778fe 100644 --- a/tests/QS/regtest-lsroks/TEST_FILES +++ b/tests/QS/regtest-lsroks/TEST_FILES @@ -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 diff --git a/tests/QS/regtest-lsroks/ch2o_hf.inp b/tests/QS/regtest-lsroks/ch2o_hf.inp new file mode 100644 index 0000000000..8637254bd4 --- /dev/null +++ b/tests/QS/regtest-lsroks/ch2o_hf.inp @@ -0,0 +1,92 @@ +&GLOBAL + PROJECT ch2o + RUN_TYPE geo_opt + PRINT_LEVEL low +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END + +&FORCE_EVAL + METHOD Quickstep + + &DFT + BASIS_SET_FILE_NAME BASIS_SET + BASIS_SET_FILE_NAME BASIS_ADMM + POTENTIAL_FILE_NAME POTENTIAL + + ROKS + MULTIPLICITY 3 + &LOW_SPIN_ROKS + ENERGY_SCALING -2.0 2.0 + SPIN_CONFIGURATION 1 1 + SPIN_CONFIGURATION 1 2 + &END + + &MGRID + CUTOFF 200 + &END MGRID + + &QS + METHOD GPW + &END QS + + &SCF + SCF_GUESS ATOMIC + EPS_SCF 5.0E-6 + MAX_SCF 20 + &OT + MINIMIZER DIIS + PRECONDITIONER FULL_SINGLE_INVERSE + STEPSIZE 0.1 + ROTATION + &END + &OUTER_SCF + EPS_SCF 5.0E-6 + MAX_SCF 5 + &END + &END SCF + + &XC +# &XC_FUNCTIONAL PBE0 +# &END XC_FUNCTIONAL + &XC_FUNCTIONAL NONE + &END XC_FUNCTIONAL + &HF + &END HF + &END XC + &POISSON + POISSON_SOLVER MT + PERIODIC NONE + &END POISSON + &END DFT + + &SUBSYS + &CELL + ABC 4.0 4.0 4.0 + PERIODIC NONE + &END CELL + + &COORD + O 0.9588431900 1.1234806613 1.8643358699 + C 1.0045827842 1.0372747429 0.6713062328 + H 1.0304990091 1.9328340670 0.0202209074 + H 1.0234151411 0.0574091066 0.1554705073 + &END COORD + + &KIND O + BASIS_SET cFIT3 + POTENTIAL GTH-PBE-q6 + &END KIND + &KIND C + BASIS_SET cFIT3 + POTENTIAL GTH-PBE-q4 + &END KIND + &KIND H + BASIS_SET cFIT3 + POTENTIAL GTH-PBE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-lsroks/ch2o_rs.inp b/tests/QS/regtest-lsroks/ch2o_rs.inp new file mode 100644 index 0000000000..16c96b6549 --- /dev/null +++ b/tests/QS/regtest-lsroks/ch2o_rs.inp @@ -0,0 +1,82 @@ +&GLOBAL + PROJECT ch2o + RUN_TYPE geo_opt + PRINT_LEVEL low +&END GLOBAL +&MOTION + &GEO_OPT + MAX_ITER 1 + &END +&END + +&FORCE_EVAL + METHOD Quickstep + + &DFT + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME POTENTIAL + + ROKS + MULTIPLICITY 3 + &LOW_SPIN_ROKS + ENERGY_SCALING -2.0 2.0 + SPIN_CONFIGURATION 1 1 + SPIN_CONFIGURATION 1 2 + &END + + &MGRID + CUTOFF 200 + &END MGRID + + &QS + METHOD GPW + &END QS + + &SCF + SCF_GUESS ATOMIC + EPS_SCF 1.0E-4 + MAX_SCF 20 + &OT + MINIMIZER DIIS + PRECONDITIONER FULL_SINGLE_INVERSE + STEPSIZE 0.1 + ROTATION + &END + &OUTER_SCF + EPS_SCF 1.0E-4 + MAX_SCF 5 + &END + &END SCF + + &XC + &XC_FUNCTIONAL PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + + &COORD + O 0.9588431900 1.1234806613 1.8643358699 + C 1.0045827842 1.0372747429 0.6713062328 + H 1.0304990091 1.9328340670 0.0202209074 + H 1.0234151411 0.0574091066 0.1554705073 + &END COORD + + &KIND O + BASIS_SET DZVP-GTH-PBE + POTENTIAL GTH-PBE-q6 + &END KIND + &KIND C + BASIS_SET DZVP-GTH-PBE + POTENTIAL GTH-PBE-q4 + &END KIND + &KIND H + BASIS_SET DZV-GTH-PBE + POTENTIAL GTH-PBE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL