q_int_sl.F90 4.27 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
MODULE m_qintsl
  USE m_juDFT
CONTAINS
  SUBROUTINE q_int_sl(isp,stars,atoms,sym, volsl,volintsl, cell,&
       z,ne,lapw, nsl,zsl,nmtsl,oneD, qintslk)          
    !     *******************************************************
    !     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
    !
    !     .. 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)
#if ( defined(CPP_INVERSION) && !defined(CPP_SOC) )
     REAL,    INTENT (IN) :: z(:,:)!(dimension%nbasfcn,dimension%neigd)
#else
    COMPLEX, INTENT (IN) :: z(:,:)
#endif
    !     ..
    !     .. Local Scalars ..
39
    REAL q1,zsl1,zsl2,qi,volsli,volintsli
40
    INTEGER i ,indp,ix1,iy1,iz1,j,n,ns,ind
41
    COMPLEX x,phase,phasep
42 43 44 45 46 47 48 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 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104
    !     ..
    !     .. Local Arrays ..
    COMPLEX, ALLOCATABLE :: stfunint(:,:),z_z(:)
#if ( defined(CPP_INVERSION) && !defined(CPP_SOC) )
     REAL,    ALLOCATABLE :: z_h(:,:)
#else
    COMPLEX, ALLOCATABLE :: z_h(:,:)
#endif
    !     ..
    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
       DO  i = 1,lapw%nv(isp)
#if ( defined(CPP_INVERSION) && !defined(CPP_SOC) )
     q1 = q1 + z(i,n)*z(i,n)
#else
          q1 = q1 + REAL(z(i,n)*CONJG(z(i,n)))
#endif
       ENDDO
       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)
#if ( defined(CPP_INVERSION) && !defined(CPP_SOC) )
105 106
             z_z(ind)  = z_z(ind)  + z(j,n)*z(i,n)*REAL(phase)
             z_z(indp) = z_z(indp) + z(i,n)*z(j,n)*REAL(phasep)
107
#else
108 109
             z_z(ind) = z_z(ind) +z(j,n)*CONJG(z(i,n))*phase     
             z_z(indp)= z_z(indp)+z(i,n)*CONJG(z(j,n))*phasep
110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
#endif
          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