vintcz.f90 8.16 KB
Newer Older
1 2 3 4 5 6 7 8 9
MODULE m_vintcz
  !     *************************************************************
  !     z-dependent part of Coulomb potential in the interstitial   *
  !     [-d/2,d/2] region            c.l.fu, r.podloucky            *
  !     *************************************************************
  !     modified for thick films to avoid underflows gb`06
  !---------------------------------------------------------------
CONTAINS
  COMPLEX FUNCTION vintcz(&
10
       &                        stars,vacuum,cell,sym,input,field,&
11 12 13 14 15 16 17 18 19 20 21 22
       &                        z,nrec2,psq,vxy,vz,rhobar,sig1dh,vz1dh,alphm)

    USE m_constants
    USE m_types
    IMPLICIT NONE
    !     ..
    !     .. Scalar Arguments ..
    TYPE(t_stars),INTENT(IN)  :: stars
    TYPE(t_vacuum),INTENT(IN) :: vacuum
    TYPE(t_cell),INTENT(IN)   :: cell
    TYPE(t_sym),INTENT(IN)    :: sym
    TYPE(t_input),INTENT(IN)  :: input
23
    TYPE(t_field),INTENT(IN)  :: field
24 25 26 27 28
    INTEGER, INTENT (IN) :: nrec2
    COMPLEX, INTENT (IN) :: rhobar
    REAL,    INTENT (IN) :: sig1dh,vz1dh  ,z
    !     ..
    !     .. 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 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
    REAL bj0,dh,fit,g,g3,q,qdh,signz,vcons2,zf
    REAL e_m,e_p,cos_q,sin_q
    INTEGER ig3n,im,iq,ivac,k1,k2,m0,nrz,nz
    !     ..
    !     .. Intrinsic Functions ..
    INTRINSIC abs,cmplx,conjg,cos,exp,sin
    !     ..
    !

    dh = cell%z1
    sumrr = (0.,0.)
    vintcz = (0.,0.)
    !--->    if z is in the vacuum, use vacuum representations (m.w.)
    IF (ABS(z).GE.cell%z1) THEN
       ivac = 1
       IF (z.LT.0.0) THEN
          ivac = 2
          IF (sym%invs .OR. sym%zrfs) ivac = 1
       END IF
       zf = (ABS(z)-cell%z1)/vacuum%delz + 1.0
       im = zf
       q = zf - im
       IF (nrec2.EQ.1) THEN
57 58 59
          fit = 0.5* (q-1.)* (q-2.)*vz(im,ivac) -&
               &            q* (q-2.)*vz(im+1,ivac) +&
               &            0.5*q* (q-1.)*vz(im+2,ivac)
60 61
          vintcz = CMPLX(fit,0.0)
       ELSE IF (im+2.LE.vacuum%nmzxy) THEN
62 63 64
          vintcz = 0.5* (q-1.)* (q-2.)*vxy(im,nrec2-1,ivac) -&
               &               q* (q-2.)*vxy(im+1,nrec2-1,ivac) +&
               &               0.5*q* (q-1.)*vxy(im+2,nrec2-1,ivac)
65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
          IF ((sym%invs.AND. (.NOT.sym%zrfs)) .AND.&
               &          z.LT.0) vintcz = CONJG(vintcz)
       END IF
       RETURN
    END IF
    !
    IF (nrec2.EQ.1) THEN
       m0 = -stars%mx3
       IF (sym%zrfs .OR. sym%invs) m0 = 0
       !     ---->    g=0 coefficient
       !           -----> v1(z)
       DO  iq = m0,stars%mx3
          IF (iq.EQ.0) CYCLE
          ig3n = stars%ig(0,0,iq)
          !     ----> use only stars within the g_max sphere (oct.97 shz)
          IF (ig3n.NE.0) THEN
             q = iq*cell%bmat(3,3)
             nz = stars%nstr(ig3n)
             sumrr = (0.0,0.0)
             !      --> sum over gz-stars
             DO  nrz = 1,nz
                signz = 3. - 2.*nrz
                q = signz*q
                qdh = q*dh
                bj0 = SIN(qdh)/qdh
90 91 92
                argr = ImagUnit*q*z
                sumrr = sumrr + (EXP(argr)-EXP(ImagUnit*qdh))/ (q*q) +&
                     &                 ImagUnit*COS(qdh)* (dh-z)/q + bj0* (z*z-dh*dh)/2.
93 94 95 96 97 98 99
             ENDDO
             vintcz = vintcz + fpi_const*psq(ig3n)*sumrr
          ENDIF
       ENDDO
       !           -----> v2(z)
       vintcz = vintcz + vz1dh - fpi_const* (dh-z)*&
            &            (sig1dh-rhobar/2.* (dh-z))
100 101
       IF (field%efield%dirichlet .AND. field%efield%vslope /= 0.0) THEN
          vintcz = vintcz + field%efield%vslope * (dh-z)
102 103 104 105 106 107 108 109 110 111 112
       END IF
       !     ---->    (g.ne.0)  coefficients
    ELSE
       m0 = -stars%mx3
       IF (sym%zrfs) m0 = 0
       k1 = stars%kv2(1,nrec2)
       k2 = stars%kv2(2,nrec2)
       DO  iq = m0,stars%mx3
          ig3n = stars%ig(k1,k2,iq)
          !     ----> use only stars within the g_max sphere (oct.97 shz)
          IF (ig3n.NE.0) THEN
113
             c_ph = stars%rgphs(k1,k2,iq)
114 115 116 117
             !           -----> v3(z)
             q = iq*cell%bmat(3,3)
             g = stars%sk2(nrec2)
             g3 = stars%sk3(ig3n)
118
             vcons1 = fpi_const*psq(ig3n)*c_ph / (g3*g3)
119 120
             nz = 1
             IF (sym%zrfs) nz = stars%nstr(ig3n)/stars%nstr2(nrec2)
121
             IF (field%efield%dirichlet) THEN
122 123
                e_m = 0.0
                e_p = 0.0
124
                vcons1  = vcons1/(g*SINH(g*2*(dh+field%efield%zsigma)))
125 126 127
                loop_vacua: DO nrz = 1,nz
                   signz = 3. - 2.*nrz
                   q = signz*q
128 129 130 131 132
                   e_m = e_m - EXP(-ImagUnit*q*dh)*(g*COSH(g*(-field%efield%zsigma))&
                        &                                  +ImagUnit*q*SINH(g*(-field%efield%zsigma)))
                   e_p = e_p + EXP(ImagUnit*q*dh)*(-g*COSH(g*(-field%efield%zsigma))&
                        &                                  +ImagUnit*q*SINH(g*(-field%efield%zsigma)))
                   sumrr = EXP(ImagUnit*q*z)
133
                   e_m = e_m + sumrr*(g*COSH(g*(z+dh+field%efield%zsigma))&
134
                        &                              -ImagUnit*q*SINH(g*(z+dh+field%efield%zsigma)))
135
                   e_p = e_p + sumrr*(g*COSH(g*(z-dh-field%efield%zsigma))&
136
                        &                              -ImagUnit*q*SINH(g*(z-dh-field%efield%zsigma)))
137 138
                END DO loop_vacua
                vintcz = vintcz&
139
                     &             + vcons1*(e_m*SINH(g*(field%efield%zsigma+dh-z))&
140
                     &                      +e_p*SINH(g*(field%efield%zsigma+dh+z)))
141 142 143 144 145 146 147 148 149 150 151 152
             ELSE
                sumrr = (0.0,0.0)
                vcons2 = - 1.0 / (2.*g)
                e_m = exp_SAVE( -g*(z+dh) ) !exp_save handles overflow
                e_p = exp_SAVE( g*(z-dh) )
                !               --> sum over gz-stars
                vacua: DO nrz = 1,nz
                   signz = 3. - 2.*nrz
                   q = signz*q
                   cos_q = COS(q*dh)
                   sin_q = SIN(q*dh)
                   sumrr = sumrr + CMPLX(COS(q*z),SIN(q*z)) + vcons2 *&
153 154
                        &                      ( (g + ImagUnit*q) * e_p * (cos_q + ImagUnit*sin_q) +&
                        &                        (g - ImagUnit*q) * e_m * (cos_q - ImagUnit*sin_q) )
155 156 157 158 159 160
                END DO vacua
                vintcz = vintcz + vcons1*sumrr
             END IF ! Neumann (vs. Dirichlet)
          END IF ! ig3d /= 0
       ENDDO
       !  ----> v4(z)
161 162 163
       IF (field%efield%dirichlet) THEN
          e_m = SINH(g*(field%efield%zsigma+dh - z))
          e_p = SINH(g*(field%efield%zsigma+dh + z))
164
          test = e_m*alphm(nrec2-1,2) + e_p*alphm(nrec2-1,1)
165
          test = fpi_const/(g*SINH(g*2*(field%efield%zsigma+dh))) * test
166 167
          IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
          vintcz = vintcz + test
168
          IF (ALLOCATED (field%efield%C1)) THEN
169 170 171
             g = stars%sk2(nrec2)
             e_m = exp_save (-g*z)
             e_p = exp_save ( g*z)
172 173
             vintcz = vintcz + field%efield%C1(nrec2-1)*e_m&
                  &                       + field%efield%C2(nrec2-1)*e_p
174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212
          END IF
       ELSE ! Neumann
          e_m = exp_save( -g*z  )
          e_p = exp_save( g*z  )
          test = e_m*alphm(nrec2-1,2) + e_p*alphm(nrec2-1,1)
          IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
          vintcz = vintcz + tpi_const/g* test
       END IF

    ENDIF

  END FUNCTION vintcz

  PURE REAL FUNCTION exp_save(x)
    ! replace exp by a function that does not under/overflow dw09
    IMPLICIT NONE
    REAL   ,INTENT(IN)     :: x
    REAL, PARAMETER ::    maxexp = LOG(2.0)*MAXEXPONENT(2.0)
    REAL, PARAMETER ::    minexp = LOG(2.0)*MINEXPONENT(2.0)

    IF ( ABS(x)>minexp .AND. ABS(x)<maxexp ) THEN
       exp_SAVE = EXP(x)
    ELSE
       IF ( x > 0 ) THEN
          IF ( x > minexp ) THEN
             exp_save = EXP(maxexp)
          ELSE
             exp_save = EXP(minexp)
          ENDIF
       ELSE
          IF ( -x > minexp ) THEN
             exp_save = EXP(-maxexp)
          ELSE
             exp_save = EXP(-minexp)
          ENDIF
       ENDIF
    ENDIF
  END FUNCTION exp_save
END MODULE m_vintcz