LCOV - code coverage report
Current view: top level - vgen - vvacxy.f90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 42.3 % 123 52
Test Date: 2025-06-14 04:34:23 Functions: 50.0 % 2 1

            Line data    Source code
       1              : module m_vvacxy
       2              :   !     **********************************************************
       3              :   !     g/=0 coefficient of vacuum coulomb potential           *
       4              :   !     due to warped vacuum charged density                     *
       5              :   !                                 c.l.fu, r.podloucky          *
       6              :   !     **********************************************************
       7              :   !     modified for thick films to avoid underflows gb`06
       8              :   !---------------------------------------------------------------
       9              : 
      10              :    use m_judft
      11              : 
      12              : contains
      13           37 :    subroutine vvacxy( stars, vacuum, cell, sym, input, field, rhtxy, vxy, alphm, l_dfptvgen )
      14              :       use m_types
      15              :       use m_constants
      16              :       use m_intgr, only: intgz1
      17              :       use m_qsf
      18              : 
      19              :       implicit none
      20              : 
      21              :       type(t_input),  intent(in)     :: input
      22              :       type(t_field),  intent(in)     :: field
      23              :       type(t_vacuum), intent(in)     :: vacuum
      24              :       type(t_sym),    intent(in)     :: sym
      25              :       type(t_stars),  intent(in)     :: stars
      26              :       type(t_cell),   intent(in)     :: cell
      27              : 
      28              :       complex,        intent(in)     :: rhtxy(vacuum%nmzxyd,stars%ng2,2)
      29              :       complex,        intent(inout)  :: vxy(vacuum%nmzxyd,stars%ng2,2)
      30              :       complex,        intent(out)    :: alphm(stars%ng2,2)
      31              :       logical,        intent(in)     :: l_dfptvgen
      32              : 
      33              :       complex                        :: alph0, alphaz, betaz, test, phas
      34              :       real                           :: g, vcons, z, e_m, e_p,arg
      35              :       integer                        :: imz, ip, irec2, ivac, ncsh,irec2r, start_star
      36           37 :       real                           :: fra(vacuum%nmzxyd), frb(vacuum%nmzxyd), fia(vacuum%nmzxyd), fib(vacuum%nmzxyd)
      37           37 :       real                           :: alpha(vacuum%nmzxyd,2,2,stars%ng2), beta(vacuum%nmzxyd,2,2,stars%ng2)
      38              : 
      39           37 :       if ( allocated( field%efield%rhoef ) .or. field%efield%dirichlet ) then
      40            0 :          ncsh = field%efield%zsigma / vacuum%delz + 1.01
      41              :          ! if nmzxy < ncsh, the inhomogenous field cannot be represented.
      42              :          ! if nmzxy > ncsh, the boundary condition is wrong - and the
      43              :          ! potential is very wavy - try it yourself, if you don't believe.
      44            0 :          if ( vacuum%nmzxyd < ncsh ) then
      45            0 :             write (oUnit,*) 'error vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '< ncsh(zsigma) = ', ncsh
      46            0 :             call judft_error( "error: vacuum%nmzxyd < ncsh", calledby="vvacxy" )
      47            0 :          else if ( vacuum%nmzxyd > ncsh ) then
      48            0 :             write (oUnit,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
      49            0 :             write (*,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
      50            0 :             call judft_warn( "nmzxyd > ncsh", calledby="vvacxy" )
      51              :          end if
      52              :       end if
      53              : 
      54           37 :       start_star = 2
      55              :       ! For q/=0 in DFPT, there is no G+q=0, so all stars are treated in the G/=0 way.
      56           37 :       if (l_dfptvgen) then
      57            0 :          if (norm2(stars%center)>1e-8) start_star = 1
      58              :       end if
      59              :       ! 2-dim star loop g/=0
      60           37 :       ip = vacuum%nmzxy + 1
      61              :       !    irec2_loop: do irec2 = 2, stars%ng2
      62              :       !      g = stars%sk2(irec2)
      63              :       !      vcons = tpi_const / g
      64           37 :       if ( field%efield%dirichlet ) then ! Dirichlet
      65            0 :          do irec2 = start_star, stars%ng2
      66            0 :             g = stars%sk2(irec2)
      67            0 :             vcons = tpi_const / g
      68            0 :             if ( allocated( field%efield%rhoef ) ) then
      69            0 :                vxy(ncsh:vacuum%nmzxy,irec2,1) = field%efield%rhoef(irec2-1, 1)
      70            0 :                if ( vacuum%nvac == 2 ) then
      71            0 :                   vxy(ncsh:vacuum%nmzxy,irec2,2) = field%efield%rhoef(irec2-1, 2)
      72              :                end if
      73              :             else
      74            0 :                vxy(ncsh:vacuum%nmzxy,irec2,1:vacuum%nvac) = 0.0
      75              :             end if
      76            0 :             vcons = 2.0 * vcons / sinh( g * 2.0 * ( field%efield%zsigma + cell%z1 ) )
      77            0 :             do ivac = 1, vacuum%nvac
      78            0 :                z = cell%z1
      79            0 :                do imz = 1, ncsh-1
      80              :                   ! as "z" > 0 in this subroutine, the integrand is the same
      81              :                   ! for both ivac -- but the integral bounds are reversed
      82            0 :                   e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
      83            0 :                   e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
      84            0 :                   fra(ncsh-imz) = real(  rhtxy(imz,irec2,ivac) ) * e_m
      85            0 :                   fia(ncsh-imz) = aimag( rhtxy(imz,irec2,ivac) ) * e_m
      86            0 :                   frb(imz)      = real(  rhtxy(imz,irec2,ivac) ) * e_p
      87            0 :                   fib(imz)      = aimag( rhtxy(imz,irec2,ivac) ) * e_p
      88            0 :                   z = z + vacuum%delz
      89              :                end do 
      90            0 :                call intgz1( fra, vacuum%delz, ncsh - 1, alpha(1,ivac,1,irec2), .false. )
      91            0 :                call intgz1( fia, vacuum%delz, ncsh - 1, alpha(1,ivac,2,irec2), .false. )
      92            0 :                call qsf( vacuum%delz, frb, beta(1,ivac,1,irec2), ncsh - 1, 1 )
      93            0 :                call qsf( vacuum%delz, fib, beta(1,ivac,2,irec2), ncsh - 1, 1 )
      94              :             end do 
      95              : 
      96            0 :             if ( ivac == 2 ) then
      97              :                ! honour reversed integral bounds
      98            0 :                alpha(:,ivac,:,irec2) = - alpha(:,ivac,:,irec2)
      99            0 :                beta(:,ivac,:,irec2)  = - beta(:,ivac,:,irec2)
     100              :             end if
     101              :          end do
     102            0 :          do irec2 = start_star, stars%ng2
     103            0 :             g = stars%sk2(irec2)
     104            0 :             vcons = tpi_const / g
     105            0 :             alphm(irec2,1) = cmplx( alpha(vacuum%nmzxy,1,1,irec2), alpha(vacuum%nmzxy,1,2,irec2) )
     106            0 :             if ( vacuum%nvac == 1 ) then
     107            0 :                call stars%map_2nd_vac(vacuum,irec2,irec2r,phas)
     108            0 :                alphm(irec2r,2) = phas * cmplx(alpha(vacuum%nmzxy,1,1,irec2),alpha(vacuum%nmzxy,1,2,irec2))
     109              :             else
     110            0 :                alphm(irec2,2) = cmplx( alpha(vacuum%nmzxy,2,1,irec2), alpha(vacuum%nmzxy,2,2,irec2) )
     111              :             end if
     112            0 :             do ivac = 1, vacuum%nvac
     113            0 :                z = cell%z1
     114            0 :                if ( ivac == 1 ) alph0 = alphm(irec2,2)
     115            0 :                if ( ivac == 2 ) alph0 = alphm(irec2,1)
     116            0 :                do imz = 1, ncsh-1
     117            0 :                   betaz  = cmplx( beta(imz,ivac,1,irec2), beta(imz,ivac,2,irec2) )
     118            0 :                   alphaz = cmplx( alpha(ncsh-imz,ivac,1,irec2), alpha(ncsh-imz,ivac,2,irec2) )
     119            0 :                   e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
     120            0 :                   e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
     121            0 :                   test = e_m * ( alph0 + betaz ) + e_p * alphaz
     122            0 :                   if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
     123            0 :                   vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + vcons * test
     124            0 :                   if ( allocated( field%efield%c1 ) ) then
     125            0 :                      e_m = exp_safe( - g * z )
     126            0 :                      e_p = exp_safe(   g * z )
     127            0 :                      if ( ivac == 1 ) then ! z > 0
     128            0 :                         vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + field%efield%c1(irec2-1) * e_p + field%efield%c2(irec2-1) * e_m
     129              :                      else ! z < 0
     130            0 :                         vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + field%efield%c1(irec2-1) * e_m + field%efield%c2(irec2-1) * e_p
     131              :                      end if
     132              :                   end if
     133            0 :                   z = z + vacuum%delz
     134              :                end do 
     135              :             end do 
     136              :          end do 
     137              :       else ! Neumann
     138        13450 :          do irec2 = start_star, stars%ng2
     139        13413 :             g = stars%sk2(irec2)
     140        13413 :             vcons = tpi_const / g
     141              :      
     142        28513 :             do ivac = 1, vacuum%nvac
     143        15063 :                z = cell%z1
     144      1521363 :                do imz = 1, vacuum%nmzxy
     145      3012600 :                   e_m = exp_safe( - g * z )
     146      3012600 :                   e_p = exp_safe(   g * z )
     147      1506300 :                   fra(ip-imz) = real(  rhtxy(imz,irec2,ivac) ) * e_m
     148      1506300 :                   fia(ip-imz) = aimag( rhtxy(imz,irec2,ivac) ) * e_m
     149      1506300 :                   frb(imz)    = real(  rhtxy(imz,irec2,ivac) ) * e_p
     150      1506300 :                   fib(imz)    = aimag( rhtxy(imz,irec2,ivac) ) * e_p
     151      1521363 :                   z = z + vacuum%delz
     152              :                end do 
     153              :                ! add external field, if segmented
     154        15063 :                if ( allocated( field%efield%rhoef ) ) then
     155            0 :                   z = cell%z1 + field%efield%zsigma
     156            0 :                   e_m = exp_safe( - g * z )
     157            0 :                   e_p = exp_safe(   g * z )
     158              :                   ! the equation has a minus sign as "rhtxy" contains the electron density
     159              :                   ! (a positive number representing a negative charge) while rhoef
     160              :                   ! specifies the charges in terms of the (positive) elementary charge "e".
     161            0 :                   fra(ip-ncsh) = fra(ip-ncsh) - real(  field%efield%rhoef(irec2-1,ivac) ) * e_m
     162            0 :                   fia(ip-ncsh) = fia(ip-ncsh) - aimag( field%efield%rhoef(irec2-1,ivac) ) * e_m
     163            0 :                   frb(ncsh)    = frb(ncsh)    - real(  field%efield%rhoef(irec2-1,ivac) ) * e_p
     164            0 :                   fib(ncsh)    = fib(ncsh)    - aimag( field%efield%rhoef(irec2-1,ivac) ) * e_p
     165              :                end if
     166        15063 :                call intgz1( fra, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,1,irec2), .true. )
     167        15063 :                call intgz1( fia, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,2,irec2), .true. )
     168        15063 :                call qsf( vacuum%delz, frb, beta(1,ivac,1,irec2), vacuum%nmzxy, 1 )
     169        28476 :                call qsf( vacuum%delz, fib, beta(1,ivac,2,irec2), vacuum%nmzxy, 1 )
     170              :             end do 
     171              :          end do
     172              :        
     173        13450 :          do irec2 = start_star, stars%ng2
     174        13413 :             alphm(irec2,1) = cmplx( alpha(vacuum%nmzxy,1,1,irec2), alpha(vacuum%nmzxy,1,2,irec2) )
     175        13450 :             if ( vacuum%nvac == 1 ) then
     176        11763 :                call stars%map_2nd_vac(vacuum,irec2,irec2r,phas)
     177        11763 :                alphm(irec2r,2) = phas * cmplx(alpha(vacuum%nmzxy,1,1,irec2),alpha(vacuum%nmzxy,1,2,irec2))
     178              :             else
     179         1650 :                alphm(irec2,2) = cmplx( alpha(vacuum%nmzxy,2,1,irec2), alpha(vacuum%nmzxy,2,2,irec2) )
     180              :             end if
     181              :          end do
     182              : 
     183        13450 :          do irec2 = start_star, stars%ng2
     184        13413 :             g = stars%sk2(irec2)
     185        13413 :             vcons = tpi_const / g
     186              :               
     187        28513 :             do ivac = 1, vacuum%nvac
     188        15063 :                z = cell%z1
     189        15063 :                if ( ivac == 1 ) alph0 = alphm(irec2,2)
     190        15063 :                if ( ivac == 2 ) alph0 = alphm(irec2,1)
     191      1534776 :                do imz = 1, vacuum%nmzxy
     192      1506300 :                   betaz  = cmplx( beta(imz,ivac,1,irec2), beta(imz,ivac,2,irec2) )
     193      1506300 :                   alphaz = cmplx( alpha(ip-imz,ivac,1,irec2), alpha(ip-imz,ivac,2,irec2) )
     194      3012600 :                   e_m = exp_safe( - g * z )
     195      3012600 :                   e_p = exp_safe(   g * z )
     196      1506300 :                   test = e_m * ( alph0 + betaz ) + e_p * alphaz
     197      1506300 :                   if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
     198      1506300 :                   vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + vcons * test
     199      1521363 :                   z = z + vacuum%delz
     200              :                end do 
     201              :             end do 
     202              :          end do
     203              :       end if
     204           37 :    end subroutine vvacxy
     205              : 
     206      6025200 :    pure real function exp_safe( x )
     207              :       ! replace exp by a function that does not under/overflow
     208              :       implicit none
     209              : 
     210              :       real, intent(in) :: x
     211              :       real, parameter  :: maxexp = log( 2.0 ) * maxexponent( 2.0 )
     212              :       real, parameter  :: minexp = log( 2.0 ) * minexponent( 2.0 )
     213              : 
     214      6025200 :       if ( abs( x ) > minexp .and. abs( x ) < maxexp ) then
     215      6025200 :          exp_safe = exp( x )
     216              :       else
     217            0 :          if ( x > 0 ) then
     218              :             if ( x > minexp ) then
     219              :                exp_safe = exp( maxexp )
     220              :             else
     221              :                exp_safe = exp( minexp )
     222              :             end if
     223              :          else
     224              :             if ( - x > minexp ) then
     225              :                exp_safe = exp( - maxexp )
     226              :             else
     227              :                exp_safe = exp( - minexp )
     228              :             end if
     229              :          end if
     230              :       end if
     231            0 :    end function exp_safe
     232              : end module m_vvacxy
        

Generated by: LCOV version 2.0-1