Add support for a separate XC section within the AS module

This commit is contained in:
Stefano Battaglia 2025-12-01 14:01:10 +01:00 committed by GitHub
parent 26c18786f4
commit 07059287c0
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 302 additions and 78 deletions

View file

@ -159,11 +159,13 @@ CONTAINS
!> \brief ...
!> \param qs_env ...
!> \param calculate_forces ...
!> \param ext_xc_section ...
! **************************************************************************************************
SUBROUTINE hfx_admm_init(qs_env, calculate_forces)
SUBROUTINE hfx_admm_init(qs_env, calculate_forces, ext_xc_section)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_admm_init'
@ -208,6 +210,7 @@ CONTAINS
IF (PRESENT(calculate_forces)) calc_forces = .TRUE.
hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
IF (PRESENT(ext_xc_section)) hfx_sections => section_vals_get_subs_vals(ext_xc_section, "HF")
CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
IF (n_rep_hf > 1) &
@ -220,6 +223,7 @@ CONTAINS
CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
CALL set_qs_env(qs_env, admm_env=admm_env)
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
IF (PRESENT(ext_xc_section)) xc_section => ext_xc_section
CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
admm_env=admm_env)
@ -1128,12 +1132,13 @@ CONTAINS
!> \param just_energy ...
!> \param v_rspace_new ...
!> \param v_tau_rspace ...
!> \param ext_xc_section ...
!> \par History
!> refactoring 03-2011 [MI]
! **************************************************************************************************
SUBROUTINE hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, &
just_energy, v_rspace_new, v_tau_rspace)
just_energy, v_rspace_new, v_tau_rspace, ext_xc_section)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks
@ -1141,6 +1146,7 @@ CONTAINS
TYPE(qs_energy_type), POINTER :: energy
LOGICAL, INTENT(in) :: calculate_forces, just_energy
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_tau_rspace
TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ks_matrix'
@ -1203,6 +1209,8 @@ CONTAINS
IF (use_virial .AND. calculate_forces) virial%pv_fock_4c = 0.0_dp
hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
IF (PRESENT(ext_xc_section)) hfx_sections => section_vals_get_subs_vals(ext_xc_section, "HF")
CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
i_rep_section=1)

View file

@ -25,6 +25,7 @@ MODULE input_cp2k_as
eri_operator_trunc, eri_operator_yukawa, gaussian, manual_selection, mao_projection, &
no_solver, qiskit_solver, wannier_projection
USE input_cp2k_loc, ONLY: create_localize_section
USE input_cp2k_xc, ONLY: create_xc_section
USE input_keyword_types, ONLY: keyword_create,&
keyword_release,&
keyword_type
@ -170,7 +171,7 @@ CONTAINS
CALL keyword_create(keyword, __LOCATION__, name="ALPHA", &
description="Fraction of new density to be mixed with previous one in SCF embedding.", &
usage="ALPHA 0.25", type_of_var=real_t, &
default_r_val=0.4_dp)
default_r_val=0.8_dp)
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)
@ -201,6 +202,10 @@ CONTAINS
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
CALL create_xc_section(subsection)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)
END SUBROUTINE create_active_space_section
! **************************************************************************************************

View file

@ -122,10 +122,12 @@ CONTAINS
section1 => section_vals_get_subs_vals(section, "XAS_TDP%KERNEL")
section2 => section_vals_get_subs_vals(section, "XAS_TDP%KERNEL%XC_FUNCTIONAL")
CALL xc_functionals_expand(section2, section1)
CALL xc_functionals_expand(section2, section1)
section1 => section_vals_get_subs_vals(sections, "PROPERTIES%RIXS%XAS_TDP%KERNEL")
section2 => section_vals_get_subs_vals(sections, "PROPERTIES%RIXS%XAS_TDP%KERNEL%XC_FUNCTIONAL")
CALL xc_functionals_expand(section2, section1)
section1 => section_vals_get_subs_vals(section, "ACTIVE_SPACE%XC")
section2 => section_vals_get_subs_vals(section, "ACTIVE_SPACE%XC%XC_FUNCTIONAL")
CALL xc_functionals_expand(section2, section1)
END DO
CALL timestop(handle)

View file

