LCOV - code coverage report
Current view: top level - vgen - psqpw.F90 (source / functions) Hit Total Coverage
Test: FLEUR test coverage Lines: 84 90 93.3 %
Date: 2024-03-29 04:21:46 Functions: 1 1 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         696 :   subroutine psqpw( fmpi, atoms, sphhar, stars, vacuum,  cell, input, sym, juphon,  &
      20         696 :        &     den, ispin, l_xyav, potdenType, psq, sigma_disc, rhoimag, stars2, iDtype, iDir, rho0, qpw0, iDir2, mat2ord )
      21             : 
      22             : #ifdef CPP_MPI
      23             :     use mpi
      24             : #endif
      25             :     use m_constants
      26             :     use m_phasy1
      27             :     use m_mpmom
      28             :     use m_sphbes
      29             :     use m_qsf
      30             :     USE m_mpi_reduce_tool
      31             :      
      32             :      
      33             :     use m_types
      34             :     use m_DoubleFactorial
      35             :     use m_SphBessel
      36             :     implicit none
      37             : 
      38             :     type(t_mpi),        intent(in)  :: fmpi
      39             :     type(t_atoms),      intent(in)  :: atoms
      40             :     type(t_sphhar),     intent(in)  :: sphhar
      41             :     type(t_stars),      intent(in)  :: stars
      42             :     type(t_vacuum),     intent(in)  :: vacuum
      43             :     type(t_cell),       intent(in)  :: cell
      44             :     type(t_input),      intent(in)  :: input
      45             :     type(t_sym),        intent(in)  :: sym
      46             :     type(t_juphon),     intent(in)  :: juphon
      47             :     type(t_potden),     intent(in)  :: den
      48             :      
      49             :     logical,            intent(in)  :: l_xyav
      50             :     !complex,            intent(in)  :: qpw(stars%ng3)
      51             :     !real,               intent(in)  :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      52             :     !real,               intent(in)  :: rht(vacuum%nmzd,2)
      53             :     integer,            intent(in)  :: potdenType, ispin
      54             :     complex,            intent(out) :: psq(stars%ng3)
      55             :     complex,            intent(out) :: sigma_disc(2)
      56             : 
      57             :     REAL, OPTIONAL, INTENT(IN)      :: rhoimag(atoms%jmtd,0:sphhar%nlhd,atoms%ntype), rho0(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      58             : 
      59             :     TYPE(t_stars), OPTIONAL, INTENT(IN) :: stars2
      60             :     COMPLEX, OPTIONAL, INTENT(IN)   :: qpw0(:)
      61             : 
      62             :     INTEGER, OPTIONAL, INTENT(IN)   :: iDtype, iDir ! DFPT: Type and direction of displaced atom
      63             :     INTEGER, OPTIONAL, INTENT(IN)   :: iDir2
      64             :     COMPLEX, OPTIONAL, INTENT(IN)   :: mat2ord(5,3,3)
      65             : 
      66             :     complex                         :: psint, sa, sl, sm, qvac, fact, fc
      67             :     real                            :: f, fpo, gz, p, rmtl, s, fJ, gr, g
      68             :     integer                         :: ivac, k, l, n, n1, nc, ncvn, lm, ll1, nd, m, nz, kStart, kEnd
      69         696 :     complex                         :: psq_local(stars%ng3)
      70         696 :     complex                         :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
      71         696 :     complex                         :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
      72         696 :     complex                         :: q2(vacuum%nmzd)
      73         696 :     real                            :: q2r(vacuum%nmzd),q2i(vacuum%nmzd)
      74         696 :     real                            :: pn(0:atoms%lmaxd,atoms%ntype)
      75        1950 :     real                            :: aj(0:atoms%lmaxd+maxval(atoms%ncv)+1)
      76         696 :     real, allocatable, dimension(:) :: il, kl
      77         696 :     real                            :: g0(atoms%ntype)
      78         696 :     complex                         :: qpw(stars%ng3)
      79         696 :     real                            :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      80         696 :     complex                         :: rht(vacuum%nmzd,2)
      81             :     LOGICAL :: l_dfptvgen ! If this is true, we handle things differently!
      82             : 
      83             : #ifdef CPP_MPI
      84             :     integer                         :: ierr
      85             : #endif
      86             : 
      87         696 :     l_dfptvgen = PRESENT(stars2)
      88     1562384 :     qpw = den%pw(:,ispin)
      89    24324510 :     rho = den%mt(:,:,:,ispin)
      90       44872 :     IF (input%film) rht = den%vac(:,1,:,ispin)
      91         696 :     sigma_disc = cmplx(0.0,0.0)
      92             : 
      93             :     ! Calculate multipole moments
      94         696 :     call timestart("mpmom")
      95         696 :     IF (.NOT.l_dfptvgen) THEN
      96         696 :         call mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell,   qpw, rho, potdenType, qlm )
      97           0 :     ELSE IF (PRESENT(iDir2)) THEN
      98             :         ! DFPT case:
      99             :         ! Additional contributions to qlm due to surface corrections.
     100             :         call mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell,   qpw, rho, potdenType, qlm, &
     101           0 :                   & rhoimag=rhoimag, stars2=stars2, iDtype=iDtype, iDir=iDir, rho0=rho0, qpw0=qpw0, iDir2=iDir2, mat2ord=mat2ord )
     102             :     ELSE
     103             :         ! DFPT case:
     104             :         ! Additional contributions to qlm due to surface corrections.
     105             :         call mpmom( input, fmpi, atoms, sphhar, stars, sym, juphon, cell,   qpw, rho, potdenType, qlm, &
     106           0 :                   & rhoimag=rhoimag, stars2=stars2, iDtype=iDtype, iDir=iDir, rho0=rho0, qpw0=qpw0 )
     107             :     END IF
     108         696 :     call timestop("mpmom")
     109             : 
     110     1562384 :     psq(:) = cmplx( 0.0, 0.0 )
     111     1562384 :     psq_local(:) = cmplx( 0.0, 0.0 )
     112             : #ifdef CPP_MPI
     113         696 :     call MPI_BCAST( qpw, size(qpw), MPI_DOUBLE_COMPLEX, 0, fmpi%mpi_comm, ierr )
     114         696 :     nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
     115         696 :     call MPI_BCAST( qlm, nd, MPI_DOUBLE_COMPLEX, 0, fmpi%MPI_COMM, ierr )
     116             : #endif
     117             : 
     118             :     ! prefactor in (A10) (Coulomb case) or (A11) (Yukawa case)
     119             :     ! nc(n) is the integer p in the paper; ncv(n) is l + p
     120             :     ! Coulomb case: pn(l,n) = (2 * l + 2 * p + 3)!! / ( (2 * l + 1)!! * R ** (ncv(n) + 1) ),
     121             :     ! Yukawa case: pn(l,n) = lambda ** (l + p + 1) / ( i_{l+p+1}(lambda * R) * (2 * l + 1)!! )
     122             :     ! g0 is the prefactor for the q=0 component in (A13)
     123       13796 :     pn = 0.
     124        1950 :     do n = 1, atoms%ntype
     125        1950 :       if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     126       12256 :         do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
     127       12256 :           pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
     128             :         end do
     129             :       else
     130         120 :         allocate( il(0:atoms%ncv(n)+1), kl(0:atoms%ncv(n)+1) )
     131          30 :         call ModSphBessel( il(0:), kl(0:), input%preconditioning_param * atoms%rmt(n), atoms%ncv(n) + 1 )
     132          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)
     133         300 :         do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
     134         300 :           pn(l, n) = input%preconditioning_param ** ( atoms%ncv(n) + 1 ) / il( atoms%ncv(n) + 1 ) / DoubleFactorial( l )
     135             :         end do
     136          30 :         deallocate( il, kl )
     137             :       end if
     138             :     end do
     139             : 
     140             :     ! q=0 term: see (A12) (Coulomb case) or (A13) (Yukawa case)
     141        2784 :     if( fmpi%irank == 0 .AND. norm2(stars%center)<=1e-8) then
     142             :        s = 0.
     143         975 :        do n = 1, atoms%ntype
     144         975 :          if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     145         612 :            s = s + atoms%neq(n) * real( qlm(0,0,n) )
     146             :          else
     147          15 :            s = s + atoms%neq(n) * real( qlm(0,0,n) ) * g0(n)
     148             :          end if
     149             :        end do
     150             :     !if( fmpi%irank == 0 ) then
     151         348 :        psq_local(1) = qpw(1) + ( sfp_const / cell%omtil ) * s
     152             :     end if
     153             : 
     154             :     ! q/=0 term: see (A10) (Coulomb case) or (A11) (Yukawa case)
     155         696 :     fpo = 1. / cell%omtil
     156             : 
     157         696 :     call timestart("loop in psqpw")
     158             : 
     159        2784 :     CALL calcIndexBounds(fmpi, MERGE(2,1,norm2(stars%center)<=1e-8), stars%ng3, kStart, kEnd)
     160             :     !$OMP parallel do default( NONE ) &
     161             :     !$OMP SHARED(atoms,stars,sym,cell,kStart,kEnd,psq_local,qpw,qlm,pn,fpo,iDir,iDir2,iDtype) &
     162         696 :     !$OMP private( pylm, sa, n, ncvn, aj, sl, l, n1, ll1, sm, m, lm )
     163             :     do k = kStart, kEnd
     164             :       call phasy1( atoms, stars, sym, cell, k, pylm )
     165             :       sa = 0.0
     166             :       do n = 1, atoms%ntype
     167             :         ncvn = atoms%ncv(n)
     168             :         call sphbes( ncvn + 1 , stars%sk3(k) * atoms%rmt(n), aj )
     169             :         sl = 0.
     170             :         do l = 0, atoms%lmax(n)
     171             :           IF(l.LT.ncvn) THEN
     172             :              n1 = ncvn - l + 1
     173             :              ll1 = l * ( l + 1 ) + 1
     174             :              sm = 0.
     175             :              do m = -l, l
     176             :                lm = ll1 + m
     177             :                sm = sm + qlm(m,l,n) * conjg( pylm(lm,n) )
     178             :              end do
     179             :           END IF
     180             :         sl = sl + pn(l,n) / ( stars%sk3(k) ** n1 ) * aj( ncvn + 1 ) * sm
     181             :         end do
     182             :         sa = sa + atoms%neq(n) * sl
     183             :         IF (PRESENT(iDir2)) THEN
     184             :           IF (iDir2==iDir) THEN
     185             :             IF (n==iDtype.OR.0==iDtype) THEN 
     186             :               sa = sa + atoms%zatom(n) * 5 * pn(2,n) * conjg(pylm(1,n)) * aj( ncvn + 1 ) / sfp_const / ( stars%sk3(k) ** (ncvn-1) )
     187             :             END IF
     188             :           END IF
     189             :         END IF
     190             :       end do
     191             :       psq_local(k) = qpw(k) + fpo * sa
     192             :     end do
     193             :     !$omp end parallel do
     194             : 
     195         696 :     CALL mpi_sum_reduce(psq_local,psq,fmpi%MPI_COMM)
     196             : 
     197         696 :     call timestop("loop in psqpw")
     198             : 
     199        2784 :     IF (.NOT.norm2(stars%center)<1e-8) RETURN ! TODO: Change this!
     200             :     !IF (l_dfptvgen) RETURN
     201             :     ! TODO: As of yet unclear if we need to do somecorrection here for
     202             :     !       DFPT, to make up for possible charge shifts from/to the vacuum.
     203             : 
     204         696 :     if ( fmpi%irank == 0 ) then
     205         348 :       if ( potdenType == POTDEN_TYPE_POTYUK ) return
     206             :       ! Check: integral of the pseudo charge density within the slab
     207         345 :       if ( input%film ) then
     208          44 :         psint = psq(1) * stars%nstr(1) * vacuum%dvac
     209      168234 :         do k = 2, stars%ng3
     210      168234 :           if ( stars%ig2(k) == 1 ) then
     211         800 :             gz = stars%kv3(3,k) * cell%bmat(3,3)
     212         800 :             f = 2. * sin( gz * cell%z1 ) / gz
     213         800 :             psint = psint + stars%nstr(k) * psq(k) * f
     214             :           end if
     215             :         end do
     216          44 :         psint = cell%area * psint
     217             :       else if ( .not. input%film ) then
     218         301 :         psint = psq(1) * stars%nstr(1) * cell%omtil
     219             :       end if
     220         345 :       write(oUnit, fmt=8000 ) psint
     221             : 8000  format (/,10x,'integral of pseudo charge density inside the slab='&
     222             :             &       ,5x,2f11.6)
     223         345 :       if ( .not. input%film .or. potdenType == POTDEN_TYPE_POTYUK ) return
     224             : 
     225             :       ! Normalized pseudo density
     226          44 :         qvac = cmplx(0.0,0.0)
     227          98 :         do ivac = 1, vacuum%nvac
     228          54 :           if (.not.l_dfptvgen) then
     229       13608 :             call qsf( vacuum%delz, real(rht(:,ivac)), q2r, vacuum%nmz, 0 )
     230       13554 :             q2 = q2r
     231             :           else
     232           0 :             call qsf( vacuum%delz, real(rht(:,ivac)), q2r, vacuum%nmz, 0 )
     233           0 :             call qsf( vacuum%delz,aimag(rht(:,ivac)), q2i, vacuum%nmz, 0 )
     234           0 :             q2 = q2r + ImagUnit*q2i
     235             :           end if
     236          54 :           q2(1) = q2(1) * cell%area
     237          98 :           qvac = qvac + q2(1) * 2. / real( vacuum%nvac )
     238             :         end do
     239             :         !TODO: reactivate electric fields
     240             :         !qvac = qvac - 2 * input%sigma
     241          44 :       if ( l_xyav ) return
     242          44 :       fact = ( qvac + psint ) / ( stars%nstr(1) * cell%vol )
     243          44 :       psq(1) = psq(1) - fact
     244             :       ! Find discontinuity at the boundaries
     245      168278 :       do k = 1, stars%ng3
     246      168278 :          if ( stars%ig2(k) == 1 ) then
     247         844 :             gz = stars%kv3(3,k) * cell%bmat(3,3)
     248         844 :             fc = cmplx(cos(gz * cell%z1),sin(gz * cell%z1))
     249         844 :             sigma_disc(1) = sigma_disc(1) - stars%nstr(k) * psq(k) * fc
     250         844 :             sigma_disc(2) = sigma_disc(2) + stars%nstr(k) * psq(k) * conjg(fc)
     251             :          end if
     252             :       end do
     253          44 :       sigma_disc(1) = sigma_disc(1) + rht(1,1)
     254          44 :       sigma_disc(2) = sigma_disc(2) - rht(1,2)
     255             : 8010  format (/,10x,'                     1000 * normalization const. ='&
     256             :             &       ,5x,2f11.6)
     257             :     end if ! fmpi%irank == 0
     258             : 
     259         696 :   end subroutine psqpw
     260             : 
     261             : end module m_psqpw

Generated by: LCOV version 1.14