vvacxy.f90 10.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
MODULE m_vvacxy
  USE m_juDFT
  !     **********************************************************
  !     g.ne.0 coefficient of vacuum coulomb potential           *
  !     due to warped vacuum charged density                     *
  !                                 c.l.fu, r.podloucky          *
  !     **********************************************************
  !     modified for thick films to avoid underflows gb`06
  !---------------------------------------------------------------
CONTAINS
  SUBROUTINE vvacxy(&
12
       &                  stars,vacuum,cell,sym,input,field,&
13 14 15 16 17 18 19 20 21 22
       &                  rhtxy,&
       &                  vxy,&
       &                  alphm)

    USE m_intgr, ONLY : intgz1
    USE m_constants
    USE m_types
    USE m_qsf
    IMPLICIT NONE
    TYPE(t_input),INTENT(IN)       :: input
23
    TYPE(t_field),INTENT(IN)       :: field
24 25 26 27 28 29 30
    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 ..
31 32
    COMPLEX, INTENT (IN)    :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
    COMPLEX, INTENT (INOUT) :: vxy(vacuum%nmzxyd,stars%ng2-1,2)
33
    COMPLEX, INTENT (OUT)   :: alphm(stars%ng2,2)
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
    !     ..
    !     .. Local Scalars ..
    COMPLEX alph0 ,alph2,alph1,alphaz,betaz,test
    REAL g,vcons,z,e_m,e_p
    INTEGER imz,ip,irec2,ivac,ncsh
    LOGICAL tail
    !     ..
    !     .. Local Arrays ..
    REAL fra(vacuum%nmzxyd),frb(vacuum%nmzxyd),fia(vacuum%nmzxyd),fib(vacuum%nmzxyd)
    REAL alpha(vacuum%nmzxyd,2,2),beta(vacuum%nmzxyd,2,2)
    REAL, ALLOCATABLE :: sig_top(:), sig_bot(:)
    !     ..
    !     .. Intrinsic Functions ..
    INTRINSIC aimag,cmplx,conjg,exp,REAL
    !     ..
49 50
    IF (ALLOCATED (field%efield%rhoEF) .OR. field%efield%dirichlet) THEN
       ncsh = field%efield%zsigma/vacuum%delz + 1.01
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
       ! If nmzxy < ncsh, the inhomogenous field cannot be represented.
       ! If nmzxy > ncsh, the boundary condition is wrong - and the
       ! potential is very wavy - try it yourself, if you don't believe.
       IF (vacuum%nmzxyd < ncsh) THEN
          WRITE (6,*) 'ERROR vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd,&
               &                '< ncsh(zsigma) = ',ncsh
          CALL juDFT_error("ERROR: vacuum%nmzxyd < ncsh",calledby ="vvacxy")
       ELSE IF (vacuum%nmzxyd > ncsh) THEN
          WRITE (6,*) 'WARNING vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd,&
               &                '> ncsh(zsigma) = ',ncsh
          WRITE (0,*) 'WARNING vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd,&
               &                '> ncsh(zsigma) = ',ncsh
          CALL juDFT_warn("nmzxyd > ncsh",calledby ="vvacxy")
       END IF
    END IF

    !     ..
    !     2-dim star loop g.ne.0
    ip = vacuum%nmzxy + 1
    irec2_loop: DO irec2 = 2,stars%ng2
       g = stars%sk2(irec2)
       vcons = tpi_const/g

       ! ********** DIRICHLET ************************************
75 76 77
       IF (field%efield%dirichlet) THEN
          IF (ALLOCATED (field%efield%rhoEF)) THEN
             vxy(ncsh:vacuum%nmzxy,irec2-1,1) = field%efield%rhoEF(irec2-1, 1)
78
             IF (vacuum%nvac == 2) THEN
79
                vxy(ncsh:vacuum%nmzxy,irec2-1,2) = field%efield%rhoEF(irec2-1, 2)
80 81
             END IF
          ELSE
82
             vxy(ncsh:vacuum%nmzxy,irec2-1,1:vacuum%nvac) = 0.0
83 84
          END IF

85
          vcons = 2.0*vcons/SINH(g*2.0*(field%efield%zsigma+cell%z1))
86 87 88 89 90 91

          ivac_loop1: DO ivac = 1,vacuum%nvac
             z = cell%z1
             imz_loop1: DO imz = 1,ncsh-1
                ! As "z" > 0 in this subroutine, the integrand is the same
                ! for both ivac -- but the integral bounds are reversed
92 93
                e_m = SINH(g*(field%efield%zsigma+cell%z1-z))
                e_p = SINH(g*(z+field%efield%zsigma+cell%z1))
94 95 96 97
                fra(ncsh-imz) = REAL(rhtxy(imz,irec2-1,ivac))* e_m
                fia(ncsh-imz) = AIMAG(rhtxy(imz,irec2-1,ivac))*e_m
                frb(imz) = REAL(rhtxy(imz,irec2-1,ivac))* e_p
                fib(imz) = AIMAG(rhtxy(imz,irec2-1,ivac))*e_p
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
                z = z + vacuum%delz
             END DO imz_loop1
             CALL intgz1(fra,vacuum%delz,ncsh-1,alpha(1,ivac,1),tail=.FALSE.)
             CALL intgz1(fia,vacuum%delz,ncsh-1,alpha(1,ivac,2),tail=.FALSE.)
             CALL qsf(vacuum%delz,frb,beta(1,ivac,1),ncsh-1,1)
             CALL qsf(vacuum%delz,fib,beta(1,ivac,2),ncsh-1,1)
          END DO ivac_loop1

          IF (ivac == 2) THEN
             ! Honour reversed integral bounds
             alpha(:,ivac,1) = -alpha(:,ivac,1)
             alpha(:,ivac,2) = -alpha(:,ivac,2)
             beta(:,ivac,1) = -beta(:,ivac,1)
             beta(:,ivac,2) = -beta(:,ivac,2)
          END IF

          alph1 = CMPLX(alpha(ncsh-1,1,1),alpha(ncsh-1,1,2))
          IF (vacuum%nvac == 1) THEN
             IF (sym%invs) THEN
                alph2 = CONJG(alph1)
             ELSE
                alph2 = alph1
             END IF
          ELSE
             alph2 = CMPLX(alpha(ncsh-1,2,1),alpha(ncsh-1,2,2))
          END IF
          ivac_loop2: DO ivac = 1,vacuum%nvac
             z = cell%z1
             IF (ivac == 1) alph0 = alph2
             IF (ivac == 2) alph0 = alph1
             imz_loop2: DO imz = 1,ncsh-1
                betaz = CMPLX(beta(imz,ivac,1),beta(imz,ivac,2))
                alphaz = CMPLX(alpha(ncsh-imz,ivac,1),&
                     &                        alpha(ncsh-imz,ivac,2))
132 133
                e_m = SINH(g*(field%efield%zsigma+cell%z1-z))
                e_p = SINH(g*(z+field%efield%zsigma+cell%z1))
134 135
                test = e_m*(alph0+betaz) + e_p*alphaz
                IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
136
                vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac)&
137
                     &                  +  vcons * test
138
                IF (ALLOCATED (field%efield%C1)) THEN
139 140 141
                   e_m = exp_save( -g*z  )
                   e_p = exp_save( g*z  )
                   IF (ivac == 1) THEN ! z > 0
142
                      vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac)&
143 144
                           &                                       + field%efield%C1(irec2-1)*e_p&
                           &                                       + field%efield%C2(irec2-1)*e_m
145
                   ELSE ! z < 0
146
                      vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac)&
147 148
                           &                                       + field%efield%C1(irec2-1)*e_m&
                           &                                       + field%efield%C2(irec2-1)*e_p
149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
                   END IF
                END IF
                z = z + vacuum%delz
             END DO imz_loop2
          END DO ivac_loop2
          alphm(irec2-1,1) = alph1
          alphm(irec2-1,2) = alph2

          ! ********** NEUMANN ************************************
       ELSE
          ivac_loop3: DO ivac = 1,vacuum%nvac
             z = cell%z1
             imz_loop3: DO imz = 1,vacuum%nmzxy
                e_m = exp_save( -g*z )
                e_p = exp_save( g*z )
164 165 166 167
                fra(ip-imz) = REAL(rhtxy(imz,irec2-1,ivac))* e_m
                fia(ip-imz) = AIMAG(rhtxy(imz,irec2-1,ivac))*e_m
                frb(imz) = REAL(rhtxy(imz,irec2-1,ivac))* e_p
                fib(imz) = AIMAG(rhtxy(imz,irec2-1,ivac))*e_p
168 169 170 171
                z = z + vacuum%delz
             END DO imz_loop3

             ! Add external field, if segmented
172 173
             IF (ALLOCATED (field%efield%rhoEF)) THEN
                z = cell%z1 + field%efield%zsigma
174 175 176 177 178 179 180
                e_m = exp_save( -g*z )
                e_p = exp_save( g*z )

                ! The equation has a minus sign as "rhtxy" contains the electron density
                ! (a positive number representing a negative charge) while rhoEF
                ! specifies the charges in terms of the (positive) elementary charge "e".
                fra(ip-ncsh) = fra(ip-ncsh)&
181
                     &                        - REAL (field%efield%rhoEF(irec2-1, ivac))*e_m
182
                fia(ip-ncsh) = fia(ip-ncsh)&
183
                     &                        - AIMAG (field%efield%rhoEF(irec2-1, ivac))*e_m
184
                frb(ncsh)    = frb(ncsh)&
185
                     &                       - REAL (field%efield%rhoEF(irec2-1, ivac))*e_p
186
                fib(ncsh)    = fib(ncsh)&
187
                     &                       - AIMAG (field%efield%rhoEF(irec2-1, ivac))*e_p
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 213 214 215
             END IF
             CALL intgz1(fra,vacuum%delz,vacuum%nmzxy,alpha(1,ivac,1),tail=.TRUE.)
             CALL intgz1(fia,vacuum%delz,vacuum%nmzxy,alpha(1,ivac,2),tail=.TRUE.)
             CALL qsf(vacuum%delz,frb,beta(1,ivac,1),vacuum%nmzxy,1)
             CALL qsf(vacuum%delz,fib,beta(1,ivac,2),vacuum%nmzxy,1)
          END DO ivac_loop3

          alph1 = CMPLX(alpha(vacuum%nmzxy,1,1),alpha(vacuum%nmzxy,1,2))
          IF (vacuum%nvac.EQ.1) THEN
             IF (sym%invs) THEN
                alph2 = CONJG(alph1)
             ELSE
                alph2 = alph1
             END IF
          ELSE
             alph2 = CMPLX(alpha(vacuum%nmzxy,2,1),alpha(vacuum%nmzxy,2,2))
          END IF
          ivac_loop4: DO ivac = 1,vacuum%nvac
             z = cell%z1
             IF (ivac.EQ.1) alph0 = alph2
             IF (ivac.EQ.2) alph0 = alph1
             imz_loop4: DO imz = 1,vacuum%nmzxy
                betaz = CMPLX(beta(imz,ivac,1),beta(imz,ivac,2))
                alphaz = CMPLX(alpha(ip-imz,ivac,1),alpha(ip-imz,ivac,2))
                e_m = exp_save( -g*z  )
                e_p = exp_save( g*z  )
                test = e_m*(alph0+betaz) + e_p*alphaz
                IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
216
                vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) +&
217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
                     &                    vcons * test
                z = z + vacuum%delz
             END DO imz_loop4
          END DO ivac_loop4
          alphm(irec2-1,1) = alph1
          alphm(irec2-1,2) = alph2
       END IF ! Neumann (vs. Dirichlet)
    END DO irec2_loop
  END SUBROUTINE vvacxy

  !------------------------------------------------------------------
  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_vvacxy