@ -13,7 +13,8 @@
! **************************************************************************************************
MODULE qs_active_space_methods
USE admm_types, ONLY: admm_type, &
get_admm_env
get_admm_env, &
admm_env_release
USE atomic_kind_types, ONLY: atomic_kind_type
USE basis_set_types, ONLY: allocate_sto_basis_set, &
create_gto_from_sto_basis, &
@ -72,6 +73,7 @@ MODULE qs_active_space_methods
default_string_length, &
dp, &
int_8
USE hfx_types, ONLY: hfx_create, hfx_release
USE machine, ONLY: m_walltime, m_flush
USE mathlib, ONLY: diamat_all
USE mathconstants, ONLY: fourpi, twopi, pi, rootpi
@ -132,9 +134,10 @@ MODULE qs_active_space_methods
set_qs_env
USE qs_integrate_potential, ONLY: integrate_v_rspace
USE qs_kind_types, ONLY: qs_kind_type
USE qs_ks_methods, ONLY: qs_ks_update_qs_env
USE qs_ks_methods, ONLY: qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
evaluate_core_matrix_traces
USE qs_ks_types, ONLY: qs_ks_did_change, &
qs_ks_env_type
qs_ks_env_type, set_ks_env
USE qs_mo_io, ONLY: write_mo_set_to_output_unit
USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
USE qs_mo_types, ONLY: allocate_mo_set, &
@ -148,6 +151,7 @@ MODULE qs_active_space_methods
qs_rho_type
USE qs_subsys_types, ONLY: qs_subsys_get, &
qs_subsys_type
USE qs_scf_post_scf, ONLY: qs_scf_compute_properties
USE scf_control_types, ONLY: scf_control_type
#ifndef __NO_SOCKETS
USE sockets_interface, ONLY: accept_socket, &
@ -212,20 +216,22 @@ CONTAINS
CHARACTER(len=10) :: cshell, lnam(5)
CHARACTER(len=default_path_length) :: qcschema_filename
CHARACTER(LEN=default_string_length) :: basis_type
INTEGER :: as_solver, eri_method, eri_operator, eri_print, group_size, handle, i, iatom, &
ishell, isp, ispin, iw, j, jm, m, max_orb_ind, mselect, n1, n2, nao, natom, nel, &
nelec_active, nelec_inactive, nelec_total, nmo, nmo_active, nmo_available, nmo_inactive, &
nmo_inactive_remaining, nmo_occ, nmo_virtual, nn1, nn2, nrow_global, nspins
INTEGER, DIMENSION(5) :: nshell
INTEGER, DIMENSION(:), POINTER :: invals
LOGICAL :: do_kpoints, ex_omega, ex_operator, &
ex_perd, ex_rcut, explicit, &
stop_after_print, store_wfn
REAL(KIND=dp) :: eri_eps_filter, eri_eps_grid, eri_eps_int, eri_gpw_cutoff, eri_op_omega, &
eri_rcut, eri_rel_cutoff, fel, focc, maxocc, nze_percentage
LOGICAL :: do_ddapc, do_kpoints, ex_omega, &
ex_operator, ex_perd, ex_rcut, &
explicit, stop_after_print, store_wfn
REAL(KIND=dp) :: alpha, eri_eps_filter, eri_eps_grid, eri_eps_int, eri_gpw_cutoff, &
eri_op_omega, eri_rcut, eri_rel_cutoff, fel, focc, maxocc, nze_percentage
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: eigenvalues
REAL(KIND=dp), DIMENSION(:), POINTER :: evals_virtual
TYPE(active_space_type), POINTER :: active_space_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_blacs_env_type), POINTER :: context
TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
@ -240,13 +246,17 @@ CONTAINS
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(mo_set_type), POINTER :: mo_set, mo_set_active, mo_set_inactive
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(preconditioner_type), POINTER :: local_preconditioner
TYPE(qcschema_type) :: qcschema_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(scf_control_type), POINTER :: scf_control
TYPE(section_vals_type), POINTER :: as_input, input, loc_print, loc_section, &
print_orb
TYPE(section_vals_type), POINTER :: adiabatic_rescaling, as_input, &
hfx_section, input, loc_print, &
loc_section, print_orb, xc_section
!--------------------------------------------------------------------------------------------!
@ -262,16 +272,45 @@ CONTAINS
IF (iw > 0) THEN
WRITE (iw, '(/,T2,A)') &
'!-----------------------------------------------------------------------------!'
WRITE (iw, '(T26,A)') "Generate Active Space Hamiltonian"
WRITE (iw, '(T27,A)') "Interface to CAS-CI and DMRG-CI"
WRITE (iw, '(T26,A)') "Active Space Embedding Module"
WRITE (iw, '(T2,A)') &
'!-----------------------------------------------------------------------------!'
END IF
! k-points?
CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
CALL get_qs_env(qs_env, do_kpoints=do_kpoints, dft_control=dft_control)
IF (do_kpoints) THEN
CALL cp_abort(__LOCATION__, "k-points not supported in active space interface")
CALL cp_abort(__LOCATION__, "k-points not supported in active space module")
END IF
! adiabatic rescaling?
adiabatic_rescaling => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
CALL section_vals_get(adiabatic_rescaling, explicit=explicit)
IF (explicit) THEN
CALL cp_abort(__LOCATION__, "Adiabatic rescaling not supported in active space module")
END IF
! Setup the possible usage of DDAPC charges
do_ddapc = dft_control%qs_control%ddapc_restraint .OR. &
qs_env%cp_ddapc_ewald%do_decoupling .OR. &
qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
qs_env%cp_ddapc_ewald%do_solvation
IF (do_ddapc) THEN
CALL cp_abort(__LOCATION__, "DDAPC charges are not supported in the active space module")
END IF
IF (dft_control%do_sccs) THEN
CALL cp_abort(__LOCATION__, "SCCS is not supported in the active space module")
END IF
IF (dft_control%correct_surf_dip) THEN
IF (dft_control%surf_dip_correct_switch) THEN
CALL cp_abort(__LOCATION__, "Surface dipole correction not supported in the AS module")
END IF
END IF
IF (dft_control%smeagol_control%smeagol_enabled) THEN
CALL cp_abort(__LOCATION__, "SMEAGOL is not supported in the active space module")
END IF
IF (dft_control%qs_control%do_kg) THEN
CALL cp_abort(__LOCATION__, "KG correction not supported in the active space module")
END IF
NULLIFY (active_space_env)
@ -305,6 +344,10 @@ CONTAINS
CPABORT("The remaining number of inactive electrons has to be even.")
END IF
CALL section_vals_val_get(as_input, "ALPHA", r_val=alpha)
IF (alpha < 0.0 .OR. alpha > 1.0) CPABORT("Specify a damping factor between 0 and 1.")
active_space_env%alpha = alpha
IF (iw > 0) THEN
WRITE (iw, '(T3,A,T70,I10)') "Total number of electrons", nelec_total
WRITE (iw, '(T3,A,T70,I10)') "Number of inactive electrons", nelec_inactive
@ -572,9 +615,9 @@ CONTAINS
ELSE
IF (dft_control%do_admm) THEN
! For admm_mo, the auxiliary density is computed from the MOs, which never change
! in the rs-dft embedding, therefore the energy is wrong as the LR HFX never changes
! in the rs-dft embedding, therefore the energy is wrong as the LR HFX never changes.
! For admm_dm, the auxiliary density is computed from the density matrix, which is
! updated at each iteration and therefore works!
! updated at each iteration and therefore works.
IF (dft_control%do_admm_mo) THEN
CPABORT("ADMM currently possible only with purification none_dm")
END IF
@ -1002,10 +1045,50 @@ CONTAINS
! compute alpha-beta overlap matrix in case of spin-polarized calculation
CALL calculate_spin_pol_overlap(active_space_env%mos_active, qs_env, active_space_env)
! figure out if we have a new xc section for the AS
xc_section => section_vals_get_subs_vals(input, "DFT%ACTIVE_SPACE%XC")
explicit = .FALSE.
IF (ASSOCIATED(xc_section)) CALL section_vals_get(xc_section, explicit=explicit)
! rebuild KS matrix if needed
IF (explicit) THEN
! release the hfx data if it was part of the SCF functional
IF (ASSOCIATED(qs_env%x_data)) CALL hfx_release(qs_env%x_data)
! also release the admm environment in case we are using admm
IF (ASSOCIATED(qs_env%admm_env)) CALL admm_env_release(qs_env%admm_env)
CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
particle_set=particle_set, cell=cell, ks_env=ks_env)
IF (dft_control%do_admm) THEN
basis_type = 'AUX_FIT'
ELSE
basis_type = 'ORB'
END IF
hfx_section => section_vals_get_subs_vals(xc_section, "HF")
CALL hfx_create(qs_env%x_data, para_env, hfx_section, atomic_kind_set, &
qs_kind_set, particle_set, dft_control, cell, orb_basis=basis_type, &
nelectron_total=nelec_total)
qs_env%requires_matrix_vxc = .TRUE. ! needs to be set only once
! a bit of a hack: this forces a new re-init of HFX
CALL set_ks_env(ks_env, s_mstruct_changed=.TRUE.)
CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., &
just_energy=.FALSE., &
ext_xc_section=xc_section)
! we need to reset it to false
CALL set_ks_env(ks_env, s_mstruct_changed=.FALSE.)
ELSE
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
END IF
! set the xc_section
active_space_env%xc_section => xc_section
CALL get_qs_env(qs_env, energy=energy)
! transform KS/Fock, Vxc and Hcore to AS MO basis
CALL calculate_operators(active_space_env%mos_active, qs_env, active_space_env)
! set the reference energy in the active space
CALL get_qs_env(qs_env, energy=energy)
active_space_env%energy_ref = energy%total
! calculate inactive energy and embedding potential
CALL subspace_fock_matrix(active_space_env)
@ -1022,6 +1105,7 @@ CONTAINS
END IF
CASE (qiskit_solver)
CALL rsdft_embedding(qs_env, active_space_env, as_input)
CALL qs_scf_compute_properties(qs_env, wf_type="MC-DFT", do_mp2=.FALSE.)
CASE DEFAULT
CPABORT("Unknown active space solver")
END SELECT
@ -1046,7 +1130,7 @@ CONTAINS
END SUBROUTINE active_space_main
! **************************************************************************************************
!> \brief computes the one-electron operators in the subspace of the provided orbital set
!> \brief computes the alpha-beta overlap within the active subspace
!> \param mos the molecular orbital set within the active subspace
!> \param qs_env ...
!> \param active_space_env ...
@ -1099,7 +1183,7 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'calculate_operators'
INTEGER :: handle, is, nmo, nspins
INTEGER :: handle, ispin, nmo, nspins
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: h_matrix, ks_matrix
@ -1110,13 +1194,10 @@ CONTAINS
! Kohn-Sham / Fock operator
CALL cp_fm_release(active_space_env%ks_sub)
CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix)
IF (SIZE(ks_matrix, 2) > 1) THEN
CPABORT("No k-points allowed at this point")
END IF
ALLOCATE (active_space_env%ks_sub(nspins))
DO is = 1, nspins
CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
CALL subspace_operator(mo_coeff, nmo, ks_matrix(is, 1)%matrix, active_space_env%ks_sub(is))
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
CALL subspace_operator(mo_coeff, nmo, ks_matrix(ispin, 1)%matrix, active_space_env%ks_sub(ispin))
END DO
! Core Hamiltonian
@ -1124,13 +1205,10 @@ CONTAINS
NULLIFY (h_matrix)
CALL get_qs_env(qs_env=qs_env, matrix_h_kp=h_matrix)
IF (SIZE(h_matrix, 2) > 1) THEN
CPABORT("No k-points allowed at this point")
END IF
ALLOCATE (active_space_env%h_sub(nspins))
DO is = 1, nspins
CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
CALL subspace_operator(mo_coeff, nmo, h_matrix(1, 1)%matrix, active_space_env%h_sub(is))
DO ispin = 1, nspins
CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
CALL subspace_operator(mo_coeff, nmo, h_matrix(1, 1)%matrix, active_space_env%h_sub(ispin))
END DO
CALL timestop(handle)
@ -1139,21 +1217,21 @@ CONTAINS
! **************************************************************************************************
!> \brief computes a one-electron operator in the subspace of the provided orbital set
!> \param orbitals the orbital coefficient matrix
!> \param nmo the number of orbitals
!> \param mo_coeff the orbital coefficient matrix
!> \param nmo the number of subspace orbitals
!> \param op_matrix operator matrix in AO basis
!> \param op_sub operator in orbital basis
!> \param orbitals_b ...
!> \param mo_coeff_b the beta orbital coefficients
!> \par History
!> 04.2016 created [JGH]
! **************************************************************************************************
SUBROUTINE subspace_operator(orbitals, nmo, op_matrix, op_sub, orbitals_b)
SUBROUTINE subspace_operator(mo_coeff, nmo, op_matrix, op_sub, mo_coeff_b)
TYPE(cp_fm_type), INTENT(IN) :: orbitals
TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
INTEGER, INTENT(IN) :: nmo
TYPE(dbcsr_type), POINTER :: op_matrix
TYPE(cp_fm_type), INTENT(INOUT) :: op_sub
TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: orbitals_b
TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: mo_coeff_b
CHARACTER(len=*), PARAMETER :: routineN = 'subspace_operator'
@ -1162,22 +1240,22 @@ CONTAINS
CALL timeset(routineN, handle)
CALL cp_fm_get_info(matrix=orbitals, ncol_global=ncol, nrow_global=nrow)
CALL cp_fm_get_info(matrix=mo_coeff, ncol_global=ncol, nrow_global=nrow)
CPASSERT(nmo <= ncol)
IF (nmo > 0) THEN
CALL cp_fm_create(vectors, orbitals%matrix_struct, "vectors")
CALL create_subspace_matrix(orbitals, op_sub, nmo)
CALL cp_fm_create(vectors, mo_coeff%matrix_struct, "vectors")
CALL create_subspace_matrix(mo_coeff, op_sub, nmo)
IF (PRESENT(orbitals_b)) THEN
IF (PRESENT(mo_coeff_b)) THEN
! if beta orbitals are present, compute the cross alpha_beta term
CALL cp_dbcsr_sm_fm_multiply(op_matrix, orbitals_b, vectors, nmo)
CALL cp_dbcsr_sm_fm_multiply(op_matrix, mo_coeff_b, vectors, nmo)
ELSE
! otherwise the same spin, whatever that is
CALL cp_dbcsr_sm_fm_multiply(op_matrix, orbitals, vectors, nmo)
CALL cp_dbcsr_sm_fm_multiply(op_matrix, mo_coeff, vectors, nmo)
END IF
CALL parallel_gemm('T', 'N', nmo, nmo, nrow, 1.0_dp, orbitals, vectors, 0.0_dp, op_sub)
CALL parallel_gemm('T', 'N', nmo, nmo, nrow, 1.0_dp, mo_coeff, vectors, 0.0_dp, op_sub)
CALL cp_fm_release(vectors)
END IF
@ -1484,7 +1562,7 @@ CONTAINS
END DO
IF (print1 .AND. iw > 0) THEN
rsize = rsize*8._dp/1000000._dp
WRITE (iw, "(T4,'ERI_GPW|',' Store active orbitals on real space grid ',T63,F12.3,' MB')") rsize
WRITE (iw, "(T2,'ERI_GPW|',' Store active orbitals on real space grid ',T66,F12.3,' MB')") rsize
END IF
ALLOCATE (wfn_a(nmo, nspins))
DO ispin = 1, nspins
@ -1495,7 +1573,7 @@ CONTAINS
CALL calculate_wavefunction(mo_coeff, iwfn, wfn_a(iwfn, ispin), rho_g, atomic_kind_set, &
qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
IF (print2 .AND. iw > 0) THEN
WRITE (iw, "(T4,'ERI_GPW|',' Orbital stored ',I4,' Spin ',i1)") iwfn, ispin
WRITE (iw, "(T2,'ERI_GPW|',' Orbital stored ',I4,' Spin ',i1)") iwfn, ispin
END IF
END DO
END DO
@ -1693,7 +1771,7 @@ CONTAINS
END DO
IF (print1 .AND. iw > 0) THEN
WRITE (iw, "(T4,'ERI_GPW|',' Number of Integrals stored ',T68,I10)") stored_integrals
WRITE (iw, "(T2,'ERI_GPW|',' Number of Integrals stored ',T71,I10)") stored_integrals
END IF
IF (eri_env%eri_gpw%store_wfn) THEN
@ -2870,7 +2948,7 @@ CONTAINS
CALL diamat_all(pmat, noon(:, ispin))
IF (iw > 0) THEN
WRITE (iw, '(T3,A,I2,A)') "Natural orbitals occupation numbers for spin ", ispin
WRITE (iw, '(/,T3,A,I2,A)') "Natural orbitals occupation numbers for spin ", ispin
DO i1 = 1, nmo_active, 8
jm = MIN(7, nmo_active - i1)
! noons are stored in ascending order, so reverse-print them
@ -2909,15 +2987,15 @@ CONTAINS
INTEGER :: iw, client_fd, socket_fd, iter, max_iter
LOGICAL :: converged, do_scf_embedding, ionode
REAL(KIND=dp) :: alpha, delta_E, energy_corr, energy_new, &
REAL(KIND=dp) :: delta_E, energy_corr, energy_new, &
energy_old, energy_scf, eps_iter, t1, t2
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(dft_control_type), POINTER :: dft_control
CALL timeset(routineN, handle)
@ -2936,10 +3014,6 @@ CONTAINS
IF (max_iter < 0) CPABORT("Specify a non-negative number of max iterations.")
CALL section_vals_val_get(as_input, "EPS_ITER", r_val=eps_iter)
IF (eps_iter < 0.0) CPABORT("Specify a non-negative convergence threshold.")
CALL section_vals_val_get(as_input, "ALPHA", r_val=alpha)
! alpha = 0.0
IF (alpha < 0.0 .OR. alpha > 1.0) CPABORT("Specify a damping factor between 0 and 1.")
active_space_env%alpha = alpha
! create the socket and wait for the client to connect
CALL initialize_socket(socket_fd, client_fd, as_input, ionode)
@ -2949,15 +3023,16 @@ CONTAINS
CALL send_eri_to_client(client_fd, active_space_env, para_env)
! get pointer to density in ao basis
CALL get_qs_env(qs_env, rho=rho, energy=energy, dft_control=dft_control)
CALL get_qs_env(qs_env, rho=rho, energy=energy, ks_env=ks_env)
CALL qs_rho_get(rho, rho_ao=rho_ao)
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,T2,A,/)") "RS-DFT SELF-CONSISTENT OPTIMIZATION"
WRITE (UNIT=iw, FMT="(/,T2,A,/)") &
"RANGE-SEPARATED DFT EMBEDDING SELF-CONSISTENT OPTIMIZATION"
WRITE (iw, '(T3,A,T68,I12)') "Max. iterations", max_iter
WRITE (iw, '(T3,A,T68,E12.4)') "Conv. threshold", eps_iter
WRITE (iw, '(T3,A,T66,F14.2)') "Density damping", alpha
WRITE (iw, '(T3,A,T66,F14.2)') "Density damping", active_space_env%alpha
WRITE (UNIT=iw, FMT="(/,T3,A,T11,A,T21,A,T34,A,T55,A,T75,A,/,T3,A)") &
"Iter", "Update", "Time", "Corr. energy", "Total energy", "Change", REPEAT("-", 78)
@ -3001,10 +3076,15 @@ CONTAINS
CALL update_density_ao(active_space_env, rho_ao) ! rho_ao is updated
! calculate F_ks in AO basis (which contains Vxc) with the new density
qs_env%requires_matrix_vxc = .TRUE.
CALL qs_rho_update_rho(rho, qs_env=qs_env) ! updates rho_r and rho_g using rho_ao
CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) ! set flags about the change
CALL qs_ks_update_qs_env(qs_env) ! this call actually calculates F_ks
! Re-evaluate the traces between the density matrix and the core Hamiltonians
CALL evaluate_core_matrix_traces(qs_env)
! the ks matrix will be rebuilt so this is fine now
! CALL set_ks_env(qs_env%ks_env, potential_changed=.FALSE.)
CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., &
just_energy=.FALSE., &
ext_xc_section=active_space_env%xc_section)
! update the reference energy
active_space_env%energy_ref = energy%total
@ -3044,6 +3124,20 @@ CONTAINS
! update qs total energy to the final rs-DFT energy
energy%total = active_space_env%energy_total
! print final energy contributions
IF (iw > 0) THEN
WRITE (UNIT=iw, FMT="(/,T3,A)") &
"Final energy contributions:"
WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
"Inactive energy:", active_space_env%energy_inactive
WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
"Active energy:", active_space_env%energy_active
WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
"Correlation energy:", energy_corr
WRITE (UNIT=iw, FMT="(T6,A,T56,F20.10)") &
"Total rs-DFT energy:", active_space_env%energy_total
END IF
! print the AS rdm and the natural orbital occupation numbers
CALL print_pmat_noon(active_space_env, iw)

