Line data Source code
1 : !!****m* ABINIT/m_gwls_GenerateEpsilon
2 : !! NAME
3 : !! m_gwls_GenerateEpsilon
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 :
23 : module m_gwls_GenerateEpsilon
24 : !----------------------------------------------------------------------------------------------------
25 : ! This module contains routines to compute and store the dielectric matrix, which plays a
26 : ! central role in the self energy computations. In particular, global arrays are used to store
27 : ! the static dielectric matrix.
28 : !----------------------------------------------------------------------------------------------------
29 : ! local modules
30 : use m_gwls_utility
31 : use m_gwls_wf
32 : use m_gwls_hamiltonian
33 : use m_gwls_lineqsolver
34 : use m_gwls_TimingLog
35 : use m_gwls_polarisability
36 : use m_gwls_model_polarisability
37 : use m_gwls_GWlanczos
38 : ! Abinit modules
39 : use m_abicore
40 : use defs_basis
41 : use m_dtset
42 :
43 : use m_io_tools, only : get_unit
44 :
45 : implicit none
46 : save
47 : private
48 : !!***
49 :
50 : ! Global arrays
51 :
52 : real(dp), public, allocatable :: epsilon_eigenvalues_0(:) ! eigenvalues of the static dielectric matrix
53 :
54 : complex(dp), public, allocatable :: epsilon_inverse_0(:,:) ! eps^{-1}-1 in diagonal basis
55 :
56 : integer, public :: kmax, nseeds, lmax
57 : integer, public :: first_seed
58 : !!***
59 :
60 : public :: driver_generate_dielectric_matrix
61 : public :: GeneratePrintDielectricEigenvalues
62 : public :: Driver_GeneratePrintDielectricEigenvalues
63 : !!***
64 :
65 : contains
66 :
67 : !!****f* m_gwls_GenerateEpsilon/driver_generate_dielectric_matrix
68 : !! NAME
69 : !! driver_generate_dielectric_matrix
70 : !!
71 : !! FUNCTION
72 : !! .
73 : !!
74 : !! INPUTS
75 : !!
76 : !! OUTPUT
77 : !!
78 : !! SOURCE
79 :
80 13 : subroutine driver_generate_dielectric_matrix(epsilon_matrix_function,nseeds,kmax,&
81 13 : epsilon_eigenvalues,Lbasis,debug)
82 : !----------------------------------------------------------------------
83 : ! This routine computes the Lanczos approximate representation of the
84 : ! implicit dielectic operator and then diagonalizes the banded
85 : ! Lanczos matrix.
86 : !----------------------------------------------------------------------
87 : interface
88 : subroutine epsilon_matrix_function(v_out,v_in,l)
89 : use defs_basis
90 :
91 : integer, intent(in) :: l
92 : complex(dp), intent(out) :: v_out(l)
93 : complex(dp), intent(in) :: v_in(l)
94 :
95 : end subroutine epsilon_matrix_function
96 : end interface
97 :
98 : integer, intent(in) :: nseeds, kmax
99 : logical, intent(in) :: debug
100 :
101 : real (dp), intent(out) :: epsilon_eigenvalues(nseeds*kmax)
102 : complex(dp), intent(out) :: Lbasis(npw_k,nseeds*kmax) ! array containing the Lanczos basis
103 :
104 :
105 : ! local variables
106 :
107 13 : complex(dp), allocatable :: seeds(:,:)
108 13 : complex(dp),allocatable :: alpha(:,:,:)
109 13 : complex(dp),allocatable :: beta (:,:,:)
110 :
111 : integer :: mpi_communicator
112 :
113 : ! *************************************************************************
114 :
115 :
116 : ! The epsilon operator will act in LA mode.
117 13 : mpi_communicator = mpi_enreg%comm_bandfft
118 :
119 :
120 : !Create seeds
121 52 : ABI_MALLOC(seeds,(npw_k,nseeds))
122 13 : call get_seeds(first_seed, nseeds, seeds)
123 :
124 : ! compute the Lanczos basis
125 65 : ABI_MALLOC(alpha,(nseeds,nseeds,kmax))
126 52 : ABI_MALLOC(beta ,(nseeds,nseeds,kmax))
127 :
128 : call block_lanczos_algorithm(mpi_communicator,epsilon_matrix_function,kmax,nseeds,npw_k, &
129 13 : seeds,alpha,beta,Lbasis)
130 :
131 : ! Diagonalize the epsilon matrix, which is banded
132 13 : call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis,epsilon_eigenvalues,debug)
133 :
134 13 : if (debug) then
135 0 : call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,nseeds*kmax,npw_k,Lbasis,epsilon_eigenvalues)
136 : end if
137 :
138 13 : ABI_FREE(seeds)
139 13 : ABI_FREE(alpha)
140 13 : ABI_FREE(beta)
141 :
142 13 : end subroutine driver_generate_dielectric_matrix
143 : !!***
144 :
145 : !!****f* m_gwls_GenerateEpsilon/GeneratePrintDielectricEigenvalues
146 : !! NAME
147 : !! GeneratePrintDielectricEigenvalues
148 : !!
149 : !! FUNCTION
150 : !! .
151 : !!
152 : !! INPUTS
153 : !!
154 : !! OUTPUT
155 : !!
156 : !! SOURCE
157 :
158 2 : subroutine GeneratePrintDielectricEigenvalues(epsilon_matrix_function,nseeds,kmax,output_filename,Lbasis,alpha,beta)
159 : !----------------------------------------------------------------------
160 : ! This routine computes the Lanczos approximate representation of the
161 : ! implicit dielectic operator and then diagonalizes the banded
162 : ! Lanczos matrix.
163 : !----------------------------------------------------------------------
164 : interface
165 : subroutine epsilon_matrix_function(v_out,v_in,l)
166 : use defs_basis
167 :
168 : integer, intent(in) :: l
169 : complex(dp), intent(out) :: v_out(l)
170 : complex(dp), intent(in) :: v_in(l)
171 :
172 : end subroutine epsilon_matrix_function
173 : end interface
174 :
175 : integer, intent(in) :: nseeds, kmax
176 :
177 : character(*), intent(in) :: output_filename
178 :
179 :
180 : complex(dp), intent(out) :: Lbasis(:,:)
181 : complex(dp), intent(out) :: alpha(:,:,:)
182 : complex(dp), intent(out) :: beta (:,:,:)
183 :
184 :
185 : ! local variables
186 :
187 :
188 2 : complex(dp),allocatable :: seeds(:,:)
189 2 : complex(dp),allocatable :: Lbasis_diag(:,:)
190 :
191 :
192 2 : real(dp), allocatable :: psik(:,:)
193 2 : real(dp), allocatable :: psir(:,:,:,:)
194 :
195 2 : real(dp), allocatable :: epsilon_eigenvalues(:)
196 :
197 :
198 : integer :: mpi_communicator
199 : integer :: io_unit
200 : integer :: lmax
201 : integer :: l
202 : integer :: ir1, ir2, ir3
203 : integer :: n1, n2, n3
204 :
205 : real(dp) :: R, G
206 : real(dp) :: sigma_R, sigma_G
207 : real(dp) :: x, y, z
208 :
209 2 : real(dp),allocatable :: G_array(:)
210 2 : real(dp),allocatable :: R_array(:,:,:)
211 :
212 : logical :: debug
213 :
214 : ! *************************************************************************
215 :
216 :
217 2 : debug = .false.
218 2 : lmax = kmax*nseeds
219 2 : mpi_communicator = mpi_enreg%comm_bandfft
220 : !Create seeds
221 8 : ABI_MALLOC(seeds,(npw_k,nseeds))
222 2 : call get_seeds(first_seed, nseeds, seeds)
223 :
224 : ! compute the Lanczos basis
225 8 : ABI_MALLOC(Lbasis_diag,(npw_k,lmax))
226 6 : ABI_MALLOC(epsilon_eigenvalues,(lmax))
227 :
228 6 : ABI_MALLOC(psik,(2,npw_k))
229 10 : ABI_MALLOC(psir,(2,n4,n5,n6))
230 6 : ABI_MALLOC(G_array,(npw_k))
231 10 : ABI_MALLOC(R_array,(n4,n5,n6))
232 :
233 39710 : psir = zero
234 13718 : R_array = zero
235 :
236 2 : n1 = n4-1
237 2 : n2 = n5-1
238 2 : n3 = n6
239 :
240 : ! Generate the Lanczos basis and banded eigenvalue representation
241 2 : call block_lanczos_algorithm(mpi_communicator, epsilon_matrix_function,kmax,nseeds,npw_k, seeds,alpha,beta,Lbasis)
242 :
243 4132 : Lbasis_diag = Lbasis
244 :
245 : ! Diagonalize the epsilon matrix, which is banded
246 2 : call diagonalize_lanczos_banded(kmax,nseeds,npw_k,alpha,beta,Lbasis_diag,epsilon_eigenvalues,debug)
247 :
248 2 : call ritz_analysis_general(mpi_communicator, epsilon_matrix_function,lmax,npw_k,Lbasis_diag,epsilon_eigenvalues)
249 :
250 2 : io_unit = get_unit()
251 2 : open(file=output_filename,status=files_status_new,unit=io_unit)
252 2 : write(io_unit,10) '#----------------------------------------------------------------------------'
253 2 : write(io_unit,10) '# '
254 2 : write(io_unit,10) '# Partial eigenvalues '
255 2 : write(io_unit,10) '# ========================================== '
256 2 : write(io_unit,10) '# '
257 2 : write(io_unit,10) '# Tabulate the eigenvalues of the dielectic matrix, as well as some '
258 2 : write(io_unit,10) '# information regarding the eigenstates. '
259 2 : write(io_unit,10) '# '
260 2 : write(io_unit,10) '# '
261 2 : write(io_unit,10) '# definitions: '
262 2 : write(io_unit,10) '# l index of the eigenvalue '
263 2 : write(io_unit,10) '# eig eigenvalue '
264 2 : write(io_unit,10) '# '
265 2 : write(io_unit,10) '# '
266 2 : write(io_unit,10) '# (the following vectors are expressed in crystal units) '
267 2 : write(io_unit,10) '# '
268 2 : write(io_unit,10) '# R = < l | |r| | l > '
269 2 : write(io_unit,10) '# sigma_R = sqrt{< l | (|r|-R)^2 | l > } '
270 2 : write(io_unit,10) '# '
271 2 : write(io_unit,10) '# G = < l | |G| | l > '
272 2 : write(io_unit,10) '# sigma_G = sqrt{< l | (|G|-G)^2 | l > } '
273 2 : write(io_unit,10) '# '
274 2 : write(io_unit,10) '# '
275 2 : write(io_unit,10) '# l eig r sigma_r G sigma_G '
276 2 : write(io_unit,10) '#----------------------------------------------------------------------------'
277 2 : flush(io_unit)
278 :
279 :
280 516 : G_array(:) = kg_k(1,:)**2+ kg_k(2,:)**2+ kg_k(3,:)**2
281 :
282 516 : G_array(:) = sqrt(G_array(:))
283 :
284 13718 : R_array = zero
285 :
286 38 : do ir3=1,n3
287 :
288 36 : if (ir3 <= n3/2 ) then
289 18 : z = (one*ir3)/(one*n3)
290 : else
291 18 : z = (one*ir3)/(one*n3)-one
292 : end if
293 :
294 686 : do ir2=1,n2
295 :
296 648 : if (ir2 <= n2/2 ) then
297 324 : y = (one*ir2)/(one*n2)
298 : else
299 324 : y = (one*ir2)/(one*n2)-one
300 : end if
301 :
302 12348 : do ir1=1,n1
303 :
304 11664 : if (ir1 <= n1/2 ) then
305 5832 : x = (one*ir1)/(one*n1)
306 : else
307 5832 : x = (one*ir1)/(one*n1)-one
308 : end if
309 :
310 :
311 12312 : R_array(ir1,ir2,ir3) = sqrt(x**2+y**2+z**2)
312 : end do
313 : end do
314 : end do
315 :
316 :
317 :
318 18 : do l=1, lmax
319 :
320 4128 : psik(1,:) = dble (Lbasis_diag(:,l))
321 4128 : psik(2,:) = dimag(Lbasis_diag(:,l))
322 :
323 16 : call g_to_r(psir ,psik)
324 :
325 :
326 4128 : G = sum(G_array(:)*(psik(1,:)**2+psik(2,:)**2))
327 109744 : R = sum(R_array(:,:,:)*(psir(1,:,:,:)**2+psir(2,:,:,:)**2) )*ucvol/nfft
328 :
329 4128 : sigma_G = sqrt(sum((G_array(:) -G)**2*(psik(1,:)**2 +psik(2,:)**2)))
330 109744 : sigma_R = sqrt(sum((R_array(:,:,:)-R)**2*(psir(1,:,:,:)**2+psir(2,:,:,:)**2))*ucvol/nfft)
331 :
332 :
333 :
334 18 : write(io_unit,20) l, epsilon_eigenvalues(l), R,sigma_R, G,sigma_G
335 :
336 : end do
337 :
338 2 : close(io_unit)
339 :
340 2 : ABI_FREE(seeds)
341 2 : ABI_FREE(Lbasis_diag)
342 2 : ABI_FREE(psik)
343 2 : ABI_FREE(psir)
344 2 : ABI_FREE(G_array)
345 2 : ABI_FREE(R_array)
346 :
347 2 : ABI_FREE(epsilon_eigenvalues)
348 :
349 :
350 : 10 format(A)
351 : 20 format(I5,ES24.16,4F12.8)
352 :
353 2 : end subroutine GeneratePrintDielectricEigenvalues
354 : !!***
355 :
356 : !!****f* m_gwls_GenerateEpsilon/Driver_GeneratePrintDielectricEigenvalues
357 : !! NAME
358 : !! Driver_GeneratePrintDielectricEigenvalues
359 : !!
360 : !! FUNCTION
361 : !! .
362 : !!
363 : !! INPUTS
364 : !!
365 : !! OUTPUT
366 : !!
367 : !! SOURCE
368 :
369 1 : subroutine Driver_GeneratePrintDielectricEigenvalues(dtset)
370 : !----------------------------------------------------------------------
371 : ! Compute the eigenvalues of the various dielectric operators
372 : !----------------------------------------------------------------------
373 : type(dataset_type),intent(in) :: dtset
374 :
375 : integer :: kmax_exact, kmax_model, kmax
376 : real(dp) :: second_model_parameter
377 :
378 :
379 : integer :: lm, k, lmax, l1, l2
380 : integer :: io_unit
381 : integer :: io_unit2
382 :
383 : real(dp) :: time1, time2, time
384 : ! local variables
385 : character(128) :: output_filename
386 : character(256) :: timing_string
387 :
388 :
389 : complex(dp), allocatable :: Lbasis_exact(:,:)
390 : complex(dp), allocatable :: Lbasis_model(:,:)
391 1 : complex(dp), allocatable :: sub_Lbasis_exact(:,:)
392 1 : complex(dp), allocatable :: sub_Lbasis_model(:,:)
393 1 : complex(dp), allocatable :: dummy(:,:)
394 1 : complex(dp), allocatable :: dummy2(:,:)
395 1 : complex(dp), allocatable :: dummy3(:,:)
396 : complex(dp), allocatable :: alpha_exact(:,:,:)
397 : complex(dp), allocatable :: beta_exact (:,:,:)
398 : complex(dp), allocatable :: alpha_model(:,:,:)
399 : complex(dp), allocatable :: beta_model (:,:,:)
400 :
401 1 : real(dp), allocatable :: eig_exact(:)
402 1 : real(dp), allocatable :: eig_model(:)
403 :
404 1 : complex(dp), allocatable :: model_epsilon_matrix(:,:)
405 1 : complex(dp), allocatable :: vector(:)
406 :
407 : real(dp) :: tr_eps_1, tr_eps_2, tr_eps_3
408 :
409 :
410 : integer :: lwork, lrwork, liwork, info
411 1 : complex(dp), allocatable :: work(:)
412 1 : real(dp) , allocatable :: rwork(:)
413 1 : integer , allocatable :: iwork(:)
414 :
415 : integer :: debug_unit
416 : character(50) :: debug_filename
417 : ! *************************************************************************
418 :
419 1 : kmax_exact = dtset%gwls_stern_kmax
420 1 : kmax_model = dtset%gwls_kmax_complement
421 :
422 : !second_model_parameter = dtset%gwls_second_model_parameter
423 1 : second_model_parameter = zero
424 :
425 :
426 : ! global stuff
427 1 : nseeds = dtset%gwls_nseeds
428 1 : first_seed = dtset%gwls_first_seed
429 1 : e = dtset%gwls_band_index
430 :
431 :
432 :
433 4 : ABI_MALLOC(Lbasis_exact,(npw_k,kmax_exact*nseeds))
434 4 : ABI_MALLOC(Lbasis_model,(npw_k,kmax_model*nseeds))
435 :
436 5 : ABI_MALLOC(alpha_exact, (nseeds,nseeds,kmax_exact))
437 4 : ABI_MALLOC(beta_exact , (nseeds,nseeds,kmax_exact))
438 5 : ABI_MALLOC(alpha_model, (nseeds,nseeds,kmax_model))
439 4 : ABI_MALLOC(beta_model , (nseeds,nseeds,kmax_model))
440 :
441 :
442 : ! set omega=0 for exact dielectric operator
443 1 : call set_dielectric_function_frequency([zero,zero])
444 :
445 :
446 :
447 1 : call cpu_time(time1)
448 1 : output_filename = 'EIGENVALUES_EXACT.dat'
449 : call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_k, nseeds, kmax_exact, &
450 1 : output_filename, Lbasis_exact, alpha_exact, beta_exact)
451 :
452 :
453 :
454 1 : call cpu_time(time2)
455 1 : time = time2-time1
456 1 : write(timing_string,'(A)') "Time to compute the EXACT Static Dielectric Matrix : "
457 1 : call write_timing_log(timing_string,time)
458 :
459 :
460 :
461 1 : call cpu_time(time1)
462 1 : call setup_Pk_model(zero,second_model_parameter)
463 1 : output_filename = 'EIGENVALUES_MODEL.dat'
464 : call GeneratePrintDielectricEigenvalues(matrix_function_epsilon_model_operator, nseeds, kmax_model, &
465 1 : &output_filename, Lbasis_model, alpha_model, beta_model)
466 1 : call cpu_time(time2)
467 1 : time = time2-time1
468 1 : write(timing_string,'(A)') "Time to compute the MODEL Static Dielectric Matrix : "
469 1 : call write_timing_log(timing_string,time)
470 :
471 :
472 1 : call cpu_time(time1)
473 1 : if (kmax_exact <= kmax_model) then
474 : kmax = kmax_exact
475 : else
476 : kmax = kmax_model
477 : end if
478 1 : lmax = nseeds*kmax
479 :
480 : ! Build model operator matrix elements in the exact basis
481 4 : ABI_MALLOC(model_epsilon_matrix, (lmax,lmax))
482 3 : ABI_MALLOC(vector, (npw_k))
483 :
484 73 : model_epsilon_matrix = cmplx_0
485 :
486 9 : do l2 =1 , lmax
487 8 : call matrix_function_epsilon_model_operator(vector ,Lbasis_exact(:,l2),npw_k)
488 :
489 :
490 73 : do l1 =1, lmax
491 :
492 72 : model_epsilon_matrix(l1, l2) = complex_vector_product(Lbasis_exact(:,l1),vector,npw_k)
493 :
494 : end do
495 :
496 :
497 : end do
498 :
499 1 : ABI_FREE(vector)
500 :
501 1 : call cpu_time(time2)
502 1 : time = time2-time1
503 1 : write(timing_string,'(A)') "Compute MODEL matrix elements in EXACT basis : "
504 1 : call write_timing_log(timing_string,time)
505 :
506 :
507 :
508 :
509 : ! Compare the traces
510 :
511 :
512 1 : io_unit = get_unit()
513 1 : open(file='DIELECTRIC_TRACE.dat',status=files_status_new,unit=io_unit)
514 1 : write(io_unit,10) '#----------------------------------------------------------------------------'
515 1 : write(io_unit,10) '# '
516 1 : write(io_unit,10) '# Partial traces '
517 1 : write(io_unit,10) '# ========================================== '
518 1 : write(io_unit,10) '# '
519 1 : write(io_unit,10) '# Tabulate the trace of various operators, as function of the number '
520 1 : write(io_unit,10) '# of lanczos steps performed. '
521 1 : write(io_unit,10) '# '
522 1 : write(io_unit,10) '# '
523 1 : write(io_unit,10) '# NOTES: '
524 1 : write(io_unit,10) '# Tr[1-eps^{-1}] is evaluated in the Lanczos basis of eps '
525 1 : write(io_unit,10) '# Tr[1-eps_m^{-1}] is evaluated in the Lanczos basis of eps_m '
526 1 : write(io_unit,10) '# Tr[eps_m^{-1}-eps^{-1}] is evaluated in the Lanczos basis of eps '
527 1 : write(io_unit,10) '# '
528 1 : write(io_unit,10) '# '
529 1 : write(io_unit,10) '# k Tr[1-eps^{-1}] Tr[1-eps_m^{-1}] Tr[eps_m^{-1}-eps^{-1}]'
530 1 : write(io_unit,10) '#----------------------------------------------------------------------------'
531 1 : flush(io_unit)
532 :
533 :
534 1 : io_unit2 = get_unit()
535 1 : open(file='RPA_ENERGY.dat',status=files_status_new,unit=io_unit2)
536 1 : write(io_unit2,10) '#----------------------------------------------------------------------------'
537 1 : write(io_unit2,10) '# '
538 1 : write(io_unit2,10) '# RPA TOTAL ENERGY '
539 1 : write(io_unit2,10) '# ========================================== '
540 1 : write(io_unit2,10) '# '
541 1 : write(io_unit2,10) '# It can be shown that the correlation energy, within the RPA, is given '
542 1 : write(io_unit2,10) '# by: '
543 1 : write(io_unit2,10) '# E_c = int_0^{infty} dw/(2pi) Tr[ ln(eps{iw)}+1-eps(iw)] '
544 1 : write(io_unit2,10) '# '
545 1 : write(io_unit2,10) '# As a gauge of what can be expected as far as convergence is concerned, '
546 1 : write(io_unit2,10) '# the following will be printed. '
547 1 : write(io_unit2,10) '# '
548 1 : write(io_unit2,10) '# I_1 = Tr[ ln(eps) + 1 - eps ] '
549 1 : write(io_unit2,10) '# I_2 = Tr[ ln(eps_m) + 1 - eps_m ] '
550 1 : write(io_unit2,10) '# I_3 = Tr[ ln(eps^{-1}_m . eps ) + eps_m - eps ] '
551 1 : write(io_unit2,10) '# '
552 1 : write(io_unit2,10) '# '
553 1 : write(io_unit2,10) '# '
554 1 : write(io_unit2,10) '# k I_1 I_2 I_3 '
555 1 : write(io_unit2,10) '#----------------------------------------------------------------------------'
556 1 : flush(io_unit2)
557 :
558 :
559 : ! Iterate every 10 values of k max, or else the linear algebra gets too expensive...
560 1 : do k = 4, kmax, 4
561 :
562 8 : ABI_MALLOC(sub_Lbasis_exact,(npw_k,k*nseeds))
563 6 : ABI_MALLOC(sub_Lbasis_model,(npw_k,k*nseeds))
564 :
565 :
566 6 : ABI_MALLOC(eig_exact,(k*nseeds))
567 4 : ABI_MALLOC(eig_model,(k*nseeds))
568 8 : ABI_MALLOC(dummy,(k*nseeds,k*nseeds))
569 6 : ABI_MALLOC(dummy2,(k*nseeds,k*nseeds))
570 :
571 3098 : sub_Lbasis_exact(:,:) = Lbasis_exact(:,1:k*nseeds)
572 3098 : sub_Lbasis_model(:,:) = Lbasis_model(:,1:k*nseeds)
573 :
574 : ! Diagonalize the epsilon matrix, which is banded
575 2 : call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_exact(:,:,1:k),beta_exact(:,:,1:k),sub_Lbasis_exact,eig_exact,.false.)
576 2 : call diagonalize_lanczos_banded(k,nseeds,npw_k,alpha_model(:,:,1:k),beta_model(:,:,1:k),sub_Lbasis_model,eig_model,.false.)
577 :
578 14 : tr_eps_1 = sum(one-one/eig_exact(:))
579 14 : tr_eps_2 = sum(one-one/eig_model(:))
580 :
581 14 : tr_eps_3 = -sum(one/eig_exact(:))
582 :
583 94 : dummy(:,:) = model_epsilon_matrix(1:k*nseeds, 1:k*nseeds)
584 :
585 2 : call driver_invert_positive_definite_hermitian_matrix(dummy,k*nseeds)
586 :
587 14 : do lm = 1, k*nseeds
588 :
589 14 : tr_eps_3 = tr_eps_3 + dble(dummy(lm,lm))
590 : end do
591 :
592 :
593 2 : write(io_unit,20) k, tr_eps_1, tr_eps_2, tr_eps_3
594 2 : flush(io_unit)
595 :
596 :
597 :
598 14 : tr_eps_1 = sum(log(eig_exact(:))+one-eig_exact(:))
599 14 : tr_eps_2 = sum(log(eig_model(:))+one-eig_model(:))
600 :
601 94 : dummy2(:,:) = zero
602 2 : tr_eps_3 = zero
603 14 : do lm = 1, k*nseeds
604 12 : dummy2(lm,lm) = eig_exact(lm)
605 14 : tr_eps_3 = tr_eps_3 + dble(model_epsilon_matrix(lm,lm)) - dble(eig_exact(lm))
606 : end do
607 :
608 8 : ABI_MALLOC(dummy3,(k*nseeds,k*nseeds))
609 : call ZGEMM( 'N', & ! Hermitian conjugate the first array
610 : 'N', & ! Leave second array as is
611 : k*nseeds, & ! the number of rows of the matrix op( A )
612 : k*nseeds, & ! the number of columns of the matrix op( B )
613 : k*nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
614 : cmplx_1, & ! alpha constant
615 : dummy2, & ! matrix A
616 : k*nseeds, & ! LDA
617 : dummy, & ! matrix B
618 : k*nseeds, & ! LDB
619 : cmplx_0, & ! beta constant
620 : dummy3, & ! matrix C
621 2 : k*nseeds) ! LDC
622 :
623 94 : dummy2(:,:) = dummy3(:,:)
624 2 : ABI_FREE(dummy3)
625 :
626 : ! find eigenvalues
627 : !call heevd(dummy2, eig_exact)
628 :
629 2 : lwork = k*nseeds+1
630 2 : lrwork = k*nseeds
631 2 : liwork = 1
632 :
633 6 : ABI_MALLOC(work,(lwork))
634 6 : ABI_MALLOC(rwork,(lrwork))
635 2 : ABI_MALLOC(iwork,(liwork))
636 :
637 2 : call zheevd('N', 'U',k*nseeds, dummy2, k*nseeds, eig_exact, work, lwork, rwork, lrwork, iwork, liwork, info)
638 2 : if ( info /= 0) then
639 0 : debug_unit = get_unit()
640 0 : write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
641 :
642 0 : open(debug_unit,file=trim(debug_filename),status='unknown')
643 :
644 0 : write(debug_unit,'(A)') '*************************************************************************************'
645 0 : write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZHEEVD(1), gwls_GenerateEpsilon'
646 0 : write(debug_unit,'(A)') '*************************************************************************************'
647 :
648 0 : close(debug_unit)
649 :
650 : end if
651 :
652 :
653 :
654 :
655 :
656 14 : tr_eps_3 = tr_eps_3 + sum(log(eig_exact))
657 :
658 2 : write(io_unit2,20) k, tr_eps_1, tr_eps_2, tr_eps_3
659 2 : flush(io_unit2)
660 :
661 :
662 :
663 2 : ABI_FREE(work)
664 2 : ABI_FREE(rwork)
665 2 : ABI_FREE(iwork)
666 :
667 :
668 :
669 2 : ABI_FREE(sub_Lbasis_exact)
670 2 : ABI_FREE(sub_Lbasis_model)
671 :
672 2 : ABI_FREE(eig_exact)
673 2 : ABI_FREE(eig_model)
674 2 : ABI_FREE(dummy)
675 2 : ABI_FREE(dummy2)
676 : end do
677 :
678 1 : close(io_unit)
679 1 : close(io_unit2)
680 :
681 1 : call cpu_time(time2)
682 1 : time = time2-time1
683 1 : write(timing_string,'(A)') "Time to compute the TRACES of the Dielectric Matrices: "
684 1 : call write_timing_log(timing_string,time)
685 :
686 1 : ABI_FREE(Lbasis_exact)
687 1 : ABI_FREE(Lbasis_model)
688 1 : ABI_FREE(alpha_exact)
689 1 : ABI_FREE(beta_exact )
690 1 : ABI_FREE(alpha_model)
691 1 : ABI_FREE(beta_model )
692 1 : ABI_FREE(model_epsilon_matrix)
693 :
694 :
695 :
696 : 10 format(A)
697 : 20 format(I5,3ES24.16)
698 :
699 1 : end subroutine Driver_GeneratePrintDielectricEigenvalues
700 : !!***
701 :
702 : end module m_gwls_GenerateEpsilon
703 : !!***
|
