LCOV - code coverage report
Current view: top level - vgen - vvacis.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 39 51 76.5 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             : module m_vvacis
       2             :   !     **********************************************************
       3             :   !     g.ne.0 coefficients of vacuum coulomb potential          *
       4             :   !     due to the interstitial charge density inside slab       *
       5             :   !                                   c.l.fu, r.podloucky        *
       6             :   !     **********************************************************
       7             :   !     modified for thick films to avoid underflows gb`06
       8             :   !---------------------------------------------------------------
       9             :   contains
      10             : 
      11           5 :   subroutine vvacis( stars, vacuum, sym, cell, psq, input, field, vxy )
      12             : 
      13             :     use m_constants
      14             :     use m_types
      15             :     implicit none
      16             : 
      17             :     type(t_input),  intent(in)  :: input
      18             :     type(t_field),  intent(in)  :: field
      19             :     type(t_vacuum), intent(in)  :: vacuum
      20             :     type(t_sym),    intent(in)  :: sym
      21             :     type(t_stars),  intent(in)  :: stars
      22             :     type(t_cell),   intent(in)  :: cell
      23             :     complex,        intent(in)  :: psq(stars%ng3)
      24             :     complex,        intent(out) :: vxy(vacuum%nmzxyd,stars%ng2-1,2)
      25             : 
      26             :     complex                     :: arg, c_ph, sumr(2)
      27             :     real                        :: dh, g, qz, sign, signz, vcons, z, e_m, arg_r, arg_i
      28             :     integer                     :: i2d, ig3n, imz, imzxy, ivac, k1, k2, kz, m0, nrec2, nrz, nz
      29             :     intrinsic exp
      30             : 
      31          15 :     vxy(:,:,:) = cmplx( 0., 0. )
      32           5 :     dh = cell%z1
      33           5 :     m0 = -stars%mx3
      34           5 :     if ( sym%zrfs ) m0 = 0
      35        1125 :     do  nrec2 = 2, stars%ng2
      36        1120 :       k1 = stars%kv2(1,nrec2)
      37        1120 :       k2 = stars%kv2(2,nrec2)
      38        1120 :       g = stars%sk2(nrec2)
      39        1120 :       if ( field%efield%dirichlet ) then
      40           0 :         vcons = 2.0 * tpi_const / ( g * sinh( g * 2.0 * ( field%efield%zsigma + dh ) ) )
      41           0 :         arg_r = g * ( dh + field%efield%zsigma + dh )
      42             :       else ! neumann
      43        1120 :         vcons = tpi_const / g
      44        1120 :         arg_r = exp_save( - 2 * dh * g )
      45             :       end if
      46        2245 :       do ivac = 1, vacuum%nvac
      47        1120 :         sumr(ivac) = ( 0.0, 0.0 )
      48        1120 :         sign = 3. - 2. * ivac
      49       43808 :         do kz = m0, stars%mx3
      50       42688 :           ig3n = stars%ig(k1,k2,kz)
      51             :           ! use only stars within the g_max sphere
      52       43808 :           if ( ig3n /= 0 ) then
      53       27688 :             c_ph = stars%rgphs(k1,k2,kz)
      54       27688 :             nz = 1
      55       27688 :             if (sym%zrfs) nz = stars%nstr(ig3n) / stars%nstr2(nrec2)
      56       27688 :             qz = kz * cell%bmat(3,3)
      57             :             ! sum over gz-stars
      58       58896 :             do  nrz = 1, nz
      59       31208 :               signz = 3. - 2. * nrz
      60       58896 :               if ( field%efield%dirichlet ) then
      61             :                 ! prefactor
      62           0 :                 arg = exp( - ImagUnit * signz * qz * dh ) / ( 2 * ( g ** 2 + qz ** 2 ) ) * psq(ig3n)
      63           0 :                 if ( ivac == 1 ) then
      64             :                   sumr(ivac) = sumr(ivac) + c_ph * exp( - arg_r ) * arg * ( &                           ! c_ph not tested in this case
      65             :                       ( - exp( 2 * g * ( field%efield%zsigma + dh ) ) + exp( 2 * ( ImagUnit * signz * qz * dh + arg_r ) ) ) * ( g - ImagUnit * signz * qz ) &
      66           0 :                     + ( - exp( 2 * g * dh ) + exp( 2 * ImagUnit * signz * qz * dh ) )                                       * ( g + ImagUnit * signz * qz ) )
      67             :                 else
      68             :                   sumr(ivac) = sumr(ivac) + c_ph * arg * ( &
      69             :                     exp(   arg_r ) * ( g + ImagUnit * signz * qz ) &
      70             :                     + exp( - arg_r ) * ( g - ImagUnit * signz * qz ) &
      71             :                     + 2 * exp( 2 * ImagUnit * signz * qz * dh ) &
      72             :                     * ( - g * cosh( g * ( - field%efield%zsigma ) ) &
      73           0 :                       + ImagUnit * signz * qz * sinh( - g * field%efield%zsigma ) ) )
      74             :                 end if
      75             :               else
      76       31208 :                 arg = g + sign * ImagUnit * signz * qz
      77       31208 :                 arg_i = sign * signz * qz * dh
      78       31208 :                 sumr(ivac) = sumr(ivac) + c_ph * psq(ig3n) * ( cos( arg_i ) * ( 1 - arg_r ) + ImagUnit * sin( arg_i ) * ( 1 + arg_r ) ) / arg
      79             :               end if
      80             :             enddo
      81             :           endif
      82             :         enddo
      83        1120 :         z = 0 ! moved cell%z1 into above equations
      84      114240 :         do imz = 1, vacuum%nmzxy
      85      112000 :           if ( field%efield%dirichlet ) then
      86           0 :             e_m = sinh( g * ( field%efield%zsigma - z ) )
      87             :           else ! neumann
      88      112000 :             e_m = exp_save( - g * z )
      89             :           end if
      90      112000 :           vxy(imz,nrec2-1,ivac) = vxy(imz,nrec2-1,ivac) + vcons * sumr(ivac) * e_m
      91      113120 :           z = z + vacuum%delz
      92             :         enddo
      93             :       enddo
      94             :     enddo
      95             : 
      96           5 :   end subroutine vvacis
      97             : 
      98             : 
      99      113120 :   pure real function exp_save(x)
     100             :     ! replace exp by a function that does not under/overflow
     101             : 
     102             :     implicit none
     103             :     real, intent(in) :: x
     104             :     real, parameter  :: maxexp = log( 2.0 ) * maxexponent( 2.0 )
     105             :     real, parameter  :: minexp = log( 2.0 ) * minexponent( 2.0 )
     106             : 
     107      113120 :     if ( abs(x) > minexp .and. abs(x) < maxexp ) then
     108      113120 :       exp_save = exp( x )
     109             :     else
     110           0 :       if ( x > 0 ) then
     111           0 :         if ( x > minexp ) then
     112             :           exp_save = exp( maxexp )
     113             :         else
     114           0 :           exp_save = exp( minexp )
     115             :         endif
     116             :       else
     117           0 :         if ( -x > minexp ) then
     118             :           exp_save = exp( -maxexp )
     119             :         else
     120           0 :           exp_save = exp( -minexp )
     121             :         endif
     122             :       endif
     123             :     endif
     124      113120 :   end function exp_save
     125             : 
     126             : end module m_vvacis

Generated by: LCOV version 1.13