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 70950 : COMPLEX FUNCTION vintcz(&
10 : & stars,vacuum,cell,sym,input,field,&
11 70950 : & z,nrec2,psq,vxy,vz,rhobar,sig1dh,vz1dh,alphm)
12 :
13 : USE m_constants
14 : USE m_types
15 : IMPLICIT NONE
16 : ! ..
17 : ! .. Scalar Arguments ..
18 : TYPE(t_stars),INTENT(IN) :: stars
19 : TYPE(t_vacuum),INTENT(IN) :: vacuum
20 : TYPE(t_cell),INTENT(IN) :: cell
21 : TYPE(t_sym),INTENT(IN) :: sym
22 : TYPE(t_input),INTENT(IN) :: input
23 : TYPE(t_field),INTENT(IN) :: field
24 : INTEGER, INTENT (IN) :: nrec2
25 : COMPLEX, INTENT (IN) :: rhobar
26 : REAL, INTENT (IN) :: sig1dh,vz1dh ,z
27 : ! ..
28 : ! .. Array Arguments ..
29 : COMPLEX, INTENT (IN) :: psq(stars%ng3),vxy(:,:,:) !(vacuum%nmzxyd,stars%ng2-1,2)
30 : COMPLEX, INTENT (IN) :: alphm(stars%ng2,2)
31 : REAL, INTENT (IN) :: vz(:,:) !(vacuum%nmzd,2,jspins)
32 :
33 : COMPLEX argr,sumrr,vcons1,test,c_ph
34 : REAL bj0,dh,fit,g,g3,q,qdh,signz,vcons2,zf
35 : REAL e_m,e_p,cos_q,sin_q
36 : INTEGER ig3n,im,iq,ivac,k1,k2,m0,nrz,nz
37 : ! ..
38 : ! .. Intrinsic Functions ..
39 : INTRINSIC abs,cmplx,conjg,cos,exp,sin
40 : ! ..
41 : !
42 :
43 70950 : dh = cell%z1
44 70950 : sumrr = (0.,0.)
45 70950 : vintcz = (0.,0.)
46 : !---> if z is in the vacuum, use vacuum representations (m.w.)
47 70950 : IF (ABS(z).GE.cell%z1) THEN
48 22945 : ivac = 1
49 22945 : IF (z.LT.0.0) THEN
50 10910 : ivac = 2
51 10910 : IF (sym%invs .OR. sym%zrfs) ivac = 1
52 : END IF
53 22945 : zf = (ABS(z)-cell%z1)/vacuum%delz + 1.0
54 22945 : im = zf
55 22945 : q = zf - im
56 22945 : IF (nrec2.EQ.1) THEN
57 : fit = 0.5* (q-1.)* (q-2.)*vz(im,ivac) -&
58 : & q* (q-2.)*vz(im+1,ivac) +&
59 97 : & 0.5*q* (q-1.)*vz(im+2,ivac)
60 97 : vintcz = CMPLX(fit,0.0)
61 22848 : ELSE IF (im+2.LE.vacuum%nmzxy) THEN
62 : vintcz = 0.5* (q-1.)* (q-2.)*vxy(im,nrec2-1,ivac) -&
63 : & q* (q-2.)*vxy(im+1,nrec2-1,ivac) +&
64 22848 : & 0.5*q* (q-1.)*vxy(im+2,nrec2-1,ivac)
65 22848 : IF ((sym%invs.AND. (.NOT.sym%zrfs)) .AND.&
66 7824 : & z.LT.0) vintcz = CONJG(vintcz)
67 : END IF
68 : RETURN
69 : END IF
70 : !
71 48005 : IF (nrec2.EQ.1) THEN
72 293 : m0 = -stars%mx3
73 293 : IF (sym%zrfs .OR. sym%invs) m0 = 0
74 : ! ----> g=0 coefficient
75 : ! -----> v1(z)
76 19187 : DO iq = m0,stars%mx3
77 9447 : IF (iq.EQ.0) CYCLE
78 9154 : ig3n = stars%ig(0,0,iq)
79 : ! ----> use only stars within the g_max sphere (oct.97 shz)
80 9447 : IF (ig3n.NE.0) THEN
81 8861 : q = iq*cell%bmat(3,3)
82 8861 : nz = stars%nstr(ig3n)
83 8861 : sumrr = (0.0,0.0)
84 : ! --> sum over gz-stars
85 26583 : DO nrz = 1,nz
86 17722 : signz = 3. - 2.*nrz
87 17722 : q = signz*q
88 17722 : qdh = q*dh
89 17722 : bj0 = SIN(qdh)/qdh
90 17722 : argr = ImagUnit*q*z
91 : sumrr = sumrr + (EXP(argr)-EXP(ImagUnit*qdh))/ (q*q) +&
92 26583 : & ImagUnit*COS(qdh)* (dh-z)/q + bj0* (z*z-dh*dh)/2.
93 : ENDDO
94 8861 : vintcz = vintcz + fpi_const*psq(ig3n)*sumrr
95 : ENDIF
96 : ENDDO
97 : ! -----> v2(z)
98 : vintcz = vintcz + vz1dh - fpi_const* (dh-z)*&
99 293 : & (sig1dh-rhobar/2.* (dh-z))
100 293 : IF (field%efield%dirichlet .AND. field%efield%vslope /= 0.0) THEN
101 0 : vintcz = vintcz + field%efield%vslope * (dh-z)
102 : END IF
103 : ! ----> (g.ne.0) coefficients
104 : ELSE
105 47712 : m0 = -stars%mx3
106 47712 : IF (sym%zrfs) m0 = 0
107 47712 : k1 = stars%kv2(1,nrec2)
108 47712 : k2 = stars%kv2(2,nrec2)
109 2284320 : DO iq = m0,stars%mx3
110 2236608 : ig3n = stars%ig(k1,k2,iq)
111 : ! ----> use only stars within the g_max sphere (oct.97 shz)
112 2284320 : IF (ig3n.NE.0) THEN
113 1469736 : c_ph = stars%rgphs(k1,k2,iq)
114 : ! -----> v3(z)
115 1469736 : q = iq*cell%bmat(3,3)
116 1469736 : g = stars%sk2(nrec2)
117 1469736 : g3 = stars%sk3(ig3n)
118 1469736 : vcons1 = fpi_const*psq(ig3n)*c_ph / (g3*g3)
119 1469736 : nz = 1
120 1469736 : IF (sym%zrfs) nz = stars%nstr(ig3n)/stars%nstr2(nrec2)
121 1469736 : IF (field%efield%dirichlet) THEN
122 0 : e_m = 0.0
123 0 : e_p = 0.0
124 0 : vcons1 = vcons1/(g*SINH(g*2*(dh+field%efield%zsigma)))
125 0 : loop_vacua: DO nrz = 1,nz
126 0 : signz = 3. - 2.*nrz
127 0 : q = signz*q
128 : e_m = e_m - EXP(-ImagUnit*q*dh)*(g*COSH(g*(-field%efield%zsigma))&
129 0 : & +ImagUnit*q*SINH(g*(-field%efield%zsigma)))
130 : e_p = e_p + EXP(ImagUnit*q*dh)*(-g*COSH(g*(-field%efield%zsigma))&
131 0 : & +ImagUnit*q*SINH(g*(-field%efield%zsigma)))
132 0 : sumrr = EXP(ImagUnit*q*z)
133 : e_m = e_m + sumrr*(g*COSH(g*(z+dh+field%efield%zsigma))&
134 0 : & -ImagUnit*q*SINH(g*(z+dh+field%efield%zsigma)))
135 : e_p = e_p + sumrr*(g*COSH(g*(z-dh-field%efield%zsigma))&
136 0 : & -ImagUnit*q*SINH(g*(z-dh-field%efield%zsigma)))
137 : END DO loop_vacua
138 : vintcz = vintcz&
139 : & + vcons1*(e_m*SINH(g*(field%efield%zsigma+dh-z))&
140 0 : & +e_p*SINH(g*(field%efield%zsigma+dh+z)))
141 : ELSE
142 1469736 : sumrr = (0.0,0.0)
143 1469736 : vcons2 = - 1.0 / (2.*g)
144 1469736 : e_m = exp_SAVE( -g*(z+dh) ) !exp_save handles overflow
145 1469736 : e_p = exp_SAVE( g*(z-dh) )
146 : ! --> sum over gz-stars
147 3055632 : vacua: DO nrz = 1,nz
148 1585896 : signz = 3. - 2.*nrz
149 1585896 : q = signz*q
150 1585896 : cos_q = COS(q*dh)
151 1585896 : sin_q = SIN(q*dh)
152 : sumrr = sumrr + CMPLX(COS(q*z),SIN(q*z)) + vcons2 *&
153 : & ( (g + ImagUnit*q) * e_p * (cos_q + ImagUnit*sin_q) +&
154 3055632 : & (g - ImagUnit*q) * e_m * (cos_q - ImagUnit*sin_q) )
155 : END DO vacua
156 1469736 : vintcz = vintcz + vcons1*sumrr
157 : END IF ! Neumann (vs. Dirichlet)
158 : END IF ! ig3d /= 0
159 : ENDDO
160 : ! ----> v4(z)
161 47712 : IF (field%efield%dirichlet) THEN
162 0 : e_m = SINH(g*(field%efield%zsigma+dh - z))
163 0 : e_p = SINH(g*(field%efield%zsigma+dh + z))
164 0 : test = e_m*alphm(nrec2-1,2) + e_p*alphm(nrec2-1,1)
165 0 : test = fpi_const/(g*SINH(g*2*(field%efield%zsigma+dh))) * test
166 0 : IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
167 0 : vintcz = vintcz + test
168 0 : IF (ALLOCATED (field%efield%C1)) THEN
169 0 : g = stars%sk2(nrec2)
170 0 : e_m = exp_save (-g*z)
171 0 : e_p = exp_save ( g*z)
172 : vintcz = vintcz + field%efield%C1(nrec2-1)*e_m&
173 0 : & + field%efield%C2(nrec2-1)*e_p
174 : END IF
175 : ELSE ! Neumann
176 47712 : e_m = exp_save( -g*z )
177 47712 : e_p = exp_save( g*z )
178 47712 : test = e_m*alphm(nrec2-1,2) + e_p*alphm(nrec2-1,1)
179 47712 : IF ( 2.0 * test == test ) test = CMPLX(0.0,0.0)
180 47712 : vintcz = vintcz + tpi_const/g* test
181 : END IF
182 :
183 : ENDIF
184 :
185 : END FUNCTION vintcz
186 :
187 3034896 : PURE REAL FUNCTION exp_save(x)
188 : ! replace exp by a function that does not under/overflow dw09
189 : IMPLICIT NONE
190 : REAL ,INTENT(IN) :: x
191 : REAL, PARAMETER :: maxexp = LOG(2.0)*MAXEXPONENT(2.0)
192 : REAL, PARAMETER :: minexp = LOG(2.0)*MINEXPONENT(2.0)
193 :
194 3034896 : IF ( ABS(x)>minexp .AND. ABS(x)<maxexp ) THEN
195 3034896 : exp_SAVE = EXP(x)
196 : ELSE
197 0 : IF ( x > 0 ) THEN
198 0 : IF ( x > minexp ) THEN
199 : exp_save = EXP(maxexp)
200 : ELSE
201 0 : exp_save = EXP(minexp)
202 : ENDIF
203 : ELSE
204 0 : IF ( -x > minexp ) THEN
205 : exp_save = EXP(-maxexp)
206 : ELSE
207 0 : exp_save = EXP(-minexp)
208 : ENDIF
209 : ENDIF
210 : ENDIF
211 3034896 : END FUNCTION exp_save
212 : END MODULE m_vintcz
|
