q_int_sl.F90 4.31 KB
Newer Older
1 2 3 4
MODULE m_qintsl
  USE m_juDFT
CONTAINS
  SUBROUTINE q_int_sl(isp,stars,atoms,sym, volsl,volintsl, cell,&
5
       ne,lapw, nsl,zsl,nmtsl,oneD, qintslk,zMat,realdata)          
6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
    !     *******************************************************
    !     calculate the charge of the En(k) state 
    !     in the interstitial region of each leyer
    !                                             Yu.M. Koroteev
    !             From pwden_old.F and pwint.F by  c.l.fu
    !     *******************************************************
#include"cpp_double.h"
    USE m_pwintsl
    USE m_types
    IMPLICIT NONE

    TYPE(t_lapw),INTENT(IN)   :: lapw
    TYPE(t_oneD),INTENT(IN)   :: oneD
    TYPE(t_sym),INTENT(IN)    :: sym
    TYPE(t_stars),INTENT(IN)  :: stars
    TYPE(t_cell),INTENT(IN)   :: cell
    TYPE(t_atoms),INTENT(IN)  :: atoms
23
    TYPE(t_zMat),INTENT(IN)   :: zMat
24 25 26 27 28 29 30 31 32
    !
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: ne,isp  ,nsl
    !     ..
    !     .. Array Arguments ..
    INTEGER, INTENT (IN) :: nmtsl(atoms%ntypd,nsl) 
    REAL,    INTENT (IN) :: volintsl(atoms%natd)  
    REAL,    INTENT (IN) :: zsl(2,atoms%natd) ,volsl(atoms%natd) 
    REAL,    INTENT (OUT):: qintslk(:,:)!(nsl,dimension%neigd)
33
    LOGICAL,OPTIONAL, INTENT (IN) :: realdata
34 35
    !     ..
    !     .. Local Scalars ..
36
    REAL q1,zsl1,zsl2,qi,volsli,volintsli
37
    INTEGER i ,indp,ix1,iy1,iz1,j,n,ns,ind
38
    COMPLEX x,phase,phasep
39 40 41
    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: stfunint(:,:),z_z(:)
42 43 44 45 46

    LOGICAL :: l_real
    IF (PRESENT(realdata)) THEN
       l_real=realdata
    ELSE
47
       l_real=zMat%l_real
48
    ENDIF
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
    !     ..
    IF (oneD%odi%d1) CALL juDFT_error("well, does not work with 1D. Not clear how to define a layer.",calledby ="q_int_sl")
    !
    !     calculate the star function expansion coefficients of
    !     the plane wave charge density for each En(k)
    !    
    !     ----> g=0 star
    !
    ALLOCATE ( stfunint(stars%n3d,nsl), z_z(stars%n3d) ) 
    !
    !  -----> calculate the integrals of star functions over
    !                     the layer interstitial
    !
    DO i = 1,nsl
       zsl1 = zsl(1,i)
       zsl2 = zsl(2,i)
       volsli = volsl(i)
       volintsli = volintsl(i)
       DO j = 1,stars%ng3
          CALL pwint_sl(stars,atoms,sym,zsl1,zsl2,&
               volsli,volintsli,cell,nmtsl(1,i),stars%kv3(1,j),x)
          stfunint(j,i) =  x*stars%nstr(j)
       ENDDO  ! over 3D stars
    ENDDO     ! over vacuum%layers
    !
    ! Here, I reordered the stuff to save memory
    !
    DO  n = 1,ne
       z_z(:) = CMPLX(0.0,0.0)
       q1 = 0.0
79 80
       IF (l_real) THEN
          DO  i = 1,lapw%nv(isp)
81
             q1 = q1 + zMat%z_r(i,n)*zMat%z_r(i,n)
82 83 84
          ENDDO
       ELSE
          DO  i = 1,lapw%nv(isp)
85
             q1 = q1 + REAL(zMat%z_c(i,n)*CONJG(zMat%z_c(i,n)))
86 87
          ENDDO
       ENDIF
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
       z_z(1) = q1/cell%omtil
       !
       !     ----> g.ne.0 stars
       !
       DO  i = 1,lapw%nv(isp)
          DO  j = 1,i-1
             ix1 = lapw%k1(j,isp) - lapw%k1(i,isp)
             iy1 = lapw%k2(j,isp) - lapw%k2(i,isp)
             iz1 = lapw%k3(j,isp) - lapw%k3(i,isp)
             IF (iabs(ix1).GT.stars%mx1) CYCLE
             IF (iabs(iy1).GT.stars%mx2) CYCLE
             IF (iabs(iz1).GT.stars%mx3) CYCLE
             ind = stars%ig(ix1,iy1,iz1)
             indp = stars%ig(-ix1,-iy1,-iz1)
             IF (ind.EQ.0 .OR. indp.EQ.0) CYCLE
             phase = stars%rgphs(ix1,iy1,iz1)/ (stars%nstr(ind)*cell%omtil)
             phasep = stars%rgphs(-ix1,-iy1,-iz1)/ (stars%nstr(indp)*cell%omtil)
105
             IF (l_real) THEN
106 107
                z_z(ind)  = z_z(ind)  + zMat%z_r(j,n)*zMat%z_r(i,n)*REAL(phase)
                z_z(indp) = z_z(indp) + zMat%z_r(i,n)*zMat%z_r(j,n)*REAL(phasep)
108
             ELSE
109 110
                z_z(ind) = z_z(ind) +zMat%z_c(j,n)*CONJG(zMat%z_c(i,n))*phase
                z_z(indp)= z_z(indp)+zMat%z_c(i,n)*CONJG(zMat%z_c(j,n))*phasep
111
             ENDIF
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
          ENDDO
       ENDDO
       ! ----> calculate a charge in the layer interstitial region of the film
       !
       DO i = 1,nsl
          qi = 0.0
          DO j = 1,stars%ng3
             qi = qi + z_z(j)*stfunint(j,i)
          ENDDO
          qintslk(i,n) = qi 
       ENDDO    ! over vacuum%layers         

    ENDDO ! over states

    DEALLOCATE ( stfunint, z_z ) 

  END SUBROUTINE q_int_sl
END MODULE m_qintsl