od_vvac.f90 3.25 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
MODULE m_od_vvac
CONTAINS
  SUBROUTINE od_vvac(&
       &     stars,vacuum,cell,&
       &     psq,rht,&
       &     vz)

    !     subroutine which calculates the non warped part of the
    !     vacuum potential (m=0,gz=0)
    !                               Y. Mokrousov
    !     the potential in this subroutine can be defined in two
    !     equivalent ways, which nevertheless give a bit defferent 
    !     results, 2nd one seems to be more precise 

    USE m_qsf
    USE m_od_cylbes
    USE m_types
    USE m_constants
    IMPLICIT NONE
    TYPE(t_vacuum),INTENT(IN)   :: vacuum
    TYPE(t_stars),INTENT(IN)    :: stars
    TYPE(t_cell),INTENT(IN)     :: cell


31
    COMPLEX, INTENT (IN) :: psq(stars%ng3)
32 33
    REAL,    INTENT (IN) :: rht(:,:) !(vacuum%nmzd,2)
    REAL,    INTENT (OUT) :: vz(:,:) !(vacuum%nmzd,2)
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

    COMPLEX  rhobar
    INTEGER  k1,k2,irec3,irec2,i,j,ivac,imz,imz1
    REAL     g2 ,a(vacuum%nmzd)
    REAL     fJ,z,zp,phi
    REAL     rht1(vacuum%nmzd)
    REAL     f2(vacuum%nmzd),f22(vacuum%nmzd)

    INTRINSIC cmplx


    DO i = 1,vacuum%nmz
       f2(i) = 0.
       f22(i) = 0.
       DO ivac = 1,vacuum%nvac
49
          vz(i,ivac) = 0.
50 51 52 53 54 55
       END DO
    END DO


    rhobar = -psq(1)

56 57
    DO  k1 = -stars%mx1,stars%mx1
       DO  k2 = -stars%mx2,stars%mx2
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
          irec3 = stars%ig(k1,k2,0)
          IF (irec3.NE.0) THEN
             irec2 = stars%ig2(irec3)
             IF (irec2.NE.1) THEN
                g2 = stars%sk2(irec2)
                phi = stars%phi2(irec2)
                CALL od_cylbes(1,cell%z1*g2,fJ)
                rhobar = rhobar - 2.*psq(irec3)*CMPLX(fJ/(g2*cell%z1),0.0)

             END IF
          END IF
       ENDDO
    ENDDO
    !----> 1st equivalent way      

    DO  i=1,vacuum%nmz
74
       rht1(i) = fpi_const*(cell%z1+(i-1)*vacuum%delz)*rht(i,1)
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
    ENDDO
    CALL qsf(vacuum%delz,rht1(1),f2(1),vacuum%nmz,1)

    DO  i = 1,vacuum%nmz
       f2(i) = tpi_const*cell%z1*cell%z1*rhobar-f2(i)
    ENDDO

    DO  i = 1,vacuum%nmz
       DO  j = 1,vacuum%nmz
          IF (j.LT.i) THEN
             f22(j) = 0.0
          ELSE
             f22(j) = f2(j)/(cell%z1+vacuum%delz*(j-1))
          END IF
       ENDDO
       CALL qsf(vacuum%delz,f22(1),a,vacuum%nmz,0)
       DO  ivac =1,vacuum%nvac
92
          vz(i,ivac) = -a(1)
93 94 95 96 97 98 99 100 101
       ENDDO
    ENDDO
    !----> 2nd equivalent way (via the Green function)

    DO imz = 1,vacuum%nmz
       z = cell%z1 + (imz-1)*vacuum%delz
       DO imz1 = 1,vacuum%nmz
          zp = cell%z1 +  (imz1-1)*vacuum%delz
          IF (imz1.LE.imz) THEN
102
             rht1(imz1) = fpi_const*LOG(z)*zp*rht(imz1,1)
103
          ELSE
104
             rht1(imz1) = fpi_const*LOG(zp)*zp*rht(imz1,1)
105 106 107 108
          END IF

       END DO
       CALL qsf(vacuum%delz,rht1,a,vacuum%nmz,0)
109
       vz(imz,1) = tpi_const*LOG(z)*(cell%z1*cell%z1)*rhobar - a(1)
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131
    END DO

    RETURN
  END SUBROUTINE od_vvac
END MODULE m_od_vvac