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

          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             : 
      14           5 :   subroutine vvacxy( stars, vacuum, cell, sym, input, field, rhtxy, vxy, alphm )
      15             : 
      16             :     use m_intgr, only: intgz1
      17             :     use m_constants
      18             :     use m_types
      19             :     use m_qsf
      20             :     implicit none
      21             : 
      22             :     type(t_input),  intent(in)     :: input
      23             :     type(t_field),  intent(in)     :: field
      24             :     type(t_vacuum), intent(in)     :: vacuum
      25             :     type(t_sym),    intent(in)     :: sym
      26             :     type(t_stars),  intent(in)     :: stars
      27             :     type(t_cell),   intent(in)     :: cell
      28             :     complex,        intent(in)     :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
      29             :     complex,        intent(inout)  :: vxy(vacuum%nmzxyd,stars%ng2-1,2)
      30             :     complex,        intent(out)    :: alphm(stars%ng2,2)
      31             : 
      32             :     complex                        :: alph0, alph2, alph1, alphaz, betaz, test
      33             :     real                           :: g, vcons, z, e_m, e_p
      34             :     integer                        :: imz, ip, irec2, ivac, ncsh
      35             :     logical                        :: tail
      36          10 :     real                           :: fra(vacuum%nmzxyd), frb(vacuum%nmzxyd), fia(vacuum%nmzxyd), fib(vacuum%nmzxyd)
      37          10 :     real                           :: alpha(vacuum%nmzxyd,2,2), beta(vacuum%nmzxyd,2,2)
      38             :     real, allocatable              :: sig_top(:), sig_bot(:)
      39             :     intrinsic aimag, cmplx, conjg, exp, real
      40             : 
      41           5 :     if ( allocated( field%efield%rhoef ) .or. field%efield%dirichlet ) then
      42           0 :       ncsh = field%efield%zsigma / vacuum%delz + 1.01
      43             :       ! if nmzxy < ncsh, the inhomogenous field cannot be represented.
      44             :       ! if nmzxy > ncsh, the boundary condition is wrong - and the
      45             :       ! potential is very wavy - try it yourself, if you don't believe.
      46           0 :       if ( vacuum%nmzxyd < ncsh ) then
      47           0 :         write (6,*) 'error vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '< ncsh(zsigma) = ', ncsh
      48           0 :         call judft_error( "error: vacuum%nmzxyd < ncsh", calledby="vvacxy" )
      49           0 :       else if ( vacuum%nmzxyd > ncsh ) then
      50           0 :         write (6,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
      51           0 :         write (0,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
      52           0 :         call judft_warn( "nmzxyd > ncsh", calledby="vvacxy" )
      53             :       end if
      54             :     end if
      55             : 
      56             :     ! 2-dim star loop g/=0
      57           5 :     ip = vacuum%nmzxy + 1
      58        1125 :     irec2_loop: do irec2 = 2, stars%ng2
      59        1120 :       g = stars%sk2(irec2)
      60        1120 :       vcons = tpi_const / g
      61             :       ! ********** dirichlet ************************************
      62        1125 :       if ( field%efield%dirichlet ) then
      63           0 :         if ( allocated( field%efield%rhoef ) ) then
      64           0 :           vxy(ncsh:vacuum%nmzxy,irec2-1,1) = field%efield%rhoef(irec2-1, 1)
      65           0 :           if ( vacuum%nvac == 2 ) then
      66           0 :             vxy(ncsh:vacuum%nmzxy,irec2-1,2) = field%efield%rhoef(irec2-1, 2)
      67             :           end if
      68             :         else
      69           0 :           vxy(ncsh:vacuum%nmzxy,irec2-1,1:vacuum%nvac) = 0.0
      70             :         end if
      71           0 :         vcons = 2.0 * vcons / sinh( g * 2.0 * ( field%efield%zsigma + cell%z1 ) )
      72           0 :         ivac_loop1: do ivac = 1, vacuum%nvac
      73           0 :           z = cell%z1
      74           0 :           imz_loop1: do imz = 1, ncsh-1
      75             :             ! as "z" > 0 in this subroutine, the integrand is the same
      76             :             ! for both ivac -- but the integral bounds are reversed
      77           0 :             e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
      78           0 :             e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
      79           0 :             fra(ncsh-imz) = real(  rhtxy(imz,irec2-1,ivac) ) * e_m
      80           0 :             fia(ncsh-imz) = aimag( rhtxy(imz,irec2-1,ivac) ) * e_m
      81           0 :             frb(imz)      = real(  rhtxy(imz,irec2-1,ivac) ) * e_p
      82           0 :             fib(imz)      = aimag( rhtxy(imz,irec2-1,ivac) ) * e_p
      83           0 :             z = z + vacuum%delz
      84             :           end do imz_loop1
      85           0 :           call intgz1( fra, vacuum%delz, ncsh - 1, alpha(1,ivac,1), .false. )
      86           0 :           call intgz1( fia, vacuum%delz, ncsh - 1, alpha(1,ivac,2), .false. )
      87           0 :           call qsf( vacuum%delz, frb, beta(1,ivac,1), ncsh - 1, 1 )
      88           0 :           call qsf( vacuum%delz, fib, beta(1,ivac,2), ncsh - 1, 1 )
      89             :         end do ivac_loop1
      90             : 
      91           0 :         if ( ivac == 2 ) then
      92             :           ! honour reversed integral bounds
      93           0 :           alpha(:,ivac,1) = - alpha(:,ivac,1)
      94           0 :           alpha(:,ivac,2) = - alpha(:,ivac,2)
      95           0 :           beta(:,ivac,1)  = - beta(:,ivac,1)
      96           0 :           beta(:,ivac,2)  = - beta(:,ivac,2)
      97             :         end if
      98             : 
      99           0 :         alph1 = cmplx( alpha(ncsh-1,1,1), alpha(ncsh-1,1,2) )
     100           0 :         if ( vacuum%nvac == 1 ) then
     101           0 :           if ( sym%invs ) then
     102           0 :             alph2 = conjg( alph1 )
     103             :           else
     104             :             alph2 = alph1
     105             :           end if
     106             :         else
     107           0 :           alph2 = cmplx( alpha(ncsh-1,2,1), alpha(ncsh-1,2,2) )
     108             :         end if
     109           0 :         ivac_loop2: do ivac = 1, vacuum%nvac
     110           0 :           z = cell%z1
     111           0 :           if ( ivac == 1 ) alph0 = alph2
     112           0 :           if ( ivac == 2 ) alph0 = alph1
     113           0 :           imz_loop2: do imz = 1, ncsh-1
     114           0 :             betaz  = cmplx( beta(imz,ivac,1), beta(imz,ivac,2) )
     115           0 :             alphaz = cmplx( alpha(ncsh-imz,ivac,1), alpha(ncsh-imz,ivac,2) )
     116           0 :             e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
     117           0 :             e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
     118           0 :             test = e_m * ( alph0 + betaz ) + e_p * alphaz
     119           0 :             if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
     120           0 :             vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) + vcons * test
     121           0 :             if ( allocated( field%efield%c1 ) ) then
     122           0 :               e_m = exp_save( - g * z )
     123           0 :               e_p = exp_save(   g * z )
     124           0 :               if ( ivac == 1 ) then ! z > 0
     125             :                 vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) &
     126             :                 + field%efield%c1(irec2-1) * e_p &
     127           0 :                 + field%efield%c2(irec2-1) * e_m
     128             :               else ! z < 0
     129             :                 vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) &
     130             :                 + field%efield%c1(irec2-1) * e_m &
     131           0 :                 + field%efield%c2(irec2-1) * e_p
     132             :               end if
     133             :             end if
     134           0 :             z = z + vacuum%delz
     135             :           end do imz_loop2
     136             :         end do ivac_loop2
     137           0 :         alphm(irec2-1,1) = alph1
     138           0 :         alphm(irec2-1,2) = alph2
     139             : 
     140             :         ! ********** neumann ************************************
     141             :       else
     142        2240 :         ivac_loop3: do ivac = 1, vacuum%nvac
     143        1120 :           z = cell%z1
     144      113120 :           imz_loop3: do imz = 1, vacuum%nmzxy
     145      112000 :             e_m = exp_save( - g * z )
     146      112000 :             e_p = exp_save(   g * z )
     147      112000 :             fra(ip-imz) = real(  rhtxy(imz,irec2-1,ivac) ) * e_m
     148      112000 :             fia(ip-imz) = aimag( rhtxy(imz,irec2-1,ivac) ) * e_m
     149      112000 :             frb(imz)    = real(  rhtxy(imz,irec2-1,ivac) ) * e_p
     150      112000 :             fib(imz)    = aimag( rhtxy(imz,irec2-1,ivac) ) * e_p
     151      113120 :             z = z + vacuum%delz
     152             :           end do imz_loop3
     153             :           ! add external field, if segmented
     154        1120 :           if ( allocated( field%efield%rhoef ) ) then
     155           0 :             z = cell%z1 + field%efield%zsigma
     156           0 :             e_m = exp_save( - g * z )
     157           0 :             e_p = exp_save(   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        1120 :           call intgz1( fra, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,1), .true. )
     167        1120 :           call intgz1( fia, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,2), .true. )
     168        1120 :           call qsf( vacuum%delz, frb, beta(1,ivac,1), vacuum%nmzxy, 1 )
     169        2240 :           call qsf( vacuum%delz, fib, beta(1,ivac,2), vacuum%nmzxy, 1 )
     170             :         end do ivac_loop3
     171        1120 :         alph1 = cmplx( alpha(vacuum%nmzxy,1,1), alpha(vacuum%nmzxy,1,2) )
     172        1120 :         if ( vacuum%nvac == 1 ) then
     173        1120 :           if ( sym%invs ) then
     174        1120 :             alph2 = conjg( alph1 )
     175             :           else
     176             :             alph2 = alph1
     177             :           end if
     178             :         else
     179           0 :           alph2 = cmplx( alpha(vacuum%nmzxy,2,1), alpha(vacuum%nmzxy,2,2) )
     180             :         end if
     181        2240 :         ivac_loop4: do ivac = 1, vacuum%nvac
     182        1120 :           z = cell%z1
     183        1120 :           if ( ivac == 1 ) alph0 = alph2
     184        1120 :           if ( ivac == 2 ) alph0 = alph1
     185      114240 :           imz_loop4: do imz = 1, vacuum%nmzxy
     186      112000 :             betaz  = cmplx( beta(imz,ivac,1), beta(imz,ivac,2) )
     187      112000 :             alphaz = cmplx( alpha(ip-imz,ivac,1), alpha(ip-imz,ivac,2) )
     188      112000 :             e_m = exp_save( - g * z )
     189      112000 :             e_p = exp_save(   g * z )
     190      112000 :             test = e_m * ( alph0 + betaz ) + e_p * alphaz
     191      112000 :             if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
     192      112000 :             vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) + vcons * test
     193      113120 :             z = z + vacuum%delz
     194             :           end do imz_loop4
     195             :         end do ivac_loop4
     196        1120 :         alphm(irec2-1,1) = alph1
     197        1120 :         alphm(irec2-1,2) = alph2
     198             :       end if
     199             :     end do irec2_loop
     200           5 :   end subroutine vvacxy
     201             : 
     202             :   !------------------------------------------------------------------
     203      448000 :   pure real function exp_save( x )
     204             :     ! replace exp by a function that does not under/overflow
     205             :     implicit none
     206             :     real, intent(in) :: x
     207             :     real, parameter  :: maxexp = log( 2.0 ) * maxexponent( 2.0 )
     208             :     real, parameter  :: minexp = log( 2.0 ) * minexponent( 2.0 )
     209             : 
     210      448000 :     if ( abs( x ) > minexp .and. abs( x ) < maxexp ) then
     211      448000 :       exp_save = exp( x )
     212             :     else
     213           0 :       if ( x > 0 ) then
     214           0 :         if ( x > minexp ) then
     215             :           exp_save = exp( maxexp )
     216             :         else
     217           0 :           exp_save = exp( minexp )
     218             :         endif
     219             :       else
     220           0 :         if ( - x > minexp ) then
     221             :           exp_save = exp( - maxexp )
     222             :         else
     223           0 :           exp_save = exp( - minexp )
     224             :         endif
     225             :       endif
     226             :     endif
     227      448000 :   end function exp_save
     228             : 
     229             : end module m_vvacxy

Generated by: LCOV version 1.13