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 37 : subroutine vvacxy( stars, vacuum, cell, sym, input, field, rhtxy, vxy, alphm, l_dfptvgen )
14 : use m_types
15 : use m_constants
16 : use m_intgr, only: intgz1
17 : use m_qsf
18 :
19 : implicit none
20 :
21 : type(t_input), intent(in) :: input
22 : type(t_field), intent(in) :: field
23 : type(t_vacuum), intent(in) :: vacuum
24 : type(t_sym), intent(in) :: sym
25 : type(t_stars), intent(in) :: stars
26 : type(t_cell), intent(in) :: cell
27 :
28 : complex, intent(in) :: rhtxy(vacuum%nmzxyd,stars%ng2,2)
29 : complex, intent(inout) :: vxy(vacuum%nmzxyd,stars%ng2,2)
30 : complex, intent(out) :: alphm(stars%ng2,2)
31 : logical, intent(in) :: l_dfptvgen
32 :
33 : complex :: alph0, alphaz, betaz, test, phas
34 : real :: g, vcons, z, e_m, e_p,arg
35 : integer :: imz, ip, irec2, ivac, ncsh,irec2r, start_star
36 37 : real :: fra(vacuum%nmzxyd), frb(vacuum%nmzxyd), fia(vacuum%nmzxyd), fib(vacuum%nmzxyd)
37 37 : real :: alpha(vacuum%nmzxyd,2,2,stars%ng2), beta(vacuum%nmzxyd,2,2,stars%ng2)
38 :
39 37 : if ( allocated( field%efield%rhoef ) .or. field%efield%dirichlet ) then
40 0 : ncsh = field%efield%zsigma / vacuum%delz + 1.01
41 : ! if nmzxy < ncsh, the inhomogenous field cannot be represented.
42 : ! if nmzxy > ncsh, the boundary condition is wrong - and the
43 : ! potential is very wavy - try it yourself, if you don't believe.
44 0 : if ( vacuum%nmzxyd < ncsh ) then
45 0 : write (oUnit,*) 'error vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '< ncsh(zsigma) = ', ncsh
46 0 : call judft_error( "error: vacuum%nmzxyd < ncsh", calledby="vvacxy" )
47 0 : else if ( vacuum%nmzxyd > ncsh ) then
48 0 : write (oUnit,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
49 0 : write (*,*) 'warning vvacxy.f: vacuum%nmzxyd =', vacuum%nmzxyd, '> ncsh(zsigma) = ', ncsh
50 0 : call judft_warn( "nmzxyd > ncsh", calledby="vvacxy" )
51 : end if
52 : end if
53 :
54 37 : start_star = 2
55 : ! For q/=0 in DFPT, there is no G+q=0, so all stars are treated in the G/=0 way.
56 37 : if (l_dfptvgen) then
57 0 : if (norm2(stars%center)>1e-8) start_star = 1
58 : end if
59 : ! 2-dim star loop g/=0
60 37 : ip = vacuum%nmzxy + 1
61 : ! irec2_loop: do irec2 = 2, stars%ng2
62 : ! g = stars%sk2(irec2)
63 : ! vcons = tpi_const / g
64 37 : if ( field%efield%dirichlet ) then ! Dirichlet
65 0 : do irec2 = start_star, stars%ng2
66 0 : g = stars%sk2(irec2)
67 0 : vcons = tpi_const / g
68 0 : if ( allocated( field%efield%rhoef ) ) then
69 0 : vxy(ncsh:vacuum%nmzxy,irec2,1) = field%efield%rhoef(irec2-1, 1)
70 0 : if ( vacuum%nvac == 2 ) then
71 0 : vxy(ncsh:vacuum%nmzxy,irec2,2) = field%efield%rhoef(irec2-1, 2)
72 : end if
73 : else
74 0 : vxy(ncsh:vacuum%nmzxy,irec2,1:vacuum%nvac) = 0.0
75 : end if
76 0 : vcons = 2.0 * vcons / sinh( g * 2.0 * ( field%efield%zsigma + cell%z1 ) )
77 0 : do ivac = 1, vacuum%nvac
78 0 : z = cell%z1
79 0 : do imz = 1, ncsh-1
80 : ! as "z" > 0 in this subroutine, the integrand is the same
81 : ! for both ivac -- but the integral bounds are reversed
82 0 : e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
83 0 : e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
84 0 : fra(ncsh-imz) = real( rhtxy(imz,irec2,ivac) ) * e_m
85 0 : fia(ncsh-imz) = aimag( rhtxy(imz,irec2,ivac) ) * e_m
86 0 : frb(imz) = real( rhtxy(imz,irec2,ivac) ) * e_p
87 0 : fib(imz) = aimag( rhtxy(imz,irec2,ivac) ) * e_p
88 0 : z = z + vacuum%delz
89 : end do
90 0 : call intgz1( fra, vacuum%delz, ncsh - 1, alpha(1,ivac,1,irec2), .false. )
91 0 : call intgz1( fia, vacuum%delz, ncsh - 1, alpha(1,ivac,2,irec2), .false. )
92 0 : call qsf( vacuum%delz, frb, beta(1,ivac,1,irec2), ncsh - 1, 1 )
93 0 : call qsf( vacuum%delz, fib, beta(1,ivac,2,irec2), ncsh - 1, 1 )
94 : end do
95 :
96 0 : if ( ivac == 2 ) then
97 : ! honour reversed integral bounds
98 0 : alpha(:,ivac,:,irec2) = - alpha(:,ivac,:,irec2)
99 0 : beta(:,ivac,:,irec2) = - beta(:,ivac,:,irec2)
100 : end if
101 : end do
102 0 : do irec2 = start_star, stars%ng2
103 0 : g = stars%sk2(irec2)
104 0 : vcons = tpi_const / g
105 0 : alphm(irec2,1) = cmplx( alpha(vacuum%nmzxy,1,1,irec2), alpha(vacuum%nmzxy,1,2,irec2) )
106 0 : if ( vacuum%nvac == 1 ) then
107 0 : call stars%map_2nd_vac(vacuum,irec2,irec2r,phas)
108 0 : alphm(irec2r,2) = phas * cmplx(alpha(vacuum%nmzxy,1,1,irec2),alpha(vacuum%nmzxy,1,2,irec2))
109 : else
110 0 : alphm(irec2,2) = cmplx( alpha(vacuum%nmzxy,2,1,irec2), alpha(vacuum%nmzxy,2,2,irec2) )
111 : end if
112 0 : do ivac = 1, vacuum%nvac
113 0 : z = cell%z1
114 0 : if ( ivac == 1 ) alph0 = alphm(irec2,2)
115 0 : if ( ivac == 2 ) alph0 = alphm(irec2,1)
116 0 : do imz = 1, ncsh-1
117 0 : betaz = cmplx( beta(imz,ivac,1,irec2), beta(imz,ivac,2,irec2) )
118 0 : alphaz = cmplx( alpha(ncsh-imz,ivac,1,irec2), alpha(ncsh-imz,ivac,2,irec2) )
119 0 : e_m = sinh( g * ( field%efield%zsigma + cell%z1 - z ) )
120 0 : e_p = sinh( g * ( field%efield%zsigma + cell%z1 + z ) )
121 0 : test = e_m * ( alph0 + betaz ) + e_p * alphaz
122 0 : if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
123 0 : vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + vcons * test
124 0 : if ( allocated( field%efield%c1 ) ) then
125 0 : e_m = exp_safe( - g * z )
126 0 : e_p = exp_safe( g * z )
127 0 : if ( ivac == 1 ) then ! z > 0
128 0 : vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + field%efield%c1(irec2-1) * e_p + field%efield%c2(irec2-1) * e_m
129 : else ! z < 0
130 0 : vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + field%efield%c1(irec2-1) * e_m + field%efield%c2(irec2-1) * e_p
131 : end if
132 : end if
133 0 : z = z + vacuum%delz
134 : end do
135 : end do
136 : end do
137 : else ! Neumann
138 13450 : do irec2 = start_star, stars%ng2
139 13413 : g = stars%sk2(irec2)
140 13413 : vcons = tpi_const / g
141 :
142 28513 : do ivac = 1, vacuum%nvac
143 15063 : z = cell%z1
144 1521363 : do imz = 1, vacuum%nmzxy
145 3012600 : e_m = exp_safe( - g * z )
146 3012600 : e_p = exp_safe( g * z )
147 1506300 : fra(ip-imz) = real( rhtxy(imz,irec2,ivac) ) * e_m
148 1506300 : fia(ip-imz) = aimag( rhtxy(imz,irec2,ivac) ) * e_m
149 1506300 : frb(imz) = real( rhtxy(imz,irec2,ivac) ) * e_p
150 1506300 : fib(imz) = aimag( rhtxy(imz,irec2,ivac) ) * e_p
151 1521363 : z = z + vacuum%delz
152 : end do
153 : ! add external field, if segmented
154 15063 : if ( allocated( field%efield%rhoef ) ) then
155 0 : z = cell%z1 + field%efield%zsigma
156 0 : e_m = exp_safe( - g * z )
157 0 : e_p = exp_safe( 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 15063 : call intgz1( fra, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,1,irec2), .true. )
167 15063 : call intgz1( fia, vacuum%delz, vacuum%nmzxy, alpha(1,ivac,2,irec2), .true. )
168 15063 : call qsf( vacuum%delz, frb, beta(1,ivac,1,irec2), vacuum%nmzxy, 1 )
169 28476 : call qsf( vacuum%delz, fib, beta(1,ivac,2,irec2), vacuum%nmzxy, 1 )
170 : end do
171 : end do
172 :
173 13450 : do irec2 = start_star, stars%ng2
174 13413 : alphm(irec2,1) = cmplx( alpha(vacuum%nmzxy,1,1,irec2), alpha(vacuum%nmzxy,1,2,irec2) )
175 13450 : if ( vacuum%nvac == 1 ) then
176 11763 : call stars%map_2nd_vac(vacuum,irec2,irec2r,phas)
177 11763 : alphm(irec2r,2) = phas * cmplx(alpha(vacuum%nmzxy,1,1,irec2),alpha(vacuum%nmzxy,1,2,irec2))
178 : else
179 1650 : alphm(irec2,2) = cmplx( alpha(vacuum%nmzxy,2,1,irec2), alpha(vacuum%nmzxy,2,2,irec2) )
180 : end if
181 : end do
182 :
183 13450 : do irec2 = start_star, stars%ng2
184 13413 : g = stars%sk2(irec2)
185 13413 : vcons = tpi_const / g
186 :
187 28513 : do ivac = 1, vacuum%nvac
188 15063 : z = cell%z1
189 15063 : if ( ivac == 1 ) alph0 = alphm(irec2,2)
190 15063 : if ( ivac == 2 ) alph0 = alphm(irec2,1)
191 1534776 : do imz = 1, vacuum%nmzxy
192 1506300 : betaz = cmplx( beta(imz,ivac,1,irec2), beta(imz,ivac,2,irec2) )
193 1506300 : alphaz = cmplx( alpha(ip-imz,ivac,1,irec2), alpha(ip-imz,ivac,2,irec2) )
194 3012600 : e_m = exp_safe( - g * z )
195 3012600 : e_p = exp_safe( g * z )
196 1506300 : test = e_m * ( alph0 + betaz ) + e_p * alphaz
197 1506300 : if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
198 1506300 : vxy(imz,irec2,ivac) = vxy(imz,irec2,ivac) + vcons * test
199 1521363 : z = z + vacuum%delz
200 : end do
201 : end do
202 : end do
203 : end if
204 37 : end subroutine vvacxy
205 :
206 6025200 : pure real function exp_safe( x )
207 : ! replace exp by a function that does not under/overflow
208 : implicit none
209 :
210 : real, intent(in) :: x
211 : real, parameter :: maxexp = log( 2.0 ) * maxexponent( 2.0 )
212 : real, parameter :: minexp = log( 2.0 ) * minexponent( 2.0 )
213 :
214 6025200 : if ( abs( x ) > minexp .and. abs( x ) < maxexp ) then
215 6025200 : exp_safe = exp( x )
216 : else
217 0 : if ( x > 0 ) then
218 : if ( x > minexp ) then
219 : exp_safe = exp( maxexp )
220 : else
221 : exp_safe = exp( minexp )
222 : end if
223 : else
224 : if ( - x > minexp ) then
225 : exp_safe = exp( - maxexp )
226 : else
227 : exp_safe = exp( - minexp )
228 : end if
229 : end if
230 : end if
231 0 : end function exp_safe
232 : end module m_vvacxy
|
