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 : !!***
|
