LCOV - code coverage report
Current view: top level - vgen - od_vvacis.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 240 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_od_vvacis
       8             : CONTAINS
       9           0 :   SUBROUTINE od_vvacis(&
      10             :        &     n2d_1,input,vacuum,nq2_1,&
      11           0 :        &     kv2_1,cell,MM,stars,&
      12             :        &     nstr2_1,&
      13             :        &     oneD,&
      14           0 :        &     rht,rhtxy,psq,vz,sym,&
      15           0 :        &     vxy,vpw)
      16             : 
      17             :     !     generates m^2+gz^2.neq.0 coefficients of the vacuum potential
      18             :     !     and fourier coefficients of the interstitial potential
      19             :     !     in the case of 1-dimensional calculations
      20             :     !     based on the using of the green's function for the 
      21             :     !     radial equation for the potential:
      22             :     !     [1/r(d/dr(r*d/dr))-m^2/r^2-gz^2]V(r) = 4pi*rho(r)
      23             :     !     when gz.neq.0 green's function is represented in
      24             :     !     the following way: 
      25             :     !          gf = 4pi*I_m(gz*r<)*K_m(gz*r>)
      26             :     !     where I and K are modified bessel funcions
      27             :     !     in the case of gz.eq.0 the green's function is:
      28             :     !            gf = 2pi*((r</r>)^m)/m
      29             :     !                 Y.Mokrousov, autumn 2002
      30             : 
      31             :     !     Fully symmetrized version, Mai 2003
      32             : 
      33             :     USE m_constants
      34             :     USE m_od_cylbes
      35             :     USE m_modcyli
      36             :     USE m_modcylk
      37             :     USE m_vacp5_0
      38             :     USE m_vacp5_z
      39             :     USE m_visp5_0
      40             :     USE m_visp5_z
      41             :     USE m_angle
      42             :     USE m_qsf
      43             :     USE m_fft2d
      44             : 
      45             : 
      46             :     USE m_types
      47             :     IMPLICIT NONE
      48             : 
      49             :     TYPE(t_input),INTENT(IN)   :: INPUT
      50             : 
      51             :     TYPE(t_oneD),INTENT(IN)   :: oneD
      52             : 
      53             :     TYPE(t_vacuum),INTENT(IN)   :: vacuum
      54             : 
      55             :     TYPE(t_sym),INTENT(IN)   :: sym
      56             : 
      57             :     TYPE(t_stars),INTENT(IN)   :: stars
      58             : 
      59             :     TYPE(t_cell),INTENT(IN)   :: cell
      60             : 
      61             :     INTEGER, INTENT (IN) :: n2d_1  ,MM 
      62             :     INTEGER, INTENT (IN) :: nq2_1  
      63             :     INTEGER, INTENT (IN) :: nstr2_1(n2d_1)
      64             :     INTEGER, INTENT (IN) :: kv2_1(2,n2d_1) 
      65             :     COMPLEX, INTENT (INOUT) :: psq(stars%ng3)
      66             :     REAL,    INTENT (IN) :: vz(vacuum%nmzd,2) 
      67             :     REAL,    INTENT (IN) :: rht(vacuum%nmzd,2)
      68             :     COMPLEX, INTENT (IN) :: rhtxy(vacuum%nmzxyd,n2d_1-1,2)
      69             :     COMPLEX, INTENT (OUT):: vxy(vacuum%nmzxyd,n2d_1-1,2)
      70             :     COMPLEX, INTENT (OUT):: vpw(stars%ng3)
      71             : 
      72             :     !     local
      73             :     INTEGER :: m
      74             :     !------> reciprocal vectors staff
      75             : 
      76             :     INTEGER irec2,irec3,k1,k2,gzi,gzi1,k3,gzmin
      77             :     REAL    gz,gx,gy,g,g2 
      78             : 
      79             :     !------> different garbage
      80             : 
      81             :     REAL, PARAMETER :: tol_21 = 1.0e-21
      82             :     INTEGER  imz,imz1,i,ivac,iirec1,m1,j,irc1,im,irc
      83             :     INTEGER ix,iy,iz,s,n
      84             :     REAL    x,y,r,z,b,q,zf
      85             :     REAL    mult
      86             :     REAL    ani1,ani2,rhti
      87             :     INTEGER ivfft1,ivfft2,ivfft2d
      88             :     COMPLEX ic
      89             : 
      90             :     !------> different staff with special functions
      91             : 
      92           0 :     REAL  IR(1:stars%mx3,0:MM)
      93           0 :     REAL, ALLOCATABLE :: II(:),KK(:)
      94             :     REAL  fJ,fJ2,fJ1,iJ2,IIIR,fJ3
      95           0 :     REAL, ALLOCATABLE :: III(:),IIII(:,:,:)
      96           0 :     REAL, ALLOCATABLE :: fJJ(:,:),iJJ(:,:)
      97             : 
      98             :     !------> used for the optimization 
      99             : 
     100             :     INTEGER irec1(4),l,lmin ,l1
     101             : 
     102             :     !------> some factors used for the construction of the ch. density
     103             : 
     104           0 :     REAL, ALLOCATABLE :: fact(:)
     105             :     COMPLEX aa,a
     106             : 
     107             :     !------> values of the potential on the vacuum boundary
     108             : 
     109             :     COMPLEX val_help
     110           0 :     COMPLEX, ALLOCATABLE :: val(:),val_m(:,:)
     111             : 
     112             :     !------> Charge density
     113             : 
     114           0 :     COMPLEX, ALLOCATABLE :: rxy(:)
     115             : 
     116             :     !---------------> POTENTIALS
     117             : 
     118             :     !-> interstitial potential by its gz,m - components and total on the 
     119             :     !-> real grid 
     120             : 
     121           0 :     COMPLEX, ALLOCATABLE :: vis(:,:,:),vis_tot(:,:)
     122             : 
     123             :     !-> potential in the vacuum caused by the vacuum charge density      
     124             : 
     125           0 :     COMPLEX, ALLOCATABLE :: pvac(:)
     126             : 
     127             :     !-> potential in the vacuum caused by the interst. density
     128             : 
     129           0 :     COMPLEX, ALLOCATABLE :: pint(:)
     130             : 
     131             :     !-> radial components of the potential caused by the vacuum  
     132             :     !-> density on the real grid  
     133             : 
     134           0 :     COMPLEX, ALLOCATABLE :: vis_help(:,:),vpw_help(:)
     135             : 
     136             :     !-> real grid vis_z,vis_0 for the fft2d
     137             : 
     138           0 :     REAL,ALLOCATABLE :: af2(:),bf2(:)
     139             : 
     140             :     !-> radial grids
     141             : 
     142           0 :     INTEGER, ALLOCATABLE :: rmap(:,:)
     143           0 :     REAL   , ALLOCATABLE :: rr(:)
     144             :     INTEGER nrd
     145             : 
     146             :     !--> time
     147             : 
     148             :     REAL    gxy0,fxy0,phi
     149           0 :     COMPLEX gxy(stars%ng2-1)
     150           0 :     COMPLEX fxy(stars%ng2-1)
     151             : 
     152             :     INTRINSIC REAL,aimag
     153             : 
     154             :     !--------- preparations ---------->
     155             : 
     156             :     ALLOCATE ( KK(vacuum%nmzxyd),II(vacuum%nmzxyd),III(9*stars%mx1*stars%mx2),&
     157             :          &     IIII(9*stars%mx1*stars%mx3,1:stars%mx3,0:MM),&
     158             :          &     rmap(0:3*stars%mx1-1,0:3*stars%mx2-1),rr(1:9*stars%mx1*stars%mx2),&
     159             :          &     fact(vacuum%nmzxyd),val(n2d_1),fJJ(0:MM+1,stars%ng2),&
     160             :          &     val_m(-stars%mx3:stars%mx3,-MM:MM),rxy(vacuum%nmzxyd),iJJ(0:MM+1,1:stars%mx3),&
     161             :          &     vis(0:3*stars%mx1-1,0:3*stars%mx2-1,n2d_1),&
     162             :          &     vis_tot(0:3*stars%mx1-1,0:3*stars%mx2-1),pvac(vacuum%nmzxyd),&
     163             :          &     pint(vacuum%nmzxyd),vis_help(0:3*stars%mx1-1,0:3*stars%mx2-1),&
     164           0 :          &     af2(0:9*stars%mx1*stars%mx2-1),bf2(0:9*stars%mx1*stars%mx2-1),vpw_help(stars%ng3) )
     165             : 
     166           0 :     ivfft2d = 9*stars%mx1*stars%mx2
     167           0 :     ic = CMPLX(0.,1.)
     168           0 :     ivfft1 = 3*stars%mx1
     169           0 :     ivfft2 = 3*stars%mx2
     170           0 :     ani1 = 1./REAL(ivfft1)
     171           0 :     ani2 = 1./REAL(ivfft2)
     172             : 
     173             :     !--------- initializations -------->
     174             : 
     175             :     !----> vpw in the '1st aproximation' (V - tilde)
     176             : 
     177           0 :     vpw(1) = CMPLX(0.,0.)
     178             : 
     179           0 :     DO irec3 = 2,stars%ng3
     180             : 
     181           0 :        g = stars%sk3(irec3)
     182             : 
     183           0 :        vpw(irec3) = fpi_const*psq(irec3)/(g*g)
     184             : 
     185             :     ENDDO
     186             : 
     187           0 :     DO irc1 = 2,nq2_1
     188           0 :        DO i = 1,vacuum%nmzxy
     189           0 :           vxy(i,irc1-1,1) = CMPLX(0.,0.)
     190             :        END DO
     191             :     END DO
     192             : 
     193             :     !----> values of the potential in the 1st approximation on the boundary
     194             :     !----> if nstr2.ne.1 then it should be changed!!!
     195             : 
     196           0 :     DO m = 0,MM+1
     197           0 :        DO k3 = 1,stars%mx3
     198           0 :           CALL modcyli(m,cell%bmat(3,3)*k3*cell%z1,iJJ(m,k3))
     199             :        END DO
     200           0 :        DO irec2 = 1,stars%ng2
     201           0 :           g2 = stars%sk2(irec2)
     202           0 :           IF (irec2.NE.0) THEN
     203           0 :              CALL od_cylbes(m,g2*cell%z1,fJJ(m,irec2))
     204             :           END IF
     205             :        END DO
     206             :     END DO
     207             : 
     208           0 :     DO irc1 = 1,nq2_1
     209           0 :        val(irc1) = CMPLX(0.,0.)
     210           0 :        m = kv2_1(2,irc1)
     211           0 :        IF (m.LT.0) THEN
     212           0 :           mult = REAL((-1)**m)
     213             :        ELSE
     214             :           mult = 1.
     215             :        END IF
     216           0 :        k3 = kv2_1(1,irc1)
     217           0 :        DO irec2 = 1,stars%ng2
     218           0 :           phi = stars%phi2(irec2)
     219           0 :           irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
     220           0 :           IF (irec3.NE.0) THEN
     221             :              val(irc1) = val(irc1) +&
     222             :                   &              (ic**m)*vpw(irec3)*EXP(-ic*&
     223             :                   &              m*phi)*fJJ(iabs(m),irec2)*&
     224           0 :                   &              stars%nstr2(irec2)*mult
     225             :           END IF
     226             :        END DO
     227             :     END DO
     228             : 
     229             :     !-----> preparing arrays for radial grids:selecting from
     230             :     !-----> x & y only the radius
     231             : 
     232             :     nrd = 0
     233             : 
     234           0 :     DO ix = 0,ivfft1 - 1
     235           0 :        iy_loop: DO  iy = 0,ivfft2 - 1
     236           0 :           x = ix*ani1
     237           0 :           IF (x.GT.0.5) x = x - 1.
     238           0 :           y = iy*ani2
     239           0 :           IF (y.GT.0.5) y = y - 1.
     240             :           r = SQRT((x*cell%amat(1,1) + y*cell%amat(1,2))**2 +&
     241           0 :                &               (x*cell%amat(2,1) + y*cell%amat(2,2))**2)
     242           0 :           DO i = 1,nrd
     243           0 :              IF (ABS(r-rr(i)).LE.1.e-6) THEN
     244           0 :                 rmap(ix,iy) = i
     245           0 :                 CYCLE iy_loop
     246             :              END IF
     247             :           END DO
     248           0 :           nrd = nrd + 1
     249           0 :           rmap(ix,iy) = nrd
     250           0 :           rr(nrd) = r
     251             :        ENDDO iy_loop
     252             :     END DO
     253             : 
     254           0 :     DO gzi = -stars%mx3,stars%mx3
     255           0 :        DO m = -MM,MM
     256           0 :           IF (m.LT.0) THEN
     257           0 :              mult = REAL((-1)**m)
     258             :           ELSE
     259             :              mult = 1.
     260             :           END IF
     261           0 :           val_m(gzi,m) = CMPLX(0.,0.)
     262           0 :           irc1 = oneD%ig1(gzi,m)
     263           0 :           IF (irc1.NE.0) THEN
     264           0 :              val_m(gzi,m) = val(irc1)
     265             :           ELSE
     266           0 :              DO irec2 = 1,stars%ng2
     267           0 :                 phi = stars%phi2(irec2)
     268           0 :                 irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),gzi)
     269           0 :                 IF (irec3.NE.0) THEN
     270             :                    val_m(gzi,m) = val_m(gzi,m) +&
     271             :                         &                    (ic**m)*vpw(irec3)*EXP(-ic*&
     272             :                         &                    m*phi)*fJJ(iabs(m),irec2)*&
     273           0 :                         &                    stars%nstr2(irec2)*mult
     274             :                 END IF
     275             :              END DO
     276             :           END IF
     277             :        END DO
     278             :     END DO
     279             : 
     280             :     !-------  During the following staff we miss the m=0,gz=0 -> irec1 = 1
     281             :     !-------  component of the potential which should be also added
     282             :     !-------  in order to get the total potential for the FFT's
     283             : 
     284           0 :     DO ix = 0,ivfft1 - 1
     285           0 :        DO iy = 0,ivfft2 - 1
     286           0 :           irc = rmap(ix,iy)
     287           0 :           r = rr(irc)
     288           0 :           IF (r.GT.cell%z1) THEN
     289           0 :              zf = (r-cell%z1)/vacuum%delz + 1.0
     290           0 :              im = zf
     291           0 :              q = zf - im
     292             :              vis(ix,iy,1) = 0.5*(q-1.)*&
     293             :                   &              (q-2.)*vz(im,1) -&
     294             :                   &              q*(q-2.)*vz(im+1,1) +&
     295           0 :                   &              0.5*q*(q-1.)*vz(im+2,1)
     296             :           ELSE
     297             :              vis(ix,iy,1) = &
     298             :                   &              vz(1,1) - val(1) + tpi_const*&
     299           0 :                   &              psq(1)*(cell%z1*cell%z1 - r*r)/2.
     300             :           END IF
     301           0 :           DO irc1 = 2,nq2_1
     302           0 :              vis(ix,iy,irc1) = CMPLX(0.,0.)
     303             :           END DO
     304             :        END DO
     305             :     END DO
     306             : 
     307           0 :     DO m = 0,MM
     308           0 :        DO k3 = 1,stars%mx3
     309           0 :           CALL modcyli(m,cell%bmat(3,3)*k3*cell%z1,IR(k3,m))
     310             :        END DO
     311             :     END DO
     312             : 
     313           0 :     DO i = 1,nrd
     314           0 :        DO m = 0,MM
     315           0 :           DO k3 = 1,stars%mx3
     316           0 :              CALL modcyli(m,cell%bmat(3,3)*k3*rr(i),IIII(i,k3,m))
     317             :           END DO
     318             :        END DO
     319             :     END DO
     320             : 
     321             :     !------- cycle by positive gz---------->
     322             : 
     323           0 :     DO gzi = 0,stars%mx3                        ! gz
     324             : 
     325             :        !------- cycle by positive m ---------->
     326             : 
     327           0 :        m_loop: DO m = 0,MM
     328             : 
     329             :           !-------------------------------------->
     330             : 
     331           0 :           IF (m.NE.0 .OR. gzi.NE.0) THEN ! m^2 + gz^2.ne.0
     332             : 
     333           0 :              irec1(1) = oneD%ig1(gzi,m)
     334           0 :              irec1(2) = oneD%ig1(gzi,-m)
     335           0 :              irec1(3) = oneD%ig1(-gzi,m)
     336           0 :              irec1(4) = oneD%ig1(-gzi,-m)
     337             : 
     338           0 :              DO l = 1,3
     339           0 :                 DO l1 = l+1,4
     340           0 :                    IF (irec1(l).EQ.irec1(l1)) irec1(l1) = 0
     341             :                 END DO
     342             :              END DO
     343             : 
     344             :              !---> if all the irec1 not equal to zero
     345             : 
     346           0 :              s = 0
     347             : 
     348           0 :              DO l = 1,4
     349           0 :                 s = s + irec1(l)*irec1(l)
     350             :              END DO
     351             : 
     352             :              !---> preparing special functions which depend only on the 
     353             :              !---> absolute value of m and gz
     354             : 
     355           0 :              IF (s.NE.0 .AND. gzi.NE.0) THEN
     356           0 :                 gz = cell%bmat(3,3)*gzi
     357           0 :                 DO  imz = 1,vacuum%nmzxy
     358           0 :                    z = cell%z1 + vacuum%delz*(imz-1)  
     359           0 :                    CALL modcylk(m,gz*z,KK(imz))
     360           0 :                    CALL modcyli(m,gz*z,II(imz))
     361             :                 END DO
     362           0 :                 IIIR = II(1)
     363           0 :                 DO irc = 1,nrd
     364           0 :                    III(irc) = IIII(irc,gzi,m)
     365             :                 END DO
     366           0 :              ELSEIF (s.EQ.0) THEN
     367             :                 CYCLE m_loop
     368           0 :              ELSEIF (s.NE.0 .AND. gzi.EQ.0) THEN  
     369           0 :                 gz = 0.
     370             :              END IF
     371             : 
     372             :              !---> now we start the cycle by +-m,gz
     373             : 
     374           0 :              DO  l = 1,4
     375             : 
     376           0 :                 IF (irec1(l).NE.0) THEN
     377             : 
     378             :                    !--------------------------------------->   
     379             : 
     380           0 :                    DO ix = 0,ivfft1-1
     381           0 :                       DO iy = 0,ivfft2-1
     382           0 :                          vis_help(ix,iy) = CMPLX(0.,0.)
     383             :                       END DO
     384             :                    END DO
     385             : 
     386           0 :                    DO i = 1,vacuum%nmzxy
     387           0 :                       fact(i) = 0.
     388           0 :                       rxy(i) = CMPLX(0.,0.)
     389           0 :                       pvac(i) = CMPLX(0.,0.)
     390           0 :                       pint(i) = CMPLX(0.,0.)
     391             :                    END DO
     392             : 
     393             :                    !-------------- gz = 0 ------------------------------------->
     394             :                    !----------------------------------------------------------->
     395             : 
     396           0 :                    IF (gzi.EQ.0) THEN
     397             : 
     398             :                       !----- this form of the density is just more easy to use
     399             : 
     400           0 :                       DO imz = 1,vacuum%nmzxy
     401           0 :                          rxy(imz) = rhtxy(imz,irec1(l)-1,1)
     402             :                       END DO
     403             : 
     404             :                       !----- vacuum potential caused by the vacuum density
     405             : 
     406             :                       CALL vacp5_0(&
     407             :                            &        vacuum%nmzxyd,vacuum%nmzxy,cell%z1,tpi_const,rxy,m,vacuum%delz,&
     408           0 :                            &        pvac,fact)
     409             : 
     410             :                       !----- vacuum potential caused by the interstitial density
     411             : 
     412           0 :                       aa = CMPLX(0.,0.)
     413           0 :                       m1 = kv2_1(2,irec1(l))
     414           0 :                       gzi1 = kv2_1(1,irec1(l))
     415             : 
     416           0 :                       IF (m1.LT.0) THEN
     417           0 :                          mult = REAL((-1)**m)
     418             :                       ELSE
     419             :                          mult = 1.
     420             :                       END IF
     421             : 
     422           0 :                       DO irec2 = 1,stars%ng2
     423           0 :                          irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),0)
     424           0 :                          IF (irec3.NE.0 .AND. irec2.NE.1) THEN
     425           0 :                             phi = stars%phi2(irec2)
     426           0 :                             g2 = stars%sk2(irec2)
     427             :                             aa = aa +    &
     428             :                                  &              (ic**m1)*psq(irec3)*&
     429             :                                  &              EXP(CMPLX(0.,-m1*phi))*&
     430           0 :                                  &              CMPLX(mult*fJJ(m+1,irec2)/g2,0.)*stars%nstr2(irec2)
     431             :                          END IF
     432             :                       END DO
     433             : 
     434             :                       !----- total vacuum potential
     435             : 
     436           0 :                       pint(:vacuum%nmzxy) =  fact(:vacuum%nmzxy)*aa 
     437             : 
     438           0 :                       vxy(:vacuum%nmzxy,irec1(l)-1,1) = pvac(:vacuum%nmzxy) + pint(:vacuum%nmzxy)
     439             : 
     440             :                       !----- array val further is a boundary values of the
     441             :                       !----- potential V- \tilde \tilde which is created to compensate 
     442             :                       !----- the influence of the outer 'noice' charge density - which 
     443             :                       !----- is just a periodical continuation of the interstitial charge
     444             :                       !----- density, V - \tilde and V - \tilde\tilde are then added in
     445             :                       !----- order to obtain the real interstitial potential           
     446             : 
     447           0 :                       val_help = vxy(1,irec1(l)-1,1) - val(irec1(l))
     448             : 
     449             :                       !----- potential \tilde\tilde{V} is a solution of the Laplase equation
     450             :                       !----- in the interstitial with the boundary conditions val_0 and val_z
     451             :                       !----- further, it is generated on the uniform grid in the unit cell
     452             :                       !----- \tilde{\Omega}, in the space between the cylindrical 
     453             :                       !----- interstitial boundary and the squre boundaries it is put to
     454             :                       !----- the vacuum potential
     455             : 
     456             :                       CALL visp5_0(&
     457             :                            &        vacuum%nmzxyd,vacuum%nmzxy,vacuum%delz,m,ivfft1,ivfft2,l,&
     458             :                            &        rxy,ani1,ani2,cell%z1,cell%amat,&
     459             :                            &        pvac,pint,tpi_const,input%jspins,val_help,&
     460           0 :                            &        vis_help)
     461             : 
     462           0 :                       DO ix = 0,ivfft1 - 1
     463           0 :                          DO iy = 0,ivfft2 - 1
     464           0 :                             vis(ix,iy,irec1(l)) = vis_help(ix,iy)
     465             :                          END DO
     466             :                       END DO
     467             : 
     468             : 
     469             :                       !------- gz.NEQ.0--------------------------------------------->
     470             :                    ELSE   
     471             :                       !------------------------------------------------------------->
     472             : 
     473           0 :                       DO  imz = 1,vacuum%nmzxy
     474           0 :                          rxy(imz) = rhtxy(imz,irec1(l)-1,1)
     475             :                       END DO
     476             : 
     477             :                       !----- vacuum potential caused by the vacuum density        
     478             : 
     479             :                       CALL vacp5_z(&
     480             :                            &        vacuum%nmzxyd,vacuum%nmzxyd,cell%z1,vacuum%delz,fpi_const,II,KK,rxy,m,&
     481           0 :                            &        pvac)
     482             : 
     483             :                       !----- vacuum potential caused by the intst. density
     484             : 
     485           0 :                       a = CMPLX(0.,0.)
     486           0 :                       m1 = kv2_1(2,irec1(l))
     487           0 :                       gzi1 = kv2_1(1,irec1(l))
     488             : 
     489           0 :                       IF (m1.LT.0) THEN
     490           0 :                          mult = REAL((-1)**m)
     491             :                       ELSE
     492             :                          mult = 1.
     493             :                       END IF
     494             : 
     495           0 :                       DO irec2 = 1,stars%ng2
     496           0 :                          irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),gzi1)
     497           0 :                          IF (irec3.NE.0) THEN
     498           0 :                             g = stars%sk3(irec3)
     499           0 :                             g2 = stars%sk2(irec2)
     500           0 :                             phi = stars%phi2(irec2)
     501             :                             b = cell%z1*( gz*iJJ(m+1,gzi)*fJJ(m,irec2) + &
     502           0 :                                  &              g2*II(1)*fJJ(m+1,irec2))/(g*g)
     503             :                             a = a +    &
     504             :                                  &              (ic**m1)*mult*EXP(-ic*m1*phi)&
     505           0 :                                  &              *psq(irec3)*b*stars%nstr2(irec2)
     506             :                          END IF
     507             :                       END DO
     508             : 
     509             :                       !----- total vacuum potential ---------------
     510             : 
     511           0 :                       DO imz = 1,vacuum%nmzxy
     512           0 :                          pint(imz) = fpi_const*a*KK(imz) 
     513           0 :                          vxy(imz,irec1(l)-1,1) =  pint(imz) + pvac(imz)  
     514             :                       END DO
     515             : 
     516           0 :                       val_help = vxy(1,irec1(l)-1,1) - val(irec1(l))
     517             : 
     518             :                       CALL visp5_z(&
     519             :                            &        vacuum%nmzxyd,vacuum%nmzxyd,vacuum%delz,m,ivfft1,ivfft2,IIIR,&
     520             :                            &        rxy,ani1,ani2,cell%z1,cell%amat,pvac,pint,tpi_const,l,&
     521             :                            &        fpi_const,val_help,III,m1,gz,rmap,rr,&
     522           0 :                            &        vis_help)
     523             : 
     524           0 :                       DO ix = 0,ivfft1 - 1
     525           0 :                          DO iy = 0,ivfft2 - 1
     526           0 :                             vis(ix,iy,irec1(l)) = vis_help(ix,iy)
     527             :                          END DO
     528             :                       END DO
     529             : 
     530             :                       !---  end of vacuum contibution to vis ---------------------->
     531             : 
     532             :                    END IF                    ! gz.atoms%neq.0
     533             : 
     534             :                    !---- > finishing the cycle by +-m,gz
     535             : 
     536             :                 END IF
     537             : 
     538             :              ENDDO
     539             :           END IF                    ! m^2+gz^2.atoms%neq.0
     540             : 
     541             :        ENDDO m_loop                  ! m
     542             : 
     543             :     ENDDO                  ! gz 
     544             : 
     545             :     !*************************************************************
     546             :     !-------> finding the Fourier coefficients of the potential
     547             :     !-------  in the interstitial
     548             :     !-------  the scheme (making the potential continuous) is as follows:
     549             :     !-------  1. transforming the \tilde{V} from pw-representation
     550             :     !-------     to the real grid, putting to zero outside the
     551             :     !-------     cylindrical interstitial
     552             :     !-------  2. Adding the \tilde\tilde{V} on the real grid 
     553             :     !-------  3. Transforming back to the pw-representation 
     554             :     !-------  Everything is done in \tilda {\Omega}
     555             : 
     556           0 :     gzmin = -stars%mx3
     557             : 
     558           0 :     IF (sym%zrfs.OR.sym%invs) gzmin = 0
     559             : 
     560           0 :     DO k3 = gzmin,stars%mx3          ! collect the Fourier components
     561             : 
     562           0 :        fxy0 = 0.
     563             : 
     564           0 :        rhti = 0.
     565             : 
     566           0 :        DO irec2 = 1,stars%ng2 - 1
     567           0 :           fxy(irec2) = CMPLX(0.,0.)
     568             :        END DO
     569             : 
     570           0 :        DO irec2 = 1,stars%ng2
     571           0 :           irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
     572           0 :           IF (irec3.NE.0) THEN
     573           0 :              IF (irec2.EQ.1) THEN
     574           0 :                 fxy0 = REAL(vpw(irec3))
     575           0 :                 rhti = AIMAG(vpw(irec3))
     576             :              ELSE
     577           0 :                 fxy(irec2-1) = vpw(irec3)
     578             :              END IF
     579             :           END IF
     580             :        END DO
     581             : 
     582           0 :        af2=0.0
     583           0 :        bf2=0.0
     584             : 
     585             :        CALL fft2d(&
     586             :             &        stars,&
     587             :             &        af2(0),bf2,&
     588             :             &        fxy0,rhti,fxy,&
     589           0 :             &        1,+1)
     590             : 
     591             : 
     592             :        !--> although some harmonics of the vacuum potential are equal to zero
     593             :        !--> the same harmonics of the interstitial potential \Tilde{V} are not,
     594             :        !--> because this potential is determined from pw coefficients of the
     595             :        !--> charge density, which has the periodicity of the rectangular 
     596             :        !--> unit cell, hence, these harmonics should be extracted from the 
     597             :        !--> final interstitial potential
     598             : 
     599           0 :        i = 0
     600             : 
     601           0 :        DO iy = 0,3*stars%mx2-1
     602           0 :           DO ix = 0,3*stars%mx1-1
     603           0 :              x = ix*ani1
     604           0 :              IF (x.GT.0.5) x = x - 1.
     605           0 :              y = iy*ani2
     606           0 :              IF (y.GT.0.5) y = y - 1.
     607             :              r = SQRT((x*cell%amat(1,1) + y*cell%amat(1,2))**2 +&
     608           0 :                   &                  (x*cell%amat(2,1) + y*cell%amat(2,2))**2)               
     609             :              phi = angle(x*cell%amat(1,1) + y*cell%amat(1,2),&
     610           0 :                   &                     x*cell%amat(2,1) + y*cell%amat(2,2))
     611           0 :              vis_tot(ix,iy) = CMPLX(0.,0.)
     612           0 :              j = rmap(ix,iy)
     613           0 :              DO m = -MM,MM
     614           0 :                 irc1 = oneD%ig1(k3,m)
     615           0 :                 IF (irc1.NE.0) THEN
     616             :                    vis_tot(ix,iy) = vis_tot(ix,iy) +&
     617           0 :                         &                    vis(ix,iy,irc1)
     618             :                 ELSE
     619           0 :                    IF (k3.EQ.0) THEN
     620           0 :                       IF (r.LE.cell%z1) THEN
     621             :                          vis_tot(ix,iy) = vis_tot(ix,iy) - &
     622             :                               &                          val_m(k3,m)*&
     623             :                               &                          EXP(CMPLX(0.,m*phi))*&
     624           0 :                               &                          (r**(iabs(m)))/(cell%z1**(iabs(m)))
     625             :                       END IF
     626             :                    ELSE
     627           0 :                       IF (r.LE.cell%z1) THEN
     628             :                          vis_tot(ix,iy) = vis_tot(ix,iy) -&
     629             :                               &                    val_m(k3,m)*IIII(j,iabs(k3),iabs(m))*&
     630           0 :                               &                    EXP(CMPLX(0.,m*phi))/IR(iabs(k3),iabs(m))
     631             :                       END IF
     632             :                    END IF
     633             :                 END IF
     634             :              END DO
     635             : 
     636           0 :              IF (r.LE.cell%z1) THEN
     637           0 :                 af2(i) = af2(i) + REAL(vis_tot(ix,iy))
     638           0 :                 bf2(i) = bf2(i) + AIMAG(vis_tot(ix,iy))
     639             :                 !$$$                 af2(i) = real(vis_tot(ix,iy))
     640             :                 !$$$                 bf2(i) = aimag(vis_tot(ix,iy))
     641             :              ELSE
     642           0 :                 af2(i) = REAL(vis_tot(ix,iy))
     643           0 :                 bf2(i) = AIMAG(vis_tot(ix,iy))
     644             :              END IF
     645           0 :              i = i+1
     646             :           END DO
     647             :        END DO
     648             : 
     649           0 :        rhti = 0.
     650             : 
     651             :        CALL fft2d(&
     652             :             &        stars,&
     653             :             &        af2(0),bf2,&
     654             :             &        gxy0,rhti,gxy,&
     655           0 :             &        1,-1) 
     656             : 
     657           0 :        DO irec2 = 1,stars%ng2
     658           0 :           irec3 = stars%ig(stars%kv2(1,irec2),stars%kv2(2,irec2),k3)
     659           0 :           IF (irec3.NE.0) THEN
     660           0 :              IF (irec2.EQ.1) THEN
     661           0 :                 vpw_help(irec3) = CMPLX(gxy0,rhti)
     662             :              ELSE
     663           0 :                 vpw_help(irec3) = gxy(irec2-1)
     664             :              END IF
     665             :           END IF
     666             :        END DO
     667             : 
     668             :     END DO                    ! gz -> Vpw(.,.,gz)
     669             : 
     670           0 :     DO irec3 = 1,stars%ng3
     671           0 :        vpw(irec3) = vpw_help(irec3)
     672             :        !$$$         vpw(irec3,1) = vpw(irec3,1) + vpw_help(irec3)
     673             :     END DO
     674             : 
     675             : 
     676           0 :     DO irc1 = 2,nq2_1
     677           0 :        DO imz = 1,vacuum%nmzxy
     678           0 :           IF (ABS(vxy(imz,irc1-1,1)).LE.tol_21)&
     679           0 :                &          vxy(imz,irc1-1,1) = CMPLX(0.,0.)
     680             :        END DO
     681             :     END DO
     682             : 
     683             : 
     684             : 
     685           0 :     DEALLOCATE ( KK,II,III,IIII,fact,val,val_m,rxy,&
     686           0 :          &     vis,vis_tot,pvac,pint,vis_help,rmap,rr,&
     687           0 :          &     af2,bf2,vpw_help,fJJ,iJJ )
     688             : 
     689           0 :     RETURN
     690             :   END SUBROUTINE od_vvacis
     691             : END MODULE m_od_vvacis
     692             : 
     693             : 
     694             : 
     695             : 
     696             : 
     697             : 
     698             : 
     699             : 
     700             : 
     701             : 
     702             : 
     703             : 

Generated by: LCOV version 1.13