psqpw.F90 7.53 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
MODULE m_psqpw
  !     ***********************************************************
  !     generates the fourier coefficients of pseudo charge density
  !                                   c.l.fu
  !         corrected april 1990   m.w.
  !
  ! cf. M.Weinert J.Math.Phys. 22(11) (1981) p.2434 eq. (10)-(15)
  !
  !
  !     parallelized 04/08 gb
  !     ***********************************************************
CONTAINS
  SUBROUTINE psqpw(mpi,&
       &                 atoms,sphhar,stars,vacuum,&
       &                 DIMENSION,cell,input,sym,oneD,&
       &                 qpw,rho,rht,l_xyav,&
       &                 psq)
    !
#include"cpp_double.h"
    USE m_constants
    USE m_phasy1
    USE m_mpmom 
    USE m_sphbes
    USE m_qsf
    USE m_od_phasy
    USE m_od_cylbes
    USE m_types

    IMPLICIT NONE
    !     ..
    !     .. Scalar Arguments ..
    TYPE(t_mpi),INTENT(IN)       :: mpi
    TYPE(t_atoms),INTENT(IN)     :: atoms
    TYPE(t_sphhar),INTENT(IN)    :: sphhar
    TYPE(t_stars),INTENT(IN)     :: stars
    TYPE(t_vacuum),INTENT(IN)    :: vacuum
    TYPE(t_dimension),INTENT(IN) :: DIMENSION
    TYPE(t_cell),INTENT(IN)      :: cell
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_sym),INTENT(IN)       :: sym
    TYPE(t_oneD),INTENT(IN)      :: oneD

    LOGICAL, INTENT (IN) :: l_xyav
    !     ..
    !     .. Array Arguments ..
46 47 48
    COMPLEX, INTENT (IN) :: qpw(stars%ng3) 
    REAL,    INTENT (IN) :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) 
    REAL,    INTENT (IN) :: rht(vacuum%nmzd,2)
49
    COMPLEX, INTENT (OUT):: psq(stars%ng3)
50 51 52 53 54 55 56 57 58
    !-odim
    !+odim
    !     ..
    !     .. Local Scalars ..
    COMPLEX psint,sa,sl,sm
    REAL f,fact,fpo,gz,p,qvac,rmtl,s,fJ,gr,g
    INTEGER ivac,k,l ,n,n1,nc,ncvn,lm,ll1 ,nd,m,nz
    !     ..
    !     .. Local Arrays ..
59 60 61
    COMPLEX pylm( (atoms%lmaxd+1)**2 ,atoms%ntype )
    COMPLEX  qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
    REAL q2(vacuum%nmzd),pn(0:atoms%lmaxd,atoms%ntype),aj(0:atoms%lmaxd+DIMENSION%ncvd+1)
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
    REAL rht1(vacuum%nmz)
#ifdef CPP_MPI
    INCLUDE 'mpif.h'
    INTEGER ierr(3)
    COMPLEX, ALLOCATABLE :: c_b(:)
#endif
    !     ..


    !
    ! Calculate multipole moments
    !
    CALL mpmom(mpi,&
         &           atoms,sphhar,stars,&
         &           sym,cell,oneD,&
         &           qpw,rho,&
         &           qlm)
#ifdef CPP_MPI
    psq(:) = CMPLX(0.0,0.0)
81
    CALL MPI_BCAST(qpw,size(qpw),CPP_MPI_COMPLEX,0,mpi%mpi_comm,ierr)
82
    nd = (2*atoms%lmaxd+1)*(atoms%lmaxd+1)*atoms%ntype
83
    CALL MPI_BCAST(qlm,nd,CPP_MPI_COMPLEX,0,mpi%MPI_COMM,ierr)
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
#endif
    !
    ! pn(l,n) = (2l + 2nc(n) + 3)!! / (2l + 1)!! R^l  ;   ncv(n)=n+l in paper
    ! 
    DO n = 1,atoms%ntype
       rmtl = 1.0
       DO l = 0,atoms%lmax(n)
          IF (l.GE.atoms%ncv(n)) THEN
             pn(l,n) = 0.0
          ELSE
             p = 1.
             DO nc = l,atoms%ncv(n)
                p = p* (2*nc+3)
             ENDDO
             pn(l,n) = p/rmtl
          END IF
          rmtl = rmtl*atoms%rmt(n)
       ENDDO
    ENDDO
    !
    ! G eq 0 term (eq.29) : \tilde \rho_s (0) = \sqrt{4 pi} / \Omega \sum_i \tilde q_{00}^i
    !
    s = 0.
    DO n = 1,atoms%ntype
       s = s + atoms%neq(n)*REAL(qlm(0,0,n))
    ENDDO
    IF (mpi%irank == 0) THEN
