From 67ae0a34c14263a1356da32c81145dce161bea2e Mon Sep 17 00:00:00 2001 From: Frederick Stein <43850145+fstein93@users.noreply.github.com> Date: Fri, 18 Jul 2025 09:53:51 +0200 Subject: [PATCH] Move calls to GreenX to a separate module (#4307) --- src/CMakeLists.txt | 1 + src/emd/rt_bse_io.F | 182 ++------------------ src/emd/rt_bse_types.F | 8 +- src/greenx_interface.F | 383 +++++++++++++++++++++++++++++++++++++++++ src/mp2_grids.F | 329 +++++++++++++++++------------------ src/rpa_main.F | 301 ++++++++++++-------------------- 6 files changed, 662 insertions(+), 542 deletions(-) create mode 100644 src/greenx_interface.F diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 907745ac45..46bec2e369 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -90,6 +90,7 @@ list( constraint_fxd.F constraint_util.F constraint_vsite.F + greenx_interface.F core_ae.F core_ppl.F core_ppnl.F diff --git a/src/emd/rt_bse_io.F b/src/emd/rt_bse_io.F index e31694e245..0f7d870839 100644 --- a/src/emd/rt_bse_io.F +++ b/src/emd/rt_bse_io.F @@ -50,12 +50,7 @@ MODULE rt_bse_io USE physcon, ONLY: femtoseconds, & evolt USE mathconstants, ONLY: twopi -#if defined (__GREENX) - USE gx_ac, ONLY: create_thiele_pade, & - evaluate_thiele_pade_at, & - free_params, & - params -#endif + USE greenx_interface, ONLY: greenx_refine_pade, greenx_output_polarizability, greenx_refine_ft #include "../base/base_uses.f90" @@ -584,11 +579,6 @@ CONTAINS ! - first dimension is 6 - 1 - real part along x, 2 - imag part along x, 3 - real part along y, ... REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: ft_full_series INTEGER :: i, n, ft_unit -#if defined (__GREENX) - COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_complex, & - moments_ft_complex - COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: moments_eval_complex -#endif logger => cp_get_default_logger() @@ -632,124 +622,16 @@ CONTAINS CALL cp_print_key_finished_output(ft_unit, logger, rtbse_env%ft_section) DEALLOCATE (ft_real_series) DEALLOCATE (ft_imag_series) -#if defined (__GREENX) IF (rtbse_env%pade_requested) THEN - ! Report Padé refinement - IF (rtbse_env%unit_nr > 0) WRITE (rtbse_env%unit_nr, '(A10,A27,E23.8E3,E20.8E3)') & - " PADE_FT| ", "Evaluation grid bounds [eV]", rtbse_env%pade_e_min, rtbse_env%pade_e_max - ALLOCATE (omega_complex(n)) - ALLOCATE (moments_ft_complex(n)) - ALLOCATE (moments_eval_complex(3, rtbse_env%pade_npoints)) - omega_complex(:) = CMPLX(omega_series(:), 0.0, kind=dp) - DO i = 1, 3 - moments_ft_complex(:) = CMPLX(ft_full_series(2*i - 1, :), & - ft_full_series(2*i, :), & - kind=dp) - ! Copy the fitting parameters - ! TODO : Optional direct setting of parameters? - CALL refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, & - omega_complex, moments_ft_complex, rtbse_env%pade_x_eval, moments_eval_complex(i, :)) - END DO - ! Write into alternative file - ft_unit = cp_print_key_unit_nr(logger, rtbse_env%ft_section, extension="_PADE.dat", & - file_form="FORMATTED", file_position="REWIND") - IF (ft_unit > 0) THEN - DO i = 1, rtbse_env%pade_npoints - WRITE (ft_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & - REAL(rtbse_env%pade_x_eval(i)), REAL(moments_eval_complex(1, i)), AIMAG(moments_eval_complex(1, i)), & - REAL(moments_eval_complex(2, i)), AIMAG(moments_eval_complex(2, i)), & - REAL(moments_eval_complex(3, i)), AIMAG(moments_eval_complex(3, i)) - END DO - END IF - CALL cp_print_key_finished_output(ft_unit, logger, rtbse_env%ft_section) - DEALLOCATE (omega_complex) - DEALLOCATE (moments_ft_complex) - DEALLOCATE (moments_eval_complex) + CALL greenx_refine_pade(rtbse_env%pade_e_min, rtbse_env%pade_e_max, rtbse_env%pade_x_eval, & + n, rtbse_env%pade_npoints, logger, rtbse_env%ft_section, & + rtbse_env%unit_nr, omega_series, ft_full_series) END IF -#endif DEALLOCATE (omega_series) DEALLOCATE (ft_full_series) ! Perform control with gx_ac END SUBROUTINE output_moments_ft ! ************************************************************************************************** -!> \brief Refines the FT grid using Padé approximants -!> \param x_fit Input x-variables -!> \param y_fit Input y-variables -!> \param x_eval Refined x-variables -!> \param y_eval Refined y-variables -! ************************************************************************************************** - SUBROUTINE refine_ft(fit_e_min, fit_e_max, x_fit, y_fit, x_eval, y_eval, n_pade_opt) - REAL(kind=dp) :: fit_e_min, & - fit_e_max - COMPLEX(kind=dp), DIMENSION(:) :: x_fit, & - y_fit, & - x_eval, & - y_eval - INTEGER, OPTIONAL :: n_pade_opt -#if defined (__GREENX) - CHARACTER(len=*), PARAMETER :: routineN = "refine_ft" - INTEGER :: handle, & - fit_start, & - fit_end, & - max_fit, & - n_fit, & - n_pade, & - n_eval, & - i - TYPE(params) :: pade_params - CALL timeset(routineN, handle) - - ! Get the sizes from arrays - max_fit = SIZE(x_fit) - n_eval = SIZE(x_eval) - - ! Search for the fit start and end indices - fit_start = -1 - fit_end = -1 - ! Search for the subset of FT points which is within energy limits given by - ! the input - ! Do not search when automatic request of highest energy is made - IF (fit_e_max < 0) fit_end = max_fit - DO i = 1, max_fit - IF (fit_start == -1 .AND. REAL(x_fit(i)) >= fit_e_min) fit_start = i - IF (fit_end == -1 .AND. REAL(x_fit(i)) > fit_e_max) fit_end = i - 1 - IF (fit_start > 0 .AND. fit_end > 0) EXIT - END DO - IF (fit_start == -1) fit_start = 1 - IF (fit_end == -1) fit_end = max_fit - n_fit = fit_end - fit_start + 1 - - n_pade = n_fit/2 - IF (PRESENT(n_pade_opt)) n_pade = n_pade_opt - - ! Warn about a large number of Padé parameters - IF (n_pade > 1000) THEN - CPWARN("More then 1000 Padé parameters requested - may reduce with FIT_E_MIN/FIT_E_MAX.") - END IF - ! TODO : Symmetry mode settable? - ! Here, we assume that ft corresponds to transform of real trace - pade_params = create_thiele_pade(n_pade, x_fit(fit_start:fit_end), y_fit(fit_start:fit_end), & - enforce_symmetry="conjugate") - - ! Check whetner the splice is needed or not - y_eval(1:n_eval) = evaluate_thiele_pade_at(pade_params, x_eval) - - CALL free_params(pade_params) - - CALL timestop(handle) -#else - ! Mark used - MARK_USED(fit_e_min) - MARK_USED(fit_e_max) - MARK_USED(x_fit) - MARK_USED(y_fit) - MARK_USED(x_eval) - MARK_USED(y_eval) - MARK_USED(n_pade_opt) - CPABORT("refine_ft called but CP2K compiled without GreenX - GreenX needed.") -#endif - END SUBROUTINE refine_ft -! ************************************************************************************************** !> \brief Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega), !> where i and j are provided by the configuration. The tensor element is energy dependent and !> has real and imaginary parts @@ -842,11 +724,11 @@ CONTAINS rtbse_env%dft_control%rtp_control%delta_pulse_direction(rtbse_env%pol_elements(k, 2))* & rtbse_env%dft_control%rtp_control%delta_pulse_scale, 0.0, kind=dp) ELSE - CALL refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, omega_complex, field_series, & - rtbse_env%pade_x_eval, field_refined) + CALL greenx_refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, omega_complex, field_series, & + rtbse_env%pade_x_eval, field_refined) END IF - CALL refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, omega_complex, moment_series, & - rtbse_env%pade_x_eval, moment_refined) + CALL greenx_refine_ft(rtbse_env%pade_fit_e_min, rtbse_env%pade_fit_e_max, omega_complex, moment_series, & + rtbse_env%pade_x_eval, moment_refined) polarizability_refined(k, :) = moment_refined(:)/field_refined(:) END IF END DO @@ -897,51 +779,11 @@ CONTAINS END DO CALL cp_print_key_finished_output(pol_unit, logger, rtbse_env%pol_section) END IF -#if defined(__GREENX) - ! Print out the refined polarizability to a file - pol_unit = cp_print_key_unit_nr(logger, rtbse_env%pol_section, extension="_PADE.dat", & - file_form="FORMATTED", file_position="REWIND") - ! Printing for both the stdout and separate file - IF (pol_unit > 0 .AND. rtbse_env%pade_requested) THEN - IF (pol_unit == rtbse_env%unit_nr) THEN - ! Print the stdout preline - WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|" - ELSE - ! Print also the energy in atomic units - WRITE (pol_unit, '(A1,A19)', advance="no") "#", "omega [a.u.]" - END IF - ! Common - print the energy in eV - WRITE (pol_unit, '(A20)', advance="no") "Energy [eV]" - ! Print a header for each polarizability element - DO k = 1, n_elems - 1 - WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)', advance="no") & - "Real pol.", rtbse_env%pol_elements(k, 1), rtbse_env%pol_elements(k, 2), & - "Imag pol.", rtbse_env%pol_elements(k, 1), rtbse_env%pol_elements(k, 2) - END DO - WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)') & - "Real pol.", rtbse_env%pol_elements(n_elems, 1), rtbse_env%pol_elements(n_elems, 2), & - "Imag pol.", rtbse_env%pol_elements(n_elems, 1), rtbse_env%pol_elements(n_elems, 2) - DO i = 1, SIZE(rtbse_env%pade_x_eval) - IF (pol_unit == rtbse_env%unit_nr) THEN - ! Print the stdout preline - WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|" - ELSE - ! omega in a.u. - WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(rtbse_env%pade_x_eval(i), kind=dp) - END IF - ! Common values - WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(rtbse_env%pade_x_eval(i), kind=dp)*evolt - DO k = 1, n_elems - 1 - WRITE (pol_unit, '(E20.8E3,E20.8E3)', advance="no") & - REAL(polarizability_refined(k, i)), AIMAG(polarizability_refined(k, i)) - END DO - ! Print the final value and advance - WRITE (pol_unit, '(E20.8E3,E20.8E3)') & - REAL(polarizability_refined(n_elems, i)), AIMAG(polarizability_refined(n_elems, i)) - END DO - CALL cp_print_key_finished_output(pol_unit, logger, rtbse_env%pol_section) + + IF (rtbse_env%pade_requested) THEN + CALL greenx_output_polarizability(logger, rtbse_env%pol_section, rtbse_env%unit_nr, & + rtbse_env%pol_elements, rtbse_env%pade_x_eval, polarizability_refined) END IF -#endif DEALLOCATE (value_series) DEALLOCATE (ft_real_series) diff --git a/src/emd/rt_bse_types.F b/src/emd/rt_bse_types.F index b7d038edfe..d0edccc2dd 100644 --- a/src/emd/rt_bse_types.F +++ b/src/emd/rt_bse_types.F @@ -468,15 +468,11 @@ CONTAINS rtbse_env%pade_npoints = INT((rtbse_env%pade_e_max - rtbse_env%pade_e_min)/rtbse_env%pade_e_step) ! Evaluate the evaluation grid IF (rtbse_env%pade_requested) THEN -#if defined(__GREENX) NULLIFY (rtbse_env%pade_x_eval) ALLOCATE (rtbse_env%pade_x_eval(rtbse_env%pade_npoints)) DO i = 1, rtbse_env%pade_npoints rtbse_env%pade_x_eval(i) = CMPLX(rtbse_env%pade_e_step*REAL(i - 1, kind=dp), 0.0, kind=dp) END DO -#else - CPABORT("Padé refinement of FT requires linking GreenX - recompile CP2K with GreenX.") -#endif END IF ! Allocate moments matrices @@ -668,9 +664,7 @@ CONTAINS DEALLOCATE (rtbse_env%time_trace%series) DEALLOCATE (rtbse_env%pol_elements) -#if defined(__GREENX) - IF (rtbse_env%pade_requested) DEALLOCATE (rtbse_env%pade_x_eval) -#endif + IF (ASSOCIATED(rtbse_env%pade_x_eval)) DEALLOCATE (rtbse_env%pade_x_eval) ! Deallocate the neighbour list that is not deallocated in gw anymore IF (ASSOCIATED(rtbse_env%bs_env%nl_3c%ij_list)) CALL neighbor_list_3c_destroy(rtbse_env%bs_env%nl_3c) diff --git a/src/greenx_interface.F b/src/greenx_interface.F new file mode 100644 index 0000000000..9ddd775633 --- /dev/null +++ b/src/greenx_interface.F @@ -0,0 +1,383 @@ +!--------------------------------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations ! +! Copyright 2000-2025 CP2K developers group ! +! ! +! SPDX-License-Identifier: GPL-2.0-or-later ! +!--------------------------------------------------------------------------------------------------! + +! ************************************************************************************************** +!> \brief Interface to the Greenx library +!> \par History +!> 07.2025 Refactored from RPA and BSE modules [Frederick Stein] +! ************************************************************************************************** +MODULE greenx_interface + USE kinds, ONLY: dp + USE cp_log_handling, ONLY: cp_logger_type, & + cp_get_default_logger + USE cp_output_handling, ONLY: cp_print_key_unit_nr, & + cp_print_key_finished_output, & + cp_print_key_generate_filename, & + low_print_level, & + medium_print_level + USE input_section_types, ONLY: section_vals_type + USE machine, ONLY: m_flush + USE physcon, ONLY: evolt +#if defined (__GREENX) + USE gx_ac, ONLY: create_thiele_pade, & + evaluate_thiele_pade_at, & + free_params, & + params + USE gx_minimax, ONLY: gx_minimax_grid +#endif + +#include "./base/base_uses.f90" + + IMPLICIT NONE + + PRIVATE + + CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'greenx_interface' + + PUBLIC :: greenx_refine_pade, greenx_output_polarizability, greenx_refine_ft, greenx_get_minimax_grid + +CONTAINS + +! ************************************************************************************************** +!> \brief Refines Pade approximants using GreenX, skips this step if GreenX is not available +!> \param e_min ... +!> \param e_max ... +!> \param x_eval ... +!> \param number_of_simulation_steps ... +!> \param number_of_pade_points ... +!> \param logger ... +!> \param ft_section ... +!> \param bse_unit ... +!> \param omega_series ... +!> \param ft_full_series ... +! ************************************************************************************************** + SUBROUTINE greenx_refine_pade(e_min, e_max, x_eval, number_of_simulation_steps, number_of_pade_points, & + logger, ft_section, bse_unit, omega_series, ft_full_series) + REAL(KIND=dp), INTENT(IN) :: e_min, e_max + COMPLEX(KIND=dp), DIMENSION(:), POINTER :: x_eval + INTEGER, INTENT(IN) :: number_of_simulation_steps, number_of_pade_points + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), POINTER :: ft_section + INTEGER, INTENT(IN) :: bse_unit + REAL(KIND=dp), DIMENSION(number_of_simulation_steps + 2), INTENT(INOUT) :: omega_series + REAL(KIND=dp), DIMENSION(6, number_of_simulation_steps + 2), INTENT(INOUT) :: ft_full_series +#if defined (__GREENX) + INTEGER :: i, ft_unit + COMPLEX(kind=dp), DIMENSION(:), ALLOCATABLE :: omega_complex, & + moments_ft_complex + COMPLEX(kind=dp), DIMENSION(:, :), ALLOCATABLE :: moments_eval_complex + + ! Report Padé refinement + IF (bse_unit > 0) WRITE (bse_unit, '(A10,A27,E23.8E3,E20.8E3)') & + " PADE_FT| ", "Evaluation grid bounds [eV]", e_min, e_max + ALLOCATE (omega_complex(number_of_simulation_steps + 2)) + ALLOCATE (moments_ft_complex(number_of_simulation_steps + 2)) + ALLOCATE (moments_eval_complex(3, number_of_pade_points)) + omega_complex(:) = CMPLX(omega_series(:), 0.0, kind=dp) + DO i = 1, 3 + moments_ft_complex(:) = CMPLX(ft_full_series(2*i - 1, :), & + ft_full_series(2*i, :), & + kind=dp) + ! Copy the fitting parameters + ! TODO : Optional direct setting of parameters? + CALL greenx_refine_ft(e_min, e_max, omega_complex, moments_ft_complex, x_eval, moments_eval_complex(i, :)) + END DO + ! Write into alternative file + ft_unit = cp_print_key_unit_nr(logger, ft_section, extension="_PADE.dat", & + file_form="FORMATTED", file_position="REWIND") + IF (ft_unit > 0) THEN + DO i = 1, number_of_pade_points + WRITE (ft_unit, '(E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3,E20.8E3)') & + REAL(x_eval(i)), REAL(moments_eval_complex(1, i)), AIMAG(moments_eval_complex(1, i)), & + REAL(moments_eval_complex(2, i)), AIMAG(moments_eval_complex(2, i)), & + REAL(moments_eval_complex(3, i)), AIMAG(moments_eval_complex(3, i)) + END DO + END IF + CALL cp_print_key_finished_output(ft_unit, logger, ft_section) + DEALLOCATE (omega_complex) + DEALLOCATE (moments_ft_complex) + DEALLOCATE (moments_eval_complex) +#else + IF (bse_unit > 0) WRITE (bse_unit, '(A10,A70)') & + " PADE_FT| ", "GreenX library is not available. Refinement is skipped" + MARK_USED(e_min) + MARK_USED(e_max) + MARK_USED(x_eval) + MARK_USED(number_of_simulation_steps) + MARK_USED(number_of_pade_points) + MARK_USED(logger) + MARK_USED(ft_section) + MARK_USED(omega_series) + MARK_USED(ft_full_series) +#endif + END SUBROUTINE greenx_refine_pade +! ************************************************************************************************** +!> \brief Outputs the isotropic polarizability tensor element alpha _ ij = mu_i(omega)/E_j(omega), +!> where i and j are provided by the configuration. The tensor element is energy dependent and +!> has real and imaginary parts +!> \param logger ... +!> \param pol_section ... +!> \param bse_unit ... +!> \param pol_elements ... +!> \param x_eval ... +!> \param polarizability_refined ... +! ************************************************************************************************** + SUBROUTINE greenx_output_polarizability(logger, pol_section, bse_unit, pol_elements, x_eval, polarizability_refined) + TYPE(cp_logger_type), POINTER :: logger + TYPE(section_vals_type), POINTER :: pol_section + INTEGER, INTENT(IN) :: bse_unit + INTEGER, DIMENSION(:, :), POINTER :: pol_elements + COMPLEX(KIND=dp), DIMENSION(:), POINTER :: x_eval + COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: polarizability_refined +#if defined(__GREENX) + INTEGER :: pol_unit, & + i, k, n_elems + + n_elems = SIZE(pol_elements, 1) + ! Print out the refined polarizability to a file + pol_unit = cp_print_key_unit_nr(logger, pol_section, extension="_PADE.dat", & + file_form="FORMATTED", file_position="REWIND") + ! Printing for both the stdout and separate file + IF (pol_unit > 0) THEN + IF (pol_unit == bse_unit) THEN + ! Print the stdout preline + WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|" + ELSE + ! Print also the energy in atomic units + WRITE (pol_unit, '(A1,A19)', advance="no") "#", "omega [a.u.]" + END IF + ! Common - print the energy in eV + WRITE (pol_unit, '(A20)', advance="no") "Energy [eV]" + ! Print a header for each polarizability element + DO k = 1, n_elems - 1 + WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)', advance="no") & + "Real pol.", pol_elements(k, 1), pol_elements(k, 2), & + "Imag pol.", pol_elements(k, 1), pol_elements(k, 2) + END DO + WRITE (pol_unit, '(A16,I2,I2,A16,I2,I2)') & + "Real pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2), & + "Imag pol.", pol_elements(n_elems, 1), pol_elements(n_elems, 2) + DO i = 1, SIZE(x_eval) + IF (pol_unit == bse_unit) THEN + ! Print the stdout preline + WRITE (pol_unit, '(A21)', advance="no") " POLARIZABILITY_PADE|" + ELSE + ! omega in a.u. + WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(x_eval(i), kind=dp) + END IF + ! Common values + WRITE (pol_unit, '(E20.8E3)', advance="no") REAL(x_eval(i), kind=dp)*evolt + DO k = 1, n_elems - 1 + WRITE (pol_unit, '(E20.8E3,E20.8E3)', advance="no") & + REAL(polarizability_refined(k, i)), AIMAG(polarizability_refined(k, i)) + END DO + ! Print the final value and advance + WRITE (pol_unit, '(E20.8E3,E20.8E3)') & + REAL(polarizability_refined(n_elems, i)), AIMAG(polarizability_refined(n_elems, i)) + END DO + CALL cp_print_key_finished_output(pol_unit, logger, pol_section) + END IF +#else + MARK_USED(logger) + MARK_USED(pol_section) + MARK_USED(bse_unit) + MARK_USED(pol_elements) + MARK_USED(x_eval) + MARK_USED(polarizability_refined) +#endif + END SUBROUTINE greenx_output_polarizability +! ************************************************************************************************** +!> \brief Refines the FT grid using Padé approximants +!> \param fit_e_min ... +!> \param fit_e_max ... +!> \param x_fit Input x-variables +!> \param y_fit Input y-variables +!> \param x_eval Refined x-variables +!> \param y_eval Refined y-variables +!> \param n_pade_opt ... +! ************************************************************************************************** + SUBROUTINE greenx_refine_ft(fit_e_min, fit_e_max, x_fit, y_fit, x_eval, y_eval, n_pade_opt) + REAL(kind=dp) :: fit_e_min, & + fit_e_max + COMPLEX(kind=dp), DIMENSION(:) :: x_fit, & + y_fit, & + x_eval, & + y_eval + INTEGER, OPTIONAL :: n_pade_opt +#if defined (__GREENX) + INTEGER :: fit_start, & + fit_end, & + max_fit, & + n_fit, & + n_pade, & + n_eval, & + i + TYPE(params) :: pade_params + + ! Get the sizes from arrays + max_fit = SIZE(x_fit) + n_eval = SIZE(x_eval) + + ! Search for the fit start and end indices + fit_start = -1 + fit_end = -1 + ! Search for the subset of FT points which is within energy limits given by + ! the input + ! Do not search when automatic request of highest energy is made + IF (fit_e_max < 0) fit_end = max_fit + DO i = 1, max_fit + IF (fit_start == -1 .AND. REAL(x_fit(i)) >= fit_e_min) fit_start = i + IF (fit_end == -1 .AND. REAL(x_fit(i)) > fit_e_max) fit_end = i - 1 + IF (fit_start > 0 .AND. fit_end > 0) EXIT + END DO + IF (fit_start == -1) fit_start = 1 + IF (fit_end == -1) fit_end = max_fit + n_fit = fit_end - fit_start + 1 + + n_pade = n_fit/2 + IF (PRESENT(n_pade_opt)) n_pade = n_pade_opt + + ! Warn about a large number of Padé parameters + IF (n_pade > 1000) THEN + CPWARN("More then 1000 Padé parameters requested - may reduce with FIT_E_MIN/FIT_E_MAX.") + END IF + ! TODO : Symmetry mode settable? + ! Here, we assume that ft corresponds to transform of real trace + pade_params = create_thiele_pade(n_pade, x_fit(fit_start:fit_end), y_fit(fit_start:fit_end), & + enforce_symmetry="conjugate") + + ! Check whetner the splice is needed or not + y_eval(1:n_eval) = evaluate_thiele_pade_at(pade_params, x_eval) + + CALL free_params(pade_params) +#else + ! Mark used + MARK_USED(fit_e_min) + MARK_USED(fit_e_max) + MARK_USED(x_fit) + MARK_USED(y_fit) + MARK_USED(x_eval) + MARK_USED(y_eval) + MARK_USED(n_pade_opt) + CPABORT("Calls to GreenX require CP2K to be compiled with support for GreenX.") +#endif + END SUBROUTINE + +! ************************************************************************************************** +!> \brief ... +!> \param unit_nr ... +!> \param num_integ_points ... +!> \param emin ... +!> \param emax ... +!> \param tau_tj ... +!> \param tau_wj ... +!> \param regularization_minimax ... +!> \param tj ... +!> \param wj ... +!> \param weights_cos_tf_t_to_w ... +!> \param weights_cos_tf_w_to_t ... +!> \param weights_sin_tf_t_to_w ... +!> \param ierr ... +! ************************************************************************************************** + SUBROUTINE greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, & + tau_tj, tau_wj, regularization_minimax, & + tj, wj, weights_cos_tf_t_to_w, & + weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr) + + INTEGER, INTENT(IN) :: unit_nr, num_integ_points + REAL(KIND=dp), INTENT(IN) :: emin, emax + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & + INTENT(OUT) :: tau_tj, tau_wj + REAL(KIND=dp), INTENT(IN) :: regularization_minimax + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & + INTENT(INOUT) :: tj, wj + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & + INTENT(OUT) :: weights_cos_tf_t_to_w, & + weights_cos_tf_w_to_t, & + weights_sin_tf_t_to_w + INTEGER, INTENT(OUT) :: ierr +#if defined (__GREENX) + INTEGER :: gi + REAL(KIND=dp) :: cosft_duality_error_greenx, & + max_errors_greenx(3) + + CALL gx_minimax_grid(num_integ_points, Emin, Emax, tau_tj, tau_wj, tj, wj, & + weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, & + max_errors_greenx, cosft_duality_error_greenx, ierr, & + bare_cos_sin_weights=.TRUE., & + regularization=regularization_minimax) + ! Factor 4 is hard-coded in the RPA weights in the internal CP2K minimax routines + wj(:) = wj(:)*4.0_dp + IF (ierr == 0) THEN + IF (unit_nr > 0) THEN + WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & + "GREENX MINIMAX_INFO| Number of integration points:", num_integ_points + WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") & + "GREENX MINIMAX_INFO| Gap (Emin):", Emin + WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") & + "GREENX MINIMAX_INFO| Maximum eigenvalue difference (Emax):", Emax + WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.4)") & + "GREENX MINIMAX_INFO| Energy range (Emax/Emin):", Emax/Emin + WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") & + "GREENX MINIMAX_INFO| Frequency grid (scaled):", "Weights", "Abscissas" + DO gi = 1, num_integ_points + WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(gi), tj(gi) + END DO + WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") & + "GREENX MINIMAX_INFO| Time grid (scaled):", "Weights", "Abscissas" + DO gi = 1, num_integ_points + WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(gi), tau_tj(gi) + END DO + CALL m_flush(unit_nr) + END IF + ELSE + IF (unit_nr > 0) THEN + WRITE (UNIT=unit_nr, FMT="(T3,A,T75)") & + "GREENX MINIMAX_INFO| Grid not available, use internal CP2K grids" + CALL m_flush(unit_nr) + END IF + IF (ALLOCATED(tau_tj)) THEN + DEALLOCATE (tau_tj) + END IF + IF (ALLOCATED(tau_wj)) THEN + DEALLOCATE (tau_wj) + END IF + IF (ALLOCATED(tj)) THEN + DEALLOCATE (tj) + END IF + IF (ALLOCATED(wj)) THEN + DEALLOCATE (wj) + END IF + IF (ALLOCATED(weights_cos_tf_t_to_w)) THEN + DEALLOCATE (weights_cos_tf_t_to_w) + END IF + IF (ALLOCATED(weights_cos_tf_w_to_t)) THEN + DEALLOCATE (weights_cos_tf_w_to_t) + END IF + IF (ALLOCATED(weights_sin_tf_t_to_w)) THEN + DEALLOCATE (weights_sin_tf_t_to_w) + END IF + END IF +#else + ierr = 1 + MARK_USED(unit_nr) + MARK_USED(num_integ_points) + MARK_USED(emin) + MARK_USED(emax) + MARK_USED(tau_tj) + MARK_USED(tau_wj) + MARK_USED(regularization_minimax) + MARK_USED(tj) + MARK_USED(wj) + MARK_USED(weights_cos_tf_t_to_w) + MARK_USED(weights_cos_tf_w_to_t) + MARK_USED(weights_sin_tf_t_to_w) +#endif + + END SUBROUTINE greenx_get_minimax_grid + +END MODULE greenx_interface diff --git a/src/mp2_grids.F b/src/mp2_grids.F index b305eab38e..5b251bc2e9 100644 --- a/src/mp2_grids.F +++ b/src/mp2_grids.F @@ -14,6 +14,7 @@ MODULE mp2_grids USE cp_fm_types, ONLY: cp_fm_get_info,& cp_fm_type + USE greenx_interface, ONLY: greenx_get_minimax_grid USE input_section_types, ONLY: section_vals_type,& section_vals_val_set USE kinds, ONLY: dp @@ -41,7 +42,7 @@ MODULE mp2_grids CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_grids' PUBLIC :: get_minimax_grid, get_clenshaw_grid, test_least_square_ft, get_l_sq_wghts_cos_tf_t_to_w, & - get_l_sq_wghts_cos_tf_w_to_t, get_l_sq_wghts_sin_tf_t_to_w, init_greenx_grids + get_l_sq_wghts_cos_tf_w_to_t, get_l_sq_wghts_sin_tf_t_to_w CONTAINS @@ -86,7 +87,7 @@ CONTAINS LOGICAL, INTENT(IN) :: do_gw_im_time, do_kpoints_cubic_RPA REAL(KIND=dp), INTENT(OUT) :: e_fermi REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & - INTENT(INOUT) :: tj, wj + INTENT(OUT) :: tj, wj REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & INTENT(OUT) :: weights_cos_tf_t_to_w, & weights_cos_tf_w_to_t, & @@ -96,18 +97,30 @@ CONTAINS CHARACTER(LEN=*), PARAMETER :: routineN = 'get_minimax_grid' INTEGER, PARAMETER :: num_points_per_magnitude = 200 - INTEGER :: handle, ierr, ispin, jquad, nspins - LOGICAL :: my_do_kpoints, my_open_shell - REAL(KIND=dp) :: Emax, Emin, max_error_min, my_E_Range, & + INTEGER :: handle, ierr, jquad, nspins + LOGICAL :: my_do_kpoints + REAL(KIND=dp) :: E_Range, Emax, Emin, max_error_min, & my_regularization, scaling REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: x_tw - TYPE(section_vals_type), POINTER :: input CALL timeset(routineN, handle) + CALL determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, & + do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi) + + CALL greenx_get_minimax_grid(unit_nr, num_integ_points, emin, emax, & + tau_tj, tau_wj, qs_env%mp2_env%ri_g0w0%regularization_minimax, & + tj, wj, weights_cos_tf_t_to_w, & + weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, ierr) + +! Shortcut if Greenx was available and successful + IF (ierr == 0) THEN + CALL timestop(handle) + RETURN + END IF + ! Test for spin unrestricted nspins = SIZE(homo) - my_open_shell = (nspins == 2) ! Test whether all necessary variables are available my_do_kpoints = .FALSE. @@ -118,177 +131,145 @@ CONTAINS my_regularization = 0.0_dp IF (PRESENT(regularization)) THEN my_regularization = regularization - END IF - IF (my_do_kpoints) THEN - CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi) - my_E_Range = Emax/Emin - ELSE - IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN - Emin = HUGE(dp) - Emax = 0.0_dp - DO ispin = 1, nspins - IF (homo(ispin) > 0) THEN - Emin = MIN(Emin, Eigenval(homo(ispin) + 1, 1, ispin) - Eigenval(homo(ispin), 1, ispin)) - Emax = MAX(Emax, MAXVAL(Eigenval(:, :, ispin)) - MINVAL(Eigenval(:, :, ispin))) - END IF - END DO - my_E_Range = Emax/Emin - qs_env%mp2_env%e_range = my_e_range - qs_env%mp2_env%e_gap = Emin - - CALL get_qs_env(qs_env, input=input) - CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=my_e_range) - CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin) - ELSE - my_E_range = qs_env%mp2_env%E_range - Emin = qs_env%mp2_env%E_gap - Emax = Emin*my_E_range + IF (num_integ_points > 20 .AND. e_range < 100.0_dp) THEN + IF (unit_nr > 0) & + CALL cp_warn(__LOCATION__, & + "You requested a large minimax grid (> 20 points) for a small minimax range R (R < 100). "// & + "That may lead to numerical "// & + "instabilities when computing minimax grid weights. You can prevent small ranges by choosing "// & + "a larger basis set with higher angular momenta or alternatively using all-electron calculations.") END IF - END IF - - IF (num_integ_points > 20 .AND. my_E_Range < 100.0_dp) THEN - IF (unit_nr > 0) & - CALL cp_warn(__LOCATION__, & - "You requested a large minimax grid (> 20 points) for a small minimax range R (R < 100). "// & - "That may lead to numerical "// & - "instabilities when computing minimax grid weights. You can prevent small ranges by choosing "// & - "a larger basis set with higher angular momenta or alternatively using all-electron calculations.") - END IF - - IF (.NOT. do_ri_sos_laplace_mp2) THEN - ALLOCATE (x_tw(2*num_integ_points)) - x_tw = 0.0_dp - ierr = 0 - IF (num_integ_points .LE. 20) THEN - CALL get_rpa_minimax_coeff(num_integ_points, my_E_Range, x_tw, ierr) - ELSE - CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, my_E_Range, x_tw) - END IF - - ALLOCATE (tj(num_integ_points)) - tj = 0.0_dp - - ALLOCATE (wj(num_integ_points)) - wj = 0.0_dp - - DO jquad = 1, num_integ_points - tj(jquad) = x_tw(jquad) - wj(jquad) = x_tw(jquad + num_integ_points) - END DO - - ! for the smaller grids, the factor of 4 is included in get_rpa_minimax_coeff for wj - IF (num_integ_points >= 26) THEN - wj(:) = wj(:)*4.0_dp - END IF - - DEALLOCATE (x_tw) - - IF (unit_nr > 0 .AND. do_print) THEN - WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & - "MINIMAX_INFO| Number of integration points:", num_integ_points - WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & - "MINIMAX_INFO| Gap for the minimax approximation:", Emin - WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & - "MINIMAX_INFO| Range for the minimax approximation:", my_E_Range - WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas" - DO jquad = 1, num_integ_points - WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(jquad), tj(jquad) - END DO - CALL m_flush(unit_nr) - END IF - - ! scale the minimax parameters - tj(:) = tj(:)*Emin - wj(:) = wj(:)*Emin - ELSE - ! When we perform SOS-MP2, we need an additional factor of 2 for the energies (compare with mp2_laplace.F) - ! We do not need weights etc. for the cosine transform - ! We do not scale Emax because it is not needed for SOS-MP2 - Emin = Emin*2.0_dp - END IF - - ! set up the minimax time grid - IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN - - ALLOCATE (x_tw(2*num_integ_points)) - x_tw = 0.0_dp - - IF (num_integ_points .LE. 20) THEN - CALL get_exp_minimax_coeff(num_integ_points, my_E_Range, x_tw) - ELSE - CALL get_exp_minimax_coeff_gw(num_integ_points, my_E_Range, x_tw) - END IF - - ! For RPA we include already a factor of two (see later steps) - scaling = 2.0_dp - IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp - - ALLOCATE (tau_tj(num_integ_points)) - tau_tj = 0.0_dp - - ALLOCATE (tau_wj(num_integ_points)) - tau_wj = 0.0_dp - - DO jquad = 1, num_integ_points - tau_tj(jquad) = x_tw(jquad)/scaling - tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling - END DO - - DEALLOCATE (x_tw) - - IF (unit_nr > 0 .AND. do_print) THEN - WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & - "MINIMAX_INFO| Range for the minimax approximation:", my_E_Range - ! For testing the gap - WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & - "MINIMAX_INFO| Gap:", Emin - WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") & - "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas" - DO jquad = 1, num_integ_points - WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad) - END DO - CALL m_flush(unit_nr) - END IF - - ! scale grid from [1,R] to [Emin,Emax] - tau_tj(:) = tau_tj(:)/Emin - tau_wj(:) = tau_wj(:)/Emin IF (.NOT. do_ri_sos_laplace_mp2) THEN - ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points)) - weights_cos_tf_t_to_w = 0.0_dp + ALLOCATE (x_tw(2*num_integ_points)) + x_tw = 0.0_dp + ierr = 0 + IF (num_integ_points .LE. 20) THEN + CALL get_rpa_minimax_coeff(num_integ_points, e_range, x_tw, ierr) + ELSE + CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, e_range, x_tw) + END IF - CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, tj, & - Emin, Emax, max_error_min, num_points_per_magnitude, & - my_regularization) + ALLOCATE (tj(num_integ_points)) + tj = 0.0_dp - ! get the weights for the cosine transform W^c(iw) -> W^c(it) - ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points)) - weights_cos_tf_w_to_t = 0.0_dp + ALLOCATE (wj(num_integ_points)) + wj = 0.0_dp - CALL get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, tj, & - Emin, Emax, max_error_min, num_points_per_magnitude, & - my_regularization) + DO jquad = 1, num_integ_points + tj(jquad) = x_tw(jquad) + wj(jquad) = x_tw(jquad + num_integ_points) + END DO - IF (do_gw_im_time) THEN + ! for the smaller grids, the factor of 4 is included in get_rpa_minimax_coeff for wj + IF (num_integ_points >= 26) THEN + wj(:) = wj(:)*4.0_dp + END IF - ! get the weights for the sine transform Sigma^sin(it) -> Sigma^sin(iw) (PRB 94, 165109 (2016), Eq. 71) - ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points)) - weights_sin_tf_t_to_w = 0.0_dp + DEALLOCATE (x_tw) - CALL get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, tj, & + IF (unit_nr > 0 .AND. do_print) THEN + WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & + "MINIMAX_INFO| Number of integration points:", num_integ_points + WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & + "MINIMAX_INFO| Gap for the minimax approximation:", Emin + WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & + "MINIMAX_INFO| Range for the minimax approximation:", e_range + WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "MINIMAX_INFO| Minimax parameters:", "Weights", "Abscissas" + DO jquad = 1, num_integ_points + WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(jquad), tj(jquad) + END DO + CALL m_flush(unit_nr) + END IF + + ! scale the minimax parameters + tj(:) = tj(:)*Emin + wj(:) = wj(:)*Emin + END IF + + ! set up the minimax time grid + IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN + + ALLOCATE (x_tw(2*num_integ_points)) + x_tw = 0.0_dp + + IF (num_integ_points .LE. 20) THEN + CALL get_exp_minimax_coeff(num_integ_points, e_range, x_tw) + ELSE + CALL get_exp_minimax_coeff_gw(num_integ_points, e_range, x_tw) + END IF + + ! For RPA we include already a factor of two (see later steps) + scaling = 2.0_dp + IF (do_ri_sos_laplace_mp2) scaling = 1.0_dp + + ALLOCATE (tau_tj(num_integ_points)) + tau_tj = 0.0_dp + + ALLOCATE (tau_wj(num_integ_points)) + tau_wj = 0.0_dp + + DO jquad = 1, num_integ_points + tau_tj(jquad) = x_tw(jquad)/scaling + tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling + END DO + + DEALLOCATE (x_tw) + + IF (unit_nr > 0 .AND. do_print) THEN + WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & + "MINIMAX_INFO| Range for the minimax approximation:", e_range + ! For testing the gap + WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") & + "MINIMAX_INFO| Gap:", Emin + WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") & + "MINIMAX_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas" + DO jquad = 1, num_integ_points + WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad) + END DO + CALL m_flush(unit_nr) + END IF + + ! scale grid from [1,R] to [Emin,Emax] + tau_tj(:) = tau_tj(:)/Emin + tau_wj(:) = tau_wj(:)/Emin + + IF (.NOT. do_ri_sos_laplace_mp2) THEN + ALLOCATE (weights_cos_tf_t_to_w(num_integ_points, num_integ_points)) + weights_cos_tf_t_to_w = 0.0_dp + + CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, tj, & Emin, Emax, max_error_min, num_points_per_magnitude, & my_regularization) - IF (unit_nr > 0) THEN - WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") & - "MINIMAX_INFO| Maximum deviation of the imag. time fit:", max_error_min + ! get the weights for the cosine transform W^c(iw) -> W^c(it) + ALLOCATE (weights_cos_tf_w_to_t(num_integ_points, num_integ_points)) + weights_cos_tf_w_to_t = 0.0_dp + + CALL get_l_sq_wghts_cos_tf_w_to_t(num_integ_points, tau_tj, weights_cos_tf_w_to_t, tj, & + Emin, Emax, max_error_min, num_points_per_magnitude, & + my_regularization) + + IF (do_gw_im_time) THEN + + ! get the weights for the sine transform Sigma^sin(it) -> Sigma^sin(iw) (PRB 94, 165109 (2016), Eq. 71) + ALLOCATE (weights_sin_tf_t_to_w(num_integ_points, num_integ_points)) + weights_sin_tf_t_to_w = 0.0_dp + + CALL get_l_sq_wghts_sin_tf_t_to_w(num_integ_points, tau_tj, weights_sin_tf_t_to_w, tj, & + Emin, Emax, max_error_min, num_points_per_magnitude, & + my_regularization) + + IF (unit_nr > 0) THEN + WRITE (UNIT=unit_nr, FMT="(T3,A,T66,ES15.2)") & + "MINIMAX_INFO| Maximum deviation of the imag. time fit:", max_error_min + END IF END IF + END IF END IF - END IF CALL timestop(handle) @@ -1538,10 +1519,11 @@ CONTAINS !> \param do_kpoints_cubic_RPA flag for cubic-scaling RPA with k-points !> \param Emin minimal eigenvalue difference (gap of the system) !> \param Emax maximal eigenvalue difference +!> \param e_range ... !> \param e_fermi Fermi level ! ************************************************************************************************** - SUBROUTINE init_greenx_grids(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, & - do_kpoints_cubic_RPA, Emin, Emax, e_fermi) + SUBROUTINE determine_energy_range(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, & + do_kpoints_cubic_RPA, Emin, Emax, e_range, e_fermi) TYPE(qs_environment_type), POINTER :: qs_env TYPE(mp_para_env_type), INTENT(IN) :: para_env @@ -1549,19 +1531,17 @@ CONTAINS REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: Eigenval LOGICAL, INTENT(IN) :: do_ri_sos_laplace_mp2, & do_kpoints_cubic_RPA - REAL(KIND=dp), INTENT(OUT) :: Emin, Emax, e_fermi + REAL(KIND=dp), INTENT(OUT) :: Emin, Emax, e_range, e_fermi - CHARACTER(LEN=*), PARAMETER :: routineN = 'init_greenx_grids' + CHARACTER(LEN=*), PARAMETER :: routineN = 'determine_energy_range' INTEGER :: handle, ispin, nspins - LOGICAL :: my_do_kpoints, my_open_shell - REAL(KIND=dp) :: my_E_Range + LOGICAL :: my_do_kpoints TYPE(section_vals_type), POINTER :: input CALL timeset(routineN, handle) ! Test for spin unrestricted nspins = SIZE(homo) - my_open_shell = (nspins == 2) ! Test whether all necessary variables are available my_do_kpoints = .FALSE. @@ -1571,7 +1551,7 @@ CONTAINS IF (my_do_kpoints) THEN CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi) - my_E_Range = Emax/Emin + E_Range = Emax/Emin ELSE IF (qs_env%mp2_env%E_range <= 1.0_dp .OR. qs_env%mp2_env%E_gap <= 0.0_dp) THEN Emin = HUGE(dp) @@ -1582,17 +1562,17 @@ CONTAINS Emax = MAX(Emax, MAXVAL(Eigenval(:, :, ispin)) - MINVAL(Eigenval(:, :, ispin))) END IF END DO - my_E_Range = Emax/Emin - qs_env%mp2_env%e_range = my_e_range + E_Range = Emax/Emin + qs_env%mp2_env%e_range = e_range qs_env%mp2_env%e_gap = Emin CALL get_qs_env(qs_env, input=input) - CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=my_e_range) + CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_RANGE", r_val=e_range) CALL section_vals_val_set(input, "DFT%XC%WF_CORRELATION%E_GAP", r_val=emin) ELSE - my_E_range = qs_env%mp2_env%E_range + E_range = qs_env%mp2_env%E_range Emin = qs_env%mp2_env%E_gap - Emax = Emin*my_E_range + Emax = Emin*E_range END IF END IF @@ -1601,6 +1581,7 @@ CONTAINS ! We do not scale Emax because it is not needed for SOS-MP2 IF (do_ri_sos_laplace_mp2) THEN Emin = Emin*2.0_dp + Emax = Emax*2.0_dp END IF CALL timestop(handle) diff --git a/src/rpa_main.F b/src/rpa_main.F index 1ddbe1d6b3..137fea15f4 100644 --- a/src/rpa_main.F +++ b/src/rpa_main.F @@ -15,110 +15,105 @@ !> 03.2019 Refactoring [Frederick Stein] ! ************************************************************************************************** MODULE rpa_main - USE bibliography, ONLY: & - Bates2013, DelBen2013, DelBen2015, Freeman1977, Gruneis2009, Ren2011, Ren2013, & - Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, cite_reference - USE bse_main, ONLY: start_bse_calculation - USE cp_blacs_env, ONLY: cp_blacs_env_create, & - cp_blacs_env_release, & - cp_blacs_env_type - USE cp_cfm_types, ONLY: cp_cfm_type - USE cp_dbcsr_api, ONLY: dbcsr_add, & - dbcsr_clear, & - dbcsr_get_info, & - dbcsr_p_type, & - dbcsr_type - USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add - 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_set_all, & - cp_fm_to_fm, & - cp_fm_type - USE dbt_api, ONLY: dbt_type - USE dgemm_counter_types, ONLY: dgemm_counter_init, & - dgemm_counter_type, & - dgemm_counter_write - USE group_dist_types, ONLY: create_group_dist, & - get_group_dist, & - group_dist_d1_type, & - maxsize, & - release_group_dist - USE hfx_types, ONLY: block_ind_type, & - hfx_compression_type - USE input_constants, ONLY: rpa_exchange_axk, & - rpa_exchange_none, & - rpa_exchange_sosex, & - sigma_none, & - wfc_mm_style_gemm - USE kinds, ONLY: dp, & - int_8 - USE kpoint_types, ONLY: kpoint_type - USE machine, ONLY: m_flush, & - m_memory - USE mathconstants, ONLY: pi, & - z_zero - USE message_passing, ONLY: mp_comm_type, & - mp_para_env_release, & - mp_para_env_type - USE minimax_exp, ONLY: check_exp_minimax_range - USE mp2_grids, ONLY: get_clenshaw_grid, & - get_minimax_grid, & - init_greenx_grids - USE mp2_laplace, ONLY: SOS_MP2_postprocessing - USE mp2_ri_grad_util, ONLY: array2fm - USE mp2_types, ONLY: mp2_type, & - three_dim_real_array, & - two_dim_int_array, & - two_dim_real_array - USE qs_environment_types, ONLY: get_qs_env, & - qs_environment_type - USE rpa_exchange, ONLY: rpa_exchange_needed_mem, & - rpa_exchange_work_type - USE rpa_grad, ONLY: rpa_grad_copy_Q, & - rpa_grad_create, & - rpa_grad_finalize, & - rpa_grad_matrix_operations, & - rpa_grad_needed_mem, & - rpa_grad_type - USE rpa_gw, ONLY: allocate_matrices_gw, & - allocate_matrices_gw_im_time, & - compute_GW_self_energy, & - compute_QP_energies, & - compute_W_cubic_GW, & - deallocate_matrices_gw, & - deallocate_matrices_gw_im_time, & - get_fermi_level_offset - USE rpa_gw_ic, ONLY: calculate_ic_correction - USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_MOs, & - invert_eps_compute_W_and_Erpa_kp - USE rpa_im_time, ONLY: compute_mat_P_omega, & - zero_mat_P_omega - USE rpa_im_time_force_methods, ONLY: calc_laplace_loop_forces, & - calc_post_loop_forces, & - calc_rpa_loop_forces, & - init_im_time_forces, & - keep_initial_quad - USE rpa_im_time_force_types, ONLY: im_time_force_release, & - im_time_force_type - USE rpa_sigma_functional, ONLY: finalize_rpa_sigma, & - rpa_sigma_create, & - rpa_sigma_matrix_spectral, & - rpa_sigma_type - USE rpa_util, ONLY: Q_trace_and_add_unit_matrix, & - alloc_im_time, & - calc_mat_Q, & - compute_Erpa_by_freq_int, & - contract_P_omega_with_mat_L, & - dealloc_im_time, & - remove_scaling_factor_rpa - USE util, ONLY: get_limit -#if defined (__GREENX) - USE gx_minimax, ONLY: gx_minimax_grid -#endif - + USE bibliography, ONLY: & + Bates2013, DelBen2013, DelBen2015, Freeman1977, Gruneis2009, Ren2011, Ren2013, & + Wilhelm2016a, Wilhelm2016b, Wilhelm2017, Wilhelm2018, cite_reference + USE bse_main, ONLY: start_bse_calculation + USE cp_blacs_env, ONLY: cp_blacs_env_create,& + cp_blacs_env_release,& + cp_blacs_env_type + USE cp_cfm_types, ONLY: cp_cfm_type + USE cp_dbcsr_api, ONLY: dbcsr_add,& + dbcsr_clear,& + dbcsr_get_info,& + dbcsr_p_type,& + dbcsr_type + USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add + 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_set_all,& + cp_fm_to_fm,& + cp_fm_type + USE dbt_api, ONLY: dbt_type + USE dgemm_counter_types, ONLY: dgemm_counter_init,& + dgemm_counter_type,& + dgemm_counter_write + USE group_dist_types, ONLY: create_group_dist,& + get_group_dist,& + group_dist_d1_type,& + maxsize,& + release_group_dist + USE hfx_types, ONLY: block_ind_type,& + hfx_compression_type + USE input_constants, ONLY: rpa_exchange_axk,& + rpa_exchange_none,& + rpa_exchange_sosex,& + sigma_none,& + wfc_mm_style_gemm + USE kinds, ONLY: dp,& + int_8 + USE kpoint_types, ONLY: kpoint_type + USE machine, ONLY: m_flush,& + m_memory + USE mathconstants, ONLY: pi,& + z_zero + USE message_passing, ONLY: mp_comm_type,& + mp_para_env_release,& + mp_para_env_type + USE minimax_exp, ONLY: check_exp_minimax_range + USE mp2_grids, ONLY: get_clenshaw_grid,& + get_minimax_grid + USE mp2_laplace, ONLY: SOS_MP2_postprocessing + USE mp2_ri_grad_util, ONLY: array2fm + USE mp2_types, ONLY: mp2_type,& + three_dim_real_array,& + two_dim_int_array,& + two_dim_real_array + USE qs_environment_types, ONLY: get_qs_env,& + qs_environment_type + USE rpa_exchange, ONLY: rpa_exchange_needed_mem,& + rpa_exchange_work_type + USE rpa_grad, ONLY: rpa_grad_copy_Q,& + rpa_grad_create,& + rpa_grad_finalize,& + rpa_grad_matrix_operations,& + rpa_grad_needed_mem,& + rpa_grad_type + USE rpa_gw, ONLY: allocate_matrices_gw,& + allocate_matrices_gw_im_time,& + compute_GW_self_energy,& + compute_QP_energies,& + compute_W_cubic_GW,& + deallocate_matrices_gw,& + deallocate_matrices_gw_im_time,& + get_fermi_level_offset + USE rpa_gw_ic, ONLY: calculate_ic_correction + USE rpa_gw_kpoints_util, ONLY: get_bandstruc_and_k_dependent_MOs,& + invert_eps_compute_W_and_Erpa_kp + USE rpa_im_time, ONLY: compute_mat_P_omega,& + zero_mat_P_omega + USE rpa_im_time_force_methods, ONLY: calc_laplace_loop_forces,& + calc_post_loop_forces,& + calc_rpa_loop_forces,& + init_im_time_forces,& + keep_initial_quad + USE rpa_im_time_force_types, ONLY: im_time_force_release,& + im_time_force_type + USE rpa_sigma_functional, ONLY: finalize_rpa_sigma,& + rpa_sigma_create,& + rpa_sigma_matrix_spectral,& + rpa_sigma_type + USE rpa_util, ONLY: Q_trace_and_add_unit_matrix,& + alloc_im_time,& + calc_mat_Q,& + compute_Erpa_by_freq_int,& + contract_P_omega_with_mat_L,& + dealloc_im_time,& + remove_scaling_factor_rpa + USE util, ONLY: get_limit #include "./base/base_uses.f90" IMPLICIT NONE @@ -1117,29 +1112,22 @@ CONTAINS COMPLEX(KIND=dp), ALLOCATABLE, & DIMENSION(:, :, :, :) :: vec_Sigma_c_gw - INTEGER :: count_ev_sc_GW, cut_memory, group_size_P, gw_corr_lev_tot, handle, handle3, & - i, ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, min_bsize, mm_style, nkp, & - nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, Pspin, Qspin, & - size_P -#if defined (__GREENX) - INTEGER :: gi, ierr -#endif + INTEGER :: count_ev_sc_GW, cut_memory, group_size_P, gw_corr_lev_tot, handle, handle3, i, & + ikp_local, ispin, iter_evGW, iter_sc_GW0, j, jquad, min_bsize, mm_style, nkp, & + nkp_self_energy, nmo, nspins, num_3c_repl, num_cells_dm, num_fit_points, Pspin, Qspin, & + size_P INTEGER(int_8) :: dbcsr_nflop INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell_3c INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cell_to_index_3c INTEGER, DIMENSION(:), POINTER :: col_blk_size, prim_blk_sizes, & RI_blk_sizes LOGICAL :: do_apply_ic_corr_to_gw, do_gw_im_time, do_ic_model, do_kpoints_cubic_RPA, & - do_periodic, do_print, do_ri_Sigma_x, exit_ev_gw, first_cycle, & - first_cycle_periodic_correction, my_open_shell, print_ic_values + do_periodic, do_print, do_ri_Sigma_x, exit_ev_gw, first_cycle, & + first_cycle_periodic_correction, my_open_shell, print_ic_values LOGICAL, ALLOCATABLE, DIMENSION(:, :, :, :, :) :: has_mat_P_blocks - REAL(KIND=dp) :: a_scaling, alpha, dbcsr_time, e_exchange, e_exchange_corr, & - eps_filter, eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, & - my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old -#if defined (__GREENX) - REAL(KIND=dp) :: cosft_duality_error_greenx, Emax, Emin, & - max_errors_greenx(3) -#endif + REAL(KIND=dp) :: a_scaling, alpha, dbcsr_time, e_exchange, e_exchange_corr, eps_filter, & + eps_filter_im_time, ext_scaling, fermi_level_offset, fermi_level_offset_input, & + my_flop_rate, omega, omega_max_fit, omega_old, tau, tau_old REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta_corr, e_fermi, tau_tj, tau_wj, tj, & trace_Qomega, vec_omega_fit_gw, wj, & wkp_W @@ -1150,8 +1138,8 @@ CONTAINS vec_Sigma_x_gw TYPE(cp_cfm_type) :: cfm_mat_Q TYPE(cp_fm_type) :: fm_mat_Q_static_bse_gemm, fm_mat_RI_global_work, fm_mat_S_ia_bse, & - fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, & - fm_scaled_dm_virt_tau + fm_mat_work, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, fm_scaled_dm_occ_tau, & + fm_scaled_dm_virt_tau TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_mat_S_gw_work, fm_mat_W, & fm_mo_coeff_occ, fm_mo_coeff_virt TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: fm_mat_L_kpoints, fm_mat_Minv_L_kpoints @@ -1288,80 +1276,12 @@ CONTAINS ALLOCATE (e_fermi(nspins)) IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN do_print = .NOT. do_ic_model -#if defined (__GREENX) - CALL init_greenx_grids(qs_env, para_env, homo, Eigenval, do_ri_sos_laplace_mp2, & - do_kpoints_cubic_RPA, Emin, Emax, e_fermi(1)) - CALL gx_minimax_grid(num_integ_points, Emin, Emax, tau_tj, tau_wj, tj, wj, & - weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, & - max_errors_greenx, cosft_duality_error_greenx, ierr, & - bare_cos_sin_weights=.TRUE., & - regularization=qs_env%mp2_env%ri_g0w0%regularization_minimax) - ! Factor 4 is hard-coded in the RPA weights in the internal CP2K minimax routines - wj(:) = wj(:)*4.0_dp - IF (ierr == 0) THEN - IF (unit_nr > 0) THEN - WRITE (UNIT=unit_nr, FMT="(T3,A,T75,i6)") & - "GREENX MINIMAX_INFO| Number of integration points:", num_integ_points - WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") & - "GREENX MINIMAX_INFO| Gap (Emin):", Emin - WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.7)") & - "GREENX MINIMAX_INFO| Maximum eigenvalue difference (Emax):", Emax - WRITE (UNIT=unit_nr, FMT="(T3,A,T61,3F20.4)") & - "GREENX MINIMAX_INFO| Energy range (Emax/Emin):", Emax/Emin - WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") & - "GREENX MINIMAX_INFO| Frequency grid (scaled):", "Weights", "Abscissas" - DO gi = 1, num_integ_points - WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") wj(gi), tj(gi) - END DO - WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") & - "GREENX MINIMAX_INFO| Time grid (scaled):", "Weights", "Abscissas" - DO gi = 1, num_integ_points - WRITE (UNIT=unit_nr, FMT="(T41,F20.10,F20.10)") tau_wj(gi), tau_tj(gi) - END DO - CALL m_flush(unit_nr) - END IF - ELSE - IF (unit_nr > 0) THEN - WRITE (UNIT=unit_nr, FMT="(T3,A,T75)") & - "GREENX MINIMAX_INFO| Grid not available, use internal CP2K grids" - CALL m_flush(unit_nr) - END IF - IF (ALLOCATED(tau_tj)) THEN - DEALLOCATE (tau_tj) - END IF - IF (ALLOCATED(tau_wj)) THEN - DEALLOCATE (tau_wj) - END IF - IF (ALLOCATED(tj)) THEN - DEALLOCATE (tj) - END IF - IF (ALLOCATED(wj)) THEN - DEALLOCATE (wj) - END IF - IF (ALLOCATED(weights_cos_tf_t_to_w)) THEN - DEALLOCATE (weights_cos_tf_t_to_w) - END IF - IF (ALLOCATED(weights_cos_tf_w_to_t)) THEN - DEALLOCATE (weights_cos_tf_w_to_t) - END IF - IF (ALLOCATED(weights_sin_tf_t_to_w)) THEN - DEALLOCATE (weights_sin_tf_t_to_w) - END IF - CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, & - do_ri_sos_laplace_mp2, do_print, & - tau_tj, tau_wj, qs_env, do_gw_im_time, & - do_kpoints_cubic_RPA, e_fermi(1), tj, wj, & - weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, & - qs_env%mp2_env%ri_g0w0%regularization_minimax) - END IF -#else CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, & do_ri_sos_laplace_mp2, do_print, & tau_tj, tau_wj, qs_env, do_gw_im_time, & do_kpoints_cubic_RPA, e_fermi(1), tj, wj, & weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, & qs_env%mp2_env%ri_g0w0%regularization_minimax) -#endif !For sos_laplace_mp2 and low-scaling RPA, potentially need to store/retrieve the initial weights IF (qs_env%mp2_env%ri_rpa_im_time%keep_quad) THEN @@ -1414,7 +1334,6 @@ CONTAINS IF (calc_forces .AND. .NOT. do_im_time) CALL rpa_grad_create(rpa_grad, fm_mat_Q(1), & fm_mat_S, homo, virtual, mp2_env, Eigenval(:, 1, :), & unit_nr, do_ri_sos_laplace_mp2) -!doubtful IF (.NOT. do_im_time .AND. .NOT. do_ri_sos_laplace_mp2) & CALL exchange_work%create(qs_env, para_env_sub, mat_munu, dimen_RI_red, & fm_mat_S, fm_mat_Q(1), fm_mat_Q_gemm(1), homo, virtual)