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 39 40 41 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 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
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 ..
    REAL phase,phasep,q1,zsl1,zsl2,qi,volsli,volintsli
    INTEGER i ,indp,ix1,iy1,iz1,j,n,ns,ind
    COMPLEX x
    !     ..
    !     .. 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) )
     z_z(ind)  = z_z(ind)  + z(j,n)*z(i,n)*phase
             z_z(indp) = z_z(indp) + z(i,n)*z(j,n)*phasep
#else
             z_z(ind) = z_z(ind) +z(j,n)*CONJG(z(i,n))*CMPLX(phase,0.0)     
             z_z(indp)= z_z(indp)+z(i,n)*CONJG(z(j,n))*CMPLX(phasep,0.0)
#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