LCOV - code coverage report
Current view: top level - vgen - psqpw.F90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 93.3 % 89 83
Test Date: 2025-06-14 04:34:23 Functions: 100.0 % 1 1

            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          682 :   subroutine psqpw( fmpi, atoms, sphhar, stars, vacuum,  cell, input, sym, juphon,  &
      20          682 :        &     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          682 :     complex                         :: psq_local(stars%ng3)
      70          682 :     complex                         :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
      71          682 :     complex                         :: qlm(-atoms%lmaxd:atoms%lmaxd,0:atoms%lmaxd,atoms%ntype)
      72          682 :     complex                         :: q2(vacuum%nmzd)
      73          682 :     real                            :: q2r(vacuum%nmzd),q2i(vacuum%nmzd)
      74          682 :     real                            :: pn(0:atoms%lmaxd,atoms%ntype)
      75         1908 :     real                            :: aj(0:atoms%lmaxd+maxval(atoms%ncv)+1)
      76          682 :     real, allocatable, dimension(:) :: il, kl
      77          682 :     real                            :: g0(atoms%ntype)
      78          682 :     complex                         :: qpw(stars%ng3)
      79          682 :     real                            :: rho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      80          682 :     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          682 :     l_dfptvgen = PRESENT(stars2)
      88      1502016 :     qpw = den%pw(:,ispin)
      89     23987348 :     rho = den%mt(:,:,:,ispin)
      90        37830 :     IF (input%film) rht = den%vac(:,1,:,ispin)
      91          682 :     sigma_disc = cmplx(0.0,0.0)
      92              : 
      93              :     ! Calculate multipole moments
      94          682 :     call timestart("mpmom")
      95          682 :     IF (.NOT.l_dfptvgen) THEN
      96          682 :         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          682 :     call timestop("mpmom")
     109              : 
     110      1502016 :     psq(:) = cmplx( 0.0, 0.0 )
     111      1502016 :     psq_local(:) = cmplx( 0.0, 0.0 )
     112              : #ifdef CPP_MPI
     113          682 :     call MPI_BCAST( qpw, size(qpw), MPI_DOUBLE_COMPLEX, 0, fmpi%mpi_comm, ierr )
     114          682 :     nd = ( 2 * atoms%lmaxd + 1 ) * ( atoms%lmaxd + 1 ) * atoms%ntype
     115          682 :     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        13558 :     pn = 0.
     124         1908 :     do n = 1, atoms%ntype
     125         1908 :       if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     126        12088 :         do l = 0, min( atoms%ncv(n) - 1, atoms%lmax(n) )
     127        12088 :           pn(l, n) = DoubleFactorial( atoms%ncv(n) + 1, l ) / ( atoms%rmt(n) ** ( atoms%ncv(n) + 1 ) )
     128              :         end do
     129              :       else
     130           30 :         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         2728 :     if( fmpi%irank == 0 .AND. norm2(stars%center)<=1e-8) then
     142              :        s = 0.
     143          954 :        do n = 1, atoms%ntype
     144          954 :          if ( potdenType /= POTDEN_TYPE_POTYUK ) then
     145          598 :            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          341 :        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          682 :     fpo = 1. / cell%omtil
     156              : 
     157          682 :     call timestart("loop in psqpw")
     158              : 
     159         2728 :     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          682 :     !$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          682 :     CALL mpi_sum_reduce(psq_local,psq,fmpi%MPI_COMM)
     196              : 
     197          682 :     call timestop("loop in psqpw")
     198              : 
     199         2728 :     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          682 :     if ( fmpi%irank == 0 ) then
     205          341 :       if ( potdenType == POTDEN_TYPE_POTYUK ) return
     206              :       ! Check: integral of the pseudo charge density within the slab
     207          338 :       if ( input%film ) then
     208           37 :         psint = psq(1) * stars%nstr(1) * vacuum%dvac
     209       138057 :         do k = 2, stars%ng3
     210       138057 :           if ( stars%ig2(k) == 1 ) then
     211          730 :             gz = stars%kv3(3,k) * cell%bmat(3,3)
     212          730 :             f = 2. * sin( gz * cell%z1 ) / gz
     213          730 :             psint = psint + stars%nstr(k) * psq(k) * f
     214              :           end if
     215              :         end do
     216           37 :         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          338 :       write(oUnit, fmt=8000 ) psint
     221              : 8000  format (/,10x,'integral of pseudo charge density inside the slab='&
     222              :             &       ,5x,2f11.6)
     223          338 :       if ( .not. input%film .or. potdenType == POTDEN_TYPE_POTYUK ) return
     224              : 
     225              :       ! Normalized pseudo density
     226           37 :         qvac = cmplx(0.0,0.0)
     227           84 :         do ivac = 1, vacuum%nvac
     228           47 :           if (.not.l_dfptvgen) then
     229        11844 :             call qsf( vacuum%delz, real(rht(:,ivac)), q2r, vacuum%nmz, 0 )
     230        11797 :             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           47 :           q2(1) = q2(1) * cell%area
     237           84 :           qvac = qvac + q2(1) * 2. / real( vacuum%nvac )
     238              :         end do
     239              :         !TODO: reactivate electric fields
     240              :         !qvac = qvac - 2 * input%sigma
     241           37 :       if ( l_xyav ) return
     242           37 :       fact = ( qvac + psint ) / ( stars%nstr(1) * cell%vol )
     243           37 :       psq(1) = psq(1) - fact
     244              :       ! Find discontinuity at the boundaries
     245       138094 :       do k = 1, stars%ng3
     246       138094 :          if ( stars%ig2(k) == 1 ) then
     247          767 :             gz = stars%kv3(3,k) * cell%bmat(3,3)
     248          767 :             fc = cmplx(cos(gz * cell%z1),sin(gz * cell%z1))
     249          767 :             sigma_disc(1) = sigma_disc(1) - stars%nstr(k) * psq(k) * fc
     250          767 :             sigma_disc(2) = sigma_disc(2) + stars%nstr(k) * psq(k) * conjg(fc)
     251              :          end if
     252              :       end do
     253              :       sigma_disc(1) = sigma_disc(1) + rht(1,1)
     254              :       sigma_disc(2) = sigma_disc(2) - rht(1,2)
     255           37 :       sigma_disc = 0.0 
     256              : 8010  format (/,10x,'                     1000 * normalization const. ='&
     257              :             &       ,5x,2f11.6)
     258              :     end if ! fmpi%irank == 0
     259              : 
     260          682 :   end subroutine psqpw
     261              : 
     262              : end module m_psqpw
        

Generated by: LCOV version 2.0-1