LCOV - code coverage report
Current view: top level - vgen - psqpw.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 85 100 85.0 %
Date: 2019-09-08 04:53:50 Functions: 2 2 100.0 %

          Line data    Source code
       1             : module m_psqpw
       2             :   !     ***********************************************************
       3             :   !     generates the fourier coefficients of pseudo charge density
       4             :   !     
       5             :   !     For yukawa_residual = .true. the subroutines calculate the
       6             :   !     pseudo charge density for the generation of the Yukawa 
       7             :   !     potential instead of the Coulomb potential. This is used in 
       8             :   !     the preconditioning of the SCF iteration for metallic systems.
       9             :   !      
      10             :   !     references:
      11             :   !     for both the Coulomb and Yukawa cases:
      12             :   !     F. Tran, P. Blaha: Phys. Rev. B 83, 235118 (2011)
      13             :   !     or see the original paper for the normal Coulomb case only:
      14             :   !     M. Weinert: J. Math. Phys. 22(11) (1981) p.2434 eq. (10)-(15)
      15             :   !     ***********************************************************
      16             : 
      17             : contains
      18             : 
      19         346 :   subroutine psqpw( mpi, atoms, sphhar, stars, vacuum,  cell, input, sym, oneD, &
      20         346 :        &     qpw, rho, rht, l_xyav, potdenType, psq )
      21             : 
      22             : #include"cpp_double.h"
      23             :     use m_constants
      24             :     use m_phasy1
      25             :     use m_mpmom 
      26             :     use m_sphbes
      27             :     use m_qsf
      28             :     use m_od_phasy
      29             :     use m_od_cylbes
      30             :     use m_types
      31             :     use m_DoubleFactorial
      32             :     use m_SphBessel
      33             :     implicit none
      34             : 
      35             :     type(t_mpi),        intent(in)  :: mpi
      36             :     type(t_atoms),      intent(in)  :: atoms
      37             :     type(t_sphhar),     intent(in)  :: sphhar
      38             :     type(t_stars),      intent(in)  :: stars
      39             :     type(t_vacuum),     intent(in)  :: vacuum
      40             :     type(t_cell),       intent(in)  :: cell
      41             :     type(t_input),      intent(in)  :: input
      42             :     type(t_sym),        intent(in)  :: sym
      43             :     type(t_oneD),       intent(in)  :: oneD
      44             :     logical,            intent(in)  :: l_xyav
      45             :     complex,            intent(in)  :: qpw(stars%ng3) 
      46             :     real,               intent(in)  :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) 
      47             :     real,               intent(in)  :: rht(vacuum%nmzd,2)
      48             :     integer,            intent(in)  :: potdenType
      49             :     complex,            intent(out) :: psq(stars%ng3)
      50             : 
      51             :     complex                         :: psint, sa, sl, sm
      52             :     real                            :: f, fact, fpo, gz, p, qvac, rmtl, s, fJ, gr, g
      53             :     integer                         :: ivac, k, l, n, n1, nc, ncvn, lm, ll1, nd, m, nz
      54         346 :     complex                         :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
      55         692 :     complex                         :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
      56         692 :     real                            :: q2(vacuum%nmzd)
      57         692 :     real                            :: pn(0:atoms%lmaxd,atoms%ntype)
      58         346 :     real                            :: aj(0:atoms%lmaxd+maxval(atoms%ncv)+1)
      59         692 :     real                            :: rht1(vacuum%nmz)
      60         346 :     real, allocatable, dimension(:) :: il, kl
      61         692 :     real                            :: g0(atoms%ntype)
      62             : #ifdef CPP_MPI
      63             :     include 'mpif.h'
      64             :     integer                         :: ierr(3)
      65         346 :     complex, allocatable            :: c_b(:)
      66             : #endif
      67             : 
      68             :     ! Calculate multipole moments
      69         346 :     call mpmom( input, mpi, atoms, sphhar, stars, sym, cell, oneD, qpw, rho, potdenType, qlm )
      70             : #ifdef CPP_MPI
      71      348950 :     psq(:) = cmplx( 0.0, 0.0 )
      72         346 :     call MPI_BCAST( qpw, size(qpw), CPP_MPI_COMPLEX, 0, mpi%mpi_comm, ierr )
      73         346 :     nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
      74         346 :     call MPI_BCAST( qlm, nd, CPP_MPI_COMPLEX, 0, mpi%MPI_COMM, ierr )
      75             : #endif
      76             : 
      77             :     ! prefactor in (A10) (Coulomb case) or (A11) (Yukawa case)
      78             :     ! nc(n) is the integer p in the paper; ncv(n) is l + p
      79             :     ! Coulomb case: pn(l,n) = (2 * l + 2 * p + 3)!! / ( (2 * l + 1)!! * R ** (ncv(n) + 1) ), 
      80             :     ! Yukawa case: pn(l,n) = lambda ** (l + p + 1) / ( i_{l+p+1}(lambda * R) * (2 * l + 1)!! )
      81             :     ! g0 is the prefactor for the q=0 component in (A13)
      82        1034 :     pn = 0.
      83        1722 :     do n = 1, atoms%ntype
      84        1034 :       if ( potdenType /= POTDEN_TYPE_POTYUK ) then
      85        6736 :         do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
      86         658 :           pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
      87             :         end do
      88             :       else
      89          30 :         allocate( il(0:atoms%ncv(n)+1), kl(0:atoms%ncv(n)+1) )
      90          30 :         call ModSphBessel( il(0:), kl(0:), input%preconditioning_param * atoms%rmt(n), atoms%ncv(n) + 1 )
      91          30 :         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)
      92         300 :         do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
      93         300 :           pn(l, n) = input%preconditioning_param ** ( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) / DoubleFactorial( l )
      94             :         end do
      95          30 :         deallocate( il, kl )
      96             :       end if
      97             :     end do
      98             : 
      99             :     ! q=0 term: see (A12) (Coulomb case) or (A13) (Yukawa case)
     100         346 :     if( mpi%irank == 0 ) then
     101         173 :     s = 0.
     102         517 :     do n = 1, atoms%ntype
     103         517 :       if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     104         329 :         s = s + atoms%neq(n) * real( qlm(0,0,n) )
     105             :       else
     106          15 :         s = s + atoms%neq(n) * real( qlm(0,0,n) ) * g0(n)
     107             :       end if
     108             :     end do
     109             :     !if( mpi%irank == 0 ) then
     110         173 :       psq(1) = qpw(1) + ( sfp_const / cell%omtil ) * s
     111             :     end if
     112             : 
     113             :     ! q/=0 term: see (A10) (Coulomb case) or (A11) (Yukawa case)
     114         346 :     fpo = 1. / cell%omtil
     115      174475 :     !$omp parallel do default( shared ) private( pylm, sa, n, ncvn, aj, sl, l, n1, ll1, sm, m, lm )
     116             :     do k = mpi%irank+2, stars%ng3, mpi%isize
     117      174129 :       if ( .not. oneD%odi%d1 ) then
     118      174129 :         call phasy1( atoms, stars, sym, cell, k, pylm )
     119             :       else
     120             :         call od_phasy( atoms%ntype, stars%ng3, atoms%nat, atoms%lmaxd, atoms%ntype, atoms%neq, &
     121           0 :              atoms%lmax, atoms%taual, cell%bmat, stars%kv3, k, oneD%odi, oneD%ods, pylm )
     122             :       end if
     123      174129 :       sa = 0.
     124      463910 :       do n = 1, atoms%ntype
     125      289781 :         ncvn = atoms%ncv(n)
     126      289781 :         call sphbes( ncvn + 1 , stars%sk3(k) * atoms%rmt(n), aj )
     127      289781 :         sl = 0.
     128     2969192 :         do l = 0, atoms%lmax(n)
     129     2679411 :           if ( l >= ncvn ) go to 60
     130     2679411 :           n1 = ncvn - l + 1
     131     2679411 :           ll1 = l * ( l + 1 ) + 1
     132     2679411 :           sm = 0.
     133    27860592 :           do m = -l, l
     134    25181181 :             lm = ll1 + m 
     135    27860592 :             sm = sm + qlm(m,l,n) * conjg( pylm(lm,n) )
     136             :           end do
     137     2969192 : 60        sl = sl + pn(l,n) / ( stars%sk3(k) ** n1 ) * aj( ncvn + 1 ) * sm
     138             :         end do
     139      463910 :         sa = sa + atoms%neq(n) * sl
     140             :       end do
     141      174821 :       psq(k) = qpw(k) + fpo * sa
     142             :     end do
     143             :     !$omp end parallel do
     144             : #ifdef CPP_MPI
     145         346 :     allocate( c_b(stars%ng3) )
     146         346 :     call MPI_REDUCE( psq, c_b, stars%ng3, CPP_MPI_COMPLEX, MPI_SUM, 0, mpi%MPI_COMM, ierr )
     147         346 :     if ( mpi%irank == 0 ) then
     148      174475 :        psq(:stars%ng3) = c_b(:stars%ng3)
     149             :     end if
     150         346 :     deallocate( c_b )
     151             : #endif
     152             : 
     153         346 :     if ( mpi%irank == 0 ) then
     154         173 :       if ( potdenType == POTDEN_TYPE_POTYUK ) return
     155             :       ! Check: integral of the pseudo charge density within the slab
     156         170 :       if ( input%film .and. .not. oneD%odi%d1 ) then
     157           5 :         psint = psq(1) * stars%nstr(1) * vacuum%dvac
     158       32588 :         do k = 2, stars%ng3
     159       16294 :           if ( stars%ig2(k) == 1 ) then
     160         125 :             gz = stars%kv3(3,k) * cell%bmat(3,3)
     161         125 :             f = 2. * sin( gz * cell%z1 ) / gz
     162         125 :             psint = psint + stars%nstr(k) * psq(k) * f
     163             :           end if
     164             :         end do
     165           5 :         psint = cell%area * psint
     166         165 :       else if ( input%film .and. oneD%odi%d1 ) then
     167           0 :         psint = (0.0, 0.0)
     168           0 :         do k = 2, stars%ng3
     169           0 :           if ( stars%kv3(3,k) == 0 ) then
     170             :             g = ( stars%kv3(1,k) * cell%bmat(1,1) + stars%kv3(2,k) * cell%bmat(2,1) ) ** 2 + &
     171           0 :               & ( stars%kv3(1,k) * cell%bmat(1,2) + stars%kv3(2,k) * cell%bmat(2,2) ) ** 2
     172           0 :             gr = sqrt( g )
     173           0 :             call od_cylbes( 1, gr * cell%z1, fJ )
     174           0 :             f = 2 * cell%vol * fJ / ( gr * cell%z1 )
     175           0 :             psint = psint + stars%nstr(k) * psq(k) * f
     176             :           end if
     177             :         end do
     178           0 :         psint = psint + psq(1) * stars%nstr(1) * cell%vol
     179         165 :       else if ( .not. input%film ) then
     180         165 :         psint = psq(1) * stars%nstr(1) * cell%omtil
     181             :       end if
     182         170 :       write( 6, fmt=8000 ) psint
     183             : 8000  format (/,10x,'integral of pseudo charge density inside the slab='&
     184             :             &       ,5x,2f11.6)
     185         170 :       if ( .not. input%film .or. potdenType == POTDEN_TYPE_POTYUK ) return
     186             : 
     187             :       ! Normalized pseudo density
     188           5 :       if ( .not. oneD%odi%d1 ) then
     189           5 :         qvac = 0.0
     190          10 :         do ivac = 1, vacuum%nvac
     191           5 :           call qsf( vacuum%delz, rht(1,ivac), q2, vacuum%nmz, 0 )
     192           5 :           q2(1) = q2(1) * cell%area
     193          10 :           qvac = qvac + q2(1) * 2. / real( vacuum%nvac )
     194             :         end do
     195           5 :         qvac = qvac - 2 * input%sigma
     196             :       else
     197           0 :         qvac = 0.0
     198           0 :         do nz = 1, vacuum%nmz
     199           0 :           rht1(nz) = ( cell%z1 + ( nz - 1 ) * vacuum%delz ) * rht(nz,vacuum%nvac)
     200             :         end do
     201           0 :         call qsf( vacuum%delz, rht1(1), q2, vacuum%nmz, 0 )
     202           0 :         qvac = cell%area * q2(1)
     203             :       end if
     204           5 :       if ( l_xyav ) return
     205           5 :       fact = ( qvac + psint ) / ( stars%nstr(1) * cell%vol )
     206           5 :       psq(1) = psq(1) - fact
     207           5 :       write( 6, fmt=8010 ) fact * 1000
     208             : 8010  format (/,10x,'                     1000 * normalization const. ='&
     209             :             &       ,5x,2f11.6)
     210             :     end if ! mpi%irank == 0 
     211             : 
     212             :   end subroutine psqpw
     213             : 
     214             : end module m_psqpw

Generated by: LCOV version 1.13