vvacis.f90 5.46 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
MODULE m_vvacis
  !     **********************************************************
  !     g.ne.0 coefficients of vacuum coulomb potential          *
  !     due to the interstitial charge density inside slab       *
  !                                   c.l.fu, r.podloucky        *
  !     **********************************************************
  !     modified for thick films to avoid underflows gb`06
  !---------------------------------------------------------------
CONTAINS
  SUBROUTINE vvacis(&
       &                  stars,vacuum,&
       &                  sym,cell,&
13
       &                  psq, input,field,&
14 15 16 17 18 19
       &                  vxy)

    USE m_constants
    USE m_types
    IMPLICIT NONE
    TYPE(t_input),INTENT(IN)   :: input
20
    TYPE(t_field),INTENT(IN)   :: field
21 22 23 24 25 26
    TYPE(t_vacuum),INTENT(IN)  :: vacuum
    TYPE(t_sym),INTENT(IN)     :: sym
    TYPE(t_stars),INTENT(IN)   :: stars
    TYPE(t_cell),INTENT(IN)    :: cell
    !     ..
    !     .. Array Arguments ..
27
    COMPLEX, INTENT (IN) :: psq(stars%ng3)
28
    COMPLEX, INTENT (OUT):: vxy(vacuum%nmzxyd,stars%ng2-1,2)
29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
    !     ..
    !     .. Local Scalars ..
    COMPLEX arg,ci
    REAL dh,g,qz,sign,signz,vcons,z,e_m
    REAL arg_r,arg_i
    INTEGER i2d,ig3n,imz,imzxy,ivac,k1,k2,kz,m0,nrec2,nrz,nz
    !     ..
    !     .. Local Arrays ..
    COMPLEX sumr(2)
    !     ..
    !     .. Intrinsic Functions ..
    INTRINSIC exp
    !     ..
    ci = CMPLX(0.0,1.0)

44
    vxy(:,:,:) = CMPLX(0.,0.)
45 46 47 48 49 50 51
    dh = cell%z1
    m0 = -stars%mx3
    IF (sym%zrfs) m0 = 0
    DO  nrec2 = 2,stars%ng2
       k1 = stars%kv2(1,nrec2)
       k2 = stars%kv2(2,nrec2)
       g = stars%sk2(nrec2)
52 53 54
       IF (field%efield%dirichlet) THEN
          vcons = 2.0*tpi_const/(g*SINH(g*2.0*(field%efield%zsigma+dh)))
          arg_r = g*(dh+field%efield%zsigma+dh)
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
       ELSE ! Neumann
          vcons = tpi_const/g
          arg_r = exp_save( - 2*dh*g )
       END IF
       DO ivac = 1,vacuum%nvac
          sumr(ivac) = (0.0,0.0)
          sign = 3. - 2.*ivac
          DO kz = m0,stars%mx3
             ig3n = stars%ig(k1,k2,kz)
             !     ----> use only stars within the g_max sphere (oct.97 shz)
             IF (ig3n.NE.0) THEN
                nz = 1
                IF (sym%zrfs) nz = stars%nstr(ig3n)/stars%nstr2(nrec2)
                qz = kz*cell%bmat(3,3)
                !     ---> sum over gz-stars
                DO  nrz = 1,nz
                   signz = 3. - 2.*nrz
72
                   IF (field%efield%dirichlet) THEN
73 74 75 76 77
                      ! prefactor
                      arg = EXP(-ci*signz*qz*dh)&
                           &                      /(2*(g**2 + qz**2)) * psq(ig3n)
                      IF (ivac == 1) THEN
                         sumr(ivac) = sumr(ivac) + EXP(-arg_r)*arg*(&
78
                              &                     (- EXP(2*g*(field%efield%zsigma+dh))&
79 80 81 82 83 84 85 86 87 88
                              &                      + EXP(2*(ci*signz*qz*dh+arg_r)))&
                              &                     *(g-ci*signz*qz)&
                              &                    +(- EXP(2*g*dh)&
                              &                      + EXP(2*ci*signz*qz*dh))&
                              &                     *(g+ci*signz*qz) )
                      ELSE
                         sumr(ivac) = sumr(ivac) + arg*(&
                              &                     EXP(arg_r)*(g+ci*signz*qz)&
                              &                     +(g-ci*signz*qz)*EXP(-arg_r)&
                              &                     +2*EXP(2*(ci*signz*qz*dh))&
89 90
                              &                      *(-g*COSH(g*(-field%efield%zsigma))&
                              &                        +ci*signz*qz*SINH(g*(-field%efield%zsigma))) )
91 92 93 94 95 96 97 98 99 100 101 102 103
                      END IF
                   ELSE
                      arg = g + sign*ci*signz*qz
                      arg_i = sign*signz*qz*dh
                      sumr(ivac) = sumr(ivac) + psq(ig3n)*(&
                           &                     COS(arg_i)*( 1 - arg_r ) +&
                           &                  ci*SIN(arg_i)*( 1 + arg_r ) ) / arg
                   END IF
                ENDDO  ! nrz 
             ENDIF
          ENDDO       ! kz 
          z = 0 ! moved cell%z1 into above equations gb`06
          DO imz = 1,vacuum%nmzxy
104 105
             IF (field%efield%dirichlet) THEN
                e_m = SINH(g*(field%efield%zsigma-z))
106 107 108
             ELSE ! NEUMANN
                e_m = exp_save( -g*z  )
             END IF
109
             vxy(imz,nrec2-1,ivac) = vxy(imz,nrec2-1,ivac) +&
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
                  &                                   vcons*sumr(ivac)*e_m
             z = z + vacuum%delz
          ENDDO  ! imz 
       ENDDO     ! ivac
    ENDDO        ! nrec2

  END SUBROUTINE vvacis
  !------------------------------------------------------------------
  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_vvacis