q_int_sl.F90 3.92 KB
Newer Older
1 2 3
MODULE m_qintsl
  USE m_juDFT
CONTAINS
4
  SUBROUTINE q_int_sl(isp,stars,atoms,sym,cell,ne,lapw,slab,oneD,zMat)          
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
    !     *******************************************************
    !     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
22
    TYPE(t_zMat),INTENT(IN)   :: zMat
23
    TYPE(t_slab),INTENT(INOUT):: slab
24 25
    !
    !     .. Scalar Arguments ..
26
    INTEGER, INTENT (IN) :: ne,isp
27 28
    !     ..
    !     .. Local Scalars ..
29
    REAL q1,zsl1,zsl2,qi,volsli,volintsli
30
    INTEGER i ,indp,ix1,iy1,iz1,j,n,ns,ind
31
    COMPLEX x,phase,phasep
32 33 34
    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: stfunint(:,:),z_z(:)
35

36 37 38 39 40 41 42 43
    !     ..
    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
    !
44
    ALLOCATE ( stfunint(stars%ng3,slab%nsl), z_z(stars%ng3) ) 
45 46 47 48
    !
    !  -----> calculate the integrals of star functions over
    !                     the layer interstitial
    !
49 50 51 52 53
    DO i = 1,slab%nsl
       zsl1 = slab%zsl(1,i)
       zsl2 = slab%zsl(2,i)
       volsli = slab%volsl(i)
       volintsli = slab%volintsl(i)
54 55
       DO j = 1,stars%ng3
          CALL pwint_sl(stars,atoms,sym,zsl1,zsl2,&
56
                        volsli,volintsli,cell,slab%nmtsl(1,i),stars%kv3(1,j),x)
57 58 59 60 61 62 63 64 65
          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
66
       IF (zmat%l_real) THEN
67
          DO  i = 1,lapw%nv(isp)
68
             q1 = q1 + zMat%z_r(i,n)*zMat%z_r(i,n)
69 70 71
          ENDDO
       ELSE
          DO  i = 1,lapw%nv(isp)
72
             q1 = q1 + REAL(zMat%z_c(i,n)*CONJG(zMat%z_c(i,n)))
73 74
          ENDDO
       ENDIF
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
       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)
92
             IF (zmat%l_real) THEN
93 94
                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)
95
             ELSE
96 97
                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
98
             ENDIF
99 100 101 102
          ENDDO
       ENDDO
       ! ----> calculate a charge in the layer interstitial region of the film
       !
103
       DO i = 1,slab%nsl
104 105 106 107
          qi = 0.0
          DO j = 1,stars%ng3
             qi = qi + z_z(j)*stfunint(j,i)
          ENDDO
108
          slab%qintsl(i,n) = qi 
109 110 111 112 113 114 115 116
       ENDDO    ! over vacuum%layers         

    ENDDO ! over states

    DEALLOCATE ( stfunint, z_z ) 

  END SUBROUTINE q_int_sl
END MODULE m_qintsl