Shifted warning into mp2_grids, removed code that was used with mixed…

… Minimax-Clenshaw grids that are not in the code any more

src/mp2_grids.F

@@ -51,8 +51,7 @@ MODULE mp2_grids
 !> \param unit_nr ...
 !> \param homo ...
 !> \param Eigenval ...
-!> \param num_time_points ...
-!> \param num_freq_points ...
+!> \param num_integ_points ...
 !> \param do_im_time ...
 !> \param do_ri_sos_laplace_mp2 ...
 !> \param do_print ...
@@ -72,18 +71,17 @@ MODULE mp2_grids
 !> \param homo_beta ...
 !> \param dimen_ia_beta ...
 !> \param Eigenval_beta ...
-!> \param do_mixed_minimax_clenshaw ...
 ! **************************************************************************************************
-   SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points, num_freq_points, &
+   SUBROUTINE 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, a_scaling, e_fermi, tj, wj, mp2_env, weights_cos_tf_t_to_w, &
                                weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
-                               homo_beta, dimen_ia_beta, Eigenval_beta, do_mixed_minimax_clenshaw)
+                               homo_beta, dimen_ia_beta, Eigenval_beta)
 
       TYPE(cp_para_env_type), POINTER                    :: para_env
       INTEGER, INTENT(IN)                                :: unit_nr, homo
       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: Eigenval
-      INTEGER, INTENT(IN)                                :: num_time_points, num_freq_points
+      INTEGER, INTENT(IN)                                :: num_integ_points
       LOGICAL, INTENT(IN)                                :: do_im_time, do_ri_sos_laplace_mp2, &
                                                             do_print
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
@@ -101,16 +99,13 @@ SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points,
                                                             weights_sin_tf_t_to_w
       INTEGER, INTENT(IN), OPTIONAL                      :: homo_beta, dimen_ia_beta
       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: Eigenval_beta
-      LOGICAL, INTENT(IN), OPTIONAL                      :: do_mixed_minimax_clenshaw
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'get_minimax_grid', &
          routineP = moduleN//':'//routineN
       INTEGER, PARAMETER                                 :: num_points_per_magnitude = 200
 
-      INTEGER                                            :: handle, ierr, jquad, num_integ_points
-      LOGICAL                                            :: my_do_kpoints, &
-                                                            my_do_mixed_minimax_clenshaw, &
-                                                            my_open_shell
+      INTEGER                                            :: handle, ierr, jquad
+      LOGICAL                                            :: my_do_kpoints, my_open_shell
       REAL(KIND=dp)                                      :: E_Range, Emax, Emax_beta, Emin, &
                                                             Emin_beta, max_error_min, scaling
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: x_tw
@@ -134,11 +129,6 @@ SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points,
          my_do_kpoints = do_kpoints_cubic_RPA
       END IF
 
-      my_do_mixed_minimax_clenshaw = .FALSE.
-      IF (PRESENT(do_mixed_minimax_clenshaw)) THEN
-         my_do_mixed_minimax_clenshaw = do_mixed_minimax_clenshaw
-      END IF
-
       IF (my_do_kpoints) THEN
          CALL gap_and_max_eig_diff_kpoints(qs_env, para_env, Emin, Emax, e_fermi)
          E_Range = Emax/Emin
@@ -156,29 +146,34 @@ SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points,
          E_Range = Emax/Emin
       END IF
 
-      IF (my_do_mixed_minimax_clenshaw) THEN
-         tj(:) = a_scaling/TAN(tj(:))
+      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
 
-      IF (.NOT. do_ri_sos_laplace_mp2 .AND. (.NOT. my_do_mixed_minimax_clenshaw)) THEN
-         ALLOCATE (x_tw(2*num_freq_points))
+      IF (.NOT. do_ri_sos_laplace_mp2) THEN
+         ALLOCATE (x_tw(2*num_integ_points))
          x_tw = 0.0_dp
          ierr = 0
