TABLE OF CONTENTS


ABINIT/m_gwls_ComputeCorrelationEnergy [ Modules ]

[ Top ] [ Modules ]

NAME

 m_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 m_gwls_ComputeCorrelationEnergy
30 
31 ! local modules
32 use m_gwls_utility
33 use m_gwls_wf
34 use m_gwls_hamiltonian
35 use m_gwls_lineqsolver
36 use m_gwls_polarisability
37 use m_gwls_model_polarisability
38 use m_gwls_DielectricArray
39 use m_gwls_ComputePoles
40 use m_gwls_Projected_AT
41 use m_gwls_Projected_BT
42 use m_gwls_GWlanczos
43 use m_gwls_GenerateEpsilon
44 use m_gwls_GWanalyticPart
45 use m_gwls_TimingLog
46 use m_gwls_LanczosBasis
47 
48 ! abinit modules
49 use defs_basis
50 use defs_datatypes
51 use defs_abitypes
52 use defs_wvltypes
53 use m_abicore
54 use m_xmpi
55 use m_pawang
56 use m_errors
57 
58 use m_time,             only : timab
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

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

PARENTS

      gwls_sternheimer

CHILDREN

SOURCE

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

PARENTS

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

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

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

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

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

m_gwls_ComputeCorrelationEnergy/output_results [ Functions ]

[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]

NAME

  output_results

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

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

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

SOURCE

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