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

            Line data    Source code
       1              : MODULE m_vintcz
       2              :   !     *************************************************************
       3              :   !     z-dependent part of Coulomb potential in the interstitial   *
       4              :   !     [-d/2,d/2] region            c.l.fu, r.podloucky            *
       5              :   !     *************************************************************
       6              :   !     modified for thick films to avoid underflows gb`06
       7              :   !---------------------------------------------------------------
       8              : CONTAINS
       9       518985 :    COMPLEX FUNCTION vintcz(stars,vacuum,cell,input,field,z,nrec2,psq,vnew,rhobar,sig1dh,vz1dh,alphm,vslope,sigma_disc,l_dfptvgen,l_corr,diff_vmz1dh,sigma_disc2)
      10              :       USE m_constants
      11              :       USE m_types
      12              : 
      13              :       IMPLICIT NONE
      14              : 
      15              :       TYPE(t_stars),  INTENT(IN) :: stars
      16              :       TYPE(t_vacuum), INTENT(IN) :: vacuum
      17              :       TYPE(t_cell),   INTENT(IN) :: cell
      18              :       TYPE(t_input),  INTENT(IN) :: input
      19              :       TYPE(t_field),  INTENT(IN) :: field
      20              : 
      21              :       INTEGER,        INTENT(IN) :: nrec2
      22              :       COMPLEX,        INTENT(IN) :: rhobar,vslope
      23              :       COMPLEX,        INTENT(IN) :: sig1dh,vz1dh
      24              :       REAL,           INTENT(IN) :: z
      25              : 
      26              :       COMPLEX,        INTENT(IN) :: psq(stars%ng3),vnew(:,:,:)!,vxy(:,:,:) !(vacuum%nmzxyd,stars%ng2-1,2)
      27              :       COMPLEX,        INTENT(IN) :: alphm(stars%ng2,2)
      28              :       complex,        intent(in) :: sigma_disc(2)
      29              :       LOGICAL,        INTENT(IN) :: l_dfptvgen, l_corr
      30              :       COMPLEX,        INTENT(IN) :: diff_vmz1dh
      31              : 
      32              :       complex, optional, intent(in) :: sigma_disc2(2)
      33              :       !REAL,           INTENT(IN) :: vz(:,:) !(vacuum%nmzd,2,jspins)
      34              : 
      35              :       COMPLEX                    :: argr,sumrr,vcons1,test,c_ph,phas
      36              :       REAL                       :: bj0,dh,fit,g,g3,q,qdh,vcons2,zf
      37              :       REAL                       :: e_m,e_p,cos_q,sin_q
      38              :       INTEGER                    :: ig3n,im,iq,ivac,k1,k2,nrec2r
      39              :       COMPLEX                    :: newdp, newdm, newdp2, newdm2
      40              :       
      41       518985 :       dh = cell%z1
      42       518985 :       sumrr = (0.,0.)
      43       518985 :       vintcz = (0.,0.)
      44              : 
      45       518985 :       newdp = CMPLX(0.0,0.0)
      46       518985 :       newdm = CMPLX(0.0,0.0)
      47       518985 :       newdp2 = CMPLX(0.0,0.0)
      48       518985 :       newdm2 = CMPLX(0.0,0.0)
      49              :       
      50              :       !--->    if z is in the vacuum, use vacuum representations (m.w.)
      51       518985 :       IF (ABS(z).GE.cell%z1) THEN
      52       177215 :          ivac = 1
      53       177215 :          IF (z.LT.0.0) THEN
      54        87212 :             ivac = 2
      55       177215 :             IF (vacuum%nvac==1) ivac = 1
      56              :          END IF
      57       177215 :          zf = (ABS(z)-cell%z1)/vacuum%delz + 1.0
      58       177215 :          im = zf
      59       177215 :          q = zf - im
      60              :          ! For q/=0 in DFPT, the is no G+q=0, so all stars are treated in the G/=0 way.
      61       708860 :          IF (nrec2.EQ.1.AND.((.NOT.l_dfptvgen).OR.norm2(stars%center)<1e-8)) THEN
      62              :             fit = 0.5* (q-1.)* (q-2.)*REAL(vnew(im,1,ivac)) -&
      63              :                &            q* (q-2.)*REAL(vnew(im+1,1,ivac)) +&
      64          532 :                &        0.5*q* (q-1.)*REAL(vnew(im+2,1,ivac))
      65          532 :             vintcz = CMPLX(fit,0.0)
      66       176683 :          ELSE IF (im+2.LE.vacuum%nmzxy) THEN
      67       176683 :             if (z<0) THEN 
      68        86955 :                call stars%map_2nd_vac(vacuum,nrec2,nrec2r,phas) ! TODO: AN TB; will this work?
      69              :             else
      70        89728 :                nrec2r=nrec2 
      71        89728 :                phas=cmplx(1.0,0.0)
      72              :             end if      
      73              :             vintcz = phas*0.5* (q-1.)* (q-2.)*vnew(im,nrec2r,ivac) -&
      74              :                                     q* (q-2.)*vnew(im+1,nrec2r,ivac) +&
      75       176683 :                                 0.5*q* (q-1.)*vnew(im+2,nrec2r,ivac)
      76              :          END IF
      77       177215 :          RETURN
      78              :       END IF
      79              : 
      80              :       ! For q/=0 in DFPT, there is no G+q=0, so all stars are treated in the G/=0 way.
      81      1367080 :       IF (nrec2==1.AND.((.NOT.l_dfptvgen).OR.norm2(stars%center)<1e-8)) THEN    !     ---->    g=0 coefficient
      82              : 
      83              :          ! New discontinuity correction
      84         1259 :          IF (l_corr) THEN
      85            0 :             newdp = - psq(1) * dh
      86            0 :             newdp2 = -psq(1) * dh**2 
      87              :          END IF
      88              : 
      89        55232 :          DO  iq = -stars%mx3,stars%mx3
      90        53973 :             IF (iq.EQ.0) CYCLE
      91        52714 :             ig3n = stars%ig(0,0,iq)
      92              :             !     ----> use only stars within the g_max sphere (oct.97 shz)
      93        53973 :             IF (ig3n.NE.0) THEN
      94        50196 :                q = iq*cell%bmat(3,3)
      95        50196 :                qdh = q*dh
      96        50196 :                bj0 = SIN(qdh)/qdh
      97        50196 :                argr = ImagUnit*q*z
      98        50196 :                sumrr = (EXP(argr)-EXP(ImagUnit*qdh))/ (q*q) + ImagUnit*COS(qdh)* (dh-z)/q + bj0* (z*z-dh*dh)/2.
      99        50196 :                vintcz = vintcz + fpi_const*psq(ig3n)*sumrr
     100              :                
     101              :                ! New discontinuity correction
     102        50196 :                IF (l_corr) THEN
     103            0 :                   newdp = newdp + ImagUnit * psq(ig3n) * cmplx(cos(qdh), sin(qdh)) * dh / qdh
     104            0 :                   newdp2 = newdp2 + psq(ig3n) * cmplx(cos(qdh), sin(qdh)) / q**2
     105              :                END IF
     106              :             END IF
     107              :          END DO
     108              :          !           -----> v2(z)
     109              :          vintcz = vintcz + vz1dh - fpi_const* (dh-z)*&
     110         1259 :             &              (sig1dh-rhobar/2.* (dh-z))
     111              :          ! Discontinuity correction
     112         1259 :          vintcz = vintcz - fpi_const * sigma_disc(1)*(dh-z)
     113              :          !IF (PRESENT(sigma_disc2)) vintcz = vintcz + fpi_const * sigma_disc2(1)
     114              : 
     115              :          ! New discontinuity correction
     116         1259 :          IF (l_dfptvgen.AND.l_corr) vintcz = vintcz - fpi_const * newdp * (dh-z) + fpi_const * newdp2
     117              : 
     118              :          ! Correct the interstitial potential by the mismatch we calculated as
     119              :          ! a boundary condition for DFPT q=0.
     120              :          !IF (l_dfptvgen) vintcz = vintcz - diff_vmz1dh*z/(2*dh) + diff_vmz1dh/2
     121              : 
     122         1259 :          IF (field%efield%dirichlet .AND. vslope /= 0.0) THEN
     123            0 :             vintcz = vintcz + vslope * (dh-z)
     124              :          END IF
     125              :          !     ---->    (g.ne.0)  coefficients
     126              :       ELSE
     127       340511 :          k1 = stars%kv2(1,nrec2)
     128       340511 :          k2 = stars%kv2(2,nrec2)
     129              : 
     130       340511 :          newdp = cmplx(0.0,0.0)
     131       340511 :          newdm = cmplx(0.0,0.0)
     132       340511 :          newdp2 = cmplx(0.0,0.0)
     133       340511 :          newdm2 = cmplx(0.0,0.0)
     134              : 
     135     10901874 :          DO  iq = -stars%mx3,stars%mx3
     136     10561363 :             ig3n = stars%ig(k1,k2,iq)
     137              :             !     ----> use only stars within the g_max sphere (oct.97 shz)
     138     10901874 :             IF (ig3n.NE.0) THEN
     139      6706139 :                c_ph = stars%rgphs(k1,k2,iq)
     140              :                !           -----> v3(z)
     141      6706139 :                q = iq*cell%bmat(3,3)
     142      6706139 :                g = stars%sk2(nrec2)
     143      6706139 :                g3 = stars%sk3(ig3n)
     144      6706139 :                vcons1 = fpi_const*psq(ig3n)*c_ph / (g3*g3)
     145      6706139 :                IF (field%efield%dirichlet) THEN
     146            0 :                   e_m = 0.0
     147            0 :                   e_p = 0.0
     148            0 :                   vcons1  = vcons1/(g*SINH(g*2*(dh+field%efield%zsigma)))
     149            0 :                   e_m = e_m - EXP(-ImagUnit*q*dh)*(g*COSH(g*(-field%efield%zsigma))+ImagUnit*q*SINH(g*(-field%efield%zsigma)))
     150            0 :                   e_p = e_p + EXP(ImagUnit*q*dh)*(-g*COSH(g*(-field%efield%zsigma))+ImagUnit*q*SINH(g*(-field%efield%zsigma)))
     151            0 :                   sumrr = EXP(ImagUnit*q*z)
     152            0 :                   e_m = e_m + sumrr*(g*COSH(g*(z+dh+field%efield%zsigma))-ImagUnit*q*SINH(g*(z+dh+field%efield%zsigma)))
     153            0 :                   e_p = e_p + sumrr*(g*COSH(g*(z-dh-field%efield%zsigma))-ImagUnit*q*SINH(g*(z-dh-field%efield%zsigma)))
     154              :                   vintcz = vintcz+ vcons1*(e_m*SINH(g*(field%efield%zsigma+dh-z))&
     155            0 :                                           +e_p*SINH(g*(field%efield%zsigma+dh+z)))
     156              :                ELSE
     157      6706139 :                   sumrr = (0.0,0.0)
     158      6706139 :                   vcons2 = - 1.0 / (2.*g)
     159     13412278 :                   e_m = exp_safe( -g*(z+dh) ) !exp_safe handles overflow
     160     13412278 :                   e_p = exp_safe( g*(z-dh) )
     161              :                   !               --> sum over gz-stars
     162      6706139 :                   cos_q = COS(q*dh)
     163      6706139 :                   sin_q = SIN(q*dh)
     164              :                   sumrr = sumrr + CMPLX(COS(q*z),SIN(q*z)) + vcons2 *&
     165              :                         &                      ( (g + ImagUnit*q) * e_p * (cos_q + ImagUnit*sin_q) +&
     166      6706139 :                         &                        (g - ImagUnit*q) * e_m * (cos_q - ImagUnit*sin_q) )
     167      6706139 :                   vintcz = vintcz + vcons1*sumrr
     168              : 
     169              :                   ! New discontinuity correction
     170      6706139 :                   IF (iq==0.AND.l_corr) newdp = newdp - psq(ig3n) * dh
     171              :                   IF (iq/=0.AND.l_corr) newdp = newdp + ImagUnit * psq(ig3n) * cmplx(cos_q, sin_q) * dh / (q*dh)
     172              :                   IF (iq==0.AND.l_corr) newdm = newdm - psq(ig3n) * dh
     173              :                   IF (iq/=0.AND.l_corr) newdm = newdm - ImagUnit * psq(ig3n) * cmplx(cos_q,-sin_q) * dh / (q*dh)
     174      6365628 :                   IF (iq==0.AND.l_corr) newdp2 = newdp2 - psq(ig3n) * dh**2
     175              :                   IF (iq/=0.AND.l_corr) newdp2 = newdp2 + psq(ig3n) * cmplx(cos_q, sin_q) / q**2
     176      6706139 :                   IF (iq==0.AND.l_corr) newdm2 = newdm2 + psq(ig3n) * dh**2
     177      6365628 :                   IF (iq/=0.AND.l_corr) newdm2 = newdm2 - psq(ig3n) * cmplx(cos_q,-sin_q) / q**2
     178              : 
     179              :                END IF ! Neumann (vs. Dirichlet)
     180              :             END IF ! ig3d /= 0
     181              :          END DO
     182              :          !  ----> v4(z)
     183       340511 :          IF (field%efield%dirichlet) THEN
     184            0 :             e_m = SINH(g*(field%efield%zsigma+dh - z))
     185            0 :             e_p = SINH(g*(field%efield%zsigma+dh + z))
     186            0 :             test = e_m*alphm(nrec2,2) + e_p*alphm(nrec2,1)
     187            0 :             test = fpi_const/(g*SINH(g*2*(field%efield%zsigma+dh))) * test
     188            0 :             IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
     189            0 :             vintcz = vintcz + test
     190            0 :             IF (ALLOCATED (field%efield%C1)) THEN
     191            0 :                g = stars%sk2(nrec2)
     192            0 :                e_m = exp_safe (-g*z)
     193            0 :                e_p = exp_safe ( g*z)
     194              :                vintcz = vintcz + field%efield%C1(nrec2-1)*e_m &
     195            0 :                   &            + field%efield%C2(nrec2-1)*e_p
     196              :             END IF
     197              :          ELSE ! Neumann
     198       681022 :             e_m = exp_safe( -g*z  )
     199       681022 :             e_p = exp_safe( g*z  )
     200       340511 :             test = e_m*alphm(nrec2,2) + e_p*alphm(nrec2,1)
     201       340511 :             IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
     202       340511 :             vintcz = vintcz + tpi_const/g* test
     203              : 
     204              :             ! New discontinuity correction
     205              :             !IF (l_dfptvgen.AND.l_corr.AND.nrec2==1) THEN
     206              :             !   e_m = exp_safe( - g * abs(z-cell%z1) )
     207              :             !   e_p = exp_safe( - g * abs(z+cell%z1) )
     208              :             !   test = e_m * newdp + e_p * newdm
     209              :             !   if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
     210              :             !   vintcz = vintcz + tpi_const/g * test
     211              :             !   !IF (abs(z-cell%z1)<1e-9) e_m = 0.0
     212              :             !   !IF (abs(z+cell%z1)<1e-9) e_p = 0.0
     213              :             !   test = e_m * newdp2 * sign(1.0,z-cell%z1)+ e_p * newdm2 * sign(1.0,z+cell%z1)
     214              :             !   if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
     215              :             !   vintcz = vintcz - tpi_const * test
     216              :             !END IF
     217              : 
     218              :          END IF
     219              :       END IF
     220              :    END FUNCTION vintcz
     221              : 
     222     14093300 :    PURE REAL FUNCTION exp_safe(x)
     223              :       ! replace exp by a function that does not under/overflow dw09
     224              :       IMPLICIT NONE
     225              : 
     226              :       REAL, INTENT(IN) :: x
     227              :       REAL, PARAMETER  :: maxexp = LOG(2.0)*MAXEXPONENT(2.0)
     228              :       REAL, PARAMETER  :: minexp = LOG(2.0)*MINEXPONENT(2.0)
     229              : 
     230     14093300 :       IF ( ABS(x)>minexp .AND. ABS(x)<maxexp ) THEN
     231     14093300 :          exp_safe = EXP(x)
     232              :       ELSE
     233            0 :          IF ( x > 0 ) THEN
     234              :             IF ( x > minexp ) THEN
     235              :                exp_safe = EXP(maxexp)
     236              :             ELSE
     237              :                exp_safe = EXP(minexp)
     238              :             END IF
     239              :          ELSE
     240              :             IF ( -x > minexp ) THEN
     241              :                exp_safe = EXP(-maxexp)
     242              :             ELSE
     243              :                exp_safe = EXP(-minexp)
     244              :             END IF
     245              :          END IF
     246              :       END IF
     247            0 :    END FUNCTION exp_safe
     248              : END MODULE m_vintcz
        

Generated by: LCOV version 2.0-1