TABLE OF CONTENTS


ABINIT/gwls_ComputeCorrelationEnergy [ Modules ]

[ Top ] [ Modules ]

NAME

 gwls_ComputeCorrelationEnergy

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 
28 
29 module gwls_ComputeCorrelationEnergy
30 
31 ! local modules
32 use m_gwls_utility
33 use gwls_wf
34 use gwls_hamiltonian
35 use gwls_lineqsolver
36 use gwls_polarisability
37 use gwls_model_polarisability
38 use gwls_DielectricArray
39 use gwls_ComputePoles
40 use gwls_Projected_AT
41 use gwls_Projected_BT
42 !use gwls_ComplementSpacePolarizability
43 use gwls_GWlanczos
44 use gwls_GenerateEpsilon
45 use gwls_GWanalyticPart
46 use gwls_TimingLog
47 use gwls_LanczosBasis
48 
49 ! abinit modules
50 use defs_basis
51 use defs_datatypes
52 use defs_abitypes
53 use defs_wvltypes
54 use m_profiling_abi
55 use m_xmpi
56 use m_pawang
57 use m_errors
58 
59 use m_io_tools,         only : get_unit,open_file
60 use m_paw_dmft,         only : paw_dmft_type
61 use m_ebands,           only : ebands_init, ebands_free
62 use m_gaussian_quadrature, only: gaussian_quadrature_gegenbauer, gaussian_quadrature_legendre
63 
64 
65 implicit none
66 save
67 private

gwls_ComputeCorrelationEnergy/compute_correlations_no_model_shift_lanczos [ Functions ]

