TABLE OF CONTENTS


ABINIT/m_gwls_model_polarisability [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_model_polarisability

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2022 ABINIT group (JLJ, BR, MC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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

m_hamiltonian/cleanup_Pk_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  cleanup_Pk_model

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

261 subroutine cleanup_Pk_model()
262 
263 implicit none
264 
265 ! *************************************************************************
266 
267 if (allocated(model_Y)) then
268   ABI_FREE(model_Y)
269 end if
270 
271 if (allocated(model_Y_LA)) then
272   ABI_FREE(model_Y_LA)
273 end if
274 
275 
276 
277 if (allocated(sqrt_density)) then
278   ABI_FREE(sqrt_density)
279 end if
280 
281 if ( allocated(psir_model)) then
282   ABI_FREE(psir_model)
283 end if
284 
285 if (allocated(psir_ext_model)) then
286   ABI_FREE(psir_ext_model)
287 end if
288 
289 end subroutine cleanup_Pk_model

m_hamiltonian/epsilon_k_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  epsilon_k_model

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

 85 subroutine epsilon_k_model(psi_out,psi_in)
 86 
 87 implicit none
 88 
 89 real(dp), intent(out) :: psi_out(2,npw_k)
 90 real(dp), intent(in)  :: psi_in(2,npw_k)
 91 
 92 real(dp) :: psik(2,npw_k)
 93 
 94 ! *************************************************************************
 95 
 96 psik = psi_in
 97 call sqrt_vc_k(psik)
 98 call Pk_model(psi_out ,psik)
 99 call sqrt_vc_k(psi_out)
100 psi_out = psi_in - psi_out
101 
102 end subroutine epsilon_k_model

m_hamiltonian/matrix_function_epsilon_model_operator [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  matrix_function_epsilon_model_operator

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

593 subroutine matrix_function_epsilon_model_operator(vector_out,vector_in,Hsize)
594 !----------------------------------------------------------------------------------------------------
595 ! This function returns the action of the operator epsilon_model on a given vector.
596 ! It is assumed that the frequency has been set in the module using setup_Pk_model.
597 !    
598 !
599 !----------------------------------------------------------------------------------------------------
600 implicit none
601 integer,      intent(in)  :: Hsize
602 complex(dpc), intent(out) :: vector_out(Hsize)
603 complex(dpc), intent(in)  :: vector_in(Hsize)
604 
605 ! local variables
606 real(dp)     :: psik (2,Hsize)
607 real(dp)     :: psik2(2,Hsize)
608 
609 ! *************************************************************************
610 
611 ! convert from one format to the other
612 psik(1,:) = dble (vector_in(:))
613 psik(2,:) = dimag(vector_in(:))
614 
615 call epsilon_k_model(psik2 ,psik)
616 
617 ! Act with  epsilon_model
618 vector_out = cmplx_1*psik2(1,:)+cmplx_i*psik2(2,:)
619 
620 end subroutine matrix_function_epsilon_model_operator

m_hamiltonian/Pk_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  Pk_model

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

305 subroutine Pk_model(psi_out,psi_in)
306 !------------------------------------------------------------------------------------------------------------------------
307 ! Returns the action of a frequency-dependent model susceptibility
308 !------------------------------------------------------------------------------------------------------------------------
309 implicit none
310 
311 real(dp), intent(out) :: psi_out(2,npw_k)
312 real(dp), intent(in)  :: psi_in(2,npw_k)
313 
314 ! *************************************************************************
315 
316 if ( dielectric_model_type == 1) then
317   call Pk_model_implementation_1(psi_out ,psi_in)
318 else if ( dielectric_model_type == 2) then
319 !  call Pk_model_implementation_2(psi_out ,psi_in)
320 end if
321 
322 end subroutine Pk_model

m_hamiltonian/Pk_model_implementation_1 [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  Pk_model_implementation_1

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

338 subroutine Pk_model_implementation_1(psi_out,psi_in)
339 !------------------------------------------------------------------------------------------------------------------------
340 ! Returns the action of a frequency-dependent model susceptibility
341 !------------------------------------------------------------------------------------------------------------------------
342 implicit none
343 
344 real(dp), intent(out) :: psi_out(2,npw_k)
345 real(dp), intent(in)  :: psi_in(2,npw_k)
346 
347 integer :: v
348 
349 real(dp), allocatable ::   psik(:,:), psik_g(:,:)
350 
351 integer, save  :: icounter = 0
352 
353 integer  :: mb, iblk
354 
355 integer  :: mpi_band_rank    
356 
357 real(dp) :: time1, time2
358 real(dp) :: total_time1, total_time2
359 
360 real(dp), save :: fft_time        = zero
361 real(dp), save :: projection_time = zero
362 real(dp), save :: Y_time          = zero
363 real(dp), save :: total_time      = zero
364 
365 
366 real(dp) :: tsec(2)
367 integer :: GWLS_TIMAB, OPTION_TIMAB
368 
369 ! *************************************************************************
370 
371 GWLS_TIMAB   = 1534
372 OPTION_TIMAB = 1
373 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
374 
375 
376 
377 call cpu_time(total_time1)
378 icounter = icounter + 1
379 
380 GWLS_TIMAB   = 1535
381 OPTION_TIMAB = 1
382 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
383 
384 
385 ABI_MALLOC(psik,           (2,npw_kb))
386 ABI_MALLOC(psik_g,         (2,npw_g))
387 
388 OPTION_TIMAB = 2
389 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
390 
391 
392 ! initialize the output to zero
393 psi_out = zero
394 
395 ! MPI information
396 mpi_band_rank = mpi_enreg%me_band
397 
398 
399 !-----------------------------------------------------------------
400 ! Put a copy of the external state on every row of FFT processors.
401 !
402 ! The, copy conjugate of initial wavefunction in local array,
403 ! and set inout array to zero.
404 !-----------------------------------------------------------------
405 call cpu_time(time1)
406 GWLS_TIMAB   = 1536
407 OPTION_TIMAB = 1
408 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
409 
410 ! fill the array psik_ext with copies of the external state
411 do mb = 1, blocksize
412 psik(:,(mb-1)*npw_k+1:mb*npw_k)   = psi_in(:,:)
413 end do
414 
415 ! change configuration of the data, from LA to FFT
416 call wf_block_distribute(psik,  psik_g,1) ! LA -> FFT
417 ! Now every row of FFT processors has a copy of the external state.
418 
419 
420 ! Store external state in real-space format, inside module-defined work array psir_ext_model
421 ! Each row of FFT processors will have a copy!
422 call g_to_r(psir_ext_model,psik_g)
423 
424 ! Conjugate the external wavefunction; result will be conjugated again later,
425 ! insuring we are in fact acting on psi, and not psi^*. This conjugation
426 ! is only done for algorithmic convenience.
427 psir_ext_model(2,:,:,:) = -psir_ext_model(2,:,:,:)
428 
429 
430 OPTION_TIMAB = 2
431 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
432 
433 
434 call cpu_time(time2)
435 fft_time = fft_time+time2-time1
436 
437 
438 ! Loop on all blocks of eigenstates
439 do iblk = 1, nbdblock
440 
441 v = (iblk-1)*blocksize + mpi_band_rank + 1 ! CAREFUL! This is a guess. Revisit this if code doesn't work as intended. 
442 
443 call cpu_time(time1)
444 GWLS_TIMAB   = 1537
445 OPTION_TIMAB = 1
446 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
447 
448 
449 ! Multiply valence state by external state, yielding  |psik_g> = | phi_v x psi_in^* > 
450 call gr_to_g(psik_g, psir_ext_model, valence_wavefunctions_FFT(:,:,iblk))
451 
452 OPTION_TIMAB = 2
453 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
454 
455 
456 call cpu_time(time2)
457 fft_time = fft_time+time2-time1
458 
459 ! Project out to conduction space
460 call cpu_time(time1)
461 GWLS_TIMAB   = 1538
462 OPTION_TIMAB = 1
463 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
464 
465 call pc_k_valence_kernel(psik_g)
466 
467 OPTION_TIMAB = 2
468 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
469 
470 
471 call cpu_time(time2)
472 projection_time = projection_time+time2-time1
473 
474 
475 ! act with model susceptibility
476 call cpu_time(time1)
477 GWLS_TIMAB   = 1539
478 OPTION_TIMAB = 1
479 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
480 
481 
482 psik_g(1,:)  = psik_g(1,:)*model_Y(:)
483 psik_g(2,:)  = psik_g(2,:)*model_Y(:)
484 
485 OPTION_TIMAB = 2
486 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
487 
488 
489 call cpu_time(time2)
490 Y_time = Y_time+time2-time1
491 
492 ! Project out to conduction space, again!
493 call cpu_time(time1)
494 GWLS_TIMAB   = 1538
495 OPTION_TIMAB = 1
496 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
497 
498 call pc_k_valence_kernel(psik_g)
499 
500 OPTION_TIMAB = 2
501 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
502 call cpu_time(time2)
503 
504 projection_time= projection_time+time2-time1
505 
506 ! Express result in real space, in module-defined work array psir_model
507 call cpu_time(time1)
508 GWLS_TIMAB   = 1536
509 OPTION_TIMAB = 1
510 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
511 
512 call g_to_r(psir_model,psik_g)
513 ! conjugate the result, cancelling the initial conjugation described earlier.
514 psir_model(2,:,:,:) = -psir_model(2,:,:,:)
515 
516 OPTION_TIMAB = 2
517 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
518 
519 call cpu_time(time2)
520 fft_time = fft_time+time2-time1
521 
522 !  Multiply by valence state in real space 
523 call cpu_time(time1)
524 GWLS_TIMAB   = 1537
525 OPTION_TIMAB = 1
526 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
527 
528 call gr_to_g(psik_g,psir_model, valence_wavefunctions_FFT(:,:,iblk))
529 
530 OPTION_TIMAB = 2
531 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
532 call cpu_time(time2)
533 fft_time = fft_time+time2-time1
534 
535 
536 
537 ! Return to linear algebra format, and add condtribution
538 GWLS_TIMAB   = 1540
539 OPTION_TIMAB = 1
540 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
541 
542 call wf_block_distribute(psik,  psik_g,2) ! FFT -> LA
543 
544 do mb = 1, blocksize
545 
546 v = (iblk-1)*blocksize + mb
547 if ( v <= nbandv) then
548   psi_out(:,:) = psi_out(:,:) + psik(:,(mb-1)*npw_k+1:mb*npw_k)
549 end if
550 end do
551 
552 OPTION_TIMAB = 2
553 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
554 
555 
556 end do ! iblk
557 
558 call cpu_time(total_time2)
559 total_time = total_time + total_time2-total_time1
560 
561 GWLS_TIMAB   = 1535
562 OPTION_TIMAB = 1
563 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
564 
565 ABI_FREE(psik)
566 ABI_FREE(psik_g)
567 
568 OPTION_TIMAB = 2
569 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
570 
571 
572 GWLS_TIMAB   = 1534
573 OPTION_TIMAB = 2
574 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
575 
576 
577 end subroutine Pk_model_implementation_1 

m_hamiltonian/setup_Pk_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  setup_Pk_model

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

118 subroutine setup_Pk_model(omega,epsilon_0)
119 !---------------------------------------------------------------
120 !
121 ! This subroutine prepares a global array in order to 
122 ! act with the model susceptibility, given by
123 !
124 ! Pk_model(r,r',i omega) = sum_v phi_v(r) Y(r-r',i omega) phi^*_v(r')
125 !
126 !
127 ! This subroutine computes the Fourier transform of Y(r,i omega),
128 ! Y(G,omega), for a given omega and epsilon_0.
129 !
130 ! It is assumed that omega and epsilon_0 >= 0. Also, the model
131 ! describes an IMAGINARY frequency i omega.
132 !---------------------------------------------------------------
133 real(dp), intent(in) :: epsilon_0, omega
134 
135 real(dp) ::  theta, R_omega
136 real(dp) ::  x, y
137 
138 
139 integer  :: ig
140 
141 real(dp),allocatable ::  G_array(:)
142 
143 ! *************************************************************************
144 
145 if (.not. allocated(psir_model)) then
146   ABI_MALLOC(psir_model, (2,n4,n5,n6))
147 end if
148 
149 if (.not. allocated(psir_ext_model)) then
150   ABI_MALLOC(psir_ext_model, (2,n4,n5,n6))
151 end if
152 
153 R_omega = 2.0_dp*sqrt(epsilon_0**2+omega**2)
154 if(abs(omega) > tol16 .or. abs(epsilon_0) > tol16) then 
155   theta   = atan2(omega,epsilon_0)
156 else
157   theta   = zero
158 end if
159 
160 x = sqrt(R_omega)*cos(0.5_dp*theta)
161 y = sqrt(R_omega)*sin(0.5_dp*theta)
162 
163 
164 if (.not. allocated(model_Y)) then
165   ABI_MALLOC(model_Y, (npw_g))
166 end if
167 
168 if (.not. allocated(model_Y_LA)) then
169   ABI_MALLOC(model_Y_LA, (npw_k))
170 end if
171 
172 !================================================================================
173 ! Compute model_Y, in FFT configuration
174 !================================================================================
175 
176 ABI_MALLOC(G_array,(npw_g))
177 G_array(:) = sqrt(2.0_dp*kinpw_gather(:))
178 
179 model_Y(:) = zero
180 do ig = 1, npw_g
181 
182 
183 if (G_array(ig) > tol12) then
184   ! G != 0.
185   model_Y(ig) = -4.0_dp/G_array(ig)*                             &
186   ((G_array(ig)+y)/((G_array(ig)+y)**2+x**2)     &
187   + (G_array(ig)-y)/((G_array(ig)-y)**2+x**2))
188 
189 
190 else
191   if ( abs(epsilon_0) < tol12 ) then
192     model_Y(ig) = zero
193   else 
194     model_Y(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2)
195   end if
196 end if
197 end do ! ig
198 
199 ABI_FREE(G_array)
200 
201 !================================================================================
202 ! Compute model_Y_LA, in LA configuration
203 !================================================================================
204 
205 ABI_MALLOC(G_array,(npw_k))
206 G_array(:) = sqrt(2.0_dp*kinpw(:))
207 
208 model_Y_LA(:) = zero
209 do ig = 1, npw_k
210 
211 if (G_array(ig) > tol12) then
212   ! G != 0.
213   model_Y_LA(ig) = -4.0_dp/G_array(ig)*                          &
214   ((G_array(ig)+y)/((G_array(ig)+y)**2+x**2)     &
215   + (G_array(ig)-y)/((G_array(ig)-y)**2+x**2))
216 
217 
218 else
219   if ( abs(epsilon_0) < tol12 ) then
220     model_Y_LA(ig) = zero
221   else 
222     model_Y_LA(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2)
223   end if
224 end if
225 end do ! ig
226 
227 ABI_FREE(G_array)
228 
229 
230 
231 if (dielectric_model_type == 2) then
232 
233   ABI_BUG('dielectric_model_type == 2 not properly implemented. Review code or input!')
234 
235   !ABI_MALLOC(sqrt_density,(2,n4,n5,n6))
236   !sqrt_density(:,:,:,:) = zero
237   !do v= 1, nbandv
238   !        sqrt_density(1,:,:,:) = sqrt_density(1,:,:,:) + valence_wfr(1,:,:,:,v)**2+valence_wfr(2,:,:,:,v)**2
239   !end do 
240   !sqrt_density(1,:,:,:) = sqrt(sqrt_density(1,:,:,:))
241 
242 end if 
243 
244 
245 end subroutine setup_Pk_model