-         IF (num_freq_points .LE. 20) THEN
-            CALL get_rpa_minimax_coeff(num_freq_points, E_Range, x_tw, ierr)
+         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_freq_points, E_Range, x_tw)
+            CALL get_rpa_minimax_coeff_larger_grid(num_integ_points, E_Range, x_tw)
          END IF
 
-         ALLOCATE (tj(num_freq_points))
+         ALLOCATE (tj(num_integ_points))
          tj = 0.0_dp
 
-         ALLOCATE (wj(num_freq_points))
+         ALLOCATE (wj(num_integ_points))
          wj = 0.0_dp
 
-         DO jquad = 1, num_freq_points
+         DO jquad = 1, num_integ_points
             tj(jquad) = x_tw(jquad)
-            wj(jquad) = x_tw(jquad + num_freq_points)
+            wj(jquad) = x_tw(jquad + num_integ_points)
          END DO
 
          DEALLOCATE (x_tw)
@@ -187,7 +182,7 @@ SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points,
             WRITE (UNIT=unit_nr, FMT="(T3,A,T66,F15.4)") &
                "INTEG_INFO| Range for the minimax approximation:", E_Range
             WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") "INTEG_INFO| Minimax parameters:", "Weights", "Abscissas"
-            DO jquad = 1, num_freq_points
+            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)
@@ -206,28 +201,28 @@ SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points,
       ! set up the minimax time grid
       IF (do_im_time .OR. do_ri_sos_laplace_mp2) THEN
 
-         ALLOCATE (x_tw(2*num_time_points))
+         ALLOCATE (x_tw(2*num_integ_points))
          x_tw = 0.0_dp
 
-         IF (num_freq_points .LE. 20) THEN
-            CALL get_exp_minimax_coeff(num_time_points, E_Range, x_tw)
+         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_time_points, E_Range, x_tw)
+            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(0:num_time_points))
+         ALLOCATE (tau_tj(0:num_integ_points))
          tau_tj = 0.0_dp
 
-         ALLOCATE (tau_wj(num_time_points))
+         ALLOCATE (tau_wj(num_integ_points))
          tau_wj = 0.0_dp
 
-         DO jquad = 1, num_time_points
+         DO jquad = 1, num_integ_points
             tau_tj(jquad) = x_tw(jquad)/scaling
-            tau_wj(jquad) = x_tw(jquad + num_time_points)/scaling
+            tau_wj(jquad) = x_tw(jquad + num_integ_points)/scaling
          END DO
 
          DEALLOCATE (x_tw)
@@ -240,7 +235,7 @@ SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points,
                "INTEG_INFO| Gap:", Emin
             WRITE (UNIT=unit_nr, FMT="(T3,A,T54,A,T72,A)") &
                "INTEG_INFO| Minimax parameters of the time grid:", "Weights", "Abscissas"
-            DO jquad = 1, num_time_points
+            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)
@@ -251,16 +246,13 @@ SUBROUTINE get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_time_points,
          tau_wj(:) = tau_wj(:)/Emin
 
          IF (.NOT. do_ri_sos_laplace_mp2) THEN
-            ALLOCATE (weights_cos_tf_t_to_w(num_freq_points, num_time_points))
+            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_freq_points, num_time_points, tau_tj, weights_cos_tf_t_to_w, tj, &
+            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)
 
-            IF (do_gw_im_time .AND. (.NOT. my_do_mixed_minimax_clenshaw)) THEN
-
-               CPASSERT(num_freq_points == num_time_points)
-               num_integ_points = num_freq_points
+            IF (do_gw_im_time) THEN
 
                ! 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))
@@ -764,8 +756,7 @@ END SUBROUTINE calculate_objfunc
 
 ! **************************************************************************************************
 !> \brief Calculate integration weights for the tau grid (in dependency of the omega node)
-!> \param num_freq_points ...
-!> \param num_time_points ...
+!> \param num_integ_points ...
 !> \param tau_tj ...
 !> \param weights_cos_tf_t_to_w ...
 !> \param omega_tj ...
@@ -774,10 +765,10 @@ END SUBROUTINE calculate_objfunc
 !> \param max_error ...
 !> \param num_points_per_magnitude ...
 ! **************************************************************************************************
-   SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_freq_points, num_time_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, &
+   SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, weights_cos_tf_t_to_w, omega_tj, &
                                            E_min, E_max, max_error, num_points_per_magnitude)
 
-      INTEGER, INTENT(IN)                                :: num_freq_points, num_time_points
+      INTEGER, INTENT(IN)                                :: num_integ_points
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
          INTENT(IN)                                      :: tau_tj
       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), &
@@ -807,38 +798,38 @@ SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_freq_points, num_time_points, tau_tj
 
       ! take at least as many x points as integration points to have clear
       ! input for the singular value decomposition
-      num_x_nodes = MAX(num_x_nodes, num_freq_points)
+      num_x_nodes = MAX(num_x_nodes, num_integ_points)
 
       ALLOCATE (x_values(num_x_nodes))
       x_values = 0.0_dp
       ALLOCATE (y_values(num_x_nodes))
       y_values = 0.0_dp
-      ALLOCATE (mat_A(num_x_nodes, num_time_points))
+      ALLOCATE (mat_A(num_x_nodes, num_integ_points))
       mat_A = 0.0_dp
-      ALLOCATE (tau_wj_work(num_time_points))
+      ALLOCATE (tau_wj_work(num_integ_points))
       tau_wj_work = 0.0_dp
-      ALLOCATE (sing_values(num_time_points))
+      ALLOCATE (sing_values(num_integ_points))
       sing_values = 0.0_dp
       ALLOCATE (mat_U(num_x_nodes, num_x_nodes))
       mat_U = 0.0_dp
-      ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_time_points))
+      ALLOCATE (mat_SinvVSinvT(num_x_nodes, num_integ_points))
 
       mat_SinvVSinvT = 0.0_dp
       ! double the value nessary for 'A' to achieve good performance
-      lwork = 8*num_time_points*num_time_points + 12*num_time_points + 2*num_x_nodes
+      lwork = 8*num_integ_points*num_integ_points + 12*num_integ_points + 2*num_x_nodes
       ALLOCATE (work(lwork))
       work = 0.0_dp
-      ALLOCATE (iwork(8*num_time_points))
+      ALLOCATE (iwork(8*num_integ_points))
       iwork = 0
-      ALLOCATE (mat_SinvVSinvSigma(num_time_points, num_x_nodes))
+      ALLOCATE (mat_SinvVSinvSigma(num_integ_points, num_x_nodes))
       mat_SinvVSinvSigma = 0.0_dp
       ALLOCATE (vec_UTy(num_x_nodes))
       vec_UTy = 0.0_dp
 
       max_error = 0.0_dp
 
       ! loop over all omega frequency points
-      DO jquad = 1, num_freq_points
+      DO jquad = 1, num_integ_points
 
          ! set the x-values logarithmically in the interval [Emin,Emax]
          multiplicator = (E_max/E_min)**(1.0_dp/(REAL(num_x_nodes, KIND=dp) - 1.0_dp))
@@ -854,22 +845,22 @@ SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_freq_points, num_time_points, tau_tj
          END DO
 
          ! calculate mat_A
-         DO jjj = 1, num_time_points
+         DO jjj = 1, num_integ_points
             DO iii = 1, num_x_nodes
                mat_A(iii, jjj) = COS(omega*tau_tj(jjj))*EXP(-x_values(iii)*tau_tj(jjj))
             END DO
          END DO
 
          ! Singular value decomposition of mat_A
-         CALL DGESDD('A', num_x_nodes, num_time_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
+         CALL DGESDD('A', num_x_nodes, num_integ_points, mat_A, num_x_nodes, sing_values, mat_U, num_x_nodes, &
                      mat_SinvVSinvT, num_x_nodes, work, lwork, iwork, info)
 
          CPASSERT(info == 0)
 
          ! integration weights = V Sigma U^T y
          ! 1) V*Sigma
-         DO jjj = 1, num_time_points
-            DO iii = 1, num_time_points
+         DO jjj = 1, num_integ_points
+            DO iii = 1, num_integ_points
                mat_SinvVSinvSigma(iii, jjj) = mat_SinvVSinvT(jjj, iii)/sing_values(jjj)
             END DO
          END DO
@@ -879,13 +870,13 @@ SUBROUTINE get_l_sq_wghts_cos_tf_t_to_w(num_freq_points, num_time_points, tau_tj
                     0.0_dp, vec_UTy, num_x_nodes)
 
          ! 3) (V*Sigma) * (U^T y)
-         CALL DGEMM('N', 'N', num_time_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_time_points, vec_UTy, &
-                    num_x_nodes, 0.0_dp, tau_wj_work, num_time_points)
+         CALL DGEMM('N', 'N', num_integ_points, 1, num_x_nodes, 1.0_dp, mat_SinvVSinvSigma, num_integ_points, vec_UTy, &
+                    num_x_nodes, 0.0_dp, tau_wj_work, num_integ_points)
 
          weights_cos_tf_t_to_w(jquad, :) = tau_wj_work(:)
 
          CALL calc_max_error_fit_tau_grid_with_cosine(max_error, omega, tau_tj, tau_wj_work, x_values, &
-                                                      y_values, num_time_points, num_x_nodes)
+                                                      y_values, num_integ_points, num_x_nodes)
 
       END DO ! jquad
 
@@ -1252,7 +1243,7 @@ SUBROUTINE test_least_square_ft(nR, iw)
                tau_wj(jquad) = x_tw(jquad + num_integ_points)/2.0_dp
             END DO
 
-            CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, num_integ_points, tau_tj, &
+            CALL get_l_sq_wghts_cos_tf_t_to_w(num_integ_points, tau_tj, &
                                               weights_cos_tf_t_to_w, tj, 1.0_dp, Rc, max_error, 200)
 
             IF (iw > 0) THEN

src/rpa_main.F

@@ -313,13 +313,6 @@ SUBROUTINE rpa_ri_compute_en(qs_env, Erpa, mp2_env, BIb_C, BIb_C_gw, BIb_C_bse_i
                                "Minimax quadrature scheme. The number of quadrature point has been reset to 30.")
                num_integ_points = 30
             END IF
-            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 for a small minimax range (< 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
          ELSE
             IF (do_minimax_quad .AND. num_integ_points > 20) THEN
                IF (unit_nr > 0) &
@@ -1512,13 +1505,13 @@ SUBROUTINE rpa_num_int(qs_env, Erpa, mp2_env, para_env, para_env_RPA, para_env_s
 
       IF (do_minimax_quad .OR. do_ri_sos_laplace_mp2) THEN
          IF (my_open_shell) THEN
-            CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, num_integ_points, do_im_time, &
+            CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, &
                                   do_ri_sos_laplace_mp2,.NOT. do_ic_model, tau_tj, tau_wj, qs_env, do_gw_im_time, &
                                   do_kpoints_cubic_RPA, a_scaling, e_fermi, tj, wj, mp2_env, &
                                   weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w, &
                                   homo_beta, dimen_ia_beta, Eigenval_beta)
          ELSE
-            CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, num_integ_points, do_im_time, &
+            CALL get_minimax_grid(para_env, unit_nr, homo, Eigenval, num_integ_points, do_im_time, &
                                   do_ri_sos_laplace_mp2,.NOT. do_ic_model, tau_tj, tau_wj, qs_env, do_gw_im_time, &
                                   do_kpoints_cubic_RPA, a_scaling, e_fermi, tj, wj, mp2_env, &
                                   weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, weights_sin_tf_t_to_w)