comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Low spin ROKS with hybrid functionals (no ADMM) + regtests (#2582)

Browse source
This commit is contained in:
Juerg Hutter 2023-02-17 12:33:58 +01:00 • committed by GitHub
parent 8c73008cef
commit 04096ff215
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
8 changed files with 340 additions and 15 deletions
Download patch file Download diff file Expand all files Collapse all files

14
src/qs_environment.F
View file

@ -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

2
src/qs_ks_methods.F
View file

@ -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

99
src/qs_ks_utils.F
View file

@ -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)

7
tests/QS/regtest-lsroks/O2.inp
View file

@ -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

54
tests/QS/regtest-lsroks/O2_pbe0.inp Normal file
View file

@ -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
tests/QS/regtest-lsroks/TEST_FILES
View file

@ -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

92
tests/QS/regtest-lsroks/ch2o_hf.inp Normal file
View file

@ -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

82
tests/QS/regtest-lsroks/ch2o_rs.inp Normal file
View file

@ -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