LCOV - code coverage report
Current view: top level - src/70_gw - m_gwls_model_polarisability.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 97.7 % 173 169
Test Date: 2026-11-10 15:53:16 Functions: 100.0 % 6 6
Branches: 68.1 % 226 154

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

Generated by: LCOV version 2.3-1