TABLE OF CONTENTS
- ABINIT/m_gwls_ComputeCorrelationEnergy
- m_gwls_ComputeCorrelationEnergy/compute_correlations_no_model_shift_lanczos
- m_gwls_ComputeCorrelationEnergy/compute_correlations_shift_lanczos
- m_gwls_ComputeCorrelationEnergy/compute_integrands_shift_lanczos
- m_gwls_ComputeCorrelationEnergy/output_epsilon_eigenvalues
- m_gwls_ComputeCorrelationEnergy/output_results
- m_gwls_ComputeCorrelationEnergy/output_Sigma_A_by_eigenvalues
ABINIT/m_gwls_ComputeCorrelationEnergy [ Modules ]
NAME
m_gwls_ComputeCorrelationEnergy
FUNCTION
.
COPYRIGHT
Copyright (C) 2009-2022 ABINIT group (JLJ, BR, MC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_gwls_ComputeCorrelationEnergy 23 24 ! local modules 25 use m_gwls_utility 26 use m_gwls_wf 27 use m_gwls_hamiltonian 28 use m_gwls_lineqsolver 29 use m_gwls_polarisability 30 use m_gwls_model_polarisability 31 use m_gwls_DielectricArray 32 use m_gwls_ComputePoles 33 use m_gwls_Projected_AT 34 use m_gwls_Projected_BT 35 use m_gwls_GWlanczos 36 use m_gwls_GenerateEpsilon 37 use m_gwls_GWanalyticPart 38 use m_gwls_TimingLog 39 use m_gwls_LanczosBasis 40 41 ! abinit modules 42 use defs_basis 43 use defs_wvltypes 44 use m_abicore 45 use m_xmpi 46 use m_pawang 47 use m_errors 48 use m_dtset 49 50 use m_time, only : timab 51 use m_io_tools, only : get_unit,open_file 52 use m_paw_dmft, only : paw_dmft_type 53 use m_ebands, only : ebands_init, ebands_free 54 use m_gaussian_quadrature, only: gaussian_quadrature_gegenbauer, gaussian_quadrature_legendre 55 56 57 implicit none 58 save 59 private
m_gwls_ComputeCorrelationEnergy/compute_correlations_no_model_shift_lanczos [ Functions ]
[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]
NAME
compute_correlations_no_model_shift_lanczos
FUNCTION
.
INPUTS
OUTPUT
SOURCE
901 subroutine compute_correlations_no_model_shift_lanczos(dtset, Sigma_x,Vxc_energy,debug) 902 !---------------------------------------------------------------------------------------------------- 903 ! 904 ! 905 ! This function computes the correlation energy Sigma_c(w) for the state we wish to correct. 906 ! 907 ! this subroutine does not rely on the use of a model dielectric operator. Thus 908 ! 909 ! Sigma^A(w) = Tr[ (eps^{-1}(0)-1). A^T(w)] 910 ! Sigma^N(w) = int dw' Tr[{(eps^{-1}(w')-1)-f(w')(eps^{-1}(0)-1)}B^T(w';w)] 911 ! 912 ! Shift lanczos is used for the resolvents. 913 !---------------------------------------------------------------------------------------------------- 914 915 type(dataset_type),intent(in) :: dtset 916 917 real(dp),intent(in) :: Sigma_x, Vxc_energy 918 logical, intent(in) :: debug 919 920 921 !Local variables 922 923 real(dp):: Sigma_x_Lanczos_projected 924 integer :: npt_gauss 925 integer :: print_debug 926 927 928 integer :: kmax_poles 929 930 integer :: lmax_model 931 932 integer :: kmax_analytic 933 integer :: kmax_numeric 934 integer :: n_ext_freq 935 936 937 real(dp) :: omega_static 938 939 real(dp) :: lorentzian 940 integer :: iw_ext 941 integer :: iw 942 943 944 945 real(dp) , allocatable :: psie_k(:,:) 946 947 real(dp), allocatable :: epsilon_eigenvalues_0(:) 948 949 950 complex(dpc), allocatable :: AT_Lanczos(:,:) 951 952 953 954 real(dp) :: time1, time2, time 955 real(dp) :: total_time1,total_time2,total_time 956 real(dp) :: setup_time1, setup_time2, setup_time 957 real(dp) :: freq_time1, freq_time2, freq_time 958 959 integer :: nfrequencies 960 real(dp),allocatable :: list_projection_frequencies(:) 961 962 963 964 complex(dpc), allocatable :: array_integrand_exact_sector(:,:) 965 complex(dpc), allocatable :: array_integrand_model_sector(:,:) 966 complex(dpc), allocatable :: tmp_dielectric_array(:,:,:) 967 968 real(dp) :: external_omega 969 970 character(256) :: timing_string 971 972 integer :: recy_line_size 973 character(128) :: recy_name 974 logical :: local_tmp_exist 975 character(500) :: msg 976 977 ! Energy contributions 978 979 real(dp) :: pole_energy 980 981 real(dp) :: sigma_A_Lanczos 982 real(dp) :: sigma_A_model_Lanczos 983 real(dp) :: sigma_B_Lanczos 984 real(dp) :: sigma_B_model_Lanczos 985 986 real(dp) :: correlations 987 real(dp) :: renormalized_energy 988 989 logical :: use_model 990 991 ! ************************************************************************* 992 993 994 !-------------------------------------------------------------------------------- 995 ! 996 ! Set up variables and allocate arrays 997 ! 998 !-------------------------------------------------------------------------------- 999 1000 !Variable allocation and initialization 1001 model_number = dtset%gwls_diel_model 1002 model_parameter = dtset%gwls_model_parameter 1003 npt_gauss = dtset%gwls_npt_gauss_quad 1004 print_debug = dtset%gwls_print_debug 1005 1006 first_seed = dtset%gwls_first_seed 1007 e = dtset%gwls_band_index 1008 1009 1010 ! set variables from gwls_GenerateEpsilon module 1011 kmax = dtset%gwls_stern_kmax 1012 nseeds = dtset%gwls_nseeds 1013 1014 kmax_poles = dtset%gwls_kmax_poles 1015 1016 kmax_poles = dtset%gwls_kmax_poles 1017 kmax_analytic = dtset%gwls_kmax_analytic 1018 kmax_numeric = dtset%gwls_kmax_numeric 1019 1020 n_ext_freq = dtset%gw_customnfreqsp 1021 1022 1023 use_model = .False. 1024 1025 1026 call cpu_time(setup_time1) 1027 !-------------------------------------------------------------------------------- 1028 ! 1029 ! Extract the frequencies at which the integrand will be evaluated 1030 ! add the value zero in the set. 1031 ! 1032 !-------------------------------------------------------------------------------- 1033 1034 call generate_frequencies_and_weights(npt_gauss) 1035 1036 !-------------------------------------------------------------------------------- 1037 ! 1038 ! 1039 ! Compute the static bases for the exact dielectric operator 1040 ! 1041 ! 1042 !-------------------------------------------------------------------------------- 1043 1044 1045 omega_static = zero 1046 ! define dimensions 1047 lmax = nseeds*kmax 1048 1049 ! Allocate arrays which will contain basis 1050 ! the 0 indicates we will not use arrays for the model dielectric operator 1051 call setup_Lanczos_basis(lmax,0) 1052 1053 ! allocate eigenvalues array 1054 ABI_MALLOC(epsilon_eigenvalues_0, (lmax)) 1055 1056 ! set omega=0 for exact dielectric operator 1057 call set_dielectric_function_frequency([0.0_dp,omega_static]) 1058 1059 ! and make note that the Sternheimer solutions must be kept (for use in the projected Sternheimer section). 1060 if(dtset%gwls_recycle == 1) then 1061 ABI_MALLOC(Sternheimer_solutions_zero,(2,npw_k,lmax,nbandv)) 1062 Sternheimer_solutions_zero = zero 1063 write_solution = .true. 1064 end if 1065 if(dtset%gwls_recycle == 2) then 1066 write(recy_name,'(A,I0.4,A)') "Recycling_",mpi_enreg%me,".dat" 1067 1068 inquire(iolength=recy_line_size) cg(:,1:npw_k) 1069 1070 inquire(file='local_tmp', exist=local_tmp_exist) 1071 if(local_tmp_exist) recy_name = 'local_tmp/' // recy_name(1:118) 1072 1073 if (open_file(file=recy_name,iomsg=msg,newunit=recy_unit,access='direct',form='unformatted',& 1074 & status='replace',recl=recy_line_size)/=0) then 1075 ABI_ERROR(msg) 1076 end if 1077 1078 write_solution = .true. 1079 end if 1080 1081 1082 call cpu_time(time1) 1083 ! Compute the Lanczos basis using Block Lanczos; store 1084 ! basis in Lbasis_lanczos 1085 call driver_generate_dielectric_matrix( matrix_function_epsilon_k, & 1086 nseeds, kmax, & 1087 epsilon_eigenvalues_0, & 1088 Lbasis_lanczos, debug) 1089 call cpu_time(time2) 1090 time = time2-time1 1091 1092 write(timing_string,'(A)') "Time to compute the EXACT Static Dielectric Matrix : " 1093 call write_timing_log(timing_string,time) 1094 1095 ! The Sternheimer solutions at $\omega = 0$ have been stored. 1096 write_solution = .false. 1097 1098 1099 call output_epsilon_eigenvalues(lmax,epsilon_eigenvalues_0,1) 1100 1101 1102 ! compute the Exchange energy, when it is projected on the Lanczos basis 1103 ! CAREFUL! This must be done BEFORE we modify the lanczos basis 1104 Sigma_x_Lanczos_projected = exchange(e, Lbasis_lanczos) 1105 1106 1107 1108 !-------------------------------------------------------------------------------- 1109 ! 1110 ! 1111 ! Prepare and compute the projection of the dielectric Sternheimer equations 1112 ! 1113 ! 1114 !-------------------------------------------------------------------------------- 1115 1116 ! Setup various arrays necessary for the Sternheimer projection scheme 1117 nfrequencies = dtset%gwls_n_proj_freq 1118 ABI_MALLOC(list_projection_frequencies,(nfrequencies)) 1119 1120 list_projection_frequencies = dtset%gwls_list_proj_freq 1121 1122 call cpu_time(time1) 1123 ! The explicit "false" as the last argument is for the optional 1124 ! variable "use_model"; we are not using a model here! 1125 1126 !call setup_projected_Sternheimer_epsilon(lmax, npt_gauss, zero, & 1127 ! list_projection_frequencies,nfrequencies,debug,.false.) 1128 1129 1130 call ProjectedSternheimerEpsilon(lmax, npt_gauss, zero, & 1131 list_projection_frequencies,nfrequencies,& 1132 epsilon_eigenvalues_0,debug,use_model) 1133 1134 1135 1136 !call cpu_time(time2) 1137 !time = time2-time1 1138 !write(timing_string,'(A)') "Time to setup the projected Sternheimer epsilon : " 1139 !call write_timing_log(timing_string,time) 1140 1141 1142 ! The Sternheimer solutions at $\omega = 0$ have been used to make the basis for the projected Sternheimer equations. 1143 if(dtset%gwls_recycle == 1) then 1144 ABI_FREE(Sternheimer_solutions_zero) 1145 end if 1146 if(dtset%gwls_recycle == 2) then 1147 close(recy_unit,status='delete') 1148 end if 1149 1150 1151 !call cpu_time(time1) 1152 !call compute_projected_Sternheimer_epsilon(lmax, npt_gauss, epsilon_eigenvalues_0,debug) 1153 1154 1155 1156 call cpu_time(time2) 1157 1158 time = time2-time1 1159 write(timing_string,'(A)') "Time to compute the projected Sternheimer epsilon : " 1160 call write_timing_log(timing_string,time) 1161 1162 1163 1164 call cpu_time(time1) 1165 call compute_eps_m1_minus_one(lmax, npt_gauss) 1166 call cpu_time(time2) 1167 time = time2-time1 1168 write(timing_string,'(A)') "Time to compute eps^{-1}-I : " 1169 call write_timing_log(timing_string,time) 1170 1171 1172 !-------------------------------------------------------------------------------- 1173 ! 1174 ! 1175 ! We no longer need the Lanczos basis in its current form! 1176 ! Modify the basis so that it now contains (V^{1/2}.l) 1177 ! 1178 ! 1179 !-------------------------------------------------------------------------------- 1180 1181 call cpu_time(time1) 1182 ABI_MALLOC(psie_k, (2,npw_k)) 1183 psie_k = cg(:,(e-1)*npw_k+1:e*npw_k) 1184 1185 lmax_model = 0 1186 call modify_Lbasis_Coulomb(psie_k, lmax, lmax_model) ! lmax_model is set to zero, such that 1187 ! the model lanczos basis (which doesn't exist 1188 ! in this case) will not be modified 1189 1190 ABI_FREE(psie_k) 1191 1192 call cpu_time(time2) 1193 time = time2-time1 1194 write(timing_string,'(A)') "Time to modify the Lanczos basis : " 1195 call write_timing_log(timing_string,time) 1196 1197 1198 call cpu_time(setup_time2) 1199 setup_time = setup_time2 - setup_time1 1200 1201 write(timing_string,'(A)') " TOTAL DIELECTRIC SETUP TIME : " 1202 call write_timing_log(timing_string,setup_time) 1203 1204 !-------------------------------------------------------------------------------- 1205 ! 1206 ! Compute the Analytic energy using Shift Lanczos 1207 ! 1208 !-------------------------------------------------------------------------------- 1209 1210 call cpu_time(time1) 1211 1212 ABI_MALLOC(AT_Lanczos,(n_ext_freq,lmax)) 1213 ! Note that the array eps^{-1}(0) - 1 is diagonal in the lanczos basis already! No need to diagonalize, so we can 1214 ! use Lbasis_lanczos directly... 1215 call compute_AT_shift_Lanczos(n_ext_freq,dtset%gw_freqsp, model_parameter, lmax, Lbasis_lanczos, kmax_analytic, AT_Lanczos) 1216 1217 call cpu_time(time2) 1218 time = time2 - time1 1219 1220 write(timing_string,'(A)') "Time to compute analytical term by SHIFT LANCZOS : " 1221 call write_timing_log(timing_string,time) 1222 1223 1224 !-------------------------------------------------------------------------------- 1225 ! 1226 ! Compute the Numeric energy using Shift Lanczos 1227 ! 1228 !-------------------------------------------------------------------------------- 1229 1230 call cpu_time(time1) 1231 1232 ABI_MALLOC(array_integrand_exact_sector,(npt_gauss+1,n_ext_freq)) 1233 ABI_MALLOC(array_integrand_model_sector,(npt_gauss+1,n_ext_freq)) 1234 1235 1236 ABI_MALLOC( tmp_dielectric_array, (lmax,lmax,npt_gauss+1)) 1237 1238 do iw = 1, npt_gauss + 1 1239 1240 lorentzian = model_parameter**2/(list_omega(iw)**2+model_parameter**2) 1241 tmp_dielectric_array(:,:,iw) = eps_m1_minus_eps_model_m1(:,:,iw)-lorentzian*eps_m1_minus_eps_model_m1(:,:,1) 1242 1243 end do 1244 1245 call compute_projected_BT_shift_Lanczos(n_ext_freq, dtset%gw_freqsp, lmax, Lbasis_lanczos, & 1246 kmax_numeric, npt_gauss, tmp_dielectric_array, array_integrand_exact_sector ) 1247 1248 array_integrand_model_sector = zero ! just a dummy array in this case 1249 1250 ABI_FREE( tmp_dielectric_array) 1251 call cpu_time(time2) 1252 time = time2 - time1 1253 write(timing_string,'(A)') "Time to compute numerical term by SHIFT LANCZOS : " 1254 call write_timing_log(timing_string,time) 1255 1256 1257 1258 !-------------------------------------------------------------------------------- 1259 ! 1260 ! set up arrays for poles 1261 ! 1262 !-------------------------------------------------------------------------------- 1263 1264 !call generate_degeneracy_table_for_poles(debug) ! so we can compute Poles contributions 1265 call generate_degeneracy_table_for_poles(.true.) ! so we can compute Poles contributions 1266 1267 !-------------------------------------------------------------------------------- 1268 ! 1269 ! Print contributions to sigma_A_Lanczos to a file 1270 ! 1271 !-------------------------------------------------------------------------------- 1272 call output_Sigma_A_by_eigenvalues(n_ext_freq,lmax,dtset%gw_freqsp,AT_Lanczos,one/epsilon_eigenvalues_0-one,1) 1273 1274 1275 !-------------------------------------------------------------------------------- 1276 ! 1277 ! Iterate on external frequencies 1278 ! 1279 !-------------------------------------------------------------------------------- 1280 1281 1282 do iw_ext = 1, dtset%gw_customnfreqsp 1283 1284 call cpu_time(freq_time1) 1285 external_omega = dtset%gw_freqsp(iw_ext) 1286 1287 write(timing_string,'(A)') "#" 1288 call write_text_block_in_Timing_log(timing_string) 1289 write(timing_string,'(A)') "#" 1290 call write_text_block_in_Timing_log(timing_string) 1291 write(timing_string,'(A,I4,A,F8.4,A)') "# Frequency # ",iw_ext," omega = ",external_omega," Ha" 1292 call write_text_block_in_Timing_log(timing_string) 1293 write(timing_string,'(A)') "#" 1294 call write_text_block_in_Timing_log(timing_string) 1295 write(timing_string,'(A)') "#" 1296 call write_text_block_in_Timing_log(timing_string) 1297 1298 1299 !-------------------------------------------------------------------------------- 1300 ! 1301 ! compute the pole term 1302 ! CAREFUL! The real valence states must still be allocated 1303 ! for the dielectric operator to work properly 1304 ! 1305 !-------------------------------------------------------------------------------- 1306 1307 call cpu_time(time1) 1308 1309 pole_energy = compute_Poles(external_omega,kmax_poles,debug) 1310 1311 call cpu_time(time2) 1312 1313 time = time2-time1 1314 write(timing_string,'(A)') "Time to compute the Poles contribution : " 1315 call write_timing_log(timing_string,time) 1316 1317 !================================================================================ 1318 ! Compute the contributions from the analytic term 1319 !================================================================================ 1320 1321 call cpu_time(time1) 1322 1323 1324 sigma_A_Lanczos = dble(sum(AT_Lanczos(iw_ext,:)*(one/epsilon_eigenvalues_0(:)-one))) 1325 1326 ABI_FREE(epsilon_eigenvalues_0) 1327 1328 call cpu_time(time2) 1329 time = time2-time1 1330 write(timing_string,'(A)') "Time for Tr[ (eps_model^{-1}-1) . AT ] AFTER SHIFT : " 1331 call write_timing_log(timing_string,time) 1332 1333 !-------------------------------------------------------------------------------- 1334 ! 1335 ! compute integrand 1336 ! 1337 !-------------------------------------------------------------------------------- 1338 1339 1340 call compute_integrands_shift_lanczos(iw_ext, n_ext_freq, npt_gauss, array_integrand_exact_sector, & 1341 array_integrand_model_sector, sigma_B_Lanczos, sigma_B_model_Lanczos) 1342 1343 1344 1345 !-------------------------------------------------------------------------------- 1346 ! 1347 ! Output results 1348 ! 1349 !-------------------------------------------------------------------------------- 1350 1351 ! just dummy variables so we can use the output_results routine 1352 sigma_A_model_Lanczos = zero 1353 sigma_B_model_Lanczos = zero 1354 call output_results(iw_ext,npt_gauss, lmax, lmax_model, model_parameter, zero, & 1355 external_omega, Sigma_x,Vxc_energy,pole_energy, & 1356 sigma_A_Lanczos,sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos, & 1357 Sigma_x_Lanczos_projected ) 1358 1359 1360 call cpu_time(freq_time2) 1361 freq_time = freq_time2-freq_time1 1362 1363 write(timing_string,'(A)') " TOTAL FREQUENCY TIME : " 1364 call write_timing_log(timing_string,freq_time) 1365 1366 1367 1368 correlations = pole_energy+sigma_A_Lanczos+sigma_B_Lanczos 1369 1370 renormalized_energy = eig(e) + Sigma_x-Vxc_energy +correlations 1371 write(std_out,10) ' ' 1372 write(std_out,14) ' For omega : ',external_omega ,' Ha = ',external_omega *Ha_eV,' eV' 1373 write(std_out,14) ' <psi_e | Sigma_c | psi_e> : ',correlations ,' Ha = ',correlations *Ha_eV,' eV' 1374 write(std_out,14) ' eps_e + <Sigma_xc - V_xc> : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV' 1375 1376 write(ab_out,10) ' ' 1377 write(ab_out,14) ' For omega : ',external_omega ,' Ha = ',external_omega *Ha_eV,' eV' 1378 write(ab_out,14) ' <psi_e | Sigma_c | psi_e>: ',correlations ,' Ha = ',correlations *Ha_eV,' eV' 1379 write(ab_out,14) ' eps_e + <Sigma_xc - V_xc> : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV' 1380 1381 1382 1383 end do 1384 1385 ABI_FREE(AT_Lanczos) 1386 ABI_FREE(array_integrand_exact_sector) 1387 ABI_FREE(array_integrand_model_sector) 1388 call clean_degeneracy_table_for_poles() 1389 call cleanup_Pk_model() 1390 call cleanup_Lanczos_basis() 1391 call cleanup_projected_Sternheimer_epsilon() 1392 1393 call cpu_time(total_time2) 1394 total_time1 = zero 1395 total_time = total_time2-total_time1 1396 write(timing_string,'(A)') " TOTAL TIME : " 1397 call write_timing_log(timing_string,total_time) 1398 1399 1400 1401 10 format(A) 1402 14 format(A,ES24.16,A,F16.8,A) 1403 1404 end subroutine compute_correlations_no_model_shift_lanczos
m_gwls_ComputeCorrelationEnergy/compute_correlations_shift_lanczos [ Functions ]
[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]
NAME
compute_correlations_shift_lanczos
FUNCTION
.
INPUTS
OUTPUT
SOURCE
81 subroutine compute_correlations_shift_lanczos(dtset, Sigma_x,Vxc_energy,debug) 82 !---------------------------------------------------------------------------------------------------- 83 ! 84 ! 85 ! This function computes the correlation energy Sigma_c(w) for the state we wish to correct. 86 ! The shift lanczos algorithm is used. 87 ! 88 !---------------------------------------------------------------------------------------------------- 89 90 type(dataset_type),intent(in) :: dtset 91 92 real(dp),intent(in) :: Sigma_x,Vxc_energy 93 logical, intent(in) :: debug 94 95 96 !Local variables 97 98 integer :: n_ext_freq 99 100 integer :: npt_gauss 101 integer :: print_debug 102 103 104 integer :: kmax_poles 105 integer :: kmax_model 106 integer :: lmax_model 107 108 109 integer :: kmax_analytic 110 integer :: kmax_numeric 111 112 real(dp) :: omega_static 113 114 real(dp) :: lorentzian 115 integer :: iw_ext 116 integer :: iw 117 118 real(dp) , allocatable :: psie_k(:,:) 119 120 121 real(dp), allocatable :: epsilon_eigenvalues_0(:) 122 real(dp), allocatable :: epsilon_model_eigenvalues_0(:) 123 124 125 complex(dpc), allocatable :: AT_Lanczos(:,:) 126 complex(dpc), allocatable :: AT_model_Lanczos(:,:) 127 128 129 ! To use lanczos instead of sqmr 130 complex(dpc), allocatable :: Lbasis_diagonalize_dielectric_terms(:,:) 131 complex(dpc), allocatable :: hermitian_static_eps_m1_minus_eps_model_m1(:,:) 132 real(dp), allocatable :: eigenvalues_static_eps_m1_minus_eps_model_m1(:) 133 real(dp), allocatable :: eigenvalues_static_eps_model_m1_minus_one(:) 134 complex(dpc), allocatable :: work(:) 135 real(dp), allocatable :: rwork(:) 136 integer :: lwork 137 integer :: info 138 integer :: l 139 140 141 integer :: debug_unit 142 character(50) :: debug_filename 143 144 real(dp) :: time1, time2, time 145 real(dp) :: total_time1,total_time2,total_time 146 real(dp) :: setup_time1, setup_time2, setup_time 147 real(dp) :: freq_time1, freq_time2, freq_time 148 149 integer :: nfrequencies 150 real(dp),allocatable :: list_projection_frequencies(:) 151 152 complex(dpc), allocatable :: array_integrand_exact_sector(:,:) 153 complex(dpc), allocatable :: array_integrand_model_sector(:,:) 154 155 156 complex(dpc), allocatable :: tmp_dielectric_array(:,:,:) 157 158 real(dp) :: external_omega 159 160 character(256) :: timing_string 161 162 integer :: recy_line_size 163 character(128) :: recy_name 164 logical :: local_tmp_exist 165 logical :: use_model 166 167 ! Energy contributions 168 169 real(dp) :: pole_energy 170 171 real(dp) :: sigma_A_Lanczos 172 real(dp) :: sigma_A_model_Lanczos 173 real(dp) :: sigma_B_Lanczos 174 real(dp) :: sigma_B_model_Lanczos 175 176 real(dp) :: correlations 177 real(dp) :: renormalized_energy 178 179 real(dp) :: second_model_parameter 180 181 real(dp):: tsec(2) 182 integer :: GWLS_TIMAB, OPTION_TIMAB 183 character(500) :: msg 184 185 ! ************************************************************************* 186 187 !-------------------------------------------------------------------------------- 188 ! 189 ! Set up variables and allocate arrays 190 ! 191 !-------------------------------------------------------------------------------- 192 193 call cpu_time(total_time1) 194 !Variable allocation and initialization 195 model_number = dtset%gwls_diel_model 196 model_parameter = dtset%gwls_model_parameter 197 npt_gauss = dtset%gwls_npt_gauss_quad 198 print_debug = dtset%gwls_print_debug 199 200 first_seed = dtset%gwls_first_seed 201 e = dtset%gwls_band_index 202 203 204 !second_model_parameter = dtset%gwls_second_model_parameter 205 second_model_parameter = zero 206 207 208 209 ! set variables from gwls_GenerateEpsilon module 210 kmax = dtset%gwls_stern_kmax 211 nseeds = dtset%gwls_nseeds 212 213 kmax_model = dtset%gwls_kmax_complement 214 kmax_poles = dtset%gwls_kmax_poles 215 kmax_analytic = dtset%gwls_kmax_analytic 216 kmax_numeric = dtset%gwls_kmax_numeric 217 218 n_ext_freq = dtset%gw_customnfreqsp 219 220 use_model = .True. 221 222 call cpu_time(setup_time1) 223 224 !-------------------------------------------------------------------------------- 225 ! 226 ! Extract the frequencies at which the integrand will be evaluated 227 ! add the value zero in the set. 228 ! 229 !-------------------------------------------------------------------------------- 230 231 call generate_frequencies_and_weights(npt_gauss) 232 233 !-------------------------------------------------------------------------------- 234 ! 235 ! 236 ! Compute the static bases for the exact and model 237 ! dielectric operator 238 ! 239 ! 240 !-------------------------------------------------------------------------------- 241 242 omega_static = zero 243 ! define dimensions 244 lmax = nseeds*kmax 245 lmax_model = nseeds*kmax_model 246 247 ! Allocate arrays which will contain basis 248 call setup_Lanczos_basis(lmax,lmax_model) 249 250 ! allocate eigenvalues array 251 ABI_MALLOC(epsilon_eigenvalues_0, (lmax)) 252 ABI_MALLOC(epsilon_model_eigenvalues_0, (lmax_model)) 253 254 ! set omega=0 for exact dielectric operator 255 call set_dielectric_function_frequency([0.0_dp,omega_static]) 256 257 ! and make note that the Sternheimer solutions must be kept (for use in the projected Sternheimer section). 258 if(dtset%gwls_recycle == 1) then 259 ABI_MALLOC(Sternheimer_solutions_zero,(2,npw_k,lmax,nbandv)) 260 Sternheimer_solutions_zero = zero 261 write_solution = .true. 262 end if 263 if(dtset%gwls_recycle == 2) then 264 write(recy_name,'(A,I0.4,A)') "Recycling_",mpi_enreg%me,".dat" 265 266 inquire(iolength=recy_line_size) cg(:,1:npw_k) 267 268 inquire(file='local_tmp', exist=local_tmp_exist) 269 if(local_tmp_exist) recy_name = 'local_tmp/' // recy_name(1:118) 270 271 if (open_file(file=recy_name,iomsg=msg,newunit=recy_unit,access='direct',form='unformatted',& 272 & status='replace',recl=recy_line_size)/=0) then 273 ABI_ERROR(msg) 274 end if 275 276 write_solution = .true. 277 end if 278 279 call cpu_time(time1) 280 ! Compute the Lanczos basis using Block Lanczos; store 281 ! basis in Lbasis_lanczos 282 GWLS_TIMAB = 1504 283 OPTION_TIMAB = 1 284 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 285 286 call driver_generate_dielectric_matrix( matrix_function_epsilon_k, & 287 nseeds, kmax, & 288 epsilon_eigenvalues_0, & 289 Lbasis_lanczos, debug) 290 OPTION_TIMAB = 2 291 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 292 293 call cpu_time(time2) 294 time = time2-time1 295 296 write(timing_string,'(A)') "Time to compute the EXACT Static Dielectric Matrix : " 297 call write_timing_log(timing_string,time) 298 299 300 call output_epsilon_eigenvalues(lmax,epsilon_eigenvalues_0,1) 301 302 303 ! The Sternheimer solutions at $\omega = 0$ have been stored. 304 write_solution = .false. 305 306 307 ! Prepare the model dielectric operator 308 call setup_Pk_model(omega_static,second_model_parameter) 309 310 call cpu_time(time1) 311 ! Compute the Lanczos basis of the model operator using Block Lanczos; store 312 ! basis in Lbasis_model_lanczos 313 GWLS_TIMAB = 1505 314 OPTION_TIMAB = 1 315 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 316 317 call driver_generate_dielectric_matrix(matrix_function_epsilon_model_operator, & 318 nseeds, kmax_model, & 319 epsilon_model_eigenvalues_0, & 320 Lbasis_model_lanczos, debug) 321 OPTION_TIMAB = 2 322 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 323 324 call cpu_time(time2) 325 time = time2-time1 326 327 write(timing_string,'(A)') "Time to compute the MODEL Static Dielectric Matrix : " 328 call write_timing_log(timing_string,time) 329 330 331 call output_epsilon_eigenvalues(lmax_model,epsilon_model_eigenvalues_0,2) 332 333 334 ABI_MALLOC(eigenvalues_static_eps_model_m1_minus_one, (lmax_model)) 335 336 do l = 1, lmax_model 337 eigenvalues_static_eps_model_m1_minus_one(l) = one/epsilon_model_eigenvalues_0(l)-one 338 end do 339 340 341 342 !-------------------------------------------------------------------------------- 343 ! 344 ! 345 ! Prepare and compute the projection of the dielectric Sternheimer equations 346 ! 347 ! 348 !-------------------------------------------------------------------------------- 349 350 ! Setup various arrays necessary for the Sternheimer projection scheme 351 nfrequencies = dtset%gwls_n_proj_freq 352 ABI_MALLOC(list_projection_frequencies,(nfrequencies)) 353 354 list_projection_frequencies = dtset%gwls_list_proj_freq 355 356 call cpu_time(time1) 357 GWLS_TIMAB = 1507 358 OPTION_TIMAB = 1 359 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 360 361 !call setup_projected_Sternheimer_epsilon(lmax, npt_gauss, second_model_parameter, & 362 ! list_projection_frequencies,nfrequencies,debug) 363 364 365 call ProjectedSternheimerEpsilon(lmax, npt_gauss, second_model_parameter, & 366 list_projection_frequencies,nfrequencies,& 367 epsilon_eigenvalues_0,debug,use_model) 368 369 370 371 372 ! The Sternheimer solutions at $\omega = 0$ have been used to make the basis for the projected Sternheimer equations. 373 if(dtset%gwls_recycle == 1) then 374 ABI_FREE(Sternheimer_solutions_zero) 375 end if 376 if(dtset%gwls_recycle == 2) then 377 close(recy_unit,status='delete') 378 end if 379 380 OPTION_TIMAB = 2 381 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 382 383 384 call cpu_time(time2) 385 time = time2-time1 386 write(timing_string,'(A)') "Time to setup and compute the projected Sternheimer epsilon : " 387 call write_timing_log(timing_string,time) 388 389 390 call cpu_time(time1) 391 392 GWLS_TIMAB = 1508 393 OPTION_TIMAB = 1 394 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 395 396 call compute_eps_m1_minus_eps_model_m1(lmax, npt_gauss) 397 398 OPTION_TIMAB = 2 399 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 400 401 call cpu_time(time2) 402 time = time2-time1 403 write(timing_string,'(A)') "Time to compute eps^{-1}-eps_model^{-1} : " 404 call write_timing_log(timing_string,time) 405 406 407 !-------------------------------------------------------------------------------- 408 ! 409 ! 410 ! Compute the model dielectric array 411 ! 412 ! 413 !-------------------------------------------------------------------------------- 414 415 call cpu_time(time1) 416 417 GWLS_TIMAB = 1509 418 OPTION_TIMAB = 1 419 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 420 421 call compute_eps_model_m1_minus_one(lmax_model, npt_gauss, second_model_parameter, epsilon_model_eigenvalues_0) 422 423 OPTION_TIMAB = 2 424 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 425 ABI_FREE(epsilon_model_eigenvalues_0) 426 427 428 call cpu_time(time2) 429 time = time2-time1 430 write(timing_string,'(A)') "Time to compute eps_model^{-1}-1 : " 431 call write_timing_log(timing_string,time) 432 433 434 !-------------------------------------------------------------------------------- 435 ! 436 ! 437 ! We no longer need the Lanczos basis in its current form! 438 ! Modify the basis so that it now contains (V^1/2.L)^*.psie 439 ! 440 ! 441 !-------------------------------------------------------------------------------- 442 443 call cpu_time(time1) 444 ABI_MALLOC(psie_k, (2,npw_k)) 445 446 GWLS_TIMAB = 1510 447 OPTION_TIMAB = 1 448 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 449 450 451 psie_k = cg(:,(e-1)*npw_k+1:e*npw_k) 452 453 call modify_Lbasis_Coulomb(psie_k, lmax, lmax_model) 454 455 ABI_FREE(psie_k) 456 457 OPTION_TIMAB = 2 458 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 459 460 461 call cpu_time(time2) 462 time = time2-time1 463 write(timing_string,'(A)') "Time to modify the Lanczos basis : " 464 call write_timing_log(timing_string,time) 465 466 !-------------------------------------------------------------------------------- 467 ! 468 ! 469 ! 470 ! Diagonalize the static array eps^{-1} - eps^{-1}_model in order to 471 ! be able to apply diagonal shift Lanczos. 472 ! 473 !-------------------------------------------------------------------------------- 474 475 call cpu_time(time1) 476 ! diagonalize the static eps^{-1} - eps^{-1}_model array, so as the use the diagonal Lanczos procedure 477 ! for the analytical term 478 GWLS_TIMAB = 1511 479 OPTION_TIMAB = 1 480 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 481 482 ABI_MALLOC(Lbasis_diagonalize_dielectric_terms, (npw_k,lmax)) 483 ABI_MALLOC(hermitian_static_eps_m1_minus_eps_model_m1, (lmax,lmax)) 484 485 ABI_MALLOC(eigenvalues_static_eps_m1_minus_eps_model_m1, (lmax)) 486 487 ABI_MALLOC(rwork, (3*lmax-2)) 488 489 490 hermitian_static_eps_m1_minus_eps_model_m1(:,:) = eps_m1_minus_eps_model_m1(:,:,1) 491 492 493 494 ! WORK QUERRY 495 lwork = -1 496 ABI_MALLOC(work, (1)) 497 call ZHEEV( 'V', & ! Compute eigenvectors and eigenvalues 498 'U', & ! use Upper triangular part 499 lmax, & ! order of matrix 500 hermitian_static_eps_m1_minus_eps_model_m1, & ! initial matrix on input; eigenvectors on output 501 lmax, & ! LDA 502 eigenvalues_static_eps_m1_minus_eps_model_m1,& ! eigenvalues 503 work, lwork, rwork, info ) ! work stuff 504 505 if ( info /= 0) then 506 debug_unit = get_unit() 507 write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log' 508 509 open(debug_unit,file=trim(debug_filename),status='unknown') 510 511 write(debug_unit,'(A)') '*********************************************************************************************' 512 write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZHEEV (1), gwls_ComputeCorrelationEnergy' 513 write(debug_unit,'(A)') '*********************************************************************************************' 514 515 close(debug_unit) 516 517 end if 518 519 ! COMPUTATION 520 lwork = nint(dble(work(1))) 521 ABI_FREE(work) 522 ABI_MALLOC(work, (lwork)) 523 call ZHEEV( 'V', & ! Compute eigenvectors and eigenvalues 524 'U', & ! use Upper triangular part 525 lmax, & ! order of matrix 526 hermitian_static_eps_m1_minus_eps_model_m1, & ! initial matrix on input; eigenvectors on output 527 lmax, & ! LDA 528 eigenvalues_static_eps_m1_minus_eps_model_m1,& ! eigenvalues 529 work, lwork, rwork, info ) ! work stuff 530 531 if ( info /= 0) then 532 debug_unit = get_unit() 533 write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log' 534 535 open(debug_unit,file=trim(debug_filename),status='unknown') 536 537 write(debug_unit,'(A)') '*********************************************************************************************' 538 write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZHEEV (2), gwls_ComputeCorrelationEnergy' 539 write(debug_unit,'(A)') '*********************************************************************************************' 540 541 close(debug_unit) 542 543 end if 544 545 546 547 ABI_FREE(work) 548 ABI_FREE(rwork) 549 550 !-------------------------------------------------------------------------------- 551 ! 552 ! update basis: L' = L . Q 553 ! 554 ! 555 ! CAREFUL!!! We must multiply by conjg(hermitian_static_eps_m1_minus_eps_model_m1), 556 ! which is the COMPLEX CONJUGATE of the eigenvectors of the matrix 557 ! eps^{-1}-eps_m^{-1}, because we have MODIFIED the basis Lbasis_lanczos 558 ! to contain the complex conjugate of the eigenvectors of eps. 559 ! This is somewhat subtle, but forgetting to do this leads to small errors 560 ! in the results... 561 !-------------------------------------------------------------------------------- 562 hermitian_static_eps_m1_minus_eps_model_m1 = conjg(hermitian_static_eps_m1_minus_eps_model_m1) 563 564 call ZGEMM('N','N',npw_k,lmax,lmax,cmplx_1,Lbasis_lanczos,npw_k, & 565 hermitian_static_eps_m1_minus_eps_model_m1, & 566 lmax,cmplx_0,Lbasis_diagonalize_dielectric_terms,npw_k) 567 568 569 OPTION_TIMAB = 2 570 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 571 572 573 call cpu_time(time2) 574 time = time2-time1 575 write(timing_string,'(A)') "Time to diagonalize eps^{-1}(0)-eps^{-1}(0)_model : " 576 call write_timing_log(timing_string,time) 577 578 ABI_FREE(hermitian_static_eps_m1_minus_eps_model_m1) 579 580 581 call cpu_time(setup_time2) 582 setup_time = setup_time2 - setup_time1 583 584 write(timing_string,'(A)') " TOTAL DIELECTRIC SETUP TIME : " 585 call write_timing_log(timing_string,setup_time) 586 587 !-------------------------------------------------------------------------------- 588 ! 589 ! Compute the Analytic energy using Shift Lanczos 590 ! 591 !-------------------------------------------------------------------------------- 592 593 call cpu_time(time1) 594 595 GWLS_TIMAB = 1512 596 OPTION_TIMAB = 1 597 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 598 599 600 ABI_MALLOC(AT_Lanczos,(n_ext_freq,lmax)) 601 call compute_AT_shift_Lanczos(n_ext_freq,dtset%gw_freqsp, model_parameter, lmax, Lbasis_diagonalize_dielectric_terms,& 602 & kmax_analytic, AT_Lanczos) 603 604 OPTION_TIMAB = 2 605 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 606 607 608 ABI_FREE(Lbasis_diagonalize_dielectric_terms) 609 610 call cpu_time(time2) 611 time = time2 - time1 612 613 write(timing_string,'(A)') "Time to compute analytical term by SHIFT LANCZOS : " 614 call write_timing_log(timing_string,time) 615 616 617 call cpu_time(time1) 618 619 GWLS_TIMAB = 1513 620 OPTION_TIMAB = 1 621 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 622 623 ABI_MALLOC(AT_model_Lanczos,(n_ext_freq,lmax_model)) 624 call compute_AT_shift_Lanczos(n_ext_freq,dtset%gw_freqsp, model_parameter, lmax_model, & 625 Lbasis_model_lanczos, kmax_analytic, AT_model_Lanczos) 626 627 OPTION_TIMAB = 2 628 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 629 630 call cpu_time(time2) 631 time = time2 - time1 632 633 write(timing_string,'(A)') "Time to compute analytical MODEL by SHIFT LANCZOS : " 634 call write_timing_log(timing_string,time) 635 636 637 !-------------------------------------------------------------------------------- 638 ! 639 ! Compute the Numeric energy using Shift Lanczos 640 ! 641 !-------------------------------------------------------------------------------- 642 643 call cpu_time(time1) 644 645 ABI_MALLOC(array_integrand_exact_sector,(npt_gauss+1,n_ext_freq)) 646 647 ABI_MALLOC( tmp_dielectric_array, (lmax,lmax,npt_gauss+1)) 648 649 do iw = 1, npt_gauss + 1 650 651 lorentzian = model_parameter**2/(list_omega(iw)**2+model_parameter**2) 652 tmp_dielectric_array(:,:,iw) = eps_m1_minus_eps_model_m1(:,:,iw)-lorentzian*eps_m1_minus_eps_model_m1(:,:,1) 653 654 end do 655 GWLS_TIMAB = 1514 656 OPTION_TIMAB = 1 657 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 658 659 call compute_projected_BT_shift_Lanczos(n_ext_freq, dtset%gw_freqsp, lmax, Lbasis_lanczos, & 660 kmax_numeric, npt_gauss, tmp_dielectric_array, array_integrand_exact_sector ) 661 662 OPTION_TIMAB = 2 663 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 664 665 ABI_FREE( tmp_dielectric_array) 666 call cpu_time(time2) 667 time = time2 - time1 668 write(timing_string,'(A)') "Time to compute numerical term by SHIFT LANCZOS : " 669 call write_timing_log(timing_string,time) 670 671 672 673 call cpu_time(time1) 674 675 ABI_MALLOC(array_integrand_model_sector,(npt_gauss+1,n_ext_freq)) 676 677 !ABI_MALLOC( tmp_dielectric_array, (lmax_model,lmax_model,npt_gauss+1)) 678 ABI_MALLOC( tmp_dielectric_array, (lmax_model,blocksize_epsilon,npt_gauss+1)) 679 680 681 do iw = 1, npt_gauss + 1 682 683 lorentzian = model_parameter**2/(list_omega(iw)**2+model_parameter**2) 684 !tmp_dielectric_array(:,:,iw) = eps_model_m1_minus_one(:,:,iw)-lorentzian*eps_model_m1_minus_one(:,:,1) 685 tmp_dielectric_array(:,:,iw) = eps_model_m1_minus_one_DISTR(:,:,iw)-lorentzian*eps_model_m1_minus_one_DISTR(:,:,1) 686 687 end do 688 689 GWLS_TIMAB = 1515 690 OPTION_TIMAB = 1 691 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 692 693 !call compute_projected_BT_shift_Lanczos(n_ext_freq , dtset%gw_freqsp, lmax_model, Lbasis_model_lanczos, & 694 ! kmax_numeric, npt_gauss, tmp_dielectric_array, array_integrand_model_sector ) 695 696 697 call compute_projected_BT_shift_Lanczos_DISTRIBUTED(n_ext_freq, dtset%gw_freqsp, lmax_model, blocksize_epsilon, & 698 model_lanczos_vector_belongs_to_this_node, model_lanczos_vector_index, & 699 Lbasis_model_lanczos, kmax_numeric, npt_gauss, tmp_dielectric_array, & 700 array_integrand_model_sector ) 701 702 OPTION_TIMAB = 2 703 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 704 705 ABI_FREE(tmp_dielectric_array) 706 call cpu_time(time2) 707 time = time2 - time1 708 write(timing_string,'(A)') "Time to compute numerical model SHIFT LANCZOS : " 709 call write_timing_log(timing_string,time) 710 711 712 !-------------------------------------------------------------------------------- 713 ! 714 ! set up arrays for poles 715 ! 716 !-------------------------------------------------------------------------------- 717 718 719 720 !call generate_degeneracy_table_for_poles(debug) ! so we can compute Poles contributions 721 call generate_degeneracy_table_for_poles(.true.) ! so we can compute Poles contributions 722 723 !-------------------------------------------------------------------------------- 724 ! 725 ! Print contributions to sigma_A_Lanczos to a file 726 ! 727 !-------------------------------------------------------------------------------- 728 call output_Sigma_A_by_eigenvalues(n_ext_freq,lmax,dtset%gw_freqsp,AT_Lanczos,eigenvalues_static_eps_m1_minus_eps_model_m1,2) 729 call output_Sigma_A_by_eigenvalues(n_ext_freq,lmax_model,dtset%gw_freqsp,AT_model_Lanczos,& 730 & eigenvalues_static_eps_model_m1_minus_one,3) 731 732 !epsilon_eigenvalues_0 733 734 ABI_FREE(epsilon_eigenvalues_0) 735 !-------------------------------------------------------------------------------- 736 ! 737 ! Iterate on external frequencies 738 ! 739 !-------------------------------------------------------------------------------- 740 741 742 do iw_ext = 1, dtset%gw_customnfreqsp 743 744 call cpu_time(freq_time1) 745 external_omega = dtset%gw_freqsp(iw_ext) 746 747 write(timing_string,'(A)') "#" 748 call write_text_block_in_Timing_log(timing_string) 749 write(timing_string,'(A)') "#" 750 call write_text_block_in_Timing_log(timing_string) 751 write(timing_string,'(A,I4,A,F8.4,A)') "# Frequency # ",iw_ext," omega = ",external_omega," Ha" 752 call write_text_block_in_Timing_log(timing_string) 753 write(timing_string,'(A)') "#" 754 call write_text_block_in_Timing_log(timing_string) 755 write(timing_string,'(A)') "#" 756 call write_text_block_in_Timing_log(timing_string) 757 758 759 !-------------------------------------------------------------------------------- 760 ! 761 ! compute the pole term 762 ! CAREFUL! The real valence states must still be allocated 763 ! for the dielectric operator to work properly 764 ! 765 !-------------------------------------------------------------------------------- 766 767 call cpu_time(time1) 768 GWLS_TIMAB = 1516 769 OPTION_TIMAB = 1 770 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 771 772 pole_energy = compute_Poles(external_omega,kmax_poles,debug) 773 774 OPTION_TIMAB = 2 775 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 776 777 call cpu_time(time2) 778 779 time = time2-time1 780 write(timing_string,'(A)') "Time to compute the Poles contribution : " 781 call write_timing_log(timing_string,time) 782 783 784 785 !================================================================================ 786 ! Compute the contributions from the analytic term 787 !================================================================================ 788 GWLS_TIMAB = 1517 789 OPTION_TIMAB = 1 790 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 791 792 793 call cpu_time(time1) 794 795 sigma_A_Lanczos = dble(sum(AT_Lanczos(iw_ext,:)*eigenvalues_static_eps_m1_minus_eps_model_m1(:))) 796 797 call cpu_time(time2) 798 time = time2-time1 799 write(timing_string,'(A)') "Time Tr[(eps^{-1}-eps_model^{-1}).AT] AFTER SHIFT : " 800 call write_timing_log(timing_string,time) 801 802 803 call cpu_time(time1) 804 805 sigma_A_model_Lanczos= dble(sum(AT_model_Lanczos(iw_ext,:)*eigenvalues_static_eps_model_m1_minus_one(:))) 806 807 call cpu_time(time2) 808 time = time2-time1 809 write(timing_string,'(A)') "Time for Tr[ (eps_model^{-1}-1) . AT ] AFTER SHIFT : " 810 call write_timing_log(timing_string,time) 811 812 OPTION_TIMAB = 2 813 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 814 !-------------------------------------------------------------------------------- 815 ! 816 ! compute integrand 817 ! 818 !-------------------------------------------------------------------------------- 819 GWLS_TIMAB = 1518 820 OPTION_TIMAB = 1 821 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 822 823 call compute_integrands_shift_lanczos(iw_ext, n_ext_freq, npt_gauss, array_integrand_exact_sector, & 824 array_integrand_model_sector, sigma_B_Lanczos, sigma_B_model_Lanczos) 825 826 OPTION_TIMAB = 2 827 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec) 828 829 !-------------------------------------------------------------------------------- 830 ! 831 ! Output results 832 ! 833 !-------------------------------------------------------------------------------- 834 835 836 call output_results(iw_ext,npt_gauss, lmax,lmax_model, model_parameter, second_model_parameter, & 837 external_omega, Sigma_x,Vxc_energy,pole_energy, & 838 sigma_A_Lanczos,sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos) 839 840 841 call cpu_time(freq_time2) 842 freq_time = freq_time2-freq_time1 843 844 write(timing_string,'(A)') " TOTAL FREQUENCY TIME : " 845 call write_timing_log(timing_string,freq_time) 846 847 848 849 correlations = pole_energy+sigma_A_Lanczos+sigma_A_model_Lanczos+sigma_B_Lanczos+sigma_B_model_Lanczos 850 851 renormalized_energy = eig(e) + Sigma_x-Vxc_energy +correlations 852 write(std_out,10) ' ' 853 write(std_out,14) ' For omega : ',external_omega ,' Ha = ',external_omega *Ha_eV,' eV' 854 write(std_out,14) ' <psi_e | Sigma_c | psi_e>: ',correlations ,' Ha = ',correlations *Ha_eV,' eV' 855 write(std_out,14) ' eps_e + <Sigma_xc - V_xc> : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV' 856 857 write(ab_out,10) ' ' 858 write(ab_out,14) ' For omega : ',external_omega ,' Ha = ',external_omega *Ha_eV,' eV' 859 write(ab_out,14) ' <psi_e | Sigma_c | psi_e>: ',correlations ,' Ha = ',correlations *Ha_eV,' eV' 860 write(ab_out,14) ' eps_e + <Sigma_xc - V_xc> : ',renormalized_energy,' Ha = ',renormalized_energy*Ha_eV,' eV' 861 862 863 end do 864 865 ABI_FREE(AT_Lanczos) 866 ABI_FREE(AT_model_Lanczos) 867 ABI_FREE(array_integrand_exact_sector) 868 ABI_FREE(array_integrand_model_sector) 869 ABI_FREE(eigenvalues_static_eps_model_m1_minus_one) 870 ABI_FREE(eigenvalues_static_eps_m1_minus_eps_model_m1) 871 call clean_degeneracy_table_for_poles() 872 call cleanup_Pk_model() 873 call cleanup_Lanczos_basis() 874 call cleanup_projected_Sternheimer_epsilon() 875 876 call cpu_time(total_time2) 877 total_time = total_time2-total_time1 878 write(timing_string,'(A)') " TOTAL TIME : " 879 call write_timing_log(timing_string,total_time) 880 881 882 10 format(A) 883 14 format(A,ES24.16,A,F16.8,A) 884 885 end subroutine compute_correlations_shift_lanczos
m_gwls_ComputeCorrelationEnergy/compute_integrands_shift_lanczos [ Functions ]
[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]
NAME
compute_integrands_shift_lanczos
FUNCTION
.
INPUTS
OUTPUT
SOURCE
1420 subroutine compute_integrands_shift_lanczos(iw_ext,n_ext_freq,npt_gauss, array_integrand_exact_sector, & 1421 array_integrand_model_sector, sigma_B_Lanczos, sigma_B_model_Lanczos) 1422 !---------------------------------------------------------------------------------------------------- 1423 ! 1424 ! This subroutine computes the integrands, assuming data was generated by shift lanczos. 1425 !---------------------------------------------------------------------------------------------------- 1426 1427 integer, intent(in) :: iw_ext, npt_gauss, n_ext_freq 1428 1429 complex(dpc), intent(in) :: array_integrand_exact_sector(npt_gauss+1,n_ext_freq) 1430 complex(dpc), intent(in) :: array_integrand_model_sector(npt_gauss+1,n_ext_freq) 1431 1432 real(dp), intent(out) :: sigma_B_Lanczos, sigma_B_model_Lanczos 1433 1434 real(dp) :: integrand_Lanczos , integrand_model_Lanczos 1435 real(dp) :: omega_prime 1436 1437 integer :: iw 1438 integer :: io_unit 1439 character(256) :: title_string 1440 character(256) :: timing_string 1441 1442 1443 real(dp) :: time1, time2, time 1444 1445 ! ************************************************************************* 1446 1447 1448 if (mpi_enreg%me == 0) then 1449 io_unit = get_unit() 1450 1451 if (iw_ext < 10) then 1452 write(title_string,'(A,I1,A)') 'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat' 1453 else if (iw_ext < 100) then 1454 write(title_string,'(A,I2,A)') 'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat' 1455 else 1456 write(title_string,'(A,I3,A)') 'APPROXIMATE_INTEGRANDS_',iw_ext,'.dat' 1457 end if 1458 1459 open(file=title_string,status=files_status_new,unit=io_unit) 1460 write(io_unit,10) '#-----------------------------------------------------------------------------------------------' 1461 write(io_unit,10) '#' 1462 write(io_unit,10) '# Approximate Integrands as a function of frequency' 1463 write(io_unit,10) '#' 1464 write(io_unit,10) '# I1 = Tr[ (eps^{-1}(iw) - eps_model^{-1}(iw)) - ' 1465 write(io_unit,10) '# f(w)(eps^{-1}(0) - eps_model^{-1}(0)) BT(w) ]' 1466 write(io_unit,10) '#' 1467 write(io_unit,10) '# I2 = Tr[ (eps_model^{-1}(iw) - 1) - f(w)(eps_model^{-1}(0)-1) BT(w)]' 1468 write(io_unit,10) '# ' 1469 write(io_unit,10) '# DIAG[I] will represent the contribution coming from taking only the diagonal elements of' 1470 write(io_unit,10) '# the arrays in the trace.' 1471 write(io_unit,10) '#' 1472 write(io_unit,10) '# omega (Ha) I1 I2 gaussian weight' 1473 write(io_unit,10) '#-----------------------------------------------------------------------------------------------' 1474 flush(io_unit) 1475 end if 1476 1477 call cpu_time(time1) 1478 1479 sigma_B_Lanczos = zero 1480 sigma_B_model_Lanczos = zero 1481 1482 do iw = 1,npt_gauss+1 1483 1484 omega_prime = list_omega(iw) 1485 1486 integrand_Lanczos = dble(array_integrand_exact_sector(iw,iw_ext)) 1487 integrand_model_Lanczos = dble(array_integrand_model_sector(iw,iw_ext)) 1488 1489 sigma_B_Lanczos = sigma_B_Lanczos + integrand_Lanczos*list_weights(iw) 1490 sigma_B_model_Lanczos = sigma_B_model_Lanczos + integrand_model_Lanczos*list_weights(iw) 1491 1492 1493 1494 if (mpi_enreg%me == 0) write(io_unit,8) omega_prime, integrand_Lanczos , integrand_model_Lanczos, list_weights(iw) 1495 1496 1497 1498 end do 1499 call cpu_time(time2) 1500 1501 if (mpi_enreg%me == 0) then 1502 write(io_unit,10) '' 1503 write(io_unit,14) '# Value of the I1 integral: ',sigma_B_Lanczos ,' Ha' 1504 write(io_unit,14) '# Value of the I2 integral: ',sigma_B_model_Lanczos ,' Ha' 1505 write(io_unit,10) '' 1506 write(io_unit,10) '' 1507 close(io_unit) 1508 end if 1509 1510 time = time2-time1 1511 write(timing_string,'(A)') "Time compute the integrands and integrals : " 1512 call write_timing_log(timing_string,time) 1513 1514 8 format(4ES24.16) 1515 10 format(A) 1516 14 format(A,ES24.16,A) 1517 1518 1519 end subroutine compute_integrands_shift_lanczos
m_gwls_ComputeCorrelationEnergy/output_epsilon_eigenvalues [ Functions ]
[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]
NAME
output_epsilon_eigenvalues
FUNCTION
.
INPUTS
OUTPUT
SOURCE
1660 subroutine output_epsilon_eigenvalues(lmax,eigenvalues,which_case) 1661 !---------------------------------------------------------------------------------------------------- 1662 ! This routine outputs the eigenvalues of the static dielectric matrix 1663 ! 1664 ! There are two cases to consider: 1665 ! 1 ) the exact dielectric matrix 1666 ! 2 ) the model dielectric matrix 1667 !---------------------------------------------------------------------------------------------------- 1668 1669 integer, intent(in) :: lmax, which_case 1670 real(dp), intent(in) :: eigenvalues(lmax) 1671 1672 1673 integer :: io_unit 1674 character(128) :: filename 1675 1676 integer :: l 1677 1678 ! ************************************************************************* 1679 1680 if (mpi_enreg%me == 0) then 1681 io_unit = get_unit() 1682 1683 if (which_case == 1) then 1684 write(filename,'(A)') "EPSILON_EIGENVALUES.dat" 1685 else if (which_case == 2) then 1686 write(filename,'(A)') "MODEL_EPSILON_EIGENVALUES.dat" 1687 end if 1688 1689 1690 open(file=filename,status=files_status_new,unit=io_unit) 1691 write(io_unit,10) '#-----------------------------------------------------------------------------------------------' 1692 write(io_unit,10) '#' 1693 write(io_unit,10) '# This file contains the computed eigenvalues of the static dielectric operator' 1694 write(io_unit,10) '# either exact or model, as indicated by the name of this file.' 1695 write(io_unit,10) '#' 1696 write(io_unit,10) '# l epsilon_l ' 1697 write(io_unit,10) '#-----------------------------------------------------------------------------------------------' 1698 1699 do l = 1, lmax 1700 1701 write(io_unit,20) l, eigenvalues(l) 1702 end do 1703 1704 1705 close(io_unit) 1706 1707 1708 end if 1709 1710 10 format(A) 1711 20 format(I7,20X,ES24.16) 1712 1713 1714 end subroutine output_epsilon_eigenvalues
m_gwls_ComputeCorrelationEnergy/output_results [ Functions ]
[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]
NAME
output_results
FUNCTION
.
INPUTS
OUTPUT
SOURCE
1535 subroutine output_results(iw_ext,npt_gauss, lmax,lmax_model, model_parameter, second_model_parameter, external_omega, & 1536 Sigma_x,Vxc_energy,pole_energy,sigma_A_Lanczos,sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos,& 1537 Sigma_x_Lanczos_projected ) 1538 !---------------------------------------------------------------------------------------------------- 1539 ! 1540 ! This subroutine computes the integrands 1541 !---------------------------------------------------------------------------------------------------- 1542 1543 integer, intent(in) :: iw_ext,lmax,lmax_model,npt_gauss 1544 real(dp), intent(in) :: model_parameter, second_model_parameter, external_omega 1545 real(dp), intent(in) :: Sigma_x,Vxc_energy,pole_energy,sigma_A_Lanczos 1546 real(dp), intent(in) :: sigma_A_model_Lanczos,sigma_B_Lanczos,sigma_B_model_Lanczos 1547 1548 real(dp), optional, intent(in) :: Sigma_x_Lanczos_projected 1549 1550 integer :: io_unit 1551 1552 character(128) :: filename 1553 1554 real(dp) :: Sigma_c 1555 1556 ! ************************************************************************* 1557 1558 1559 if (mpi_enreg%me == 0) then 1560 io_unit = get_unit() 1561 if (iw_ext < 10) then 1562 write(filename,'(A,I1,A)') 'ALL_ENERGY_',iw_ext,'.dat' 1563 else if (iw_ext < 100) then 1564 write(filename,'(A,I2,A)') 'ALL_ENERGY_',iw_ext,'.dat' 1565 else 1566 write(filename,'(A,I3,A)') 'ALL_ENERGY_',iw_ext,'.dat' 1567 end if 1568 1569 1570 open(file=filename,status=files_status_new,unit=io_unit) 1571 write(io_unit,10) '#-----------------------------------------------------------------------------------------------' 1572 write(io_unit,10) '#' 1573 write(io_unit,10) '# This file contains the results of the Correlation energy calculation.' 1574 write(io_unit,10) '# ' 1575 write(io_unit,10) '# Definitions:' 1576 write(io_unit,10) '# ' 1577 write(io_unit,10) '# eps_e = Bare DFT energy of the state ' 1578 write(io_unit,10) '# ' 1579 write(io_unit,10) '# Sigma_A_1 = Tr[(eps^{-1}(0) - eps_model^{-1}(0)) AT(W) ] ' 1580 write(io_unit,10) '# ' 1581 write(io_unit,10) '# Sigma_A_2 = Tr[(eps_model^{-1}(0) - 1) AT(W) ] ' 1582 write(io_unit,10) '# ' 1583 write(io_unit,10) '# Sigma_B_1 = Int dw I1(w) ' 1584 write(io_unit,10) '# ' 1585 write(io_unit,10) '# Sigma_B_2 = Int dw I2(w) ' 1586 write(io_unit,10) '# ' 1587 write(io_unit,10) '# I1(w) = Tr[ (eps^{-1}(iw) - eps_model^{-1}(iw)) - ' 1588 write(io_unit,10) '# f(w)(eps^{-1}(0) - eps_model^{-1}(0)) BT(w,W) ]' 1589 write(io_unit,10) '#' 1590 write(io_unit,10) '# I2(w) = Tr[ (eps_model^{-1}(iw) - 1) - f(w)(eps_model^{-1}(0)-1) BT(w,W)]' 1591 write(io_unit,10) '# ' 1592 write(io_unit,10) '# Parameters: ' 1593 write(io_unit,10) '# ' 1594 write(io_unit,25) '# lmax = ', lmax 1595 write(io_unit,25) '# lmax_model = ', lmax_model 1596 write(io_unit,25) '# npt_gauss = ', npt_gauss 1597 write(io_unit,12) '# omega0 = ', model_parameter,' Ha' 1598 write(io_unit,12) '# epsilon0 = ', second_model_parameter,' Ha' 1599 write(io_unit,14) '# omega_ext (W) = ', external_omega,' Ha' 1600 write(io_unit,10) '# ' 1601 write(io_unit,10) '# ' 1602 write(io_unit,10) '# NOTE: if lmax_model = 0, then eps_model = I, the identity. ' 1603 write(io_unit,10) '# ' 1604 write(io_unit,10) '# ' 1605 write(io_unit,10) '#-----------------------------------------------------------------------------------------------' 1606 write(io_unit,10) ' ' 1607 write(io_unit,30) ' eps_e (Ha) : ', eig(e) 1608 write(io_unit,10) ' ' 1609 write(io_unit,30) ' Sigma_x (Ha) : ', Sigma_x 1610 1611 if (present(Sigma_x_Lanczos_projected) ) then 1612 write(io_unit,30) ' Sigma_x_PROJECTED (Ha) : ', Sigma_x_Lanczos_projected 1613 end if 1614 1615 1616 1617 write(io_unit,30) ' < V_xc >_e (Ha) : ', Vxc_energy 1618 write(io_unit,10) ' ' 1619 write(io_unit,30) ' poles (Ha) : ', pole_energy 1620 write(io_unit,30) ' Sigma_A_1 (Ha) : ', sigma_A_Lanczos 1621 write(io_unit,30) ' Sigma_A_2 (Ha) : ', sigma_A_model_Lanczos 1622 write(io_unit,30) ' Sigma_B_1 (Ha) : ', sigma_B_Lanczos 1623 write(io_unit,30) ' Sigma_B_2 (Ha) : ', sigma_B_model_Lanczos 1624 write(io_unit,10) ' ' 1625 1626 Sigma_c = pole_energy+sigma_A_Lanczos+sigma_A_model_Lanczos+sigma_B_Lanczos+sigma_B_model_Lanczos 1627 1628 write(io_unit,30) ' Sigma_c (Ha) : ',Sigma_c 1629 write(io_unit,10) ' ' 1630 write(io_unit,30) ' E_e (Ha) : ', eig(e)+Sigma_x-Vxc_energy+Sigma_c 1631 1632 1633 close(io_unit) 1634 end if 1635 1636 1637 1638 10 format(A) 1639 12 format(A,ES10.3,A) 1640 14 format(A,ES24.16,A) 1641 25 format(A,I5) 1642 30 format(A,ES24.16) 1643 1644 end subroutine output_results
m_gwls_ComputeCorrelationEnergy/output_Sigma_A_by_eigenvalues [ Functions ]
[ Top ] [ m_gwls_ComputeCorrelationEnergy ] [ Functions ]
NAME
output_Sigma_A_by_eigenvalues
FUNCTION
.
INPUTS
OUTPUT
SOURCE
1731 subroutine output_Sigma_A_by_eigenvalues(n_ext_freq,lmax,external_frequencies,AT_Lanczos,eigenvalues_array,which_case) 1732 !---------------------------------------------------------------------------------------------------- 1733 ! This routine outputs the eigenvalues of the static dielectric matrix, as well as the 1734 ! contributions to Sigma_A, decomposed by eigenvalues. 1735 ! 1736 ! There are three cases to consider: 1737 ! 1 ) no model is being used; we are printing A 1738 ! 2 ) a model is being used; we are printing A1 1739 ! 2 ) a model is being used; we are printing A2 1740 !---------------------------------------------------------------------------------------------------- 1741 1742 integer, intent(in) :: n_ext_freq, lmax, which_case 1743 complex(dpc), intent(in) :: AT_Lanczos(n_ext_freq,lmax) 1744 real(dp), intent(in) :: eigenvalues_array(lmax) 1745 real(dp), intent(in) :: external_frequencies(n_ext_freq) 1746 1747 integer :: iw_ext, l 1748 real(dp) :: external_omega 1749 1750 1751 complex(dpc) :: matrix, eig 1752 1753 integer :: io_unit 1754 character(128) :: filename 1755 1756 ! ************************************************************************* 1757 1758 1759 if (mpi_enreg%me == 0) then 1760 io_unit = get_unit() 1761 1762 1763 do iw_ext = 1, n_ext_freq 1764 1765 if (which_case == 1) then 1766 write(filename,'(A,I0.4,A)') "SIGMA_A_BY_EIGENVALUES_",iw_ext,".dat" 1767 else if (which_case == 2) then 1768 write(filename,'(A,I0.4,A)') "SIGMA_A1_BY_EIGENVALUES_",iw_ext,".dat" 1769 else if (which_case == 3) then 1770 write(filename,'(A,I0.4,A)') "SIGMA_A2_BY_EIGENVALUES_",iw_ext,".dat" 1771 end if 1772 1773 1774 external_omega = external_frequencies(iw_ext) 1775 1776 open(file=filename,status=files_status_new,unit=io_unit) 1777 write(io_unit,10) '#-----------------------------------------------------------------------------------------------' 1778 write(io_unit,10) '#' 1779 write(io_unit,10) '# This file contains the contributions to Sigma_A, the analytical self-energy term,' 1780 write(io_unit,10) '# as a function of the eigenvalue of the dielectric operator.' 1781 write(io_unit,10) '#' 1782 write(io_unit,10) '# Definitions:' 1783 write(io_unit,10) '# ' 1784 write(io_unit,10) '# ' 1785 1786 if (which_case == 1) then 1787 write(io_unit,10) '# Sigma_A = Tr[(eps^{-1}(0) -I ) AT(W) ] ' 1788 write(io_unit,10) '# = sum_{l} ( 1/eps_l -1 ) < V_l | AT(W) | V_l > ' 1789 else if (which_case == 2) then 1790 write(io_unit,10) '# Sigma_A1= Tr[(eps^{-1}(0) - eps_model^{-1}(0) ) A1T(W) ] ' 1791 write(io_unit,10) '# = sum_{l} ( LBDA_l ) < V_l | A1T(W) | V_l > ' 1792 write(io_unit,10) '# where LBDA_l are the eigenvalues of eps^{-1}(0)-eps_model^{-1}(0) ' 1793 else if (which_case == 3) then 1794 write(io_unit,10) '# Sigma_A2= Tr[(eps_model^{-1}(0)- I ) A2T(W) ] ' 1795 write(io_unit,10) '# = sum_{l} (1/eps^{model}_l -1 ) < V_l | A2T(W) | V_l > ' 1796 end if 1797 1798 1799 write(io_unit,10) '# ' 1800 write(io_unit,10) '# ' 1801 write(io_unit,10) '# Parameters: ' 1802 write(io_unit,10) '# ' 1803 1804 if (which_case == 1 .or. which_case == 2) then 1805 write(io_unit,25) '# lmax = ', lmax 1806 else if (which_case == 3) then 1807 write(io_unit,25) '# lmax_model = ', lmax 1808 end if 1809 1810 write(io_unit,25) '# n_ext_freq = ', n_ext_freq 1811 write(io_unit,25) '# iw_ext = ', iw_ext 1812 write(io_unit,14) '# omega_ext (W) = ', external_omega,' Ha' 1813 write(io_unit,10) '# ' 1814 write(io_unit,10) '# ' 1815 if (which_case == 1) then 1816 write(io_unit,10) '# (1/eps_l -1) < V_l | AT(W) | V_l > (Ha) ( 1/eps_l -1 ) '& 1817 & //'< V_l | AT(W) | V_l > (Ha)' 1818 else if (which_case == 2) then 1819 write(io_unit,10) '# LBDA_1 < V_l | A1T(W) | V_l > (Ha) LBDA_l '& 1820 & //'< V_l | A1T(W) | V_l > (Ha)' 1821 else if (which_case == 3) then 1822 write(io_unit,10) '# (1/eps_m_l -1) < V_l | A2T(W) | V_l > (Ha) ( 1/eps_m_l -1 ) '& 1823 & //'< V_l | A2T(W) | V_l > (Ha)' 1824 end if 1825 1826 write(io_unit,10) '# real imaginary real '& 1827 & //'imaginary' 1828 write(io_unit,10) '#---------------------------------------------------------------------------------------------------------'& 1829 & //'----------------' 1830 1831 do l = 1, lmax 1832 1833 matrix = AT_Lanczos(iw_ext,l) 1834 eig = eigenvalues_array(l) 1835 1836 1837 write(io_unit,20) real(eig), matrix, eig*matrix 1838 end do 1839 1840 1841 close(io_unit) 1842 1843 end do 1844 1845 end if 1846 1847 10 format(A) 1848 14 format(A,ES24.16,A) 1849 20 format(ES24.16,2ES24.16,2X,2ES24.16) 1850 25 format(A,I5) 1851 1852 1853 1854 end subroutine output_Sigma_A_by_eigenvalues