LCOV - code coverage report
Current view: top level - init - efield.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 50 374 13.4 %
Date: 2019-09-08 04:53:50 Functions: 3 12 25.0 %

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

Generated by: LCOV version 1.13