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

          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       70950 :   COMPLEX FUNCTION vintcz(&
      10             :        &                        stars,vacuum,cell,sym,input,field,&
      11       70950 :        &                        z,nrec2,psq,vxy,vz,rhobar,sig1dh,vz1dh,alphm)
      12             : 
      13             :     USE m_constants
      14             :     USE m_types
      15             :     IMPLICIT NONE
      16             :     !     ..
      17             :     !     .. Scalar Arguments ..
      18             :     TYPE(t_stars),INTENT(IN)  :: stars
      19             :     TYPE(t_vacuum),INTENT(IN) :: vacuum
      20             :     TYPE(t_cell),INTENT(IN)   :: cell
      21             :     TYPE(t_sym),INTENT(IN)    :: sym
      22             :     TYPE(t_input),INTENT(IN)  :: input
      23             :     TYPE(t_field),INTENT(IN)  :: field
      24             :     INTEGER, INTENT (IN) :: nrec2
      25             :     COMPLEX, INTENT (IN) :: rhobar
      26             :     REAL,    INTENT (IN) :: sig1dh,vz1dh  ,z
      27             :     !     ..
      28             :     !     .. Array Arguments ..
      29             :     COMPLEX, INTENT (IN) :: psq(stars%ng3),vxy(:,:,:) !(vacuum%nmzxyd,stars%ng2-1,2)
      30             :     COMPLEX, INTENT (IN) :: alphm(stars%ng2,2)
      31             :     REAL,    INTENT (IN) :: vz(:,:) !(vacuum%nmzd,2,jspins)  
      32             : 
      33             :     COMPLEX argr,sumrr,vcons1,test,c_ph
      34             :     REAL bj0,dh,fit,g,g3,q,qdh,signz,vcons2,zf
      35             :     REAL e_m,e_p,cos_q,sin_q
      36             :     INTEGER ig3n,im,iq,ivac,k1,k2,m0,nrz,nz
      37             :     !     ..
      38             :     !     .. Intrinsic Functions ..
      39             :     INTRINSIC abs,cmplx,conjg,cos,exp,sin
      40             :     !     ..
      41             :     !
      42             : 
      43       70950 :     dh = cell%z1
      44       70950 :     sumrr = (0.,0.)
      45       70950 :     vintcz = (0.,0.)
      46             :     !--->    if z is in the vacuum, use vacuum representations (m.w.)
      47       70950 :     IF (ABS(z).GE.cell%z1) THEN
      48       22945 :        ivac = 1
      49       22945 :        IF (z.LT.0.0) THEN
      50       10910 :           ivac = 2
      51       10910 :           IF (sym%invs .OR. sym%zrfs) ivac = 1
      52             :        END IF
      53       22945 :        zf = (ABS(z)-cell%z1)/vacuum%delz + 1.0
      54       22945 :        im = zf
      55       22945 :        q = zf - im
      56       22945 :        IF (nrec2.EQ.1) THEN
      57             :           fit = 0.5* (q-1.)* (q-2.)*vz(im,ivac) -&
      58             :                &            q* (q-2.)*vz(im+1,ivac) +&
      59          97 :                &            0.5*q* (q-1.)*vz(im+2,ivac)
      60          97 :           vintcz = CMPLX(fit,0.0)
      61       22848 :        ELSE IF (im+2.LE.vacuum%nmzxy) THEN
      62             :           vintcz = 0.5* (q-1.)* (q-2.)*vxy(im,nrec2-1,ivac) -&
      63             :                &               q* (q-2.)*vxy(im+1,nrec2-1,ivac) +&
      64       22848 :                &               0.5*q* (q-1.)*vxy(im+2,nrec2-1,ivac)
      65       22848 :           IF ((sym%invs.AND. (.NOT.sym%zrfs)) .AND.&
      66        7824 :                &          z.LT.0) vintcz = CONJG(vintcz)
      67             :        END IF
      68             :        RETURN
      69             :     END IF
      70             :     !
      71       48005 :     IF (nrec2.EQ.1) THEN
      72         293 :        m0 = -stars%mx3
      73         293 :        IF (sym%zrfs .OR. sym%invs) m0 = 0
      74             :        !     ---->    g=0 coefficient
      75             :        !           -----> v1(z)
      76       19187 :        DO  iq = m0,stars%mx3
      77        9447 :           IF (iq.EQ.0) CYCLE
      78        9154 :           ig3n = stars%ig(0,0,iq)
      79             :           !     ----> use only stars within the g_max sphere (oct.97 shz)
      80        9447 :           IF (ig3n.NE.0) THEN
      81        8861 :              q = iq*cell%bmat(3,3)
      82        8861 :              nz = stars%nstr(ig3n)
      83        8861 :              sumrr = (0.0,0.0)
      84             :              !      --> sum over gz-stars
      85       26583 :              DO  nrz = 1,nz
      86       17722 :                 signz = 3. - 2.*nrz
      87       17722 :                 q = signz*q
      88       17722 :                 qdh = q*dh
      89       17722 :                 bj0 = SIN(qdh)/qdh
      90       17722 :                 argr = ImagUnit*q*z
      91             :                 sumrr = sumrr + (EXP(argr)-EXP(ImagUnit*qdh))/ (q*q) +&
      92       26583 :                      &                 ImagUnit*COS(qdh)* (dh-z)/q + bj0* (z*z-dh*dh)/2.
      93             :              ENDDO
      94        8861 :              vintcz = vintcz + fpi_const*psq(ig3n)*sumrr
      95             :           ENDIF
      96             :        ENDDO
      97             :        !           -----> v2(z)
      98             :        vintcz = vintcz + vz1dh - fpi_const* (dh-z)*&
      99         293 :             &            (sig1dh-rhobar/2.* (dh-z))
     100         293 :        IF (field%efield%dirichlet .AND. field%efield%vslope /= 0.0) THEN
     101           0 :           vintcz = vintcz + field%efield%vslope * (dh-z)
     102             :        END IF
     103             :        !     ---->    (g.ne.0)  coefficients
     104             :     ELSE
     105       47712 :        m0 = -stars%mx3
     106       47712 :        IF (sym%zrfs) m0 = 0
     107       47712 :        k1 = stars%kv2(1,nrec2)
     108       47712 :        k2 = stars%kv2(2,nrec2)
     109     2284320 :        DO  iq = m0,stars%mx3
     110     2236608 :           ig3n = stars%ig(k1,k2,iq)
     111             :           !     ----> use only stars within the g_max sphere (oct.97 shz)
     112     2284320 :           IF (ig3n.NE.0) THEN
     113     1469736 :              c_ph = stars%rgphs(k1,k2,iq)
     114             :              !           -----> v3(z)
     115     1469736 :              q = iq*cell%bmat(3,3)
     116     1469736 :              g = stars%sk2(nrec2)
     117     1469736 :              g3 = stars%sk3(ig3n)
     118     1469736 :              vcons1 = fpi_const*psq(ig3n)*c_ph / (g3*g3)
     119     1469736 :              nz = 1
     120     1469736 :              IF (sym%zrfs) nz = stars%nstr(ig3n)/stars%nstr2(nrec2)
     121     1469736 :              IF (field%efield%dirichlet) THEN
     122           0 :                 e_m = 0.0
     123           0 :                 e_p = 0.0
     124           0 :                 vcons1  = vcons1/(g*SINH(g*2*(dh+field%efield%zsigma)))
     125           0 :                 loop_vacua: DO nrz = 1,nz
     126           0 :                    signz = 3. - 2.*nrz
     127           0 :                    q = signz*q
     128             :                    e_m = e_m - EXP(-ImagUnit*q*dh)*(g*COSH(g*(-field%efield%zsigma))&
     129           0 :                         &                                  +ImagUnit*q*SINH(g*(-field%efield%zsigma)))
     130             :                    e_p = e_p + EXP(ImagUnit*q*dh)*(-g*COSH(g*(-field%efield%zsigma))&
     131           0 :                         &                                  +ImagUnit*q*SINH(g*(-field%efield%zsigma)))
     132           0 :                    sumrr = EXP(ImagUnit*q*z)
     133             :                    e_m = e_m + sumrr*(g*COSH(g*(z+dh+field%efield%zsigma))&
     134           0 :                         &                              -ImagUnit*q*SINH(g*(z+dh+field%efield%zsigma)))
     135             :                    e_p = e_p + sumrr*(g*COSH(g*(z-dh-field%efield%zsigma))&
     136           0 :                         &                              -ImagUnit*q*SINH(g*(z-dh-field%efield%zsigma)))
     137             :                 END DO loop_vacua
     138             :                 vintcz = vintcz&
     139             :                      &             + vcons1*(e_m*SINH(g*(field%efield%zsigma+dh-z))&
     140           0 :                      &                      +e_p*SINH(g*(field%efield%zsigma+dh+z)))
     141             :              ELSE
     142     1469736 :                 sumrr = (0.0,0.0)
     143     1469736 :                 vcons2 = - 1.0 / (2.*g)
     144     1469736 :                 e_m = exp_SAVE( -g*(z+dh) ) !exp_save handles overflow
     145     1469736 :                 e_p = exp_SAVE( g*(z-dh) )
     146             :                 !               --> sum over gz-stars
     147     3055632 :                 vacua: DO nrz = 1,nz
     148     1585896 :                    signz = 3. - 2.*nrz
     149     1585896 :                    q = signz*q
     150     1585896 :                    cos_q = COS(q*dh)
     151     1585896 :                    sin_q = SIN(q*dh)
     152             :                    sumrr = sumrr + CMPLX(COS(q*z),SIN(q*z)) + vcons2 *&
     153             :                         &                      ( (g + ImagUnit*q) * e_p * (cos_q + ImagUnit*sin_q) +&
     154     3055632 :                         &                        (g - ImagUnit*q) * e_m * (cos_q - ImagUnit*sin_q) )
     155             :                 END DO vacua
     156     1469736 :                 vintcz = vintcz + vcons1*sumrr
     157             :              END IF ! Neumann (vs. Dirichlet)
     158             :           END IF ! ig3d /= 0
     159             :        ENDDO
     160             :        !  ----> v4(z)
     161       47712 :        IF (field%efield%dirichlet) THEN
     162           0 :           e_m = SINH(g*(field%efield%zsigma+dh - z))
     163           0 :           e_p = SINH(g*(field%efield%zsigma+dh + z))
     164           0 :           test = e_m*alphm(nrec2-1,2) + e_p*alphm(nrec2-1,1)
     165           0 :           test = fpi_const/(g*SINH(g*2*(field%efield%zsigma+dh))) * test
     166           0 :           IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
     167           0 :           vintcz = vintcz + test
     168           0 :           IF (ALLOCATED (field%efield%C1)) THEN
     169           0 :              g = stars%sk2(nrec2)
     170           0 :              e_m = exp_save (-g*z)
     171           0 :              e_p = exp_save ( g*z)
     172             :              vintcz = vintcz + field%efield%C1(nrec2-1)*e_m&
     173           0 :                   &                       + field%efield%C2(nrec2-1)*e_p
     174             :           END IF
     175             :        ELSE ! Neumann
     176       47712 :           e_m = exp_save( -g*z  )
     177       47712 :           e_p = exp_save( g*z  )
     178       47712 :           test = e_m*alphm(nrec2-1,2) + e_p*alphm(nrec2-1,1)
     179       47712 :           IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
     180       47712 :           vintcz = vintcz + tpi_const/g* test
     181             :        END IF
     182             : 
     183             :     ENDIF
     184             : 
     185             :   END FUNCTION vintcz
     186             : 
     187     3034896 :   PURE REAL FUNCTION exp_save(x)
     188             :     ! replace exp by a function that does not under/overflow dw09
     189             :     IMPLICIT NONE
     190             :     REAL   ,INTENT(IN)     :: x
     191             :     REAL, PARAMETER ::    maxexp = LOG(2.0)*MAXEXPONENT(2.0)
     192             :     REAL, PARAMETER ::    minexp = LOG(2.0)*MINEXPONENT(2.0)
     193             : 
     194     3034896 :     IF ( ABS(x)>minexp .AND. ABS(x)<maxexp ) THEN
     195     3034896 :        exp_SAVE = EXP(x)
     196             :     ELSE
     197           0 :        IF ( x > 0 ) THEN
     198           0 :           IF ( x > minexp ) THEN
     199             :              exp_save = EXP(maxexp)
     200             :           ELSE
     201           0 :              exp_save = EXP(minexp)
     202             :           ENDIF
     203             :        ELSE
     204           0 :           IF ( -x > minexp ) THEN
     205             :              exp_save = EXP(-maxexp)
     206             :           ELSE
     207           0 :              exp_save = EXP(-minexp)
     208             :           ENDIF
     209             :        ENDIF
     210             :     ENDIF
     211     3034896 :   END FUNCTION exp_save
     212             : END MODULE m_vintcz

Generated by: LCOV version 1.13