q_int_sl.F90 4.15 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)          
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
    !
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: ne,isp  ,nsl
    !     ..
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
29 30 31
    INTEGER, INTENT (IN) :: nmtsl(atoms%ntype,nsl) 
    REAL,    INTENT (IN) :: volintsl(atoms%nat)  
    REAL,    INTENT (IN) :: zsl(2,atoms%nat) ,volsl(atoms%nat) 
32 33 34
    REAL,    INTENT (OUT):: qintslk(:,:)!(nsl,dimension%neigd)
    !     ..
    !     .. Local Scalars ..
35
    REAL q1,zsl1,zsl2,qi,volsli,volintsli
36
    INTEGER i ,indp,ix1,iy1,iz1,j,n,ns,ind
37
    COMPLEX x,phase,phasep
38 39 40
    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: stfunint(:,:),z_z(:)
41

42 43 44 45 46 47 48 49
    !     ..
    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
    !
50
    ALLOCATE ( stfunint(stars%ng3,nsl), z_z(stars%ng3) ) 
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
    !
    !  -----> 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
72
       IF (zmat%l_real) THEN
73
          DO  i = 1,lapw%nv(isp)
74
             q1 = q1 + zMat%z_r(i,n)*zMat%z_r(i,n)
75 76 77
          ENDDO
       ELSE
          DO  i = 1,lapw%nv(isp)
78
             q1 = q1 + REAL(zMat%z_c(i,n)*CONJG(zMat%z_c(i,n)))
79 80
          ENDDO
       ENDIF
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
       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)
98
             IF (zmat%l_real) THEN
99 100
                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)
101
             ELSE
102 103
                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
104
             ENDIF
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
          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