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