comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-28 14:15:19 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Move calls to GreenX to a separate module (#4307)

Browse source
This commit is contained in:
Frederick Stein 2025-07-18 09:53:51 +02:00 • committed by GitHub
parent 858b7a110b
commit 67ae0a34c1
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
6 changed files with 662 additions and 542 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/CMakeLists.txt
View file

@ -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

182
src/emd/rt_bse_io.F
View file

@ -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)

8
src/emd/rt_bse_types.F
View file

@ -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)

383
src/greenx_interface.F Normal file
View file

@ -0,0 +1,383 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! 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

329
src/mp2_grids.F
View file

@ -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)

301
src/rpa_main.F
View file

@ -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)