111
       psq(1) = qpw(1) + (sfp_const/cell%omtil)*s
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
    ENDIF
    !
    ! G ne 0 term (eq.28) : \tilde \rho_s (K) = 4 pi / \Omega \sum_{lmi} (-i)^l \exp{-iK\xi_i}
    !                    (2n+3)!!/(2l+1)!! * 1/R_i^l * j_{n+1}(KR_i)/(KR_i)^{n-l+1} Y_{lm} (K)
    !
    fpo = 1./cell%omtil
    !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(pylm,sa,n,ncvn,aj,&
    !$OMP&  sl,l,n1,ll1,sm,m,lm)
    DO k = mpi%irank+2, stars%ng3, mpi%isize
       IF (.NOT.oneD%odi%d1) THEN
          CALL phasy1(&
               &               atoms,stars,sym,&
               &               cell,k,&
               &               pylm)
       ELSE
          !-odim
          CALL od_phasy(&
129
               &        atoms%ntype,stars%ng3,atoms%nat,atoms%lmaxd,atoms%ntype,atoms%neq,atoms%lmax,&
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152
               &        atoms%taual,cell%bmat,stars%kv3,k,oneD%odi,oneD%ods,&
               &        pylm)
          !+odim
       END IF
       !
       sa = 0.
       DO n = 1,atoms%ntype
          ncvn = atoms%ncv(n)
          CALL sphbes(ncvn+1,stars%sk3(k)*atoms%rmt(n),aj)
          sl = 0.
          DO l = 0,atoms%lmax(n)
             IF (l.GE.ncvn) GO TO 60
             n1 = ncvn - l + 1
             ll1 = l*(l+1) + 1
             sm = 0.
             DO m = -l,l
                lm = ll1 + m 
                sm = sm + qlm(m,l,n)*CONJG(pylm(lm,n))
             ENDDO
60           sl = sl + pn(l,n)/ ((stars%sk3(k)*atoms%rmt(n))**n1)*aj(ncvn+1)*sm
          ENDDO
          sa = sa + atoms%neq(n)*sl
       ENDDO
153
       psq(k) = qpw(k) + fpo*sa
154 155 156
    ENDDO
    !$OMP END PARALLEL DO
#ifdef CPP_MPI
157 158
    ALLOCATE(c_b(stars%ng3))
    CALL MPI_REDUCE(psq,c_b,stars%ng3,CPP_MPI_COMPLEX,MPI_SUM,0,mpi%MPI_COMM,ierr)
159
    IF (mpi%irank.EQ.0) THEN
160
       psq(:stars%ng3)=c_b(:stars%ng3)
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
    ENDIF
    DEALLOCATE (c_b)
#endif
    IF (mpi%irank == 0) THEN
       !
       ! Check: integral of the pseudo charge density within the slab
       !
       IF (input%film .AND. .NOT.oneD%odi%d1) THEN
          psint = psq(1)*stars%nstr(1)*vacuum%dvac
          DO k = 2,stars%ng3
             IF (stars%ig2(k).EQ.1) THEN
                gz = stars%kv3(3,k)*cell%bmat(3,3)
                f = 2.*SIN(gz*cell%z1)/gz
                psint = psint + stars%nstr(k)*psq(k)*f
             END IF
          ENDDO
          psint = cell%area*psint
       ELSEIF (input%film .AND. oneD%odi%d1) THEN
          !-odim
          psint = (0.0,0.0)
          DO k = 2,stars%ng3
             IF (stars%kv3(3,k).EQ.0) THEN
                g = (stars%kv3(1,k)*cell%bmat(1,1) + stars%kv3(2,k)*cell%bmat(2,1))**2 +&
                     &             (stars%kv3(1,k)*cell%bmat(1,2) + stars%kv3(2,k)*cell%bmat(2,2))**2
                gr = SQRT(g)
                CALL od_cylbes(1,gr*cell%z1,fJ)
                f = 2*cell%vol*fJ/(gr*cell%z1)
                psint = psint + stars%nstr(k)*psq(k)*f
             END IF
          ENDDO
          psint = psint + psq(1)*stars%nstr(1)*cell%vol
          !+odim
       ELSEIF (.NOT.input%film) THEN
          psint = psq(1)*stars%nstr(1)*cell%omtil
       ENDIF
       WRITE (6,FMT=8000) psint
       WRITE (16,FMT=8000) psint
8000   FORMAT (/,10x,'integral of pseudo charge density inside the slab='&
            &       ,5x,2f11.6)
       IF (.NOT.input%film) RETURN
       !
       ! Normalized pseudo density
       !
       IF (.NOT.oneD%odi%d1) THEN
          qvac = 0.0
          DO ivac = 1,vacuum%nvac
207
             CALL qsf(vacuum%delz,rht(1,ivac),q2,vacuum%nmz,0)
208 209 210
             q2(1) = q2(1)*cell%area
             qvac = qvac + q2(1)*2./REAL(vacuum%nvac)
          ENDDO
211
          qvac = qvac - 2*input%sigma
212 213 214 215 216
       ELSE
          !-odim
          qvac = 0.0
          DO nz = 1,vacuum%nmz
             rht1(nz) = (cell%z1+(nz-1)*vacuum%delz)*&
217
                  &           rht(nz,vacuum%nvac)
218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238
          ENDDO
          CALL qsf(vacuum%delz,rht1(1),q2,vacuum%nmz,0)
          qvac = cell%area*q2(1)
          !+odim
       END IF
       !      fact = abs(qvac/psint)
       !      DO k = 1,nq3
       !         psq(k) = fact*psq(k)
       !      ENDDO
       IF (l_xyav) RETURN
       fact = (qvac + psint)/(stars%nstr(1)*cell%vol)
       psq(1) = psq(1) - fact
       !-gu
       WRITE (6,FMT=8010) fact*1000
       WRITE (16,FMT=8010) fact*1000
8010   FORMAT (/,10x,'                     1000 * normalization const. ='&
            &       ,5x,2f11.6)
       !
    ENDIF ! mpi%irank == 0 
  END SUBROUTINE psqpw
END MODULE m_psqpw