TABLE OF CONTENTS
- ABINIT/gwls_ComputeCorrelationEnergy
- gwls_ComputeCorrelationEnergy/compute_correlations_no_model_shift_lanczos
- gwls_ComputeCorrelationEnergy/compute_correlations_shift_lanczos
- gwls_ComputeCorrelationEnergy/compute_integrands_shift_lanczos
- gwls_ComputeCorrelationEnergy/output_epsilon_eigenvalues
- gwls_ComputeCorrelationEnergy/output_results
- gwls_ComputeCorrelationEnergy/output_Sigma_A_by_eigenvalues
ABINIT/gwls_ComputeCorrelationEnergy [ 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