Line data Source code
1 : MODULE m_vintcz
2 : ! *************************************************************
3 : ! z-dependent part of Coulomb potential in the interstitial *
4 : ! [-d/2,d/2] region c.l.fu, r.podloucky *
5 : ! *************************************************************
6 : ! modified for thick films to avoid underflows gb`06
7 : !---------------------------------------------------------------
8 : CONTAINS
9 648576 : COMPLEX FUNCTION vintcz(stars,vacuum,cell,input,field,z,nrec2,psq,vnew,rhobar,sig1dh,vz1dh,alphm,vslope,sigma_disc,l_dfptvgen,l_corr,diff_vmz1dh,sigma_disc2)
10 : USE m_constants
11 : USE m_types
12 :
13 : IMPLICIT NONE
14 :
15 : TYPE(t_stars), INTENT(IN) :: stars
16 : TYPE(t_vacuum), INTENT(IN) :: vacuum
17 : TYPE(t_cell), INTENT(IN) :: cell
18 : TYPE(t_input), INTENT(IN) :: input
19 : TYPE(t_field), INTENT(IN) :: field
20 :
21 : INTEGER, INTENT(IN) :: nrec2
22 : COMPLEX, INTENT(IN) :: rhobar,vslope
23 : COMPLEX, INTENT(IN) :: sig1dh,vz1dh
24 : REAL, INTENT(IN) :: z
25 :
26 : COMPLEX, INTENT(IN) :: psq(stars%ng3),vnew(:,:,:)!,vxy(:,:,:) !(vacuum%nmzxyd,stars%ng2-1,2)
27 : COMPLEX, INTENT(IN) :: alphm(stars%ng2,2)
28 : complex, intent(in) :: sigma_disc(2)
29 : LOGICAL, INTENT(IN) :: l_dfptvgen, l_corr
30 : COMPLEX, INTENT(IN) :: diff_vmz1dh
31 :
32 : complex, optional, intent(in) :: sigma_disc2(2)
33 : !REAL, INTENT(IN) :: vz(:,:) !(vacuum%nmzd,2,jspins)
34 :
35 : COMPLEX :: argr,sumrr,vcons1,test,c_ph,phas
36 : REAL :: bj0,dh,fit,g,g3,q,qdh,vcons2,zf
37 : REAL :: e_m,e_p,cos_q,sin_q
38 : INTEGER :: ig3n,im,iq,ivac,k1,k2,nrec2r
39 : COMPLEX :: newdp, newdm, newdp2, newdm2
40 :
41 648576 : dh = cell%z1
42 648576 : sumrr = (0.,0.)
43 648576 : vintcz = (0.,0.)
44 :
45 648576 : newdp = CMPLX(0.0,0.0)
46 648576 : newdm = CMPLX(0.0,0.0)
47 648576 : newdp2 = CMPLX(0.0,0.0)
48 648576 : newdm2 = CMPLX(0.0,0.0)
49 :
50 : !---> if z is in the vacuum, use vacuum representations (m.w.)
51 648576 : IF (ABS(z).GE.cell%z1) THEN
52 224339 : ivac = 1
53 224339 : IF (z.LT.0.0) THEN
54 110774 : ivac = 2
55 224339 : IF (vacuum%nvac==1) ivac = 1
56 : END IF
57 224339 : zf = (ABS(z)-cell%z1)/vacuum%delz + 1.0
58 224339 : im = zf
59 224339 : q = zf - im
60 : ! For q/=0 in DFPT, the is no G+q=0, so all stars are treated in the G/=0 way.
61 897356 : IF (nrec2.EQ.1.AND.((.NOT.l_dfptvgen).OR.norm2(stars%center)<1e-8)) THEN
62 : fit = 0.5* (q-1.)* (q-2.)*REAL(vnew(im,1,ivac)) -&
63 : & q* (q-2.)*REAL(vnew(im+1,1,ivac)) +&
64 616 : & 0.5*q* (q-1.)*REAL(vnew(im+2,1,ivac))
65 616 : vintcz = CMPLX(fit,0.0)
66 223723 : ELSE IF (im+2.LE.vacuum%nmzxy) THEN
67 223723 : if (z<0) THEN
68 110475 : call stars%map_2nd_vac(vacuum,nrec2,nrec2r,phas) ! TODO: AN TB; will this work?
69 : else
70 113248 : nrec2r=nrec2
71 113248 : phas=cmplx(1.0,0.0)
72 : end if
73 : vintcz = phas*0.5* (q-1.)* (q-2.)*vnew(im,nrec2r,ivac) -&
74 : q* (q-2.)*vnew(im+1,nrec2r,ivac) +&
75 223723 : 0.5*q* (q-1.)*vnew(im+2,nrec2r,ivac)
76 : END IF
77 224339 : RETURN
78 : END IF
79 :
80 : ! For q/=0 in DFPT, there is no G+q=0, so all stars are treated in the G/=0 way.
81 1696948 : IF (nrec2==1.AND.((.NOT.l_dfptvgen).OR.norm2(stars%center)<1e-8)) THEN ! ----> g=0 coefficient
82 :
83 : ! New discontinuity correction
84 1406 : IF (l_corr) THEN
85 0 : newdp = - psq(1) * dh
86 0 : newdp2 = -psq(1) * dh**2
87 : END IF
88 :
89 58760 : DO iq = -stars%mx3,stars%mx3
90 57354 : IF (iq.EQ.0) CYCLE
91 55948 : ig3n = stars%ig(0,0,iq)
92 : ! ----> use only stars within the g_max sphere (oct.97 shz)
93 57354 : IF (ig3n.NE.0) THEN
94 53136 : q = iq*cell%bmat(3,3)
95 53136 : qdh = q*dh
96 53136 : bj0 = SIN(qdh)/qdh
97 53136 : argr = ImagUnit*q*z
98 53136 : sumrr = (EXP(argr)-EXP(ImagUnit*qdh))/ (q*q) + ImagUnit*COS(qdh)* (dh-z)/q + bj0* (z*z-dh*dh)/2.
99 53136 : vintcz = vintcz + fpi_const*psq(ig3n)*sumrr
100 :
101 : ! New discontinuity correction
102 53136 : IF (l_corr) THEN
103 0 : newdp = newdp + ImagUnit * psq(ig3n) * cmplx(cos(qdh), sin(qdh)) * dh / qdh
104 0 : newdp2 = newdp2 + psq(ig3n) * cmplx(cos(qdh), sin(qdh)) / q**2
105 : END IF
106 : END IF
107 : END DO
108 : ! -----> v2(z)
109 : vintcz = vintcz + vz1dh - fpi_const* (dh-z)*&
110 1406 : & (sig1dh-rhobar/2.* (dh-z))
111 : ! Discontinuity correction
112 1406 : vintcz = vintcz - fpi_const * sigma_disc(1)*(dh-z)
113 : !IF (PRESENT(sigma_disc2)) vintcz = vintcz + fpi_const * sigma_disc2(1)
114 :
115 : ! New discontinuity correction
116 1406 : IF (l_dfptvgen.AND.l_corr) vintcz = vintcz - fpi_const * newdp * (dh-z) + fpi_const * newdp2
117 :
118 : ! Correct the interstitial potential by the mismatch we calculated as
119 : ! a boundary condition for DFPT q=0.
120 : !IF (l_dfptvgen) vintcz = vintcz - diff_vmz1dh*z/(2*dh) + diff_vmz1dh/2
121 :
122 1406 : IF (field%efield%dirichlet .AND. vslope /= 0.0) THEN
123 0 : vintcz = vintcz + vslope * (dh-z)
124 : END IF
125 : ! ----> (g.ne.0) coefficients
126 : ELSE
127 422831 : k1 = stars%kv2(1,nrec2)
128 422831 : k2 = stars%kv2(2,nrec2)
129 :
130 422831 : newdp = cmplx(0.0,0.0)
131 422831 : newdm = cmplx(0.0,0.0)
132 422831 : newdp2 = cmplx(0.0,0.0)
133 422831 : newdm2 = cmplx(0.0,0.0)
134 :
135 12877554 : DO iq = -stars%mx3,stars%mx3
136 12454723 : ig3n = stars%ig(k1,k2,iq)
137 : ! ----> use only stars within the g_max sphere (oct.97 shz)
138 12877554 : IF (ig3n.NE.0) THEN
139 7888019 : c_ph = stars%rgphs(k1,k2,iq)
140 : ! -----> v3(z)
141 7888019 : q = iq*cell%bmat(3,3)
142 7888019 : g = stars%sk2(nrec2)
143 7888019 : g3 = stars%sk3(ig3n)
144 7888019 : vcons1 = fpi_const*psq(ig3n)*c_ph / (g3*g3)
145 7888019 : IF (field%efield%dirichlet) THEN
146 0 : e_m = 0.0
147 0 : e_p = 0.0
148 0 : vcons1 = vcons1/(g*SINH(g*2*(dh+field%efield%zsigma)))
149 0 : e_m = e_m - EXP(-ImagUnit*q*dh)*(g*COSH(g*(-field%efield%zsigma))+ImagUnit*q*SINH(g*(-field%efield%zsigma)))
150 0 : e_p = e_p + EXP(ImagUnit*q*dh)*(-g*COSH(g*(-field%efield%zsigma))+ImagUnit*q*SINH(g*(-field%efield%zsigma)))
151 0 : sumrr = EXP(ImagUnit*q*z)
152 0 : e_m = e_m + sumrr*(g*COSH(g*(z+dh+field%efield%zsigma))-ImagUnit*q*SINH(g*(z+dh+field%efield%zsigma)))
153 0 : e_p = e_p + sumrr*(g*COSH(g*(z-dh-field%efield%zsigma))-ImagUnit*q*SINH(g*(z-dh-field%efield%zsigma)))
154 : vintcz = vintcz+ vcons1*(e_m*SINH(g*(field%efield%zsigma+dh-z))&
155 0 : +e_p*SINH(g*(field%efield%zsigma+dh+z)))
156 : ELSE
157 7888019 : sumrr = (0.0,0.0)
158 7888019 : vcons2 = - 1.0 / (2.*g)
159 7888019 : e_m = exp_safe( -g*(z+dh) ) !exp_safe handles overflow
160 7888019 : e_p = exp_safe( g*(z-dh) )
161 : ! --> sum over gz-stars
162 7888019 : cos_q = COS(q*dh)
163 7888019 : sin_q = SIN(q*dh)
164 : sumrr = sumrr + CMPLX(COS(q*z),SIN(q*z)) + vcons2 *&
165 : & ( (g + ImagUnit*q) * e_p * (cos_q + ImagUnit*sin_q) +&
166 7888019 : & (g - ImagUnit*q) * e_m * (cos_q - ImagUnit*sin_q) )
167 7888019 : vintcz = vintcz + vcons1*sumrr
168 :
169 : ! New discontinuity correction
170 7888019 : IF (iq==0.AND.l_corr) newdp = newdp - psq(ig3n) * dh
171 : IF (iq/=0.AND.l_corr) newdp = newdp + ImagUnit * psq(ig3n) * cmplx(cos_q, sin_q) * dh / (q*dh)
172 : IF (iq==0.AND.l_corr) newdm = newdm - psq(ig3n) * dh
173 422831 : IF (iq/=0.AND.l_corr) newdm = newdm - ImagUnit * psq(ig3n) * cmplx(cos_q,-sin_q) * dh / (q*dh)
174 7465188 : IF (iq==0.AND.l_corr) newdp2 = newdp2 - psq(ig3n) * dh**2
175 422831 : IF (iq/=0.AND.l_corr) newdp2 = newdp2 + psq(ig3n) * cmplx(cos_q, sin_q) / q**2
176 7465188 : IF (iq==0.AND.l_corr) newdm2 = newdm2 + psq(ig3n) * dh**2
177 12454723 : IF (iq/=0.AND.l_corr) newdm2 = newdm2 - psq(ig3n) * cmplx(cos_q,-sin_q) / q**2
178 :
179 : END IF ! Neumann (vs. Dirichlet)
180 : END IF ! ig3d /= 0
181 : END DO
182 : ! ----> v4(z)
183 422831 : IF (field%efield%dirichlet) THEN
184 0 : e_m = SINH(g*(field%efield%zsigma+dh - z))
185 0 : e_p = SINH(g*(field%efield%zsigma+dh + z))
186 0 : test = e_m*alphm(nrec2,2) + e_p*alphm(nrec2,1)
187 0 : test = fpi_const/(g*SINH(g*2*(field%efield%zsigma+dh))) * test
188 0 : IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
189 0 : vintcz = vintcz + test
190 0 : IF (ALLOCATED (field%efield%C1)) THEN
191 0 : g = stars%sk2(nrec2)
192 0 : e_m = exp_safe (-g*z)
193 0 : e_p = exp_safe ( g*z)
194 : vintcz = vintcz + field%efield%C1(nrec2-1)*e_m &
195 0 : & + field%efield%C2(nrec2-1)*e_p
196 : END IF
197 : ELSE ! Neumann
198 422831 : e_m = exp_safe( -g*z )
199 422831 : e_p = exp_safe( g*z )
200 422831 : test = e_m*alphm(nrec2,2) + e_p*alphm(nrec2,1)
201 422831 : IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
202 422831 : vintcz = vintcz + tpi_const/g* test
203 :
204 : ! New discontinuity correction
205 : !IF (l_dfptvgen.AND.l_corr.AND.nrec2==1) THEN
206 : ! e_m = exp_safe( - g * abs(z-cell%z1) )
207 : ! e_p = exp_safe( - g * abs(z+cell%z1) )
208 : ! test = e_m * newdp + e_p * newdm
209 : ! if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
210 : ! vintcz = vintcz + tpi_const/g * test
211 : ! !IF (abs(z-cell%z1)<1e-9) e_m = 0.0
212 : ! !IF (abs(z+cell%z1)<1e-9) e_p = 0.0
213 : ! test = e_m * newdp2 * sign(1.0,z-cell%z1)+ e_p * newdm2 * sign(1.0,z+cell%z1)
214 : ! if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
215 : ! vintcz = vintcz - tpi_const * test
216 : !END IF
217 :
218 : END IF
219 : END IF
220 : END FUNCTION vintcz
221 :
222 16621700 : PURE REAL FUNCTION exp_safe(x)
223 : ! replace exp by a function that does not under/overflow dw09
224 : IMPLICIT NONE
225 :
226 : REAL, INTENT(IN) :: x
227 : REAL, PARAMETER :: maxexp = LOG(2.0)*MAXEXPONENT(2.0)
228 : REAL, PARAMETER :: minexp = LOG(2.0)*MINEXPONENT(2.0)
229 :
230 16621700 : IF ( ABS(x)>minexp .AND. ABS(x)<maxexp ) THEN
231 16621700 : exp_safe = EXP(x)
232 : ELSE
233 0 : IF ( x > 0 ) THEN
234 0 : IF ( x > minexp ) THEN
235 : exp_safe = EXP(maxexp)
236 : ELSE
237 0 : exp_safe = EXP(minexp)
238 : END IF
239 : ELSE
240 0 : IF ( -x > minexp ) THEN
241 : exp_safe = EXP(-maxexp)
242 : ELSE
243 0 : exp_safe = EXP(-minexp)
244 : END IF
245 : END IF
246 : END IF
247 16621700 : END FUNCTION exp_safe
248 : END MODULE m_vintcz
|
