TABLE OF CONTENTS


ABINIT/m_gwls_model_polarisability [ Modules ]

[ Top ] [ Modules ]

NAME

 m_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 m_gwls_model_polarisability
29 
30 ! local modules
31 use m_gwls_utility
32 use m_gwls_wf
33 use m_gwls_valenceWavefunctions
34 use m_gwls_hamiltonian
35 use m_gwls_lineqsolver
36 
37 ! abinit modules
38 use defs_basis
39 use m_abicore
40 use m_bandfft_kpt
41 use m_errors
42 
43 use m_time,      only : timab
44 
45 implicit none
46 save
47 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

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

m_hamiltonian/epsilon_k_model [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  epsilon_k_model

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      m_gwls_model_polarisability

CHILDREN

      epsilon_k_model

SOURCE

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

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

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

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

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