Line data Source code
1 : MODULE m_efield
2 : USE m_juDFT
3 : USE m_constants
4 : IMPLICIT NONE
5 : PRIVATE
6 : PUBLIC :: e_field
7 : CONTAINS
8 80 : SUBROUTINE e_field(atoms, stars, sym, vacuum, cell, input,efield)
9 : !
10 : !*********************************************************************
11 : ! sets the values of the sheets of charge for external electric
12 : ! fields by requiring charge neutrality.
13 : ! the position of the sheets of charge relative to the vacuum
14 : ! boundary is fixed (10 a.u.), but can be set to a different
15 : ! value in the file apwefl.
16 : !
17 : ! modified and fixed 10-99
18 : !*********************************************************************
19 :
20 : USE m_types
21 : !USE m_setcor, ONLY: setcor
22 : IMPLICIT NONE
23 : ! ..
24 : ! .. Scalar Arguments ..
25 : TYPE(t_atoms), INTENT (IN) :: atoms
26 : Type(t_stars),INTENT(IN) :: stars
27 : TYPE(t_sym),INTENT(IN) :: sym
28 : TYPE(t_vacuum),INTENT(IN) :: vacuum
29 : TYPE(t_cell),INTENT(IN) :: cell
30 : TYPE(t_input),INTENT(IN) :: input
31 : TYPE(t_efield),INTENT(INOUT) :: efield
32 : ! ..
33 : ! ..
34 : ! .. Local Scalars ..
35 : REAL, parameter :: eps = 1.0e-7
36 : REAL qn,qe,bmu
37 : INTEGER n,iwd,nst,nc
38 : ! ..
39 : ! ..
40 : ! .. Local Parameters ..
41 : INTEGER, PARAMETER :: pTOP = 1, pBOT = 2, pTOPBOT = 3
42 : INTEGER, PARAMETER :: pADD = 1, pREPLACE = 2
43 : INTEGER, PARAMETER :: pALL = 1, pZERO = 2, pNONZERO = 3
44 : !
45 : !---> obtain total nuclear charge
46 : !
47 80 : qn=0.0
48 217 : DO n=1,atoms%ntype
49 217 : qn=qn+atoms%neq(n)*atoms%zatom(n)
50 : ENDDO
51 : !
52 : !---> obtain total electronic charge (in electron units)
53 : !
54 80 : qe=0.0
55 : !---> core electrons
56 217 : DO n = 1,atoms%ntype
57 217 : IF (atoms%zatom(n).GE.1.0) THEN
58 137 : qe=qe+atoms%neq(n)*atoms%econf(n)%core_electrons
59 137 : WRITE (oUnit,*) 'neq= ',atoms%neq(n),' ncore= ',qe
60 : ENDIF
61 : ENDDO
62 : !---> semi-core and valence electrons
63 80 : qe=qe+input%zelec
64 80 : WRITE (oUnit,"(A, F16.8)") 'zelec= ',input%zelec
65 :
66 80 : WRITE (oUnit, '(/,/,a)') ' parameters for external electric field:'
67 80 : WRITE (oUnit, '(3x,a,f12.5)') 'total electronic charge =', qe
68 80 : WRITE (oUnit, '(3x,a,f12.5)') 'total nuclear charge =', qn
69 :
70 80 : CALL read_efield (efield, stars%mx1, stars%mx2, vacuum%nvac, cell%area)
71 :
72 : ! Sign convention of electric field: E>0 repels electrons,
73 : ! consistent with conventional choice that current is
74 : ! opposite to electron flow. With this choice, qe < qn
75 : ! corresponds to E>0
76 : ! In case of Dirichlet boundary conditions, we ignore the
77 : ! value of "sigma" and take the surplus charge into account
78 : ! in vvac.
79 80 : if (efield%autocomp .or. efield%dirichlet) efield%sigma = 0.5*(qe-qn)
80 :
81 80 : CALL print_efield (oUnit,efield, cell%area, vacuum%nvac, cell%amat,vacuum%dvac)
82 :
83 : IF (.NOT. efield%dirichlet&
84 240 : & .AND. ABS (SUM (efield%sig_b) + 2*efield%sigma - (qe-qn)) > eps) THEN
85 0 : IF (ABS (SUM (efield%sig_b) - (qe-qn)) < eps) THEN
86 : CALL juDFT_error&
87 : & ("E field: top+bottom+film charge does not add up to "&
88 : & //"zero.",&
89 : & hint="Consider setting 'autocomp = false' in apwefl. "&
90 : & //"(By default, the number of electrons in the energy "&
91 : & //"window is automatically compensated via the charge "&
92 : & //"sheets.)",&
93 0 : & calledby ="efield")
94 : END IF
95 : CALL juDFT_error&
96 : & ("E field: top+bottom+film charge does not add up to zero"&
97 0 : & ,calledby ="efield")
98 : ENDIF
99 :
100 240 : IF (ABS (efield%sigma) > 0.49 .OR. ANY (ABS (efield%sig_b) > 0.49)) THEN
101 0 : WRITE (oUnit,*) 'If you really want to calculate an e-field this'
102 0 : WRITE (oUnit,*) 'big, uncomment STOP in efield.f !'
103 : CALL juDFT_error("E-field too big or No. of e- not correct"&
104 0 : & ,calledby ="efield")
105 : ENDIF
106 :
107 80 : IF (efield%l_segmented) THEN
108 : CALL V_seg_EF(&
109 : & efield,&
110 0 : & vacuum, stars)
111 :
112 0 : IF (efield%plot_rho)&
113 : & CALL print_rhoEF(&
114 : & efield, stars%mx1, stars%mx2, vacuum%nvac, stars%ng2, sym%nop, sym%nop2,&
115 : & stars%ng2, stars%kv2, sym%mrot, sym%symor, sym%tau, sym%invtab,&
116 0 : & cell%area, stars%nstr2, cell%amat)
117 : END IF
118 :
119 80 : IF (ALLOCATED (efield%sigEF)) DEALLOCATE (efield%sigEF)
120 :
121 80 : IF (efield%dirichlet .AND. ALLOCATED (efield%rhoEF))&
122 0 : & call set_dirchlet_coeff (efield, vacuum%dvac, stars%ng2, stars%sk2, vacuum%nvac)
123 :
124 : CONTAINS
125 :
126 0 : SUBROUTINE set_dirchlet_coeff (E, dvac, nq2, sk2, nvac)
127 : TYPE(t_efield), INTENT(INOUT), TARGET :: E
128 : REAL, INTENT(IN) :: dvac
129 : INTEGER, INTENT(IN) :: nq2, nvac
130 : REAL, INTENT(IN) :: sk2(:)
131 :
132 0 : COMPLEX, POINTER :: V(:,:)
133 : REAL :: x, y, g
134 : INTEGER :: ig, ivac1, ivac2
135 :
136 0 : ALLOCATE (E%C1(nq2-1), E%C2(nq2-1))
137 0 : E%l_dirichlet_coeff = .true.
138 :
139 0 : V => E%rhoEF
140 0 : ivac1 = 1
141 0 : ivac2 = nvac
142 :
143 0 : DO ig = 1, nq2-1
144 0 : g = sk2(ig+1)
145 0 : x = EXP ( g*(E%zsigma+Dvac/2.0))
146 0 : y = EXP (-g*(E%zsigma+Dvac/2.0))
147 0 : E%C1(ig) = (V(ig, ivac1)*x - V(ig, ivac2)*y) / (x**2-y**2)
148 0 : E%C2(ig) = (V(ig, ivac2)*x - V(ig, ivac1)*y) / (x**2-y**2)
149 : END DO
150 0 : END SUBROUTINE set_dirchlet_coeff
151 :
152 : ! Read the electric field parameters from the 'apwefl' file
153 : ! (except for sigma, which is determined by the charge in the film).
154 : !
155 : ! The position of the charged sheet can be given (in a.u.) via
156 : ! zsigma, an additional, in x-y homogeneous field can be given via
157 : ! sig_b(1) and sig_b(2). NOTE: The position zsigma counts from the
158 : ! interstitial-vacuum boundary ("z1") and not from the center (z=0).
159 : !
160 : ! There are two input formats. The old format is
161 : ! zsigma
162 : ! sig_b(nvac)
163 : ! where sig_b is optional. zsigma is 10.0 by default while b_sig is zero.
164 : !
165 : ! In the new format, lines starting with "!" or "#" are treated as comments
166 : ! and all keywords are case insensitive.
167 : !
168 : ! Namelist (optionally, defaults to the shown values):
169 : ! &efield zsigma=10.0, sig_b=0.0,0.0, plot_charge=f, plot_rho=f,
170 : ! autocomp=t, dirichlet=f, eV=f /
171 : !
172 : ! By default, Neumann boundary conditions are assumed and additional charge
173 : ! is placed as surface charge density at the charge sheets located at zsigma.
174 : ! If Dirichlet is true, a (segmented) metal plate at zsigma is assumed,
175 : ! which has a charge of zero - or the value given by sig_b or via adding
176 : ! different shapes. In case of Dirichlet boundary conditions, the value is
177 : ! regarded as potential in Hartree (or [electron] Volt, if "eV" is set to
178 : ! true.)
179 : !
180 : ! A fancy potential can now be created by placing charge in a rectangular
181 : ! or circular region via the "rect" and "circ" directive. Both directives
182 : ! are followed by "top, "bot" or "TopBot"/"BotTop" and then a point (x,y),
183 : ! which denotes for "rect" the bottom-left coordinate and for "circ" the
184 : ! centre. Next parameter is the width/height for "rect" and the diameter
185 : ! radius for "circ". Followed by a parameter for the charge (or the
186 : ! potential in case of Dirichlet boundary conditions). Optional
187 : ! arguments: The charge can be added ("add", default) to the previous
188 : ! charge in this sheet - or it replaces ("replace") the charge. The
189 : ! charge can be placed in the whole area ("all", default) or only where
190 : ! previously no charge ("zero") or a "nonzero" charge was.
191 : ! The sig_b charge is added after all the shapes and is thus treated like
192 : ! a "rect 0,0 top/bot 1.0,1.0 b_sig(1)/b_sig(2)" would be.
193 : !
194 : ! Note: All coordinates are relative coordinates. The regions can exceed the
195 : ! x-y extend of the film; e.g. using
196 : ! circ top 0,0 0.25 0.5
197 : ! one places 1/2 an electron in a quarter circle with origin (0,0).
198 : !
199 : ! WARNING: "circ" creates a perfect circle on the grid, however, it only
200 : ! matches a circle and not an ellipse if the k1d/k2d ratio matches the
201 : ! crystal a/b ratio.
202 : ! TODO: Support input in absolute values instead of relative ones.
203 : ! (Could be an option to the namelist or another (optional) tag.)
204 : !
205 : ! rectSinX: Create a sinusoidal potential in x direction (constant in y
206 : ! direction for any x value), i.e.
207 : ! V(x,y) = A * sin(2*pi*n*x + 2*pi*delta)
208 : ! A is the amplitude; however, the argument in apwefl is not A directly
209 : ! but A' = A * Lz, where Lz is the number of points in z direction.
210 : ! Contrary to "circ" and "rect", charges are mask out without being
211 : ! redistributed to non-masked positions. Note:
212 : ! Int_0^(2pi) A*|sin(x)| dx = 4*A)
213 : ! n is the order and delta the offset. "x" is relative to the rectangle
214 : ! and thus between 0 and 1. The syntax is:
215 : ! rectSinX top/bot x,y w,h, A', n, delta [options]
216 : !
217 : ! datafile top/bot filename [zero_avg] [option]
218 : ! Read data points from "filename"; if zero_avg is present, average the
219 : ! read data to zero; replace/add is supported, but zero/not_zero is not.
220 : !
221 : ! polygon: Create a polygon-shaped charge distribution; note that the
222 : ! currently used algorithm does not always give the perfactly shaped
223 : ! polygon - and the edge points are not always included in the polygon.
224 : ! Syntax:
225 : ! polygon top/bot n_points, x1,y1, ..., xn,yn charge [options]
226 : !
227 : !
228 : ! Example 1: To have two top plates (segments):
229 : ! rect top 0, 0 0.5,1.0 0.2
230 : ! rect top 0.5,0 0.5,1.0 -0.2
231 : !
232 : ! Example 2: To have a charged ring with 0.2e and -0.2e of charge
233 : ! evenly distributed outside this ring:
234 : ! circ topBot 0.5,0.50.2 1 ! Create temporary an inner ring
235 : ! circ topBot 0.5,0.50.3 0.2 zero ! Create outer ring
236 : ! circ topBot 0.5,0.50.2 0 replace ! Delete inner ring
237 : ! rect topBot 0,0 1,1 -0.2 zero ! Place smeared opposite charge
238 : !
239 : !
240 : ! If autocomp is .true., the extra charge in the film (Window) is distributed
241 : ! over both external sheets; if it is .false. one needs to compensate it
242 : ! manually.
243 : !
244 : ! Note: The namelist has to be placed before the shapes when Dirichlet
245 : ! boundary conditions are used as otherwise the potential is divided by
246 : ! the number of grid points...
247 :
248 80 : SUBROUTINE read_efield (E, k1d, k2d, nvac, area)
249 : TYPE(t_efield), INTENT(INOUT) :: E
250 : INTEGER, INTENT(IN) :: k1d, k2d, nvac
251 : REAL, INTENT(IN) :: area
252 :
253 : REAL :: tmp
254 : INTEGER :: i
255 :
256 : ! New format
257 400 : ALLOCATE(E%sigEF(3*k1d, 3*k2d, nvac))
258 207681 : E%sigEF = 0.0
259 80 : if (allocated(e%shapes)) then
260 80 : DO i=1,SIZE(e%shapes)
261 80 : CALL read_shape (E, e%shapes(i), nvac)
262 : END DO
263 : endif
264 80 : IF (e%l_eV) THEN
265 0 : E%sig_b(:) = E%sig_b/hartree_to_ev_const
266 0 : E%sigEF(:,:,:) = E%sigEF/hartree_to_ev_const
267 : END IF
268 : ! Save average sigEF potential in sig_b and remove it from
269 : ! sigEF to avoid double counting; i.e. make g_|| = 0 of sigEF == 0.
270 80 : IF (E%dirichlet) THEN
271 0 : tmp = sum(E%sigEF(:,:,1))/SIZE(E%sigEF)*nvac
272 0 : E%sig_b(1) = E%sig_b(1) + tmp
273 0 : E%sigEF(:,:,1) = E%sigEF(:,:,1) - tmp
274 : ELSE
275 111020 : tmp = sum(E%sigEF(:,:,1))
276 80 : E%sig_b(1) = E%sig_b(1) + tmp
277 111260 : E%sigEF(:,:,1) = E%sigEF(:,:,1) - tmp/SIZE(E%sigEF)*nvac
278 : END IF
279 80 : IF (nvac > 1) THEN
280 71 : IF (E%dirichlet) THEN
281 0 : tmp = sum(E%sigEF(:,:,2))/SIZE(E%sigEF)*nvac
282 0 : E%sig_b(2) = E%sig_b(2) + tmp
283 0 : E%sigEF(:,:,2) = E%sigEF(:,:,2) - tmp
284 : ELSE
285 96581 : tmp = sum(E%sigEF(:,:,2))
286 71 : E%sig_b(2) = E%sig_b(2) + tmp
287 96794 : E%sigEF(:,:,2) = E%sigEF(:,:,2) - tmp/SIZE(E%sigEF)*nvac
288 : END IF
289 9 : ELSE IF (E%sig_b(2) /= 0.0) THEN
290 : CALL juDFT_error&
291 : & ("z-mirror/inversion symmetry but second e-field"//&
292 0 : & "sheet specified",calledby ="efield")
293 : END IF
294 :
295 207681 : IF (ALL (E%sigEF == 0.0)) THEN
296 80 : DEALLOCATE (E%sigEF)
297 : ELSE
298 0 : E%l_segmented = .TRUE.
299 : END IF
300 80 : END SUBROUTINE read_efield
301 :
302 0 : SUBROUTINE read_shape (E, orig_str, nvac)
303 : USE m_constants, ONLY : pimach
304 : TYPE(t_efield), INTENT(INOUT) :: E
305 : CHARACTER(*), INTENT(IN) :: orig_str
306 : INTEGER, INTENT(IN) :: nvac
307 :
308 : REAL :: TWO_PI
309 : INTEGER :: ipos, action, iopt, ivac, ix, iy, nx, ny, cnt, i, j
310 : INTEGER :: np
311 : CHARACTER(10) :: tag, pos
312 : CHARACTER(200) :: str
313 : REAL :: x, y, h, w, radius, charge, order, shift, tmp
314 0 : REAL, ALLOCATABLE :: poly(:,:)
315 0 : REAL :: data(UBOUND(E%sigEF, dim=1), UBOUND(E%sigEF, dim=2))
316 0 : LOGICAL :: mask(UBOUND(E%sigEF, dim=1), UBOUND(E%sigEF, dim=2))
317 :
318 0 : TWO_PI = 2.0 * pimach()
319 0 : str = lower_case (orig_str)
320 :
321 0 : IF (str(1:8) == 'rectsinx') THEN
322 0 : READ (str, *) tag, pos, x, y, w, h, charge, order, shift
323 0 : IF (tag /= 'rectsinx')&
324 : & CALL juDFT_error("Internal error in read_shape (rectSinX)"&
325 0 : & ,calledby ="efield")
326 0 : ELSE IF (str(1:4) == 'circ') THEN
327 0 : READ (str, *) tag, pos, x, y, radius, charge
328 0 : IF (tag /= 'circ') CALL juDFT_error&
329 : & ("Internal error in read_shape (circ)",calledby ="efield"&
330 0 : & )
331 0 : ELSE IF (str(1:4) == 'rect') THEN
332 0 : READ (str, *) tag, pos, x, y, w, h, charge
333 0 : IF (tag /= 'rect') CALL juDFT_error&
334 : & ("Internal error in read_shape (rect)",calledby ="efield"&
335 0 : & )
336 0 : ELSE IF (str(1:7) == 'polygon') THEN
337 0 : READ (str, *) tag, pos, np
338 0 : ALLOCATE (poly(2,np))
339 0 : READ (str, *) tag, pos, np, poly, charge
340 0 : IF (tag /= 'polygon')&
341 : & CALL juDFT_error("Internal error in read_shape (polygon)"&
342 0 : & ,calledby ="efield")
343 0 : ELSE IF (str(1:8) == "datafile") THEN
344 0 : READ (str, *) tag, pos
345 0 : IF (tag /= 'datafile') CALL juDFT_error&
346 0 : & ("Internal error in read_shape",calledby ="efield")
347 : ELSE
348 0 : CALL juDFT_error("ERROR reading ",calledby ="efield")
349 : END IF
350 :
351 0 : IF (pos == 'top') THEN
352 0 : ipos = pTOP
353 0 : ELSEIF (pos == 'bot') THEN
354 0 : ipos = pBOT
355 0 : ELSEIF (pos == 'topbot' .OR. pos == 'bottop') THEN
356 0 : ipos = pTOPBOT
357 : ELSE
358 0 : CALL juDFT_error("ERROR reading ",calledby ="efield")
359 : END IF
360 :
361 0 : IF (nvac == 1 .AND. ipos /= pTOP)&
362 0 : & CALL juDFT_error("ERROR reading ",calledby ="efield")
363 :
364 0 : action = pADD
365 0 : IF (INDEX (str, 'replace') > 0) action = pREPLACE
366 0 : iopt = pALL
367 0 : IF (INDEX (str, ' zero') > 0) iopt = pZERO
368 0 : IF (INDEX (str, 'nonzero') > 0) iopt = pNONZERO
369 :
370 0 : mask = .false.
371 0 : nx = UBOUND (mask, dim=1)
372 0 : ny = UBOUND (mask, dim=2)
373 :
374 0 : IF (tag == 'rect') THEN
375 : mask(MAX (FLOOR(x*nx-0.5)+2,1) : MIN (FLOOR((x+w)*nx+0.5),nx),&
376 : & MAX (FLOOR(y*ny-0.5)+2,1) : MIN (FLOOR((y+h)*ny+0.5),ny))&
377 0 : & = .true.
378 0 : WRITE (oUnit, *) tag, pos2str (ipos), x, y, h, w, charge,&
379 0 : & action2str (action), opt2str (iopt)
380 0 : ELSE IF (tag == 'rectsinx') THEN
381 : mask(MAX (FLOOR(x*nx-0.5)+2,1) : MIN (FLOOR((x+w)*nx+0.5),nx),&
382 : & MAX (FLOOR(y*ny-0.5)+2,1) : MIN (FLOOR((y+h)*ny+0.5),ny))&
383 0 : & = .true.
384 0 : WRITE (oUnit, *) tag, pos2str (ipos), x, y, h, w, charge, order,&
385 0 : & shift, action2str (action), opt2str (iopt)
386 0 : ELSE IF (tag == 'circ') THEN
387 0 : DO iy = 1, ny
388 0 : DO ix = 1, nx
389 0 : IF (SQRT ((REAL (ix)/nx - x)**2 + (REAL(iy)/ny - y)**2)&
390 : & <= radius)&
391 0 : & mask(ix, iy) = .true.
392 : END DO
393 : END DO
394 0 : WRITE (oUnit, *) tag, pos2str (ipos), x, y, radius, charge,&
395 0 : & action2str (action), opt2str (iopt)
396 0 : ELSE IF (tag == 'datafile') THEN
397 0 : WRITE (oUnit, '(1x,a,1x,a)', advance="no") tag, pos2str (ipos)
398 0 : call readDataFile (str, data)
399 0 : WRITE (oUnit, *) action2str (action), opt2str (iopt)
400 0 : mask = .true.
401 0 : ELSE IF (tag == 'polygon') THEN
402 0 : DO iy = 1, ny
403 0 : DO ix = 1, nx
404 : mask(ix, iy) = in_polygon (np, poly, (/real(ix-1)/nx,&
405 0 : & real(iy-1)/ny/))
406 : END DO
407 : END DO
408 : ELSE
409 : CALL juDFT_error("Internal error in read_shape",calledby&
410 0 : & ="efield")
411 : END IF
412 :
413 0 : IF (action == pADD) THEN
414 : ! All relevant vacua
415 0 : DO ivac = 2 - MOD (ipos, 2), MIN (ipos, 2)
416 0 : IF (iopt == pZERO) THEN
417 0 : mask = mask .and. (E%sigEF(:,:,ivac) == 0.0)
418 0 : ELSE IF (iopt == pNONZERO) THEN
419 0 : mask = mask .and. (E%sigEF(:,:,ivac) /= 0.0)
420 : END IF
421 :
422 0 : IF (tag == 'rectsinx') THEN
423 : ! number of grid points in y direction
424 0 : IF (E%dirichlet) THEN
425 : cnt = 1
426 : ELSE
427 0 : cnt = MAXVAL (COUNT (mask, dim=2))
428 : END IF
429 0 : i = MAX (FLOOR(x*nx-0.5)+2,1)
430 0 : j = MIN (FLOOR((x+w)*nx+0.5),nx)
431 0 : DO ix = i, j
432 0 : tmp = REAL (ix-1)/REAL (j) ! Range: [0,1)
433 0 : DO iy = 1, ny
434 0 : IF (.NOT. mask (ix, iy)) CYCLE
435 : E%sigEF(ix,iy,ivac) = E%sigEF(ix,iy,ivac)&
436 0 : & + charge/cnt * SIN (TWO_PI*(order*tmp+shift))
437 : END DO
438 : END DO
439 0 : ELSE IF (tag == 'datafile') THEN
440 0 : WHERE (mask) E%sigEF(:,:,ivac) = E%sigEF(:,:,ivac) + data
441 : ELSE ! circ, rect
442 0 : IF (E%dirichlet) THEN
443 : cnt = 1
444 : ELSE
445 0 : cnt = COUNT (mask)
446 : END IF
447 : WHERE (mask) E%sigEF(:,:,ivac) = E%sigEF(:,:,ivac)&
448 0 : & + charge/cnt
449 : END IF
450 : END DO ! ivac
451 : ELSE ! pREPLACE
452 : ! All relevant vacua
453 0 : DO ivac = 2 - MOD (ipos, 2), MIN (ipos, 2)
454 0 : IF (iopt == pZERO) THEN
455 0 : mask = mask .and. (E%sigEF(:,:,ivac) == 0.0)
456 0 : ELSE IF (iopt == pNONZERO) THEN
457 0 : mask = mask .and. (E%sigEF(:,:,ivac) /= 0.0)
458 : END IF
459 :
460 0 : IF (tag == 'rectsinx') THEN
461 : ! number of grid points in y direction
462 0 : IF (E%dirichlet) THEN
463 : cnt = 1
464 : ELSE
465 0 : cnt = MAXVAL (COUNT (mask, dim=2))
466 : END IF
467 0 : i = MAX (FLOOR(x*nx-0.5)+2,1)
468 0 : j = MIN (FLOOR((x+w)*nx+0.5),nx)
469 0 : DO ix = i, j
470 0 : tmp = REAL (ix-1)/REAL (j) ! Range: [0,1)
471 0 : DO iy = 1, ny
472 0 : IF (.NOT. mask (ix, iy)) CYCLE
473 : E%sigEF(ix,iy,ivac) =&
474 0 : & + charge/cnt * SIN (TWO_PI*(order*tmp+shift))
475 : END DO
476 : END DO
477 0 : ELSE IF (tag == 'datafile') THEN
478 0 : WHERE (mask) E%sigEF(:,:,ivac) = data
479 : ELSE ! circ, rect
480 0 : IF (E%dirichlet) THEN
481 : cnt = 1
482 : ELSE
483 0 : cnt = COUNT (mask)
484 : END IF
485 0 : WHERE (mask) E%sigEF(:,:,ivac) = charge/cnt
486 : END IF
487 : END DO ! ivac
488 : END IF
489 0 : END SUBROUTINE read_shape
490 :
491 : ! Read electric-field data points from "file"
492 : ! Format: First line, number of x and y points
493 : ! Second line: charge(x=1,y=1)
494 : ! Third line: charge(x=1,y=2)
495 : ! etc.
496 0 : SUBROUTINE readDataFile (str, data)
497 : CHARACTER(*), intent(in) :: str
498 : REAL, intent(out) :: data(:,:)
499 :
500 : integer :: i, j, nx, ny
501 : logical :: l_exist
502 0 : character(len (str)) :: file, dummy
503 :
504 0 : read(str, *) dummy, dummy, file
505 :
506 0 : nx = ubound(data, dim=1)
507 0 : ny = ubound(data, dim=2)
508 :
509 0 : INQUIRE (file=file, exist=l_exist)
510 0 : IF (.not. l_exist) GOTO 110
511 :
512 0 : OPEN (1243, file=file, status='old')
513 0 : READ (1243, *, end=110) i, j ! nx ny
514 0 : IF (i /= nx .or. j /= ny) THEN
515 0 : WRITE (*,'(3a)') "ERROR: Invalid number of points in '",&
516 0 : & trim(file), "' during electrical field dataFile"
517 0 : WRITE (*,"(a,i0,a,i0,a)") "Expected: ", nx, ",", ny
518 0 : WRITE (*,"(a,i0,a,i0,a)") "Read : ", i, ", ", j
519 : CALL juDFT_error("ERROR: dataFile of electric fiel&
520 0 : &d",calledby ="efield")
521 : END IF
522 :
523 0 : DO j = 1, ny
524 0 : DO i = 1, nx
525 0 : READ (1243, *) data(i,j)
526 : END DO
527 : END DO
528 0 : CLOSE (1243)
529 :
530 0 : WRITE (oUnit, '(1x,3a,1x)', advance="no") '"', trim (file), '"'
531 0 : IF (INDEX (lower_case (str), 'zero_avg') > 0) THEN
532 0 : data(:,:) = data - sum(data)/(nx*ny)
533 0 : WRITE (oUnit, '(a,1x)', advance="no") "zero_avg"
534 : END IF
535 :
536 0 : RETURN ! No error
537 :
538 : ! ERROR
539 0 : 110 WRITE (*,'(3a)') "ERROR: While reading '",trim(file),&
540 0 : & "' during electrical field dataFile"
541 0 : WRITE (*,"(a,i0,a,i0,a)") "Expecting ", nx, ",", ny, " points"
542 : CALL juDFT_error("ERROR: dataFile of electric field",calledby&
543 0 : & ="efield")
544 : END SUBROUTINE readDataFile
545 :
546 :
547 : ! Give a polygon and a trial point. Now a line is taken from the edge
548 : ! of the grid and it is counted how many times an edges of the polygon
549 : ! is crossed. If it was an odd time until one reaches the point, it is
550 : ! in the polygon - otherwise not.
551 : ! Note: The technique does not always find the optimal polygon.
552 : !
553 : ! Cf. http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
554 : !
555 0 : PURE FUNCTION in_polygon (n, poly, point)
556 : INTEGER, intent(in) :: n
557 : REAL, intent(in) :: poly(2,n), point(2)
558 : LOGICAL :: in_polygon
559 : INTEGER :: i, j
560 :
561 0 : i = 1
562 0 : j = n
563 0 : in_polygon = .FALSE.
564 0 : DO
565 0 : IF (i > n) EXIT
566 0 : IF ((poly(2,i) > point(2)) .NEQV. (poly(2,j) > point(2))) THEN
567 0 : IF (point(1) < (poly(1,j) - poly(1,i))&
568 : & * (point(2)-poly(2,i)) / (poly(2,j)-poly(2,i))&
569 : & + poly(1,i))&
570 0 : & in_polygon = .NOT. in_polygon
571 : END IF
572 0 : j = i
573 0 : i = i+1
574 : END DO
575 0 : END FUNCTION in_polygon
576 :
577 0 : PURE FUNCTION pos2str (ipos)
578 : INTEGER, INTENT(IN) :: ipos
579 : CHARACTER(6) :: pos2str
580 0 : IF (ipos == pTOP) THEN
581 0 : pos2str = 'Top'
582 0 : ELSE IF (ipos == pBOT) THEN
583 0 : pos2str = 'Bot'
584 : ELSE
585 0 : pos2str = 'TopBot'
586 : END IF
587 0 : END FUNCTION pos2str
588 :
589 :
590 0 : PURE FUNCTION action2str (action)
591 : INTEGER, INTENT(IN) :: action
592 : CHARACTER(7) :: action2str
593 0 : IF (action == pADD) THEN
594 0 : action2str = 'Add'
595 : ELSE
596 0 : action2str = 'Replace'
597 : END IF
598 0 : END FUNCTION action2str
599 :
600 :
601 0 : PURE FUNCTION opt2str (iopt)
602 : INTEGER, INTENT(IN) :: iopt
603 : CHARACTER(7) :: opt2str
604 0 : IF (iopt == pALL) THEN
605 0 : opt2str = 'All'
606 0 : ELSE IF (iopt == pZERO) THEN
607 0 : opt2str = 'Zero'
608 : ELSE
609 0 : opt2str = 'NonZero'
610 : END IF
611 0 : END FUNCTION opt2str
612 :
613 :
614 80 : SUBROUTINE print_efield (unit, E, area, nvac, amat,dvac)
615 : INTEGER, INTENT(IN) :: unit, nvac
616 : TYPE(t_efield), INTENT(IN) :: E
617 : REAL, INTENT(IN) :: area, amat(3,3),dvac
618 :
619 : ! electric constant
620 : REAL, PARAMETER :: epsilon0 = 8.854187817620e-12 ! F/m
621 : ! Bohr radius
622 : REAL, PARAMETER :: a0 = 0.52917720859e-10 ! m
623 : ! elemental charge
624 : REAL, PARAMETER :: ec = 1.602176487e-19 ! C
625 :
626 : ! htr -> electron Volt
627 : REAL, PARAMETER :: htr_eV = 27.21138386 ! eV
628 :
629 : ! Conversion of surface charge sigma to an electric field
630 : ! near a plate: E = -sigma[1/m^2]*e/eps0
631 : ! = -sigma[1/a0^2]*e/(eps0 * a0^2)
632 : ! Assuming field in V/cm and charge in atomic units
633 : ! The 10^-2 is due to centimetres
634 : !
635 : ! NOTE: This equation assumes that one has two plates, one
636 : ! charged "sigma" and one "-sigma"; then E = sigma/eps0.
637 : ! If one plate were 0 and the other sigma, one had to use
638 : ! E = sigma/(2*eps0)
639 : ! Cf. http://hyperphysics.phy-astr.gsu.edu/hbase/electric/elesht.html
640 : !
641 : ! BE CAREFUL TO AVOID DOUBLE COUNTING
642 : REAL, PARAMETER :: sig_to_E = -ec/(epsilon0*a0**2)*1e-2
643 :
644 : REAL :: tmp, pt_rel(3), pt_abs(3)
645 : INTEGER :: ivac, i, j, nx, ny
646 :
647 : IF (SUM (ABS (E%sig_b)) < 1e-15&
648 : &.AND. ABS (E%sigma) < 1e-15&
649 : &.AND. .NOT. ALLOCATED (E%sigEF)&
650 240 : & .AND. .NOT. E%dirichlet) RETURN
651 :
652 0 : WRITE (unit,'(3x,a,2(f12.5,a))')'z-z1 of external sheet =',&
653 0 : & E%zsigma, ' a0 = ', E%zsigma*a0*1e10, ' A'
654 :
655 :
656 0 : IF (E%dirichlet) THEN
657 :
658 0 : IF (E%sigma > 0)&
659 0 : & WRITE (unit,'(3x,a,f12.5)') 'Surplus charge: ', 2.0*E%sigma
660 :
661 0 : IF (ALLOCATED(E%sigEF))&
662 0 : & WRITE (unit,'(3x,a)') 'Average potential:'
663 0 : WRITE (unit, '(3x,a,f12.5,a, f12.5,a)') 'on sheet 1: ',&
664 0 : & E%sig_b(1),' htr = ',E%sig_b(1)*hartree_to_ev_const,' V'
665 0 : IF (nvac > 1) THEN
666 0 : WRITE (unit, '(3x,a,f12.5,a, f12.5,a)') 'on sheet 2: ',&
667 0 : & E%sig_b(2),' htr = ',E%sig_b(2)*hartree_to_ev_const,' V'
668 :
669 : WRITE (unit,'(3x,a,f14.5,a)')&
670 0 : & 'Average field (plate to plate):',&
671 0 : & (E%sig_b(2)-E%sig_b(1))/((2*E%zsigma+Dvac)*a0*100),&
672 0 : & ' V/cm'
673 : END IF ! nvac > 1
674 :
675 : ELSE ! Neumann
676 :
677 0 : IF (ALLOCATED(E%sigEF))&
678 0 : & WRITE (unit,'(3x,a)') 'Average charges:'
679 0 : tmp = E%sigma+E%sig_b(1)
680 : WRITE (unit, '((3x,a,f12.5,5x,a,f12.5,a))')&
681 0 : & 'charge on external sheet 1 =', tmp,&
682 0 : & '(surface density=', tmp/area, ' e/a.u.**2)'
683 0 : tmp = tmp/area*sig_to_E
684 : WRITE (unit, '(3x,a,e13.5,5x,a)')&
685 0 : & 'external field on sheet 1 =', tmp, 'V/cm'
686 :
687 0 : IF (nvac > 1) THEN
688 0 : tmp = E%sigma+E%sig_b(2)
689 : WRITE (unit, '((3x,a,f12.5,5x,a,f12.5,a))')&
690 0 : & 'charge on external sheet 2 =', tmp,&
691 0 : & '(surface density=', tmp/area, ' e/a.u.**2)'
692 0 : tmp = tmp/area*sig_to_E
693 : WRITE (unit, '(3x,a,e13.5,5x,a)')&
694 0 : & 'external field on sheet 2 =', tmp, 'V/cm'
695 :
696 : ! Cf. comment above, where "sig_to_E" is defined
697 : WRITE (unit, '(a,/,a)') 'NOTE: The equation for the E field'&
698 0 : & // ' assumes two oppositely charged plates.',&
699 : & ' You may need to divide by two before '&
700 0 : & //'summing the external fields to avoid double counting'
701 : ELSE
702 : WRITE (unit, '(a)') 'NOTE: The equation for the E field '&
703 : & // 'already assumes two oppositely charged plates - '&
704 0 : & // 'be careful to avoid double counting'
705 : END IF ! nvac > 1
706 :
707 : END IF ! Neumann (vs. Dirichlet)
708 :
709 0 : IF (ALLOCATED (E%sigEF) .AND. E%plot_charge) THEN
710 0 : nx = UBOUND (E%sigEF, dim=1) - 1
711 0 : ny = UBOUND (E%sigEF, dim=2) - 1
712 0 : DO ivac = 1, nvac
713 0 : IF (ivac == 1) THEN
714 0 : OPEN (748, file='efield-1.dat')
715 : ELSE
716 0 : OPEN (748, file='efield-2.dat')
717 : END IF
718 0 : IF (E%dirichlet) THEN
719 0 : WRITE (748, '(3a)') '# X[a0] Y[a0] X[rel] ',&
720 0 : & 'Y[rel] potential[htr] potential[V]'
721 : ELSE ! Neumann
722 0 : WRITE (748, '(3a)') '# X[a0] Y[a0] X[rel] ',&
723 0 : & 'Y[rel] charge per grid point',&
724 0 : & ' E_field (V/cm) in vicinity of the grid point'
725 : END IF
726 0 : pt_rel(3) = 0.0
727 0 : DO i = 1, nx+1
728 0 : pt_rel(1) = REAL(i-1)/nx
729 0 : DO j = 1, ny+1
730 0 : pt_rel(2) = REAL(j-1)/ny
731 0 : pt_abs=matmul(amat,pt_rel)
732 0 : IF (E%dirichlet) THEN
733 : WRITE (748, '(4f12.5,2g16.5)')&
734 0 : & pt_abs(1:2), pt_rel(1:2),&
735 : & E%sigEF(i,j,ivac)&
736 0 : & + E%sig_b(ivac),&
737 : & (E%sigEF(i,j,ivac)&
738 0 : & + E%sig_b(ivac))*hartree_to_ev_const
739 : ELSE ! Neumann
740 : WRITE (748, '(4f12.5,2g16.5)')&
741 0 : & pt_abs(1:2), pt_rel(1:2),&
742 : & E%sigEF(i,j,ivac)&
743 0 : & + (E%sigma + E%sig_b(ivac))/((nx+1)*(ny+1)),&
744 : & (E%sigEF(i,j,ivac)*((nx+1)*(ny+1))&
745 0 : & + E%sigma + E%sig_b(ivac))/area*sig_to_E
746 : END IF
747 : END DO
748 0 : WRITE (748, *)
749 : END DO
750 0 : CLOSE (748)
751 : END DO
752 : END IF
753 : END SUBROUTINE print_efield
754 :
755 0 : SUBROUTINE V_seg_EF(&
756 : & efield,&
757 : & vacuum, stars)
758 : USE m_fft2d
759 : use m_types
760 : ! Dummy variables:
761 : TYPE(t_efield), INTENT(INOUT) :: efield
762 : TYPE(t_vacuum), INTENT(IN) :: vacuum
763 : TYPE(t_stars), INTENT(IN) :: stars
764 :
765 : ! Local variables:
766 : INTEGER :: i, ivac
767 : REAL :: fg, fgi
768 0 : COMPLEX :: fgfull(stars%ng2)
769 0 : REAL :: rhoRS(3*stars%mx1,3*stars%mx2), rhoRSimag(3*stars%mx1,3*stars%mx2) ! Real space density
770 :
771 0 : ALLOCATE (efield%rhoEF (stars%ng2-1, vacuum%nvac))
772 :
773 0 : DO ivac = 1, vacuum%nvac
774 0 : rhoRSimag = 0.0 ! Required in the loop. fft2d overrides this
775 : ! The fft2d algorithm puts the normalization to the isn=-1
776 : ! (r->k) transformation, thus we need to multiply by
777 : ! ifft2 = 3*k1d*3*k2d to compensate.
778 0 : IF (efield%dirichlet) THEN
779 0 : rhoRS(:,:) = efield%sigEF(:,:,ivac)
780 : ELSE
781 0 : rhoRS(:,:) = efield%sigEF(:,:,ivac)*(3*stars%mx1*3*stars%mx2)
782 : END IF
783 :
784 : ! FFT rhoRS(r_2d) -> rhoEF(g_2d)
785 : CALL fft2d (&
786 : & stars,&
787 : & rhoRS, rhoRSimag,&
788 0 : & fgfull, -1)
789 0 : efield%rhoEF(:,ivac) = fgfull(2:)
790 : ! FFT gives the the average charge per grid point
791 : ! while sig_b stores the (total) charge per sheet
792 0 : IF (efield%dirichlet .and. ABS (REAL(fgfull(1))) > 1.0e-15) THEN
793 0 : PRINT *, 'INFO: Difference of average potential: fg=',&
794 0 : & REAL(fgfull(1)),', sig_b=', efield%sig_b(ivac),&
795 0 : & ", ivac=", ivac
796 0 : ELSE IF (ABS (REAL(fgfull(1))/(3*stars%mx1*3*stars%mx2)) > 1.0e-15) THEN
797 0 : PRINT *, 'INFO: Difference of average potential: fg=',&
798 0 : & REAL(fgfull(1))/(3*stars%mx1*3*stars%mx2),', sig_b=', efield%sig_b(ivac),&
799 0 : & ", ivac=", ivac
800 : END IF
801 : END DO ! ivac
802 0 : END SUBROUTINE V_seg_EF
803 :
804 :
805 0 : SUBROUTINE print_rhoEF(&
806 : & efield, k1d, k2d, nvac, n2d, nop, nop2,&
807 0 : & nq2, kv2, mrot, symor, tau, invtab, area,&
808 0 : & nstr2, amat)
809 : USE m_starf, ONLY: starf2
810 : ! Arguments
811 : TYPE(t_efield), INTENT(IN) :: efield
812 : INTEGER, INTENT(IN) :: k1d, k2d, nvac, n2d, nq2, nop, nop2
813 : REAL, INTENT (IN) :: area
814 : LOGICAL, INTENT (IN) :: symor
815 : INTEGER, INTENT (IN) :: nstr2(n2d), kv2(2,n2d), mrot(3,3,nop),&
816 : & invtab(nop)
817 : REAL, INTENT (IN) :: tau(3,nop), amat(3,3)
818 :
819 : ! Local variables
820 : real :: pt_rel(3), pt_abs(3), rcc(3)
821 : integer :: ix, iy, ivac, k, nix, niy
822 : real :: rhoOut, rhoOutIm
823 0 : complex :: sf2(n2d)
824 :
825 0 : nix = 3*k1d-1
826 0 : niy = 3*k2d-1
827 :
828 0 : DO ivac = 1,nvac
829 0 : IF (ivac == 1) THEN
830 0 : OPEN (754,file='efield-fft-1.dat')
831 : ELSE
832 0 : OPEN (754,file='efield-fft-2.dat')
833 : END IF
834 :
835 0 : IF (efield%dirichlet) THEN
836 0 : WRITE (754, '(2a)') '# X[a0] Y[a0] X[rel] Y[rel] ',&
837 0 : & 'potential per grid point (Re, Im, Abs), grid (ix, iy)'
838 : ELSE ! Neumann
839 0 : WRITE (754, '(2a)') '# X[a0] Y[a0] X[rel] Y[rel] ',&
840 0 : & 'charge per grid point (Re, Im, Abs), grid (ix, iy)'
841 : END IF
842 0 : pt_rel(3) = 0.0
843 0 : DO ix = 0, nix
844 0 : pt_rel(1) = REAL (ix)/REAL (nix)
845 0 : DO iy = 0, niy
846 0 : pt_rel(2) = REAL (iy)/REAL (niy)
847 : CALL starf2 (&
848 : & nop2,nq2,kv2,mrot,symor,tau,pt_rel,invtab,&
849 0 : & sf2)
850 0 : IF (efield%dirichlet) THEN
851 0 : rhoOut = efield%sig_b(ivac)
852 : ELSE
853 0 : rhoOut = (efield%sigma + efield%sig_b(ivac))/area
854 : END IF
855 0 : rhoOutIm = 0.0
856 : ! As we moved the normalization to g-space while the
857 : ! fft2d routine works in r space, we need to compensate
858 : ! by dividing by ifft2 = 3*k1d*3*k2d
859 0 : DO k = 2, nq2
860 0 : IF (efield%dirichlet) THEN
861 : rhoOut = rhoOut&
862 : & + REAL (efield%rhoEF(k-1,ivac)*sf2(k))&
863 0 : & * nstr2(k)
864 : rhoOutIm = rhoOutIm&
865 : & + AIMAG (efield%rhoEF(k-1,ivac)*sf2(k))&
866 0 : & * nstr2(k)
867 : ELSE
868 : rhoOut = rhoOut&
869 : & + REAL (efield%rhoEF(k-1,ivac)*sf2(k))&
870 0 : & * nstr2(k)/(3*k1d*3*k2d)
871 : rhoOutIm = rhoOutIm&
872 : & + AIMAG (efield%rhoEF(k-1,ivac)*sf2(k))&
873 0 : & * nstr2(k)/(3*k1d*3*k2d)
874 : END IF
875 : END DO ! k
876 0 : pt_abs=matmul(amat,pt_rel)
877 : WRITE (754,'(4f14.7,3g16.7,2i6)')&
878 0 : & pt_abs(1:2),pt_rel(1:2),&
879 0 : & rhoOut, rhoOutIm, SQRT (rhoOut**2+rhoOutIm**2), ix, iy
880 : END DO ! iy
881 0 : WRITE (754,*)
882 : END DO ! ix
883 0 : CLOSE (754)
884 : END DO ! ivac
885 0 : END SUBROUTINE print_rhoEF
886 : END SUBROUTINE e_field
887 :
888 : SUBROUTINE read_namelist (iou, E, eV)
889 : USE m_types, only: t_efield
890 : INTEGER, INTENT(IN) :: iou
891 : TYPE(t_efield), INTENT(INOUT) :: E
892 : LOGICAL, INTENT(INOUT) :: eV
893 :
894 : REAL :: zsigma, sig_b(2)
895 : LOGICAL :: plot_charge
896 : LOGICAL :: plot_rho
897 : LOGICAL :: autocomp
898 : LOGICAL :: dirichlet
899 : NAMELIST /efield/ zsigma, sig_b, plot_charge, plot_rho,&
900 : & autocomp, dirichlet, eV
901 :
902 : zsigma = E%zsigma
903 : plot_charge = E%plot_charge
904 : plot_rho = E%plot_rho
905 : autocomp = E%autocomp
906 : dirichlet = E%dirichlet
907 : sig_b = E%sig_b
908 :
909 : READ (iou,nml=efield)
910 :
911 : E%zsigma = zsigma
912 : E%plot_charge = plot_charge
913 : E%plot_rho = plot_rho
914 : E%autocomp = autocomp
915 : E%dirichlet = dirichlet
916 : E%sig_b = sig_b
917 : END SUBROUTINE read_namelist
918 :
919 :
920 0 : ELEMENTAL FUNCTION lower_case(string)
921 : CHARACTER(len=*), INTENT(IN) :: string
922 : CHARACTER(len=len(string)) :: lower_case
923 :
924 : INTEGER :: i
925 :
926 0 : DO i = 1, LEN (string)
927 : IF ( IACHAR ('A') <= IACHAR (string(i:i)) &
928 0 : .and. IACHAR ('Z') >= IACHAR (string(i:i))) THEN
929 : lower_case(i:i) = ACHAR (IACHAR (string(i:i))&
930 0 : + IACHAR ('a') - IACHAR ('A'))
931 : ELSE
932 0 : lower_case(i:i) = string(i:i)
933 : END IF
934 : END do
935 0 : END FUNCTION lower_case
936 : END MODULE m_efield
|
