Line data Source code
1 : module m_vvacxy
2 : ! **********************************************************
3 : ! g/=0 coefficient of vacuum coulomb potential *
4 : ! due to warped vacuum charged density *
5 : ! c.l.fu, r.podloucky *
6 : ! **********************************************************
7 : ! modified for thick films to avoid underflows gb`06
8 : !---------------------------------------------------------------
9 :
10 : use m_judft
11 :
12 : contains
13 :
14 5 : subroutine vvacxy( stars, vacuum, cell, sym, input, field, rhtxy, vxy, alphm )
15 :
16 : use m_intgr, only: intgz1
17 : use m_constants
18 : use m_types
19 : use m_qsf
20 : implicit none
21 :
22 : type(t_input), intent(in) :: input
23 : type(t_field), intent(in) :: field
24 : type(t_vacuum), intent(in) :: vacuum
25 : type(t_sym), intent(in) :: sym
26 : type(t_stars), intent(in) :: stars
27 : type(t_cell), intent(in) :: cell
28 : complex, intent(in) :: rhtxy(vacuum%nmzxyd,stars%ng2-1,2)
29 : complex, intent(inout) :: vxy(vacuum%nmzxyd,stars%ng2-1,2)
30 : complex, intent(out) :: alphm(stars%ng2,2)
31 :
32 : complex :: alph0, alph2, alph1, alphaz, betaz, test
33 : real :: g, vcons, z, e_m, e_p
34 : integer :: imz, ip, irec2, ivac, ncsh
35 : logical :: tail
36 10 : real :: fra(vacuum%nmzxyd), frb(vacuum%nmzxyd), fia(vacuum%nmzxyd), fib(vacuum%nmzxyd)
37 10 : real :: alpha(vacuum%nmzxyd,2,2), beta(vacuum%nmzxyd,2,2)
38 : real, allocatable :: sig_top(:), sig_bot(:)
39 : intrinsic aimag, cmplx, conjg, exp, real
40 :
41 5 : if ( allocated( field%efield%rhoef ) .or. field%efield%dirichlet ) then
42 0 : ncsh = field%efield%zsigma / vacuum%delz + 1.01
43 : ! if nmzxy < ncsh, the inhomogenous field cannot be represented.
44 : ! if nmzxy > ncsh, the boundary condition is wrong - and the
45 : ! potential is very wavy - try it yourself, if you don't believe.
46 0 : if ( vacuum%nmzxyd < ncsh ) then
47 0 : write (6,*) 'error vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '< ncsh(zsigma) = ', ncsh
48 0 : call judft_error( "error: vacuum%nmzxyd < ncsh", calledby="vvacxy" )
49 0 : else if ( vacuum%nmzxyd > ncsh ) then
50 0 : write (6,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
51 0 : write (0,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
52 0 : call judft_warn( "nmzxyd > ncsh", calledby="vvacxy" )
53 : end if
54 : end if
55 :
56 : ! 2-dim star loop g/=0
57 5 : ip = vacuum%nmzxy + 1
58 1125 : irec2_loop: do irec2 = 2, stars%ng2
59 1120 : g = stars%sk2(irec2)
60 1120 : vcons = tpi_const / g
61 : ! ********** dirichlet ************************************
62 1125 : if ( field%efield%dirichlet ) then
63 0 : if ( allocated( field%efield%rhoef ) ) then
64 0 : vxy(ncsh:vacuum%nmzxy,irec2-1,1) = field%efield%rhoef(irec2-1, 1)
65 0 : if ( vacuum%nvac == 2 ) then
66 0 : vxy(ncsh:vacuum%nmzxy,irec2-1,2) = field%efield%rhoef(irec2-1, 2)
67 : end if
68 : else
69 0 : vxy(ncsh:vacuum%nmzxy,irec2-1,1:vacuum%nvac) = 0.0
70 : end if
71 0 : vcons = 2.0 * vcons / sinh( g * 2.0 * ( field%efield%zsigma + cell%z1 ) )
72 0 : ivac_loop1: do ivac = 1, vacuum%nvac
73 0 : z = cell%z1
74 0 : imz_loop1: do imz = 1, ncsh-1
75 : ! as "z" > 0 in this subroutine, the integrand is the same
76 : ! for both ivac -- but the integral bounds are reversed
77 0 : e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
78 0 : e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
79 0 : fra(ncsh-imz) = real( rhtxy(imz,irec2-1,ivac) ) * e_m
80 0 : fia(ncsh-imz) = aimag( rhtxy(imz,irec2-1,ivac) ) * e_m
81 0 : frb(imz) = real( rhtxy(imz,irec2-1,ivac) ) * e_p
82 0 : fib(imz) = aimag( rhtxy(imz,irec2-1,ivac) ) * e_p
83 0 : z = z + vacuum%delz
84 : end do imz_loop1
85 0 : call intgz1( fra, vacuum%delz, ncsh - 1, alpha(1,ivac,1), .false. )
86 0 : call intgz1( fia, vacuum%delz, ncsh - 1, alpha(1,ivac,2), .false. )
87 0 : call qsf( vacuum%delz, frb, beta(1,ivac,1), ncsh - 1, 1 )
88 0 : call qsf( vacuum%delz, fib, beta(1,ivac,2), ncsh - 1, 1 )
89 : end do ivac_loop1
90 :
91 0 : if ( ivac == 2 ) then
92 : ! honour reversed integral bounds
93 0 : alpha(:,ivac,1) = - alpha(:,ivac,1)
94 0 : alpha(:,ivac,2) = - alpha(:,ivac,2)
95 0 : beta(:,ivac,1) = - beta(:,ivac,1)
96 0 : beta(:,ivac,2) = - beta(:,ivac,2)
97 : end if
98 :
99 0 : alph1 = cmplx( alpha(ncsh-1,1,1), alpha(ncsh-1,1,2) )
100 0 : if ( vacuum%nvac == 1 ) then
101 0 : if ( sym%invs ) then
102 0 : alph2 = conjg( alph1 )
103 : else
104 : alph2 = alph1
105 : end if
106 : else
107 0 : alph2 = cmplx( alpha(ncsh-1,2,1), alpha(ncsh-1,2,2) )
108 : end if
109 0 : ivac_loop2: do ivac = 1, vacuum%nvac
110 0 : z = cell%z1
111 0 : if ( ivac == 1 ) alph0 = alph2
112 0 : if ( ivac == 2 ) alph0 = alph1
113 0 : imz_loop2: do imz = 1, ncsh-1
114 0 : betaz = cmplx( beta(imz,ivac,1), beta(imz,ivac,2) )
115 0 : alphaz = cmplx( alpha(ncsh-imz,ivac,1), alpha(ncsh-imz,ivac,2) )
116 0 : e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
117 0 : e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
118 0 : test = e_m * ( alph0 + betaz ) + e_p * alphaz
119 0 : if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
120 0 : vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) + vcons * test
121 0 : if ( allocated( field%efield%c1 ) ) then
122 0 : e_m = exp_save( - g * z )
123 0 : e_p = exp_save( g * z )
124 0 : if ( ivac == 1 ) then ! z > 0
125 : vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) &
126 : + field%efield%c1(irec2-1) * e_p &
127 0 : + field%efield%c2(irec2-1) * e_m
128 : else ! z < 0
129 : vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) &
130 : + field%efield%c1(irec2-1) * e_m &
131 0 : + field%efield%c2(irec2-1) * e_p
132 : end if
133 : end if
134 0 : z = z + vacuum%delz
135 : end do imz_loop2
136 : end do ivac_loop2
137 0 : alphm(irec2-1,1) = alph1
138 0 : alphm(irec2-1,2) = alph2
139 :
140 : ! ********** neumann ************************************
141 : else
142 2240 : ivac_loop3: do ivac = 1, vacuum%nvac
143 1120 : z = cell%z1
144 113120 : imz_loop3: do imz = 1, vacuum%nmzxy
145 112000 : e_m = exp_save( - g * z )
146 112000 : e_p = exp_save( g * z )
147 112000 : fra(ip-imz) = real( rhtxy(imz,irec2-1,ivac) ) * e_m
148 112000 : fia(ip-imz) = aimag( rhtxy(imz,irec2-1,ivac) ) * e_m
149 112000 : frb(imz) = real( rhtxy(imz,irec2-1,ivac) ) * e_p
150 112000 : fib(imz) = aimag( rhtxy(imz,irec2-1,ivac) ) * e_p
151 113120 : z = z + vacuum%delz
152 : end do imz_loop3
153 : ! add external field, if segmented
154 1120 : if ( allocated( field%efield%rhoef ) ) then
155 0 : z = cell%z1 + field%efield%zsigma
156 0 : e_m = exp_save( - g * z )
157 0 : e_p = exp_save( g * z )
158 : ! the equation has a minus sign as "rhtxy" contains the electron density
159 : ! (a positive number representing a negative charge) while rhoef
160 : ! specifies the charges in terms of the (positive) elementary charge "e".
161 0 : fra(ip-ncsh) = fra(ip-ncsh) - real( field%efield%rhoef(irec2-1,ivac) ) * e_m
162 0 : fia(ip-ncsh) = fia(ip-ncsh) - aimag( field%efield%rhoef(irec2-1,ivac) ) * e_m
163 0 : frb(ncsh) = frb(ncsh) - real( field%efield%rhoef(irec2-1,ivac) ) * e_p
164 0 : fib(ncsh) = fib(ncsh) - aimag( field%efield%rhoef(irec2-1,ivac) ) * e_p
165 : end if
166 1120 : call intgz1( fra, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,1), .true. )
167 1120 : call intgz1( fia, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,2), .true. )
168 1120 : call qsf( vacuum%delz, frb, beta(1,ivac,1), vacuum%nmzxy, 1 )
169 2240 : call qsf( vacuum%delz, fib, beta(1,ivac,2), vacuum%nmzxy, 1 )
170 : end do ivac_loop3
171 1120 : alph1 = cmplx( alpha(vacuum%nmzxy,1,1), alpha(vacuum%nmzxy,1,2) )
172 1120 : if ( vacuum%nvac == 1 ) then
173 1120 : if ( sym%invs ) then
174 1120 : alph2 = conjg( alph1 )
175 : else
176 : alph2 = alph1
177 : end if
178 : else
179 0 : alph2 = cmplx( alpha(vacuum%nmzxy,2,1), alpha(vacuum%nmzxy,2,2) )
180 : end if
181 2240 : ivac_loop4: do ivac = 1, vacuum%nvac
182 1120 : z = cell%z1
183 1120 : if ( ivac == 1 ) alph0 = alph2
184 1120 : if ( ivac == 2 ) alph0 = alph1
185 114240 : imz_loop4: do imz = 1, vacuum%nmzxy
186 112000 : betaz = cmplx( beta(imz,ivac,1), beta(imz,ivac,2) )
187 112000 : alphaz = cmplx( alpha(ip-imz,ivac,1), alpha(ip-imz,ivac,2) )
188 112000 : e_m = exp_save( - g * z )
189 112000 : e_p = exp_save( g * z )
190 112000 : test = e_m * ( alph0 + betaz ) + e_p * alphaz
191 112000 : if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
192 112000 : vxy(imz,irec2-1,ivac) = vxy(imz,irec2-1,ivac) + vcons * test
193 113120 : z = z + vacuum%delz
194 : end do imz_loop4
195 : end do ivac_loop4
196 1120 : alphm(irec2-1,1) = alph1
197 1120 : alphm(irec2-1,2) = alph2
198 : end if
199 : end do irec2_loop
200 5 : end subroutine vvacxy
201 :
202 : !------------------------------------------------------------------
203 448000 : pure real function exp_save( x )
204 : ! replace exp by a function that does not under/overflow
205 : implicit none
206 : real, intent(in) :: x
207 : real, parameter :: maxexp = log( 2.0 ) * maxexponent( 2.0 )
208 : real, parameter :: minexp = log( 2.0 ) * minexponent( 2.0 )
209 :
210 448000 : if ( abs( x ) > minexp .and. abs( x ) < maxexp ) then
211 448000 : exp_save = exp( x )
212 : else
213 0 : if ( x > 0 ) then
214 0 : if ( x > minexp ) then
215 : exp_save = exp( maxexp )
216 : else
217 0 : exp_save = exp( minexp )
218 : endif
219 : else
220 0 : if ( - x > minexp ) then
221 : exp_save = exp( - maxexp )
222 : else
223 0 : exp_save = exp( - minexp )
224 : endif
225 : endif
226 : endif
227 448000 : end function exp_save
228 :
229 : end module m_vvacxy
|
