From d41ea5a552677d864facb5fae100b4982af8c9ae Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Fri, 15 Nov 2024 14:32:17 +0100 Subject: [PATCH] Harris and EHT methods: enable LS solver (#3780) --- src/CMakeLists.txt | 1 + src/dm_ls_scf.F | 117 +++++++++++++-------- src/dm_ls_scf_methods.F | 51 ++++++---- src/dm_ls_scf_qs.F | 129 +++++++++++++++++++++-- src/iterate_matrix.F | 35 ++++--- src/pao_methods.F | 2 +- src/qs_harris_utils.F | 1 + src/qs_nonscf.F | 133 ++---------------------- src/qs_nonscf_utils.F | 152 ++++++++++++++++++++++++++++ tests/QS/regtest-harris/TEST_FILES | 1 + tests/QS/regtest-harris/fl1.inp | 79 +++++++++++++++ tests/xTB/regtest-gfn0/TEST_FILES | 2 + tests/xTB/regtest-gfn0/ch2o_ls1.inp | 38 +++++++ tests/xTB/regtest-gfn0/ch2o_ls2.inp | 47 +++++++++ 14 files changed, 581 insertions(+), 207 deletions(-) create mode 100644 src/qs_nonscf_utils.F create mode 100644 tests/QS/regtest-harris/fl1.inp create mode 100644 tests/xTB/regtest-gfn0/ch2o_ls1.inp create mode 100644 tests/xTB/regtest-gfn0/ch2o_ls2.inp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 455c25b5e4..42af3e0f23 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -649,6 +649,7 @@ list( qs_neighbor_list_types.F qs_nl_hash_table_types.F qs_nonscf.F + qs_nonscf_utils.F qs_o3c_methods.F qs_o3c_types.F qs_oce_methods.F diff --git a/src/dm_ls_scf.F b/src/dm_ls_scf.F index 86c6f8a9c3..948b7f5195 100644 --- a/src/dm_ls_scf.F +++ b/src/dm_ls_scf.F @@ -37,13 +37,9 @@ MODULE dm_ls_scf density_matrix_tc2,& density_matrix_trs4,& ls_scf_init_matrix_S - USE dm_ls_scf_qs, ONLY: ls_scf_dm_to_ks,& - ls_scf_init_qs,& - ls_scf_qs_atomic_guess,& - matrix_ls_create,& - matrix_ls_to_qs,& - matrix_qs_to_ls,& - rho_mixing_ls_init + USE dm_ls_scf_qs, ONLY: & + ls_nonscf_energy, ls_nonscf_ks, ls_scf_dm_to_ks, ls_scf_init_qs, ls_scf_qs_atomic_guess, & + matrix_ls_create, matrix_ls_to_qs, matrix_qs_to_ls, rho_mixing_ls_init USE dm_ls_scf_types, ONLY: ls_scf_env_type USE ec_env_types, ONLY: energy_correction_type USE input_constants, ONLY: ls_cluster_atomic,& @@ -78,6 +74,7 @@ MODULE dm_ls_scf USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_ks_types, ONLY: qs_ks_env_type + USE qs_nonscf_utils, ONLY: qs_nonscf_print_summary USE qs_scf_post_gpw, ONLY: qs_scf_post_moments,& write_mo_free_results USE qs_scf_post_tb, ONLY: scf_post_calculation_tb @@ -100,33 +97,43 @@ MODULE dm_ls_scf !> \brief perform an linear scaling scf procedure: entry point !> !> \param qs_env ... +!> \param nonscf ... !> \par History !> 2010.10 created [Joost VandeVondele] !> \author Joost VandeVondele ! ************************************************************************************************** - SUBROUTINE ls_scf(qs_env) + SUBROUTINE ls_scf(qs_env, nonscf) TYPE(qs_environment_type), POINTER :: qs_env + LOGICAL, INTENT(IN), OPTIONAL :: nonscf CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf' INTEGER :: handle - LOGICAL :: pao_is_done + LOGICAL :: do_scf, pao_is_done TYPE(ls_scf_env_type), POINTER :: ls_scf_env CALL timeset(routineN, handle) + do_scf = .TRUE. + IF (PRESENT(nonscf)) do_scf = .NOT. nonscf + CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env) - CALL pao_optimization_start(qs_env, ls_scf_env) - - pao_is_done = .FALSE. - DO WHILE (.NOT. pao_is_done) - CALL ls_scf_init_scf(qs_env, ls_scf_env) - CALL pao_update(qs_env, ls_scf_env, pao_is_done) - CALL ls_scf_main(qs_env, ls_scf_env) - CALL pao_post_scf(qs_env, ls_scf_env, pao_is_done) - CALL ls_scf_post(qs_env, ls_scf_env) - END DO - CALL pao_optimization_end(ls_scf_env) + IF (do_scf) THEN + CALL pao_optimization_start(qs_env, ls_scf_env) + pao_is_done = .FALSE. + DO WHILE (.NOT. pao_is_done) + CALL ls_scf_init_scf(qs_env, ls_scf_env, .FALSE.) + CALL pao_update(qs_env, ls_scf_env, pao_is_done) + CALL ls_scf_main(qs_env, ls_scf_env, .FALSE.) + CALL pao_post_scf(qs_env, ls_scf_env, pao_is_done) + CALL ls_scf_post(qs_env, ls_scf_env) + END DO + CALL pao_optimization_end(ls_scf_env) + ELSE + CALL ls_scf_init_scf(qs_env, ls_scf_env, .TRUE.) + CALL ls_scf_main(qs_env, ls_scf_env, .TRUE.) + CALL ls_scf_post(qs_env, ls_scf_env) + END IF CALL timestop(handle) @@ -136,13 +143,15 @@ END SUBROUTINE ls_scf !> \brief initialization needed for scf !> \param qs_env ... !> \param ls_scf_env ... +!> \param nonscf ... !> \par History !> 2010.10 created [Joost VandeVondele] !> \author Joost VandeVondele ! ************************************************************************************************** - SUBROUTINE ls_scf_init_scf(qs_env, ls_scf_env) + SUBROUTINE ls_scf_init_scf(qs_env, ls_scf_env, nonscf) TYPE(qs_environment_type), POINTER :: qs_env TYPE(ls_scf_env_type) :: ls_scf_env + LOGICAL, INTENT(IN) :: nonscf CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_init_scf' @@ -210,7 +219,7 @@ SUBROUTINE ls_scf_init_scf(qs_env, ls_scf_env) CALL ls_scf_init_matrix_s(matrix_s(1)%matrix, ls_scf_env) ! get the initial guess for the SCF - CALL ls_scf_initial_guess(qs_env, ls_scf_env) + CALL ls_scf_initial_guess(qs_env, ls_scf_env, nonscf) IF (ls_scf_env%do_rho_mixing) THEN CALL rho_mixing_ls_init(qs_env, ls_scf_env) @@ -232,13 +241,15 @@ END SUBROUTINE ls_scf_init_scf !> \brief deal with the scf initial guess !> \param qs_env ... !> \param ls_scf_env ... +!> \param nonscf ... !> \par History !> 2012.11 created [Joost VandeVondele] !> \author Joost VandeVondele ! ************************************************************************************************** - SUBROUTINE ls_scf_initial_guess(qs_env, ls_scf_env) + SUBROUTINE ls_scf_initial_guess(qs_env, ls_scf_env, nonscf) TYPE(qs_environment_type), POINTER :: qs_env TYPE(ls_scf_env_type) :: ls_scf_env + LOGICAL, INTENT(IN) :: nonscf CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_initial_guess' INTEGER, PARAMETER :: aspc_guess = 2, atomic_guess = 1, & @@ -279,7 +290,8 @@ SUBROUTINE ls_scf_initial_guess(qs_env, ls_scf_env) ! how to get the initial guess SELECT CASE (initial_guess_type) CASE (atomic_guess) - CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env%energy_init) + CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env, ls_scf_env%energy_init, nonscf) + IF (unit_nr > 0) WRITE (unit_nr, '()') CASE (restart_guess) project_name = logger%iter_info%project_name DO ispin = 1, SIZE(ls_scf_env%matrix_p) @@ -293,7 +305,11 @@ SUBROUTINE ls_scf_initial_guess(qs_env, ls_scf_env) END DO ! directly go to computing the corresponding energy and ks matrix - CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0) + IF (nonscf) THEN + CALL ls_nonscf_ks(qs_env, ls_scf_env, ls_scf_env%energy_init) + ELSE + CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0) + END IF CASE (aspc_guess) CALL cite_reference(Kolafa2004) CALL cite_reference(Kuhne2007) @@ -346,7 +362,11 @@ SUBROUTINE ls_scf_initial_guess(qs_env, ls_scf_env) END IF ! compute corresponding energy and ks matrix - CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0) + IF (nonscf) THEN + CALL ls_nonscf_ks(qs_env, ls_scf_env, ls_scf_env%energy_init) + ELSE + CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0) + END IF END SELECT IF (unit_nr > 0) THEN @@ -440,13 +460,15 @@ END SUBROUTINE ls_scf_store_result !> \brief Main SCF routine. Can we keep it clean ? !> \param qs_env ... !> \param ls_scf_env ... +!> \param nonscf ... !> \par History !> 2010.10 created [Joost VandeVondele] !> \author Joost VandeVondele ! ************************************************************************************************** - SUBROUTINE ls_scf_main(qs_env, ls_scf_env) + SUBROUTINE ls_scf_main(qs_env, ls_scf_env, nonscf) TYPE(qs_environment_type), POINTER :: qs_env TYPE(ls_scf_env_type) :: ls_scf_env + LOGICAL, INTENT(IN), OPTIONAL :: nonscf CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_main' @@ -456,7 +478,7 @@ SUBROUTINE ls_scf_main(qs_env, ls_scf_env) LOGICAL :: check_convergence, diis_step, do_transport, extra_scf, maxscf_reached, & scf_converged, should_stop, transm_maxscf_reached, transm_scf_converged REAL(KIND=dp) :: energy_diff, energy_new, energy_old, & - eps_diis, t1, t2 + eps_diis, t1, t2, tdiag TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_ks_deviation, matrix_mixing_old @@ -560,7 +582,9 @@ SUBROUTINE ls_scf_main(qs_env, ls_scf_env) ELSE ! turn the KS matrix in a density matrix DO ispin = 1, nspin - IF (ls_scf_env%do_rho_mixing) THEN + IF (nonscf) THEN + CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin)) + ELSEIF (ls_scf_env%do_rho_mixing) THEN CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin)) ELSE IF (iscf == 1) THEN @@ -652,14 +676,16 @@ SUBROUTINE ls_scf_main(qs_env, ls_scf_env) CALL density_matrix_tc2(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, & nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), & ls_scf_env%lumo_spin(ispin), non_monotonic=ls_scf_env%non_monotonic, & - eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos) + eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, & + iounit=-1) CASE (ls_scf_trs4) CALL density_matrix_trs4(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, & nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), & ls_scf_env%lumo_spin(ispin), ls_scf_env%mu_spin(ispin), & dynamic_threshold=ls_scf_env%dynamic_threshold, & matrix_ks_deviation=matrix_ks_deviation(ispin), & - eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos) + eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, & + iounit=-1) CASE (ls_scf_pexsi) IF (ls_scf_env%has_s_preconditioner) & CPABORT("S preconditioning not implemented in combination with the PEXSI library. ") @@ -684,22 +710,31 @@ SUBROUTINE ls_scf_main(qs_env, ls_scf_env) END IF ! compute the corresponding new energy KS matrix and new energy - CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf) + IF (nonscf) THEN + CALL ls_nonscf_energy(qs_env, ls_scf_env) + ELSE + CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf) + END IF IF (ls_scf_env%do_pexsi) THEN CALL pexsi_to_qs(ls_scf_env, qs_env, kTS=ls_scf_env%pexsi%kTS) END IF - ! report current SCF loop - energy_diff = energy_new - energy_old - energy_old = energy_new - t2 = m_walltime() - IF (unit_nr > 0) THEN - WRITE (unit_nr, *) - WRITE (unit_nr, '(T2,A,I6,F20.9,F20.9,F12.6)') "SCF", iscf, energy_new, energy_diff, t2 - t1 - WRITE (unit_nr, *) - CALL m_flush(unit_nr) + IF (nonscf) THEN + tdiag = t2 - t1 + CALL qs_nonscf_print_summary(qs_env, tdiag, ls_scf_env%nelectron_total, unit_nr) + EXIT + ELSE + ! report current SCF loop + energy_diff = energy_new - energy_old + energy_old = energy_new + IF (unit_nr > 0) THEN + WRITE (unit_nr, *) + WRITE (unit_nr, '(T2,A,I6,F20.9,F20.9,F12.6)') "SCF", iscf, energy_new, energy_diff, t2 - t1 + WRITE (unit_nr, *) + CALL m_flush(unit_nr) + END IF END IF IF (do_transport) THEN diff --git a/src/dm_ls_scf_methods.F b/src/dm_ls_scf_methods.F index 9aaa5554bd..45c6a51fc0 100644 --- a/src/dm_ls_scf_methods.F +++ b/src/dm_ls_scf_methods.F @@ -137,7 +137,8 @@ SUBROUTINE ls_scf_init_matrix_S(matrix_s, ls_scf_env) CALL matrix_sqrt_Newton_Schulz(ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_s_sqrt_inv, & ls_scf_env%matrix_s, ls_scf_env%eps_filter, & ls_scf_env%s_sqrt_order, & - ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos) + ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos, & + iounit=-1) CASE DEFAULT CPABORT("Unknown sqrt method.") END SELECT @@ -211,7 +212,8 @@ END SUBROUTINE ls_scf_init_matrix_s !> \author Joost VandeVondele ! ************************************************************************************************** SUBROUTINE compute_matrix_preconditioner(matrix_s, preconditioner_type, ls_mstruct, & - matrix_bs_sqrt, matrix_bs_sqrt_inv, threshold, order, eps_lanczos, max_iter_lanczos) + matrix_bs_sqrt, matrix_bs_sqrt_inv, threshold, order, & + eps_lanczos, max_iter_lanczos) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s INTEGER :: preconditioner_type @@ -291,7 +293,8 @@ SUBROUTINE compute_matrix_preconditioner(matrix_s, preconditioner_type, ls_mstru ! XXXXXXXXXXX CALL matrix_sqrt_Newton_Schulz(matrix_bs_sqrt, matrix_bs_sqrt_inv, matrix_bs, & threshold=MIN(threshold, 1.0E-10_dp), order=order, & - eps_lanczos=eps_lanczos, max_iter_lanczos=max_iter_lanczos) + eps_lanczos=eps_lanczos, max_iter_lanczos=max_iter_lanczos, & + iounit=-1) END SELECT CALL dbcsr_release(matrix_bs) @@ -526,7 +529,7 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o SELECT CASE (sign_method) CASE (ls_scf_sign_ns) - CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order) + CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, iounit=-1) CASE (ls_scf_sign_proot) CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order) CASE (ls_scf_sign_submatrix) @@ -547,7 +550,7 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o ! compute sign(inv(S)*H-I*mu) SELECT CASE (sign_method) CASE (ls_scf_sign_ns) - CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order) + CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order, iounit=-1) CASE (ls_scf_sign_proot) CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order) CASE (ls_scf_sign_submatrix) @@ -574,7 +577,8 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o 0.0_dp, matrix_tmp, filter_eps=threshold) CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp) frob_matrix = dbcsr_frobenius_norm(matrix_tmp) - IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix + IF (unit_nr > 0 .AND. frob_matrix > 0.001_dp) & + WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix IF (sign_symmetric) THEN CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, & @@ -684,7 +688,8 @@ SUBROUTINE density_matrix_sign_internal_mu(matrix_p, trace, mu, sign_method, mat 0.0_dp, matrix_tmp, filter_eps=threshold) CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp) frob_matrix = dbcsr_frobenius_norm(matrix_tmp) - IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix + IF (unit_nr > 0 .AND. frob_matrix > 0.001_dp) & + WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, & 0.0_dp, matrix_tmp, filter_eps=threshold) @@ -712,6 +717,7 @@ END SUBROUTINE density_matrix_sign_internal_mu !> \param max_iter_lanczos ... !> \param eps_lanczos ... !> \param converged ... +!> \param iounit ... !> \par History !> 2012.06 created [Florian Thoele] !> \author Florian Thoele @@ -719,7 +725,7 @@ END SUBROUTINE density_matrix_sign_internal_mu SUBROUTINE density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, & nelectron, threshold, e_homo, e_lumo, e_mu, & dynamic_threshold, matrix_ks_deviation, & - max_iter_lanczos, eps_lanczos, converged) + max_iter_lanczos, eps_lanczos, converged, iounit) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv @@ -731,6 +737,7 @@ SUBROUTINE density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, & INTEGER, INTENT(IN) :: max_iter_lanczos REAL(KIND=dp), INTENT(IN) :: eps_lanczos LOGICAL, INTENT(OUT), OPTIONAL :: converged + INTEGER, INTENT(IN), OPTIONAL :: iounit CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_trs4' INTEGER, PARAMETER :: max_iter = 100 @@ -756,11 +763,15 @@ SUBROUTINE density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, & CALL timeset(routineN, handle) - logger => cp_get_default_logger() - IF (logger%para_env%is_source()) THEN - unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + IF (PRESENT(iounit)) THEN + unit_nr = iounit ELSE - unit_nr = -1 + 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 END IF do_dyn_threshold = .FALSE. @@ -1010,11 +1021,12 @@ END SUBROUTINE density_matrix_trs4 !> \param non_monotonic ... !> \param eps_lanczos ... !> \param max_iter_lanczos ... +!> \param iounit ... !> \author Jonathan Mullin ! ************************************************************************************************** SUBROUTINE density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, & nelectron, threshold, e_homo, e_lumo, & - non_monotonic, eps_lanczos, max_iter_lanczos) + non_monotonic, eps_lanczos, max_iter_lanczos, iounit) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv @@ -1024,6 +1036,7 @@ SUBROUTINE density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, & LOGICAL, INTENT(IN), OPTIONAL :: non_monotonic REAL(KIND=dp), INTENT(IN) :: eps_lanczos INTEGER, INTENT(IN) :: max_iter_lanczos + INTEGER, INTENT(IN), OPTIONAL :: iounit CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_tc2' INTEGER, PARAMETER :: max_iter = 100 @@ -1040,11 +1053,15 @@ SUBROUTINE density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, & CALL timeset(routineN, handle) - logger => cp_get_default_logger() - IF (logger%para_env%is_source()) THEN - unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + IF (PRESENT(iounit)) THEN + unit_nr = iounit ELSE - unit_nr = -1 + 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 END IF do_non_monotonic = .FALSE. diff --git a/src/dm_ls_scf_qs.F b/src/dm_ls_scf_qs.F index 2f369046b8..7caee96f92 100644 --- a/src/dm_ls_scf_qs.F +++ b/src/dm_ls_scf_qs.F @@ -47,12 +47,15 @@ MODULE dm_ls_scf_qs pw_r3d_rs_type USE qs_atomic_block, ONLY: calculate_atomic_block_dm USE qs_collocate_density, ONLY: calculate_rho_elec + USE qs_core_energies, ONLY: calculate_ptrace USE qs_density_mixing_types, ONLY: direct_mixing_nr,& gspace_mixing_nr USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type USE qs_gspace_mixing, ONLY: gspace_mixing + USE qs_harris_types, ONLY: harris_type + USE qs_harris_utils, ONLY: harris_density_update USE qs_initial_guess, ONLY: calculate_mopac_dm USE qs_kind_types, ONLY: qs_kind_type USE qs_ks_methods, ONLY: qs_ks_update_qs_env @@ -78,8 +81,8 @@ MODULE dm_ls_scf_qs CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_qs' PUBLIC :: matrix_ls_create, matrix_qs_to_ls, matrix_ls_to_qs, ls_scf_init_qs, & - ls_scf_dm_to_ks, ls_scf_qs_atomic_guess, write_matrix_to_cube, rho_mixing_ls_init, & - matrix_decluster + ls_nonscf_ks, ls_nonscf_energy, ls_scf_dm_to_ks, ls_scf_qs_atomic_guess, & + write_matrix_to_cube, rho_mixing_ls_init, matrix_decluster CONTAINS @@ -448,20 +451,24 @@ END SUBROUTINE ls_scf_init_qs ! ************************************************************************************************** !> \brief get an atomic initial guess !> \param qs_env ... +!> \param ls_scf_env ... !> \param energy ... +!> \param nonscf ... !> \par History !> 2012.11 created [Joost VandeVondele] !> \author Joost VandeVondele ! ************************************************************************************************** - SUBROUTINE ls_scf_qs_atomic_guess(qs_env, energy) + SUBROUTINE ls_scf_qs_atomic_guess(qs_env, ls_scf_env, energy, nonscf) TYPE(qs_environment_type), POINTER :: qs_env + TYPE(ls_scf_env_type) :: ls_scf_env REAL(KIND=dp) :: energy + LOGICAL, INTENT(IN), OPTIONAL :: nonscf CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_qs_atomic_guess' INTEGER :: handle, nspin, unit_nr INTEGER, DIMENSION(2) :: nelectron_spin - LOGICAL :: has_unit_metric + LOGICAL :: do_scf, has_unit_metric TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set TYPE(cp_logger_type), POINTER :: logger TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao @@ -513,11 +520,16 @@ SUBROUTINE ls_scf_qs_atomic_guess(qs_env, energy) nspin, nelectron_spin, unit_nr, para_env) END IF - CALL qs_rho_update_rho(rho, qs_env=qs_env) - CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) - CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.) - - energy = qs_energy%total + do_scf = .TRUE. + IF (PRESENT(nonscf)) do_scf = .NOT. nonscf + IF (do_scf) THEN + CALL qs_rho_update_rho(rho, qs_env=qs_env) + CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) + CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.) + energy = qs_energy%total + ELSE + CALL ls_nonscf_ks(qs_env, ls_scf_env, energy) + END IF CALL timestop(handle) @@ -604,6 +616,105 @@ SUBROUTINE ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf) END SUBROUTINE ls_scf_dm_to_ks +! ************************************************************************************************** +!> \brief use the external density in ls_scf_env to compute the new KS matrix +!> \param qs_env ... +!> \param ls_scf_env ... +!> \param energy_new ... +! ************************************************************************************************** + SUBROUTINE ls_nonscf_ks(qs_env, ls_scf_env, energy_new) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(ls_scf_env_type) :: ls_scf_env + REAL(KIND=dp) :: energy_new + + CHARACTER(len=*), PARAMETER :: routineN = 'ls_nonscf_ks' + + INTEGER :: handle, ispin, nspin + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao + TYPE(harris_type), POINTER :: harris_env + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(qs_energy_type), POINTER :: energy + TYPE(qs_rho_type), POINTER :: rho + + NULLIFY (energy, rho, rho_ao) + CALL timeset(routineN, handle) + + nspin = ls_scf_env%nspins + CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho) + CALL qs_rho_get(rho, rho_ao=rho_ao) + + ! set the new density matrix + DO ispin = 1, nspin + CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), & + ls_scf_env%ls_mstruct, covariant=.FALSE.) + END DO + + IF (qs_env%harris_method) THEN + CALL get_qs_env(qs_env, harris_env=harris_env) + CALL harris_density_update(qs_env, harris_env) + END IF + ! compute the corresponding KS matrix and new energy + CALL qs_rho_update_rho(rho, qs_env=qs_env) + IF (ls_scf_env%do_rho_mixing) THEN + CPABORT("P mixing not implemented in linear scaling NONSCF. ") + END IF + + CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) + CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., & + just_energy=.FALSE., print_active=.TRUE.) + energy_new = energy%total + + CALL timestop(handle) + + END SUBROUTINE ls_nonscf_ks + +! ************************************************************************************************** +!> \brief use the new density matrix in ls_scf_env to compute the new energy +!> \param qs_env ... +!> \param ls_scf_env ... +! ************************************************************************************************** + SUBROUTINE ls_nonscf_energy(qs_env, ls_scf_env) + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(ls_scf_env_type) :: ls_scf_env + + CHARACTER(len=*), PARAMETER :: routineN = 'ls_nonscf_energy' + + INTEGER :: handle, ispin, nspin + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_ks, rho_ao + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(qs_energy_type), POINTER :: energy + TYPE(qs_rho_type), POINTER :: rho + + NULLIFY (energy, rho, rho_ao) + CALL timeset(routineN, handle) + IF (qs_env%qmmm) THEN + CPABORT("NYA") + END IF + + nspin = ls_scf_env%nspins + CALL get_qs_env(qs_env, para_env=para_env, energy=energy, rho=rho) + CALL qs_rho_get(rho, rho_ao=rho_ao) + + ! set the new density matrix + DO ispin = 1, nspin + CALL matrix_ls_to_qs(rho_ao(ispin)%matrix, ls_scf_env%matrix_p(ispin), & + ls_scf_env%ls_mstruct, covariant=.FALSE.) + END DO + + CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) + + ! band energy : Tr(PH) + CALL get_qs_env(qs_env, matrix_ks=matrix_ks) + CALL calculate_ptrace(matrix_ks, rho_ao, energy%band, nspin) + ! core energy : Tr(Ph) + energy%total = energy%total - energy%core + CALL get_qs_env(qs_env, matrix_h=matrix_h) + CALL calculate_ptrace(matrix_h, rho_ao, energy%core, nspin) + + CALL timestop(handle) + + END SUBROUTINE ls_nonscf_energy + ! ************************************************************************************************** !> \brief ... !> \param qs_env ... diff --git a/src/iterate_matrix.F b/src/iterate_matrix.F index 8df188e448..3f58b5232a 100644 --- a/src/iterate_matrix.F +++ b/src/iterate_matrix.F @@ -672,16 +672,17 @@ END SUBROUTINE invert_Hotelling !> \param matrix ... !> \param threshold ... !> \param sign_order ... +!> \param iounit ... !> \par History !> 2010.10 created [Joost VandeVondele] !> 2019.05 extended to order byxond 2 [Robert Schade] !> \author Joost VandeVondele, Robert Schade ! ************************************************************************************************** - SUBROUTINE matrix_sign_Newton_Schulz(matrix_sign, matrix, threshold, sign_order) + SUBROUTINE matrix_sign_Newton_Schulz(matrix_sign, matrix, threshold, sign_order, iounit) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sign, matrix REAL(KIND=dp), INTENT(IN) :: threshold - INTEGER, INTENT(IN), OPTIONAL :: sign_order + INTEGER, INTENT(IN), OPTIONAL :: sign_order, iounit CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_sign_Newton_Schulz' @@ -696,11 +697,15 @@ SUBROUTINE matrix_sign_Newton_Schulz(matrix_sign, matrix, threshold, sign_order) CALL timeset(routineN, handle) - logger => cp_get_default_logger() - IF (logger%para_env%is_source()) THEN - unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + IF (PRESENT(iounit)) THEN + unit_nr = iounit ELSE - unit_nr = -1 + 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 END IF IF (PRESENT(sign_order)) THEN @@ -1059,8 +1064,6 @@ SUBROUTINE matrix_sign_proot(matrix_sign, matrix, threshold, sign_order) symmetrize = .FALSE. CALL matrix_sqrt_proot(matrix_sqrt, matrix_sqrt_inv, matrix2, threshold, order, & 0.01_dp, 100, symmetrize, converged) -! call matrix_sqrt_Newton_Schulz(matrix_sqrt, matrix_sqrt_inv, matrix2, threshold, order, & -! 0.01_dp,100, symmetrize,converged) CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_sqrt_inv, 0.0_dp, matrix_sign, & filter_eps=threshold, flop=flop1) @@ -1615,18 +1618,20 @@ END SUBROUTINE matrix_sign_submatrix_mu_adjust !> \param max_iter_lanczos ... !> \param symmetrize ... !> \param converged ... +!> \param iounit ... !> \par History !> 2010.10 created [Joost VandeVondele] !> \author Joost VandeVondele ! ************************************************************************************************** SUBROUTINE matrix_sqrt_Newton_Schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, & - eps_lanczos, max_iter_lanczos, symmetrize, converged) + eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit) TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sqrt, matrix_sqrt_inv, matrix REAL(KIND=dp), INTENT(IN) :: threshold INTEGER, INTENT(IN) :: order REAL(KIND=dp), INTENT(IN) :: eps_lanczos INTEGER, INTENT(IN) :: max_iter_lanczos LOGICAL, OPTIONAL :: symmetrize, converged + INTEGER, INTENT(IN), OPTIONAL :: iounit CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_sqrt_Newton_Schulz' @@ -1642,11 +1647,15 @@ SUBROUTINE matrix_sqrt_Newton_Schulz(matrix_sqrt, matrix_sqrt_inv, matrix, thres CALL timeset(routineN, handle) - logger => cp_get_default_logger() - IF (logger%para_env%is_source()) THEN - unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.) + IF (PRESENT(iounit)) THEN + unit_nr = iounit ELSE - unit_nr = -1 + 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 END IF IF (PRESENT(converged)) converged = .FALSE. diff --git a/src/pao_methods.F b/src/pao_methods.F index c694b1df27..83cddc8fb2 100644 --- a/src/pao_methods.F +++ b/src/pao_methods.F @@ -878,7 +878,7 @@ SUBROUTINE pao_guess_initial_P(pao, qs_env, ls_scf_env) pao%need_initial_scf = .FALSE. pao%preopt_dm_file = "" ! load only for first MD step ELSE - CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env%energy_init) + CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env, ls_scf_env%energy_init) IF (pao%iw > 0) WRITE (pao%iw, '(A,F20.9)') & " PAO| Energy from initial atomic guess:", ls_scf_env%energy_init pao%need_initial_scf = .TRUE. diff --git a/src/qs_harris_utils.F b/src/qs_harris_utils.F index 88184bdc62..2fbf51963a 100644 --- a/src/qs_harris_utils.F +++ b/src/qs_harris_utils.F @@ -282,6 +282,7 @@ SUBROUTINE harris_density_update(qs_env, harris_env) END IF DEALLOCATE (density, coef) END DO + harris_env%rhoin%frozen = .TRUE. END IF CASE DEFAULT CPABORT("Illeagal value of harris_env%density_source") diff --git a/src/qs_nonscf.F b/src/qs_nonscf.F index 71814bf6a2..576709b634 100644 --- a/src/qs_nonscf.F +++ b/src/qs_nonscf.F @@ -19,14 +19,13 @@ MODULE qs_nonscf USE cp_log_handling, ONLY: cp_get_default_logger,& cp_logger_get_default_io_unit,& cp_logger_type + USE dm_ls_scf, ONLY: ls_scf USE input_section_types, ONLY: section_vals_get_subs_vals,& section_vals_type USE kinds, ONLY: dp USE kpoint_types, ONLY: kpoint_type - USE machine, ONLY: m_flush,& - m_walltime + USE machine, ONLY: m_walltime USE message_passing, ONLY: mp_para_env_type - USE qs_charges_types, ONLY: qs_charges_type USE qs_core_energies, ONLY: calculate_ptrace USE qs_energy_types, ONLY: qs_energy_type USE qs_environment_types, ONLY: get_qs_env,& @@ -36,6 +35,7 @@ MODULE qs_nonscf USE qs_ks_types, ONLY: qs_ks_did_change,& qs_ks_env_type USE qs_mo_types, ONLY: mo_set_type + USE qs_nonscf_utils, ONLY: qs_nonscf_print_summary USE qs_rho_types, ONLY: qs_rho_get,& qs_rho_type USE qs_scf, ONLY: init_scf_loop @@ -75,18 +75,17 @@ SUBROUTINE nonscf(qs_env) TYPE(qs_scf_env_type), POINTER :: scf_env TYPE(scf_control_type), POINTER :: scf_control - CALL get_qs_env(qs_env, scf_env=scf_env, scf_control=scf_control, & - dft_control=dft_control) + CALL get_qs_env(qs_env, dft_control=dft_control) IF (dft_control%qs_control%do_ls_scf) THEN ! Density matrix based solver - CPABORT("NOT AVAILABLE") - !CALL ls_scf(qs_env=qs_env) + CALL ls_scf(qs_env, nonscf=.TRUE.) ELSE ! Wavefunction based solver + CALL get_qs_env(qs_env, scf_env=scf_env, scf_control=scf_control) IF (.NOT. ASSOCIATED(scf_env)) THEN CALL qs_scf_env_initialize(qs_env, scf_env) CALL set_qs_env(qs_env, scf_env=scf_env) @@ -219,128 +218,10 @@ SUBROUTINE do_nonscf(qs_env, scf_env, scf_control) t2 = m_walltime() tdiag = t2 - t1 - CALL qs_nonscf_print_summary(qs_env, tdiag, iounit) + CALL qs_nonscf_print_summary(qs_env, tdiag, scf_env%nelectron, iounit) CALL timestop(handle) END SUBROUTINE do_nonscf -! ************************************************************************************************** -!> \brief writes a summary of information after diagonalization -!> \param qs_env ... -!> \param tdiag ... -!> \param iounit ... -! ************************************************************************************************** - SUBROUTINE qs_nonscf_print_summary(qs_env, tdiag, iounit) - TYPE(qs_environment_type), POINTER :: qs_env - REAL(KIND=dp), INTENT(IN) :: tdiag - INTEGER, INTENT(IN) :: iounit - - INTEGER :: nelectron_total - REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r - TYPE(dft_control_type), POINTER :: dft_control - TYPE(qs_charges_type), POINTER :: qs_charges - TYPE(qs_energy_type), POINTER :: energy - TYPE(qs_rho_type), POINTER :: rho - TYPE(qs_scf_env_type), POINTER :: scf_env - - IF (iounit > 0) THEN - CALL get_qs_env(qs_env=qs_env, energy=energy, dft_control=dft_control) - IF (qs_env%harris_method) THEN - CPASSERT(.NOT. dft_control%qs_control%gapw) - CPASSERT(.NOT. dft_control%qs_control%gapw_xc) - - CALL get_qs_env(qs_env=qs_env, rho=rho, scf_env=scf_env, qs_charges=qs_charges) - nelectron_total = scf_env%nelectron - CALL qs_rho_get(rho, tot_rho_r=tot_rho_r) - WRITE (UNIT=iounit, FMT="(/,(T3,A,T41,2F20.10))") & - "Electronic density on regular grids: ", & - SUM(tot_rho_r), & - SUM(tot_rho_r) + nelectron_total, & - "Core density on regular grids:", & - qs_charges%total_rho_core_rspace, & - qs_charges%total_rho_core_rspace - REAL(nelectron_total + dft_control%charge, dp) - WRITE (UNIT=iounit, FMT="(T3,A,T41,F20.10)") & - "Total charge density on r-space grids: ", & - SUM(tot_rho_r) + & - qs_charges%total_rho_core_rspace, & - "Total charge density g-space grids: ", & - qs_charges%total_rho_gspace - - WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & - "Diagonalization", "Time:", tdiag, energy%band - - WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & - "Core Hamiltonian energy: ", energy%core, & - "Overlap energy of the core charge distribution:", energy%core_overlap, & - "Self energy of the core charge distribution: ", energy%core_self, & - "Hartree energy: ", energy%hartree, & - "Exchange-correlation energy: ", energy%exc - IF (energy%dispersion /= 0.0_dp) & - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "Dispersion energy: ", energy%dispersion - IF (energy%gcp /= 0.0_dp) & - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "gCP energy: ", energy%gcp - IF (energy%efield /= 0.0_dp) & - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "Electric field interaction energy: ", energy%efield - - ELSEIF (dft_control%qs_control%semi_empirical) THEN - CPABORT("NONSCF not available") - ELSEIF (dft_control%qs_control%dftb) THEN - CPASSERT(energy%dftb3 == 0.0_dp) - energy%total = energy%total + energy%band + energy%qmmm_el - WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & - "Diagonalization", "Time:", tdiag, energy%total - WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & - "Core Hamiltonian energy: ", energy%core, & - "Repulsive potential energy: ", energy%repulsive, & - "Dispersion energy: ", energy%dispersion - IF (energy%efield /= 0.0_dp) & - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "Electric field interaction energy: ", energy%efield - ELSEIF (dft_control%qs_control%xtb) THEN - energy%total = energy%total + energy%band + energy%qmmm_el - WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & - "Diagonalization", "Time:", tdiag, energy%total - CPASSERT(dft_control%qs_control%xtb_control%gfn_type == 0) - WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & - "Core Hamiltonian energy: ", energy%core, & - "Repulsive potential energy: ", energy%repulsive, & - "SRB Correction energy: ", energy%srb, & - "Charge equilibration energy: ", energy%eeq, & - "Dispersion energy: ", energy%dispersion - IF (dft_control%qs_control%xtb_control%do_nonbonded) & - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "Correction for nonbonded interactions: ", energy%xtb_nonbonded - IF (energy%efield /= 0.0_dp) & - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "Electric field interaction energy: ", energy%efield - ELSE - CPABORT("NONSCF not available") - END IF - IF (dft_control%smear) THEN - WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & - "Electronic entropic energy:", energy%kTS - WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & - "Fermi energy:", energy%efermi - END IF - IF (energy%qmmm_el /= 0.0_dp) THEN - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "QM/MM Electrostatic energy: ", energy%qmmm_el - IF (qs_env%qmmm_env_qm%image_charge) THEN - WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & - "QM/MM image charge energy: ", energy%image_charge - END IF - END IF - - WRITE (UNIT=iounit, FMT="(/,(T3,A,T56,F25.14))") & - "Total energy: ", energy%total - - CALL m_flush(iounit) - END IF - - END SUBROUTINE qs_nonscf_print_summary - END MODULE qs_nonscf diff --git a/src/qs_nonscf_utils.F b/src/qs_nonscf_utils.F new file mode 100644 index 0000000000..b8b0f3bed8 --- /dev/null +++ b/src/qs_nonscf_utils.F @@ -0,0 +1,152 @@ +!--------------------------------------------------------------------------------------------------! +! 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 Quickstep NON-SCF run. +!> \par History +!> - initial setup [JGH, 2024] +!> \author JGH (13.05.2024) +! ************************************************************************************************** +MODULE qs_nonscf_utils + USE cp_control_types, ONLY: dft_control_type + USE kinds, ONLY: dp + USE machine, ONLY: m_flush + USE qs_charges_types, ONLY: qs_charges_type + USE qs_energy_types, ONLY: qs_energy_type + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE qs_rho_types, ONLY: qs_rho_get,& + qs_rho_type +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_nonscf_utils' + + PUBLIC :: qs_nonscf_print_summary + +CONTAINS + +! ************************************************************************************************** +!> \brief writes a summary of information after diagonalization +!> \param qs_env ... +!> \param tdiag ... +!> \param nelectron_total ... +!> \param iounit ... +! ************************************************************************************************** + SUBROUTINE qs_nonscf_print_summary(qs_env, tdiag, nelectron_total, iounit) + TYPE(qs_environment_type), POINTER :: qs_env + REAL(KIND=dp), INTENT(IN) :: tdiag + INTEGER, INTENT(IN) :: nelectron_total, iounit + + REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r + TYPE(dft_control_type), POINTER :: dft_control + TYPE(qs_charges_type), POINTER :: qs_charges + TYPE(qs_energy_type), POINTER :: energy + TYPE(qs_rho_type), POINTER :: rho + + IF (iounit > 0) THEN + CALL get_qs_env(qs_env=qs_env, energy=energy, dft_control=dft_control) + IF (qs_env%harris_method) THEN + CPASSERT(.NOT. dft_control%qs_control%gapw) + CPASSERT(.NOT. dft_control%qs_control%gapw_xc) + + CALL get_qs_env(qs_env=qs_env, rho=rho, qs_charges=qs_charges) + CALL qs_rho_get(rho, tot_rho_r=tot_rho_r) + WRITE (UNIT=iounit, FMT="(/,(T3,A,T41,2F20.10))") & + "Electronic density on regular grids: ", & + SUM(tot_rho_r), & + SUM(tot_rho_r) + nelectron_total, & + "Core density on regular grids:", & + qs_charges%total_rho_core_rspace, & + qs_charges%total_rho_core_rspace - REAL(nelectron_total + dft_control%charge, dp) + WRITE (UNIT=iounit, FMT="(T3,A,T41,F20.10)") & + "Total charge density on r-space grids: ", & + SUM(tot_rho_r) + & + qs_charges%total_rho_core_rspace, & + "Total charge density g-space grids: ", & + qs_charges%total_rho_gspace + + WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & + "Diagonalization", "Time:", tdiag, energy%band + + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Core Hamiltonian energy: ", energy%core, & + "Overlap energy of the core charge distribution:", energy%core_overlap, & + "Self energy of the core charge distribution: ", energy%core_self, & + "Hartree energy: ", energy%hartree, & + "Exchange-correlation energy: ", energy%exc + IF (energy%dispersion /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Dispersion energy: ", energy%dispersion + IF (energy%gcp /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "gCP energy: ", energy%gcp + IF (energy%efield /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Electric field interaction energy: ", energy%efield + + ELSEIF (dft_control%qs_control%semi_empirical) THEN + CPABORT("NONSCF not available") + ELSEIF (dft_control%qs_control%dftb) THEN + CPASSERT(energy%dftb3 == 0.0_dp) + energy%total = energy%total + energy%band + energy%qmmm_el + WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & + "Diagonalization", "Time:", tdiag, energy%total + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Core Hamiltonian energy: ", energy%core, & + "Repulsive potential energy: ", energy%repulsive, & + "Dispersion energy: ", energy%dispersion + IF (energy%efield /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Electric field interaction energy: ", energy%efield + ELSEIF (dft_control%qs_control%xtb) THEN + energy%total = energy%total + energy%band + energy%qmmm_el + WRITE (UNIT=iounit, FMT="(/,T2,A,T40,A,F10.2,T61,F20.10)") & + "Diagonalization", "Time:", tdiag, energy%total + CPASSERT(dft_control%qs_control%xtb_control%gfn_type == 0) + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Core Hamiltonian energy: ", energy%core, & + "Repulsive potential energy: ", energy%repulsive, & + "SRB Correction energy: ", energy%srb, & + "Charge equilibration energy: ", energy%eeq, & + "Dispersion energy: ", energy%dispersion + IF (dft_control%qs_control%xtb_control%do_nonbonded) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Correction for nonbonded interactions: ", energy%xtb_nonbonded + IF (energy%efield /= 0.0_dp) & + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "Electric field interaction energy: ", energy%efield + ELSE + CPABORT("NONSCF not available") + END IF + IF (dft_control%smear) THEN + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Electronic entropic energy:", energy%kTS + WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") & + "Fermi energy:", energy%efermi + END IF + IF (energy%qmmm_el /= 0.0_dp) THEN + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "QM/MM Electrostatic energy: ", energy%qmmm_el + IF (qs_env%qmmm_env_qm%image_charge) THEN + WRITE (UNIT=iounit, FMT="(T3,A,T56,F25.14)") & + "QM/MM image charge energy: ", energy%image_charge + END IF + END IF + + WRITE (UNIT=iounit, FMT="(/,(T3,A,T56,F25.14))") & + "Total energy: ", energy%total + + CALL m_flush(iounit) + END IF + + END SUBROUTINE qs_nonscf_print_summary + +END MODULE qs_nonscf_utils diff --git a/tests/QS/regtest-harris/TEST_FILES b/tests/QS/regtest-harris/TEST_FILES index 66feaddf79..7050f0d0a0 100644 --- a/tests/QS/regtest-harris/TEST_FILES +++ b/tests/QS/regtest-harris/TEST_FILES @@ -3,4 +3,5 @@ # see regtest/TEST_FILES 2H2O_t01.inp 11 1e-010 -34.570515851298339 ft1.inp 0 +fl1.inp 0 #EOF diff --git a/tests/QS/regtest-harris/fl1.inp b/tests/QS/regtest-harris/fl1.inp new file mode 100644 index 0000000000..b4832515a3 --- /dev/null +++ b/tests/QS/regtest-harris/fl1.inp @@ -0,0 +1,79 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT cc02 + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + CHECK_ATOM_FORCE 1 xyz + DE 0.001 + DEBUG_DIPOLE .FALSE. + DEBUG_FORCES .TRUE. + DEBUG_POLARIZABILITY .FALSE. + DEBUG_STRESS_TENSOR .FALSE. + EPS_NO_ERROR_CHECK 1.E-4 + MAX_RELATIVE_ERROR 0.05 + STOP_ON_MISMATCH T +&END DEBUG + +&FORCE_EVAL + METHOD Quickstep + &DFT + BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH + POTENTIAL_FILE_NAME POTENTIAL_UZH + &HARRIS_METHOD ON + &END HARRIS_METHOD + &LS_SCF + EPS_FILTER 1.0E-7 + MATRIX_CLUSTER_TYPE ATOMIC + MU -0.20 + NON_MONOTONIC T + PURIFICATION_METHOD TC2 + REPORT_ALL_SPARSITIES OFF + S_PRECONDITIONER MOLECULAR + &END LS_SCF + &MGRID + CUTOFF 300 + REL_CUTOFF 60 + &END MGRID + &QS + EPS_DEFAULT 1.E-14 + LS_SCF + &END QS + &XC + &VDW_POTENTIAL + DISPERSION_FUNCTIONAL PAIR_POTENTIAL + &PAIR_POTENTIAL + PARAMETER_FILE_NAME dftd3.dat + REFERENCE_FUNCTIONAL PBE + TYPE DFTD3(BJ) + &END PAIR_POTENTIAL + &END VDW_POTENTIAL + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &END DFT + &SUBSYS + &CELL + ABC 6.0 6.0 6.0 + &END CELL + &COORD + O 0.094933 -0.000368 0.895642 + C -0.031077 -0.000121 -0.307326 + H -0.090437 0.947608 -0.895642 + H -0.094933 -0.947608 -0.895562 + &END COORD + &KIND H + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q1 + POTENTIAL GTH-GGA-q1 + &END KIND + &KIND C + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q4 + POTENTIAL GTH-GGA-q4 + &END KIND + &KIND O + BASIS_SET ORB DZVP-MOLOPT-GGA-GTH-q6 + POTENTIAL GTH-GGA-q6 + &END KIND + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/xTB/regtest-gfn0/TEST_FILES b/tests/xTB/regtest-gfn0/TEST_FILES index 76d427399a..8c4ca70330 100644 --- a/tests/xTB/regtest-gfn0/TEST_FILES +++ b/tests/xTB/regtest-gfn0/TEST_FILES @@ -19,4 +19,6 @@ SiC-stress.inp 0 ch2o_eeq.inp 1 1.0E-11 -7.84498782090161 ch2o_eloc.inp 0 ch2o_eper.inp 0 +ch2o_ls1.inp 1 1.0E-12 -6.46515972056537 +ch2o_ls2.inp 0 #EOF diff --git a/tests/xTB/regtest-gfn0/ch2o_ls1.inp b/tests/xTB/regtest-gfn0/ch2o_ls1.inp new file mode 100644 index 0000000000..708db197e4 --- /dev/null +++ b/tests/xTB/regtest-gfn0/ch2o_ls1.inp @@ -0,0 +1,38 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT CH2O + RUN_TYPE ENERGY +&END GLOBAL + +&FORCE_EVAL + &DFT + &LS_SCF + EPS_FILTER 1.0E-7 + MATRIX_CLUSTER_TYPE ATOMIC + MU -0.20 + NON_MONOTONIC T + PURIFICATION_METHOD TC2 + REPORT_ALL_SPARSITIES OFF + S_PRECONDITIONER MOLECULAR + &END LS_SCF + &QS + LS_SCF + METHOD xTB + &XTB + GFN_TYPE 0 + &END XTB + &END QS + &END DFT + &SUBSYS + &CELL + ABC 20.0 20.0 20.0 + PERIODIC NONE + &END CELL + &COORD + O 0.051368 0.000000 0.000000 + C 1.278612 0.000000 0.000000 + H 1.870460 0.939607 0.000000 + H 1.870460 -0.939607 0.000000 + &END COORD + &END SUBSYS +&END FORCE_EVAL diff --git a/tests/xTB/regtest-gfn0/ch2o_ls2.inp b/tests/xTB/regtest-gfn0/ch2o_ls2.inp new file mode 100644 index 0000000000..6acecdba69 --- /dev/null +++ b/tests/xTB/regtest-gfn0/ch2o_ls2.inp @@ -0,0 +1,47 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT CH2O + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DEBUG_FORCES T + DEBUG_STRESS_TENSOR F + DX 0.0001 + EPS_NO_ERROR_CHECK 0.000001 + STOP_ON_MISMATCH T +&END DEBUG + +&FORCE_EVAL + &DFT + &LS_SCF + EPS_FILTER 1.0E-7 + MATRIX_CLUSTER_TYPE ATOMIC + MU -0.20 + NON_MONOTONIC T + PURIFICATION_METHOD TC2 + REPORT_ALL_SPARSITIES OFF + S_PRECONDITIONER MOLECULAR + &END LS_SCF + &QS + LS_SCF + METHOD xTB + &XTB + GFN_TYPE 0 + VDW_POTENTIAL DFTD3 + &END XTB + &END QS + &END DFT + &SUBSYS + &CELL + ABC 20.0 20.0 20.0 + PERIODIC NONE + &END CELL + &COORD + O 0.051368 0.000000 0.000000 + C 1.278612 0.000000 0.000000 + H 1.870460 0.939607 0.000000 + H 1.870460 -0.939607 0.000000 + &END COORD + &END SUBSYS +&END FORCE_EVAL