TABLE OF CONTENTS


ABINIT/m_gwls_ComputeCorrelationEnergy [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_ComputeCorrelationEnergy

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2022 ABINIT group (JLJ, BR, MC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_gwls_ComputeCorrelationEnergy
23 
24 ! local modules
25 use m_gwls_utility
26 use m_gwls_wf
27 use m_gwls_hamiltonian
28 use m_gwls_lineqsolver
29 use m_gwls_polarisability
30 use m_gwls_model_polarisability
31 use m_gwls_DielectricArray
32 use m_gwls_ComputePoles
33 use m_gwls_Projected_AT
34 use m_gwls_Projected_BT
35 use m_gwls_GWlanczos
36 use m_gwls_GenerateEpsilon
37 use m_gwls_GWanalyticPart
38 use m_gwls_TimingLog
39 use m_gwls_LanczosBasis
40 
41 ! abinit modules
42 use defs_basis
43 use defs_wvltypes
44 use m_abicore
45 use m_xmpi
46 use m_pawang
47 use m_errors
48 use m_dtset
49 
50 use m_time,             only : timab
51 use m_io_tools,         only : get_unit,open_file
52 use m_paw_dmft,         only : paw_dmft_type
53 use m_ebands,           only : ebands_init, ebands_free
54 use m_gaussian_quadrature, only: gaussian_quadrature_gegenbauer, gaussian_quadrature_legendre
55 
56 
57 implicit none
58 save
59 private

m_gwls_ComputeCorrelationEnergy/compute_correlations_no_model_shift_lanczos [ Functions ]

[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  compute_correlations_no_model_shift_lanczos

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

 901 subroutine compute_correlations_no_model_shift_lanczos(dtset, Sigma_x,Vxc_energy,debug)
 902 !----------------------------------------------------------------------------------------------------
 903 !
 904 !
 905 ! This function computes the correlation energy Sigma_c(w) for the state we wish to correct.
 906 !
 907 ! this subroutine does not rely on the use of a model dielectric operator. Thus
 908 !
 909 !        Sigma^A(w) = Tr[ (eps^{-1}(0)-1). A^T(w)]
 910 !        Sigma^N(w) = int dw' Tr[{(eps^{-1}(w')-1)-f(w')(eps^{-1}(0)-1)}B^T(w';w)]
 911 !
 912 ! Shift lanczos is used for the resolvents.
 913 !----------------------------------------------------------------------------------------------------
 914 
 915 type(dataset_type),intent(in) :: dtset
 916 
 917 real(dp),intent(in)  :: Sigma_x, Vxc_energy
 918 logical, intent(in)  :: debug
 919 
 920 
 921 !Local variables
 922 
 923 real(dp):: Sigma_x_Lanczos_projected
 924 integer :: npt_gauss
 925 integer :: print_debug
 926 
 927 
 928 integer :: kmax_poles
 929 
 930 integer :: lmax_model
 931 
 932 integer :: kmax_analytic
 933 integer :: kmax_numeric
 934 integer :: n_ext_freq
 935 
 936 
 937 real(dp) :: omega_static
 938 
 939 real(dp) :: lorentzian
 940 integer  :: iw_ext
 941 integer  :: iw
 942 
 943 
 944 
 945 real(dp) , allocatable :: psie_k(:,:)
 946 
 947 real(dp),  allocatable :: epsilon_eigenvalues_0(:)
 948 
 949 
 950 complex(dpc), allocatable   :: AT_Lanczos(:,:)
 951 
 952 
 953 
 954 real(dp)       :: time1, time2, time
 955 real(dp)       :: total_time1,total_time2,total_time
 956 real(dp)       :: setup_time1, setup_time2, setup_time
 957 real(dp)       :: freq_time1, freq_time2, freq_time
 958 
 959 integer              :: nfrequencies
 960 real(dp),allocatable :: list_projection_frequencies(:)
 961 
 962 
 963 
 964 complex(dpc), allocatable   :: array_integrand_exact_sector(:,:)
 965 complex(dpc), allocatable   :: array_integrand_model_sector(:,:)
 966 complex(dpc), allocatable   :: tmp_dielectric_array(:,:,:)
 967 
 968 real(dp)        :: external_omega
 969 
 970 character(256)  :: timing_string
 971 
 972 integer :: recy_line_size
 973 character(128) :: recy_name
 974 logical :: local_tmp_exist
 975 character(500) :: msg
 976 
 977 ! Energy contributions
 978 
 979 real(dp)       :: pole_energy
 980 
 981 real(dp)       :: sigma_A_Lanczos
 982 real(dp)       :: sigma_A_model_Lanczos
 983 real(dp)       :: sigma_B_Lanczos
 984 real(dp)       :: sigma_B_model_Lanczos
 985 
 986 real(dp)       :: correlations
 987 real(dp)       :: renormalized_energy
 988 
 989 logical :: use_model
 990 
 991 ! *************************************************************************
 992 
 993 
 994 !--------------------------------------------------------------------------------
 995 !
 996 ! Set up variables and allocate arrays
 997 !
 998 !--------------------------------------------------------------------------------
 999 