View file

@ -20,6 +20,7 @@ MODULE qs_active_space_types
USE cp_fm_types, ONLY: cp_fm_release,&
cp_fm_type
USE input_constants, ONLY: eri_method_gpw_ht
USE input_section_types, ONLY: section_vals_type
USE kinds, ONLY: default_path_length,&
dp
USE message_passing, ONLY: mp_comm_type
@ -109,6 +110,7 @@ MODULE qs_active_space_types
TYPE(cp_fm_type), DIMENSION(:), POINTER :: fock_sub => NULL()
TYPE(cp_fm_type), DIMENSION(:), POINTER :: sab_sub => NULL()
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmat_inactive => NULL()
TYPE(section_vals_type), POINTER :: xc_section => NULL()
END TYPE active_space_type
ABSTRACT INTERFACE

View file

@ -59,8 +59,7 @@ MODULE qs_ks_methods
smeagol_runtype_emtransport
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
section_vals_type
USE kg_correction, ONLY: kg_ekin_subset
USE kinds, ONLY: default_string_length,&
dp
@ -146,7 +145,7 @@ MODULE qs_ks_methods
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_methods'
PUBLIC :: calc_rho_tot_gspace, qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
qs_ks_allocate_basics, evaluate_core_matrix_p_mix_new
qs_ks_allocate_basics, evaluate_core_matrix_p_mix_new, evaluate_core_matrix_traces
CONTAINS
@ -160,6 +159,7 @@ CONTAINS
!> ks matrix. Defaults to false
!> \param print_active ...
!> \param ext_ks_matrix ...
!> \param ext_xc_section ...
!> \par History
!> 06.2002 moved from qs_scf to qs_ks_methods, use of ks_env
!> new did_change scheme [fawzi]
@ -175,12 +175,13 @@ CONTAINS
!> to qs_env components?
! **************************************************************************************************
SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
print_active, ext_ks_matrix)
print_active, ext_ks_matrix, ext_xc_section)
TYPE(qs_environment_type), POINTER :: qs_env
LOGICAL, INTENT(in) :: calculate_forces, just_energy
LOGICAL, INTENT(IN), OPTIONAL :: print_active
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
POINTER :: ext_ks_matrix
TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_build_kohn_sham_matrix'
@ -188,7 +189,7 @@ CONTAINS
INTEGER :: handle, iatom, img, ispin, nimages, &
nspins
LOGICAL :: do_adiabatic_rescaling, do_ddapc, do_hfx, do_ppl, dokp, gapw, gapw_xc, &
hfx_treat_lsd_in_core, just_energy_xc, lrigpw, my_print, rigpw, use_virial
just_energy_xc, lrigpw, my_print, rigpw, use_virial
REAL(KIND=dp) :: ecore_ppl, edisp, ee_ener, ekin_mol, &
mulliken_order_p, vscale
REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc
@ -272,12 +273,6 @@ CONTAINS
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
CALL section_vals_get(hfx_sections, explicit=do_hfx)
IF (do_hfx) THEN
CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
i_rep_section=1)
END IF
adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
just_energy_xc = just_energy
@ -509,7 +504,11 @@ CONTAINS
! calculate the density matrix for the fitted mo_coeffs
IF (dft_control%do_admm) THEN
CALL hfx_admm_init(qs_env, calculate_forces)
IF (PRESENT(ext_xc_section)) THEN
CALL hfx_admm_init(qs_env, calculate_forces, ext_xc_section)
ELSE
CALL hfx_admm_init(qs_env, calculate_forces)
END IF
IF (dft_control%do_admm_mo) THEN
IF (qs_env%run_rtp) THEN
@ -565,6 +564,8 @@ CONTAINS
xc_section => admm_env%xc_section_primary
ELSE
xc_section => section_vals_get_subs_vals(input, "DFT%XC")
! build ks matrix with an xc section potentially different from the one defined in input
IF (PRESENT(ext_xc_section)) xc_section => ext_xc_section
END IF
IF (gapw_xc) THEN
@ -621,14 +622,16 @@ CONTAINS
END IF
! *** Add Hartree-Fock contribution if required ***
hfx_sections => section_vals_get_subs_vals(xc_section, "HF")
CALL section_vals_get(hfx_sections, explicit=do_hfx)
IF (do_hfx) THEN
IF (dokp) THEN
CALL hfx_ks_matrix_kp(qs_env, ks_matrix, energy, calculate_forces)
ELSE
! ext_xc_section may contain a hfx section
CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, &
just_energy, v_rspace_new, v_tau_rspace)
just_energy, v_rspace_new, v_tau_rspace, ext_xc_section=xc_section)
END IF
!! Adiabatic rescaling only if do_hfx; right?????
END IF !do_hfx
IF (do_ppl .AND. calculate_forces) THEN
@ -745,7 +748,7 @@ CONTAINS
ELSE
IF (do_hfx) THEN
IF (.FALSE.) THEN
CPWARN("KS matrix not longer correct. Check possible problems with property calculations!")
CPWARN("KS matrix no longer correct. Check possible problems with property calculations!")
END IF
END IF
END IF ! .NOT. just energy

