TABLE OF CONTENTS

ABINIT/gwls_model_polarisability [ Modules ]

[ Top ] [ Modules ]

NAME

 gwls_model_polarisability

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 
28 module gwls_model_polarisability
29 ! local modules
30 use m_gwls_utility
31 use gwls_wf
32 use gwls_valenceWavefunctions
33 use gwls_hamiltonian
34 use gwls_lineqsolver
35 
36 ! abinit modules
37 use defs_basis
38 use m_profiling_abi
39 use m_bandfft_kpt
40 use m_errors
41 
42 implicit none
43 save
44 private

m_hamiltonian/cleanup_Pk_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  cleanup_Pk_model

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy

CHILDREN

      epsilon_k_model

SOURCE

295 subroutine cleanup_Pk_model()
296 
297 
298 !This section has been created automatically by the script Abilint (TD).
299 !Do not modify the following lines by hand.
300 #undef ABI_FUNC
301 #define ABI_FUNC 'cleanup_Pk_model'
302 !End of the abilint section
303 
304 implicit none
305 
306 ! *************************************************************************
307 
308 if (allocated(model_Y)) then
309   ABI_DEALLOCATE(model_Y)
310 end if
311 
312 if (allocated(model_Y_LA)) then
313   ABI_DEALLOCATE(model_Y_LA)
314 end if
315 
316 
317 
318 if (allocated(sqrt_density)) then
319   ABI_DEALLOCATE(sqrt_density)
320 end if
321 
322 if ( allocated(psir_model)) then
323   ABI_DEALLOCATE(psir_model)
324 end if
325 
326 if (allocated(psir_ext_model)) then
327   ABI_DEALLOCATE(psir_ext_model)
328 end if
329 
330 end subroutine cleanup_Pk_model

m_hamiltonian/epsilon_k_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  epsilon_k_model

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_model_polarisability

CHILDREN

      epsilon_k_model

SOURCE

 93 subroutine epsilon_k_model(psi_out,psi_in)
 94 
 95 
 96 !This section has been created automatically by the script Abilint (TD).
 97 !Do not modify the following lines by hand.
 98 #undef ABI_FUNC
 99 #define ABI_FUNC 'epsilon_k_model'
100 !End of the abilint section
101 
102 implicit none
103 
104 real(dp), intent(out) :: psi_out(2,npw_k)
105 real(dp), intent(in)  :: psi_in(2,npw_k)
106 
107 real(dp) :: psik(2,npw_k)
108 
109 ! *************************************************************************
110 
111 psik = psi_in
112 call sqrt_vc_k(psik)
113 call Pk_model(psi_out ,psik)
114 call sqrt_vc_k(psi_out)
115 psi_out = psi_in - psi_out
116 
117 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

PARENTS

      gwls_GenerateEpsilon

CHILDREN

      epsilon_k_model

SOURCE

667 subroutine matrix_function_epsilon_model_operator(vector_out,vector_in,Hsize)
668 !----------------------------------------------------------------------------------------------------
669 ! This function returns the action of the operator epsilon_model on a given vector.
670 ! It is assumed that the frequency has been set in the module using setup_Pk_model.
671 !    
672 !
673 !----------------------------------------------------------------------------------------------------
674 
675 !This section has been created automatically by the script Abilint (TD).
676 !Do not modify the following lines by hand.
677 #undef ABI_FUNC
678 #define ABI_FUNC 'matrix_function_epsilon_model_operator'
679 !End of the abilint section
680 
681 implicit none
682 integer,      intent(in)  :: Hsize
683 complex(dpc), intent(out) :: vector_out(Hsize)
684 complex(dpc), intent(in)  :: vector_in(Hsize)
685 
686 ! local variables
687 real(dp)     :: psik (2,Hsize)
688 real(dp)     :: psik2(2,Hsize)
689 
690 ! *************************************************************************
691 
692 ! convert from one format to the other
693 psik(1,:) = dble (vector_in(:))
694 psik(2,:) = dimag(vector_in(:))
695 
696 call epsilon_k_model(psik2 ,psik)
697 
698 ! Act with  epsilon_model
699 vector_out = cmplx_1*psik2(1,:)+cmplx_i*psik2(2,:)
700 
701 end subroutine matrix_function_epsilon_model_operator

m_hamiltonian/Pk_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  Pk_model

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_model_polarisability

CHILDREN

      epsilon_k_model

SOURCE