1000 !Variable allocation and initialization
1001 model_number        = dtset%gwls_diel_model
1002 model_parameter     = dtset%gwls_model_parameter
1003 npt_gauss           = dtset%gwls_npt_gauss_quad
1004 print_debug         = dtset%gwls_print_debug
1005 
1006 first_seed          = dtset%gwls_first_seed
1007 e                   = dtset%gwls_band_index
1008 
1009 
1010 ! set variables from gwls_GenerateEpsilon module
1011 kmax   = dtset%gwls_stern_kmax
1012 nseeds = dtset%gwls_nseeds
1013 
1014 kmax_poles  = dtset%gwls_kmax_poles
1015 
1016 kmax_poles    = dtset%gwls_kmax_poles
1017 kmax_analytic = dtset%gwls_kmax_analytic
1018 kmax_numeric  = dtset%gwls_kmax_numeric
1019 
1020 n_ext_freq    = dtset%gw_customnfreqsp
1021 
1022 
1023 use_model     = .False.
1024 
1025 
1026 call cpu_time(setup_time1)
1027 !--------------------------------------------------------------------------------
1028 !
1029 ! Extract the frequencies at which the integrand will be evaluated
1030 ! add the value zero in the set.
1031 !
1032 !--------------------------------------------------------------------------------
1033 
1034 call generate_frequencies_and_weights(npt_gauss)
1035 
1036 !--------------------------------------------------------------------------------
1037 !
1038 !
1039 ! Compute the static bases for the exact dielectric operator
1040 !
1041 !
1042 !--------------------------------------------------------------------------------
1043 
1044 
1045 omega_static = zero
1046 ! define dimensions
1047 lmax         = nseeds*kmax
1048 
1049 ! Allocate arrays which will contain basis
1050 ! the 0 indicates we will not use arrays for the model dielectric operator
1051 call setup_Lanczos_basis(lmax,0)
1052 
1053 ! allocate eigenvalues array
1054 ABI_MALLOC(epsilon_eigenvalues_0, (lmax))
1055 
1056 ! set omega=0 for exact dielectric operator
1057 call set_dielectric_function_frequency([0.0_dp,omega_static])
1058 
1059 ! and make note that the Sternheimer solutions must be kept (for use in the projected Sternheimer section).
1060 if(dtset%gwls_recycle == 1) then
1061   ABI_MALLOC(Sternheimer_solutions_zero,(2,npw_k,lmax,nbandv))
1062   Sternheimer_solutions_zero = zero
1063   write_solution = .true.
1064 end if
1065 if(dtset%gwls_recycle == 2) then
1066   write(recy_name,'(A,I0.4,A)') "Recycling_",mpi_enreg%me,".dat"
1067 
1068   inquire(iolength=recy_line_size) cg(:,1:npw_k)
1069 
1070   inquire(file='local_tmp', exist=local_tmp_exist)
1071   if(local_tmp_exist) recy_name = 'local_tmp/' // recy_name(1:118)
1072 
1073   if (open_file(file=recy_name,iomsg=msg,newunit=recy_unit,access='direct',form='unformatted',&
1074 &               status='replace',recl=recy_line_size)/=0) then
1075     ABI_ERROR(msg)
1076   end if
1077 
1078   write_solution = .true.
1079 end if
1080 
1081 
1082 call cpu_time(time1)
1083 ! Compute the Lanczos basis using Block Lanczos; store
1084 ! basis in Lbasis_lanczos
1085 call driver_generate_dielectric_matrix( matrix_function_epsilon_k, &
1086 nseeds, kmax, &
1087 epsilon_eigenvalues_0, &
1088 Lbasis_lanczos, debug)
1089 call cpu_time(time2)
1090 time = time2-time1
1091 
1092 write(timing_string,'(A)')  "Time to compute the EXACT Static Dielectric Matrix  :   "
1093 call write_timing_log(timing_string,time)
1094 
1095 ! The Sternheimer solutions at $\omega = 0$ have been stored.
1096 write_solution = .false.
1097 
1098 
1099 call output_epsilon_eigenvalues(lmax,epsilon_eigenvalues_0,1)
1100 
1101 
1102 ! compute the Exchange energy, when it is projected on the Lanczos basis
1103 ! CAREFUL! This must be done BEFORE we modify the lanczos basis
1104 Sigma_x_Lanczos_projected =  exchange(e, Lbasis_lanczos)
1105 
1106 
1107 
1108 !--------------------------------------------------------------------------------
1109 !
1110 !
1111 ! Prepare and compute the projection of the dielectric Sternheimer equations
1112 !
1113 !
1114 !--------------------------------------------------------------------------------
1115 
1116 !  Setup various arrays necessary for the Sternheimer projection scheme
1117 nfrequencies = dtset%gwls_n_proj_freq
1118 ABI_MALLOC(list_projection_frequencies,(nfrequencies))
1119 
1120 list_projection_frequencies = dtset%gwls_list_proj_freq
1121 
1122 call cpu_time(time1)
1123 ! The explicit "false" as the last argument is for the optional
1124 ! variable "use_model"; we are not using a model here!
1125 
1126 !call setup_projected_Sternheimer_epsilon(lmax, npt_gauss, zero, &
1127 !                list_projection_frequencies,nfrequencies,debug,.false.)
1128 
1129 
1130 call ProjectedSternheimerEpsilon(lmax, npt_gauss, zero, &
1131 list_projection_frequencies,nfrequencies,&
1132 epsilon_eigenvalues_0,debug,use_model)
1133 
1134 
1135 
1136 !call cpu_time(time2)
1137 !time = time2-time1
1138 !write(timing_string,'(A)')  "Time to setup the projected Sternheimer epsilon     :   "
1139 !call write_timing_log(timing_string,time)
1140 
1141 
1142 ! The Sternheimer solutions at $\omega = 0$ have been used to make the basis for the projected Sternheimer equations.
1143 if(dtset%gwls_recycle == 1) then
1144   ABI_FREE(Sternheimer_solutions_zero)
1145 end if
1146 if(dtset%gwls_recycle == 2) then
1147   close(recy_unit,status='delete')
1148 end if
1149 
1150 
1151 !call cpu_time(time1)
1152 !call compute_projected_Sternheimer_epsilon(lmax, npt_gauss, epsilon_eigenvalues_0,debug)
1153 
1154 
1155 
1156 call cpu_time(time2)
1157 
1158 time = time2-time1
1159 write(timing_string,'(A)')  "Time to compute the projected Sternheimer epsilon   :   "
1160 call write_timing_log(timing_string,time)
1161 
1162 
1163 
1164 call cpu_time(time1)
1165 call compute_eps_m1_minus_one(lmax, npt_gauss)
1166 call cpu_time(time2)
1167 time = time2-time1
1168 write(timing_string,'(A)')  "Time to compute eps^{-1}-I                          :   "
1169 call write_timing_log(timing_string,time)
1170 
1171 
1172 !--------------------------------------------------------------------------------
1173 !
1174 !
1175 ! We no longer need the Lanczos basis in its current form!
1176 ! Modify the basis so that it now contains (V^{1/2}.l)
1177 !
1178 !
1179 !--------------------------------------------------------------------------------
1180 
1181 call cpu_time(time1)
1182 ABI_MALLOC(psie_k, (2,npw_k))
1183 psie_k = cg(:,(e-1)*npw_k+1:e*npw_k)
1184 
1185 lmax_model = 0
1186 call modify_Lbasis_Coulomb(psie_k, lmax, lmax_model) ! lmax_model is set to zero, such that
1187 ! the model lanczos basis (which doesn't exist
1188 ! in this case) will not be modified
1189 
1190 ABI_FREE(psie_k)
1191 
1192 call cpu_time(time2)
1193 time = time2-time1
1194 write(timing_string,'(A)')  "Time to modify the Lanczos basis                    :   "
1195 call write_timing_log(timing_string,time)
1196 
1197 
1198 call cpu_time(setup_time2)
1199 setup_time = setup_time2 - setup_time1
1200 
1201 write(timing_string,'(A)')  "       TOTAL DIELECTRIC SETUP TIME                  :   "
1202 call write_timing_log(timing_string,setup_time)
1203 
1204 !--------------------------------------------------------------------------------
1205 !
1206 ! Compute the Analytic energy using Shift Lanczos
1207 !
1208 !--------------------------------------------------------------------------------
1209 
1210 call cpu_time(time1)
1211 
1212 ABI_MALLOC(AT_Lanczos,(n_ext_freq,lmax))
1213 ! Note that the array eps^{-1}(0) - 1 is diagonal in the lanczos basis already! No need to diagonalize, so we can
1214 ! use Lbasis_lanczos directly...
1215 call compute_AT_shift_Lanczos(n_ext_freq,dtset%gw_freqsp, model_parameter, lmax, Lbasis_lanczos, kmax_analytic, AT_Lanczos)
1216 
1217 call cpu_time(time2)
1218 time = time2 - time1
1219 
1220 write(timing_string,'(A)')  "Time to compute analytical term by SHIFT LANCZOS    :   "
1221 call write_timing_log(timing_string,time)
1222 
1223 
1224 !--------------------------------------------------------------------------------
1225 !
1226 ! Compute the Numeric energy using Shift Lanczos
1227 !
1228 !--------------------------------------------------------------------------------
1229 
1230 call cpu_time(time1)
1231 
1232 ABI_MALLOC(array_integrand_exact_sector,(npt_gauss+1,n_ext_freq))
1233 ABI_MALLOC(array_integrand_model_sector,(npt_gauss+1,n_ext_freq))
1234 
1235 
1236 ABI_MALLOC( tmp_dielectric_array, (lmax,lmax,npt_gauss+1))
1237 
1238 do iw = 1, npt_gauss + 1
1239 
1240 lorentzian      = model_parameter**2/(list_omega(iw)**2+model_parameter**2)
1241 tmp_dielectric_array(:,:,iw) = eps_m1_minus_eps_model_m1(:,:,iw)-lorentzian*eps_m1_minus_eps_model_m1(:,:,1)
1242 
1243 end do
1244 
1245 call compute_projected_BT_shift_Lanczos(n_ext_freq, dtset%gw_freqsp, lmax, Lbasis_lanczos,         &
1246 kmax_numeric, npt_gauss, tmp_dielectric_array, array_integrand_exact_sector )
1247 
1248 array_integrand_model_sector = zero ! just a dummy array in this case
1249 
1250 ABI_FREE( tmp_dielectric_array)
1251 call cpu_time(time2)
1252 time = time2 - time1
1253 write(timing_string,'(A)')  "Time to compute numerical term by SHIFT LANCZOS     :   "
1254 call write_timing_log(timing_string,time)
1255 
1256 
1257 
1258 !--------------------------------------------------------------------------------
1259 !
1260 ! set up arrays for poles
1261 !
1262 !--------------------------------------------------------------------------------
1263 
1264 !call generate_degeneracy_table_for_poles(debug) ! so we can compute Poles contributions
1265 call generate_degeneracy_table_for_poles(.true.) ! so we can compute Poles contributions
1266 
1267 !--------------------------------------------------------------------------------
1268 !
1269 ! Print contributions to sigma_A_Lanczos  to a file
1270 !
1271 !--------------------------------------------------------------------------------
1272 call output_Sigma_A_by_eigenvalues(n_ext_freq,lmax,dtset%gw_freqsp,AT_Lanczos,one/epsilon_eigenvalues_0-one,1)
1273 
1274 
1275 !--------------------------------------------------------------------------------
1276 !
1277 ! Iterate on external frequencies
1278 !
1279 !--------------------------------------------------------------------------------
1280 
1281 
1282 do iw_ext = 1, dtset%gw_customnfreqsp
1283 
1284 call cpu_time(freq_time1)
1285 external_omega = dtset%gw_freqsp(iw_ext)
1286 
1287 write(timing_string,'(A)')  "#"
1288 call write_text_block_in_Timing_log(timing_string)
1289 write(timing_string,'(A)')  "#"
1290 call write_text_block_in_Timing_log(timing_string)
1291 write(timing_string,'(A,I4,A,F8.4,A)')  "#  Frequency # ",iw_ext," omega = ",external_omega," Ha"
1292 call write_text_block_in_Timing_log(timing_string)
1293 write(timing_string,'(A)')  "#"
1294 call write_text_block_in_Timing_log(timing_string)
1295 write(timing_string,'(A)')  "#"
1296 call write_text_block_in_Timing_log(timing_string)
1297 
1298 
1299 !--------------------------------------------------------------------------------
1300 !
1301 ! compute the pole term
1302 !                        CAREFUL! The real valence states must still be allocated
1303 !                        for the dielectric operator to work properly
1304 !
1305 !--------------------------------------------------------------------------------
1306 
1307 call cpu_time(time1)
1308 
1309 pole_energy    = compute_Poles(external_omega,kmax_poles,debug)
1310 
1311 call cpu_time(time2)
1312 
1313 time = time2-time1
1314 write(timing_string,'(A)')  "Time to compute the Poles contribution              :   "
1315 call write_timing_log(timing_string,time)
1316 
1317 !================================================================================
1318 ! Compute the contributions from the analytic term
1319 !================================================================================
1320 
1321 call cpu_time(time1)
1322 
1323 
1324 sigma_A_Lanczos      = dble(sum(AT_Lanczos(iw_ext,:)*(one/epsilon_eigenvalues_0(:)-one)))
1325 
1326 ABI_FREE(epsilon_eigenvalues_0)
1327 
1328 call cpu_time(time2)
1329 time = time2-time1
1330 write(timing_string,'(A)')  "Time for Tr[ (eps_model^{-1}-1) . AT ] AFTER SHIFT  :   "
1331 call write_timing_log(timing_string,time)
1332 
1333 !--------------------------------------------------------------------------------
1334 !
1335 ! compute integrand
1336 !
1337 !--------------------------------------------------------------------------------
1338 
1339 
1340 call compute_integrands_shift_lanczos(iw_ext, n_ext_freq, npt_gauss, array_integrand_exact_sector, &
1341 array_integrand_model_sector, sigma_B_Lanczos, sigma_B_model_Lanczos)
1342 
1343 
1344 
1345 !--------------------------------------------------------------------------------
1346 !
1347 ! Output results
1348 !
1349 !--------------------------------------------------------------------------------
1350 
1351 ! just dummy variables so we can use the output_results routine
1352 sigma_A_model_Lanczos = zero
1353 sigma_B_model_Lanczos = zero
1354 call output_results(iw_ext,npt_gauss, lmax, lmax_model, model_parameter, zero,  &
1355 external_omega, Sigma_x,Vxc_energy,pole_energy,        &
1356 sigma_A_Lanczos,sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos, &
1357 Sigma_x_Lanczos_projected )
1358 
1359 
1360 call cpu_time(freq_time2)
1361 freq_time = freq_time2-freq_time1
1362 
1363 write(timing_string,'(A)')  "               TOTAL FREQUENCY TIME                 :   "
1364 call write_timing_log(timing_string,freq_time)
1365 
1366 
1367 
1368 correlations        = pole_energy+sigma_A_Lanczos+sigma_B_Lanczos
1369 
1370 renormalized_energy = eig(e) + Sigma_x-Vxc_energy +correlations
1371 write(std_out,10) '                               '
1372 write(std_out,14) ' For omega                   : ',external_omega     ,' Ha = ',external_omega     *Ha_eV,' eV'
1373 write(std_out,14) '  <psi_e | Sigma_c  | psi_e> : ',correlations       ,' Ha = ',correlations       *Ha_eV,' eV'
1374 write(std_out,14) '  eps_e + <Sigma_xc - V_xc>  : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV'
1375 
1376 write(ab_out,10) '                               '
1377 write(ab_out,14) ' For omega                   : ',external_omega     ,' Ha = ',external_omega     *Ha_eV,' eV'
1378 write(ab_out,14) '   <psi_e | Sigma_c  | psi_e>: ',correlations       ,' Ha = ',correlations       *Ha_eV,' eV'
1379 write(ab_out,14) '   eps_e + <Sigma_xc - V_xc> : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV'
1380 
1381 
1382 
1383 end do
1384 
1385 ABI_FREE(AT_Lanczos)
1386 ABI_FREE(array_integrand_exact_sector)
1387 ABI_FREE(array_integrand_model_sector)
1388 call clean_degeneracy_table_for_poles()
1389 call cleanup_Pk_model()
1390 call cleanup_Lanczos_basis()
1391 call cleanup_projected_Sternheimer_epsilon()
1392 
1393 call cpu_time(total_time2)
1394 total_time1 = zero
1395 total_time = total_time2-total_time1
1396 write(timing_string,'(A)')  "               TOTAL TIME                           :   "
1397 call write_timing_log(timing_string,total_time)
1398 
1399 
1400 
1401 10 format(A)
1402 14 format(A,ES24.16,A,F16.8,A)
1403 
1404 end subroutine compute_correlations_no_model_shift_lanczos

