diff --git a/src/hfx_admm_utils.F b/src/hfx_admm_utils.F index 999fa728bc..5c2f06e503 100644 --- a/src/hfx_admm_utils.F +++ b/src/hfx_admm_utils.F @@ -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) diff --git a/src/input_cp2k_as.F b/src/input_cp2k_as.F index 70651aa91c..6abb1ded91 100644 --- a/src/input_cp2k_as.F +++ b/src/input_cp2k_as.F @@ -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 ! ************************************************************************************************** diff --git a/src/input_cp2k_check.F b/src/input_cp2k_check.F index ef90ebec53..8103ca1305 100644 --- a/src/input_cp2k_check.F +++ b/src/input_cp2k_check.F @@ -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) diff --git a/src/qs_active_space_methods.F b/src/qs_active_space_methods.F index 7cbf1814b1..a201c69d8e 100644 --- a/src/qs_active_space_methods.F +++ b/src/qs_active_space_methods.F @@ -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) diff --git a/src/qs_active_space_types.F b/src/qs_active_space_types.F index cf683b554a..fdc4a4b552 100644 --- a/src/qs_active_space_types.F +++ b/src/qs_active_space_types.F @@ -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 diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index febf57e2ab..28ea4d672e 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -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 diff --git a/tests/QS/regtest-as-dft/TEST_FILES.toml b/tests/QS/regtest-as-dft/TEST_FILES.toml index 2b97ddc3d1..fa31065e69 100644 --- a/tests/QS/regtest-as-dft/TEST_FILES.toml +++ b/tests/QS/regtest-as-dft/TEST_FILES.toml @@ -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 diff --git a/tests/QS/regtest-as-dft/h2_as_xc_section.inp b/tests/QS/regtest-as-dft/h2_as_xc_section.inp new file mode 100644 index 0000000000..c2077553fa --- /dev/null +++ b/tests/QS/regtest-as-dft/h2_as_xc_section.inp @@ -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