psqpw.F90 8.54 KB
Newer Older
1
module m_psqpw
2 3
  !     ***********************************************************
  !     generates the fourier coefficients of pseudo charge density
4 5 6 7 8 9 10 11 12 13 14
  !     
  !     For yukawa_residual = .true. the subroutines calculate the
  !     pseudo charge density for the generation of the Yukawa 
  !     potential instead of the Coulomb potential. This is used in 
  !     the preconditioning of the SCF iteration for metallic systems.
  !      
  !     references:
  !     for both the Coulomb and Yukawa cases:
  !     F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
  !     or see the original paper for the normal Coulomb case only:
  !     M. Weinert: J. Math. Phys. 22(11) (1981) p.2434 eq. (10)-(15)
15
  !     ***********************************************************
16

17
contains
18

19
  subroutine psqpw( mpi, atoms, sphhar, stars, vacuum, dimension, cell, input, sym, oneD, &
20
       &     qpw, rho, rht, l_xyav, potdenType, psq )
21

22
#include"cpp_double.h"
23 24 25 26 27 28 29 30
    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
31 32
    use m_DoubleFactorial
    use m_SphBessel
33
    implicit none
34

35 36 37 38 39 40 41 42 43 44 45 46 47 48
    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
    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
    integer,            intent(in)  :: potdenType
50
    complex,            intent(out) :: psq(stars%ng3)
51

52 53 54 55 56 57
    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
    complex                         :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
    complex                         :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
    real                            :: q2(vacuum%nmzd)
58 59
    real                            :: pn(0:atoms%lmaxd,atoms%ntype)
    real                            :: aj(0:atoms%lmaxd+DIMENSION%ncvd+1)
60 61 62
    real                            :: rht1(vacuum%nmz)
    real, allocatable, dimension(:) :: il, kl
    real                            :: g0(atoms%ntype)
63
#ifdef CPP_MPI
64 65 66
    include 'mpif.h'
    integer                         :: ierr(3)
    complex, allocatable            :: c_b(:)
67 68 69
#endif

    ! Calculate multipole moments
70
    call mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, potdenType, qlm )
71
#ifdef CPP_MPI
72 73 74 75
    psq(:) = cmplx( 0.0, 0.0 )
    call MPI_BCAST( qpw, size(qpw), CPP_MPI_COMPLEX, 0, mpi%mpi_comm, ierr )
    nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
    call MPI_BCAST( qlm, nd, CPP_MPI_COMPLEX, 0, mpi%MPI_COMM, ierr )
76
#endif
77 78 79 80 81 82 83 84

    ! prefactor in (A10) (Coulomb case) or (A11) (Yukawa case)
    ! nc(n) is the integer p in the paper; ncv(n) is l + p
    ! Coulomb case: pn(l,n) = (2 * l + 2 * p + 3)!! / ( (2 * l + 1)!! * R ** (ncv(n) + 1) ), 
    ! Yukawa case: pn(l,n) = lambda ** (l + p + 1) / ( i_{l+p+1}(lambda * R) * (2 * l + 1)!! )
    ! g0 is the prefactor for the q=0 component in (A13)
    pn = 0.
    do n = 1, atoms%ntype
85
      if ( potdenType /= POTDEN_TYPE_POTYUK ) then
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
        do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
          pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
        end do
      else
        allocate( il(0:atoms%ncv(n)+1), kl(0:atoms%ncv(n)+1) )
        call ModSphBessel( il(0:), kl(0:), input%preconditioning_param * atoms%rmt(n), atoms%ncv(n) + 1 )
        g0(n) = ( input%preconditioning_param * atoms%rmt(n) ) ** ( atoms%ncv(n) + 1 ) / DoubleFactorial( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) !p=ncv(n)
        do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
          pn(l, n) = input%preconditioning_param ** ( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) / DoubleFactorial( l )
        end do
        deallocate( il, kl )
      end if
    end do

    ! q=0 term: see (A12) (Coulomb case) or (A13) (Yukawa case)
101
    if( mpi%irank == 0 ) then
102
    s = 0.
103
    do n = 1, atoms%ntype
104
      if ( potdenType /= POTDEN_TYPE_POTYUK ) then
105 106 107 108 109
        s = s + atoms%neq(n) * real( qlm(0,0,n) )
      else
        s = s + atoms%neq(n) * real( qlm(0,0,n) ) * g0(n)
      end if
    end do
110
    !if( mpi%irank == 0 ) then
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
      psq(1) = qpw(1) + ( sfp_const / cell%omtil ) * s
    end if

    ! q/=0 term: see (A10) (Coulomb case) or (A11) (Yukawa case)
    fpo = 1. / cell%omtil
    !$omp parallel do default( shared ) private( pylm, sa, n, ncvn, aj, 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
        call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, atoms%neq, &
             atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm )
      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 >= 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) )
          end do
60        sl = sl + pn(l,n) / ( stars%sk3(k) ** n1 ) * aj( ncvn + 1 ) * sm
        end do
        sa = sa + atoms%neq(n) * sl
      end do
      psq(k) = qpw(k) + fpo * sa
    end do
    !$omp end parallel do
145
#ifdef CPP_MPI
146 147 148 149 150 151
    allocate( c_b(stars%ng3) )
    call MPI_REDUCE( psq, c_b, stars%ng3, CPP_MPI_COMPLEX, MPI_SUM, 0, mpi%MPI_COMM, ierr )
    if ( mpi%irank == 0 ) then
       psq(:stars%ng3) = c_b(:stars%ng3)
    end if
    deallocate( c_b )
152
#endif
153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183

    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) == 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
        end do
        psint = cell%area * psint
      else if ( input%film .and. oneD%odi%d1 ) then
        psint = (0.0, 0.0)
        do k = 2, stars%ng3
          if ( stars%kv3(3,k) == 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
        end do
        psint = psint + psq(1) * stars%nstr(1) * cell%vol
      else if ( .not. input%film ) then
        psint = psq(1) * stars%nstr(1) * cell%omtil
      end if
      write( 6, fmt=8000 ) psint
      write( 16, fmt=8000 ) psint
184
8000  format (/,10x,'integral of pseudo charge density inside the slab='&
185
            &       ,5x,2f11.6)
186
      if ( .not. input%film .or. potdenType == POTDEN_TYPE_POTYUK ) return
187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209

      ! Normalized pseudo density
      if ( .not. oneD%odi%d1 ) then
        qvac = 0.0
        do ivac = 1, vacuum%nvac
          call qsf( vacuum%delz, rht(1,ivac), q2, vacuum%nmz, 0 )
          q2(1) = q2(1) * cell%area
          qvac = qvac + q2(1) * 2. / real( vacuum%nvac )
        end do
        qvac = qvac - 2 * input%sigma
      else
        qvac = 0.0
        do nz = 1, vacuum%nmz
          rht1(nz) = ( cell%z1 + ( nz - 1 ) * vacuum%delz ) * rht(nz,vacuum%nvac)
        end do
        call qsf( vacuum%delz, rht1(1), q2, vacuum%nmz, 0 )
        qvac = cell%area * q2(1)
      end if
      if ( l_xyav ) return
      fact = ( qvac + psint ) / ( stars%nstr(1) * cell%vol )
      psq(1) = psq(1) - fact
      write( 6, fmt=8010 ) fact * 1000
      write( 16, fmt=8010 ) fact * 1000
210
8010  format (/,10x,'                     1000 * normalization const. ='&
211
            &       ,5x,2f11.6)
212
    end if ! mpi%irank == 0 
213

214
  end subroutine psqpw
215

216
end module m_psqpw