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