352 subroutine Pk_model(psi_out,psi_in)
353 !------------------------------------------------------------------------------------------------------------------------
354 ! Returns the action of a frequency-dependent model susceptibility
355 !------------------------------------------------------------------------------------------------------------------------
356 
357 !This section has been created automatically by the script Abilint (TD).
358 !Do not modify the following lines by hand.
359 #undef ABI_FUNC
360 #define ABI_FUNC 'Pk_model'
361 !End of the abilint section
362 
363 implicit none
364 
365 real(dp), intent(out) :: psi_out(2,npw_k)
366 real(dp), intent(in)  :: psi_in(2,npw_k)
367 
368 ! *************************************************************************
369 
370 if ( dielectric_model_type == 1) then
371   call Pk_model_implementation_1(psi_out ,psi_in)
372 else if ( dielectric_model_type == 2) then
373 !  call Pk_model_implementation_2(psi_out ,psi_in)
374 end if
375 
376 end subroutine Pk_model

m_hamiltonian/Pk_model_implementation_1 [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  Pk_model_implementation_1

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_model_polarisability

CHILDREN

      epsilon_k_model

SOURCE

398 subroutine Pk_model_implementation_1(psi_out,psi_in)
399 !------------------------------------------------------------------------------------------------------------------------
400 ! Returns the action of a frequency-dependent model susceptibility
401 !------------------------------------------------------------------------------------------------------------------------
402 
403 !This section has been created automatically by the script Abilint (TD).
404 !Do not modify the following lines by hand.
405 #undef ABI_FUNC
406 #define ABI_FUNC 'Pk_model_implementation_1'
407  use interfaces_18_timing
408 !End of the abilint section
409 
410 implicit none
411 
412 real(dp), intent(out) :: psi_out(2,npw_k)
413 real(dp), intent(in)  :: psi_in(2,npw_k)
414 
415 integer :: v
416 
417 real(dp), allocatable ::   psik(:,:), psik_g(:,:)
418 
419 integer, save  :: icounter = 0
420 
421 integer  :: mb, iblk
422 
423 integer  :: mpi_band_rank    
424 
425 real(dp) :: time1, time2
426 real(dp) :: total_time1, total_time2
427 
428 real(dp), save :: fft_time        = zero
429 real(dp), save :: projection_time = zero
430 real(dp), save :: Y_time          = zero
431 real(dp), save :: total_time      = zero
432 
433 
434 real(dp) :: tsec(2)
435 integer :: GWLS_TIMAB, OPTION_TIMAB
436 
437 ! *************************************************************************
438 
439 GWLS_TIMAB   = 1534
440 OPTION_TIMAB = 1
441 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
442 
443 
444 
445 call cpu_time(total_time1)
446 icounter = icounter + 1
447 
448 GWLS_TIMAB   = 1535
449 OPTION_TIMAB = 1
450 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
451 
452 
453 ABI_ALLOCATE(psik,           (2,npw_kb))
454 ABI_ALLOCATE(psik_g,         (2,npw_g))
455 
456 OPTION_TIMAB = 2
457 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
458 
459 
460 ! initialize the output to zero
461 psi_out = zero
462 
463 ! MPI information
464 mpi_band_rank = mpi_enreg%me_band
465 
466 
467 !-----------------------------------------------------------------
468 ! Put a copy of the external state on every row of FFT processors.
469 !
470 ! The, copy conjugate of initial wavefunction in local array,
471 ! and set inout array to zero.
472 !-----------------------------------------------------------------
473 call cpu_time(time1)
474 GWLS_TIMAB   = 1536
475 OPTION_TIMAB = 1
476 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
477 
478 ! fill the array psik_ext with copies of the external state
479 do mb = 1, blocksize
480 psik(:,(mb-1)*npw_k+1:mb*npw_k)   = psi_in(:,:)
481 end do
482 
483 ! change configuration of the data, from LA to FFT
484 call wf_block_distribute(psik,  psik_g,1) ! LA -> FFT
485 ! Now every row of FFT processors has a copy of the external state.
486 
487 
488 ! Store external state in real-space format, inside module-defined work array psir_ext_model
489 ! Each row of FFT processors will have a copy!
490 call g_to_r(psir_ext_model,psik_g)
491 
492 ! Conjugate the external wavefunction; result will be conjugated again later,
493 ! insuring we are in fact acting on psi, and not psi^*. This conjugation
494 ! is only done for algorithmic convenience.
495 psir_ext_model(2,:,:,:) = -psir_ext_model(2,:,:,:)
496 
497 
498 OPTION_TIMAB = 2
499 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
500 
501 
502 call cpu_time(time2)
503 fft_time = fft_time+time2-time1
504 
505 
506 ! Loop on all blocks of eigenstates
507 do iblk = 1, nbdblock
508 
509 v = (iblk-1)*blocksize + mpi_band_rank + 1 ! CAREFUL! This is a guess. Revisit this if code doesn't work as intended. 
510 
511 call cpu_time(time1)
512 GWLS_TIMAB   = 1537
513 OPTION_TIMAB = 1
514 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
515 
516 
517 ! Multiply valence state by external state, yielding  |psik_g> = | phi_v x psi_in^* > 
518 call gr_to_g(psik_g, psir_ext_model, valence_wavefunctions_FFT(:,:,iblk))
519 
520 OPTION_TIMAB = 2
521 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
522 
523 
524 call cpu_time(time2)
525 fft_time = fft_time+time2-time1
526 
527 ! Project out to conduction space
528 call cpu_time(time1)
529 GWLS_TIMAB   = 1538
530 OPTION_TIMAB = 1
531 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
532 
533 call pc_k_valence_kernel(psik_g)
534 
535 OPTION_TIMAB = 2
536 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
537 
538 
539 call cpu_time(time2)
540 projection_time = projection_time+time2-time1
541 
542 
543 ! act with model susceptibility
544 call cpu_time(time1)
545 GWLS_TIMAB   = 1539
546 OPTION_TIMAB = 1
547 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
548 
549 
550 psik_g(1,:)  = psik_g(1,:)*model_Y(:)
551 psik_g(2,:)  = psik_g(2,:)*model_Y(:)
552 
553 OPTION_TIMAB = 2
554 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
555 
556 
557 call cpu_time(time2)
558 Y_time = Y_time+time2-time1
559 
560 ! Project out to conduction space, again!
561 call cpu_time(time1)
562 GWLS_TIMAB   = 1538
563 OPTION_TIMAB = 1
564 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
565 
566 call pc_k_valence_kernel(psik_g)
567 
568 OPTION_TIMAB = 2
569 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
570 call cpu_time(time2)
571 
572 projection_time= projection_time+time2-time1
573 
574 ! Express result in real space, in module-defined work array psir_model
575 call cpu_time(time1)
576 GWLS_TIMAB   = 1536
577 OPTION_TIMAB = 1
578 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
579 
580 call g_to_r(psir_model,psik_g)
581 ! conjugate the result, cancelling the initial conjugation described earlier.
582 psir_model(2,:,:,:) = -psir_model(2,:,:,:)
583 
584 OPTION_TIMAB = 2
585 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
586 
587 call cpu_time(time2)
588 fft_time = fft_time+time2-time1
589 
590 !  Multiply by valence state in real space 
591 call cpu_time(time1)
592 GWLS_TIMAB   = 1537
593 OPTION_TIMAB = 1
594 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
595 
596 call gr_to_g(psik_g,psir_model, valence_wavefunctions_FFT(:,:,iblk))
597 
598 OPTION_TIMAB = 2
599 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
600 call cpu_time(time2)
601 fft_time = fft_time+time2-time1
602 
603 
604 
605 ! Return to linear algebra format, and add condtribution
606 GWLS_TIMAB   = 1540
607 OPTION_TIMAB = 1
608 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
609 
610 call wf_block_distribute(psik,  psik_g,2) ! FFT -> LA
611 
612 do mb = 1, blocksize
613 
614 v = (iblk-1)*blocksize + mb
615 if ( v <= nbandv) then
616   psi_out(:,:) = psi_out(:,:) + psik(:,(mb-1)*npw_k+1:mb*npw_k)
617 end if
618 end do
619 
620 OPTION_TIMAB = 2
621 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
622 
623 
624 end do ! iblk
625 
626 call cpu_time(total_time2)
627 total_time = total_time + total_time2-total_time1
628 
629 GWLS_TIMAB   = 1535
630 OPTION_TIMAB = 1
631 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
632 
633 ABI_DEALLOCATE(psik)
634 ABI_DEALLOCATE(psik_g)
635 
636 OPTION_TIMAB = 2
637 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
638 
639 
640 GWLS_TIMAB   = 1534
641 OPTION_TIMAB = 2
642 call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
643 
644 
645 end subroutine Pk_model_implementation_1

m_hamiltonian/setup_Pk_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  setup_Pk_model

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputeCorrelationEnergy,gwls_DielectricArray,gwls_GenerateEpsilon

CHILDREN

      epsilon_k_model

SOURCE

139 subroutine setup_Pk_model(omega,epsilon_0)
140 !---------------------------------------------------------------
141 !
142 ! This subroutine prepares a global array in order to 
143 ! act with the model susceptibility, given by
144 !
145 ! Pk_model(r,r',i omega) = sum_v phi_v(r) Y(r-r',i omega) phi^*_v(r')
146 !
147 !
148 ! This subroutine computes the Fourier transform of Y(r,i omega),
149 ! Y(G,omega), for a given omega and epsilon_0.
150 !
151 ! It is assumed that omega and epsilon_0 >= 0. Also, the model
152 ! describes an IMAGINARY frequency i omega.
153 !---------------------------------------------------------------
154 
155 !This section has been created automatically by the script Abilint (TD).
156 !Do not modify the following lines by hand.
157 #undef ABI_FUNC
158 #define ABI_FUNC 'setup_Pk_model'
159 !End of the abilint section
160 
161 real(dp), intent(in) :: epsilon_0, omega
162 
163 real(dp) ::  theta, R_omega
164 real(dp) ::  x, y
165 
166 
167 integer  :: ig
168 
169 real(dp),allocatable ::  G_array(:)
170 
171 ! *************************************************************************
172 
173 if (.not. allocated(psir_model)) then
174   ABI_ALLOCATE(psir_model, (2,n4,n5,n6))
175 end if
176 
177 if (.not. allocated(psir_ext_model)) then
178   ABI_ALLOCATE(psir_ext_model, (2,n4,n5,n6))
179 end if
180 
181 R_omega = 2.0_dp*sqrt(epsilon_0**2+omega**2)
182 if(abs(omega) > tol16 .or. abs(epsilon_0) > tol16) then 
183   theta   = atan2(omega,epsilon_0)
184 else
185   theta   = zero
186 end if
187 
188 x = sqrt(R_omega)*cos(0.5_dp*theta)
189 y = sqrt(R_omega)*sin(0.5_dp*theta)
190 
191 
192 if (.not. allocated(model_Y)) then
193   ABI_ALLOCATE(model_Y, (npw_g))
194 end if
195 
196 if (.not. allocated(model_Y_LA)) then
197   ABI_ALLOCATE(model_Y_LA, (npw_k))
198 end if
199 
200 !================================================================================
201 ! Compute model_Y, in FFT configuration
202 !================================================================================
203 
204 ABI_ALLOCATE(G_array,(npw_g))
205 G_array(:) = sqrt(2.0_dp*kinpw_gather(:))
206 
207 model_Y(:) = zero
208 do ig = 1, npw_g
209 
210 
211 if (G_array(ig) > tol12) then
212   ! G != 0.
213   model_Y(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(ig) = zero
221   else 
222     model_Y(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2)
223   end if
224 end if
225 end do ! ig
226 
227 ABI_DEALLOCATE(G_array)
228 
229 !================================================================================
230 ! Compute model_Y_LA, in LA configuration
231 !================================================================================
232 
233 ABI_ALLOCATE(G_array,(npw_k))
234 G_array(:) = sqrt(2.0_dp*kinpw(:))
235 
236 model_Y_LA(:) = zero
237 do ig = 1, npw_k
238 
239 if (G_array(ig) > tol12) then
240   ! G != 0.
241   model_Y_LA(ig) = -4.0_dp/G_array(ig)*                          &
242   ((G_array(ig)+y)/((G_array(ig)+y)**2+x**2)     &
243   + (G_array(ig)-y)/((G_array(ig)-y)**2+x**2))
244 
245 
246 else
247   if ( abs(epsilon_0) < tol12 ) then
248     model_Y_LA(ig) = zero
249   else 
250     model_Y_LA(ig) = -4.0_dp*epsilon_0/(epsilon_0**2+omega**2)
251   end if
252 end if
253 end do ! ig
254 
255 ABI_DEALLOCATE(G_array)
256 
257 
258 
259 if (dielectric_model_type == 2) then
260 
261   MSG_BUG('dielectric_model_type == 2 not properly implemented. Review code or input!')
262 
263   !ABI_ALLOCATE(sqrt_density,(2,n4,n5,n6))
264   !sqrt_density(:,:,:,:) = zero
265   !do v= 1, nbandv
266   !        sqrt_density(1,:,:,:) = sqrt_density(1,:,:,:) + valence_wfr(1,:,:,:,v)**2+valence_wfr(2,:,:,:,v)**2
267   !end do 
268   !sqrt_density(1,:,:,:) = sqrt(sqrt_density(1,:,:,:))
269 
270 end if 
271 
272 
273 end subroutine setup_Pk_model