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 4303119 : COMPLEX FUNCTION vintcz(stars,vacuum,cell,input,field,z,nrec2,psq,vnew,rhobar,sig1dh,vz1dh,alphm,vslope,l_dfptvgen)
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 4303119 : dh = cell%z1
42 4303119 : sumrr = (0.,0.)
43 4303119 : vintcz = (0.,0.)
44 :
45 : ! newdp = CMPLX(0.0,0.0)
46 : ! newdm = CMPLX(0.0,0.0)
47 : ! newdp2 = CMPLX(0.0,0.0)
48 : ! newdm2 = CMPLX(0.0,0.0)
49 :
50 : !---> if z is in the vacuum, use vacuum representations (m.w.)
51 4303119 : IF (ABS(z).GE.cell%z1) THEN
52 1967916 : ivac = 1
53 1967916 : IF (z.LT.0.0) THEN
54 972672 : ivac = 2
55 1967916 : IF (vacuum%nvac==1) ivac = 1
56 : END IF
57 1967916 : zf = (ABS(z)-cell%z1)/vacuum%delz + 1.0
58 1967916 : im = zf
59 1967916 : 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 7871664 : 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 8147 : & 0.5*q* (q-1.)*REAL(vnew(im+2,1,ivac))
65 8147 : vintcz = CMPLX(fit,0.0)
66 1959769 : ELSE IF (im+2.LE.vacuum%nmzxy) THEN
67 1959769 : if (z<0) THEN
68 968641 : call stars%map_2nd_vac(vacuum,nrec2,nrec2r,phas) ! TODO: AN TB; will this work?
69 : else
70 991128 : nrec2r=nrec2
71 991128 : 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 1959769 : 0.5*q* (q-1.)*vnew(im+2,nrec2r,ivac)
76 : END IF
77 1967916 : 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 9340812 : IF (nrec2==1.AND.((.NOT.l_dfptvgen).OR.norm2(stars%center)<1e-8)) THEN ! ----> g=0 coefficient
82 :
83 : ! ! New discontinuity correction
84 : ! IF (l_corr) THEN
85 : ! newdp = - psq(1) * dh
86 : ! newdp2 = -psq(1) * dh**2
87 : ! END IF
88 :
89 517698 : DO iq = -stars%mx3,stars%mx3
90 508076 : IF (iq.EQ.0) CYCLE
91 498454 : ig3n = stars%ig(0,0,iq)
92 : ! ----> use only stars within the g_max sphere (oct.97 shz)
93 508076 : IF (ig3n.NE.0) THEN
94 479210 : q = iq*cell%bmat(3,3)
95 479210 : qdh = q*dh
96 479210 : bj0 = SIN(qdh)/qdh
97 479210 : argr = ImagUnit*q*z
98 479210 : sumrr = (EXP(argr)-EXP(ImagUnit*qdh))/ (q*q) + ImagUnit*COS(qdh)* (dh-z)/q + bj0* (z*z-dh*dh)/2.
99 479210 : vintcz = vintcz + fpi_const*psq(ig3n)*sumrr
100 :
101 : ! ! New discontinuity correction
102 : ! IF (l_corr) THEN
103 : ! newdp = newdp + ImagUnit * psq(ig3n) * cmplx(cos(qdh), sin(qdh)) * dh / qdh
104 : ! 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 9622 : & (sig1dh-rhobar/2.* (dh-z))
111 :
112 : ! ! Discontinuity correction
113 : ! vintcz = vintcz - fpi_const * sigma_disc(1)*(dh-z)
114 : ! !IF (PRESENT(sigma_disc2)) vintcz = vintcz + fpi_const * sigma_disc2(1)
115 :
116 : ! ! New discontinuity correction
117 : ! IF (l_dfptvgen.AND.l_corr) vintcz = vintcz - fpi_const * newdp * (dh-z) + fpi_const * newdp2
118 :
119 : ! Correct the interstitial potential by the mismatch we calculated as
120 : ! a boundary condition for DFPT q=0.
121 : !IF (l_dfptvgen) vintcz = vintcz - diff_vmz1dh*z/(2*dh) + diff_vmz1dh/2
122 :
123 9622 : IF (field%efield%dirichlet .AND. vslope /= 0.0) THEN
124 0 : vintcz = vintcz + vslope * (dh-z)
125 : END IF
126 : ! ----> (g.ne.0) coefficients
127 : ELSE
128 2325581 : k1 = stars%kv2(1,nrec2)
129 2325581 : k2 = stars%kv2(2,nrec2)
130 :
131 : !newdp = cmplx(0.0,0.0)
132 : !newdm = cmplx(0.0,0.0)
133 : !newdp2 = cmplx(0.0,0.0)
134 : !newdm2 = cmplx(0.0,0.0)
135 :
136 118105238 : DO iq = -stars%mx3,stars%mx3
137 115779657 : ig3n = stars%ig(k1,k2,iq)
138 : ! ----> use only stars within the g_max sphere (oct.97 shz)
139 118105238 : IF (ig3n.NE.0) THEN
140 105476161 : c_ph = stars%rgphs(k1,k2,iq)
141 : ! -----> v3(z)
142 105476161 : q = iq*cell%bmat(3,3)
143 105476161 : g = stars%sk2(nrec2)
144 105476161 : g3 = stars%sk3(ig3n)
145 105476161 : vcons1 = fpi_const*psq(ig3n)*c_ph / (g3*g3)
146 105476161 : IF (field%efield%dirichlet) THEN
147 0 : e_m = 0.0
148 0 : e_p = 0.0
149 0 : vcons1 = vcons1/(g*SINH(g*2*(dh+field%efield%zsigma)))
150 0 : e_m = e_m - EXP(-ImagUnit*q*dh)*(g*COSH(g*(-field%efield%zsigma))+ImagUnit*q*SINH(g*(-field%efield%zsigma)))
151 0 : e_p = e_p + EXP(ImagUnit*q*dh)*(-g*COSH(g*(-field%efield%zsigma))+ImagUnit*q*SINH(g*(-field%efield%zsigma)))
152 0 : sumrr = EXP(ImagUnit*q*z)
153 0 : e_m = e_m + sumrr*(g*COSH(g*(z+dh+field%efield%zsigma))-ImagUnit*q*SINH(g*(z+dh+field%efield%zsigma)))
154 0 : e_p = e_p + sumrr*(g*COSH(g*(z-dh-field%efield%zsigma))-ImagUnit*q*SINH(g*(z-dh-field%efield%zsigma)))
155 : vintcz = vintcz+ vcons1*(e_m*SINH(g*(field%efield%zsigma+dh-z))&
156 0 : +e_p*SINH(g*(field%efield%zsigma+dh+z)))
157 : ELSE
158 105476161 : sumrr = (0.0,0.0)
159 105476161 : vcons2 = - 1.0 / (2.*g)
160 210952322 : e_m = exp_safe( -g*(z+dh) ) !exp_safe handles overflow
161 210952322 : e_p = exp_safe( g*(z-dh) )
162 : ! --> sum over gz-stars
163 105476161 : cos_q = COS(q*dh)
164 105476161 : sin_q = SIN(q*dh)
165 : sumrr = sumrr + CMPLX(COS(q*z),SIN(q*z)) + vcons2 *&
166 : & ( (g + ImagUnit*q) * e_p * (cos_q + ImagUnit*sin_q) +&
167 105476161 : & (g - ImagUnit*q) * e_m * (cos_q - ImagUnit*sin_q) )
168 105476161 : vintcz = vintcz + vcons1*sumrr
169 :
170 : ! ! New discontinuity correction
171 : ! IF (iq==0.AND.l_corr) newdp = newdp - psq(ig3n) * dh
172 : ! IF (iq/=0.AND.l_corr) newdp = newdp + ImagUnit * psq(ig3n) * cmplx(cos_q, sin_q) * dh / (q*dh)
173 : ! IF (iq==0.AND.l_corr) newdm = newdm - psq(ig3n) * dh
174 : ! IF (iq/=0.AND.l_corr) newdm = newdm - ImagUnit * psq(ig3n) * cmplx(cos_q,-sin_q) * dh / (q*dh)
175 : ! IF (iq==0.AND.l_corr) newdp2 = newdp2 - psq(ig3n) * dh**2
176 : ! IF (iq/=0.AND.l_corr) newdp2 = newdp2 + psq(ig3n) * cmplx(cos_q, sin_q) / q**2
177 : ! IF (iq==0.AND.l_corr) newdm2 = newdm2 + psq(ig3n) * dh**2
178 : ! IF (iq/=0.AND.l_corr) newdm2 = newdm2 - psq(ig3n) * cmplx(cos_q,-sin_q) / q**2
179 :
180 : END IF ! Neumann (vs. Dirichlet)
181 : END IF ! ig3d /= 0
182 : END DO
183 : ! ----> v4(z)
184 2325581 : IF (field%efield%dirichlet) THEN
185 0 : e_m = SINH(g*(field%efield%zsigma+dh - z))
186 0 : e_p = SINH(g*(field%efield%zsigma+dh + z))
187 0 : test = e_m*alphm(nrec2,2) + e_p*alphm(nrec2,1)
188 0 : test = fpi_const/(g*SINH(g*2*(field%efield%zsigma+dh))) * test
189 0 : IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
190 0 : vintcz = vintcz + test
191 0 : IF (ALLOCATED (field%efield%C1)) THEN
192 0 : g = stars%sk2(nrec2)
193 0 : e_m = exp_safe (-g*z)
194 0 : e_p = exp_safe ( g*z)
195 : vintcz = vintcz + field%efield%C1(nrec2-1)*e_m &
196 0 : & + field%efield%C2(nrec2-1)*e_p
197 : END IF
198 : ELSE ! Neumann
199 4651162 : e_m = exp_safe( -g*z )
200 4651162 : e_p = exp_safe( g*z )
201 2325581 : test = e_m*alphm(nrec2,2) + e_p*alphm(nrec2,1)
202 2325581 : IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
203 2325581 : vintcz = vintcz + tpi_const/g* test
204 :
205 : ! New discontinuity correction
206 : !IF (l_dfptvgen.AND.l_corr.AND.nrec2==1) THEN
207 : ! e_m = exp_safe( - g * abs(z-cell%z1) )
208 : ! e_p = exp_safe( - g * abs(z+cell%z1) )
209 : ! test = e_m * newdp + e_p * newdm
210 : ! if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
211 : ! vintcz = vintcz + tpi_const/g * test
212 : ! !IF (abs(z-cell%z1)<1e-9) e_m = 0.0
213 : ! !IF (abs(z+cell%z1)<1e-9) e_p = 0.0
214 : ! test = e_m * newdp2 * sign(1.0,z-cell%z1)+ e_p * newdm2 * sign(1.0,z+cell%z1)
215 : ! if ( 2.0 * test == test ) test = cmplx( 0.0, 0.0 )
216 : ! vintcz = vintcz - tpi_const * test
217 : !END IF
218 :
219 : END IF
220 : END IF
221 : END FUNCTION vintcz
222 :
223 215603484 : PURE REAL FUNCTION exp_safe(x)
224 : ! replace exp by a function that does not under/overflow dw09
225 : IMPLICIT NONE
226 :
227 : REAL, INTENT(IN) :: x
228 : REAL, PARAMETER :: maxexp = LOG(2.0)*MAXEXPONENT(2.0)
229 : REAL, PARAMETER :: minexp = LOG(2.0)*MINEXPONENT(2.0)
230 :
231 215603484 : IF ( ABS(x)>minexp .AND. ABS(x)<maxexp ) THEN
232 215603484 : exp_safe = EXP(x)
233 : ELSE
234 0 : IF ( x > 0 ) THEN
235 : IF ( x > minexp ) THEN
236 : exp_safe = EXP(maxexp)
237 : ELSE
238 : exp_safe = EXP(minexp)
239 : END IF
240 : ELSE
241 : IF ( -x > minexp ) THEN
242 : exp_safe = EXP(-maxexp)
243 : ELSE
244 : exp_safe = EXP(-minexp)
245 : END IF
246 : END IF
247 : END IF
248 0 : END FUNCTION exp_safe
249 : END MODULE m_vintcz
|