m_gwls_ComputeCorrelationEnergy/compute_correlations_shift_lanczos [ Functions ]

[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  compute_correlations_shift_lanczos

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

 81 subroutine compute_correlations_shift_lanczos(dtset, Sigma_x,Vxc_energy,debug)
 82 !----------------------------------------------------------------------------------------------------
 83 !
 84 !
 85 ! This function computes the correlation energy Sigma_c(w) for the state we wish to correct.
 86 ! The shift lanczos algorithm is used.
 87 !
 88 !----------------------------------------------------------------------------------------------------
 89 
 90 type(dataset_type),intent(in) :: dtset
 91 
 92 real(dp),intent(in)  :: Sigma_x,Vxc_energy
 93 logical, intent(in)  :: debug
 94 
 95 
 96 !Local variables
 97 
 98 integer :: n_ext_freq
 99 
100 integer :: npt_gauss
101 integer :: print_debug
102 
103 
104 integer :: kmax_poles
105 integer :: kmax_model
106 integer :: lmax_model
107 
108 
109 integer :: kmax_analytic
110 integer :: kmax_numeric
111 
112 real(dp) :: omega_static
113 
114 real(dp) :: lorentzian
115 integer  :: iw_ext
116 integer  :: iw
117 
118 real(dp) , allocatable :: psie_k(:,:)
119 
120 
121 real(dp),  allocatable :: epsilon_eigenvalues_0(:)
122 real(dp),  allocatable :: epsilon_model_eigenvalues_0(:)
123 
124 
125 complex(dpc), allocatable   :: AT_Lanczos(:,:)
126 complex(dpc), allocatable   :: AT_model_Lanczos(:,:)
127 
128 
129 ! To use lanczos instead of sqmr
130 complex(dpc), allocatable   :: Lbasis_diagonalize_dielectric_terms(:,:)
131 complex(dpc), allocatable   :: hermitian_static_eps_m1_minus_eps_model_m1(:,:)
132 real(dp), allocatable       :: eigenvalues_static_eps_m1_minus_eps_model_m1(:)
133 real(dp), allocatable       :: eigenvalues_static_eps_model_m1_minus_one(:)
134 complex(dpc), allocatable   :: work(:)
135 real(dp), allocatable       :: rwork(:)
136 integer                     :: lwork
137 integer                     :: info
138 integer                     :: l
139 
140 
141 integer        :: debug_unit
142 character(50)  :: debug_filename
143 
144 real(dp)       :: time1, time2, time
145 real(dp)       :: total_time1,total_time2,total_time
146 real(dp)       :: setup_time1, setup_time2, setup_time
147 real(dp)       :: freq_time1, freq_time2, freq_time
148 
149 integer              :: nfrequencies
150 real(dp),allocatable :: list_projection_frequencies(:)
151 
152 complex(dpc), allocatable   :: array_integrand_exact_sector(:,:)
153 complex(dpc), allocatable   :: array_integrand_model_sector(:,:)
154 
155 
156 complex(dpc), allocatable :: tmp_dielectric_array(:,:,:)
157 
158 real(dp)        :: external_omega
159 
160 character(256)  :: timing_string
161 
162 integer :: recy_line_size
163 character(128) :: recy_name
164 logical :: local_tmp_exist
165 logical :: use_model
166 
167 ! Energy contributions
168 
169 real(dp)       :: pole_energy
170 
171 real(dp)       :: sigma_A_Lanczos
172 real(dp)       :: sigma_A_model_Lanczos
173 real(dp)       :: sigma_B_Lanczos
174 real(dp)       :: sigma_B_model_Lanczos
175 
176 real(dp)       :: correlations
177 real(dp)       :: renormalized_energy
178 
179 real(dp)       :: second_model_parameter
180 
181 real(dp):: tsec(2)
182 integer :: GWLS_TIMAB, OPTION_TIMAB
183 character(500) :: msg
184 
185 ! *************************************************************************
186 
187 !--------------------------------------------------------------------------------
188 !
189 ! Set up variables and allocate arrays
190 !
191 !--------------------------------------------------------------------------------
192 
193 call cpu_time(total_time1)
194 !Variable allocation and initialization
195 model_number        = dtset%gwls_diel_model
196 model_parameter     = dtset%gwls_model_parameter
197 npt_gauss           = dtset%gwls_npt_gauss_quad
198 print_debug         = dtset%gwls_print_debug
199 
200 first_seed          = dtset%gwls_first_seed
201 e                   = dtset%gwls_band_index
202 
203 
204 !second_model_parameter  = dtset%gwls_second_model_parameter
205 second_model_parameter  = zero
206 
207 
208 
209 ! set variables from gwls_GenerateEpsilon module
210 kmax   = dtset%gwls_stern_kmax
211 nseeds = dtset%gwls_nseeds
212 
213 kmax_model    = dtset%gwls_kmax_complement
214 kmax_poles    = dtset%gwls_kmax_poles
215 kmax_analytic = dtset%gwls_kmax_analytic
216 kmax_numeric  = dtset%gwls_kmax_numeric
217 
218 n_ext_freq    = dtset%gw_customnfreqsp
219 
220 use_model     = .True.
221 
222 call cpu_time(setup_time1)
223 
224 !--------------------------------------------------------------------------------
225 !
226 ! Extract the frequencies at which the integrand will be evaluated
227 ! add the value zero in the set.
228 !
229 !--------------------------------------------------------------------------------
230 
231 call generate_frequencies_and_weights(npt_gauss)
232 
233 !--------------------------------------------------------------------------------
234 !
235 !
236 ! Compute the static bases for the exact and model
237 ! dielectric operator
238 !
239 !
240 !--------------------------------------------------------------------------------
241 
242 omega_static = zero
243 ! define dimensions
244 lmax         = nseeds*kmax
245 lmax_model   = nseeds*kmax_model
246 
247 ! Allocate arrays which will contain basis
248 call setup_Lanczos_basis(lmax,lmax_model)
249 
250 ! allocate eigenvalues array
251 ABI_MALLOC(epsilon_eigenvalues_0, (lmax))
252 ABI_MALLOC(epsilon_model_eigenvalues_0, (lmax_model))
253 
254 ! set omega=0 for exact dielectric operator
255 call set_dielectric_function_frequency([0.0_dp,omega_static])
256 
257 ! and make note that the Sternheimer solutions must be kept (for use in the projected Sternheimer section).
258 if(dtset%gwls_recycle == 1) then
259   ABI_MALLOC(Sternheimer_solutions_zero,(2,npw_k,lmax,nbandv))
260   Sternheimer_solutions_zero = zero
261   write_solution = .true.
262 end if
263 if(dtset%gwls_recycle == 2) then
264   write(recy_name,'(A,I0.4,A)') "Recycling_",mpi_enreg%me,".dat"
265 
266   inquire(iolength=recy_line_size) cg(:,1:npw_k)
267 
268   inquire(file='local_tmp', exist=local_tmp_exist)
269   if(local_tmp_exist) recy_name = 'local_tmp/' // recy_name(1:118)
270 
271   if (open_file(file=recy_name,iomsg=msg,newunit=recy_unit,access='direct',form='unformatted',&
272 &               status='replace',recl=recy_line_size)/=0) then
273     ABI_ERROR(msg)
274   end if
275 
276   write_solution = .true.
277 end if
278 
279 call cpu_time(time1)
280 ! Compute the Lanczos basis using Block Lanczos; store
281 ! basis in Lbasis_lanczos
282 GWLS_TIMAB   = 1504
283 OPTION_TIMAB = 1
284 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
285 
286 call driver_generate_dielectric_matrix( matrix_function_epsilon_k, &
287 nseeds, kmax, &
288 epsilon_eigenvalues_0, &
289 Lbasis_lanczos, debug)
290 OPTION_TIMAB = 2
291 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
292 
293 call cpu_time(time2)
294 time = time2-time1
295 
296 write(timing_string,'(A)')  "Time to compute the EXACT Static Dielectric Matrix  :   "
297 call write_timing_log(timing_string,time)
298 
299 
300 call output_epsilon_eigenvalues(lmax,epsilon_eigenvalues_0,1)
301 
302 
303 ! The Sternheimer solutions at $\omega = 0$ have been stored.
304 write_solution = .false.
305 
306 
307 ! Prepare the model dielectric operator
308 call setup_Pk_model(omega_static,second_model_parameter)
309 
310 call cpu_time(time1)
311 ! Compute the Lanczos basis of the model operator using Block Lanczos; store
312 ! basis in Lbasis_model_lanczos
313 GWLS_TIMAB   = 1505
314 OPTION_TIMAB = 1
315 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
316 
317 call driver_generate_dielectric_matrix(matrix_function_epsilon_model_operator, &
318 nseeds, kmax_model, &
319 epsilon_model_eigenvalues_0, &
320 Lbasis_model_lanczos, debug)
321 OPTION_TIMAB = 2
322 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
323 
324 call cpu_time(time2)
325 time = time2-time1
326 
327 write(timing_string,'(A)')  "Time to compute the MODEL Static Dielectric Matrix  :   "
328 call write_timing_log(timing_string,time)
329 
330 
331 call output_epsilon_eigenvalues(lmax_model,epsilon_model_eigenvalues_0,2)
332 
333 
334 ABI_MALLOC(eigenvalues_static_eps_model_m1_minus_one, (lmax_model))
335 
336 do l = 1, lmax_model
337 eigenvalues_static_eps_model_m1_minus_one(l) = one/epsilon_model_eigenvalues_0(l)-one
338 end do
339 
340 
341 
342 !--------------------------------------------------------------------------------
343 !
344 !
345 ! Prepare and compute the projection of the dielectric Sternheimer equations
346 !
347 !
348 !--------------------------------------------------------------------------------
349 
350 !  Setup various arrays necessary for the Sternheimer projection scheme
351 nfrequencies = dtset%gwls_n_proj_freq
352 ABI_MALLOC(list_projection_frequencies,(nfrequencies))
353 
354 list_projection_frequencies = dtset%gwls_list_proj_freq
355 
356 call cpu_time(time1)
357 GWLS_TIMAB   = 1507
358 OPTION_TIMAB = 1
359 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
360 
361 !call setup_projected_Sternheimer_epsilon(lmax, npt_gauss, second_model_parameter, &
362 !                                        list_projection_frequencies,nfrequencies,debug)
363 
364 
365 call ProjectedSternheimerEpsilon(lmax, npt_gauss, second_model_parameter, &
366 list_projection_frequencies,nfrequencies,&
367 epsilon_eigenvalues_0,debug,use_model)
368 
369 
370 
371 
372 ! The Sternheimer solutions at $\omega = 0$ have been used to make the basis for the projected Sternheimer equations.
373 if(dtset%gwls_recycle == 1) then
374   ABI_FREE(Sternheimer_solutions_zero)
375 end if
376 if(dtset%gwls_recycle == 2) then
377   close(recy_unit,status='delete')
378 end if
379 
380 OPTION_TIMAB = 2
381 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
382 
383 
384 call cpu_time(time2)
385 time = time2-time1
386 write(timing_string,'(A)')  "Time to setup and compute the projected Sternheimer epsilon     :   "
387 call write_timing_log(timing_string,time)
388 
389 
390 call cpu_time(time1)
391 
392 GWLS_TIMAB   = 1508
393 OPTION_TIMAB = 1
394 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
395 
396 call compute_eps_m1_minus_eps_model_m1(lmax, npt_gauss)
397 
398 OPTION_TIMAB = 2
399 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
400 
401 call cpu_time(time2)
402 time = time2-time1
403 write(timing_string,'(A)')  "Time to compute eps^{-1}-eps_model^{-1}             :   "
404 call write_timing_log(timing_string,time)
405 
406 
407 !--------------------------------------------------------------------------------
408 !
409 !
410 ! Compute the model dielectric array
411 !
412 !
413 !--------------------------------------------------------------------------------
414 
415 call cpu_time(time1)
416 
417 GWLS_TIMAB   = 1509
418 OPTION_TIMAB = 1
419 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
420 
421 call compute_eps_model_m1_minus_one(lmax_model, npt_gauss, second_model_parameter, epsilon_model_eigenvalues_0)
422 
423 OPTION_TIMAB = 2
424 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
425 ABI_FREE(epsilon_model_eigenvalues_0)
426 
427 
428 call cpu_time(time2)
429 time = time2-time1
430 write(timing_string,'(A)')  "Time to compute eps_model^{-1}-1                    :   "
431 call write_timing_log(timing_string,time)
432 
433 
434 !--------------------------------------------------------------------------------
435 !
436 !
437 ! We no longer need the Lanczos basis in its current form!
438 ! Modify the basis so that it now contains (V^1/2.L)^*.psie
439 !
440 !
441 !--------------------------------------------------------------------------------
442 
443 call cpu_time(time1)
444 ABI_MALLOC(psie_k, (2,npw_k))
445 
446 GWLS_TIMAB   = 1510
447 OPTION_TIMAB = 1
448 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
449 
450 
451 psie_k = cg(:,(e-1)*npw_k+1:e*npw_k)
452 
453 call modify_Lbasis_Coulomb(psie_k, lmax, lmax_model)
454 
455 ABI_FREE(psie_k)
456 
457 OPTION_TIMAB = 2
458 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
459 
460 
461 call cpu_time(time2)
462 time = time2-time1
463 write(timing_string,'(A)')  "Time to modify the Lanczos basis                    :   "
464 call write_timing_log(timing_string,time)
465 
466 !--------------------------------------------------------------------------------
467 !
468 !
469 !
470 ! Diagonalize the static array eps^{-1} - eps^{-1}_model in order to
471 ! be able to apply diagonal shift Lanczos.
472 !
473 !--------------------------------------------------------------------------------
474 
475 call cpu_time(time1)
476 ! diagonalize the static eps^{-1} - eps^{-1}_model array, so as the use the diagonal Lanczos procedure
477 ! for the analytical term
478 GWLS_TIMAB   = 1511
479 OPTION_TIMAB = 1
480 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
481 
482 ABI_MALLOC(Lbasis_diagonalize_dielectric_terms, (npw_k,lmax))
483 ABI_MALLOC(hermitian_static_eps_m1_minus_eps_model_m1, (lmax,lmax))
484 
485 ABI_MALLOC(eigenvalues_static_eps_m1_minus_eps_model_m1, (lmax))
486 
487 ABI_MALLOC(rwork, (3*lmax-2))
488 
489 
490 hermitian_static_eps_m1_minus_eps_model_m1(:,:) = eps_m1_minus_eps_model_m1(:,:,1)
491 
492 
493 
494 ! WORK QUERRY
495 lwork = -1
496 ABI_MALLOC(work, (1))
497 call ZHEEV( 'V',        & ! Compute eigenvectors and eigenvalues
498 'U',        & ! use Upper triangular part
499 lmax,        & ! order of matrix
500 hermitian_static_eps_m1_minus_eps_model_m1,  & ! initial matrix on input; eigenvectors on output
501 lmax,        & ! LDA
502 eigenvalues_static_eps_m1_minus_eps_model_m1,& ! eigenvalues
503 work, lwork, rwork, info )  ! work stuff
504 
505 if ( info /= 0) then
506   debug_unit = get_unit()
507   write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
508 
509   open(debug_unit,file=trim(debug_filename),status='unknown')
510 
511   write(debug_unit,'(A)')      '*********************************************************************************************'
512   write(debug_unit,'(A,I4,A)') '*      ERROR: info = ',info,' in ZHEEV (1), gwls_ComputeCorrelationEnergy'
513   write(debug_unit,'(A)')      '*********************************************************************************************'
514 
515   close(debug_unit)
516 
517 end if
518 
519 ! COMPUTATION
520 lwork = nint(dble(work(1)))
521 ABI_FREE(work)
522 ABI_MALLOC(work, (lwork))
523 call ZHEEV( 'V',        & ! Compute eigenvectors and eigenvalues
524 'U',        & ! use Upper triangular part
525 lmax,        & ! order of matrix
526 hermitian_static_eps_m1_minus_eps_model_m1,  & ! initial matrix on input; eigenvectors on output
527 lmax,        & ! LDA
528 eigenvalues_static_eps_m1_minus_eps_model_m1,& ! eigenvalues
529 work, lwork, rwork, info )  ! work stuff
530 
531 if ( info /= 0) then
532   debug_unit = get_unit()
533   write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
534 
535   open(debug_unit,file=trim(debug_filename),status='unknown')
536 
537   write(debug_unit,'(A)')      '*********************************************************************************************'
538   write(debug_unit,'(A,I4,A)') '*      ERROR: info = ',info,' in ZHEEV (2), gwls_ComputeCorrelationEnergy'
539   write(debug_unit,'(A)')      '*********************************************************************************************'
540 
541   close(debug_unit)
542 
543 end if
544 
545 
546 
547 ABI_FREE(work)
548 ABI_FREE(rwork)
549 
550 !--------------------------------------------------------------------------------
551 !
552 ! update basis: L' = L . Q
553 !
554 !
555 ! CAREFUL!!! We must multiply by conjg(hermitian_static_eps_m1_minus_eps_model_m1),
556 !            which is the COMPLEX CONJUGATE of the eigenvectors of the matrix
557 !            eps^{-1}-eps_m^{-1}, because we have MODIFIED the basis Lbasis_lanczos
558 !            to contain the complex conjugate of the eigenvectors of eps.
559 !            This is somewhat subtle, but forgetting to do this leads to small errors
560 !            in the results...
561 !--------------------------------------------------------------------------------
562 hermitian_static_eps_m1_minus_eps_model_m1 = conjg(hermitian_static_eps_m1_minus_eps_model_m1)
563 
564 call ZGEMM('N','N',npw_k,lmax,lmax,cmplx_1,Lbasis_lanczos,npw_k,  &
565 hermitian_static_eps_m1_minus_eps_model_m1, &
566 lmax,cmplx_0,Lbasis_diagonalize_dielectric_terms,npw_k)
567 
568 
569 OPTION_TIMAB = 2
570 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
571 
572 
573 call cpu_time(time2)
574 time = time2-time1
575 write(timing_string,'(A)')  "Time to diagonalize eps^{-1}(0)-eps^{-1}(0)_model   :   "
576 call write_timing_log(timing_string,time)
577 
578 ABI_FREE(hermitian_static_eps_m1_minus_eps_model_m1)
579 
580 
581 call cpu_time(setup_time2)
582 setup_time = setup_time2 - setup_time1
583 
584 write(timing_string,'(A)')  "       TOTAL DIELECTRIC SETUP TIME                  :   "
585 call write_timing_log(timing_string,setup_time)
586 
587 !--------------------------------------------------------------------------------
588 !
589 ! Compute the Analytic energy using Shift Lanczos
590 !
591 !--------------------------------------------------------------------------------
592 
593 call cpu_time(time1)
594 
595 GWLS_TIMAB   = 1512
596 OPTION_TIMAB = 1
597 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
598 
599 
600 ABI_MALLOC(AT_Lanczos,(n_ext_freq,lmax))
601 call compute_AT_shift_Lanczos(n_ext_freq,dtset%gw_freqsp, model_parameter, lmax, Lbasis_diagonalize_dielectric_terms,&
602 &                             kmax_analytic, AT_Lanczos)
603 
604 OPTION_TIMAB = 2
605 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
606 
607 
608 ABI_FREE(Lbasis_diagonalize_dielectric_terms)
609 
610 call cpu_time(time2)
611 time = time2 - time1
612 
613 write(timing_string,'(A)')  "Time to compute analytical term by SHIFT LANCZOS    :   "
614 call write_timing_log(timing_string,time)
615 
616 
617 call cpu_time(time1)
618 
619 GWLS_TIMAB   = 1513
620 OPTION_TIMAB = 1
621 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
622 
623 ABI_MALLOC(AT_model_Lanczos,(n_ext_freq,lmax_model))
624 call compute_AT_shift_Lanczos(n_ext_freq,dtset%gw_freqsp, model_parameter, lmax_model,  &
625 Lbasis_model_lanczos, kmax_analytic, AT_model_Lanczos)
626 
627 OPTION_TIMAB = 2
628 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
629 
630 call cpu_time(time2)
631 time = time2 - time1
632 
633 write(timing_string,'(A)')  "Time to compute analytical MODEL by SHIFT LANCZOS   :  "
634 call write_timing_log(timing_string,time)
635 
636 
637 !--------------------------------------------------------------------------------
638 !
639 ! Compute the Numeric energy using Shift Lanczos
640 !
641 !--------------------------------------------------------------------------------
642 
643 call cpu_time(time1)
644 
645 ABI_MALLOC(array_integrand_exact_sector,(npt_gauss+1,n_ext_freq))
646 
647 ABI_MALLOC( tmp_dielectric_array, (lmax,lmax,npt_gauss+1))
648 
649 do iw = 1, npt_gauss + 1
650 
651 lorentzian      = model_parameter**2/(list_omega(iw)**2+model_parameter**2)
652 tmp_dielectric_array(:,:,iw) = eps_m1_minus_eps_model_m1(:,:,iw)-lorentzian*eps_m1_minus_eps_model_m1(:,:,1)
653 
654 end do
655 GWLS_TIMAB   = 1514
656 OPTION_TIMAB = 1
657 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
658 
659 call compute_projected_BT_shift_Lanczos(n_ext_freq, dtset%gw_freqsp, lmax, Lbasis_lanczos,         &
660 kmax_numeric, npt_gauss, tmp_dielectric_array, array_integrand_exact_sector )
661 
662 OPTION_TIMAB = 2
663 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
664 
665 ABI_FREE( tmp_dielectric_array)
666 call cpu_time(time2)
667 time = time2 - time1
668 write(timing_string,'(A)')  "Time to compute numerical term by SHIFT LANCZOS     :   "
669 call write_timing_log(timing_string,time)
670 
671 
672 
673 call cpu_time(time1)
674 
675 ABI_MALLOC(array_integrand_model_sector,(npt_gauss+1,n_ext_freq))
676 
677 !ABI_MALLOC( tmp_dielectric_array, (lmax_model,lmax_model,npt_gauss+1))
678 ABI_MALLOC( tmp_dielectric_array, (lmax_model,blocksize_epsilon,npt_gauss+1))
679 
680 
681 do iw = 1, npt_gauss + 1
682 
683 lorentzian      = model_parameter**2/(list_omega(iw)**2+model_parameter**2)
684 !tmp_dielectric_array(:,:,iw) = eps_model_m1_minus_one(:,:,iw)-lorentzian*eps_model_m1_minus_one(:,:,1)
685 tmp_dielectric_array(:,:,iw) = eps_model_m1_minus_one_DISTR(:,:,iw)-lorentzian*eps_model_m1_minus_one_DISTR(:,:,1)
686 
687 end do
688 
689 GWLS_TIMAB   = 1515
690 OPTION_TIMAB = 1
691 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
692 
693 !call compute_projected_BT_shift_Lanczos(n_ext_freq , dtset%gw_freqsp, lmax_model, Lbasis_model_lanczos,         &
694 !                                        kmax_numeric, npt_gauss, tmp_dielectric_array, array_integrand_model_sector )
695 
696 
697 call compute_projected_BT_shift_Lanczos_DISTRIBUTED(n_ext_freq, dtset%gw_freqsp, lmax_model, blocksize_epsilon, &
698 model_lanczos_vector_belongs_to_this_node, model_lanczos_vector_index,  &
699 Lbasis_model_lanczos, kmax_numeric, npt_gauss, tmp_dielectric_array, &
700 array_integrand_model_sector )
701 
702 OPTION_TIMAB = 2
703 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
704 
705 ABI_FREE(tmp_dielectric_array)
706 call cpu_time(time2)
707 time = time2 - time1
708 write(timing_string,'(A)')  "Time to compute numerical model   SHIFT LANCZOS     :   "
709 call write_timing_log(timing_string,time)
710 
711 
712 !--------------------------------------------------------------------------------
713 !
714 ! set up arrays for poles
715 !
716 !--------------------------------------------------------------------------------
717 
718 
719 
720 !call generate_degeneracy_table_for_poles(debug) ! so we can compute Poles contributions
721 call generate_degeneracy_table_for_poles(.true.) ! so we can compute Poles contributions
722 
723 !--------------------------------------------------------------------------------
724 !
725 ! Print contributions to sigma_A_Lanczos  to a file
726 !
727 !--------------------------------------------------------------------------------
728 call output_Sigma_A_by_eigenvalues(n_ext_freq,lmax,dtset%gw_freqsp,AT_Lanczos,eigenvalues_static_eps_m1_minus_eps_model_m1,2)
729 call output_Sigma_A_by_eigenvalues(n_ext_freq,lmax_model,dtset%gw_freqsp,AT_model_Lanczos,&
730 &                                  eigenvalues_static_eps_model_m1_minus_one,3)
731 
732 !epsilon_eigenvalues_0
733 
734 ABI_FREE(epsilon_eigenvalues_0)
735 !--------------------------------------------------------------------------------
736 !
737 ! Iterate on external frequencies
738 !
739 !--------------------------------------------------------------------------------
740 
741 
742 do iw_ext = 1, dtset%gw_customnfreqsp
743 
744 call cpu_time(freq_time1)
745 external_omega = dtset%gw_freqsp(iw_ext)
746 
747 write(timing_string,'(A)')  "#"
748 call write_text_block_in_Timing_log(timing_string)
749 write(timing_string,'(A)')  "#"
750 call write_text_block_in_Timing_log(timing_string)
751 write(timing_string,'(A,I4,A,F8.4,A)')  "#  Frequency # ",iw_ext," omega = ",external_omega," Ha"
752 call write_text_block_in_Timing_log(timing_string)
753 write(timing_string,'(A)')  "#"
754 call write_text_block_in_Timing_log(timing_string)
755 write(timing_string,'(A)')  "#"
756 call write_text_block_in_Timing_log(timing_string)
757 
758 
759 !--------------------------------------------------------------------------------
760 !
761 ! compute the pole term
762 !                        CAREFUL! The real valence states must still be allocated
763 !                        for the dielectric operator to work properly
764 !
765 !--------------------------------------------------------------------------------
766 
767 call cpu_time(time1)
768 GWLS_TIMAB   = 1516
769 OPTION_TIMAB = 1
770 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
771 
772 pole_energy    = compute_Poles(external_omega,kmax_poles,debug)
773 
774 OPTION_TIMAB = 2
775 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
776 
777 call cpu_time(time2)
778 
779 time = time2-time1
780 write(timing_string,'(A)')  "Time to compute the Poles contribution              :   "
781 call write_timing_log(timing_string,time)
782 
783 
784 
785 !================================================================================
786 ! Compute the contributions from the analytic term
787 !================================================================================
788 GWLS_TIMAB   = 1517
789 OPTION_TIMAB = 1
790 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
791 
792 
793 call cpu_time(time1)
794 
795 sigma_A_Lanczos      = dble(sum(AT_Lanczos(iw_ext,:)*eigenvalues_static_eps_m1_minus_eps_model_m1(:)))
796 
797 call cpu_time(time2)
798 time = time2-time1
799 write(timing_string,'(A)')  "Time  Tr[(eps^{-1}-eps_model^{-1}).AT] AFTER SHIFT  :   "
800 call write_timing_log(timing_string,time)
801 
802 
803 call cpu_time(time1)
804 
805 sigma_A_model_Lanczos= dble(sum(AT_model_Lanczos(iw_ext,:)*eigenvalues_static_eps_model_m1_minus_one(:)))
806 
807 call cpu_time(time2)
808 time = time2-time1
809 write(timing_string,'(A)')  "Time for Tr[ (eps_model^{-1}-1) . AT ] AFTER SHIFT  :   "
810 call write_timing_log(timing_string,time)
811 
812 OPTION_TIMAB = 2
813 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
814 !--------------------------------------------------------------------------------
815 !
816 ! compute integrand
817 !
818 !--------------------------------------------------------------------------------
819 GWLS_TIMAB   = 1518
820 OPTION_TIMAB = 1
821 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
822 
823 call compute_integrands_shift_lanczos(iw_ext, n_ext_freq, npt_gauss, array_integrand_exact_sector, &
824 array_integrand_model_sector, sigma_B_Lanczos, sigma_B_model_Lanczos)
825 
826 OPTION_TIMAB = 2
827 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
828 
829 !--------------------------------------------------------------------------------
830 !
831 ! Output results
832 !
833 !--------------------------------------------------------------------------------
834 
835 
836 call output_results(iw_ext,npt_gauss, lmax,lmax_model, model_parameter, second_model_parameter,  &
837 external_omega, Sigma_x,Vxc_energy,pole_energy,                          &
838 sigma_A_Lanczos,sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos)
839 
840 
841 call cpu_time(freq_time2)
842 freq_time = freq_time2-freq_time1
843 
844 write(timing_string,'(A)')  "               TOTAL FREQUENCY TIME                 :   "
845 call write_timing_log(timing_string,freq_time)
846 
847 
848 
849 correlations        = pole_energy+sigma_A_Lanczos+sigma_A_model_Lanczos+sigma_B_Lanczos+sigma_B_model_Lanczos
850 
851 renormalized_energy = eig(e) + Sigma_x-Vxc_energy +correlations
852 write(std_out,10) '                               '
853 write(std_out,14) ' For omega                   : ',external_omega     ,' Ha = ',external_omega     *Ha_eV,' eV'
854 write(std_out,14) '   <psi_e | Sigma_c  | psi_e>: ',correlations       ,' Ha = ',correlations       *Ha_eV,' eV'
855 write(std_out,14) '   eps_e + <Sigma_xc - V_xc> : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV'
856 
857 write(ab_out,10) '                               '
858 write(ab_out,14) ' For omega                   : ',external_omega     ,' Ha = ',external_omega     *Ha_eV,' eV'
859 write(ab_out,14) '   <psi_e | Sigma_c  | psi_e>: ',correlations       ,' Ha = ',correlations       *Ha_eV,' eV'
860 write(ab_out,14) '   eps_e + <Sigma_xc - V_xc> : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV'
861 
862 
863 end do
864 
865 ABI_FREE(AT_Lanczos)
866 ABI_FREE(AT_model_Lanczos)
867 ABI_FREE(array_integrand_exact_sector)
868 ABI_FREE(array_integrand_model_sector)
869 ABI_FREE(eigenvalues_static_eps_model_m1_minus_one)
870 ABI_FREE(eigenvalues_static_eps_m1_minus_eps_model_m1)
871 call clean_degeneracy_table_for_poles()
872 call cleanup_Pk_model()
873 call cleanup_Lanczos_basis()
874 call cleanup_projected_Sternheimer_epsilon()
875 
876 call cpu_time(total_time2)
877 total_time = total_time2-total_time1
878 write(timing_string,'(A)')  "               TOTAL TIME                           :   "
879 call write_timing_log(timing_string,total_time)
880 
881 
882 10 format(A)
883 14 format(A,ES24.16,A,F16.8,A)
884 
885 end subroutine compute_correlations_shift_lanczos

m_gwls_ComputeCorrelationEnergy/compute_integrands_shift_lanczos [ Functions ]

[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  compute_integrands_shift_lanczos

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

1420 subroutine compute_integrands_shift_lanczos(iw_ext,n_ext_freq,npt_gauss, array_integrand_exact_sector, &
1421 array_integrand_model_sector, sigma_B_Lanczos, sigma_B_model_Lanczos)
1422 !----------------------------------------------------------------------------------------------------
1423 !
1424 ! This subroutine computes the integrands, assuming data was generated by shift lanczos.
1425 !----------------------------------------------------------------------------------------------------
1426 
1427 integer,  intent(in)   :: iw_ext, npt_gauss, n_ext_freq
1428 
1429 complex(dpc),  intent(in)   :: array_integrand_exact_sector(npt_gauss+1,n_ext_freq)
1430 complex(dpc),  intent(in)   :: array_integrand_model_sector(npt_gauss+1,n_ext_freq)
1431 
1432 real(dp), intent(out)  :: sigma_B_Lanczos, sigma_B_model_Lanczos
1433 
1434 real(dp) :: integrand_Lanczos , integrand_model_Lanczos
1435 real(dp) :: omega_prime
1436 
1437 integer         :: iw
1438 integer         :: io_unit
1439 character(256)  :: title_string
1440 character(256)  :: timing_string
1441 
1442 
1443 real(dp)       :: time1, time2, time
1444 
1445 ! *************************************************************************
1446 
1447 
1448 if (mpi_enreg%me == 0) then
1449   io_unit = get_unit()
1450 
1451   if (iw_ext < 10) then
1452     write(title_string,'(A,I1,A)')  'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat'
1453   else if (iw_ext < 100) then
1454     write(title_string,'(A,I2,A)')  'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat'
1455   else
1456     write(title_string,'(A,I3,A)')  'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat'
1457   end if
1458 
1459   open(file=title_string,status=files_status_new,unit=io_unit)
1460   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1461   write(io_unit,10) '#'
1462   write(io_unit,10) '#               Approximate Integrands as a function of frequency'
1463   write(io_unit,10) '#'
1464   write(io_unit,10) '#        I1 = Tr[ (eps^{-1}(iw) - eps_model^{-1}(iw)) -                      '
1465   write(io_unit,10) '#                          f(w)(eps^{-1}(0) - eps_model^{-1}(0)) BT(w) ]'
1466   write(io_unit,10) '#'
1467   write(io_unit,10) '#        I2 = Tr[ (eps_model^{-1}(iw) - 1) - f(w)(eps_model^{-1}(0)-1) BT(w)]'
1468   write(io_unit,10) '#                                                                            '
1469   write(io_unit,10) '#         DIAG[I] will represent the contribution coming from taking only the diagonal elements of'
1470   write(io_unit,10) '#         the arrays in the trace.'
1471   write(io_unit,10) '#'
1472   write(io_unit,10) '#      omega (Ha)                  I1                      I2                 gaussian weight'
1473   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1474   flush(io_unit)
1475 end if
1476 
1477 call cpu_time(time1)
1478 
1479 sigma_B_Lanczos              = zero
1480 sigma_B_model_Lanczos        = zero
1481 
1482 do iw = 1,npt_gauss+1
1483 
1484 omega_prime     = list_omega(iw)
1485 
1486 integrand_Lanczos         =  dble(array_integrand_exact_sector(iw,iw_ext))
1487 integrand_model_Lanczos   =  dble(array_integrand_model_sector(iw,iw_ext))
1488 
1489 sigma_B_Lanczos       = sigma_B_Lanczos       + integrand_Lanczos*list_weights(iw)
1490 sigma_B_model_Lanczos = sigma_B_model_Lanczos + integrand_model_Lanczos*list_weights(iw)
1491 
1492 
1493 
1494 if (mpi_enreg%me == 0) write(io_unit,8) omega_prime, integrand_Lanczos , integrand_model_Lanczos, list_weights(iw)
1495 
1496 
1497 
1498 end do
1499 call cpu_time(time2)
1500 
1501 if (mpi_enreg%me == 0) then
1502   write(io_unit,10) ''
1503   write(io_unit,14) '# Value of the I1 integral: ',sigma_B_Lanczos ,' Ha'
1504   write(io_unit,14) '# Value of the I2 integral: ',sigma_B_model_Lanczos ,' Ha'
1505   write(io_unit,10) ''
1506   write(io_unit,10) ''
1507   close(io_unit)
1508 end if
1509 
1510 time = time2-time1
1511 write(timing_string,'(A)')  "Time compute the integrands and integrals           :   "
1512 call write_timing_log(timing_string,time)
1513 
1514 8  format(4ES24.16)
1515 10 format(A)
1516 14 format(A,ES24.16,A)
1517 
1518 
1519 end subroutine compute_integrands_shift_lanczos

m_gwls_ComputeCorrelationEnergy/output_epsilon_eigenvalues [ Functions ]

[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  output_epsilon_eigenvalues

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

1660 subroutine output_epsilon_eigenvalues(lmax,eigenvalues,which_case)
1661 !----------------------------------------------------------------------------------------------------
1662 !       This routine outputs the eigenvalues of the static dielectric matrix
1663 !
1664 !       There are two cases to consider:
1665 !                               1 )   the exact dielectric matrix
1666 !                               2 )   the model dielectric matrix
1667 !----------------------------------------------------------------------------------------------------
1668 
1669 integer,      intent(in) ::  lmax, which_case
1670 real(dp),     intent(in) :: eigenvalues(lmax)
1671 
1672 
1673 integer         :: io_unit
1674 character(128) :: filename
1675 
1676 integer         :: l
1677 
1678 ! *************************************************************************
1679 
1680 if (mpi_enreg%me == 0) then
1681   io_unit = get_unit()
1682 
1683   if (which_case == 1) then
1684     write(filename,'(A)') "EPSILON_EIGENVALUES.dat"
1685   else if (which_case == 2) then
1686     write(filename,'(A)') "MODEL_EPSILON_EIGENVALUES.dat"
1687   end if
1688 
1689 
1690   open(file=filename,status=files_status_new,unit=io_unit)
1691   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1692   write(io_unit,10) '#'
1693   write(io_unit,10) '#  This file contains the computed eigenvalues of the static dielectric operator'
1694   write(io_unit,10) '#  either exact or model, as indicated by the name of this file.'
1695   write(io_unit,10) '#'
1696   write(io_unit,10) '#     l                           epsilon_l                                                     '
1697   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1698 
1699   do l = 1, lmax
1700 
1701   write(io_unit,20) l, eigenvalues(l)
1702   end do
1703 
1704 
1705   close(io_unit)
1706 
1707 
1708 end if
1709 
1710 10 format(A)
1711 20 format(I7,20X,ES24.16)
1712 
1713 
1714 end subroutine output_epsilon_eigenvalues

m_gwls_ComputeCorrelationEnergy/output_results [ Functions ]

[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  output_results

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

1535 subroutine output_results(iw_ext,npt_gauss, lmax,lmax_model, model_parameter, second_model_parameter,  external_omega, &
1536 Sigma_x,Vxc_energy,pole_energy,sigma_A_Lanczos,sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos,&
1537 Sigma_x_Lanczos_projected )
1538 !----------------------------------------------------------------------------------------------------
1539 !
1540 ! This subroutine computes the integrands
1541 !----------------------------------------------------------------------------------------------------
1542 
1543 integer,  intent(in)   :: iw_ext,lmax,lmax_model,npt_gauss
1544 real(dp), intent(in)   :: model_parameter, second_model_parameter,  external_omega
1545 real(dp), intent(in)   :: Sigma_x,Vxc_energy,pole_energy,sigma_A_Lanczos
1546 real(dp), intent(in)   :: sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos
1547 
1548 real(dp), optional, intent(in)   :: Sigma_x_Lanczos_projected
1549 
1550 integer         :: io_unit
1551 
1552 character(128) :: filename
1553 
1554 real(dp)       :: Sigma_c
1555 
1556 ! *************************************************************************
1557 
1558 
1559 if (mpi_enreg%me == 0) then
1560   io_unit = get_unit()
1561   if (iw_ext < 10) then
1562     write(filename,'(A,I1,A)')  'ALL_ENERGY_',iw_ext,'.dat'
1563   else if (iw_ext < 100) then
1564     write(filename,'(A,I2,A)')  'ALL_ENERGY_',iw_ext,'.dat'
1565   else
1566     write(filename,'(A,I3,A)')  'ALL_ENERGY_',iw_ext,'.dat'
1567   end if
1568 
1569 
1570   open(file=filename,status=files_status_new,unit=io_unit)
1571   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1572   write(io_unit,10) '#'
1573   write(io_unit,10) '#               This file contains the results of the Correlation energy calculation.'
1574   write(io_unit,10) '# '
1575   write(io_unit,10) '# Definitions:'
1576   write(io_unit,10) '# '
1577   write(io_unit,10) '#                eps_e     = Bare DFT energy of the state                                           '
1578   write(io_unit,10) '# '
1579   write(io_unit,10) '#                Sigma_A_1 = Tr[(eps^{-1}(0) - eps_model^{-1}(0)) AT(W) ]                           '
1580   write(io_unit,10) '# '
1581   write(io_unit,10) '#                Sigma_A_2 = Tr[(eps_model^{-1}(0) - 1) AT(W) ]                           '
1582   write(io_unit,10) '# '
1583   write(io_unit,10) '#                Sigma_B_1 = Int dw I1(w)                                                '
1584   write(io_unit,10) '# '
1585   write(io_unit,10) '#                Sigma_B_2 = Int dw I2(w)                                                '
1586   write(io_unit,10) '# '
1587   write(io_unit,10) '#        I1(w) = Tr[ (eps^{-1}(iw) - eps_model^{-1}(iw)) -                      '
1588   write(io_unit,10) '#                          f(w)(eps^{-1}(0) - eps_model^{-1}(0)) BT(w,W) ]'
1589   write(io_unit,10) '#'
1590   write(io_unit,10) '#        I2(w) = Tr[ (eps_model^{-1}(iw) - 1) - f(w)(eps_model^{-1}(0)-1) BT(w,W)]'
1591   write(io_unit,10) '#                                                                            '
1592   write(io_unit,10) '#  Parameters:                                                               '
1593   write(io_unit,10) '#                                                                            '
1594   write(io_unit,25) '#                 lmax          = ', lmax
1595   write(io_unit,25) '#                 lmax_model    = ', lmax_model
1596   write(io_unit,25) '#                 npt_gauss     = ', npt_gauss
1597   write(io_unit,12) '#                 omega0        = ', model_parameter,'  Ha'
1598   write(io_unit,12) '#                 epsilon0      = ', second_model_parameter,'  Ha'
1599   write(io_unit,14) '#                 omega_ext (W) = ', external_omega,'  Ha'
1600   write(io_unit,10) '#                                                                            '
1601   write(io_unit,10) '#                                                                            '
1602   write(io_unit,10) '#   NOTE: if lmax_model = 0, then eps_model = I, the identity.               '
1603   write(io_unit,10) '#                                                                            '
1604   write(io_unit,10) '#                                                                            '
1605   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1606   write(io_unit,10) '                                                                             '
1607   write(io_unit,30) '       eps_e    (Ha)      :     ', eig(e)
1608   write(io_unit,10) '                                                                             '
1609   write(io_unit,30) '      Sigma_x   (Ha)      :     ', Sigma_x
1610 
1611   if (present(Sigma_x_Lanczos_projected) ) then
1612     write(io_unit,30) '   Sigma_x_PROJECTED (Ha) :     ', Sigma_x_Lanczos_projected
1613   end if
1614 
1615 
1616 
1617   write(io_unit,30) '    < V_xc >_e  (Ha)      :     ', Vxc_energy
1618   write(io_unit,10) '                                                                             '
1619   write(io_unit,30) '       poles    (Ha)      :     ', pole_energy
1620   write(io_unit,30) '     Sigma_A_1  (Ha)      :     ', sigma_A_Lanczos
1621   write(io_unit,30) '     Sigma_A_2  (Ha)      :     ', sigma_A_model_Lanczos
1622   write(io_unit,30) '     Sigma_B_1  (Ha)      :     ', sigma_B_Lanczos
1623   write(io_unit,30) '     Sigma_B_2  (Ha)      :     ', sigma_B_model_Lanczos
1624   write(io_unit,10) '                                                                             '
1625 
1626   Sigma_c = pole_energy+sigma_A_Lanczos+sigma_A_model_Lanczos+sigma_B_Lanczos+sigma_B_model_Lanczos
1627 
1628   write(io_unit,30) '      Sigma_c   (Ha)      :     ',Sigma_c
1629   write(io_unit,10) '                                                                             '
1630   write(io_unit,30) '       E_e      (Ha)      :     ', eig(e)+Sigma_x-Vxc_energy+Sigma_c
1631 
1632 
1633   close(io_unit)
1634 end if
1635 
1636 
1637 
1638 10 format(A)
1639 12 format(A,ES10.3,A)
1640 14 format(A,ES24.16,A)
1641 25 format(A,I5)
1642 30 format(A,ES24.16)
1643 
1644 end subroutine output_results

m_gwls_ComputeCorrelationEnergy/output_Sigma_A_by_eigenvalues [ Functions ]

[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  output_Sigma_A_by_eigenvalues

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

1731 subroutine output_Sigma_A_by_eigenvalues(n_ext_freq,lmax,external_frequencies,AT_Lanczos,eigenvalues_array,which_case)
1732 !----------------------------------------------------------------------------------------------------
1733 !       This routine outputs the eigenvalues of the static dielectric matrix, as well as the
1734 !       contributions to Sigma_A, decomposed by eigenvalues.
1735 !
1736 !       There are three cases to consider:
1737 !                               1 )   no model is being used; we are printing A
1738 !                               2 )   a model is being used; we are printing A1
1739 !                               2 )   a model is being used; we are printing A2
1740 !----------------------------------------------------------------------------------------------------
1741 
1742 integer,      intent(in) :: n_ext_freq, lmax, which_case
1743 complex(dpc), intent(in) :: AT_Lanczos(n_ext_freq,lmax)
1744 real(dp),     intent(in) :: eigenvalues_array(lmax)
1745 real(dp),     intent(in) :: external_frequencies(n_ext_freq)
1746 
1747 integer   :: iw_ext, l
1748 real(dp)  :: external_omega
1749 
1750 
1751 complex(dpc)  :: matrix, eig
1752 
1753 integer         :: io_unit
1754 character(128) :: filename
1755 
1756 ! *************************************************************************
1757 
1758 
1759 if (mpi_enreg%me == 0) then
1760   io_unit = get_unit()
1761 
1762 
1763   do iw_ext = 1, n_ext_freq
1764 
1765   if (which_case == 1) then
1766     write(filename,'(A,I0.4,A)') "SIGMA_A_BY_EIGENVALUES_",iw_ext,".dat"
1767   else if (which_case == 2) then
1768     write(filename,'(A,I0.4,A)') "SIGMA_A1_BY_EIGENVALUES_",iw_ext,".dat"
1769   else if (which_case == 3) then
1770     write(filename,'(A,I0.4,A)') "SIGMA_A2_BY_EIGENVALUES_",iw_ext,".dat"
1771   end if
1772 
1773 
1774   external_omega = external_frequencies(iw_ext)
1775 
1776   open(file=filename,status=files_status_new,unit=io_unit)
1777   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1778   write(io_unit,10) '#'
1779   write(io_unit,10) '#  This file contains the contributions to Sigma_A, the analytical self-energy term,'
1780   write(io_unit,10) '#  as a function of the eigenvalue of the dielectric operator.'
1781   write(io_unit,10) '#'
1782   write(io_unit,10) '# Definitions:'
1783   write(io_unit,10) '# '
1784   write(io_unit,10) '# '
1785 
1786   if (which_case == 1) then
1787     write(io_unit,10) '#            Sigma_A = Tr[(eps^{-1}(0) -I ) AT(W) ]                           '
1788     write(io_unit,10) '#                    = sum_{l} ( 1/eps_l -1 ) < V_l | AT(W) | V_l >           '
1789   else if (which_case == 2) then
1790     write(io_unit,10) '#            Sigma_A1= Tr[(eps^{-1}(0) - eps_model^{-1}(0) ) A1T(W) ]         '
1791     write(io_unit,10) '#                    = sum_{l} ( LBDA_l ) < V_l | A1T(W) | V_l >          '
1792     write(io_unit,10) '#                    where LBDA_l are the eigenvalues of eps^{-1}(0)-eps_model^{-1}(0) '
1793   else if (which_case == 3) then
1794     write(io_unit,10) '#            Sigma_A2= Tr[(eps_model^{-1}(0)- I ) A2T(W) ]         '
1795     write(io_unit,10) '#                    = sum_{l} (1/eps^{model}_l -1 ) < V_l | A2T(W) | V_l >          '
1796   end if
1797 
1798 
1799   write(io_unit,10) '# '
1800   write(io_unit,10) '#                                                                            '
1801   write(io_unit,10) '#  Parameters:                                                               '
1802   write(io_unit,10) '#                                                                            '
1803 
1804   if (which_case == 1 .or. which_case == 2) then
1805     write(io_unit,25) '#                 lmax          = ', lmax
1806   else if (which_case == 3) then
1807     write(io_unit,25) '#                 lmax_model    = ', lmax
1808   end if
1809 
1810   write(io_unit,25) '#                 n_ext_freq    = ', n_ext_freq
1811   write(io_unit,25) '#                 iw_ext        = ', iw_ext
1812   write(io_unit,14) '#                 omega_ext (W) = ', external_omega,'  Ha'
1813   write(io_unit,10) '#                                                                            '
1814   write(io_unit,10) '#                                                                            '
1815   if (which_case == 1) then
1816     write(io_unit,10) '#     (1/eps_l -1)                  < V_l | AT(W) | V_l >  (Ha)                ( 1/eps_l -1 ) '&
1817     &                     //'< V_l | AT(W) | V_l > (Ha)'
1818   else if (which_case == 2) then
1819     write(io_unit,10) '#        LBDA_1                     < V_l | A1T(W) | V_l >  (Ha)                   LBDA_l     '&
1820     &                     //'< V_l | A1T(W) | V_l > (Ha)'
1821   else if (which_case == 3) then
1822     write(io_unit,10) '#     (1/eps_m_l -1)                < V_l | A2T(W) | V_l >  (Ha)               ( 1/eps_m_l -1 ) '&
1823     &                     //'< V_l | A2T(W) | V_l > (Ha)'
1824   end if
1825 
1826   write(io_unit,10) '#                                 real                 imaginary                    real                  '&
1827   &                   //'imaginary'
1828   write(io_unit,10) '#---------------------------------------------------------------------------------------------------------'&
1829   &                   //'----------------'
1830 
1831   do l = 1, lmax
1832 
1833   matrix = AT_Lanczos(iw_ext,l)
1834   eig    = eigenvalues_array(l)
1835 
1836 
1837   write(io_unit,20) real(eig), matrix, eig*matrix
1838   end do
1839 
1840 
1841   close(io_unit)
1842 
1843   end do
1844 
1845 end if
1846 
1847 10 format(A)
1848 14 format(A,ES24.16,A)
1849 20 format(ES24.16,2ES24.16,2X,2ES24.16)
1850 25 format(A,I5)
1851 
1852 
1853 
1854 end subroutine output_Sigma_A_by_eigenvalues