Branch data Line data Source code
1 : : !!****m* ABINIT/m_gwls_DielectricArray
2 : : !! NAME
3 : : !! m_gwls_DielectricArray
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_DielectricArray
24 : : !----------------------------------------------------------------------------------------------------
25 : : ! This module generates and stores the arrays
26 : : !
27 : : ! { eps^{-1}(iw)-eps_model^{-1}(iw) } in Lanczos basis
28 : : ! { eps_model^{-1}(iw) - 1 } in model Lanczos basis
29 : : !
30 : : ! It makes sense to build these only once, as they do not depend on the external frequency.
31 : : !
32 : : !----------------------------------------------------------------------------------------------------
33 : :
34 : : ! local modules
35 : : use m_gwls_utility
36 : : use m_gwls_wf
37 : : use m_gwls_valenceWavefunctions
38 : : use m_gwls_hamiltonian
39 : : use m_gwls_lineqsolver
40 : : use m_gwls_polarisability
41 : : use m_gwls_model_polarisability
42 : : use m_gwls_GenerateEpsilon
43 : : use m_gwls_TimingLog
44 : : use m_gwls_QR_factorization
45 : : use m_gwls_LanczosBasis
46 : :
47 : : ! abinit modules
48 : : use defs_basis
49 : : use m_abicore
50 : : use m_xmpi
51 : : use m_cgtools
52 : :
53 : : use m_time, only : timab
54 : : use m_io_tools, only: get_unit
55 : : use m_gaussian_quadrature, only: get_frequencies_and_weights_legendre
56 : :
57 : :
58 : : implicit none
59 : : save
60 : :
61 : : private
62 : : !!***
63 : :
64 : : ! Frequencies and weights for Legendre integration
65 : : real(dp), public, allocatable :: list_omega(:)
66 : : real(dp), public, allocatable :: list_weights(:)
67 : :
68 : : ! Arrays to store the combinations of dielectric operators
69 : : complex(dp), public, allocatable :: model_dielectric_Lanczos_basis(:,:,:)
70 : : complex(dp), public, allocatable :: projected_dielectric_Lanczos_basis(:,:,:)
71 : : complex(dp), public, allocatable :: eps_m1_minus_eps_model_m1(:,:,:)
72 : :
73 : : complex(dp), public, allocatable :: eps_model_m1_minus_one(:,:,:)
74 : : complex(dp),public, allocatable :: eps_model_m1_minus_one_DISTR(:,:,:)
75 : :
76 : :
77 : : ! dimensions of blocks in the model dielectric matrix
78 : : integer,public :: nbdblock_epsilon
79 : : integer,public :: blocksize_epsilon
80 : : logical,public, allocatable :: model_lanczos_vector_belongs_to_this_node(:)
81 : : integer,public, allocatable :: model_lanczos_vector_index(:)
82 : :
83 : :
84 : : ! Arrays necessary to project the Sternheimer equation within
85 : : ! the computation of the dielectric matrix.
86 : : complex(dp), public, allocatable :: projected_epsilon_M_matrix(:,:,:)
87 : : complex(dp), public, allocatable :: projected_epsilon_B_matrix(:,:,:)
88 : : complex(dp), public, allocatable :: projected_epsilon_G_matrix(:,:,:)
89 : :
90 : : integer,public, allocatable :: list_lsolutions_EpsilonProjected(:)
91 : : !!***
92 : :
93 : : public :: generate_frequencies_and_weights
94 : : public :: cleanup_projected_Sternheimer_epsilon
95 : : public :: compute_eps_m1_minus_eps_model_m1
96 : : public :: compute_eps_m1_minus_one
97 : : public :: compute_eps_model_m1_minus_one
98 : : public :: ProjectedSternheimerEpsilon
99 : : !!***
100 : :
101 : : contains
102 : :
103 : : !!****f* m_hamiltonian/generate_frequencies_and_weights
104 : : !! NAME
105 : : !! generate_frequencies_and_weights
106 : : !!
107 : : !! FUNCTION
108 : : !! .
109 : : !!
110 : : !! INPUTS
111 : : !!
112 : : !! OUTPUT
113 : : !!
114 : : !! SOURCE
115 : :
116 : 7 : subroutine generate_frequencies_and_weights(npt_gauss)
117 : : !--------------------------------------------------------------------------------
118 : : !
119 : : ! This subroutine computes the frequencies and weights necessary for Gauss-Legendre
120 : : ! quadrature, and stores the results in module arrays.
121 : : !
122 : : !--------------------------------------------------------------------------------
123 : : integer, intent(in) :: npt_gauss
124 : :
125 : :
126 : 7 : real(dp), allocatable :: list_omega_tmp(:)
127 : 7 : real(dp), allocatable :: list_weights_tmp(:)
128 : :
129 : : integer :: i
130 : :
131 : : ! *************************************************************************
132 : :
133 [ - + + - : 21 : ABI_MALLOC(list_omega_tmp, (npt_gauss))
+ - + - -
+ - + ]
134 [ - + + - : 14 : ABI_MALLOC(list_weights_tmp, (npt_gauss))
- + - + ]
135 : :
136 : 7 : call get_frequencies_and_weights_legendre(npt_gauss,list_omega_tmp,list_weights_tmp)
137 : :
138 : :
139 [ + - - + : 21 : ABI_MALLOC(list_omega, (npt_gauss+1))
+ - + - +
- - + - +
- + ]
140 [ - + + - : 14 : ABI_MALLOC(list_weights, (npt_gauss+1))
- + - + -
+ ]
141 : :
142 : : ! make sure the first frequency in zero!
143 : 7 : list_omega(1) = zero
144 : 7 : list_weights(1) = zero
145 : :
146 [ + + ]: 77 : do i = 1,npt_gauss
147 : :
148 : : ! inverse the order of the frequency points, as they come out
149 : : ! in reverse order from the generating subroutine
150 : 70 : list_omega (i+1) = list_omega_tmp (npt_gauss+1-i)
151 : 77 : list_weights(i+1) = list_weights_tmp(npt_gauss+1-i)
152 : :
153 : : end do
154 : :
155 : :
156 [ - + ]: 7 : ABI_FREE(list_omega_tmp)
157 [ - + ]: 7 : ABI_FREE(list_weights_tmp)
158 : :
159 : :
160 : 7 : end subroutine generate_frequencies_and_weights
161 : : !!***
162 : :
163 : : !!****f* m_hamiltonian/compute_eps_m1_minus_eps_model_m1
164 : : !! NAME
165 : : !! compute_eps_m1_minus_eps_model_m1
166 : : !!
167 : : !! FUNCTION
168 : : !! .
169 : : !!
170 : : !! INPUTS
171 : : !!
172 : : !! OUTPUT
173 : : !!
174 : : !! SOURCE
175 : :
176 : 6 : subroutine compute_eps_m1_minus_eps_model_m1(lmax, npt_gauss)
177 : : !----------------------------------------------------------------------------------------------------
178 : : !
179 : : ! This subroutine computes the array
180 : : !
181 : : ! eps^{-1}(iw) - eps_model^{-1}(iw),
182 : : !
183 : : ! for all relevant frequencies in the Lanczos basis.
184 : : !----------------------------------------------------------------------------------------------------
185 : : integer , intent(in) :: lmax, npt_gauss
186 : :
187 : : character(256) :: timing_string
188 : : real(dp) :: time1, time2
189 : : real(dp) :: time
190 : :
191 : : integer :: iw, l
192 : 6 : complex(dp),allocatable :: dummy_matrix(:,:)
193 : 6 : complex(dp),allocatable :: iden(:,:)
194 : : ! *************************************************************************
195 : :
196 : 6 : timing_string = "#"
197 : 6 : call write_text_block_in_Timing_log(timing_string)
198 : 6 : timing_string = "# Computing eps^{-1}(iw) - eps_model^{-1}(iw) "
199 : 6 : call write_text_block_in_Timing_log(timing_string)
200 : 6 : timing_string = "#"
201 : 6 : call write_text_block_in_Timing_log(timing_string)
202 : :
203 : :
204 : 6 : call cpu_time(time1)
205 : : ! Allocate the module array
206 [ - + + - : 30 : ABI_MALLOC(eps_m1_minus_eps_model_m1, (lmax,lmax,npt_gauss+1))
- + + - +
- - + + -
+ - + - -
+ - + -
+ ]
207 [ - + + - : 24 : ABI_MALLOC(dummy_matrix, (lmax,lmax))
- + + - +
- + - - +
- + ]
208 [ - + + - : 18 : ABI_MALLOC(iden, (lmax,lmax))
- + + - -
+ + - ]
209 : :
210 : :
211 [ + + + + ]: 438 : iden = cmplx_0
212 : :
213 [ + + ]: 54 : do l = 1, lmax
214 : 54 : iden(l,l) = cmplx_1
215 : : end do
216 : :
217 [ + + ]: 72 : do iw = 1, npt_gauss + 1
218 : :
219 : :
220 [ + + + + ]: 4818 : dummy_matrix(:,:) = projected_dielectric_Lanczos_basis(:,:,iw)
221 : :
222 : 66 : call driver_invert_positive_definite_hermitian_matrix(dummy_matrix,lmax)
223 : :
224 : :
225 [ + + + + ]: 4818 : eps_m1_minus_eps_model_m1(:,:,iw) = dummy_matrix(:,:)
226 : :
227 [ + + + + ]: 4818 : dummy_matrix(:,:) = model_dielectric_Lanczos_basis(:,:,iw)
228 : :
229 : :
230 : :
231 : 66 : call driver_invert_positive_definite_hermitian_matrix(dummy_matrix,lmax)
232 : :
233 [ + + + + ]: 4824 : eps_m1_minus_eps_model_m1(:,:,iw) = eps_m1_minus_eps_model_m1(:,:,iw) - dummy_matrix(:,:)
234 : :
235 : :
236 : : end do
237 : :
238 : :
239 : :
240 : :
241 : : ! Deallocate the arrays which are no longer needed
242 [ - + ]: 6 : ABI_FREE(model_dielectric_Lanczos_basis)
243 [ - + ]: 6 : ABI_FREE(projected_dielectric_Lanczos_basis)
244 : :
245 : :
246 : :
247 [ - + ]: 6 : ABI_FREE(dummy_matrix)
248 [ - + ]: 6 : ABI_FREE(iden)
249 : :
250 : :
251 : :
252 : 6 : call cpu_time(time2)
253 : 6 : time = time2-time1
254 : :
255 : 6 : timing_string = "# Total time : "
256 : 6 : call write_timing_log(timing_string,time)
257 : :
258 : :
259 : :
260 : :
261 : 6 : end subroutine compute_eps_m1_minus_eps_model_m1
262 : : !!***
263 : :
264 : : !!****f* m_hamiltonian/compute_eps_m1_minus_one
265 : : !! NAME
266 : : !! compute_eps_m1_minus_one
267 : : !!
268 : : !! FUNCTION
269 : : !! .
270 : : !!
271 : : !! INPUTS
272 : : !!
273 : : !! OUTPUT
274 : : !!
275 : : !! SOURCE
276 : :
277 : 1 : subroutine compute_eps_m1_minus_one(lmax, npt_gauss)
278 : : !----------------------------------------------------------------------------------------------------
279 : : !
280 : : ! This subroutine computes the array
281 : : !
282 : : ! eps^{-1}(iw) - I
283 : : !
284 : : ! for all relevant frequencies in the Lanczos basis.
285 : : !----------------------------------------------------------------------------------------------------
286 : : integer , intent(in) :: lmax, npt_gauss
287 : :
288 : : character(256) :: timing_string
289 : : real(dp) :: time1, time2
290 : : real(dp) :: time
291 : :
292 : : integer :: iw, l
293 : 1 : complex(dp),allocatable :: dummy_matrix(:,:)
294 : 1 : complex(dp),allocatable :: iden(:,:)
295 : : ! *************************************************************************
296 : :
297 : 1 : timing_string = "#"
298 : 1 : call write_text_block_in_Timing_log(timing_string)
299 : 1 : timing_string = "# Computing eps^{-1}(iw) - I "
300 : 1 : call write_text_block_in_Timing_log(timing_string)
301 : 1 : timing_string = "#"
302 : 1 : call write_text_block_in_Timing_log(timing_string)
303 : :
304 : 1 : call cpu_time(time1)
305 : : ! Allocate the module array
306 : :
307 : : ! The array eps_m1_minus_eps_model_m1 will be used to store
308 : : ! eps^{-1}-1; we can think of eps_model = I in this case.
309 [ - + + - : 5 : ABI_MALLOC(eps_m1_minus_eps_model_m1, (lmax,lmax,npt_gauss+1))
- + + - +
- - + + -
+ - + - -
+ - + -
+ ]
310 : :
311 [ - + + - : 4 : ABI_MALLOC(dummy_matrix, (lmax,lmax))
- + + - +
- + - - +
- + ]
312 [ - + + - : 3 : ABI_MALLOC(iden, (lmax,lmax))
- + + - -
+ + - ]
313 : :
314 [ + + + + ]: 73 : iden = cmplx_0
315 : :
316 [ + + ]: 9 : do l = 1, lmax
317 : 9 : iden(l,l) = cmplx_1
318 : : end do
319 : :
320 : :
321 [ + + ]: 12 : do iw = 1, npt_gauss + 1
322 : :
323 [ + + + + ]: 803 : dummy_matrix(:,:) = projected_dielectric_Lanczos_basis(:,:,iw)
324 : :
325 : 11 : call driver_invert_positive_definite_hermitian_matrix(dummy_matrix,lmax)
326 : :
327 [ + + + + ]: 804 : eps_m1_minus_eps_model_m1(:,:,iw) = dummy_matrix(:,:)-iden(:,:)
328 : :
329 : : end do
330 : :
331 : :
332 : : ! Deallocate the arrays which are no longer needed
333 [ - + ]: 1 : ABI_FREE(projected_dielectric_Lanczos_basis)
334 : :
335 [ - + ]: 1 : ABI_FREE(dummy_matrix)
336 [ - + ]: 1 : ABI_FREE(iden)
337 : :
338 : 1 : call cpu_time(time2)
339 : 1 : time = time2-time1
340 : :
341 : 1 : timing_string = "# Total time : "
342 : 1 : call write_timing_log(timing_string,time)
343 : :
344 : 1 : end subroutine compute_eps_m1_minus_one
345 : : !!***
346 : :
347 : : !!****f* m_hamiltonian/compute_eps_model_m1_minus_one
348 : : !! NAME
349 : : !! compute_eps_model_m1_minus_one
350 : : !!
351 : : !! FUNCTION
352 : : !! .
353 : : !!
354 : : !! INPUTS
355 : : !!
356 : : !! OUTPUT
357 : : !!
358 : : !! SOURCE
359 : :
360 : 6 : subroutine compute_eps_model_m1_minus_one(lmax_model, npt_gauss, second_model_parameter, epsilon_model_eigenvalues_0)
361 : : !----------------------------------------------------------------------------------------------------
362 : : !
363 : : ! This subroutine computes the array
364 : : !
365 : : ! eps_model^{-1}(iw) - 1
366 : : !
367 : : ! for all relevant frequencies in the model Lanczos basis.
368 : : !
369 : : ! This array can potentially get very large, as the complementary basis gets big to achieve
370 : : ! convergence. Correspondingly, it makes sense to DISTRIBUTE this array on all processors.
371 : : !
372 : : ! The algorithm will go as follows:
373 : : !
374 : : ! eps_m = 1 - V . P . V, V = sqrt{vc}
375 : : !
376 : : ! I ) compute VPV, storing blocks on different processors:
377 : : !
378 : : ! VPV = [ ------- -------- ] = VPV[lc, nB, nW]
379 : : ! | [| block | block | ]
380 : : ! lc [| 1 | 2 | ... ]
381 : : ! | [| | | ]
382 : : ! | [| | | ]
383 : : ! | [ ------- --------- ]
384 : : ! bsize
385 : : !
386 : : ! This construction *does* involve a fair bit of communications, but it takes
387 : : ! a lot less RAM!
388 : : !
389 : : ! II) once VPV is constructed, do, one frequency at a time:
390 : : ! - import all blocks to the HEAD processor
391 : : ! - compute eps_m = 1- VPV
392 : : ! - invert eps_m^{-1}
393 : : ! - subtract identity eps_m^{-1} - 1
394 : : ! - redistribute, block by block
395 : : !
396 : : ! Doing this frequency by frequency will reduce the RAM weight on the head node.
397 : : !
398 : : !
399 : : !----------------------------------------------------------------------------------------------------
400 : : integer, intent(in) :: lmax_model, npt_gauss
401 : : real(dp), intent(in) :: second_model_parameter
402 : : real(dp), intent(in) :: epsilon_model_eigenvalues_0(lmax_model)
403 : :
404 : : integer :: l, l1, l2
405 : : integer :: iw
406 : : integer :: v
407 : : integer :: ierr
408 : :
409 : 6 : real(dp), allocatable :: psikg_valence(:,:)
410 : 6 : real(dp), allocatable :: psir_valence(:,:,:,:)
411 : :
412 : :
413 : 6 : real(dp), allocatable :: psik_wrk(:,:)
414 : 6 : real(dp), allocatable :: psikb_wrk(:,:)
415 : 6 : real(dp), allocatable :: psikg_wrk(:,:)
416 : 6 : real(dp), allocatable :: psikg_tmp(:,:)
417 : :
418 : 6 : complex(dp),allocatable :: local_Lbasis_conjugated(:,:)
419 : :
420 : :
421 : 6 : complex(dp),allocatable :: VPV(:,:,:)
422 : 6 : real(dp), allocatable :: re_buffer(:,:), im_buffer(:,:)
423 : :
424 : 6 : complex(dp),allocatable :: vpv_row(:)
425 : :
426 : :
427 : 6 : complex(dp),allocatable :: epsilon_head(:,:)
428 : :
429 : 6 : real(dp), allocatable :: re_BUFFER_head(:,:)
430 : 6 : real(dp), allocatable :: im_BUFFER_head(:,:)
431 : :
432 : :
433 : :
434 : 6 : complex(dp),allocatable :: YL(:)
435 : :
436 : : character(256) :: timing_string
437 : : real(dp) :: time1, time2, time
438 : : real(dp) :: fft_time1, fft_time2, fft_time
439 : : real(dp) :: prod_time1, prod_time2, prod_time
440 : :
441 : : complex(dp) :: z
442 : :
443 : :
444 : : integer :: iblk_lanczos, nbdblock_lanczos
445 : : integer :: mb
446 : : integer :: lb
447 : : integer :: ik
448 : :
449 : : integer :: mpi_communicator
450 : : integer :: mpi_nproc
451 : : integer :: mpi_rank
452 : : integer :: mpi_head_rank
453 : :
454 : :
455 : 6 : integer,allocatable :: sendcounts(:), displs(:)
456 : :
457 : : integer :: sendcount, recvcount
458 : :
459 : :
460 : : logical :: head
461 : :
462 : :
463 : : ! timing
464 : : real(dp) :: tsec(2)
465 : : integer :: GWLS_TIMAB, OPTION_TIMAB
466 : :
467 : : ! *************************************************************************
468 : :
469 : :
470 : 6 : GWLS_TIMAB = 1541
471 : 6 : OPTION_TIMAB = 1
472 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
473 : :
474 : :
475 : :
476 : 6 : timing_string = "#"
477 : 6 : call write_text_block_in_Timing_log(timing_string)
478 : 6 : timing_string = "# computing eps_model_m1_minus_one"
479 : 6 : call write_text_block_in_Timing_log(timing_string)
480 : 6 : timing_string = "#"
481 : 6 : call write_text_block_in_Timing_log(timing_string)
482 : :
483 : :
484 : : ! Number of blocks of lanczos vectors (blocksize = npband)
485 : 6 : nbdblock_lanczos = lmax_model/blocksize
486 [ - + ]: 6 : if (modulo(lmax_model,blocksize) /= 0) nbdblock_lanczos = nbdblock_lanczos + 1
487 : :
488 : :
489 : : ! communicator
490 : 6 : mpi_communicator = mpi_enreg%comm_bandfft
491 : :
492 : : ! total number of processors in the communicator
493 : 6 : mpi_nproc = xmpi_comm_size(mpi_communicator )
494 : :
495 : : ! what is the rank of this processor?
496 : 6 : mpi_rank = xmpi_comm_rank(mpi_communicator )
497 : :
498 : : ! rank of the "head" processor
499 : 6 : mpi_head_rank = 0
500 : :
501 : :
502 : : ! number of blocks in the model dielectric matrix, which is equal to the number of processors
503 : 6 : nbdblock_epsilon = mpi_nproc
504 : :
505 : :
506 : : ! number of vectors in every block
507 : 6 : blocksize_epsilon = lmax_model/mpi_nproc
508 [ - + ]: 6 : if (modulo(lmax_model,mpi_nproc) /= 0) blocksize_epsilon = blocksize_epsilon + 1
509 : :
510 : : ! attribute blocks to every nodes, and tabulate ownership in logical array
511 : : ! This is not *the most efficient* implementation possible, but it is convenient
512 [ - + + - : 18 : ABI_MALLOC( model_lanczos_vector_belongs_to_this_node, (lmax_model))
+ - + - -
+ - + -
+ ]
513 [ - + + - : 12 : ABI_MALLOC( model_lanczos_vector_index, (lmax_model))
- + - + -
+ ]
514 : :
515 : :
516 [ + + ]: 54 : model_lanczos_vector_index = 0
517 [ + + ]: 54 : model_lanczos_vector_belongs_to_this_node = .false.
518 : :
519 [ + + ]: 54 : do l =1, lmax_model
520 [ + + ]: 54 : if (mpi_rank == (l-1)/blocksize_epsilon) then
521 : 24 : model_lanczos_vector_belongs_to_this_node(l) = .true.
522 : 24 : model_lanczos_vector_index(l) = l-mpi_rank*blocksize_epsilon
523 : : end if
524 : : end do
525 : :
526 : : !write(100+mpi_rank,*) "model_lanczos_vector_belongs_to_this_node = ",model_lanczos_vector_belongs_to_this_node(:)
527 : : !write(100+mpi_rank,*) "model_lanczos_vector_index = ",model_lanczos_vector_index(:)
528 : : !flush(100+mpi_rank)
529 : :
530 : : ! Prepare the array that will contain the matrix elements of the model operator
531 [ - + + - : 30 : ABI_MALLOC(VPV, (lmax_model,blocksize_epsilon,npt_gauss+1))
- + + - +
- - + + -
+ - + - -
+ - + ]
532 : :
533 [ + + + + : 2448 : VPV(:,:,:) = cmplx_0
+ + ]
534 : :
535 [ - + + - : 18 : ABI_MALLOC(vpv_row, (lmax_model))
+ - + - -
+ - + ]
536 : :
537 : :
538 : : ! various working arrays
539 : 6 : GWLS_TIMAB = 1542
540 : : OPTION_TIMAB = 1
541 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
542 : :
543 [ - + + - : 18 : ABI_MALLOC(psikg_valence,(2,npw_g))
+ - + - -
+ - + ]
544 [ - + + - : 30 : ABI_MALLOC(psir_valence ,(2,n4,n5,n6))
- + + - -
+ + - + -
+ - - + -
+ ]
545 : :
546 : :
547 [ - + + - : 18 : ABI_MALLOC(psik_wrk, (2,npw_k))
+ - + - -
+ - + ]
548 [ - + + - : 18 : ABI_MALLOC(psikb_wrk, (2,npw_kb))
+ - + - -
+ - + ]
549 [ - + + - : 12 : ABI_MALLOC(psikg_wrk, (2,npw_g))
- + - + ]
550 [ - + + - : 12 : ABI_MALLOC(psikg_tmp, (2,npw_g))
- + - + ]
551 : :
552 [ - + + - : 24 : ABI_MALLOC(local_Lbasis_conjugated,(npw_k,lmax_model))
- + + - +
- + - - +
- + ]
553 [ - + + - : 18 : ABI_MALLOC(YL,(npw_k))
+ - + - -
+ - + ]
554 : :
555 : 6 : OPTION_TIMAB = 2
556 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
557 : :
558 : :
559 : :
560 : 6 : fft_time = zero
561 : 6 : prod_time = zero
562 : :
563 : :
564 : 6 : call cpu_time(time1)
565 : : ! loop on all valence bands
566 : :
567 [ + + ]: 30 : do v = 1, nbandv
568 : :
569 : :
570 : : ! copy pre-calculated valence state in this covenient local array
571 [ + + + + ]: 12360 : psikg_valence(:,:) = kernel_wavefunctions_FFT(:,:,v)
572 : :
573 : : ! compute fourier transform of valence state, and conjugate
574 : 24 : call g_to_r(psir_valence,psikg_valence)
575 [ + + + + : 164616 : psir_valence(2,:,:,:) = -psir_valence(2,:,:,:)
+ + ]
576 : :
577 : : !--------------------------------------------------------------------------
578 : : !
579 : : ! Step 1: build the modified basis Pc . [ (V^{1/2} l) . psi_v^*].
580 : : !
581 : : !--------------------------------------------------------------------------
582 : 24 : GWLS_TIMAB = 1543
583 : 24 : OPTION_TIMAB = 1
584 : 24 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
585 : :
586 : : ! loop on all blocks of lanczos vectors
587 [ + + ]: 152 : do iblk_lanczos = 1, nbdblock_lanczos
588 : : ! loop on all states within this block
589 [ + + ]: 320 : do mb = 1, blocksize
590 : :
591 : : ! Determine the index of the Lanczos vector
592 : 192 : l = (iblk_lanczos-1)*blocksize + mb
593 : :
594 [ - + ]: 192 : if ( l <= lmax_model) then
595 [ + + ]: 24864 : psik_wrk(1,:) = dble (Lbasis_model_lanczos(:,l))
596 [ + + ]: 24864 : psik_wrk(2,:) = dimag(Lbasis_model_lanczos(:,l))
597 : : else
598 [ # # # # ]: 0 : psik_wrk(:,:) = zero
599 : : end if
600 : :
601 : : ! apply Coulomb potential
602 : 192 : call sqrt_vc_k(psik_wrk)
603 : :
604 : : ! Store in array of blocks of wavefunctions
605 [ + + + + ]: 74336 : psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:)
606 : :
607 : : end do ! mb
608 : :
609 : 128 : call cpu_time(fft_time1)
610 : :
611 : : ! Transform to FFT representation
612 : 128 : call wf_block_distribute(psikb_wrk, psikg_wrk,1) ! LA -> FFT
613 : :
614 : :
615 : : ! generate the vector Pc [ (sqrt_V_c.l) psi_v^*]
616 : :
617 : : ! Compute the real space product, and return to k space, in FFT configuration
618 : 128 : call gr_to_g(psikg_tmp,psir_valence, psikg_wrk)
619 : :
620 : :
621 : 128 : call cpu_time(fft_time2)
622 : 128 : fft_time = fft_time + fft_time2-fft_time1
623 : :
624 : : ! project
625 : 128 : call pc_k_valence_kernel(psikg_tmp)
626 : :
627 : : ! Return to LA configuration
628 : :
629 : : ! Transform to FFT representation
630 : 128 : call wf_block_distribute(psikb_wrk, psikg_tmp, 2) ! FFT -> LA
631 : :
632 [ + + ]: 344 : do mb = 1, blocksize
633 : :
634 : : ! Determine the index of the Lanczos vector
635 : 192 : l = (iblk_lanczos-1)*blocksize + mb
636 : :
637 : :
638 [ + - ]: 320 : if ( l <= lmax_model) then
639 [ + + + + ]: 74208 : psik_wrk(:,:) = psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k)
640 [ + + ]: 24864 : local_Lbasis_conjugated(:,l) = cmplx_1*psik_wrk(1,:)+cmplx_i*psik_wrk(2,:)
641 : : end if
642 : :
643 : :
644 : : end do ! mb
645 : :
646 : : end do !iblk_lanczos
647 : 24 : OPTION_TIMAB = 2
648 : 24 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
649 : :
650 : :
651 : : !--------------------------------------------------------------------------
652 : : ! Step 2: Now that we have the modified basis, compute the matrix
653 : : ! elements of the model dielectric operator
654 : : !
655 : : !
656 : : !--------------------------------------------------------------------------
657 : 24 : GWLS_TIMAB = 1544
658 : 24 : OPTION_TIMAB = 1
659 : 24 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
660 : :
661 [ + + ]: 264 : do iw = 2, npt_gauss+1
662 : :
663 : 240 : call setup_Pk_model(list_omega(iw),second_model_parameter)
664 : :
665 : 240 : call cpu_time(prod_time1)
666 [ + + ]: 2160 : do l1 = 1, lmax_model
667 : :
668 : : ! Apply core function Y to left-vector; conjugate
669 [ + + ]: 248640 : do ik = 1, npw_k
670 : 248640 : YL(ik) = model_Y_LA(ik)*conjg(local_Lbasis_conjugated(ik,l1))
671 : : end do
672 : :
673 : : ! Only compute lower diagonal part of matrix; epsilon is hermitian conjugate!
674 [ + + ]: 17280 : vpv_row = cmplx_0
675 [ + + ]: 10560 : do l2 = 1, l1
676 [ + + ]: 1120800 : do ik = 1, npw_k
677 : 1118880 : vpv_row(l2) = vpv_row(l2) + YL(ik)*local_Lbasis_conjugated(ik,l2)
678 : : end do
679 : : end do
680 : :
681 : : ! the code below is twice as long!
682 : : !call ZGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY )
683 : : !call ZGEMV ( 'T', npw_k, lmax_model, cmplx_1, local_Lbasis_conjugated, npw_k, YL, 1, cmplx_0, vpv_row, 1)
684 : :
685 : : !do l2 = 1, l1
686 : : ! eps_model_m1_minus_one(l1,l2,iw) = eps_model_m1_minus_one(l1,l2,iw) &
687 : : ! -complex_vector_product(YL, local_Lbasis_conjugated(:,l2),npw_k)
688 : : !end do
689 : :
690 : : ! Sum on all processors, making sure all processors have the total vpv_row
691 : 1920 : call xmpi_sum(vpv_row, mpi_communicator, ierr) ! sum on all processors for LA configuration
692 : :
693 : : ! Each processor takes its slice!
694 [ + + ]: 10800 : do l2 =1, l1
695 [ + + ]: 10560 : if ( model_lanczos_vector_belongs_to_this_node(l2) ) then
696 : 4320 : lb = model_lanczos_vector_index(l2)
697 : 4320 : VPV(l1,lb,iw) = VPV(l1,lb,iw) + vpv_row(l2)
698 : : end if
699 : : end do
700 : :
701 : : end do
702 : 240 : call cpu_time(prod_time2)
703 : :
704 : 264 : prod_time = prod_time + prod_time2-prod_time1
705 : :
706 : : end do ! iw
707 : 24 : OPTION_TIMAB = 2
708 : 30 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
709 : :
710 : : end do ! v
711 : :
712 : 6 : call cpu_time(time2)
713 : :
714 : 6 : time = time2-time1
715 : 6 : timing_string = "# computing VPV : "
716 : 6 : call write_timing_log(timing_string,time)
717 : :
718 : 6 : timing_string = "# --- of which is FFT transforms : "
719 : 6 : call write_timing_log(timing_string,fft_time)
720 : :
721 : 6 : timing_string = "# --- of which is products : "
722 : 6 : call write_timing_log(timing_string,prod_time)
723 : :
724 : :
725 : :
726 : :
727 : :
728 : 6 : GWLS_TIMAB = 1542
729 : 6 : OPTION_TIMAB = 1
730 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
731 [ - + ]: 6 : ABI_FREE(local_Lbasis_conjugated)
732 [ - + ]: 6 : ABI_FREE(YL)
733 [ - + ]: 6 : ABI_FREE(psikg_valence)
734 [ - + ]: 6 : ABI_FREE(psir_valence)
735 [ - + ]: 6 : ABI_FREE(psik_wrk)
736 [ - + ]: 6 : ABI_FREE(psikb_wrk)
737 [ - + ]: 6 : ABI_FREE(psikg_wrk)
738 [ - + ]: 6 : ABI_FREE(psikg_tmp)
739 [ - + ]: 6 : ABI_FREE(vpv_row)
740 : :
741 : 6 : OPTION_TIMAB = 2
742 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
743 : :
744 : :
745 : :
746 : 6 : call cpu_time(time1)
747 : 6 : GWLS_TIMAB = 1547
748 : 6 : OPTION_TIMAB = 1
749 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
750 : :
751 : :
752 : : !--------------------------------------------------------------------------------
753 : : !
754 : : !
755 : : ! Gather dielectric matrix on head node, invert, and re-distribute. This
756 : : ! Saves a lot of RAM, without needing the full machinery of ScaLAPACK.
757 : : !
758 : : !--------------------------------------------------------------------------------
759 : :
760 : :
761 [ - + + - : 30 : ABI_MALLOC(eps_model_m1_minus_one_DISTR, (lmax_model,blocksize_epsilon,npt_gauss+1))
- + + - -
+ + - + -
+ - - + -
+ - + ]
762 [ - + + - : 24 : ABI_MALLOC(re_buffer, (lmax_model,blocksize_epsilon))
- + + - +
- + - - +
- + ]
763 [ - + + - : 18 : ABI_MALLOC(im_buffer, (lmax_model,blocksize_epsilon))
- + + - -
+ - + ]
764 : :
765 [ + + + + : 2448 : eps_model_m1_minus_one_DISTR(:,:,:) = cmplx_0
+ + ]
766 : :
767 : : ! Define the head node, which will invert the dielectric matrices
768 : 6 : head = mpi_rank == mpi_head_rank
769 : :
770 : : ! Amount of data received and sent
771 : 6 : sendcount = lmax_model*blocksize_epsilon
772 : 6 : recvcount = lmax_model*blocksize_epsilon
773 : :
774 [ - + + - : 18 : ABI_MALLOC(sendcounts,(mpi_nproc))
+ - + - -
+ - + ]
775 [ - + + - : 12 : ABI_MALLOC(displs ,(mpi_nproc))
- + - + ]
776 [ + + ]: 24 : sendcounts(:) = sendcount
777 : :
778 [ + + ]: 24 : do l =1, mpi_nproc
779 : 24 : displs(l) = (l-1)*sendcount
780 : : end do
781 : :
782 [ + + ]: 6 : if (head) then
783 : : ! build and invert the dielectric array
784 : : ! careful! lmax_model not necessarily equal to blocksize_epsilon*nbdblock_epsilon
785 [ - + + - : 12 : ABI_MALLOC(re_BUFFER_head, (lmax_model, blocksize_epsilon*nbdblock_epsilon))
- + + - +
- + - - +
- + ]
786 [ - + + - : 9 : ABI_MALLOC(im_BUFFER_head, (lmax_model, blocksize_epsilon*nbdblock_epsilon))
- + + - -
+ - + ]
787 [ - + + - : 12 : ABI_MALLOC(epsilon_head, (lmax_model, lmax_model))
- + + - +
- + - - +
- + ]
788 : : else
789 : : !This looks superfluous and it is on large number of systems, but sending these
790 : : !unallocated in xmpi_scatterv caused 'cannot allocate memory' cryptic errors on
791 : : !several parallel test farm computers (cronos_gcc46_paral, petrus_nag, inca_gcc44_sdebug)
792 [ - + ]: 3 : ABI_MALLOC(re_BUFFER_head, (1,1))
793 [ - + ]: 3 : ABI_MALLOC(im_BUFFER_head, (1,1))
794 [ - + ]: 3 : ABI_MALLOC(epsilon_head, (1,1))
795 : : end if
796 : :
797 : : ! Do one frequency at a time, to avoid overflowing the RAM
798 [ + + ]: 72 : do iw = 1, npt_gauss+1
799 : : ! Gather, except for static case
800 [ + + ]: 66 : if ( iw /=1 ) then
801 : : ! Gather VPV on head node, for this frequency
802 [ + - + + : 2220 : call xmpi_gather(dble(VPV(:,:,iw)), sendcount , re_BUFFER_head, recvcount, mpi_head_rank, mpi_communicator,ierr)
+ + ]
803 [ + + + + ]: 2220 : call xmpi_gather(dimag(VPV(:,:,iw)), sendcount , im_BUFFER_head, recvcount, mpi_head_rank, mpi_communicator,ierr)
804 : : end if
805 : :
806 [ + + ]: 66 : if ( head ) then
807 : :
808 : : ! fill the dielectric matrix
809 : :
810 [ + + + + ]: 2409 : epsilon_head(:,:) = cmplx_0
811 [ + + ]: 33 : if (iw ==1) then
812 : : ! STATIC CASE, diagonal matrix
813 [ + + ]: 27 : do l= 1, lmax_model
814 : 27 : epsilon_head(l,l) = cmplx_1/epsilon_model_eigenvalues_0(l)-cmplx_1
815 : : end do
816 : :
817 : : else
818 : : ! DYNAMIC CASE, compute
819 [ + + ]: 270 : do l1 =1, lmax_model
820 [ + + ]: 1320 : do l2 =1, l1
821 : 1080 : z = -cmplx_1*re_BUFFER_head(l1,l2)-cmplx_i*im_BUFFER_head(l1,l2)
822 : 1080 : epsilon_head(l1,l2) = z
823 : 1320 : epsilon_head(l2,l1) = conjg(z)
824 : : end do
825 : 270 : epsilon_head(l1,l1) = epsilon_head(l1,l1) + cmplx_1
826 : : end do
827 : :
828 : : ! invert the matrix
829 : 30 : call driver_invert_positive_definite_hermitian_matrix(epsilon_head,lmax_model)
830 : :
831 : : ! subtract identity
832 [ + + ]: 270 : do l =1, lmax_model
833 : 270 : epsilon_head(l,l) = epsilon_head(l,l) - cmplx_1
834 : : end do
835 : : end if
836 : :
837 : : ! copy in head buffer
838 [ + + + + ]: 2409 : re_BUFFER_head(:,:) = zero
839 [ + + + + ]: 2409 : im_BUFFER_head(:,:) = zero
840 [ + + ]: 297 : do l1 =1, lmax_model
841 [ + + ]: 2409 : do l2 =1, lmax_model
842 : 2112 : z = epsilon_head(l1,l2)
843 : 2112 : re_BUFFER_head(l1,l2) = dble(z)
844 : 2376 : im_BUFFER_head(l1,l2) = dimag(z)
845 : : end do
846 : : end do
847 : :
848 : : end if
849 : :
850 : : ! Scatter back the data on the head to all processors
851 : 66 : call xmpi_scatterv(re_BUFFER_head, sendcounts, displs, re_buffer, recvcount, mpi_head_rank, mpi_communicator, ierr)
852 : 66 : call xmpi_scatterv(im_BUFFER_head, sendcounts, displs, im_buffer, recvcount, mpi_head_rank, mpi_communicator, ierr)
853 : :
854 [ + + + + ]: 2448 : eps_model_m1_minus_one_DISTR(:,:,iw) = cmplx_1*re_buffer(:,:) + cmplx_i*im_buffer(:,:)
855 : :
856 : : end do
857 : :
858 [ - + ]: 6 : ABI_FREE(re_BUFFER_head)
859 [ - + ]: 6 : ABI_FREE(im_BUFFER_head)
860 [ - + ]: 6 : ABI_FREE(epsilon_head)
861 : :
862 : :
863 : :
864 : : !================================================================================
865 : : !
866 : : ! For debugging purposes, store distributed dielectric matrix back in the
867 : : ! complete local copies, to insure the rest of the code works.
868 : : !
869 : : !================================================================================
870 : :
871 : : if (.false.) then
872 : : ! Prepare the array that will contain the matrix elements of the model operator
873 : : ! THIS IS ONLY FOR THE REST OF THE CODE TO WORK; WE WILL REMOVE THIS
874 : : ! TO SAVE RAM LATER
875 : : ABI_MALLOC(eps_model_m1_minus_one, (lmax_model,lmax_model,npt_gauss+1))
876 : :
877 : : ! initialize the array with zeros
878 : : eps_model_m1_minus_one = cmplx_0
879 : :
880 : : ! Amount of data received and sent
881 : : sendcount = lmax_model*blocksize_epsilon
882 : : recvcount = lmax_model*blocksize_epsilon*nbdblock_epsilon
883 : :
884 : :
885 : : ABI_MALLOC(re_BUFFER_head, (lmax_model,blocksize_epsilon*nbdblock_epsilon))
886 : : ABI_MALLOC(im_BUFFER_head, (lmax_model,blocksize_epsilon*nbdblock_epsilon))
887 : :
888 : : do iw = 1, npt_gauss+1
889 : :
890 : : re_buffer(:,:) = dble( eps_model_m1_minus_one_DISTR(:,:,iw))
891 : : im_buffer(:,:) = dimag(eps_model_m1_minus_one_DISTR(:,:,iw))
892 : :
893 : : call xmpi_allgather(re_buffer,sendcount,re_BUFFER_head,mpi_communicator,ierr)
894 : : call xmpi_allgather(im_buffer,sendcount,im_BUFFER_head,mpi_communicator,ierr)
895 : :
896 : : do l = 1, lmax_model
897 : : eps_model_m1_minus_one(:,l,iw) = cmplx_1*re_BUFFER_head(:,l)+ cmplx_i*im_BUFFER_head(:,l)
898 : : end do
899 : : end do
900 : :
901 : : ABI_FREE(re_BUFFER_head)
902 : : ABI_FREE(im_BUFFER_head)
903 : : end if
904 : :
905 : 6 : call cpu_time(time2)
906 : 6 : OPTION_TIMAB = 2
907 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
908 : :
909 : : !================================================================================
910 : : !
911 : : !================================================================================
912 : :
913 : :
914 : 6 : time = time2-time1
915 : 6 : timing_string = "# inverting / distributing : "
916 : :
917 : 6 : call write_timing_log(timing_string,time)
918 : :
919 : :
920 : :
921 [ - + ]: 6 : ABI_FREE(re_buffer )
922 [ - + ]: 6 : ABI_FREE(im_buffer )
923 [ - + ]: 6 : ABI_FREE(sendcounts)
924 [ - + ]: 6 : ABI_FREE(displs )
925 [ - + ]: 6 : ABI_FREE(VPV )
926 : :
927 : :
928 : :
929 : 6 : GWLS_TIMAB = 1541
930 : : OPTION_TIMAB = 2
931 : 6 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
932 : :
933 : :
934 : 6 : end subroutine compute_eps_model_m1_minus_one
935 : : !!***
936 : :
937 : : !!****f* m_hamiltonian/cleanup_projected_Sternheimer_epsilon
938 : : !! NAME
939 : : !! cleanup_projected_Sternheimer_epsilon
940 : : !!
941 : : !! FUNCTION
942 : : !! .
943 : : !!
944 : : !! INPUTS
945 : : !!
946 : : !! OUTPUT
947 : : !!
948 : : !! SOURCE
949 : :
950 : 7 : subroutine cleanup_projected_Sternheimer_epsilon
951 : :
952 : : ! *************************************************************************
953 : :
954 [ - + ]: 7 : ABI_SFREE(projected_epsilon_M_matrix)
955 [ - + ]: 7 : ABI_SFREE(projected_epsilon_B_matrix)
956 [ - + ]: 7 : ABI_SFREE(projected_epsilon_G_matrix)
957 [ + - ]: 7 : ABI_SFREE(eps_m1_minus_eps_model_m1)
958 [ + - ]: 7 : ABI_SFREE(list_omega)
959 [ + - ]: 7 : ABI_SFREE(list_weights)
960 [ + + ]: 7 : ABI_SFREE(eps_model_m1_minus_one_DISTR)
961 [ + + ]: 7 : ABI_SFREE(model_lanczos_vector_belongs_to_this_node)
962 [ + + ]: 7 : ABI_SFREE(model_lanczos_vector_index)
963 : :
964 : 7 : end subroutine cleanup_projected_Sternheimer_epsilon
965 : : !!***
966 : :
967 : :
968 : : !!****f* m_hamiltonian/ProjectedSternheimerEpsilon
969 : : !! NAME
970 : : !! ProjectedSternheimerEpsilon
971 : : !!
972 : : !! FUNCTION
973 : : !! .
974 : : !!
975 : : !! INPUTS
976 : : !!
977 : : !! OUTPUT
978 : : !!
979 : : !! SOURCE
980 : :
981 : 7 : subroutine ProjectedSternheimerEpsilon(lmax, npt_gauss, second_model_parameter, &
982 : 7 : list_projection_frequencies,nfrequencies,&
983 : 7 : epsilon_eigenvalues_0,debug,use_model)
984 : : !----------------------------------------------------------------------------------------------------
985 : : ! This subroutine combines in a single subprogram the jobs of previous routines
986 : : !
987 : : ! - setup_projected_Sternheimer_epsilon
988 : : ! - compute_projected_Sternheimer_epsilon
989 : : !
990 : : ! The purpose of this combination is to avoid independent loops on nbandv, requiring the
991 : : ! arrays
992 : : ! projected_epsilon_M_matrix
993 : : ! projected_epsilon_B_matrix
994 : : ! projected_epsilon_G_matrix
995 : : !
996 : : ! from scaling like N^3, which grows very large with problem size.
997 : : !
998 : : ! Thus, this routine:
999 : : !
1000 : : ! - Computes the frequency-dependent dielectric matrix in the Lanczos basis, using
1001 : : ! the projected Sternheimer equations.
1002 : : !
1003 : : ! - Computes the frequency-dependent MODEL dielectric matrix in the complementary Lanczos basis.
1004 : : !
1005 : : !
1006 : : ! This routine will be verbose and write log files; indeed, large jobs crash in here, it will
1007 : : ! be important to know where/why!
1008 : : !
1009 : : ! The subroutine also computes the matrix elements on epsilon_model(iw) in the Lanczos basis;
1010 : : ! this is done here to avoid preforming direct products with the valence states again later.
1011 : : !----------------------------------------------------------------------------------------------------
1012 : : real(dp), parameter :: svd_tolerance = 1.0e-16_dp
1013 : :
1014 : : integer, intent(in) :: lmax, npt_gauss
1015 : : integer, intent(in) :: nfrequencies
1016 : : real(dp), intent(in) :: list_projection_frequencies(nfrequencies)
1017 : : logical, intent(in) :: debug
1018 : : real(dp), intent(in) :: epsilon_eigenvalues_0(lmax)
1019 : :
1020 : : logical,optional,intent(in) :: use_model
1021 : :
1022 : : real(dp), intent(in) :: second_model_parameter
1023 : :
1024 : :
1025 : : integer :: l, l1, l2
1026 : : integer :: i, iw, v
1027 : : integer :: recy_i
1028 : : integer :: lsolutions_max, lsolutions, ls
1029 : : integer :: projection
1030 : :
1031 : 7 : complex(dp), allocatable :: sternheimer_A0(:,:)
1032 : 7 : complex(dp), allocatable :: sternheimer_A(:,:)
1033 : 7 : complex(dp), allocatable :: sternheimer_B(:,:)
1034 : 7 : complex(dp), allocatable :: sternheimer_X(:,:)
1035 : 7 : complex(dp), allocatable :: sternheimer_G(:,:)
1036 : :
1037 : :
1038 : 7 : complex(dp), allocatable :: dummy_tmp_1(:,:)
1039 : 7 : complex(dp), allocatable :: dummy_tmp_2(:,:)
1040 : :
1041 : 7 : integer, allocatable :: ipiv(:)
1042 : :
1043 : :
1044 : :
1045 : 7 : complex(dp),allocatable :: local_Lbasis(:,:)
1046 : 7 : complex(dp),allocatable :: local_Lbasis_conjugated(:,:)
1047 : 7 : complex(dp),allocatable :: YL(:)
1048 : :
1049 : 7 : real(dp), allocatable :: psikg_in(:,:), psikg_out(:,:)
1050 : :
1051 : 7 : real(dp), allocatable :: psik_wrk(:,:), psikg_wrk(:,:), psikb_wrk(:,:)
1052 : 7 : real(dp), allocatable :: psi_gamma_l1(:,:), psi_gamma_l2(:,:)
1053 : :
1054 : 7 : real(dp), allocatable :: psikg_valence(:,:)
1055 : 7 : real(dp), allocatable :: psir_valence(:,:,:,:)
1056 : :
1057 : 7 : real(dp), allocatable :: psi_rhs(:,:,:)
1058 : :
1059 : 7 : real(dp), allocatable :: psikg_VL(:,:)
1060 : :
1061 : :
1062 : 7 : complex(dp), allocatable :: check_matrix(:,:), check_matrix2(:,:)
1063 : 7 : complex(dp), allocatable :: c_sternheimer_solutions(:,:)
1064 : 7 : complex(dp), allocatable :: QR_orthonormal_basis(:,:)
1065 : :
1066 : 7 : complex(dp), allocatable :: svd_matrix(:,:)
1067 : 7 : real (dp ), allocatable :: svd_values(:)
1068 : :
1069 : :
1070 : : integer :: iblk_lanczos, nbdblock_lanczos
1071 : : integer :: iblk_solutions, nbdblock_solutions
1072 : : integer :: mb
1073 : :
1074 : : character(128) :: filename
1075 : : logical :: file_exists
1076 : : integer :: io_unit
1077 : :
1078 : : character(128) :: filename_log
1079 : : integer :: io_unit_log
1080 : :
1081 : :
1082 : : real(dp) :: omega
1083 : :
1084 : : character(256) :: timing_string
1085 : : real(dp) :: time1, time2
1086 : : real(dp) :: time_exact
1087 : :
1088 : :
1089 : : integer :: info
1090 : : integer :: ierr
1091 : :
1092 : : real(dp) :: z(2)
1093 : :
1094 : :
1095 : : logical :: omega_is_imaginary
1096 : : real(dp) :: omega0
1097 : :
1098 : : logical :: model
1099 : : logical :: write_debug
1100 : :
1101 : :
1102 : : integer :: mpi_communicator, mpi_rank, mpi_group
1103 : :
1104 : : ! *************************************************************************
1105 : :
1106 : :
1107 : : !================================================================================
1108 : : ! Prepare MPI information
1109 : : !================================================================================
1110 : :
1111 : : ! for LA configuration ,The processors communicate over band+FFT
1112 : 7 : mpi_communicator = mpi_enreg%comm_bandfft
1113 : :
1114 : : ! what is the rank of this processor, within its group?
1115 : 7 : mpi_rank = mpi_enreg%me_fft
1116 : :
1117 : : ! Which group does this processor belong to, given the communicator?
1118 : 7 : mpi_group = mpi_enreg%me_band
1119 : :
1120 : :
1121 : :
1122 : :
1123 : : !================================================================================
1124 : : ! Setup a log file, to keep track of the algorithm
1125 : : !================================================================================
1126 : :
1127 : :
1128 : 7 : write(filename_log,'(A,I4.4,A)') 'ProjectedSternheimerEpsilon_PROC=',mpi_enreg%me,'.log'
1129 : :
1130 : 7 : io_unit_log = get_unit()
1131 : 7 : open(io_unit_log,file=filename_log,status=files_status_new)
1132 : 7 : write(io_unit_log,10) ''
1133 : 7 : write(io_unit_log,10) '#===================================================================================================='
1134 : 7 : write(io_unit_log,10) "# ProjectedSternheimerEpsilon: log file "
1135 : 7 : write(io_unit_log,10) "# ------------------------------------------------------------------- "
1136 : 7 : write(io_unit_log,10) "# "
1137 : 7 : write(io_unit_log,10) "# This file tracks the algorithm in the routine ProjectedSternheimerEpsilon. The goal is to "
1138 : 7 : write(io_unit_log,10) "# establish where the algorithm crashes if it does, and/or to track are far along the code is. "
1139 : 7 : write(io_unit_log,10) '#'
1140 : 7 : write(io_unit_log,10) '# MPI data for this process:'
1141 : 7 : write(io_unit_log,10) '#'
1142 : 7 : write(io_unit_log,22) '# mpi_rank :',mpi_rank,' (rank of this processor in its band group)'
1143 : 7 : write(io_unit_log,22) '# mpi_group:',mpi_group,' (band group to which this processor belongs)'
1144 : 7 : write(io_unit_log,10) '#===================================================================================================='
1145 : 7 : flush(io_unit_log)
1146 : :
1147 : :
1148 : : !================================================================================
1149 : : ! Setup timing; prepare arrays
1150 : : !================================================================================
1151 : :
1152 : 7 : write(io_unit_log,10) " - Preparing and allocating arrays ...."
1153 : 7 : flush(io_unit_log)
1154 : :
1155 : :
1156 : :
1157 : 7 : timing_string = "#"
1158 : 7 : call write_text_block_in_Timing_log(timing_string)
1159 : 7 : timing_string = "# ProjectedSternheimerEpsilon "
1160 : 7 : call write_text_block_in_Timing_log(timing_string)
1161 : 7 : timing_string = "#"
1162 : 7 : call write_text_block_in_Timing_log(timing_string)
1163 : :
1164 : : ! Allocate the module array
1165 [ - + + - : 35 : ABI_MALLOC(projected_dielectric_Lanczos_basis, (lmax,lmax,npt_gauss+1))
- + + - +
- - + + -
+ - + - -
+ - + -
+ ]
1166 : :
1167 [ + + + + : 5628 : projected_dielectric_Lanczos_basis(:,:,:) = cmplx_0
+ + ]
1168 : :
1169 : : ! initialize zero frequency with exact solution
1170 [ + + ]: 63 : do l = 1, lmax
1171 : 63 : projected_dielectric_Lanczos_basis(l,l,1) = cmplx_1*epsilon_eigenvalues_0(l)
1172 : : end do
1173 : :
1174 : :
1175 : : ! initialize other frequencies with the identity
1176 [ + + ]: 77 : do iw = 2, npt_gauss + 1
1177 [ + + ]: 637 : do l = 1, lmax
1178 : 630 : projected_dielectric_Lanczos_basis(l,l,iw) = cmplx_1
1179 : : end do
1180 : : end do
1181 : :
1182 : 7 : time_exact = zero
1183 : :
1184 : :
1185 : : !================================================================================
1186 : : ! Parallelisation of the code is subtle; Hamiltonian must act over
1187 : : ! FFT rows, we must be careful with memory, etc...
1188 : : !
1189 : : ! We will parallelise in block of Lanczos vectors, not over bands.
1190 : : !================================================================================
1191 : :
1192 : : ! Number of blocks of lanczos vectors
1193 : 7 : nbdblock_lanczos = lmax/blocksize
1194 [ - + ]: 7 : if (modulo(lmax,blocksize) /= 0) nbdblock_lanczos = nbdblock_lanczos + 1
1195 : :
1196 : :
1197 [ + - ]: 7 : if (present(use_model)) then
1198 : 7 : model = use_model
1199 : : else
1200 : : model = .true.
1201 : : end if
1202 : :
1203 [ + + ]: 7 : if (model) then
1204 : : ! Prepare the array that will contain the matrix elements of the model operator
1205 [ - + + - : 24 : ABI_MALLOC(model_dielectric_Lanczos_basis, (lmax,lmax,npt_gauss+1))
- + + - -
+ + - - +
- + - + ]
1206 : :
1207 : : ! initialize with zero. NOT with the identity, in order to avoid extra communications
1208 : : ! (see below)
1209 [ + + + + : 4824 : model_dielectric_Lanczos_basis(:,:,:) = cmplx_0
+ + ]
1210 : :
1211 : : end if
1212 : :
1213 : : ! various working arrays
1214 : :
1215 [ - + + - : 21 : ABI_MALLOC(psikg_valence ,(2,npw_g))
+ - + - -
+ - + ]
1216 [ - + + - : 35 : ABI_MALLOC(psir_valence ,(2,n4,n5,n6))
- + + - -
+ + - + -
+ - - + -
+ ]
1217 : :
1218 : :
1219 [ - + + - : 14 : ABI_MALLOC(psikg_VL ,(2,npw_g))
- + - + ]
1220 : :
1221 : :
1222 [ - + + - : 21 : ABI_MALLOC(psik_wrk ,(2,npw_k))
+ - + - -
+ - + ]
1223 [ - + + - : 21 : ABI_MALLOC(psikb_wrk ,(2,npw_kb))
+ - + - -
+ - + ]
1224 [ - + + - : 14 : ABI_MALLOC(psikg_wrk ,(2,npw_g))
- + - + ]
1225 : :
1226 : :
1227 [ - + + - : 14 : ABI_MALLOC(psi_gamma_l1 ,(2,npw_k))
- + - + ]
1228 [ - + + - : 14 : ABI_MALLOC(psi_gamma_l2 ,(2,npw_k))
- + - + ]
1229 : :
1230 [ - + + - : 28 : ABI_MALLOC(psi_rhs ,(2,npw_k,lmax))
- + + - +
- + - - +
- + ]
1231 : :
1232 : :
1233 [ - + + - : 14 : ABI_MALLOC(psikg_in ,(2,npw_g))
- + - + ]
1234 [ - + + - : 14 : ABI_MALLOC(psikg_out ,(2,npw_g))
- + - + ]
1235 : :
1236 : :
1237 : : ! maximal possible dimension of the solution space
1238 : : ! +1 because the solutions at $\omega=\infty$ are free.
1239 : : ! +1 if recycling is activated, because the solutions at $\omega=0$ are then available.
1240 : 7 : i=1
1241 [ + - ]: 7 : if(dtset%gwls_recycle == 1 .or. dtset%gwls_recycle == 2) then
1242 : 7 : i=2
1243 : : end if
1244 : 7 : lsolutions_max = lmax*(nfrequencies+i)
1245 : :
1246 : :
1247 [ - + + - : 28 : ABI_MALLOC(local_Lbasis, (npw_k,lmax))
- + + - +
- + - - +
- + ]
1248 [ - + + - : 21 : ABI_MALLOC(local_Lbasis_conjugated,(npw_k,lmax))
- + + - -
+ - + ]
1249 [ - + + - : 21 : ABI_MALLOC(YL,(npw_k))
+ - + - -
+ - + ]
1250 : :
1251 [ - + + - : 28 : ABI_MALLOC(c_sternheimer_solutions,(npw_k,lsolutions_max))
- + + - +
- + - - +
- + ]
1252 [ - + + - : 21 : ABI_MALLOC(QR_orthonormal_basis ,(npw_k,lsolutions_max))
- + + - -
+ - + ]
1253 : :
1254 : :
1255 : 7 : omega_is_imaginary = .true.
1256 : :
1257 : :
1258 [ - + + - : 21 : ABI_MALLOC(svd_matrix,(npw_k,lsolutions_max))
- + + - -
+ - + ]
1259 [ - + + - : 21 : ABI_MALLOC(svd_values,(lsolutions_max))
+ - + - -
+ - + ]
1260 : :
1261 : :
1262 : : ! Prepare files for writing
1263 [ - + - - ]: 7 : write_debug = debug .and. mpi_enreg%me == 0
1264 : :
1265 : : if ( write_debug ) then
1266 : :
1267 : 0 : write(filename,'(A)') "ProjectedSternheimerEpsilon.log"
1268 : 0 : inquire(file=filename,exist=file_exists)
1269 : :
1270 : 0 : i = 0
1271 [ # # ]: 0 : do while (file_exists)
1272 : 0 : i = i+1
1273 : 0 : write (filename,'(A,I0.4,A)') "ProjectedSternheimerEpsilon_",i,".log"
1274 : 0 : inquire(file=filename,exist=file_exists)
1275 : : end do
1276 : :
1277 : :
1278 : 0 : io_unit = get_unit()
1279 : 0 : open(io_unit,file=filename,status=files_status_new)
1280 : 0 : write(io_unit,10) ''
1281 : 0 : write(io_unit,10) '#===================================================================================================='
1282 : 0 : write(io_unit,10) "# Building the dielectic matrix using projected Sternheimer equation "
1283 : 0 : write(io_unit,10) "# ------------------------------------------------------------------- "
1284 : 0 : write(io_unit,10) "# "
1285 : 0 : write(io_unit,10) '# This file contains some tests to check if the various elements entering the projected '
1286 : 0 : write(io_unit,10) '# dielectric matrix have the right properties. At this point, this is mostly for debugging.'
1287 : 0 : write(io_unit,10) '# The wavefunctions and other related arrays are stored in reciprocal space.'
1288 : 0 : write(io_unit,10) '#'
1289 : 0 : write(io_unit,10) '#===================================================================================================='
1290 : 0 : write(io_unit,10) ''
1291 : 0 : flush(io_unit)
1292 : : end if
1293 : :
1294 [ - + ]: 7 : if (debug) then
1295 [ # # # # : 0 : ABI_MALLOC(check_matrix ,(lsolutions_max,lsolutions_max))
# # # # #
# # # # #
# # ]
1296 [ # # # # : 0 : ABI_MALLOC(check_matrix2,(lsolutions_max,lsolutions_max))
# # # # #
# # # ]
1297 : : end if
1298 : :
1299 : : !================================================================================
1300 : : ! Loop on all valence bands
1301 : : !================================================================================
1302 : :
1303 : :
1304 : :
1305 : :
1306 [ + + ]: 35 : do v = 1, nbandv
1307 : 28 : write(io_unit_log,10) '#===================================================================================================='
1308 : 28 : write(io_unit_log,20) '# valence band index:', v
1309 : 28 : write(io_unit_log,10) '#===================================================================================================='
1310 : 28 : flush(io_unit_log)
1311 : :
1312 : :
1313 : :
1314 : :
1315 [ - + ]: 28 : if ( write_debug ) then
1316 : 0 : write(io_unit,10) '#===================================================================================================='
1317 : 0 : write(io_unit,20) '# valence band index:', v
1318 : 0 : write(io_unit,10) '#===================================================================================================='
1319 : 0 : flush(io_unit)
1320 : : end if
1321 : :
1322 : 28 : write(io_unit_log,10) ' - Fourier transform valence state ...'
1323 : 28 : flush(io_unit_log)
1324 : :
1325 : :
1326 : : ! copy pre-calculated valence state in this covenient local array
1327 [ + + + + ]: 15448 : psikg_valence(:,:) = kernel_wavefunctions_FFT(:,:,v)
1328 : :
1329 : : ! compute fourier transform of valence state, and conjugate
1330 : 28 : call g_to_r(psir_valence,psikg_valence)
1331 [ + + + + : 192052 : psir_valence(2,:,:,:) = -psir_valence(2,:,:,:)
+ + ]
1332 : :
1333 : : ! loop on all blocks of lanczos vectors
1334 : 28 : write(io_unit_log,10) ' - Loop on all lanczos blocks to generate modified basis and Sternheimer RHS:'
1335 : 28 : flush(io_unit_log)
1336 [ + + ]: 188 : do iblk_lanczos = 1, nbdblock_lanczos
1337 : : !--------------------------------------------------------------------------
1338 : : ! Below, we build the modified basis, [ (V^{1/2} l)^* . psi_v],
1339 : : ! as well as conjugated, projected form, Pc . [ (V^{1/2} l) . psi_v^*].
1340 : : !
1341 : : ! It is very irritating to have to do it this way, but I don't
1342 : : ! see an alternative; see discussion below.
1343 : : !
1344 : : !--------------------------------------------------------------------------
1345 : :
1346 : 160 : write(io_unit_log,23) ' iblk_lanczos = ',iblk_lanczos," / ",nbdblock_lanczos
1347 : :
1348 : :
1349 : 160 : write(io_unit_log,10) ' -- Prepare modified basis computation...'
1350 : 160 : flush(io_unit_log)
1351 : :
1352 : :
1353 : :
1354 : : ! loop on all states within this block
1355 [ + + ]: 384 : do mb = 1, blocksize
1356 : :
1357 : : ! Determine the index of the Lanczos vector
1358 : 224 : l = (iblk_lanczos-1)*blocksize + mb
1359 : :
1360 : : ! take a single lanczos vector
1361 : :
1362 [ + - ]: 224 : if ( l <= lmax) then
1363 [ + + ]: 33120 : psik_wrk(1,:) = dble (Lbasis_lanczos(:,l))
1364 [ + + ]: 33120 : psik_wrk(2,:) = dimag(Lbasis_lanczos(:,l))
1365 : : else
1366 [ # # # # ]: 0 : psik_wrk(:,:) = zero
1367 : : end if
1368 : :
1369 : :
1370 : : ! Apply coulomb potential
1371 : 224 : call sqrt_vc_k(psik_wrk)
1372 : :
1373 : : ! Store in array of blocks of wavefunctions
1374 [ + + + + ]: 99072 : psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:)
1375 : :
1376 : : end do ! mb
1377 : :
1378 : : ! Transform to FFT representation
1379 : 160 : call wf_block_distribute(psikb_wrk, psikg_VL,1) ! LA -> FFT
1380 : :
1381 : : ! psikg_VL now contains | V^1/2 . l >, in FFT configuration
1382 : :
1383 : :
1384 : : !----------------------------------------------------------
1385 : : ! i) Compute the modified basis
1386 : : !----------------------------------------------------------
1387 : 160 : write(io_unit_log,10) ' -- compute modified basis ...'
1388 : 160 : flush(io_unit_log)
1389 : :
1390 : :
1391 : :
1392 : : ! Fourier transform to real space, and conjugate (psir1 is a global work array)
1393 : 160 : call g_to_r(psir1,psikg_VL)
1394 [ + + + + : 1097440 : psir1(2,:,:,:) = -psir1(2,:,:,:) ! IS THIS STACK-DANGEROUS?
+ + ]
1395 : :
1396 : : ! Compute the real space product, and return to k space, in FFT configuration
1397 : 160 : call gr_to_g(psikg_wrk,psir1,psikg_valence)
1398 : :
1399 : : ! psikg_wrk contains | (V^1/2 . l)^* phi_v >, in FFT configuration
1400 : :
1401 : : ! return to LA representation
1402 : 160 : call wf_block_distribute(psikb_wrk, psikg_wrk,2) ! FFT -> LA
1403 : :
1404 : : ! store data, in LA representation
1405 [ + + ]: 384 : do mb = 1, blocksize
1406 : 224 : l = (iblk_lanczos-1)*blocksize + mb
1407 : :
1408 [ + - ]: 384 : if ( l <= lmax) then
1409 : : ! local_Lbasis
1410 [ + + + + ]: 98912 : psik_wrk(:,:) = psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k)
1411 [ + + ]: 33120 : local_Lbasis(:,l) = cmplx_1*psik_wrk(1,:)+cmplx_i*psik_wrk(2,:)
1412 : : end if
1413 : :
1414 : : end do !mb
1415 : :
1416 : : !--------------------------------------------------------------------------
1417 : : ! ii) Build the Sternheimer coefficients, at various frequencies,
1418 : : ! to define the Sternheimer basis.
1419 : : !--------------------------------------------------------------------------
1420 : 160 : write(io_unit_log,20) ' -- Compute Sternheimer RHS...'
1421 : 160 : flush(io_unit_log)
1422 : :
1423 : :
1424 : : ! psikg_wrk still contains | (V^1/2 . l)^* phi_v >, in FFT configuration
1425 : :
1426 : : ! Create right-hand-side of Sternheimer equation, in FFT configuration
1427 : 160 : call pc_k_valence_kernel(psikg_wrk)
1428 : 160 : call Hpsik(psikg_in,psikg_wrk,eig(v))
1429 : 160 : call pc_k_valence_kernel(psikg_in)
1430 [ + + + + ]: 98848 : psikg_in(:,:) = -psikg_in(:,:) ! IS THIS STACK-DANGEROUS?
1431 : :
1432 : : ! return RHS to LA representation, for explicit storage
1433 : 160 : call wf_block_distribute(psikb_wrk, psikg_in,2) ! FFT -> LA
1434 : :
1435 : : ! store data, in LA representation
1436 [ + + ]: 384 : do mb = 1, blocksize
1437 : 224 : l = (iblk_lanczos-1)*blocksize + mb
1438 : :
1439 [ + - ]: 384 : if ( l <= lmax) then
1440 : : ! psi_rhs
1441 [ + + + + ]: 98912 : psi_rhs(:,:,l) = psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k)
1442 : : end if
1443 : :
1444 : : end do !mb
1445 : :
1446 : : !----------------------------------------------------------
1447 : : ! iii) extract solutions for all projection frequencies
1448 : : !----------------------------------------------------------
1449 : 160 : write(io_unit_log,20) ' -- Extract solutions for all projection frequencies...'
1450 : 160 : flush(io_unit_log)
1451 : :
1452 : :
1453 : :
1454 [ - + ]: 160 : do iw = 1, nfrequencies
1455 : :
1456 : 0 : omega0 = list_projection_frequencies(iw)
1457 : : ! Solve Sternheimer equation
1458 : :
1459 : 0 : projection = 0
1460 [ # # ]: 0 : if(omega0 < 1d-12) projection=1
1461 : :
1462 : : ! solve A x = b, over the whole lanczos block
1463 : 0 : call sqmr(psikg_in, psikg_out, eig(v), projection, omega0, omega_is_imaginary)
1464 : :
1465 : : ! return LA representation, for explicit storage
1466 : 0 : call wf_block_distribute(psikb_wrk, psikg_out, 2) ! FFT -> LA
1467 : :
1468 [ - - ]: 160 : do mb = 1, blocksize
1469 : 0 : l = (iblk_lanczos-1)*blocksize + mb
1470 : :
1471 [ # # ]: 0 : if ( l <= lmax) then
1472 : 0 : ls = (l-1)*nfrequencies+iw
1473 : :
1474 [ # # # # ]: 0 : psik_wrk(:,:) = psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k)
1475 : :
1476 [ # # ]: 0 : c_sternheimer_solutions(:,ls)= cmplx_1*psik_wrk(1,:)+cmplx_i*psik_wrk(2,:)
1477 : : end if
1478 : :
1479 : : end do ! mb
1480 : :
1481 : : end do ! iw
1482 : :
1483 : : !----------------------------------------------------------
1484 : : ! iv) Compute the conjugated, projected modified basis
1485 : : !----------------------------------------------------------
1486 : :
1487 : 160 : write(io_unit_log,20) ' -- compute the conjugated, projected modified basis...'
1488 : 160 : flush(io_unit_log)
1489 : :
1490 : :
1491 : :
1492 : :
1493 : : ! Compute the real space product, | (V^1/2. l) phi_v^* > and return to k space, in FFT configuration
1494 : 160 : call gr_to_g(psikg_wrk, psir_valence, psikg_VL)
1495 : :
1496 : : ! project on conduction states
1497 : 160 : call pc_k_valence_kernel(psikg_wrk)
1498 : :
1499 : : ! return to LA representation
1500 : 160 : call wf_block_distribute(psikb_wrk, psikg_wrk,2) ! FFT -> LA
1501 : :
1502 : : ! store back, in LA configuration
1503 [ + + ]: 412 : do mb = 1, blocksize
1504 : 224 : l = (iblk_lanczos-1)*blocksize + mb
1505 : :
1506 [ + - ]: 384 : if ( l <= lmax) then
1507 [ + + + + ]: 98912 : psik_wrk(:,:) = psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k)
1508 [ + + ]: 33120 : local_Lbasis_conjugated(:,l) = cmplx_1*psik_wrk(1,:)+cmplx_i*psik_wrk(2,:)
1509 : : end if
1510 : :
1511 : : end do !mb
1512 : :
1513 : : end do ! iblk_lanczos
1514 : :
1515 : :
1516 : :
1517 : 28 : write(io_unit_log,10) ' - Read in solutions at w=0 and/or w = infinity, if appropriate...'
1518 : 28 : flush(io_unit_log)
1519 : :
1520 : : ! Begin with the storage of the solutions at $\omega = 0.$, which are free.
1521 [ + + ]: 28 : if(dtset%gwls_recycle == 1) then
1522 : : c_sternheimer_solutions(:,lsolutions_max-2*lmax+1:lsolutions_max-lmax) = cmplx_1*Sternheimer_solutions_zero(1,:,:,v) + &
1523 [ + + + + ]: 8368 : & cmplx_i*Sternheimer_solutions_zero(2,:,:,v)
1524 : : end if
1525 [ + + ]: 28 : if(dtset%gwls_recycle == 2) then
1526 [ + + ]: 108 : do i=1,lmax
1527 : 96 : recy_i = (i-1)*nbandv + v
1528 : : !BUG : On petrus, NAG 5.3.1 + OpenMPI 1.6.2 cause read(...,rec=i) to read the data written by write(...,rec=i+1).
1529 : 96 : read(recy_unit,rec=recy_i) psik_wrk
1530 [ + + ]: 24780 : c_sternheimer_solutions(:,lsolutions_max-2*lmax+i) = cmplx_1*psik_wrk(1,:) + cmplx_i*psik_wrk(2,:)
1531 : : end do
1532 : : end if
1533 : :
1534 : : ! and then continue with the storage of the vectors on which the Sternheimer solutions will be projected.
1535 [ + + + + ]: 33148 : c_sternheimer_solutions(:,lsolutions_max-lmax+1:lsolutions_max) = cmplx_1*psi_rhs(1,:,:) + cmplx_i*psi_rhs(2,:,:)
1536 : : ! Previously was = local_Lbasis; but analysis in the Lanczos article reveals psi_rhs should be better.
1537 : : ! Furthermore, tests show that, with psi_rhs, silane@1Ha has Sigma_c 0.01mHa away from the result with
1538 : : ! gwls_list_proj_freq 0.0 1.0, in contrast with local_Lbasis, which has Sigma_c 0.3mHa away from the same result.
1539 : :
1540 [ + + ]: 28 : if ( model ) then
1541 : :
1542 : 24 : write(io_unit_log,10) ' - USE MODEL: model = .true., hence compute model model dielectric matrix...'
1543 : 24 : flush(io_unit_log)
1544 : :
1545 : :
1546 : :
1547 : : !--------------------------------------------------------------------------
1548 : : ! Now that we have the modified basis, compute the matrix
1549 : : ! elements of the model dielectric operator
1550 : : !
1551 : : ! CAREFUL!
1552 : : !
1553 : : ! The model is given by
1554 : : !
1555 : : ! P_model(iw) = sum_{v} phi_v(r) P_c.Y(iw).P_c phi_v^*(r')
1556 : : !
1557 : : ! such that
1558 : : !
1559 : : ! <l1 | eps_model(iw) | l2 > = delta_{l1,l2}
1560 : : ! - sum_{v} < (V^{1/2}.l1).phi_v^*| Pc . Y . Pc | (V^{1/2}.l2).phi_v^* >
1561 : : !
1562 : : ! But local_Lbasis defined above corresponds to
1563 : : ! Pc | (V^{1/2} .l )^* phi_v >.
1564 : : !
1565 : : ! This is why we must define local_Lbasis_conjugated, of the form
1566 : : ! Pc | (V^{1/2} .l ) phi_v^* >.
1567 : : !
1568 : : !--------------------------------------------------------------------------
1569 : :
1570 [ + + ]: 288 : do iw = 1, npt_gauss+1
1571 : :
1572 : 264 : call setup_Pk_model(list_omega(iw),second_model_parameter)
1573 : :
1574 : : ! Only build the lower triangular part; the upper triangular part is obtained from the Hermitian conjugate
1575 [ + + ]: 2400 : do l1 = 1, lmax
1576 : :
1577 [ + + ]: 273504 : YL(:) = model_Y_LA(:)*local_Lbasis_conjugated(:,l1)
1578 [ + + ]: 11880 : do l2 = 1, l1
1579 : : model_dielectric_Lanczos_basis(l1,l2,iw) = model_dielectric_Lanczos_basis(l1,l2,iw) &
1580 : 11616 : -complex_vector_product(YL, local_Lbasis_conjugated(:,l2),npw_k)
1581 : :
1582 : : end do
1583 : : end do
1584 : :
1585 : : end do ! iw
1586 : :
1587 : :
1588 : : end if
1589 : :
1590 : :
1591 : : !--------------------------------------------------------------------------
1592 : : ! Check explicitly that solutions satisfy the Sternheimer equations
1593 : : !--------------------------------------------------------------------------
1594 : :
1595 [ - + ]: 28 : if ( debug ) then
1596 : :
1597 [ # # ]: 0 : if (write_debug) then
1598 : 0 : write(io_unit,10) "#--------------------------------------------------------------------------------"
1599 : 0 : write(io_unit,10) "# Check explicitly that solutions satisfy the Sternheimer equation. "
1600 : 0 : write(io_unit,10) "# "
1601 : 0 : write(io_unit,10) "# Define: "
1602 : 0 : write(io_unit,10) "# E_l = || (omega^2+[H-Ev]^2) |phi_l> + Pc.[H-Ev].Pc |(V^1/2.q_l)^*.phi_v > || "
1603 : 0 : write(io_unit,10) "#--------------------------------------------------------------------------------"
1604 : 0 : write(io_unit,10) '# l Im[omega] (Ha) E_l'
1605 : 0 : write(io_unit,10) "#--------------------------------------------------------------------------------"
1606 : 0 : flush(io_unit)
1607 : : end if
1608 : :
1609 : :
1610 [ # # ]: 0 : do iblk_lanczos = 1, nbdblock_lanczos
1611 : :
1612 : : ! loop on all states within this block
1613 [ # # ]: 0 : do mb = 1, blocksize
1614 : :
1615 : 0 : l = (iblk_lanczos-1)*blocksize + mb
1616 : :
1617 : :
1618 [ # # ]: 0 : if ( l <= lmax) then
1619 : : ! psik_wrk = | (V^{1/2} .l )^* phi_v >
1620 [ # # ]: 0 : psik_wrk(1,:) = dble (local_Lbasis(:,l))
1621 [ # # ]: 0 : psik_wrk(2,:) = dimag(local_Lbasis(:,l))
1622 : : else
1623 [ # # # # ]: 0 : psik_wrk(:,:) = zero
1624 : : end if
1625 : :
1626 : : ! Store in array of blocks of wavefunctions
1627 [ # # # # ]: 0 : psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:)
1628 : : end do ! mb
1629 : :
1630 : : ! Transform to FFT representation
1631 : 0 : call wf_block_distribute(psikb_wrk, psikg_wrk, 1) ! LA -> FFT
1632 : :
1633 : : ! Create right-hand-side of Sternheimer equation
1634 : 0 : call pc_k_valence_kernel(psikg_wrk)
1635 : 0 : call Hpsik(psikg_in,psikg_wrk,eig(v))
1636 : 0 : call pc_k_valence_kernel(psikg_in)
1637 [ # # # # ]: 0 : psikg_in(:,:) = -psikg_in(:,:) ! IS THIS STACK-DANGEROUS?
1638 : :
1639 [ # # ]: 0 : do iw = 1, nfrequencies
1640 : : ! loop on all states within this block
1641 [ # # ]: 0 : do mb = 1, blocksize
1642 : :
1643 : 0 : l = (iblk_lanczos-1)*blocksize + mb
1644 : :
1645 : 0 : ls = (l-1)*nfrequencies+iw
1646 : :
1647 [ # # ]: 0 : if ( l <= lmax) then
1648 [ # # ]: 0 : psik_wrk(1,:) = dble (c_sternheimer_solutions(:,ls))
1649 [ # # ]: 0 : psik_wrk(2,:) = dimag(c_sternheimer_solutions(:,ls))
1650 : : else
1651 [ # # # # ]: 0 : psik_wrk(:,:) = zero
1652 : : end if
1653 : :
1654 : : ! Store in array of blocks of wavefunctions
1655 [ # # # # ]: 0 : psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:)
1656 : : end do
1657 : :
1658 : : ! Transform to FFT representation
1659 : 0 : call wf_block_distribute(psikb_wrk, psikg_wrk, 1) ! LA -> FFT
1660 : :
1661 : :
1662 : 0 : omega0 = list_projection_frequencies(iw)
1663 : :
1664 [ # # # # ]: 0 : psikg_out(:,:) = omega0**2*psikg_wrk(:,:)
1665 : :
1666 : 0 : call Hpsik(psikg_wrk,cte=eig(v))
1667 : 0 : call Hpsik(psikg_wrk,cte=eig(v))
1668 : :
1669 : :
1670 [ # # # # ]: 0 : psikg_out(:,:) = psikg_out(:,:) + psikg_wrk(:,:)-psikg_in(:,:)
1671 : : ! psikg_out now contains [ w0^2 + [H-epsilon_v]^2 ] | x > - |RHS>, in FFT configuration.
1672 : :
1673 : : ! bring it back to LA configuration
1674 : :
1675 : : ! Transform to FFT representation
1676 : 0 : call wf_block_distribute(psikb_wrk, psikg_out, 2) ! FFT -> LA
1677 : :
1678 [ # # ]: 0 : do mb = 1, blocksize
1679 : 0 : l = (iblk_lanczos-1)*blocksize + mb
1680 : :
1681 [ # # ]: 0 : if ( l <= lmax) then
1682 [ # # # # ]: 0 : psik_wrk(:,:) = psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k)
1683 : :
1684 : 0 : z(:) = cg_zdotc(npw_k ,psik_wrk, psik_wrk)
1685 : :
1686 : 0 : call xmpi_sum(z, mpi_communicator,ierr) ! sum on all processors working on FFT!
1687 [ # # ]: 0 : if (write_debug) write(io_unit,21) l, omega0, sqrt(z(1))
1688 : : end if
1689 : : end do ! mb
1690 : :
1691 : : end do ! iw
1692 : : end do ! iblk_lanczos
1693 : :
1694 : : end if
1695 : :
1696 : : !--------------------------------------------------------------------------
1697 : : ! Step 5: Perform a singular value decomposition to extract a
1698 : : ! linearly independent basis for the solution space.
1699 : : !
1700 : : !--------------------------------------------------------------------------
1701 : 28 : write(io_unit_log,10) ' - Perform SVD to extract linearly independent basis to Sternheimer equation...'
1702 : 28 : flush(io_unit_log)
1703 : :
1704 : :
1705 : :
1706 : :
1707 [ + + + + ]: 66268 : svd_matrix(:,:) = c_sternheimer_solutions(:,:)
1708 : :
1709 : 28 : call extract_SVD(mpi_communicator, npw_k,lsolutions_max,svd_matrix,svd_values)
1710 : :
1711 [ - + ]: 28 : if ( write_debug ) then
1712 : 0 : write(io_unit,10) "#--------------------------------------------------------------------------------"
1713 : 0 : write(io_unit,10) "# Check the singular value decomposition of the arrays"
1714 : 0 : write(io_unit,10) '# l svd'
1715 : 0 : write(io_unit,10) "#--------------------------------------------------------------------------------"
1716 : 0 : flush(io_unit)
1717 : : end if
1718 : :
1719 : 28 : lsolutions = 0
1720 [ + + + + ]: 66268 : QR_orthonormal_basis(:,:) = cmplx_0
1721 [ + + ]: 476 : do l=1, lsolutions_max
1722 [ + - ]: 476 : if (svd_values(l) > svd_tolerance ) then
1723 : 448 : lsolutions = lsolutions + 1
1724 : :
1725 [ - + ]: 448 : if ( write_debug ) then
1726 : 0 : write(io_unit,14) l,svd_values(l)
1727 : 0 : flush(io_unit)
1728 : : end if
1729 [ + + ]: 66240 : QR_orthonormal_basis(:,l) = svd_matrix(:,l)
1730 : :
1731 : : else
1732 [ # # ]: 0 : if ( write_debug ) then
1733 : 0 : write(io_unit,15) l,svd_values(l),' SVD value too small! Vector to be discarded!'
1734 : 0 : flush(io_unit)
1735 : : end if
1736 : : end if
1737 : : end do
1738 : :
1739 : : !--------------------------------------------------------------------------
1740 : : ! Step 6: project all relevant arrays onto the newly defined orthonormal
1741 : : ! basis.
1742 : : !--------------------------------------------------------------------------
1743 : 28 : write(io_unit_log,10) ' - Compute the B matrix...'
1744 : 28 : flush(io_unit_log)
1745 : :
1746 [ - + ]: 28 : if (debug) then
1747 [ # # # # ]: 0 : check_matrix(:,:) = cmplx_0
1748 [ # # ]: 0 : do l = 1, lsolutions
1749 : 0 : check_matrix(l,l) = -cmplx_1
1750 : : end do
1751 : : end if
1752 : :
1753 [ - + + - : 84 : ABI_MALLOC(ipiv ,(lsolutions))
+ - + - -
+ - + ]
1754 [ - + + - : 112 : ABI_MALLOC(sternheimer_A0 ,(lsolutions,lsolutions))
- + + - +
- + - - +
- + ]
1755 [ - + + - : 112 : ABI_MALLOC(sternheimer_B ,(lsolutions,lmax))
- + + - +
- + - - +
- + ]
1756 [ - + + - : 112 : ABI_MALLOC(sternheimer_G ,(lmax,lsolutions))
- + + - +
- + - - +
- + ]
1757 : :
1758 : : ! Compute the X matrix and the check_matrix
1759 [ + + ]: 476 : do l1 = 1, lsolutions
1760 : :
1761 [ + + ]: 66240 : psi_gamma_l1(1,:) = real (QR_orthonormal_basis(:,l1))
1762 [ + + ]: 66240 : psi_gamma_l1(2,:) = dimag(QR_orthonormal_basis(:,l1))
1763 : :
1764 [ + + ]: 7644 : do l2 = 1, lsolutions
1765 : :
1766 [ + + ]: 7168 : if (l2 <= lmax) then
1767 : 3584 : z(:) = cg_zdotc(npw_k, psi_gamma_l1, psi_rhs(:,:,l2))
1768 : 3584 : call xmpi_sum(z, mpi_communicator,ierr) ! sum on all processors for LA configuration
1769 : :
1770 : 3584 : sternheimer_B(l1,l2) = cmplx_1*z(1)+cmplx_i*z(2)
1771 : : end if
1772 : :
1773 [ - + ]: 7616 : if (debug) then
1774 [ # # ]: 0 : psi_gamma_l2(1,:) = dble (QR_orthonormal_basis(:,l2))
1775 [ # # ]: 0 : psi_gamma_l2(2,:) = dimag(QR_orthonormal_basis(:,l2))
1776 : :
1777 : 0 : z(:) = cg_zdotc(npw_k, psi_gamma_l1, psi_gamma_l2)
1778 : 0 : call xmpi_sum(z, mpi_communicator,ierr) ! sum on all processors
1779 : :
1780 : 0 : check_matrix(l1,l2) = check_matrix(l1,l2) + cmplx_1*z(1)+cmplx_i*z(2)
1781 : : end if
1782 : :
1783 : : end do ! l2
1784 : : end do ! l1
1785 : :
1786 : :
1787 : : ! Number of blocks of solution vectors
1788 : 28 : nbdblock_solutions = lsolutions/blocksize
1789 : :
1790 [ - + ]: 28 : if (modulo(lsolutions,blocksize) /= 0) nbdblock_solutions = nbdblock_solutions + 1
1791 : :
1792 : :
1793 : 28 : write(io_unit_log,10) ' - Compute the A0 matrix...'
1794 : 28 : flush(io_unit_log)
1795 : :
1796 : : ! Compute the A matrix
1797 [ + + ]: 348 : do iblk_solutions =1, nbdblock_solutions
1798 : :
1799 [ + + ]: 768 : do mb = 1, blocksize
1800 : 448 : l2 = (iblk_solutions-1)*blocksize + mb
1801 : :
1802 [ + - ]: 448 : if ( l2 <= lsolutions) then
1803 [ + + ]: 66240 : psik_wrk(1,:) = dble (QR_orthonormal_basis(:,l2))
1804 [ + + ]: 66240 : psik_wrk(2,:) = dimag(QR_orthonormal_basis(:,l2))
1805 : : else
1806 [ # # # # ]: 0 : psik_wrk(:,:) = zero
1807 : : end if
1808 : :
1809 : : ! Store in array of blocks of wavefunctions
1810 [ + + + + ]: 198144 : psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k) = psik_wrk(:,:)
1811 : : end do
1812 : :
1813 : : ! Transform to FFT representation
1814 : 320 : call wf_block_distribute(psikb_wrk, psikg_wrk,1) ! LA -> FFT
1815 : :
1816 : : ! act twice with the Hamiltonian operator
1817 : 320 : call Hpsik(psikg_out,psikg_wrk,eig(v))
1818 : 320 : call Hpsik(psikg_out,cte=eig(v))
1819 : :
1820 : : ! return to LA representation
1821 : 320 : call wf_block_distribute(psikb_wrk, psikg_out,2) ! FFT -> LA
1822 : :
1823 [ + + ]: 796 : do mb = 1, blocksize
1824 : 448 : l2 = (iblk_solutions-1)*blocksize + mb
1825 : :
1826 [ + - ]: 448 : if ( l2 <= lsolutions) then
1827 [ + + + + ]: 197824 : psik_wrk(:,:) = psikb_wrk(:,(mb-1)*npw_k+1:mb*npw_k)
1828 : : else
1829 [ # # # # ]: 0 : psik_wrk(:,:) = zero
1830 : : end if
1831 : :
1832 [ + + ]: 7936 : do l1 = 1, lsolutions
1833 : :
1834 [ + + ]: 1059840 : psi_gamma_l1(1,:) = real (QR_orthonormal_basis(:,l1))
1835 [ + + ]: 1059840 : psi_gamma_l1(2,:) = dimag(QR_orthonormal_basis(:,l1))
1836 : :
1837 : 7168 : z(:) = cg_zdotc(npw_k, psi_gamma_l1, psik_wrk)
1838 : 7168 : call xmpi_sum(z,mpi_communicator,ierr) ! sum on all processors working on FFT!
1839 : :
1840 : :
1841 [ + - ]: 7616 : if ( l2 <= lsolutions ) then
1842 : 7168 : sternheimer_A0(l1,l2) = cmplx_1*z(1)+cmplx_i*z(2)
1843 : : end if
1844 : :
1845 : :
1846 : : end do ! l1
1847 : : end do ! mb
1848 : : end do ! iblk_solutions
1849 : :
1850 : :
1851 [ - + ]: 28 : if (debug) then
1852 : : ! HERE we use a dummy variable to avoid operations which might blow the stack!
1853 : : ! Stack overflows lead to hard-to-find bugs; let's avoid putting them in here.
1854 [ # # # # : 0 : ABI_MALLOC(dummy_tmp_1,(lsolutions,lsolutions))
# # # # #
# # # ]
1855 [ # # # # : 0 : ABI_MALLOC(dummy_tmp_2,(lsolutions,lsolutions))
# # # # #
# # # ]
1856 : :
1857 : :
1858 [ # # # # ]: 0 : dummy_tmp_1(:,:)= transpose(sternheimer_A0(:,:))
1859 [ # # # # ]: 0 : dummy_tmp_2(:,:)= conjg(dummy_tmp_1(:,:))
1860 : :
1861 [ # # # # ]: 0 : check_matrix2(:,:) = sternheimer_A0(:,:)-dummy_tmp_2(:,:)
1862 : :
1863 [ # # ]: 0 : ABI_FREE(dummy_tmp_1)
1864 [ # # ]: 0 : ABI_FREE(dummy_tmp_2)
1865 : : end if
1866 : :
1867 : 28 : write(io_unit_log,10) ' - Compute the GAMMA matrix...'
1868 : 28 : flush(io_unit_log)
1869 : :
1870 : :
1871 : : ! Compute the GAMMA matrices
1872 [ + + ]: 252 : do l1 = 1, lmax
1873 : :
1874 : : ! psik_wrk = | (V^{1/2} . l )^* phi_v >
1875 [ + + ]: 33120 : psik_wrk(1,:) = dble (local_Lbasis(:,l1) )
1876 [ + + ]: 33120 : psik_wrk(2,:) = dimag(local_Lbasis(:,l1))
1877 : :
1878 : :
1879 [ + + ]: 3836 : do l2 = 1, lsolutions
1880 : :
1881 [ + + ]: 529920 : psi_gamma_l2(1,:) = real (QR_orthonormal_basis(:,l2))
1882 [ + + ]: 529920 : psi_gamma_l2(2,:) = dimag(QR_orthonormal_basis(:,l2))
1883 : :
1884 : :
1885 : : ! Note that G_{lJ} = < l | Vc^{1/2}.(gamma_J^*.phi_v)>
1886 : : ! = < gamma_J | (Vc^{1/2}.l^*).phi_v >
1887 : :
1888 : 3584 : z(:) = cg_zdotc(npw_k,psi_gamma_l2, psik_wrk)
1889 : 3584 : call xmpi_sum(z,mpi_communicator,ierr) ! sum on all processors working on FFT!
1890 : 3808 : sternheimer_G(l1,l2) = cmplx_1*z(1)+cmplx_i*z(2)
1891 : :
1892 : : end do ! l1
1893 : : end do ! l2
1894 : :
1895 : :
1896 : :
1897 [ - + ]: 28 : if ( write_debug ) then
1898 [ # # # # ]: 0 : write(io_unit,19) '<gamma| gamma>', sqrt(sum(abs(check_matrix(:,:))**2))
1899 [ # # # # ]: 0 : write(io_unit,19) ' M hermitian ',sqrt(sum(abs(check_matrix2(:,:))**2))
1900 : 0 : write(io_unit,10) " "
1901 : 0 : write(io_unit,10) "# GAMMA Matrix:"
1902 : 0 : write(io_unit,10) " "
1903 : :
1904 [ # # ]: 0 : do l1 = 1, lmax
1905 : 0 : write(io_unit,30) sternheimer_G(l1,:)
1906 : : end do
1907 : 0 : flush(io_unit)
1908 : : end if
1909 : :
1910 : : !--------------------------------------------------------------------------
1911 : : ! Step 7: Compute the solutions
1912 : : !--------------------------------------------------------------------------
1913 : :
1914 : 28 : write(io_unit_log,10) ' - Compute the Projected Sternheimer solutions and build approximate dielectric operator...'
1915 : 28 : flush(io_unit_log)
1916 : :
1917 : :
1918 : :
1919 [ - + + - : 84 : ABI_MALLOC(sternheimer_A ,(lsolutions,lsolutions))
- + + - -
+ - + ]
1920 [ - + + - : 84 : ABI_MALLOC(sternheimer_X ,(lsolutions,lmax))
- + + - -
+ - + ]
1921 [ + + ]: 308 : do iw = 2, npt_gauss + 1
1922 : :
1923 : 280 : write(io_unit_log,23) ' -- iw = ',iw,' / ',npt_gauss+1
1924 : 280 : flush(io_unit_log)
1925 : :
1926 : :
1927 : 280 : omega = list_omega(iw)
1928 : :
1929 [ + + + + ]: 76440 : sternheimer_A(:,:) = sternheimer_A0(:,:)
1930 : :
1931 [ + + ]: 4760 : do l = 1, lsolutions
1932 : 4760 : sternheimer_A(l,l) = sternheimer_A(l,l) + omega**2
1933 : : end do
1934 : :
1935 [ + + + + ]: 38360 : sternheimer_X(:,:) = sternheimer_B(:,:)
1936 : :
1937 : : !--------------------------------------------------------------------------
1938 : : ! Step 2: solve A*X = B, a projected form of the Sternheimer equation
1939 : : !--------------------------------------------------------------------------
1940 : 280 : write(io_unit_log,10) ' -- Solve A*X = B'
1941 : 280 : flush(io_unit_log)
1942 : :
1943 : :
1944 : :
1945 : 280 : call cpu_time(time1)
1946 : : call zgesv(lsolutions, & ! number of rows of A matrix
1947 : : lmax, & ! number of columns of B matrix
1948 : : sternheimer_A, & ! The A matrix on input, the LU factorization on output
1949 : : lsolutions, & ! leading dimension of A
1950 : : ipiv, & ! array of pivots
1951 : : sternheimer_X, & ! B matrix on input, solution X on output
1952 : : lsolutions, & ! leading dimension of B
1953 : 280 : info )
1954 : 280 : call cpu_time(time2)
1955 : :
1956 : 280 : time_exact = time_exact + time2-time1
1957 : :
1958 : : !--------------------------------------------------------------------------
1959 : : ! Step 3: Add contribution to projected epsilon
1960 : : !--------------------------------------------------------------------------
1961 : : ! perform E = E -4 sum_l Gamma_v*X^*_v
1962 : :
1963 : 280 : write(io_unit_log,10) ' -- add contribution to dielectric matrix at this frequency'
1964 : 280 : flush(io_unit_log)
1965 : :
1966 : :
1967 : :
1968 [ - + + - : 1120 : ABI_MALLOC(dummy_tmp_1,(lsolutions,lmax))
- + + - +
- + - - +
- + ]
1969 [ + + + + ]: 38360 : dummy_tmp_1(:,:) = conjg(sternheimer_X) ! DO THIS to avoid potential stack problems
1970 : :
1971 : 280 : call cpu_time(time1)
1972 : : call zgemm( 'N', & ! A matrix is in normal order
1973 : : 'N', & ! B matrix is in normal order
1974 : : lmax, & ! number of rows of A
1975 : : lmax, & ! number of columns of B
1976 : : lsolutions, & ! number of columns of A
1977 : : -4.0_dp*cmplx_1, & ! premultiply A*B by this scalar
1978 : : sternheimer_G, & ! GAMMA matrix
1979 : : lmax, & ! leading dimension of A
1980 : : dummy_tmp_1, & ! B matrix
1981 : : lsolutions, & ! leading dimension of B
1982 : : cmplx_1, & ! beta is one
1983 : : projected_dielectric_Lanczos_basis(:,:,iw), & ! C matrix
1984 : 280 : lmax) ! leading dimension of C
1985 : :
1986 [ - + ]: 280 : ABI_FREE(dummy_tmp_1)
1987 : :
1988 : 280 : call cpu_time(time2)
1989 : 308 : time_exact = time_exact + time2-time1
1990 : :
1991 : : end do ! iw
1992 : :
1993 : 28 : write(io_unit_log,10) ' - Deallocate tmp arrays...'
1994 : 28 : flush(io_unit_log)
1995 : :
1996 : :
1997 : :
1998 : :
1999 [ - + ]: 28 : ABI_FREE(ipiv )
2000 [ - + ]: 28 : ABI_FREE(sternheimer_A )
2001 [ - + ]: 28 : ABI_FREE(sternheimer_A0 )
2002 [ - + ]: 28 : ABI_FREE(sternheimer_B )
2003 [ - + ]: 28 : ABI_FREE(sternheimer_X )
2004 [ - + ]: 35 : ABI_FREE(sternheimer_G )
2005 : :
2006 : : end do ! v
2007 : :
2008 : : !--------------------------------------------------------------------------
2009 : : ! Finalize, post v loop
2010 : : !--------------------------------------------------------------------------
2011 : 7 : write(io_unit_log,10) " - Finalize, after band iterations...."
2012 : 7 : flush(io_unit_log)
2013 : :
2014 : :
2015 : :
2016 : :
2017 : :
2018 : 7 : timing_string = "# Exact Sector : "
2019 : 7 : call write_timing_log(timing_string,time_exact)
2020 : :
2021 : :
2022 [ - + + - : 28 : ABI_MALLOC(dummy_tmp_1,(lmax,lmax))
- + + - +
- + - - +
- + ]
2023 [ - + + - : 21 : ABI_MALLOC(dummy_tmp_2,(lmax,lmax))
- + + - -
+ - + ]
2024 : :
2025 [ + + ]: 77 : do iw = 2, npt_gauss+1
2026 : : ! finally, make sure matrix is hermitian
2027 : :
2028 : : ! play this little game to avoid STACK problems
2029 [ + + + + ]: 5110 : dummy_tmp_1(:,:) = 0.5_dp*transpose(projected_dielectric_Lanczos_basis(:,:,iw))
2030 [ + + + + ]: 5110 : dummy_tmp_2(:,:) = conjg(dummy_tmp_1(:,:))
2031 [ + + + + ]: 5110 : dummy_tmp_1(:,:) = dummy_tmp_2(:,:) +0.5_dp*projected_dielectric_Lanczos_basis(:,:,iw)
2032 : :
2033 [ + + + + ]: 5117 : projected_dielectric_Lanczos_basis(:,:,iw) = dummy_tmp_1(:,:)
2034 : : end do
2035 : :
2036 [ - + ]: 7 : ABI_FREE(dummy_tmp_1)
2037 [ - + ]: 7 : ABI_FREE(dummy_tmp_2)
2038 : :
2039 : :
2040 : :
2041 [ - + ]: 7 : if ( write_debug ) then
2042 : : !write some results to a file
2043 : :
2044 : 0 : write(io_unit,10) '#===================================================================================================='
2045 : 0 : write(io_unit,10) "# Projected dielectric matrices "
2046 : 0 : write(io_unit,10) "# ------------------------------------------------------- "
2047 : 0 : write(io_unit,10) "# This file contains the various projected dielectric matrices as a function of frequency. "
2048 : 0 : write(io_unit,10) '#===================================================================================================='
2049 : 0 : write(io_unit,10) ''
2050 : 0 : flush(io_unit)
2051 : :
2052 [ # # ]: 0 : do iw = 1, npt_gauss+1
2053 : 0 : write(io_unit,10) "#"
2054 : 0 : write(io_unit,12) "# omega = ",list_omega(iw), " i Ha"
2055 : 0 : write(io_unit,10) "#"
2056 : :
2057 [ # # ]: 0 : do l =1, lmax
2058 : 0 : write(io_unit,30) projected_dielectric_Lanczos_basis(l,:,iw)
2059 : : end do
2060 : : end do
2061 : :
2062 : : end if
2063 : :
2064 : :
2065 [ + + ]: 7 : if ( model ) then
2066 : :
2067 : :
2068 : : ! Add all components on the processors
2069 : 6 : call xmpi_sum(model_dielectric_Lanczos_basis,mpi_communicator,ierr) ! sum on all processors
2070 : :
2071 : : ! add the identity
2072 [ + + ]: 72 : do iw = 1, npt_gauss+1
2073 [ + + ]: 600 : do l= 1, lmax
2074 : 594 : model_dielectric_Lanczos_basis(l,l,iw) = model_dielectric_Lanczos_basis(l,l,iw) + cmplx_1
2075 : : end do
2076 : : end do
2077 : :
2078 : : ! hermitian the operator
2079 [ + + ]: 72 : do iw = 1, npt_gauss+1
2080 : :
2081 [ + + ]: 600 : do l1 = 1, lmax
2082 [ + + ]: 2970 : do l2 = 1, l1
2083 : : ! operator is hermitian
2084 : 2904 : model_dielectric_Lanczos_basis(l2,l1,iw) = conjg(model_dielectric_Lanczos_basis(l1,l2,iw))
2085 : : end do
2086 : : end do
2087 : :
2088 : : end do ! iw
2089 : :
2090 : : end if
2091 : :
2092 : :
2093 : :
2094 [ - + ]: 7 : if ( write_debug ) then
2095 : 0 : close(io_unit)
2096 : : end if
2097 : :
2098 : 7 : write(io_unit_log,10) " - Deallocate and exit...."
2099 : 7 : flush(io_unit_log)
2100 : :
2101 : :
2102 : :
2103 [ - + ]: 7 : if (debug) then
2104 [ # # ]: 0 : ABI_FREE(check_matrix)
2105 [ # # ]: 0 : ABI_FREE(check_matrix2)
2106 : : end if
2107 : :
2108 : :
2109 [ - + ]: 7 : ABI_FREE(psikg_valence)
2110 [ - + ]: 7 : ABI_FREE(psir_valence)
2111 [ - + ]: 7 : ABI_FREE(psikg_VL)
2112 : :
2113 : :
2114 [ - + ]: 7 : ABI_FREE(YL)
2115 : :
2116 [ - + ]: 7 : ABI_FREE(local_Lbasis_conjugated)
2117 [ - + ]: 7 : ABI_FREE(local_Lbasis)
2118 : :
2119 : :
2120 [ - + ]: 7 : ABI_FREE(svd_matrix)
2121 [ - + ]: 7 : ABI_FREE(svd_values)
2122 [ - + ]: 7 : ABI_FREE(psi_gamma_l1)
2123 [ - + ]: 7 : ABI_FREE(psi_gamma_l2)
2124 : :
2125 [ - + ]: 7 : ABI_FREE(psikg_in)
2126 [ - + ]: 7 : ABI_FREE(psikg_out)
2127 : :
2128 [ - + ]: 7 : ABI_FREE(psi_rhs)
2129 : :
2130 [ - + ]: 7 : ABI_FREE(c_sternheimer_solutions)
2131 [ - + ]: 7 : ABI_FREE(QR_orthonormal_basis)
2132 : :
2133 [ - + ]: 7 : ABI_FREE(psik_wrk)
2134 [ - + ]: 7 : ABI_FREE(psikb_wrk)
2135 [ - + ]: 7 : ABI_FREE(psikg_wrk)
2136 : :
2137 : :
2138 : 7 : close(io_unit_log)
2139 : :
2140 : :
2141 : : 10 format(A)
2142 : : 12 format(A,F12.8,A)
2143 : : 14 format(I5,10X,ES24.12)
2144 : : 15 format(I5,10X,ES24.12,10X,A)
2145 : :
2146 : : 19 format(20X,A,15X,E24.16)
2147 : :
2148 : : 20 format(A,I5)
2149 : : 21 format(I5,10X,F8.4,15X,ES24.12)
2150 : : 22 format(A,I5,A)
2151 : : 23 format(A,I5,A,I5)
2152 : : 30 format(2X,1000(ES12.4,2X,ES12.4,5X))
2153 : :
2154 [ - + - + ]: 7 : end subroutine ProjectedSternheimerEpsilon
2155 : : !!***
2156 : :
2157 : :
2158 : : end module m_gwls_DielectricArray
2159 : : !!***
|