View file

@ -28,4 +28,5 @@
#h2_manual_UKS_1+2.inp 92 1e-8 8.52511231
"h2_manual_UKS_1+2.inp" = [{matcher="M092", tol=1e-8, ref=8.52386602}]
"h2_lr_trunc.inp" = [{matcher="M092", tol=1e-8, ref=6.67878293}]
"h2_as_xc_section.inp" = [{matcher="M092", tol=1e-8, ref=6.69976941}]
#EOF

View file

@ -0,0 +1,109 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT h2
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
&ACTIVE_SPACE
ACTIVE_ELECTRONS 2
ACTIVE_ORBITALS 4
&ERI
CUTOFF_RADIUS 4.5
EPS_INTEGRAL 1E-10
METHOD FULL_GPW
OMEGA 0.5
OPERATOR LR_TRUNC
PERIODICITY 1 1 1
&END ERI
&ERI_GPW
CUTOFF 250
&END ERI_GPW
&FCIDUMP
FILENAME __STD_OUT__
&END FCIDUMP
&XC
DENSITY_CUTOFF 1E-008
GRADIENT_CUTOFF 1E-008
&HF
FRACTION 1.0
&INTERACTION_POTENTIAL
CUTOFF_RADIUS 4.5
OMEGA 0.5
POTENTIAL_TYPE MIX_CL_TRUNC
SCALE_COULOMB 0.25
SCALE_LONGRANGE 0.75
&END INTERACTION_POTENTIAL
&END HF
&XC_FUNCTIONAL
&GGA_C_PBE_ERF_GWS
_OMEGA 0.5
&END GGA_C_PBE_ERF_GWS
&GGA_X_PBE_ERF_GWS
_OMEGA 0.5
SCALE 0.75
&END GGA_X_PBE_ERF_GWS
&END XC_FUNCTIONAL
&END XC
&END ACTIVE_SPACE
&MGRID
CUTOFF 300
&END MGRID
&POISSON
PERIODIC XYZ
POISSON_SOLVER PERIODIC
&END POISSON
&QS
METHOD GAPW
&END QS
&SCF
ADDED_MOS 3
EPS_DIIS 1E-003
EPS_EIGVAL 1E-009
EPS_SCF 1E-009
MAX_SCF 20
&PRINT
&RESTART OFF
&END RESTART
&END PRINT
&END SCF
&XC
DENSITY_CUTOFF 1E-008
GRADIENT_CUTOFF 1E-008
&HF
FRACTION 0.25
&INTERACTION_POTENTIAL
OMEGA 0.5
POTENTIAL_TYPE SHORTRANGE
&END INTERACTION_POTENTIAL
&END HF
&XC_FUNCTIONAL
&GGA_C_PBE_ERF_GWS
_OMEGA 0.5
&END GGA_C_PBE_ERF_GWS
&GGA_X_PBE
&END GGA_X_PBE
&GGA_X_PBE_ERF_GWS
_OMEGA 0.5
SCALE -0.25
&END GGA_X_PBE_ERF_GWS
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC XYZ
&END CELL
&COORD
H 0.0 0.0 5.356
H 0.0 0.0 4.644
&END COORD
&KIND H
BASIS_SET 6-31G*
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL