Print Sigma(i*omega) and corresponding fit in case &PRINT SELF_ENERGY…

… is activated
cp2k · Nov 15, 2020 · 76181cb · 76181cb
1 parent 97eec91
commit 76181cb
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 4 deletions.
diff --git a/src/input_cp2k_mp2.F b/src/input_cp2k_mp2.F
@@ -836,7 +836,8 @@ SUBROUTINE create_print_section(section)
       CALL keyword_create(keyword, __LOCATION__, name="", &
                           variants=(/"SELF_ENERGY"/), &
                           description="If true, print the self-energy for all levels for real energy "// &
-                          "together with the straight line to see the quasiparticle energy as intersection.", &
+                          "together with the straight line to see the quasiparticle energy as intersection. "// &
+                          "In addition, prints the self-energy for imaginary frequencies together with the Pade fit.", &
                           usage="SELF_ENERGY TRUE", &
                           default_l_val=.FALSE., &
                           lone_keyword_l_val=.TRUE.)

diff --git a/src/rpa_gw.F b/src/rpa_gw.F
@@ -3432,7 +3432,8 @@ SUBROUTINE continuation_pade(vec_gw_energ, vec_omega_fit_gw, &
       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'continuation_pade'
 
       CHARACTER(len=default_path_length)                 :: filename
-      COMPLEX(KIND=dp)                                   :: im_unit, re_unit, sigma_c_pade
+      COMPLEX(KIND=dp)                                   :: im_unit, re_unit, sigma_c_pade, &
+                                                            sigma_c_pade_im_freq
       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:)        :: coeff_pade, omega_points_pade, &
                                                             Sigma_c_gw_reorder
       INTEGER                                            :: handle, i_omega, iunit, jquad, &
@@ -3595,6 +3596,23 @@ SUBROUTINE continuation_pade(vec_gw_energ, vec_omega_fit_gw, &
 
          END DO
 
+         WRITE (iunit, "(A51,A39)") " w (eV)  Re(Sigma(i*w)) (eV)   Im(Sigma(i*w)) (eV) ", &
+            "  Re(Fit(i*w)) (eV)    Im(Fit(iw)) (eV)"
+
+         DO jquad = 1, num_fit_points
+
+            CALL evaluate_pade_function(vec_omega_fit_gw_sign_reorder(jquad), &
+                                        nparam_pade, omega_points_pade, &
+                                        coeff_pade, sigma_c_pade_im_freq, do_imag_freq=.TRUE.)
+
+            WRITE (iunit, "(F12.2,4F17.5)") vec_omega_fit_gw_sign_reorder(jquad)*evolt, &
+               REAL(Sigma_c_gw_reorder(jquad)*evolt), &
+               AIMAG(Sigma_c_gw_reorder(jquad)*evolt), &
+               REAL(sigma_c_pade_im_freq*evolt), &
+               AIMAG(sigma_c_pade_im_freq*evolt)
+
+         END DO
+
          CALL close_file(iunit)
 
       END IF
@@ -3680,27 +3698,37 @@ END SUBROUTINE get_pade_parameters
 !> \param xpoints selection of points of the original complex function, i.e. here of Sigma_c(iomega)
 !> \param coeff pade coefficients
 !> \param func_val function value
+!> \param do_imag_freq ...
 ! **************************************************************************************************
-   SUBROUTINE evaluate_pade_function(x_val, nparam, xpoints, coeff, func_val)
+   SUBROUTINE evaluate_pade_function(x_val, nparam, xpoints, coeff, func_val, do_imag_freq)
 
       REAL(KIND=dp), INTENT(IN)                          :: x_val
       INTEGER, INTENT(IN)                                :: nparam
       COMPLEX(KIND=dp), DIMENSION(:), INTENT(IN)         :: xpoints, coeff
       COMPLEX(KIND=dp), INTENT(OUT)                      :: func_val
+      LOGICAL, INTENT(IN), OPTIONAL                      :: do_imag_freq
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_pade_function'
 
       COMPLEX(KIND=dp)                                   :: im_unit, re_unit
       INTEGER                                            :: handle, iparam
+      LOGICAL                                            :: my_do_imag_freq
 
       CALL timeset(routineN, handle)
 
+      my_do_imag_freq = .FALSE.
+      IF (PRESENT(do_imag_freq)) my_do_imag_freq = do_imag_freq
+
       im_unit = (0.0_dp, 1.0_dp)
       re_unit = (1.0_dp, 0.0_dp)
 
       func_val = re_unit
       DO iparam = nparam, 2, -1
-         func_val = re_unit + coeff(iparam)*(re_unit*x_val - xpoints(iparam - 1))/func_val
+         IF (my_do_imag_freq) THEN
+            func_val = re_unit + coeff(iparam)*(im_unit*x_val - xpoints(iparam - 1))/func_val
+         ELSE
+            func_val = re_unit + coeff(iparam)*(re_unit*x_val - xpoints(iparam - 1))/func_val
+         END IF
       ENDDO
 
       func_val = coeff(1)/func_val