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