From 69e24ce58a054f0f96f90f31367ce63ed6de5b55 Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Mon, 14 Oct 2024 09:33:58 +0200 Subject: [PATCH] Framework to calculate response forces for external energy expressions based on KS orbitals (#3721) --- src/CMakeLists.txt | 2 +- src/ec_env_types.F | 9 + src/ec_environment.F | 23 +- src/ec_external.F | 813 ++++++++++++++++++++++++++++++++ src/energy_corrections.F | 44 +- src/input_constants.F | 3 +- src/input_cp2k_ec.F | 27 +- src/qs_basis_gradient.F | 4 +- src/qs_energy.F | 4 +- src/qs_energy_matrix_w.F | 172 ------- src/qs_local_properties.F | 4 +- src/qs_matrix_w.F | 8 +- src/qs_tddfpt2_fprint.F | 4 +- src/response_solver.F | 55 ++- src/rpa_im_time_force_methods.F | 4 +- tests/QS/regtest-ec/EXT1.inp | 63 +++ tests/QS/regtest-ec/EXT2.inp | 60 +++ tests/QS/regtest-ec/TEST_FILES | 3 + tests/TEST_TYPES | 3 +- 19 files changed, 1067 insertions(+), 238 deletions(-) create mode 100644 src/ec_external.F delete mode 100644 src/qs_energy_matrix_w.F create mode 100644 tests/QS/regtest-ec/EXT1.inp create mode 100644 tests/QS/regtest-ec/EXT2.inp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d70e2e5cdf..52df62caae 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -134,6 +134,7 @@ list( ec_efield_local.F ec_environment.F ec_env_types.F + ec_external.F ec_methods.F ec_orth_solver.F ed_analysis.F @@ -546,7 +547,6 @@ list( qs_elf_methods.F qs_energy.F qs_energy_init.F - qs_energy_matrix_w.F qs_energy_types.F qs_energy_utils.F qs_energy_window.F diff --git a/src/ec_env_types.F b/src/ec_env_types.F index e3e44bd816..ee604e1489 100644 --- a/src/ec_env_types.F +++ b/src/ec_env_types.F @@ -14,6 +14,8 @@ MODULE ec_env_types USE cp_dbcsr_api, ONLY: dbcsr_p_type USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set + USE cp_fm_types, ONLY: cp_fm_release,& + cp_fm_type USE dm_ls_scf_types, ONLY: ls_scf_env_type,& ls_scf_release USE hfx_types, ONLY: hfx_release,& @@ -66,6 +68,7 @@ MODULE ec_env_types ! debug LOGICAL :: debug_forces = .FALSE. LOGICAL :: debug_stress = .FALSE. + LOGICAL :: debug_external = .FALSE. ! basis set CHARACTER(len=20) :: basis = "" LOGICAL :: mao = .FALSE. @@ -101,6 +104,9 @@ MODULE ec_env_types matrix_w => Null() ! reduce basis TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef => Null() + ! external energy calclulation + TYPE(cp_fm_type), DIMENSION(:), POINTER :: mo_occ => NULL() + TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos => NULL() ! CP equations TYPE(qs_p_env_type), POINTER :: p_env => Null() TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz => Null(), matrix_wz => Null(), & @@ -206,6 +212,9 @@ SUBROUTINE ec_env_release(ec_env) IF (ASSOCIATED(ec_env%xc_section_aux)) CALL section_vals_release(ec_env%xc_section_aux) IF (ASSOCIATED(ec_env%xc_section_primary)) CALL section_vals_release(ec_env%xc_section_primary) + CALL cp_fm_release(ec_env%mo_occ) + CALL cp_fm_release(ec_env%cpmos) + DEALLOCATE (ec_env) END IF diff --git a/src/ec_environment.F b/src/ec_environment.F index 0cbcceb799..0d1ed263dc 100644 --- a/src/ec_environment.F +++ b/src/ec_environment.F @@ -28,12 +28,12 @@ MODULE ec_environment USE dm_ls_scf_types, ONLY: ls_scf_env_type USE ec_env_types, ONLY: energy_correction_type USE input_constants, ONLY: & - ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, & - ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, kg_cholesky, ls_cluster_atomic, & - ls_cluster_molecular, ls_s_inversion_hotelling, ls_s_inversion_none, & - ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, & - ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, xc_vdw_fun_nonloc, & - xc_vdw_fun_pairpot + ec_diagonalization, ec_functional_dc, ec_functional_ext, ec_functional_harris, & + ec_matrix_sign, ec_matrix_tc2, ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, & + kg_cholesky, ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, & + ls_s_inversion_none, ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, & + ls_s_preconditioner_molecular, ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, & + xc_vdw_fun_nonloc, xc_vdw_fun_pairpot USE input_cp2k_check, ONLY: xc_functionals_expand USE input_section_types, ONLY: section_get_ival,& section_vals_get,& @@ -198,6 +198,8 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section) l_val=ec_env%debug_forces) CALL section_vals_val_get(ec_section, "DEBUG_STRESS", & l_val=ec_env%debug_stress) + CALL section_vals_val_get(ec_section, "DEBUG_EXTERNAL_METHOD", & + l_val=ec_env%debug_external) ! ADMM CALL section_vals_val_get(ec_section, "ADMM", l_val=ec_env%do_ec_admm) @@ -281,6 +283,11 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section) "DC-DFT: Correction and ground state need to use the same basis. "// & "Checked by comparing basis set names only.") END IF + IF (ec_env%energy_functional == ec_functional_ext .AND. ec_env%basis_inconsistent) THEN + CALL cp_abort(__LOCATION__, & + "Exteranl Energy: Correction and ground state need to use the same basis. "// & + "Checked by comparing basis set names only.") + END IF ! ! set functional SELECT CASE (ec_env%energy_functional) @@ -288,6 +295,8 @@ SUBROUTINE init_ec_env(qs_env, ec_env, dft_section, ec_section) ec_env%ec_name = "Harris" CASE (ec_functional_dc) ec_env%ec_name = "DC-DFT" + CASE (ec_functional_ext) + ec_env%ec_name = "External Energy" CASE DEFAULT CPABORT("unknown energy correction") END SELECT @@ -502,6 +511,8 @@ SUBROUTINE ec_write_input(ec_env) WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "HARRIS FUNCTIONAL" CASE (ec_functional_dc) WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "DC-DFT" + CASE (ec_functional_ext) + WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Correction: ", "External" END SELECT WRITE (unit_nr, '()') diff --git a/src/ec_external.F b/src/ec_external.F new file mode 100644 index 0000000000..a89fbb8c3e --- /dev/null +++ b/src/ec_external.F @@ -0,0 +1,813 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2024 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Routines for an external energy correction on top of a Kohn-Sham calculation +!> \par History +!> 10.2024 created +!> \author JGH +! ************************************************************************************************** +MODULE ec_external + USE atomic_kind_types, ONLY: atomic_kind_type + USE cell_types, ONLY: cell_type + USE core_ppl, ONLY: build_core_ppl + USE core_ppnl, ONLY: build_core_ppnl + USE cp_control_types, ONLY: dft_control_type + USE cp_dbcsr_api, ONLY: dbcsr_add,& + dbcsr_copy,& + dbcsr_create,& + dbcsr_p_type,& + dbcsr_set,& + dbcsr_type,& + dbcsr_type_symmetric + USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl + USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,& + cp_dbcsr_sm_fm_multiply,& + dbcsr_allocate_matrix_set,& + dbcsr_deallocate_matrix_set + USE cp_fm_struct, ONLY: cp_fm_struct_create,& + cp_fm_struct_release,& + cp_fm_struct_type + USE cp_fm_types, ONLY: cp_fm_create,& + cp_fm_get_diag,& + cp_fm_get_info,& + cp_fm_maxabsval,& + cp_fm_release,& + cp_fm_set_all,& + cp_fm_to_fm,& + cp_fm_type + USE cp_log_handling, ONLY: cp_get_default_logger,& + cp_logger_get_default_unit_nr,& + cp_logger_type + USE ec_env_types, ONLY: energy_correction_type + USE kinds, ONLY: default_string_length,& + dp + USE mathlib, ONLY: det_3x3 + USE message_passing, ONLY: mp_para_env_type + USE parallel_gemm_api, ONLY: parallel_gemm + USE particle_types, ONLY: particle_type + USE physcon, ONLY: pascal + USE qs_core_energies, ONLY: calculate_ptrace + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_force_types, ONLY: qs_force_type + USE qs_kind_types, ONLY: qs_kind_type + USE qs_kinetic, ONLY: build_kinetic_matrix + USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_mo_types, ONLY: get_mo_set,& + mo_set_type + USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type + USE qs_overlap, ONLY: build_overlap_matrix + USE qs_rho_types, ONLY: qs_rho_get,& + qs_rho_type + USE virial_methods, ONLY: one_third_sum_diag + USE virial_types, ONLY: virial_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + +! *** Global parameters *** + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_external' + + PUBLIC :: ec_ext_energy + +CONTAINS + +! ************************************************************************************************** +!> \brief External energy method +!> \param qs_env ... +!> \param ec_env ... +!> \param calculate_forces ... +!> \par History +!> 10.2024 created +!> \author JGH +! ************************************************************************************************** + SUBROUTINE ec_ext_energy(qs_env, ec_env, calculate_forces) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(energy_correction_type), POINTER :: ec_env + LOGICAL, INTENT(IN) :: calculate_forces + + CHARACTER(len=*), PARAMETER :: routineN = 'ec_ext_energy' + + INTEGER :: handle, ispin, nocc, nspins, unit_nr + REAL(KIND=dp) :: focc + TYPE(cp_fm_struct_type), POINTER :: fm_struct + TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, mo_occ, mo_ref + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(cp_logger_type), POINTER :: logger + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, dft_control=dft_control) + nspins = dft_control%nspins + + logger => cp_get_default_logger() + IF (logger%para_env%is_source()) THEN + unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + ELSE + unit_nr = -1 + END IF + + CALL get_qs_env(qs_env, mos=mos) + ALLOCATE (cpmos(nspins), mo_ref(nspins), mo_occ(nspins)) + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) + NULLIFY (fm_struct) + CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, & + template_fmstruct=mo_coeff%matrix_struct) + CALL cp_fm_create(cpmos(ispin), fm_struct) + CALL cp_fm_set_all(cpmos(ispin), 0.0_dp) + CALL cp_fm_create(mo_ref(ispin), fm_struct) + CALL cp_fm_set_all(mo_ref(ispin), 0.0_dp) + CALL cp_fm_create(mo_occ(ispin), fm_struct) + CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc) + CALL cp_fm_struct_release(fm_struct) + END DO + + CALL cp_fm_release(ec_env%mo_occ) + CALL cp_fm_release(ec_env%cpmos) + IF (ec_env%debug_external) THEN + ! + ec_env%mo_occ => mo_ref + CALL ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr) + ! + IF (calculate_forces) THEN + focc = 2.0_dp + IF (nspins == 1) focc = 4.0_dp + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) + CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_h(1, 1)%matrix, ec_env%mo_occ(ispin), & + cpmos(ispin), nocc, & + alpha=focc, beta=0.0_dp) + END DO + END IF + ec_env%cpmos => cpmos + ELSE + ! get external information +!deb CALL ec_ext_interface(qs_env, ec_env%ex, mo_ref, cpmos, calculate_forces, unit_nr) + ec_env%mo_occ => mo_ref + ec_env%cpmos => cpmos + CPABORT("EC EXT NYA") + END IF + + IF (calculate_forces) THEN + ! check for orbital rotations + CALL get_qs_env(qs_env, matrix_s=matrix_s) + DO ispin = 1, nspins + CALL align_vectors(ec_env%cpmos(ispin), ec_env%mo_occ(ispin), mo_occ(ispin), & + matrix_s(1)%matrix, unit_nr) + END DO + ! set up matrices for response + CALL ec_ext_setup(qs_env, ec_env, .TRUE., unit_nr) + ! orthogonality force + CALL matrix_r_forces(qs_env, ec_env%cpmos, ec_env%mo_occ, & + ec_env%matrix_w(1, 1)%matrix, unit_nr) + ELSE + CALL ec_ext_setup(qs_env, ec_env, .FALSE., unit_nr) + END IF + + CALL cp_fm_release(mo_occ) + + CALL timestop(handle) + + END SUBROUTINE ec_ext_energy + +! ************************************************************************************************** + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param ec_env ... +!> \param calculate_forces ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE ec_ext_debug(qs_env, ec_env, calculate_forces, unit_nr) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(energy_correction_type), POINTER :: ec_env + LOGICAL, INTENT(IN) :: calculate_forces + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_debug' + + CHARACTER(LEN=default_string_length) :: headline + INTEGER :: handle, ispin, nocc, nspins + TYPE(cp_fm_type), POINTER :: mo_coeff + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, matrix_s + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mo_set_type), DIMENSION(:), POINTER :: mos + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb + TYPE(qs_rho_type), POINTER :: rho + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env=qs_env, & + dft_control=dft_control, & + sab_orb=sab_orb, & + rho=rho, & + matrix_s_kp=matrix_s, & + matrix_h_kp=matrix_h) + + nspins = dft_control%nspins + CALL get_qs_env(qs_env, mos=mos) + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) + CALL cp_fm_to_fm(mo_coeff, ec_env%mo_occ(ispin), nocc) + END DO + + ! Core Hamiltonian matrix + IF (ASSOCIATED(ec_env%matrix_h)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_h) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_h, 1, 1) + headline = "CORE HAMILTONIAN MATRIX" + ALLOCATE (ec_env%matrix_h(1, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_h(1, 1)%matrix, name=TRIM(headline), & + template=matrix_h(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_h(1, 1)%matrix, sab_orb) + CALL dbcsr_copy(ec_env%matrix_h(1, 1)%matrix, matrix_h(1, 1)%matrix) + + ! Get density matrix of reference calculation + CALL qs_rho_get(rho, rho_ao_kp=matrix_p) + ! Use Core energy as model energy + CALL calculate_ptrace(ec_env%matrix_h, matrix_p, ec_env%ex, nspins) + + IF (calculate_forces) THEN + ! force of model energy + CALL ec_debug_force(qs_env, matrix_p, unit_nr) + END IF + + CALL timestop(handle) + + END SUBROUTINE ec_ext_debug + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param matrix_p ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE ec_debug_force(qs_env, matrix_p, unit_nr) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_debug_force' + + INTEGER :: handle, iounit, nder, nimages + INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index + LOGICAL :: calculate_forces, debug_forces, & + debug_stress, use_virial + REAL(KIND=dp) :: eps_ppnl, fconv + REAL(KIND=dp), DIMENSION(3) :: fodeb + REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(cell_type), POINTER :: cell + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb, sac_ppl, sap_ppnl + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(virial_type), POINTER :: virial + + CALL timeset(routineN, handle) + + debug_forces = .TRUE. + debug_stress = .TRUE. + iounit = unit_nr + + calculate_forces = .TRUE. + + ! no k-points possible + NULLIFY (cell, dft_control, force, ks_env, para_env, virial) + CALL get_qs_env(qs_env=qs_env, & + cell=cell, & + dft_control=dft_control, & + force=force, & + ks_env=ks_env, & + para_env=para_env, & + virial=virial) + nimages = dft_control%nimages + IF (nimages /= 1) THEN + CPABORT("K-points not implemented") + END IF + + ! check for virial + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + + fconv = 1.0E-9_dp*pascal/cell%deth + IF (debug_stress .AND. use_virial) THEN + sttot = virial%pv_virial + END IF + + ! check for GAPW/GAPW_XC + IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN + CPABORT("GAPW not implemented") + END IF + + ! get neighbor lists + NULLIFY (sab_orb, sac_ppl, sap_ppnl) + CALL get_qs_env(qs_env=qs_env, & + sab_orb=sab_orb, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl) + + ! initialize src matrix + NULLIFY (scrm) + CALL dbcsr_allocate_matrix_set(scrm, 1, 1) + ALLOCATE (scrm(1, 1)%matrix) + CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_p(1, 1)%matrix) + CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb) + + nder = 1 + IF (SIZE(matrix_p, 1) == 2) THEN + CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & + alpha_scalar=1.0_dp, beta_scalar=1.0_dp) + END IF + + ! kinetic energy + IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic + CALL build_kinetic_matrix(ks_env, matrixkp_t=scrm, & + matrix_name="KINETIC ENERGY MATRIX", & + basis_type="ORB", & + sab_nl=sab_orb, calculate_forces=.TRUE., & + matrixkp_p=matrix_p) + IF (debug_forces) THEN + fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dT ", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_ekinetic - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Pout*dT', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + IF (SIZE(matrix_p, 1) == 2) THEN + CALL dbcsr_add(matrix_p(1, 1)%matrix, matrix_p(2, 1)%matrix, & + alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) + END IF + + ! compute the ppl contribution to the core hamiltonian + NULLIFY (atomic_kind_set, particle_set, qs_kind_set) + CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, & + atomic_kind_set=atomic_kind_set) + + IF (ASSOCIATED(sac_ppl)) THEN + IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl + CALL build_core_ppl(scrm, matrix_p, force, & + virial, calculate_forces, use_virial, nder, & + qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, & + nimages, cell_to_index, "ORB") + IF (calculate_forces .AND. debug_forces) THEN + fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPL ", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_ppl - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Pout*dH_PPL', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + END IF + + ! compute the ppnl contribution to the core hamiltonian *** + eps_ppnl = dft_control%qs_control%eps_ppnl + IF (ASSOCIATED(sap_ppnl)) THEN + IF (calculate_forces .AND. debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl + CALL build_core_ppnl(scrm, matrix_p, force, & + virial, calculate_forces, use_virial, nder, & + qs_kind_set, atomic_kind_set, particle_set, & + sab_orb, sap_ppnl, eps_ppnl, & + nimages, cell_to_index, "ORB") + IF (calculate_forces .AND. debug_forces) THEN + fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pout*dH_PPNL", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_ppnl - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Pout*dH_PPNL', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + END IF + + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_virial - sttot) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Stress Pout*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") ' ' + END IF + + ! delete scr matrix + CALL dbcsr_deallocate_matrix_set(scrm) + + CALL timestop(handle) + + END SUBROUTINE ec_debug_force + +! ************************************************************************************************** + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param ec_env ... +!> \param calc_forces ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE ec_ext_setup(qs_env, ec_env, calc_forces, unit_nr) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(energy_correction_type), POINTER :: ec_env + LOGICAL, INTENT(IN) :: calc_forces + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'ec_ext_setup' + + CHARACTER(LEN=default_string_length) :: headline + INTEGER :: handle, ispin, nao, nocc, nspins + REAL(KIND=dp) :: a_max, c_max + TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct + TYPE(cp_fm_type) :: emat, ksmo, smo + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_ks, matrix_p, matrix_s + TYPE(dft_control_type), POINTER :: dft_control + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb + TYPE(qs_rho_type), POINTER :: rho + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env=qs_env, & + dft_control=dft_control, & + sab_orb=sab_orb, & + rho=rho, & + matrix_s_kp=matrix_s, & + matrix_ks_kp=matrix_ks, & + matrix_h_kp=matrix_h) + + nspins = dft_control%nspins + + ! KS Hamiltonian matrix + IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_ks, nspins, 1) + headline = "HAMILTONIAN MATRIX" + DO ispin = 1, nspins + ALLOCATE (ec_env%matrix_ks(ispin, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_ks(ispin, 1)%matrix, name=TRIM(headline), & + template=matrix_ks(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_ks(ispin, 1)%matrix, sab_orb) + CALL dbcsr_copy(ec_env%matrix_ks(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix) + END DO + + ! Overlap matrix + IF (ASSOCIATED(ec_env%matrix_s)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_s) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_s, 1, 1) + headline = "OVERLAP MATRIX" + ALLOCATE (ec_env%matrix_s(1, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_s(1, 1)%matrix, name=TRIM(headline), & + template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_s(1, 1)%matrix, sab_orb) + CALL dbcsr_copy(ec_env%matrix_s(1, 1)%matrix, matrix_s(1, 1)%matrix) + + ! density matrix + ! Get density matrix of reference calculation + CALL qs_rho_get(rho, rho_ao_kp=matrix_p) + IF (ASSOCIATED(ec_env%matrix_p)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_p) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_p, nspins, 1) + headline = "DENSITY MATRIX" + DO ispin = 1, nspins + ALLOCATE (ec_env%matrix_p(ispin, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_p(ispin, 1)%matrix, name=TRIM(headline), & + template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_p(ispin, 1)%matrix, sab_orb) + CALL dbcsr_copy(ec_env%matrix_p(ispin, 1)%matrix, matrix_p(ispin, 1)%matrix) + END DO + + IF (calc_forces) THEN + ! energy weighted density matrix + ! for security, we recalculate W here (stored in qs_env) + IF (ASSOCIATED(ec_env%matrix_w)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_w) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_w, nspins, 1) + headline = "ENERGY WEIGHTED DENSITY MATRIX" + DO ispin = 1, nspins + ALLOCATE (ec_env%matrix_w(ispin, 1)%matrix) + CALL dbcsr_create(ec_env%matrix_w(ispin, 1)%matrix, name=TRIM(headline), & + template=matrix_p(ispin, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_w(ispin, 1)%matrix, sab_orb) + CALL dbcsr_set(ec_env%matrix_w(ispin, 1)%matrix, 0.0_dp) + END DO + + ! hz matrix + IF (ASSOCIATED(ec_env%matrix_hz)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_hz) + CALL dbcsr_allocate_matrix_set(ec_env%matrix_hz, nspins) + headline = "Hz MATRIX" + DO ispin = 1, nspins + ALLOCATE (ec_env%matrix_hz(ispin)%matrix) + CALL dbcsr_create(ec_env%matrix_hz(ispin)%matrix, name=TRIM(headline), & + template=matrix_s(1, 1)%matrix, matrix_type=dbcsr_type_symmetric) + CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_hz(ispin)%matrix, sab_orb) + CALL dbcsr_set(ec_env%matrix_hz(ispin)%matrix, 0.0_dp) + END DO + + ! Test for consistency of orbitals and KS matrix + DO ispin = 1, nspins + mat_struct => ec_env%mo_occ(ispin)%matrix_struct + CALL cp_fm_create(ksmo, mat_struct) + CALL cp_fm_get_info(ksmo, nrow_global=nao, ncol_global=nocc) + CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_ks(ispin, 1)%matrix, ec_env%mo_occ(ispin), & + ksmo, nocc, alpha=1.0_dp, beta=0.0_dp) + CALL cp_fm_create(smo, mat_struct) + CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_s(1, 1)%matrix, ec_env%mo_occ(ispin), & + smo, nocc, alpha=1.0_dp, beta=0.0_dp) + CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, & + template_fmstruct=mat_struct) + CALL cp_fm_create(emat, fm_struct) + CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, ec_env%mo_occ(ispin), ksmo, 0.0_dp, emat) + CALL parallel_gemm('N', 'N', nao, nocc, nocc, -1.0_dp, smo, emat, 1.0_dp, ksmo) + CALL cp_fm_maxabsval(ksmo, a_max) + CALL cp_fm_struct_release(fm_struct) + CALL cp_fm_release(smo) + CALL cp_fm_release(ksmo) + CALL cp_fm_release(emat) + CALL cp_fm_maxabsval(ec_env%mo_occ(ispin), c_max) + IF (unit_nr > 0) THEN + WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO coeficients", ispin, c_max + WRITE (unit_nr, "(T3,A,T50,I2,T61,F20.12)") "External:: Max value of MO gradients", ispin, a_max + END IF + END DO + END IF + + CALL timestop(handle) + + END SUBROUTINE ec_ext_setup + +! ************************************************************************************************** +!> \brief ... +!> \param cpmos ... +!> \param mo_ref ... +!> \param mo_occ ... +!> \param matrix_s ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE align_vectors(cpmos, mo_ref, mo_occ, matrix_s, unit_nr) + TYPE(cp_fm_type), INTENT(IN) :: cpmos, mo_ref, mo_occ + TYPE(dbcsr_type), POINTER :: matrix_s + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'align_vectors' + + INTEGER :: handle, i, nao, nocc, scg + REAL(KIND=dp) :: a_max + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag + TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct + TYPE(cp_fm_type) :: emat, smo + + CALL timeset(routineN, handle) + + mat_struct => cpmos%matrix_struct + CALL cp_fm_create(smo, mat_struct) + CALL cp_fm_get_info(smo, nrow_global=nao, ncol_global=nocc) + CALL cp_dbcsr_sm_fm_multiply(matrix_s, mo_occ, smo, nocc, alpha=1.0_dp, beta=0.0_dp) + CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, nrow_global=nocc, & + template_fmstruct=mat_struct) + CALL cp_fm_create(emat, fm_struct) + CALL parallel_gemm('T', 'N', nocc, nocc, nao, 1.0_dp, mo_ref, smo, 0.0_dp, emat) + CALL parallel_gemm('N', 'N', nao, nocc, nocc, 1.0_dp, cpmos, emat, 0.0_dp, smo) + CALL cp_fm_to_fm(smo, cpmos) + CALL cp_fm_to_fm(mo_occ, mo_ref) + ! + ALLOCATE (diag(nocc)) + CALL cp_fm_get_diag(emat, diag) + a_max = nocc - SUM(diag) + scg = 0 + DO i = 1, nocc + IF (ABS(diag(i) + 1.0_dp) < 0.001) scg = scg + 1 + END DO + IF (unit_nr > 0) THEN + WRITE (unit_nr, "(T3,A,T61,F20.8)") "External:: Orbital rotation index", a_max + WRITE (unit_nr, "(T3,A,T71,I10)") "External:: Number of orbital phase changes", scg + END IF + + DEALLOCATE (diag) + CALL cp_fm_struct_release(fm_struct) + CALL cp_fm_release(smo) + CALL cp_fm_release(emat) + + CALL timestop(handle) + + END SUBROUTINE align_vectors + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param matrix_w ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE matrix_w_forces(qs_env, matrix_w, unit_nr) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_w + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_w_forces' + + INTEGER :: handle, iounit, nder, nimages + LOGICAL :: debug_forces, debug_stress, use_virial + REAL(KIND=dp) :: fconv + REAL(KIND=dp), DIMENSION(3) :: fodeb + REAL(KIND=dp), DIMENSION(3, 3) :: stdeb, sttot + TYPE(cell_type), POINTER :: cell + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: scrm + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(virial_type), POINTER :: virial + + CALL timeset(routineN, handle) + + debug_forces = .TRUE. + debug_stress = .TRUE. + + iounit = unit_nr + + ! no k-points possible + CALL get_qs_env(qs_env=qs_env, & + cell=cell, & + dft_control=dft_control, & + force=force, & + ks_env=ks_env, & + sab_orb=sab_orb, & + para_env=para_env, & + virial=virial) + nimages = dft_control%nimages + IF (nimages /= 1) THEN + CPABORT("K-points not implemented") + END IF + + ! check for virial + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + + fconv = 1.0E-9_dp*pascal/cell%deth + IF (debug_stress .AND. use_virial) THEN + sttot = virial%pv_virial + END IF + + ! initialize src matrix + NULLIFY (scrm) + CALL dbcsr_allocate_matrix_set(scrm, 1, 1) + ALLOCATE (scrm(1, 1)%matrix) + CALL dbcsr_create(scrm(1, 1)%matrix, template=matrix_w(1, 1)%matrix) + CALL cp_dbcsr_alloc_block_from_nbl(scrm(1, 1)%matrix, sab_orb) + + nder = 1 + IF (SIZE(matrix_w, 1) == 2) THEN + CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, & + alpha_scalar=1.0_dp, beta_scalar=1.0_dp) + END IF + + ! Overlap and kinetic energy matrices + IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap + CALL build_overlap_matrix(ks_env, matrixkp_s=scrm, & + matrix_name="OVERLAP MATRIX", & + basis_type_a="ORB", & + basis_type_b="ORB", & + sab_nl=sab_orb, calculate_forces=.TRUE., & + matrixkp_p=matrix_w) + + IF (debug_forces) THEN + fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wout*dS ", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_overlap - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Wout*dS', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + IF (SIZE(matrix_w, 1) == 2) THEN + CALL dbcsr_add(matrix_w(1, 1)%matrix, matrix_w(2, 1)%matrix, & + alpha_scalar=1.0_dp, beta_scalar=-1.0_dp) + END IF + + ! delete scrm matrix + CALL dbcsr_deallocate_matrix_set(scrm) + + CALL timestop(handle) + + END SUBROUTINE matrix_w_forces + +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param cpmos ... +!> \param mo_occ ... +!> \param matrix_r ... +!> \param unit_nr ... +! ************************************************************************************************** + SUBROUTINE matrix_r_forces(qs_env, cpmos, mo_occ, matrix_r, unit_nr) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, mo_occ + TYPE(dbcsr_type), POINTER :: matrix_r + INTEGER, INTENT(IN) :: unit_nr + + CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_r_forces' + + INTEGER :: handle, iounit, ispin, nao, nocc, nspins + LOGICAL :: debug_forces, debug_stress, use_virial + REAL(KIND=dp) :: fconv, focc + REAL(KIND=dp), DIMENSION(3) :: fodeb + REAL(KIND=dp), DIMENSION(3, 3) :: stdeb + TYPE(cell_type), POINTER :: cell + TYPE(cp_fm_struct_type), POINTER :: fm_struct, mat_struct + TYPE(cp_fm_type) :: chcmat, rcvec + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: scrm + TYPE(dft_control_type), POINTER :: dft_control + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(neighbor_list_set_p_type), DIMENSION(:), & + POINTER :: sab_orb + TYPE(qs_force_type), DIMENSION(:), POINTER :: force + TYPE(qs_ks_env_type), POINTER :: ks_env + TYPE(virial_type), POINTER :: virial + + CALL timeset(routineN, handle) + + debug_forces = .TRUE. + debug_stress = .TRUE. + iounit = unit_nr + + nspins = SIZE(mo_occ) + focc = 1.0_dp + IF (nspins == 1) focc = 2.0_dp + focc = 0.25_dp*focc + + CALL dbcsr_set(matrix_r, 0.0_dp) + DO ispin = 1, nspins + CALL cp_fm_get_info(cpmos(ispin), matrix_struct=fm_struct, nrow_global=nao, ncol_global=nocc) + CALL cp_fm_create(rcvec, fm_struct) + CALL cp_fm_struct_create(mat_struct, nrow_global=nocc, ncol_global=nocc, template_fmstruct=fm_struct) + CALL cp_fm_create(chcmat, mat_struct) + CALL parallel_gemm("T", "N", nocc, nocc, nao, 1.0_dp, mo_occ(ispin), cpmos(ispin), 0.0_dp, chcmat) + CALL parallel_gemm("N", "N", nao, nocc, nocc, 1.0_dp, mo_occ(ispin), chcmat, 0.0_dp, rcvec) + CALL cp_dbcsr_plus_fm_fm_t(matrix_r, matrix_v=rcvec, matrix_g=mo_occ(ispin), ncol=nocc, alpha=focc) + CALL cp_fm_struct_release(mat_struct) + CALL cp_fm_release(rcvec) + CALL cp_fm_release(chcmat) + END DO + + CALL get_qs_env(qs_env=qs_env, & + cell=cell, & + dft_control=dft_control, & + force=force, & + ks_env=ks_env, & + sab_orb=sab_orb, & + para_env=para_env, & + virial=virial) + ! check for virial + use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) + fconv = 1.0E-9_dp*pascal/cell%deth + + IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1) + IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap + NULLIFY (scrm) + CALL build_overlap_matrix(ks_env, matrix_s=scrm, & + matrix_name="OVERLAP MATRIX", & + basis_type_a="ORB", basis_type_b="ORB", & + sab_nl=sab_orb, calculate_forces=.TRUE., & + matrix_p=matrix_r) + IF (debug_forces) THEN + fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3) + CALL para_env%sum(fodeb) + IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb + END IF + IF (debug_stress .AND. use_virial) THEN + stdeb = fconv*(virial%pv_overlap - stdeb) + CALL para_env%sum(stdeb) + IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") & + 'STRESS| Wz ', one_third_sum_diag(stdeb), det_3x3(stdeb) + END IF + CALL dbcsr_deallocate_matrix_set(scrm) + + CALL timestop(handle) + + END SUBROUTINE matrix_r_forces + +END MODULE ec_external diff --git a/src/energy_corrections.F b/src/energy_corrections.F index 7035763346..d5cdb5d659 100644 --- a/src/energy_corrections.F +++ b/src/energy_corrections.F @@ -12,6 +12,7 @@ !> 09.2019 Moved from KG to Kohn-Sham !> 08.2022 Add Density-Corrected DFT methods !> 04.2023 Add hybrid functionals for DC-DFT +!> 10.2024 Add external energy method !> \author JGH ! ************************************************************************************************** MODULE energy_corrections @@ -80,6 +81,7 @@ MODULE energy_corrections USE ec_efield_local, ONLY: ec_efield_integrals,& ec_efield_local_operator USE ec_env_types, ONLY: energy_correction_type + USE ec_external, ONLY: ec_ext_energy USE ec_methods, ONLY: create_kernel,& ec_mos_init USE external_potential_types, ONLY: get_potential,& @@ -88,9 +90,10 @@ MODULE energy_corrections USE hfx_exx, ONLY: add_exx_to_rhs,& calculate_exx USE input_constants, ONLY: & - ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_matrix_sign, ec_matrix_tc2, & - ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, ot_precond_full_single_inverse, & - ot_precond_solver_default, vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, xc_vdw_fun_pairpot + ec_diagonalization, ec_functional_dc, ec_functional_ext, ec_functional_harris, & + ec_matrix_sign, ec_matrix_tc2, ec_matrix_trs4, ec_ot_atomic, ec_ot_diag, ec_ot_gs, & + ot_precond_full_single_inverse, ot_precond_solver_default, vdw_pairpot_dftd3, & + vdw_pairpot_dftd3bj, xc_vdw_fun_pairpot USE input_section_types, ONLY: section_get_ival,& section_get_lval,& section_vals_get,& @@ -200,6 +203,7 @@ MODULE energy_corrections !> \brief Energy Correction to a Kohn-Sham simulation !> Available energy corrections: (1) Harris energy functional !> (2) Density-corrected DFT +!> (3) Energy from external source !> !> \param qs_env ... !> \param ec_init ... @@ -370,6 +374,10 @@ SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr) ! Calculate Hartree and XC related energies here CALL ec_build_ks_matrix(qs_env, ec_env) + CASE (ec_functional_ext) + + CALL ec_ext_energy(qs_env, ec_env, calculate_forces=.FALSE.) + CASE DEFAULT CPABORT("unknown energy correction") END SELECT @@ -408,6 +416,10 @@ SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr) ec_env%matrix_w) CALL ec_dc_build_ks_matrix_force(qs_env, ec_env) + CASE (ec_functional_ext) + + CALL ec_ext_energy(qs_env, ec_env, calculate_forces=.TRUE.) + CASE DEFAULT CPABORT("unknown energy correction") END SELECT @@ -445,11 +457,18 @@ SUBROUTINE energy_correction_low(qs_env, ec_env, calculate_forces, unit_nr) CALL ec_properties(qs_env, ec_env) ! Deallocate Harris density and response density on grid - DO ispin = 1, nspins - CALL auxbas_pw_pool%give_back_pw(ec_env%rhoout_r(ispin)) - CALL auxbas_pw_pool%give_back_pw(ec_env%rhoz_r(ispin)) - END DO - DEALLOCATE (ec_env%rhoout_r, ec_env%rhoz_r) + IF (ASSOCIATED(ec_env%rhoout_r)) THEN + DO ispin = 1, nspins + CALL auxbas_pw_pool%give_back_pw(ec_env%rhoout_r(ispin)) + END DO + DEALLOCATE (ec_env%rhoout_r) + END IF + IF (ASSOCIATED(ec_env%rhoz_r)) THEN + DO ispin = 1, nspins + CALL auxbas_pw_pool%give_back_pw(ec_env%rhoz_r(ispin)) + END DO + DEALLOCATE (ec_env%rhoz_r) + END IF ! Deallocate matrices IF (ASSOCIATED(ec_env%matrix_ks)) CALL dbcsr_deallocate_matrix_set(ec_env%matrix_ks) @@ -2734,7 +2753,7 @@ SUBROUTINE ec_energy(ec_env, unit_nr) CALL timeset(routineN, handle) - nspins = SIZE(ec_env%matrix_ks, 1) + nspins = SIZE(ec_env%matrix_p, 1) DO ispin = 1, nspins CALL dbcsr_dot(ec_env%matrix_p(ispin, 1)%matrix, ec_env%matrix_s(1, 1)%matrix, trace) IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T65,F16.10)') 'Tr[PS] ', trace @@ -2783,6 +2802,13 @@ SUBROUTINE ec_energy(ec_env, unit_nr) WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Energy Functional ", ec_env%etotal END IF + CASE (ec_functional_ext) + + ec_env%etotal = ec_env%ex + IF (unit_nr > 0) THEN + WRITE (unit_nr, '(T3,A,T56,F25.15)') "Etotal Energy Functional ", ec_env%etotal + END IF + CASE DEFAULT CPASSERT(.FALSE.) diff --git a/src/input_constants.F b/src/input_constants.F index f1a6bab8d2..ff73a20dab 100644 --- a/src/input_constants.F +++ b/src/input_constants.F @@ -1159,7 +1159,8 @@ MODULE input_constants ! non-scf energy corrections INTEGER, PARAMETER, PUBLIC :: ec_functional_harris = 2001, & - ec_functional_dc = 2002 + ec_functional_dc = 2002, & + ec_functional_ext = 2003 ! Energy correction solver INTEGER, PARAMETER, PUBLIC :: ec_diagonalization = 1001, & diff --git a/src/input_cp2k_ec.F b/src/input_cp2k_ec.F index 660c70b28f..771a32be96 100644 --- a/src/input_cp2k_ec.F +++ b/src/input_cp2k_ec.F @@ -20,13 +20,14 @@ MODULE input_cp2k_ec high_print_level USE input_constants, ONLY: & bqb_opt_exhaustive, bqb_opt_normal, bqb_opt_off, bqb_opt_patient, bqb_opt_quick, & - ec_diagonalization, ec_functional_dc, ec_functional_harris, ec_ls_solver, ec_matrix_sign, & - ec_matrix_tc2, ec_matrix_trs4, ec_mo_solver, ec_ot_atomic, ec_ot_diag, ec_ot_gs, & - kg_cholesky, ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, & - ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, & - ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_sign_ns, & - ls_scf_sign_proot, ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, & - ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, precond_mlp + ec_diagonalization, ec_functional_dc, ec_functional_ext, ec_functional_harris, & + ec_ls_solver, ec_matrix_sign, ec_matrix_tc2, ec_matrix_trs4, ec_mo_solver, ec_ot_atomic, & + ec_ot_diag, ec_ot_gs, kg_cholesky, ls_cluster_atomic, ls_cluster_molecular, & + ls_s_inversion_hotelling, ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, & + ls_s_preconditioner_molecular, ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, & + ls_scf_sign_ns, ls_scf_sign_proot, ot_precond_full_all, ot_precond_full_kinetic, & + ot_precond_full_single, ot_precond_full_single_inverse, ot_precond_none, & + ot_precond_s_inverse, precond_mlp USE input_cp2k_mm, ONLY: create_dipoles_section USE input_cp2k_voronoi, ONLY: create_print_voronoi_section USE input_cp2k_xc, ONLY: create_xc_section @@ -101,10 +102,11 @@ SUBROUTINE create_ec_section(section) description="Functional used in energy correction", & usage="ENERGY_FUNCTIONAL HARRIS", & default_i_val=ec_functional_harris, & - enum_c_vals=s2a("HARRIS", "DCDFT"), & + enum_c_vals=s2a("HARRIS", "DCDFT", "EXTERNAL"), & enum_desc=s2a("Harris functional", & - "Density-corrected DFT"), & - enum_i_vals=(/ec_functional_harris, ec_functional_dc/)) + "Density-corrected DFT", & + "External calculated energy"), & + enum_i_vals=(/ec_functional_harris, ec_functional_dc, ec_functional_ext/)) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) @@ -128,6 +130,11 @@ SUBROUTINE create_ec_section(section) usage="DEBUG_STRESS T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, __LOCATION__, name="DEBUG_EXTERNAL_METHOD", & + description="Uses an internal pseudo-energy to test EXTERNAL energy method.", & + usage="DEBUG_EXTERNAL_METHOD T", default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) CALL keyword_create(keyword, __LOCATION__, name="SKIP_EC", & description="Skip EC calculation if ground-state calculation has not converged.", & diff --git a/src/qs_basis_gradient.F b/src/qs_basis_gradient.F index 3408c81dc8..5740eeb54c 100644 --- a/src/qs_basis_gradient.F +++ b/src/qs_basis_gradient.F @@ -29,7 +29,6 @@ MODULE qs_basis_gradient USE particle_types, ONLY: particle_type USE qs_core_hamiltonian, ONLY: build_core_hamiltonian_matrix USE qs_density_matrices, ONLY: calculate_density_matrix - USE qs_energy_matrix_w, ONLY: qs_energies_compute_w USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_force_types, ONLY: allocate_qs_force,& @@ -44,6 +43,7 @@ MODULE qs_basis_gradient USE qs_ks_types, ONLY: get_ks_env,& qs_ks_env_type,& set_ks_env + USE qs_matrix_w, ONLY: compute_matrix_w USE qs_mixing_utils, ONLY: mixing_allocate USE qs_mo_methods, ONLY: make_basis_sm USE qs_mo_types, ONLY: get_mo_set,& @@ -160,7 +160,7 @@ SUBROUTINE qs_basis_center_gradient(qs_env) END IF END IF ! time to compute the w matrix - CALL qs_energies_compute_w(qs_env, .TRUE.) + CALL compute_matrix_w(qs_env, .TRUE.) ! core hamiltonian forces CALL build_core_hamiltonian_matrix(qs_env=qs_env, calculate_forces=.TRUE.) diff --git a/src/qs_energy.F b/src/qs_energy.F index 7fd1278fd1..7dae680295 100644 --- a/src/qs_energy.F +++ b/src/qs_energy.F @@ -33,7 +33,7 @@ MODULE qs_energy qs_environment_type USE qs_harris_utils, ONLY: harris_energy_correction USE qs_ks_methods, ONLY: qs_ks_update_qs_env - USE qs_matrix_w, ONLY: qs_energies_compute_matrix_w + USE qs_matrix_w, ONLY: compute_matrix_w USE qs_nonscf, ONLY: nonscf USE qs_scf, ONLY: scf USE scf_control_types, ONLY: scf_control_type @@ -129,7 +129,7 @@ SUBROUTINE qs_energies(qs_env, consistent_energies, calc_forces) IF (.NOT. ASSOCIATED(qs_env%mp2_env)) THEN ! if calculate forces, time to compute the w matrix - CALL qs_energies_compute_matrix_w(qs_env, my_calc_forces) + CALL compute_matrix_w(qs_env, my_calc_forces) END IF END IF diff --git a/src/qs_energy_matrix_w.F b/src/qs_energy_matrix_w.F deleted file mode 100644 index 59717e569e..0000000000 --- a/src/qs_energy_matrix_w.F +++ /dev/null @@ -1,172 +0,0 @@ -!--------------------------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright 2000-2024 CP2K developers group ! -! ! -! SPDX-License-Identifier: GPL-2.0-or-later ! -!--------------------------------------------------------------------------------------------------! - -! ***************************************************************************** -!> \brief Utility subroutine for qs energy calculation -!> \par History -!> none -!> \author MK (29.10.2002) -! ***************************************************************************** -MODULE qs_energy_matrix_w - USE cp_control_types, ONLY: dft_control_type - USE cp_dbcsr_api, ONLY: dbcsr_p_type,& - dbcsr_set - USE cp_fm_struct, ONLY: cp_fm_struct_create,& - cp_fm_struct_release,& - cp_fm_struct_type - USE cp_fm_types, ONLY: cp_fm_create,& - cp_fm_release,& - cp_fm_type - USE kinds, ONLY: dp - USE kpoint_methods, ONLY: kpoint_density_matrices,& - kpoint_density_transform - USE kpoint_types, ONLY: kpoint_type - USE qs_density_matrices, ONLY: calculate_w_matrix,& - calculate_w_matrix_ot - USE qs_environment_types, ONLY: get_qs_env,& - qs_environment_type - USE qs_mo_types, ONLY: get_mo_set,& - mo_set_type - USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type - USE qs_rho_types, ONLY: qs_rho_get,& - qs_rho_type - USE scf_control_types, ONLY: scf_control_type -#include "./base/base_uses.f90" - - IMPLICIT NONE - - PRIVATE - -! *** Global parameters *** - - CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_matrix_w' - - PUBLIC :: qs_energies_compute_w - -CONTAINS - -! ***************************************************************************** -!> \brief Refactoring of qs_energies_scf. Moves computation of matrix_w -!> into separate subroutine -!> \param qs_env ... -!> \param calc_forces ... -!> \par History -!> 05.2013 created [Florian Schiffmann] -! ************************************************************************************************** - - SUBROUTINE qs_energies_compute_w(qs_env, calc_forces) - TYPE(qs_environment_type), POINTER :: qs_env - LOGICAL, INTENT(IN) :: calc_forces - - CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_compute_w' - - INTEGER :: handle, is, ispin, nao, nspin - LOGICAL :: do_kpoints, has_unit_metric - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_w, & - mo_derivs, rho_ao - TYPE(dft_control_type), POINTER :: dft_control - TYPE(mo_set_type), DIMENSION(:), POINTER :: mos - TYPE(mo_set_type), POINTER :: mo_set - TYPE(qs_rho_type), POINTER :: rho - TYPE(scf_control_type), POINTER :: scf_control - - CALL timeset(routineN, handle) - - ! if calculate forces, time to compute the w matrix - CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric) - - IF (calc_forces .AND. .NOT. has_unit_metric) THEN - CALL get_qs_env(qs_env, do_kpoints=do_kpoints) - - IF (do_kpoints) THEN - BLOCK - TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct - TYPE(cp_fm_type), POINTER :: mo_coeff - TYPE(cp_fm_type), DIMENSION(2) :: fmwork - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrix_w_kp - TYPE(kpoint_type), POINTER :: kpoints - TYPE(neighbor_list_set_p_type), DIMENSION(:), & - POINTER :: sab_nl - - CALL get_qs_env(qs_env, & - matrix_w_kp=matrix_w_kp, & - matrix_s_kp=matrix_s_kp, & - sab_orb=sab_nl, & - mos=mos, & - kpoints=kpoints) - - CALL get_mo_set(mos(1), mo_coeff=mo_coeff, nao=nao) - CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, & - template_fmstruct=mo_coeff%matrix_struct) - - DO is = 1, SIZE(fmwork) - CALL cp_fm_create(fmwork(is), matrix_struct=ao_ao_fmstruct) - END DO - CALL cp_fm_struct_release(ao_ao_fmstruct) - - ! energy weighted density matrices in k-space - CALL kpoint_density_matrices(kpoints, energy_weighted=.TRUE.) - ! energy weighted density matrices in real space - CALL kpoint_density_transform(kpoints, matrix_w_kp, .TRUE., & - matrix_s_kp(1, 1)%matrix, sab_nl, fmwork) - - DO is = 1, SIZE(fmwork) - CALL cp_fm_release(fmwork(is)) - END DO - - END BLOCK - ELSE - - NULLIFY (dft_control, rho_ao) - CALL get_qs_env(qs_env, & - matrix_w=matrix_w, & - matrix_ks=matrix_ks, & - matrix_s=matrix_s, & - mo_derivs=mo_derivs, & - scf_control=scf_control, & - mos=mos, & - rho=rho, & - dft_control=dft_control) - - CALL qs_rho_get(rho, rho_ao=rho_ao) - - nspin = SIZE(mos) - DO ispin = 1, nspin - mo_set => mos(ispin) - IF (dft_control%roks) THEN - IF (scf_control%use_ot) THEN - IF (ispin > 1) THEN - ! not very elegant, indeed ... - CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp) - ELSE - CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, & - matrix_w(ispin)%matrix, matrix_s(1)%matrix) - END IF - ELSE - CALL calculate_w_matrix(mo_set=mo_set, & - matrix_ks=matrix_ks(ispin)%matrix, & - matrix_p=rho_ao(ispin)%matrix, & - matrix_w=matrix_w(ispin)%matrix) - END IF - ELSE - IF (scf_control%use_ot) THEN - CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, & - matrix_w(ispin)%matrix, matrix_s(1)%matrix) - ELSE - CALL calculate_w_matrix(mo_set, matrix_w(ispin)%matrix) - END IF - END IF - END DO - - END IF - - END IF - CALL timestop(handle) - - END SUBROUTINE qs_energies_compute_w - -END MODULE qs_energy_matrix_w diff --git a/src/qs_local_properties.F b/src/qs_local_properties.F index 40e8ac45e3..4441eec1e3 100644 --- a/src/qs_local_properties.F +++ b/src/qs_local_properties.F @@ -43,13 +43,13 @@ MODULE qs_local_properties pw_r3d_rs_type USE qs_collocate_density, ONLY: calculate_rho_elec USE qs_core_energies, ONLY: calculate_ptrace - USE qs_energy_matrix_w, ONLY: qs_energies_compute_w USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_ks_methods, ONLY: calc_rho_tot_gspace USE qs_ks_types, ONLY: qs_ks_env_type,& set_ks_env + USE qs_matrix_w, ONLY: compute_matrix_w USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type USE qs_vxc, ONLY: qs_xc_density @@ -145,7 +145,7 @@ SUBROUTINE qs_local_energy(qs_env, energy_density) CALL set_ks_env(ks_env, matrix_w_kp=matrix_w) END IF ! band structure energy density - CALL qs_energies_compute_w(qs_env, .TRUE.) + CALL compute_matrix_w(qs_env, .TRUE.) CALL get_qs_env(qs_env, ks_env=ks_env, matrix_w_kp=matrix_w) CALL auxbas_pw_pool%create_pw(edens_r) CALL auxbas_pw_pool%create_pw(edens_g) diff --git a/src/qs_matrix_w.F b/src/qs_matrix_w.F index 22cee54af5..f442a22409 100644 --- a/src/qs_matrix_w.F +++ b/src/qs_matrix_w.F @@ -45,7 +45,7 @@ MODULE qs_matrix_w CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_matrix_w' - PUBLIC :: qs_energies_compute_matrix_w + PUBLIC :: compute_matrix_w CONTAINS @@ -58,11 +58,11 @@ MODULE qs_matrix_w !> 05.2013 created [Florian Schiffmann] ! ************************************************************************************************** - SUBROUTINE qs_energies_compute_matrix_w(qs_env, calc_forces) + SUBROUTINE compute_matrix_w(qs_env, calc_forces) TYPE(qs_environment_type), POINTER :: qs_env LOGICAL, INTENT(IN) :: calc_forces - CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_compute_matrix_w' + CHARACTER(len=*), PARAMETER :: routineN = 'compute_matrix_w' INTEGER :: handle, is, ispin, nao, nspin LOGICAL :: do_kpoints, has_unit_metric @@ -168,6 +168,6 @@ SUBROUTINE qs_energies_compute_matrix_w(qs_env, calc_forces) CALL timestop(handle) - END SUBROUTINE qs_energies_compute_matrix_w + END SUBROUTINE compute_matrix_w END MODULE qs_matrix_w diff --git a/src/qs_tddfpt2_fprint.F b/src/qs_tddfpt2_fprint.F index 465a11a6f8..e2376950ce 100644 --- a/src/qs_tddfpt2_fprint.F +++ b/src/qs_tddfpt2_fprint.F @@ -54,7 +54,7 @@ MODULE qs_tddfpt2_fprint USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix USE qs_ks_types, ONLY: qs_ks_env_type,& set_ks_env - USE qs_matrix_w, ONLY: qs_energies_compute_matrix_w + USE qs_matrix_w, ONLY: compute_matrix_w USE qs_p_env_types, ONLY: p_env_release,& qs_p_env_type USE qs_scf, ONLY: scf @@ -323,7 +323,7 @@ SUBROUTINE gs_forces(qs_env, gs_force) CALL qs_env_rebuild_pw_env(qs_env) CALL qs_energies_init(qs_env, .TRUE.) CALL scf(qs_env) - CALL qs_energies_compute_matrix_w(qs_env, .TRUE.) + CALL compute_matrix_w(qs_env, .TRUE.) NULLIFY (para_env) CALL get_qs_env(qs_env, para_env=para_env) diff --git a/src/response_solver.F b/src/response_solver.F index e00384fd50..2a819c424d 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -351,31 +351,39 @@ SUBROUTINE response_calculation(qs_env, ec_env) ! CPKS vector cpmos - RHS of response equation as Ax + b = 0 (sign of b) ! Sign is changed in linres_solver! ! Projector Q applied in linres_solver! - CALL get_qs_env(qs_env, mos=mos) - ALLOCATE (cpmos(nspins), mo_occ(nspins)) - DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) - NULLIFY (fm_struct) - CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, & - template_fmstruct=mo_coeff%matrix_struct) - CALL cp_fm_create(cpmos(ispin), fm_struct) - CALL cp_fm_set_all(cpmos(ispin), 0.0_dp) - CALL cp_fm_create(mo_occ(ispin), fm_struct) - CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc) - CALL cp_fm_struct_release(fm_struct) - END DO + IF (ASSOCIATED(ec_env%cpmos)) THEN - focc = 2.0_dp - IF (nspins == 1) focc = 4.0_dp - DO ispin = 1, nspins - CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) - CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_hz(ispin)%matrix, mo_occ(ispin), & - cpmos(ispin), nocc, & - alpha=focc, beta=0.0_dp) - END DO - CALL cp_fm_release(mo_occ) + CALL response_equation_new(qs_env, p_env, ec_env%cpmos, unit_nr) + + ELSE + CALL get_qs_env(qs_env, mos=mos) + ALLOCATE (cpmos(nspins), mo_occ(nspins)) + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) + NULLIFY (fm_struct) + CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, & + template_fmstruct=mo_coeff%matrix_struct) + CALL cp_fm_create(cpmos(ispin), fm_struct) + CALL cp_fm_set_all(cpmos(ispin), 0.0_dp) + CALL cp_fm_create(mo_occ(ispin), fm_struct) + CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc) + CALL cp_fm_struct_release(fm_struct) + END DO - CALL response_equation_new(qs_env, p_env, cpmos, unit_nr) + focc = 2.0_dp + IF (nspins == 1) focc = 4.0_dp + DO ispin = 1, nspins + CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc) + CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_hz(ispin)%matrix, mo_occ(ispin), & + cpmos(ispin), nocc, & + alpha=focc, beta=0.0_dp) + END DO + CALL cp_fm_release(mo_occ) + + CALL response_equation_new(qs_env, p_env, cpmos, unit_nr) + + CALL cp_fm_release(cpmos) + END IF ! Get the response density matrix, ! and energy-weighted response density matrix @@ -383,7 +391,6 @@ SUBROUTINE response_calculation(qs_env, ec_env) CALL dbcsr_copy(ec_env%matrix_z(ispin)%matrix, p_env%p1(ispin)%matrix) CALL dbcsr_copy(ec_env%matrix_wz(ispin)%matrix, p_env%w1(ispin)%matrix) END DO - CALL cp_fm_release(cpmos) CASE (ec_ls_solver) diff --git a/src/rpa_im_time_force_methods.F b/src/rpa_im_time_force_methods.F index 0006294a0a..9a75b5f215 100644 --- a/src/rpa_im_time_force_methods.F +++ b/src/rpa_im_time_force_methods.F @@ -117,7 +117,6 @@ MODULE rpa_im_time_force_methods USE qs_collocate_density, ONLY: calculate_rho_elec,& collocate_function USE qs_density_matrices, ONLY: calculate_whz_matrix - USE qs_energy_matrix_w, ONLY: qs_energies_compute_w USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type,& set_qs_env @@ -133,6 +132,7 @@ MODULE rpa_im_time_force_methods USE qs_ks_reference, ONLY: ks_ref_potential USE qs_ks_types, ONLY: set_ks_env USE qs_linres_types, ONLY: linres_control_type + USE qs_matrix_w, ONLY: compute_matrix_w USE qs_mo_types, ONLY: get_mo_set,& mo_set_type USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,& @@ -3775,7 +3775,7 @@ SUBROUTINE update_im_time_forces(p_env, matrix_hz, matrix_p_F, matrix_p_F_admm, !Add to it the SCF W matrix, except if EXX (because taken care of by HF response) IF (.NOT. do_exx) THEN - CALL qs_energies_compute_w(qs_env, calc_forces=.TRUE.) + CALL compute_matrix_w(qs_env, calc_forces=.TRUE.) CALL get_qs_env(qs_env, matrix_w=matrix_w) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, 1.0_dp) IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, 1.0_dp) diff --git a/tests/QS/regtest-ec/EXT1.inp b/tests/QS/regtest-ec/EXT1.inp new file mode 100644 index 0000000000..43eb2da60e --- /dev/null +++ b/tests/QS/regtest-ec/EXT1.inp @@ -0,0 +1,63 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT H2 + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + CHECK_ATOM_FORCE 1 z + DEBUG_FORCES T + DEBUG_STRESS_TENSOR F + DX 0.001 + STOP_ON_MISMATCH F +&END DEBUG + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &ENERGY_CORRECTION + DEBUG_EXTERNAL_METHOD + DEBUG_FORCES + ENERGY_FUNCTIONAL EXTERNAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END RESPONSE_SOLVER + &END ENERGY_CORRECTION + &MGRID + CUTOFF 300 + &END MGRID + &PRINT + &DERIVATIVES + &END DERIVATIVES + &END PRINT + &QS + EPS_DEFAULT 1.E-10 + &END QS + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END SCF + &XC + &XC_FUNCTIONAL + &PADE + &END PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + H 0.000000 0.000000 0.360000 + H 0.000000 0.000000 -0.360000 + &END COORD + &KIND H + BASIS_SET ORB DZVP-GTH-BLYP + POTENTIAL GTH-PADE-q1 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-ec/EXT2.inp b/tests/QS/regtest-ec/EXT2.inp new file mode 100644 index 0000000000..c38d52e105 --- /dev/null +++ b/tests/QS/regtest-ec/EXT2.inp @@ -0,0 +1,60 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT N2 + RUN_TYPE ENERGY_FORCE +&END GLOBAL + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_SET + POTENTIAL_FILE_NAME GTH_POTENTIALS + &ENERGY_CORRECTION + DEBUG_EXTERNAL_METHOD + DEBUG_FORCES + ENERGY_FUNCTIONAL EXTERNAL + &RESPONSE_SOLVER + METHOD MO_SOLVER + PRECONDITIONER FULL_SINGLE_INVERSE + &END RESPONSE_SOLVER + &END ENERGY_CORRECTION + &MGRID + CUTOFF 300 + &END MGRID + &PRINT + &DERIVATIVES + &END DERIVATIVES + &END PRINT + &QS + EPS_DEFAULT 1.E-14 + &END QS + &SCF + EPS_SCF 1.0E-7 + SCF_GUESS ATOMIC + &END SCF + &XC + &XC_FUNCTIONAL + &PADE + &END PADE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.000000 0.000000 0.000000 + H 0.000000 -0.757136 0.504545 + H 0.000000 0.757136 0.504545 + &END COORD + &KIND H + BASIS_SET ORB DZVP-GTH-BLYP + POTENTIAL GTH-PADE-q1 + &END KIND + &KIND O + BASIS_SET ORB DZVP-GTH-BLYP + POTENTIAL GTH-PADE-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/QS/regtest-ec/TEST_FILES b/tests/QS/regtest-ec/TEST_FILES index 8b18c15f42..c3b9a8029e 100644 --- a/tests/QS/regtest-ec/TEST_FILES +++ b/tests/QS/regtest-ec/TEST_FILES @@ -61,4 +61,7 @@ N2_ec-hfx.inp 11 1e-10 N2_ec-hfx-admm.inp 11 1e-10 -19.8434174923 # harris.inp 0 +# +EXT1.inp 0 +EXT2.inp 119 1e-08 1.12631842 #EOF diff --git a/tests/TEST_TYPES b/tests/TEST_TYPES index eb7bc9f421..663109814a 100644 --- a/tests/TEST_TYPES +++ b/tests/TEST_TYPES @@ -1,4 +1,4 @@ -118 +119 Total energy:!3 MD| Potential energy!5 Total energy \[eV\]:!4 @@ -117,6 +117,7 @@ MOMENTS_TRACE_RE| 0.10000000E+000!3 MOMENTS_TRACE_IM| 0.10000000E+000!3 MOMENTS_TRACE_RE| 0.20000000E+000!3 POLARIZABILITY| 0.32436607E+002!4 +DEBUG:: Total Force !6 # # these are the tests the can be selected for regtesting. # do regtest will grep for test_grep (first column) and look if the numeric value