[ Top ] [ gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  compute_correlations_no_model_shift_lanczos

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_sternheimer

CHILDREN

SOURCE

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

gwls_ComputeCorrelationEnergy/compute_correlations_shift_lanczos [ Functions ]

[ Top ] [ gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  compute_correlations_shift_lanczos

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_sternheimer

CHILDREN

SOURCE

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

gwls_ComputeCorrelationEnergy/compute_integrands_shift_lanczos [ Functions ]

[ Top ] [ gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  compute_integrands_shift_lanczos

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

1462 subroutine compute_integrands_shift_lanczos(iw_ext,n_ext_freq,npt_gauss, array_integrand_exact_sector, &
1463 array_integrand_model_sector, sigma_B_Lanczos, sigma_B_model_Lanczos)
1464 !----------------------------------------------------------------------------------------------------
1465 !
1466 ! This subroutine computes the integrands, assuming data was generated by shift lanczos.
1467 !----------------------------------------------------------------------------------------------------
1468 
1469 !This section has been created automatically by the script Abilint (TD).
1470 !Do not modify the following lines by hand.
1471 #undef ABI_FUNC
1472 #define ABI_FUNC 'compute_integrands_shift_lanczos'
1473 !End of the abilint section
1474 
1475 implicit none
1476 
1477 integer,  intent(in)   :: iw_ext, npt_gauss, n_ext_freq
1478 
1479 complex(dpc),  intent(in)   :: array_integrand_exact_sector(npt_gauss+1,n_ext_freq)
1480 complex(dpc),  intent(in)   :: array_integrand_model_sector(npt_gauss+1,n_ext_freq)
1481 
1482 real(dp), intent(out)  :: sigma_B_Lanczos, sigma_B_model_Lanczos
1483 
1484 real(dp) :: integrand_Lanczos , integrand_model_Lanczos 
1485 real(dp) :: omega_prime
1486 
1487 integer         :: iw 
1488 integer         :: io_unit
1489 character(256)  :: title_string
1490 character(256)  :: timing_string
1491 
1492 
1493 real(dp)       :: time1, time2, time
1494 
1495 ! *************************************************************************
1496 
1497 
1498 if (mpi_enreg%me == 0) then
1499   io_unit = get_unit()
1500 
1501   if (iw_ext < 10) then
1502     write(title_string,'(A,I1,A)')  'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat'
1503   else if (iw_ext < 100) then
1504     write(title_string,'(A,I2,A)')  'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat'
1505   else 
1506     write(title_string,'(A,I3,A)')  'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat'
1507   end if
1508 
1509   open(file=title_string,status=files_status_new,unit=io_unit)
1510   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1511   write(io_unit,10) '#'
1512   write(io_unit,10) '#               Approximate Integrands as a function of frequency'
1513   write(io_unit,10) '#'
1514   write(io_unit,10) '#        I1 = Tr[ (eps^{-1}(iw) - eps_model^{-1}(iw)) -                      '
1515   write(io_unit,10) '#                          f(w)(eps^{-1}(0) - eps_model^{-1}(0)) BT(w) ]'
1516   write(io_unit,10) '#'
1517   write(io_unit,10) '#        I2 = Tr[ (eps_model^{-1}(iw) - 1) - f(w)(eps_model^{-1}(0)-1) BT(w)]'
1518   write(io_unit,10) '#                                                                            '
1519   write(io_unit,10) '#         DIAG[I] will represent the contribution coming from taking only the diagonal elements of'
1520   write(io_unit,10) '#         the arrays in the trace.'
1521   write(io_unit,10) '#'
1522   write(io_unit,10) '#      omega (Ha)                  I1                      I2                 gaussian weight'
1523   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1524   flush(io_unit)
1525 end if
1526 
1527 call cpu_time(time1)
1528 
1529 sigma_B_Lanczos              = zero
1530 sigma_B_model_Lanczos        = zero
1531 
1532 do iw = 1,npt_gauss+1
1533 
1534 omega_prime     = list_omega(iw)
1535 
1536 integrand_Lanczos         =  dble(array_integrand_exact_sector(iw,iw_ext))
1537 integrand_model_Lanczos   =  dble(array_integrand_model_sector(iw,iw_ext))
1538 
1539 sigma_B_Lanczos       = sigma_B_Lanczos       + integrand_Lanczos*list_weights(iw)
1540 sigma_B_model_Lanczos = sigma_B_model_Lanczos + integrand_model_Lanczos*list_weights(iw)
1541 
1542 
1543 
1544 if (mpi_enreg%me == 0) write(io_unit,8) omega_prime, integrand_Lanczos , integrand_model_Lanczos, list_weights(iw)
1545 
1546 
1547 
1548 end do
1549 call cpu_time(time2)
1550 
1551 if (mpi_enreg%me == 0) then
1552   write(io_unit,10) ''
1553   write(io_unit,14) '# Value of the I1 integral: ',sigma_B_Lanczos ,' Ha'
1554   write(io_unit,14) '# Value of the I2 integral: ',sigma_B_model_Lanczos ,' Ha'
1555   write(io_unit,10) ''
1556   write(io_unit,10) ''
1557   close(io_unit)
1558 end if
1559 
1560 time = time2-time1
1561 write(timing_string,'(A)')  "Time compute the integrands and integrals           :   "
1562 call write_timing_log(timing_string,time)
1563 
1564 8  format(4ES24.16)
1565 10 format(A)
1566 14 format(A,ES24.16,A)
1567 
1568 
1569 end subroutine compute_integrands_shift_lanczos

gwls_ComputeCorrelationEnergy/output_epsilon_eigenvalues [ Functions ]

[ Top ] [ gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  output_epsilon_eigenvalues

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

1728 subroutine output_epsilon_eigenvalues(lmax,eigenvalues,which_case)
1729 !----------------------------------------------------------------------------------------------------
1730 !       This routine outputs the eigenvalues of the static dielectric matrix
1731 !
1732 !       There are two cases to consider:
1733 !                               1 )   the exact dielectric matrix
1734 !                               2 )   the model dielectric matrix
1735 !----------------------------------------------------------------------------------------------------
1736 
1737 !This section has been created automatically by the script Abilint (TD).
1738 !Do not modify the following lines by hand.
1739 #undef ABI_FUNC
1740 #define ABI_FUNC 'output_epsilon_eigenvalues'
1741 !End of the abilint section
1742 
1743 implicit none
1744 
1745 integer,      intent(in) ::  lmax, which_case
1746 real(dp),     intent(in) :: eigenvalues(lmax)
1747 
1748 
1749 integer         :: io_unit
1750 character(128) :: filename
1751 
1752 integer         :: l
1753 
1754 ! *************************************************************************
1755 
1756 if (mpi_enreg%me == 0) then
1757   io_unit = get_unit()
1758 
1759   if (which_case == 1) then
1760     write(filename,'(A)') "EPSILON_EIGENVALUES.dat"
1761   else if (which_case == 2) then
1762     write(filename,'(A)') "MODEL_EPSILON_EIGENVALUES.dat"
1763   end if
1764 
1765 
1766   open(file=filename,status=files_status_new,unit=io_unit)
1767   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1768   write(io_unit,10) '#'
1769   write(io_unit,10) '#  This file contains the computed eigenvalues of the static dielectric operator'
1770   write(io_unit,10) '#  either exact or model, as indicated by the name of this file.'
1771   write(io_unit,10) '#'
1772   write(io_unit,10) '#     l                           epsilon_l                                                     '
1773   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1774 
1775   do l = 1, lmax
1776 
1777   write(io_unit,20) l, eigenvalues(l)
1778   end do 
1779 
1780 
1781   close(io_unit)
1782 
1783 
1784 end if
1785 
1786 10 format(A)
1787 20 format(I7,20X,ES24.16)
1788 
1789 
1790 end subroutine output_epsilon_eigenvalues

gwls_ComputeCorrelationEnergy/output_results [ Functions ]

[ Top ] [ gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  output_results

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

1590 subroutine output_results(iw_ext,npt_gauss, lmax,lmax_model, model_parameter, second_model_parameter,  external_omega, &
1591 Sigma_x,Vxc_energy,pole_energy,sigma_A_Lanczos,sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos,&
1592 Sigma_x_Lanczos_projected )
1593 !----------------------------------------------------------------------------------------------------
1594 !
1595 ! This subroutine computes the integrands
1596 !----------------------------------------------------------------------------------------------------
1597 
1598 !This section has been created automatically by the script Abilint (TD).
1599 !Do not modify the following lines by hand.
1600 #undef ABI_FUNC
1601 #define ABI_FUNC 'output_results'
1602 !End of the abilint section
1603 
1604 implicit none
1605 
1606 integer,  intent(in)   :: iw_ext,lmax,lmax_model,npt_gauss
1607 real(dp), intent(in)   :: model_parameter, second_model_parameter,  external_omega
1608 real(dp), intent(in)   :: Sigma_x,Vxc_energy,pole_energy,sigma_A_Lanczos
1609 real(dp), intent(in)   :: sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos
1610 
1611 real(dp), optional, intent(in)   :: Sigma_x_Lanczos_projected 
1612 
1613 integer         :: io_unit
1614 
1615 character(128) :: filename
1616 
1617 real(dp)       :: Sigma_c 
1618 
1619 ! *************************************************************************
1620 
1621 
1622 if (mpi_enreg%me == 0) then
1623   io_unit = get_unit()
1624   if (iw_ext < 10) then
1625     write(filename,'(A,I1,A)')  'ALL_ENERGY_',iw_ext,'.dat'
1626   else if (iw_ext < 100) then
1627     write(filename,'(A,I2,A)')  'ALL_ENERGY_',iw_ext,'.dat'
1628   else 
1629     write(filename,'(A,I3,A)')  'ALL_ENERGY_',iw_ext,'.dat'
1630   end if
1631 
1632 
1633   open(file=filename,status=files_status_new,unit=io_unit)
1634   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1635   write(io_unit,10) '#'
1636   write(io_unit,10) '#               This file contains the results of the Correlation energy calculation.'
1637   write(io_unit,10) '# '
1638   write(io_unit,10) '# Definitions:'
1639   write(io_unit,10) '# '
1640   write(io_unit,10) '#                eps_e     = Bare DFT energy of the state                                           '
1641   write(io_unit,10) '# '
1642   write(io_unit,10) '#                Sigma_A_1 = Tr[(eps^{-1}(0) - eps_model^{-1}(0)) AT(W) ]                           '
1643   write(io_unit,10) '# '
1644   write(io_unit,10) '#                Sigma_A_2 = Tr[(eps_model^{-1}(0) - 1) AT(W) ]                           '
1645   write(io_unit,10) '# '
1646   write(io_unit,10) '#                Sigma_B_1 = Int dw I1(w)                                                '
1647   write(io_unit,10) '# '
1648   write(io_unit,10) '#                Sigma_B_2 = Int dw I2(w)                                                '
1649   write(io_unit,10) '# '
1650   write(io_unit,10) '#        I1(w) = Tr[ (eps^{-1}(iw) - eps_model^{-1}(iw)) -                      '
1651   write(io_unit,10) '#                          f(w)(eps^{-1}(0) - eps_model^{-1}(0)) BT(w,W) ]'
1652   write(io_unit,10) '#'
1653   write(io_unit,10) '#        I2(w) = Tr[ (eps_model^{-1}(iw) - 1) - f(w)(eps_model^{-1}(0)-1) BT(w,W)]'
1654   write(io_unit,10) '#                                                                            '
1655   write(io_unit,10) '#  Parameters:                                                               '
1656   write(io_unit,10) '#                                                                            '
1657   write(io_unit,25) '#                 lmax          = ', lmax
1658   write(io_unit,25) '#                 lmax_model    = ', lmax_model
1659   write(io_unit,25) '#                 npt_gauss     = ', npt_gauss
1660   write(io_unit,12) '#                 omega0        = ', model_parameter,'  Ha'
1661   write(io_unit,12) '#                 epsilon0      = ', second_model_parameter,'  Ha'
1662   write(io_unit,14) '#                 omega_ext (W) = ', external_omega,'  Ha'
1663   write(io_unit,10) '#                                                                            '
1664   write(io_unit,10) '#                                                                            '
1665   write(io_unit,10) '#   NOTE: if lmax_model = 0, then eps_model = I, the identity.               '
1666   write(io_unit,10) '#                                                                            '
1667   write(io_unit,10) '#                                                                            '
1668   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1669   write(io_unit,10) '                                                                             '
1670   write(io_unit,30) '       eps_e    (Ha)      :     ', eig(e)
1671   write(io_unit,10) '                                                                             '
1672   write(io_unit,30) '      Sigma_x   (Ha)      :     ', Sigma_x
1673 
1674   if (present(Sigma_x_Lanczos_projected) ) then
1675     write(io_unit,30) '   Sigma_x_PROJECTED (Ha) :     ', Sigma_x_Lanczos_projected 
1676   end if 
1677 
1678 
1679 
1680   write(io_unit,30) '    < V_xc >_e  (Ha)      :     ', Vxc_energy 
1681   write(io_unit,10) '                                                                             '
1682   write(io_unit,30) '       poles    (Ha)      :     ', pole_energy
1683   write(io_unit,30) '     Sigma_A_1  (Ha)      :     ', sigma_A_Lanczos
1684   write(io_unit,30) '     Sigma_A_2  (Ha)      :     ', sigma_A_model_Lanczos
1685   write(io_unit,30) '     Sigma_B_1  (Ha)      :     ', sigma_B_Lanczos
1686   write(io_unit,30) '     Sigma_B_2  (Ha)      :     ', sigma_B_model_Lanczos
1687   write(io_unit,10) '                                                                             '
1688 
1689   Sigma_c = pole_energy+sigma_A_Lanczos+sigma_A_model_Lanczos+sigma_B_Lanczos+sigma_B_model_Lanczos
1690 
1691   write(io_unit,30) '      Sigma_c   (Ha)      :     ',Sigma_c 
1692   write(io_unit,10) '                                                                             '
1693   write(io_unit,30) '       E_e      (Ha)      :     ', eig(e)+Sigma_x-Vxc_energy+Sigma_c
1694 
1695 
1696   close(io_unit)
1697 end if
1698 
1699 
1700 
1701 10 format(A)
1702 12 format(A,ES10.3,A)
1703 14 format(A,ES24.16,A)
1704 25 format(A,I5)
1705 30 format(A,ES24.16)
1706 
1707 end subroutine output_results

gwls_ComputeCorrelationEnergy/output_Sigma_A_by_eigenvalues [ Functions ]

[ Top ] [ gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  output_Sigma_A_by_eigenvalues

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

1814 subroutine output_Sigma_A_by_eigenvalues(n_ext_freq,lmax,external_frequencies,AT_Lanczos,eigenvalues_array,which_case)
1815 !----------------------------------------------------------------------------------------------------
1816 !       This routine outputs the eigenvalues of the static dielectric matrix, as well as the 
1817 !       contributions to Sigma_A, decomposed by eigenvalues.
1818 !
1819 !       There are three cases to consider:
1820 !                               1 )   no model is being used; we are printing A
1821 !                               2 )   a model is being used; we are printing A1
1822 !                               2 )   a model is being used; we are printing A2
1823 !----------------------------------------------------------------------------------------------------
1824 
1825 !This section has been created automatically by the script Abilint (TD).
1826 !Do not modify the following lines by hand.
1827 #undef ABI_FUNC
1828 #define ABI_FUNC 'output_Sigma_A_by_eigenvalues'
1829 !End of the abilint section
1830 
1831 implicit none
1832 
1833 integer,      intent(in) :: n_ext_freq, lmax, which_case
1834 complex(dpc), intent(in) :: AT_Lanczos(n_ext_freq,lmax)
1835 real(dp),     intent(in) :: eigenvalues_array(lmax)
1836 real(dp),     intent(in) :: external_frequencies(n_ext_freq)
1837 
1838 integer   :: iw_ext, l
1839 real(dp)  :: external_omega
1840 
1841 
1842 complex(dpc)  :: matrix, eig
1843 
1844 integer         :: io_unit
1845 character(128) :: filename
1846 
1847 ! *************************************************************************
1848 
1849 
1850 if (mpi_enreg%me == 0) then
1851   io_unit = get_unit()
1852 
1853 
1854   do iw_ext = 1, n_ext_freq
1855 
1856   if (which_case == 1) then
1857     write(filename,'(A,I0.4,A)') "SIGMA_A_BY_EIGENVALUES_",iw_ext,".dat"
1858   else if (which_case == 2) then
1859     write(filename,'(A,I0.4,A)') "SIGMA_A1_BY_EIGENVALUES_",iw_ext,".dat"
1860   else if (which_case == 3) then
1861     write(filename,'(A,I0.4,A)') "SIGMA_A2_BY_EIGENVALUES_",iw_ext,".dat"
1862   end if
1863 
1864 
1865   external_omega = external_frequencies(iw_ext)
1866 
1867   open(file=filename,status=files_status_new,unit=io_unit)
1868   write(io_unit,10) '#-----------------------------------------------------------------------------------------------'
1869   write(io_unit,10) '#'
1870   write(io_unit,10) '#  This file contains the contributions to Sigma_A, the analytical self-energy term,' 
1871   write(io_unit,10) '#  as a function of the eigenvalue of the dielectric operator.'
1872   write(io_unit,10) '#'
1873   write(io_unit,10) '# Definitions:'
1874   write(io_unit,10) '# '
1875   write(io_unit,10) '# '
1876 
1877   if (which_case == 1) then
1878     write(io_unit,10) '#            Sigma_A = Tr[(eps^{-1}(0) -I ) AT(W) ]                           '
1879     write(io_unit,10) '#                    = sum_{l} ( 1/eps_l -1 ) < V_l | AT(W) | V_l >           '
1880   else if (which_case == 2) then
1881     write(io_unit,10) '#            Sigma_A1= Tr[(eps^{-1}(0) - eps_model^{-1}(0) ) A1T(W) ]         '
1882     write(io_unit,10) '#                    = sum_{l} ( LBDA_l ) < V_l | A1T(W) | V_l >          '
1883     write(io_unit,10) '#                    where LBDA_l are the eigenvalues of eps^{-1}(0)-eps_model^{-1}(0) '
1884   else if (which_case == 3) then
1885     write(io_unit,10) '#            Sigma_A2= Tr[(eps_model^{-1}(0)- I ) A2T(W) ]         '
1886     write(io_unit,10) '#                    = sum_{l} (1/eps^{model}_l -1 ) < V_l | A2T(W) | V_l >          '
1887   end if
1888 
1889 
1890   write(io_unit,10) '# '
1891   write(io_unit,10) '#                                                                            '
1892   write(io_unit,10) '#  Parameters:                                                               '
1893   write(io_unit,10) '#                                                                            '
1894 
1895   if (which_case == 1 .or. which_case == 2) then
1896     write(io_unit,25) '#                 lmax          = ', lmax
1897   else if (which_case == 3) then
1898     write(io_unit,25) '#                 lmax_model    = ', lmax
1899   end if
1900 
1901   write(io_unit,25) '#                 n_ext_freq    = ', n_ext_freq
1902   write(io_unit,25) '#                 iw_ext        = ', iw_ext
1903   write(io_unit,14) '#                 omega_ext (W) = ', external_omega,'  Ha'
1904   write(io_unit,10) '#                                                                            '
1905   write(io_unit,10) '#                                                                            '
1906   if (which_case == 1) then
1907     write(io_unit,10) '#     (1/eps_l -1)                  < V_l | AT(W) | V_l >  (Ha)                ( 1/eps_l -1 ) '&
1908     &                     //'< V_l | AT(W) | V_l > (Ha)'
1909   else if (which_case == 2) then
1910     write(io_unit,10) '#        LBDA_1                     < V_l | A1T(W) | V_l >  (Ha)                   LBDA_l     '&
1911     &                     //'< V_l | A1T(W) | V_l > (Ha)'
1912   else if (which_case == 3) then
1913     write(io_unit,10) '#     (1/eps_m_l -1)                < V_l | A2T(W) | V_l >  (Ha)               ( 1/eps_m_l -1 ) '&
1914     &                     //'< V_l | A2T(W) | V_l > (Ha)'
1915   end if
1916 
1917   write(io_unit,10) '#                                 real                 imaginary                    real                  '&
1918   &                   //'imaginary'
1919   write(io_unit,10) '#---------------------------------------------------------------------------------------------------------'&
1920   &                   //'----------------'
1921 
1922   do l = 1, lmax
1923 
1924   matrix = AT_Lanczos(iw_ext,l)
1925   eig    = eigenvalues_array(l)
1926 
1927 
1928   write(io_unit,20) real(eig), matrix, eig*matrix
1929   end do 
1930 
1931 
1932   close(io_unit)
1933 
1934   end do
1935 
1936 end if
1937 
1938 10 format(A)
1939 14 format(A,ES24.16,A)
1940 20 format(ES24.16,2ES24.16,2X,2ES24.16)
1941 25 format(A,I5)
1942 
1943 
1944 
1945 end subroutine output_Sigma_A_by_eigenvalues