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