Line data Source code
1 : !!****m* ABINIT/m_gwls_polarisability
2 : !! NAME
3 : !! m_gwls_polarisability
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_polarisability
24 : ! local modules
25 : use m_gwls_utility
26 : use m_gwls_wf
27 : use m_gwls_valenceWavefunctions
28 : use m_gwls_hamiltonian
29 : use m_gwls_lineqsolver
30 :
31 : ! abinit modules
32 : use defs_basis
33 : use m_errors
34 : use m_abicore
35 : use m_bandfft_kpt
36 :
37 : use m_time, only : timab
38 :
39 : implicit none
40 : save
41 : private
42 : !!***
43 :
44 : real(dp),public :: matrix_function_omega(2)
45 :
46 :
47 : ! Some timing variables
48 : integer, public :: counter_fft = 0, counter_sqmr = 0, counter_rprod = 0 , counter_proj = 0, counter_H = 0
49 :
50 : real(dp), public :: time1 = zero, time2 = zero, time_fft = zero
51 : real(dp), public :: time_sqmr = zero, time_rprod = zero, time_proj = zero, time_H = zero
52 :
53 : real(dp), allocatable, public :: Sternheimer_solutions_zero(:,:,:,:)
54 : integer, public :: index_solution=0
55 : integer, public :: recy_unit
56 : logical, public :: write_solution=.false.
57 :
58 : !integer :: io_unit
59 : !!***
60 :
61 : public :: Pk, epsilon_k
62 : public :: matrix_function_epsilon_k
63 : public :: set_dielectric_function_frequency
64 : !!***
65 :
66 : contains
67 :
68 : !!****f* m_hamiltonian/Pk
69 : !! NAME
70 : !! Pk
71 : !!
72 : !! FUNCTION
73 : !! .
74 : !!
75 : !! INPUTS
76 : !!
77 : !! OUTPUT
78 : !!
79 : !! SOURCE
80 :
81 184 : subroutine Pk(psi_inout,omega)
82 : !===============================================================================
83 : !
84 : ! This routine applies the polarizability operator to an arbitrary state psi_inout.
85 : !
86 : !
87 : ! A note about parallelism:
88 : ! -------------------------
89 : ! The input/output is in "linear algebra" configuration, which is to say
90 : ! that ALL processors have a fraction of the G-vectors for this state. Internally,
91 : ! this routine will parallelise over bands and FFT, thus using the "FFT"
92 : ! configuration of the data. For more information, see the lobpcgwf.F90, and
93 : ! the article by F. Bottin et al.
94 : !===============================================================================
95 :
96 : real(dp), intent(inout) :: psi_inout(2,npw_k)
97 : real(dp), intent(in) :: omega(2)
98 :
99 : logical :: omega_imaginary
100 :
101 : integer :: v, mb, iblk
102 :
103 : real(dp):: norm_omega
104 184 : real(dp), allocatable :: psik(:,:)
105 184 : real(dp), allocatable :: psik_alltoall(:,:), psik_wrk_alltoall(:,:)
106 184 : real(dp), allocatable :: psik_in_alltoall(:,:), psik_tmp_alltoall(:,:)
107 :
108 184 : real(dp), allocatable :: psik_ext(:,:), psik_ext_alltoall(:,:)
109 :
110 :
111 184 : real(dp), allocatable :: psir(:,:,:,:), psir_ext(:,:,:,:)
112 :
113 : integer :: cplex
114 : integer :: recy_i
115 :
116 :
117 : integer :: mpi_band_rank
118 :
119 : real(dp) :: list_SQMR_frequencies(2)
120 : real(dp) :: list_QMR_frequencies(2,2)
121 :
122 :
123 : integer, save :: icounter = 0
124 : real(dp), save :: total_time1 = zero, total_time2 = zero, total_time = zero
125 :
126 :
127 : integer :: num_op_v, i_op_v, case_op_v
128 : real(dp):: factor_op_v
129 : real(dp):: lbda
130 : real(dp):: zz(2)
131 :
132 : character(len=500) :: message
133 :
134 : real(dp) :: tsec(2)
135 : integer :: GWLS_TIMAB, OPTION_TIMAB
136 :
137 : ! *************************************************************************
138 :
139 184 : GWLS_TIMAB = 1524
140 184 : OPTION_TIMAB = 1
141 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
142 :
143 :
144 :
145 184 : icounter = icounter + 1
146 184 : call cpu_time(total_time1)
147 :
148 :
149 : !========================================
150 : ! Allocate work arrays and define
151 : ! important parameters
152 : !========================================
153 184 : GWLS_TIMAB = 1525
154 : OPTION_TIMAB = 1
155 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
156 :
157 :
158 184 : cplex = 2 ! complex potential
159 :
160 184 : mpi_band_rank = mpi_enreg%me_band
161 :
162 552 : ABI_MALLOC(psik, (2,npw_kb))
163 552 : ABI_MALLOC(psik_alltoall, (2,npw_g))
164 368 : ABI_MALLOC(psik_wrk_alltoall, (2,npw_g))
165 368 : ABI_MALLOC(psik_tmp_alltoall, (2,npw_g))
166 368 : ABI_MALLOC(psik_in_alltoall, (2,npw_g))
167 :
168 920 : ABI_MALLOC(psir, (2,n4,n5,n6))
169 736 : ABI_MALLOC(psir_ext,(2,n4,n5,n6))
170 :
171 :
172 184 : OPTION_TIMAB = 2
173 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
174 :
175 :
176 : !--------------------------------------------------------------------------------
177 : !
178 : ! The polarizability acting on a state | PSI > is given by
179 : !
180 : ! Pk | PSI > = 2 sum_{v} | h_v^* phi_v > ( factor of 2 comes from sum on spin)
181 : !
182 : ! | h_v > = Pc . OPERATOR_v . Pc | PSI^* phi_v >
183 : !
184 : ! There are multiple cases to consider:
185 : !
186 : ! CASE:
187 : ! 1) |omega| = 0
188 : ! OPERATOR_v = 1/[H-ev]
189 : ! prefactor = -4
190 : !
191 : ! 2) omega = i lbda (imaginary)
192 : ! OPERATOR_v = [H-ev]/(lbda^2+[H-ev]^2)
193 : ! prefactor = -4
194 : !
195 : ! 3) omega = lbda (real) (USE SQMR)
196 : ! OPERATOR_v = { 1/[H-ev+lbda]+ 1/[H-ev-lbda] }
197 : ! prefactor = -2
198 : !
199 : ! 4) omega = lbda (real) (USE QMR)
200 : ! OPERATOR_v = { 1/[H-ev+lbda]+ 1/[H-ev-lbda] }
201 : ! prefactor = -2
202 : !
203 : ! It simplifies the code below to systematize the algorithm.
204 : !
205 : !--------------------------------------------------------------------------------
206 :
207 : ! Check which part of omega is non-zero. Default is that omega is real.
208 184 : if (abs(omega(2)) < 1.0d-12) then
209 184 : omega_imaginary=.false.
210 184 : norm_omega = omega(1)
211 :
212 0 : elseif (abs(omega(1)) < 1.0d-12 .and. abs(omega(2)) > 1.0d-12) then
213 0 : omega_imaginary = .true.
214 0 : norm_omega = omega(2)
215 :
216 : else
217 : write(message,"(a,es16.8,3a)")&
218 0 : "omega=",omega,",",ch10,&
219 0 : "but either it's real or imaginary part need to be 0 for the polarisability routine to work."
220 0 : ABI_ERROR(message)
221 : end if
222 :
223 :
224 : !-----------------------------------------------------------------
225 : ! I) Prepare global values depending on the CASE being considered
226 : !-----------------------------------------------------------------
227 :
228 184 : if (norm_omega < 1.0D-12 ) then
229 72 : case_op_v = 1
230 72 : num_op_v = 1
231 72 : factor_op_v = -4.0_dp
232 :
233 112 : else if (omega_imaginary) then
234 0 : case_op_v = 2
235 0 : num_op_v = 1
236 0 : factor_op_v = -4.0_dp
237 :
238 112 : else if( .not. activate_inf_shift_poles) then
239 96 : case_op_v = 3
240 96 : num_op_v = 2
241 96 : factor_op_v = -2.0_dp
242 :
243 : else
244 16 : case_op_v = 4
245 16 : num_op_v = 2
246 16 : factor_op_v = -2.0_dp
247 :
248 16 : inf_shift_poles = dtset%zcut
249 :
250 :
251 16 : write(message,*) " inf_shift_poles = ",inf_shift_poles
252 16 : call wrtout(std_out,message,'COLL')
253 : end if
254 :
255 :
256 184 : write(message,10)" "
257 184 : call wrtout(std_out,message,'COLL')
258 184 : write(message,10) " Pk: applying the polarizability on states"
259 184 : call wrtout(std_out,message,'COLL')
260 184 : write(message,10) " =============================================================="
261 184 : call wrtout(std_out,message,'COLL')
262 184 : write(message,12) " CASE : ",case_op_v
263 184 : call wrtout(std_out,message,'COLL')
264 :
265 :
266 :
267 : !-----------------------------------------------------------------
268 : ! II) copy conjugate of initial wavefunction in local array,
269 : ! and set inout array to zero. Each FFT row of processors
270 : ! must have a copy of the initial wavefunction in FFT
271 : ! configuration!
272 : !-----------------------------------------------------------------
273 184 : call cpu_time(time1)
274 :
275 184 : GWLS_TIMAB = 1526
276 184 : OPTION_TIMAB = 1
277 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
278 :
279 :
280 552 : ABI_MALLOC(psik_ext,(2,npw_kb))
281 552 : ABI_MALLOC(psik_ext_alltoall,(2,npw_g))
282 :
283 : ! fill the array psik_ext with copies of the external state
284 464 : do mb = 1, blocksize
285 105320 : psik_ext(:,(mb-1)*npw_k+1:mb*npw_k) = psi_inout(:,:)
286 : end do
287 :
288 : ! change configuration of the data, from LA to FFT
289 184 : call wf_block_distribute(psik_ext, psik_ext_alltoall,1) ! LA -> FFT
290 : ! Now every row of FFT processors has a copy of the external state.
291 :
292 : ! Copy the external state to the real space format, appropriate for real space products to be
293 : ! used later.
294 :
295 184 : call g_to_r(psir_ext,psik_ext_alltoall)
296 1262056 : psir_ext(2,:,:,:) = -psir_ext(2,:,:,:)
297 :
298 :
299 : ! Don't need these arrays anymore...
300 184 : ABI_FREE(psik_ext)
301 184 : ABI_FREE(psik_ext_alltoall)
302 :
303 184 : OPTION_TIMAB = 2
304 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
305 :
306 184 : call cpu_time(time2)
307 184 : time_fft = time_fft + time2-time1
308 184 : counter_fft = counter_fft+1
309 :
310 : ! set external state to zero, ready to start cumulating the answer
311 86536 : psi_inout = zero
312 :
313 : !-----------------------------------------------------------------
314 : ! III) Iterate on all valence bands
315 : !-----------------------------------------------------------------
316 :
317 : ! Loop on all blocks of eigenstates
318 1544 : do iblk = 1, nbdblock
319 :
320 :
321 : ! What is the valence band index for this block and this row of FFT processors? It is not clear from the
322 : ! code in lobpcgwf.F90; I'm going to *guess*.
323 :
324 1360 : v = (iblk-1)*blocksize + mpi_band_rank + 1 ! CAREFUL! This is a guess. Revisit this if code doesn't work as intended.
325 :
326 :
327 : !Solving of Sternheiner equation
328 1360 : write(message,12) " band :", v
329 1360 : call wrtout(std_out,message,'COLL')
330 1360 : write(message,14) " eigenvalue (Ha) : ",eig(v)
331 1360 : call wrtout(std_out,message,'COLL')
332 1360 : write(message,14) " Re[omega] (Ha) : ",omega(1)
333 1360 : call wrtout(std_out,message,'COLL')
334 1360 : write(message,14) " Im[omega] (Ha) : ",omega(2)
335 1360 : call wrtout(std_out,message,'COLL')
336 :
337 : !-----------------------------------------------------------------
338 : ! IV) prepare some arrays, if they are needed
339 : !-----------------------------------------------------------------
340 1360 : if (case_op_v == 3) then
341 :
342 640 : list_SQMR_frequencies(1) = eig(v) - norm_omega
343 640 : list_SQMR_frequencies(2) = eig(v) + norm_omega
344 :
345 720 : else if (case_op_v == 4) then
346 480 : list_QMR_frequencies(:,1) = (/eig(v)-norm_omega,-inf_shift_poles/)
347 480 : list_QMR_frequencies(:,2) = (/eig(v)+norm_omega, inf_shift_poles/)
348 : end if
349 :
350 :
351 : !-----------------------------------------------------------------
352 : ! V) Compute the real-space product of the input wavefunction
353 : ! with the valence wavefunction.
354 : !-----------------------------------------------------------------
355 : ! Unfortunately, the input wavefunctions will be FFT-transformed k-> r nbandv times
356 : ! because fourwf cannot multiply a real-space potential with a real-space wavefunction!
357 1360 : call cpu_time(time1)
358 :
359 1360 : GWLS_TIMAB = 1527
360 1360 : OPTION_TIMAB = 1
361 1360 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
362 :
363 1360 : call gr_to_g(psik_in_alltoall,psir_ext,valence_wavefunctions_FFT(:,:,iblk))
364 :
365 1360 : OPTION_TIMAB = 2
366 1360 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
367 :
368 :
369 1360 : call cpu_time(time2)
370 1360 : time_fft = time_fft + time2-time1
371 1360 : counter_fft = counter_fft+1
372 :
373 : !-----------------------------------------------------------------
374 : ! VI) Project out to conduction space
375 : !-----------------------------------------------------------------
376 1360 : call cpu_time(time1)
377 :
378 1360 : GWLS_TIMAB = 1528
379 1360 : OPTION_TIMAB = 1
380 1360 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
381 :
382 : !call pc_k(psik_in)
383 1360 : call pc_k_valence_kernel(psik_in_alltoall)
384 :
385 1360 : OPTION_TIMAB = 2
386 1360 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
387 :
388 1360 : call cpu_time(time2)
389 1360 : time_proj = time_proj + time2-time1
390 1360 : counter_proj= counter_proj+1
391 :
392 :
393 : !-----------------------------------------------------------------
394 : ! VII) Loop on potential valence operators,
395 : !-----------------------------------------------------------------
396 3704 : do i_op_v = 1, num_op_v
397 :
398 2160 : if (case_op_v == 1) then
399 :
400 : ! frequency is zero, operator to apply is 1/[H-ev]
401 560 : call cpu_time(time1)
402 560 : GWLS_TIMAB = 1529
403 560 : OPTION_TIMAB = 1
404 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
405 :
406 560 : call sqmr(psik_in_alltoall,psik_tmp_alltoall,eig(v),1)
407 :
408 560 : OPTION_TIMAB = 2
409 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
410 :
411 :
412 560 : call cpu_time(time2)
413 560 : time_sqmr = time_sqmr+ time2-time1
414 560 : counter_sqmr= counter_sqmr+1
415 :
416 : ! If we are constructing the $\hat \epsilon(i\omega = 0)$ matrix (and the Lanczos basis at the same time),
417 : ! keep the Sternheimer solutions for use in the projected Sternheimer section (in LA configuration).
418 560 : if(write_solution .and. ((iblk-1)*blocksize < nbandv)) then
419 160 : call wf_block_distribute(psik, psik_tmp_alltoall, 2) ! FFT -> LA
420 384 : do mb=1,blocksize
421 224 : v = (iblk-1)*blocksize+mb
422 384 : if(v <= nbandv) then
423 224 : if(dtset%gwls_recycle == 1) then
424 24800 : Sternheimer_solutions_zero(:,:,index_solution,v) = psik(:,(mb-1)*npw_k+1:mb*npw_k)
425 : end if
426 224 : if(dtset%gwls_recycle == 2) then
427 96 : recy_i = (index_solution-1)*nbandv + v
428 : !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).
429 : !Workaround compatible only with nag : write(recy_unit,rec=recy_i+1).
430 96 : write(recy_unit,rec=recy_i) psik(:,(mb-1)*npw_k+1:mb*npw_k)
431 : end if
432 : end if
433 : end do
434 : end if
435 :
436 1600 : else if (case_op_v == 2) then
437 :
438 : ! frequency purely imaginary, operator to apply is
439 : ! [H-ev]/(lbda^2+[H-ev]^2)
440 0 : call cpu_time(time1)
441 :
442 0 : GWLS_TIMAB = 1533
443 0 : OPTION_TIMAB = 1
444 0 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
445 :
446 0 : psik_wrk_alltoall = psik_in_alltoall
447 0 : call Hpsik(psik_wrk_alltoall,eig(v))
448 :
449 0 : OPTION_TIMAB = 2
450 0 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
451 :
452 :
453 0 : call cpu_time(time2)
454 0 : time_H = time_H + time2-time1
455 0 : counter_H = counter_H+1
456 :
457 0 : call cpu_time(time1)
458 :
459 0 : GWLS_TIMAB = 1530
460 0 : OPTION_TIMAB = 1
461 0 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
462 :
463 0 : call sqmr(psik_wrk_alltoall,psik_tmp_alltoall,eig(v),0,norm_omega,omega_imaginary)
464 :
465 0 : OPTION_TIMAB = 2
466 0 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
467 :
468 :
469 0 : call cpu_time(time2)
470 0 : time_sqmr = time_sqmr+ time2-time1
471 0 : counter_sqmr= counter_sqmr+1
472 :
473 1600 : else if (case_op_v == 3) then
474 :
475 : ! frequency purely real, operator to apply is
476 : ! 1/[H-ev +/- lbda]
477 : ! TREATED WITH SQMR
478 :
479 :
480 1280 : lbda = list_SQMR_frequencies(i_op_v)
481 1280 : call cpu_time(time1)
482 :
483 1280 : GWLS_TIMAB = 1531
484 1280 : OPTION_TIMAB = 1
485 1280 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
486 :
487 1280 : call sqmr(psik_in_alltoall,psik_tmp_alltoall,lbda,1)
488 :
489 1280 : OPTION_TIMAB = 2
490 1280 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
491 :
492 :
493 :
494 1280 : call cpu_time(time2)
495 1280 : time_sqmr = time_sqmr+ time2-time1
496 1280 : counter_sqmr= counter_sqmr+1
497 :
498 320 : else if (case_op_v == 4) then
499 : ! frequency purely real, operator to apply is
500 : ! 1/[H-ev +/- lbda]
501 : ! TREATED WITH QMR
502 :
503 960 : zz(:) = list_QMR_frequencies(:,i_op_v)
504 :
505 320 : call cpu_time(time1)
506 320 : GWLS_TIMAB = 1532
507 320 : OPTION_TIMAB = 1
508 320 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
509 :
510 320 : call qmr(psik_in_alltoall,psik_tmp_alltoall,zz) !,0)
511 :
512 320 : OPTION_TIMAB = 2
513 320 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
514 :
515 :
516 320 : call cpu_time(time2)
517 320 : time_sqmr = time_sqmr+ time2-time1
518 320 : counter_sqmr= counter_sqmr+1
519 :
520 : end if
521 : !-----------------------------------------------------------------
522 : ! VIII) Project on conduction states
523 : !-----------------------------------------------------------------
524 2160 : call cpu_time(time1)
525 2160 : GWLS_TIMAB = 1528
526 2160 : OPTION_TIMAB = 1
527 2160 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
528 :
529 2160 : call pc_k_valence_kernel(psik_tmp_alltoall)
530 :
531 2160 : OPTION_TIMAB = 2
532 2160 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
533 :
534 :
535 2160 : call cpu_time(time2)
536 2160 : time_proj = time_proj + time2-time1
537 2160 : counter_proj= counter_proj+1
538 :
539 : !-----------------------------------------------------------------
540 : ! IX) Conjugate result, and express in denpot format
541 : !-----------------------------------------------------------------
542 :
543 2160 : call cpu_time(time1)
544 :
545 2160 : GWLS_TIMAB = 1526
546 2160 : OPTION_TIMAB = 1
547 2160 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
548 :
549 :
550 2160 : call g_to_r(psir,psik_tmp_alltoall)
551 14815440 : psir(2,:,:,:) = -psir(2,:,:,:)
552 :
553 :
554 2160 : OPTION_TIMAB = 2
555 2160 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
556 :
557 2160 : call cpu_time(time2)
558 2160 : time_fft = time_fft + time2-time1
559 2160 : counter_fft = counter_fft+1
560 :
561 : !-----------------------------------------------------------------
562 : ! X) Multiply by valence state in real space
563 : !-----------------------------------------------------------------
564 2160 : call cpu_time(time1)
565 2160 : GWLS_TIMAB = 1527
566 2160 : OPTION_TIMAB = 1
567 2160 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
568 :
569 2160 : call gr_to_g(psik_alltoall,psir,valence_wavefunctions_FFT(:,:,iblk))
570 :
571 :
572 2160 : OPTION_TIMAB = 2
573 2160 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
574 :
575 :
576 2160 : call cpu_time(time2)
577 2160 : time_fft = time_fft + time2-time1
578 2160 : counter_fft = counter_fft+1
579 :
580 : !-----------------------------------------------------------------
581 : ! XI) Return to LA configuration, and cumulate the sum
582 : !-----------------------------------------------------------------
583 2160 : call wf_block_distribute(psik, psik_alltoall,2) ! FFT -> LA
584 :
585 6480 : do mb = 1, blocksize
586 :
587 2960 : v = (iblk-1)*blocksize + mb
588 :
589 5120 : if (v <= nbandv) then
590 : ! only add contributions from valence
591 543968 : psi_inout(:,:) = psi_inout(:,:) + psik(:,(mb-1)*npw_k+1:mb*npw_k)
592 : end if
593 :
594 : end do
595 :
596 :
597 : end do ! i_op_v
598 :
599 : end do ! iblk
600 :
601 :
602 : !-----------------------------------------------------------------
603 : ! XII) account for prefactor
604 : !-----------------------------------------------------------------
605 86536 : psi_inout = factor_op_v*psi_inout
606 :
607 :
608 184 : GWLS_TIMAB = 1525
609 184 : OPTION_TIMAB = 1
610 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
611 :
612 :
613 184 : ABI_FREE(psik)
614 184 : ABI_FREE(psik_alltoall)
615 184 : ABI_FREE(psik_wrk_alltoall)
616 184 : ABI_FREE(psik_tmp_alltoall)
617 184 : ABI_FREE(psik_in_alltoall)
618 :
619 184 : ABI_FREE(psir)
620 184 : ABI_FREE(psir_ext)
621 :
622 184 : OPTION_TIMAB = 2
623 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
624 :
625 :
626 184 : call cpu_time(total_time2)
627 :
628 184 : total_time = total_time + total_time2 - total_time1
629 :
630 :
631 184 : GWLS_TIMAB = 1524
632 : OPTION_TIMAB = 2
633 184 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
634 :
635 : 10 format(A)
636 : 12 format(A,I5)
637 : 14 format(A,F24.16)
638 :
639 184 : end subroutine Pk
640 : !!***
641 :
642 : !!****f* m_hamiltonian/epsilon_k
643 : !! NAME
644 : !! epsilon_k
645 : !!
646 : !! FUNCTION
647 : !! .
648 : !!
649 : !! INPUTS
650 : !!
651 : !! OUTPUT
652 : !!
653 : !! SOURCE
654 :
655 184 : subroutine epsilon_k(psi_out,psi_in,omega)
656 :
657 : real(dp), intent(out) :: psi_out(2,npw_k)
658 : real(dp), intent(in) :: psi_in(2,npw_k), omega(2)
659 :
660 : ! *************************************************************************
661 :
662 86536 : psi_out = psi_in
663 184 : call sqrt_vc_k(psi_out)
664 184 : call Pk(psi_out,omega)
665 184 : call sqrt_vc_k(psi_out)
666 :
667 86536 : psi_out = psi_in - psi_out
668 184 : end subroutine epsilon_k
669 : !!***
670 :
671 : !!****f* m_hamiltonian/set_dielectric_function_frequency
672 : !! NAME
673 : !! set_dielectric_function_frequency
674 : !!
675 : !! FUNCTION
676 : !! .
677 : !!
678 : !! INPUTS
679 : !!
680 : !! OUTPUT
681 : !!
682 : !! SOURCE
683 :
684 22 : subroutine set_dielectric_function_frequency(omega)
685 : !----------------------------------------------------------------------------------------------------
686 : ! This routine sets the value of the module's frequency.
687 : !----------------------------------------------------------------------------------------------------
688 : real(dp), intent(in) :: omega(2)
689 :
690 : ! *************************************************************************
691 :
692 22 : matrix_function_omega(:) = omega(:)
693 :
694 22 : end subroutine set_dielectric_function_frequency
695 : !!***
696 :
697 : !!****f* m_hamiltonian/matrix_function_epsilon_k
698 : !! NAME
699 : !! matrix_function_epsilon_k
700 : !!
701 : !! FUNCTION
702 : !! .
703 : !!
704 : !! INPUTS
705 : !!
706 : !! OUTPUT
707 : !!
708 : !! SOURCE
709 :
710 184 : subroutine matrix_function_epsilon_k(vector_out,vector_in,Hsize)
711 : !----------------------------------------------------------------------------------------------------
712 : ! This function is a simple wrapper around epsilon_k to be fed to the Lanczos
713 : ! algorithm.
714 : !----------------------------------------------------------------------------------------------------
715 : integer, intent(in) :: Hsize
716 : complex(dp), intent(out) :: vector_out(Hsize)
717 : complex(dp), intent(in) :: vector_in(Hsize)
718 :
719 :
720 : ! local variables
721 368 : real(dp) :: psik (2,npw_k)
722 184 : real(dp) :: psik2(2,npw_k)
723 :
724 : ! *************************************************************************
725 :
726 : ! convert from one format to the other
727 28968 : psik(1,:) = dble (vector_in(:))
728 28968 : psik(2,:) = dimag(vector_in(:))
729 :
730 : ! act on vector
731 :
732 :
733 184 : call epsilon_k(psik2 ,psik, matrix_function_omega)
734 :
735 : ! convert back
736 28968 : vector_out = cmplx_1*psik2(1,:)+cmplx_i*psik2(2,:)
737 :
738 184 : end subroutine matrix_function_epsilon_k
739 : !!***
740 :
741 : end module m_gwls_polarisability
742 : !!***
|
