Line data Source code
1 : !!****m* ABINIT/m_gwls_model_polarisability
2 : !! NAME
3 : !! m_gwls_model_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_model_polarisability
24 :
25 : ! local modules
26 : use m_gwls_utility
27 : use m_gwls_wf
28 : use m_gwls_valenceWavefunctions
29 : use m_gwls_hamiltonian
30 : use m_gwls_lineqsolver
31 :
32 : ! abinit modules
33 : use defs_basis
34 : use m_abicore
35 : use m_bandfft_kpt
36 : use m_errors
37 :
38 : use m_time, only : timab
39 :
40 : implicit none
41 : save
42 : private
43 : !!***
44 :
45 : real(dp), public :: model_polarizability_epsilon_0 ! model parameter
46 : !integer :: model_polarizability_model_type ! how is epsilon_0 used to model?
47 :
48 :
49 : real(dp), allocatable, private :: psir_model(:,:,:,:), psir_ext_model(:,:,:,:)
50 :
51 :
52 : integer, public :: dielectric_model_type = 1
53 :
54 : real(dp),public, allocatable :: model_Y(:) ! model susceptibility, distributed according to FFT configuration
55 : real(dp),public, allocatable :: model_Y_LA(:) ! model susceptibility, distributed according to LA configuration
56 :
57 : real(dp),public, allocatable :: sqrt_density(:,:,:,:) ! average valence wave function...
58 : !!***
59 :
60 : public :: Pk_model
61 :
62 :
63 : public :: epsilon_k_model
64 : public :: setup_Pk_model
65 : public :: cleanup_Pk_model
66 :
67 : public :: matrix_function_epsilon_model_operator
68 : !!***
69 :
70 : contains
71 :
72 : !!****f* m_hamiltonian/epsilon_k_model
73 : !! NAME
74 : !! epsilon_k_model
75 : !!
76 : !! FUNCTION
77 : !! .
78 : !!
79 : !! INPUTS
80 : !!
81 : !! OUTPUT
82 : !!
83 : !! SOURCE
84 :
85 72 : subroutine epsilon_k_model(psi_out,psi_in)
86 :
87 : real(dp), intent(out) :: psi_out(2,npw_k)
88 : real(dp), intent(in) :: psi_in(2,npw_k)
89 :
90 72 : real(dp) :: psik(2,npw_k)
91 :
92 : ! *************************************************************************
93 :
94 37080 : psik = psi_in
95 72 : call sqrt_vc_k(psik)
96 72 : call Pk_model(psi_out ,psik)
97 72 : call sqrt_vc_k(psi_out)
98 37080 : psi_out = psi_in - psi_out
99 :
100 72 : end subroutine epsilon_k_model
101 : !!***
102 :
103 : !!****f* m_hamiltonian/setup_Pk_model
104 : !! NAME
105 : !! setup_Pk_model
106 : !!
107 : !! FUNCTION
108 : !! .
109 : !!
110 : !! INPUTS
111 : !!
112 : !! OUTPUT
113 : !!
114 : !! SOURCE
115 :
116 511 : subroutine setup_Pk_model(omega,epsilon_0)
117 : !---------------------------------------------------------------
118 : !
119 : ! This subroutine prepares a global array in order to
120 : ! act with the model susceptibility, given by
121 : !
122 : ! Pk_model(r,r',i omega) = sum_v phi_v(r) Y(r-r',i omega) phi^*_v(r')
123 : !
124 : !
125 : ! This subroutine computes the Fourier transform of Y(r,i omega),
126 : ! Y(G,omega), for a given omega and epsilon_0.
127 : !
128 : ! It is assumed that omega and epsilon_0 >= 0. Also, the model
129 : ! describes an IMAGINARY frequency i omega.
130 : !---------------------------------------------------------------
131 : real(dp), intent(in) :: epsilon_0, omega
132 :
133 : real(dp) :: theta, R_omega
134 : real(dp) :: x, y
135 :
136 :
137 : integer :: ig
138 :
139 511 : real(dp),allocatable :: G_array(:)
140 :
141 : ! *************************************************************************
142 :
143 511 : if (.not. allocated(psir_model)) then
144 30 : ABI_MALLOC(psir_model, (2,n4,n5,n6))
145 : end if
146 :
147 511 : if (.not. allocated(psir_ext_model)) then
148 30 : ABI_MALLOC(psir_ext_model, (2,n4,n5,n6))
149 : end if
150 :
151 511 : R_omega = 2.0_dp*sqrt(epsilon_0**2+omega**2)
152 511 : if(abs(omega) > tol16 .or. abs(epsilon_0) > tol16) then
153 480 : theta = atan2(omega,epsilon_0)
154 : else
155 : theta = zero
156 : end if
157 :
158 511 : x = sqrt(R_omega)*cos(0.5_dp*theta)
159 511 : y = sqrt(R_omega)*sin(0.5_dp*theta)
160 :
161 :
162 511 : if (.not. allocated(model_Y)) then
163 18 : ABI_MALLOC(model_Y, (npw_g))
164 : end if
165 :
166 511 : if (.not. allocated(model_Y_LA)) then
167 18 : ABI_MALLOC(model_Y_LA, (npw_k))
168 : end if
169 :
170 : !================================================================================
171 : ! Compute model_Y, in FFT configuration
172 : !================================================================================
173 :
174 1533 : ABI_MALLOC(G_array,(npw_g))
175 88148 : G_array(:) = sqrt(2.0_dp*kinpw_gather(:))
176 :
177 88148 : model_Y(:) = zero
178 88148 : do ig = 1, npw_g
179 :
180 :
181 88148 : if (G_array(ig) > tol12) then
182 : ! G != 0.
183 : model_Y(ig) = -4.0_dp/G_array(ig)* &
184 : ((G_array(ig)+y)/((G_array(ig)+y)**2+x**2) &
185 87296 : + (G_array(ig)-y)/((G_array(ig)-y)**2+x**2))
186 :
187 :
188 : else
189 341 : if ( abs(epsilon_0) < tol12 ) then
190 341 : model_Y(ig) = zero
191 : else
192 0 : model_Y(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2)
193 : end if
194 : end if
195 : end do ! ig
196 :
197 511 : ABI_FREE(G_array)
198 :
199 : !================================================================================
200 : ! Compute model_Y_LA, in LA configuration
201 : !================================================================================
202 :
203 1533 : ABI_MALLOC(G_array,(npw_k))
204 66303 : G_array(:) = sqrt(2.0_dp*kinpw(:))
205 :
206 66303 : model_Y_LA(:) = zero
207 66303 : do ig = 1, npw_k
208 :
209 66303 : if (G_array(ig) > tol12) then
210 : ! G != 0.
211 : model_Y_LA(ig) = -4.0_dp/G_array(ig)* &
212 : ((G_array(ig)+y)/((G_array(ig)+y)**2+x**2) &
213 65536 : + (G_array(ig)-y)/((G_array(ig)-y)**2+x**2))
214 :
215 :
216 : else
217 256 : if ( abs(epsilon_0) < tol12 ) then
218 256 : model_Y_LA(ig) = zero
219 : else
220 0 : model_Y_LA(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2)
221 : end if
222 : end if
223 : end do ! ig
224 :
225 511 : ABI_FREE(G_array)
226 :
227 :
228 :
229 511 : if (dielectric_model_type == 2) then
230 :
231 0 : ABI_BUG('dielectric_model_type == 2 not properly implemented. Review code or input!')
232 :
233 : !ABI_MALLOC(sqrt_density,(2,n4,n5,n6))
234 : !sqrt_density(:,:,:,:) = zero
235 : !do v= 1, nbandv
236 : ! sqrt_density(1,:,:,:) = sqrt_density(1,:,:,:) + valence_wfr(1,:,:,:,v)**2+valence_wfr(2,:,:,:,v)**2
237 : !end do
238 : !sqrt_density(1,:,:,:) = sqrt(sqrt_density(1,:,:,:))
239 :
240 : end if
241 :
242 :
243 511 : end subroutine setup_Pk_model
244 : !!***
245 :
246 : !!****f* m_hamiltonian/cleanup_Pk_model
247 : !! NAME
248 : !! cleanup_Pk_model
249 : !!
250 : !! FUNCTION
251 : !! .
252 : !!
253 : !! INPUTS
254 : !!
255 : !! OUTPUT
256 : !!
257 : !! SOURCE
258 :
259 7 : subroutine cleanup_Pk_model()
260 :
261 : ! *************************************************************************
262 :
263 7 : ABI_SFREE(model_Y)
264 7 : ABI_SFREE(model_Y_LA)
265 7 : ABI_SFREE(sqrt_density)
266 7 : ABI_SFREE(psir_model)
267 7 : ABI_SFREE(psir_ext_model)
268 :
269 7 : end subroutine cleanup_Pk_model
270 : !!***
271 :
272 : !!****f* m_hamiltonian/Pk_model
273 : !! NAME
274 : !! Pk_model
275 : !!
276 : !! FUNCTION
277 : !! .
278 : !!
279 : !! INPUTS
280 : !!
281 : !! OUTPUT
282 : !!
283 : !! SOURCE
284 :
285 72 : subroutine Pk_model(psi_out,psi_in)
286 : !------------------------------------------------------------------------------------------------------------------------
287 : ! Returns the action of a frequency-dependent model susceptibility
288 : !------------------------------------------------------------------------------------------------------------------------
289 :
290 : real(dp), intent(out) :: psi_out(2,npw_k)
291 : real(dp), intent(in) :: psi_in(2,npw_k)
292 :
293 : ! *************************************************************************
294 :
295 72 : if ( dielectric_model_type == 1) then
296 72 : call Pk_model_implementation_1(psi_out ,psi_in)
297 : else if ( dielectric_model_type == 2) then
298 : ! call Pk_model_implementation_2(psi_out ,psi_in)
299 : end if
300 :
301 0 : end subroutine Pk_model
302 : !!***
303 :
304 : !!****f* m_hamiltonian/Pk_model_implementation_1
305 : !! NAME
306 : !! Pk_model_implementation_1
307 : !!
308 : !! FUNCTION
309 : !! .
310 : !!
311 : !! INPUTS
312 : !!
313 : !! OUTPUT
314 : !!
315 : !! SOURCE
316 :
317 72 : subroutine Pk_model_implementation_1(psi_out,psi_in)
318 : !------------------------------------------------------------------------------------------------------------------------
319 : ! Returns the action of a frequency-dependent model susceptibility
320 : !------------------------------------------------------------------------------------------------------------------------
321 :
322 : real(dp), intent(out) :: psi_out(2,npw_k)
323 : real(dp), intent(in) :: psi_in(2,npw_k)
324 :
325 : integer :: v
326 :
327 72 : real(dp), allocatable :: psik(:,:), psik_g(:,:)
328 :
329 : integer, save :: icounter = 0
330 :
331 : integer :: mb, iblk
332 :
333 : integer :: mpi_band_rank
334 :
335 : real(dp) :: time1, time2
336 : real(dp) :: total_time1, total_time2
337 :
338 : real(dp), save :: fft_time = zero
339 : real(dp), save :: projection_time = zero
340 : real(dp), save :: Y_time = zero
341 : real(dp), save :: total_time = zero
342 :
343 :
344 : real(dp) :: tsec(2)
345 : integer :: GWLS_TIMAB, OPTION_TIMAB
346 :
347 : ! *************************************************************************
348 :
349 72 : GWLS_TIMAB = 1534
350 72 : OPTION_TIMAB = 1
351 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
352 :
353 :
354 :
355 72 : call cpu_time(total_time1)
356 72 : icounter = icounter + 1
357 :
358 72 : GWLS_TIMAB = 1535
359 : OPTION_TIMAB = 1
360 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
361 :
362 :
363 216 : ABI_MALLOC(psik, (2,npw_kb))
364 216 : ABI_MALLOC(psik_g, (2,npw_g))
365 :
366 72 : OPTION_TIMAB = 2
367 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
368 :
369 :
370 : ! initialize the output to zero
371 37080 : psi_out = zero
372 :
373 : ! MPI information
374 72 : mpi_band_rank = mpi_enreg%me_band
375 :
376 :
377 : !-----------------------------------------------------------------
378 : ! Put a copy of the external state on every row of FFT processors.
379 : !
380 : ! The, copy conjugate of initial wavefunction in local array,
381 : ! and set inout array to zero.
382 : !-----------------------------------------------------------------
383 72 : call cpu_time(time1)
384 72 : GWLS_TIMAB = 1536
385 72 : OPTION_TIMAB = 1
386 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
387 :
388 : ! fill the array psik_ext with copies of the external state
389 176 : do mb = 1, blocksize
390 43352 : psik(:,(mb-1)*npw_k+1:mb*npw_k) = psi_in(:,:)
391 : end do
392 :
393 : ! change configuration of the data, from LA to FFT
394 72 : call wf_block_distribute(psik, psik_g,1) ! LA -> FFT
395 : ! Now every row of FFT processors has a copy of the external state.
396 :
397 :
398 : ! Store external state in real-space format, inside module-defined work array psir_ext_model
399 : ! Each row of FFT processors will have a copy!
400 72 : call g_to_r(psir_ext_model,psik_g)
401 :
402 : ! Conjugate the external wavefunction; result will be conjugated again later,
403 : ! insuring we are in fact acting on psi, and not psi^*. This conjugation
404 : ! is only done for algorithmic convenience.
405 493848 : psir_ext_model(2,:,:,:) = -psir_ext_model(2,:,:,:)
406 :
407 :
408 72 : OPTION_TIMAB = 2
409 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
410 :
411 :
412 72 : call cpu_time(time2)
413 72 : fft_time = fft_time+time2-time1
414 :
415 :
416 : ! Loop on all blocks of eigenstates
417 632 : do iblk = 1, nbdblock
418 :
419 560 : v = (iblk-1)*blocksize + mpi_band_rank + 1 ! CAREFUL! This is a guess. Revisit this if code doesn't work as intended.
420 :
421 560 : call cpu_time(time1)
422 560 : GWLS_TIMAB = 1537
423 560 : OPTION_TIMAB = 1
424 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
425 :
426 :
427 : ! Multiply valence state by external state, yielding |psik_g> = | phi_v x psi_in^* >
428 560 : call gr_to_g(psik_g, psir_ext_model, valence_wavefunctions_FFT(:,:,iblk))
429 :
430 560 : OPTION_TIMAB = 2
431 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
432 :
433 :
434 560 : call cpu_time(time2)
435 560 : fft_time = fft_time+time2-time1
436 :
437 : ! Project out to conduction space
438 560 : call cpu_time(time1)
439 560 : GWLS_TIMAB = 1538
440 560 : OPTION_TIMAB = 1
441 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
442 :
443 560 : call pc_k_valence_kernel(psik_g)
444 :
445 560 : OPTION_TIMAB = 2
446 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
447 :
448 :
449 560 : call cpu_time(time2)
450 560 : projection_time = projection_time+time2-time1
451 :
452 :
453 : ! act with model susceptibility
454 560 : call cpu_time(time1)
455 560 : GWLS_TIMAB = 1539
456 560 : OPTION_TIMAB = 1
457 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
458 :
459 :
460 123920 : psik_g(1,:) = psik_g(1,:)*model_Y(:)
461 123920 : psik_g(2,:) = psik_g(2,:)*model_Y(:)
462 :
463 560 : OPTION_TIMAB = 2
464 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
465 :
466 :
467 560 : call cpu_time(time2)
468 560 : Y_time = Y_time+time2-time1
469 :
470 : ! Project out to conduction space, again!
471 560 : call cpu_time(time1)
472 560 : GWLS_TIMAB = 1538
473 560 : OPTION_TIMAB = 1
474 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
475 :
476 560 : call pc_k_valence_kernel(psik_g)
477 :
478 560 : OPTION_TIMAB = 2
479 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
480 560 : call cpu_time(time2)
481 :
482 560 : projection_time= projection_time+time2-time1
483 :
484 : ! Express result in real space, in module-defined work array psir_model
485 560 : call cpu_time(time1)
486 560 : GWLS_TIMAB = 1536
487 560 : OPTION_TIMAB = 1
488 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
489 :
490 560 : call g_to_r(psir_model,psik_g)
491 : ! conjugate the result, cancelling the initial conjugation described earlier.
492 3841040 : psir_model(2,:,:,:) = -psir_model(2,:,:,:)
493 :
494 560 : OPTION_TIMAB = 2
495 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
496 :
497 560 : call cpu_time(time2)
498 560 : fft_time = fft_time+time2-time1
499 :
500 : ! Multiply by valence state in real space
501 560 : call cpu_time(time1)
502 560 : GWLS_TIMAB = 1537
503 560 : OPTION_TIMAB = 1
504 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
505 :
506 560 : call gr_to_g(psik_g,psir_model, valence_wavefunctions_FFT(:,:,iblk))
507 :
508 560 : OPTION_TIMAB = 2
509 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
510 560 : call cpu_time(time2)
511 560 : fft_time = fft_time+time2-time1
512 :
513 :
514 :
515 : ! Return to linear algebra format, and add condtribution
516 560 : GWLS_TIMAB = 1540
517 560 : OPTION_TIMAB = 1
518 560 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
519 :
520 560 : call wf_block_distribute(psik, psik_g,2) ! FFT -> LA
521 :
522 1280 : do mb = 1, blocksize
523 :
524 720 : v = (iblk-1)*blocksize + mb
525 1280 : if ( v <= nbandv) then
526 148320 : psi_out(:,:) = psi_out(:,:) + psik(:,(mb-1)*npw_k+1:mb*npw_k)
527 : end if
528 : end do
529 :
530 560 : OPTION_TIMAB = 2
531 632 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
532 :
533 :
534 : end do ! iblk
535 :
536 72 : call cpu_time(total_time2)
537 72 : total_time = total_time + total_time2-total_time1
538 :
539 72 : GWLS_TIMAB = 1535
540 72 : OPTION_TIMAB = 1
541 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
542 :
543 72 : ABI_FREE(psik)
544 72 : ABI_FREE(psik_g)
545 :
546 72 : OPTION_TIMAB = 2
547 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
548 :
549 :
550 72 : GWLS_TIMAB = 1534
551 : OPTION_TIMAB = 2
552 72 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
553 :
554 :
555 72 : end subroutine Pk_model_implementation_1
556 : !!***
557 :
558 : !!****f* m_hamiltonian/matrix_function_epsilon_model_operator
559 : !! NAME
560 : !! matrix_function_epsilon_model_operator
561 : !!
562 : !! FUNCTION
563 : !! .
564 : !!
565 : !! INPUTS
566 : !!
567 : !! OUTPUT
568 : !!
569 : !! SOURCE
570 :
571 72 : subroutine matrix_function_epsilon_model_operator(vector_out,vector_in,Hsize)
572 : !----------------------------------------------------------------------------------------------------
573 : ! This function returns the action of the operator epsilon_model on a given vector.
574 : ! It is assumed that the frequency has been set in the module using setup_Pk_model.
575 : !
576 : !
577 : !----------------------------------------------------------------------------------------------------
578 : integer, intent(in) :: Hsize
579 : complex(dp), intent(out) :: vector_out(Hsize)
580 : complex(dp), intent(in) :: vector_in(Hsize)
581 :
582 : ! local variables
583 144 : real(dp) :: psik (2,Hsize)
584 72 : real(dp) :: psik2(2,Hsize)
585 : ! *************************************************************************
586 :
587 : ! convert from one format to the other
588 12408 : psik(1,:) = dble (vector_in(:))
589 12408 : psik(2,:) = dimag(vector_in(:))
590 :
591 72 : call epsilon_k_model(psik2 ,psik)
592 :
593 : ! Act with epsilon_model
594 12408 : vector_out = cmplx_1*psik2(1,:)+cmplx_i*psik2(2,:)
595 :
596 72 : end subroutine matrix_function_epsilon_model_operator
597 : !!***
598 :
599 :
600 : end module m_gwls_model_polarisability
601 : !!***
|
