Line data Source code
1 : module m_vvacis
2 : ! **********************************************************
3 : ! g.ne.0 coefficients of vacuum coulomb potential *
4 : ! due to the interstitial charge density inside slab *
5 : ! c.l.fu, r.podloucky *
6 : ! **********************************************************
7 : ! modified for thick films to avoid underflows gb`06
8 : !---------------------------------------------------------------
9 : contains
10 44 : subroutine vvacis( stars, vacuum, cell, psq, input, field, vxy, l_dfptvgen, l_corr )
11 : use m_constants
12 : use m_types
13 :
14 : implicit none
15 :
16 : type(t_input), intent(in) :: input
17 : type(t_field), intent(in) :: field
18 : type(t_vacuum), intent(in) :: vacuum
19 : type(t_stars), intent(in) :: stars
20 : type(t_cell), intent(in) :: cell
21 :
22 : complex, intent(in) :: psq(stars%ng3)
23 : complex, intent(inout) :: vxy(vacuum%nmzxyd,stars%ng2,2)
24 : logical, intent(in) :: l_dfptvgen, l_corr
25 :
26 : complex :: arg, c_ph, sumr(2)
27 : real :: g, qz, signu, vcons, z, e_m, arg_r, arg_i, e_p
28 : integer :: ig3n, imz, ivac, k1, k2, kz, nrec2, start_star
29 : complex :: newdp, newdm, newdp2, newdm2, test
30 :
31 44 : newdp = cmplx(0.0,0.0)
32 44 : newdm = cmplx(0.0,0.0)
33 44 : newdp2 = cmplx(0.0,0.0)
34 44 : newdm2 =cmplx(0.0,0.0)
35 :
36 44 : start_star = 2
37 : ! For q/=0 in DFPT, there is no G+q=0, so all stars are treated in the G/=0 way.
38 44 : if (l_dfptvgen) then
39 0 : if (norm2(stars%center)>1e-8) start_star = 1
40 : end if
41 :
42 3501398 : vxy(:,start_star:,:) = cmplx( 0., 0. )
43 :
44 17377 : do nrec2 = start_star, stars%ng2
45 17333 : k1 = stars%kv2(1,nrec2)
46 17333 : k2 = stars%kv2(2,nrec2)
47 17333 : g = stars%sk2(nrec2)
48 17333 : if ( field%efield%dirichlet ) then
49 0 : vcons = 2.0 * tpi_const / ( g * sinh( g * 2.0 * ( field%efield%zsigma + cell%z1 ) ) )
50 0 : arg_r = g * ( cell%z1 + field%efield%zsigma + cell%z1 )
51 : else ! neumann
52 17333 : vcons = tpi_const / g
53 17333 : arg_r = exp_safe( - 2 * cell%z1 * g )
54 : end if
55 36360 : do ivac = 1, vacuum%nvac
56 18983 : sumr(ivac) = ( 0.0, 0.0 )
57 :
58 : IF (ivac==1) newdp = cmplx(0.0,0.0)
59 : IF (ivac==2) newdm = cmplx(0.0,0.0)
60 : IF (ivac==1) newdp2 = cmplx(0.0,0.0)
61 : IF (ivac==2) newdm2 = cmplx(0.0,0.0)
62 :
63 18983 : signu = 3. - 2. * ivac
64 528402 : do kz = -stars%mx3, stars%mx3
65 509419 : ig3n = stars%ig(k1,k2,kz)
66 509419 : c_ph = stars%rgphs(k1,k2,kz)
67 : ! use only stars within the g_max sphere
68 528402 : if ( ig3n /= 0 ) then
69 319979 : qz = kz * cell%bmat(3,3)
70 : ! sum over gz-stars
71 319979 : if ( field%efield%dirichlet ) then
72 : ! prefactor
73 0 : arg = exp( - ImagUnit * qz * cell%z1 ) / ( 2 * ( g ** 2 + qz ** 2 ) ) * psq(ig3n)
74 0 : if ( ivac == 1 ) then
75 : sumr(ivac) = sumr(ivac) + c_ph * exp( - arg_r ) * arg * ( & ! c_ph not tested in this case
76 : ( - exp( 2 * g * ( field%efield%zsigma + cell%z1 ) ) + exp( 2 * ( ImagUnit * qz * cell%z1 + arg_r ) ) ) * ( g - ImagUnit * qz ) &
77 0 : + ( - exp( 2 * g * cell%z1 ) + exp( 2 * ImagUnit * qz * cell%z1 ) ) * ( g + ImagUnit * qz ) )
78 : else
79 : sumr(ivac) = sumr(ivac) + c_ph * arg * ( &
80 : exp( arg_r ) * ( g + ImagUnit * qz ) &
81 : + exp( - arg_r ) * ( g - ImagUnit * qz ) &
82 : + 2 * exp( 2 * ImagUnit * qz * cell%z1 ) &
83 : * ( - g * cosh( g * ( - field%efield%zsigma ) ) &
84 0 : + ImagUnit * qz * sinh( - g * field%efield%zsigma ) ) )
85 : end if
86 : else
87 319979 : arg = g + signu * ImagUnit * qz
88 319979 : arg_i = signu * qz * cell%z1
89 319979 : sumr(ivac) = sumr(ivac) + c_ph * psq(ig3n) * ( cos( arg_i ) * ( 1 - arg_r ) + ImagUnit * sin( arg_i ) * ( 1 + arg_r ) ) / arg
90 :
91 : ! New discontinuity correction
92 : IF (l_corr) THEN
93 : IF (kz==0.AND.ivac==1) newdp = newdp - psq(ig3n) * cell%z1
94 : IF (kz/=0.AND.ivac==1) newdp = newdp + ImagUnit * psq(ig3n) * cmplx(cos(qz * cell%z1), sin(qz * cell%z1)) / qz
95 : IF (kz==0.AND.ivac==2) newdm = newdm - psq(ig3n) * cell%z1
96 : IF (kz/=0.AND.ivac==2) newdm = newdm - ImagUnit * psq(ig3n) * cmplx(cos(qz * cell%z1), sin(qz * cell%z1)) / qz
97 : IF (kz==0.AND.ivac==1) newdp2 = newdp2 - psq(ig3n) * cell%z1**2
98 : IF (kz/=0.AND.ivac==1) newdp2 = newdp2 + psq(ig3n) * cmplx(qz * cell%z1, qz * cell%z1) / qz**2
99 : IF (kz==0.AND.ivac==2) newdm2 = newdm2 + psq(ig3n) * cell%z1**2
100 : IF (kz/=0.AND.ivac==2) newdm2 = newdm2 - psq(ig3n) * cmplx(qz * cell%z1,-qz * cell%z1) / qz**2
101 : END IF
102 : end if
103 : end if
104 : end do
105 18983 : z = 0 ! moved cell%z1 into above equations
106 1934616 : do imz = 1, vacuum%nmzxy
107 1898300 : if ( field%efield%dirichlet ) then
108 0 : e_m = sinh( g * ( field%efield%zsigma - z ) )
109 : else ! neumann
110 1898300 : e_m = exp_safe( - g * z )
111 : end if
112 1898300 : vxy(imz,nrec2,ivac) = vxy(imz,nrec2,ivac) + vcons * sumr(ivac) * e_m
113 1917283 : z = z + vacuum%delz
114 : end do
115 : end do
116 : ! New discontinuity correction
117 : !z = 0
118 : !do ivac = 1, vacuum%nvac
119 : ! do imz = 1, vacuum%nmzxy
120 : ! IF (l_dfptvgen.AND.l_corr.AND.nrec2==1) THEN
121 : ! e_m = exp_safe( - g * abs(z-cell%z1) )
122 : ! e_p = exp_safe( - g * abs(z+cell%z1) )
123 : ! test = e_m * newdp + e_p * newdm
124 : ! if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
125 : ! vxy(imz,nrec2,ivac) = vxy(imz,nrec2,ivac) + tpi_const/g * test
126 : ! IF (abs(z-cell%z1)<1e-9) e_m = 0.0
127 : ! IF (abs(z+cell%z1)<1e-9) e_p = 0.0
128 : ! test = e_m * newdp2 * sign(1.0,z-cell%z1)+ e_p * newdm2 * sign(1.0,z+cell%z1)
129 : ! if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
130 : ! vxy(imz,nrec2,ivac) = vxy(imz,nrec2,ivac) - tpi_const * test
131 : ! END IF
132 : ! z = z + vacuum%delz
133 : ! end do
134 : !end do
135 : end do
136 44 : end subroutine vvacis
137 :
138 1915633 : pure real function exp_safe(x)
139 : ! replace exp by a function that does not under/overflow
140 :
141 : implicit none
142 :
143 : real, intent(in) :: x
144 : real, parameter :: maxexp = log( 2.0 ) * maxexponent( 2.0 )
145 : real, parameter :: minexp = log( 2.0 ) * minexponent( 2.0 )
146 :
147 1915633 : if ( abs(x) > minexp .and. abs(x) < maxexp ) then
148 1915633 : exp_safe = exp( x )
149 : else
150 0 : if ( x > 0 ) then
151 0 : if ( x > minexp ) then
152 : exp_safe = exp( maxexp )
153 : else
154 0 : exp_safe = exp( minexp )
155 : end if
156 : else
157 0 : if ( -x > minexp ) then
158 : exp_safe = exp( -maxexp )
159 : else
160 0 : exp_safe = exp( -minexp )
161 : end if
162 : end if
163 : end if
164 1915633 : end function exp_safe
165 : end module m_vvacis
|
