LCOV - code coverage report
Current view: top level - vgen - vmts.F90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 69.7 % 89 62
Test Date: 2025-06-10 04:26:50 Functions: 100.0 % 1 1

            Line data    Source code
       1              : module m_vmts
       2              : #ifdef CPP_MPI
       3              :   use mpi
       4              : #endif
       5              : contains
       6              : 
       7          682 :   subroutine vmts( input, fmpi, stars, sphhar, atoms, sym, cell, juphon, dosf, vpw, rho, potdenType, vr, rhoIm, vrIm, iDtype, iDir, iDir2, mat2ord )
       8              : 
       9              :   !-------------------------------------------------------------------------
      10              :   ! This subroutine calculates the lattice harmonics expansion coefficients
      11              :   ! of the Coulomb / Yukawa potential for all atom types.
      12              :   !
      13              :   ! In more detail:
      14              :   ! the radius-dependent coefficient of the potential is
      15              :   ! V_{lm}(r) = \int_0^R G_l(r,r') \rho_{lm}(r') r'^2 dr' + V_{lm}(R) P_l(r)
      16              :   ! where
      17              :   ! [sphere boundary contribution - first part of code:]
      18              :   ! V_{lm}(R) is the spherical harmonics expansion of the interstitial
      19              :   ! potential,
      20              :   ! P_l(r) is the derivative of the spherical Green's function on the sphere
      21              :   ! boundary (times a constant),
      22              :   ! [sphere interior contribution - second part of code:]
      23              :   ! G_l(r,r') ) = green_factor u_1(r_<) u_2(r_>) is the spherical Green's
      24              :   ! function with homogeneous solutions u_1 and u_2 and
      25              :   ! r_< = min(r,r') and r_> = max(r,r')
      26              :   ! The integral is split in a part where r'=r_< and a part where r'=r_> and
      27              :   ! the integral from r to R is split in \int_0^R - \int_0^r. Resulting
      28              :   ! terms depending solely on R (not on r) are summarised in the variable
      29              :   ! termsR.
      30              :   !
      31              :   ! More information and equations can be found in
      32              :   ! F. Tran, P. Blaha: Phys. Rev. B 83, 235118(2011)
      33              :   !-------------------------------------------------------------------------
      34              : 
      35              :     use m_constants
      36              :     use m_types
      37              :     use m_mpi_reduce_tool
      38              :     use m_intgr, only : intgr2, sfint
      39              :     use m_phasy1
      40              :     use m_sphbes
      41              :      
      42              :     use m_SphBessel
      43              :     !$ use omp_lib
      44              :     implicit none
      45              : 
      46              :     type(t_input),  intent(in)        :: input
      47              :     type(t_mpi),    intent(in)        :: fmpi
      48              :     type(t_stars),  intent(in)        :: stars
      49              :     type(t_sphhar), intent(in)        :: sphhar
      50              :     type(t_atoms),  intent(in)        :: atoms
      51              :     type(t_sym),    intent(in)        :: sym
      52              :     type(t_cell),   intent(in)        :: cell
      53              :     type(t_juphon), intent(in)        :: juphon
      54              :      
      55              :     LOGICAL,        INTENT(IN)        :: dosf
      56              :     complex,        intent(in)        :: vpw(:)!(stars%ng3,input%jspins)
      57              :     real,           intent(in)        :: rho(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      58              :     integer,        intent(in)        :: potdenType
      59              :     real,           intent(out)       :: vr(:,0:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
      60              :     REAL,    OPTIONAL, INTENT(IN)     :: rhoIm(:,0:,:)
      61              :     REAL,    OPTIONAL, INTENT(OUT)    :: vrIm(:,0:,:)
      62              :     INTEGER, OPTIONAL, INTENT(IN)     :: iDtype, iDir
      63              :     INTEGER, OPTIONAL, INTENT(IN)     :: iDir2
      64              :     COMPLEX, OPTIONAL, INTENT(IN)     :: mat2ord(5,3,3)
      65              : 
      66              :     complex                           :: cp, sm
      67              :     integer                           :: i, jm, k, l, lh, n, nd, lm, m, imax, lmax, iMem, ptsym
      68              :     integer                           :: maxBunchSize, numBunches, iBunch, firstStar, iTempArray
      69              :     complex                           :: temp
      70          682 :     complex                           :: vtl(0:sphhar%nlhd, atoms%ntype)
      71          682 :     complex                           :: pylm(( atoms%lmaxd + 1 ) ** 2, atoms%ntype)
      72              :     real                              :: green_factor, termsR, pref
      73          682 :     real                              :: green_1    (1:atoms%jmtd), green_2    (1:atoms%jmtd)
      74          682 :     real                              :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
      75          682 :     real                              :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)!, integral_3 (1:atoms%jmtd)
      76          682 :     real                              :: sbf(0:atoms%lmaxd)
      77          682 :     real, allocatable, dimension(:,:) :: il, kl
      78              :     LOGICAL                           :: l_dfptvgen
      79          682 :     COMPLEX, ALLOCATABLE              :: vtlStars(:,:,:), vtlLocal(:,:)
      80              :     TYPE(t_parallelLoop)              :: mpiLoop, ompLoop
      81              : 
      82          682 :     l_dfptvgen = PRESENT(rhoIm)
      83              : 
      84              :     ! SPHERE BOUNDARY CONTRIBUTION to the coefficients calculated from the values
      85              :     ! of the interstitial Coulomb / Yukawa potential on the sphere boundary
      86              : 
      87         1364 :     ALLOCATE (vtlLocal(0:sphhar%nlhd,atoms%ntype))
      88        34318 :     vtlLocal(:,:) = cmplx(0.0,0.0)
      89              : 
      90         2728 :     firstStar = MERGE(2,1,norm2(stars%center)<=1e-8)
      91          682 :     maxBunchSize = 2*getNumberOfThreads() ! This bunch size is kind of a magic number determined from some
      92              :                                           ! naive performance tests for a 64 atom unit cell
      93          682 :     CALL calcNumberComputationBunches(firstStar, stars%ng3, maxBunchSize, numBunches)
      94              : 
      95          682 :     CALL mpiLoop%init(fmpi%irank,fmpi%isize,0,numBunches-1)
      96              : 
      97         2046 :     ALLOCATE(vtlStars(0:sphhar%nlhd,atoms%ntype,maxBunchSize))
      98       137954 :     vtlStars = CMPLX(0.0,0.0)
      99              : 
     100              :     ! q/=0 components
     101       188327 :     DO iBunch = mpiLoop%bunchMinIndex, mpiLoop%bunchMaxIndex
     102       187645 :        CALL ompLoop%init(iBunch,numBunches,firstStar,stars%ng3)
     103              :        !$OMP parallel do default( NONE ) &
     104              :        !$OMP SHARED(ompLoop,atoms,stars,sym,cell,sphhar,vpw,vtlStars,potdenType) &
     105       188327 :        !$OMP private(iTempArray,cp,pylm,n,sbf,nd,lh,l,sm,jm,m,lm)
     106              :        do k = ompLoop%bunchMinIndex, ompLoop%bunchMaxIndex
     107              :           iTempArray = k - ompLoop%bunchMinIndex + 1
     108              :           cp = vpw(k) * stars%nstr(k)
     109              :           call phasy1( atoms, stars, sym, cell, k, pylm )
     110              :           do n = 1, atoms%ntype
     111              :              call sphbes( atoms%lmax(n), stars%sk3(k) * atoms%rmt(n), sbf )
     112              :              nd = sym%ntypsy(atoms%firstAtom(n))
     113              :              do lh = 0, sphhar%nlh(nd)
     114              :                 l = sphhar%llh(lh,nd)
     115              :                 sm = (0.0,0.0)
     116              :                 do jm = 1, sphhar%nmem(lh,nd)
     117              :                    m = sphhar%mlh(jm,lh,nd)
     118              :                    lm = l * ( l + 1 ) + m + 1
     119              :                    sm = sm + conjg( sphhar%clnu(jm,lh,nd) ) * pylm(lm,n)
     120              :                 end do
     121              :                 vtlStars(lh,n,iTempArray) = vtlStars(lh,n,iTempArray) + cp * sbf(l) * sm
     122              :              end do
     123              :           end do
     124              :        end do
     125              :        !$OMP END PARALLEL DO
     126              :     END DO
     127              :     
     128              :     !$OMP parallel do default( NONE ) &
     129              :     !$OMP SHARED(atoms,sym,sphhar,maxBunchSize,vtlStars,vtlLocal) &
     130          682 :     !$OMP private(nd,lh,temp,iTempArray)
     131              :     DO n = 1, atoms%ntype
     132              :        nd = sym%ntypsy(atoms%firstAtom(n))
     133              :        do lh = 0, sphhar%nlh(nd)
     134              :           temp = CMPLX(0.0,0.0)
     135              :           DO iTempArray = 1, maxBunchSize
     136              :              temp = temp + vtlStars(lh,n,iTempArray)
     137              :           END DO
     138              :           vtlLocal(lh,n) = temp
     139              :        END DO
     140              :     END DO
     141              :     !$OMP END PARALLEL DO
     142              :     
     143              :     ! q=0 component
     144         2728 :     if ( fmpi%irank == 0 .AND. norm2(stars%center)<=1e-8 ) then
     145          954 :        DO n = 1, atoms%ntype
     146          954 :           vtlLocal(0,n) = vtlLocal(0,n) + sfp_const * vpw(1)
     147              :        END DO
     148              :     end if
     149              : 
     150          682 :     CALL mpi_sum_reduce(vtlLocal, vtl,fmpi%mpi_comm)
     151              : 
     152              :     ! SPHERE INTERIOR CONTRIBUTION to the coefficients calculated from the
     153              :     ! values of the sphere Coulomb/Yukawa potential on the sphere boundary
     154              : 
     155          682 :     if( fmpi%irank == 0 ) then
     156          341 :     if ( potdenType == POTDEN_TYPE_POTYUK ) then
     157            6 :       allocate( il(0:atoms%lmaxd, 1:atoms%jmtd), kl(0:atoms%lmaxd, 1:atoms%jmtd) )
     158              :     end if
     159              : 
     160          954 :     do n = 1, atoms%ntype
     161          613 :       nd = sym%ntypsy(atoms%firstAtom(n))
     162          613 :       imax = atoms%jri(n)
     163          613 :       lmax = sphhar%llh(sphhar%nlh(nd), nd)
     164          613 :       if ( potdenType == POTDEN_TYPE_POTYUK ) then
     165              :         !do concurrent (i = 1:imax)
     166        11370 :         do i = 1,imax
     167        11370 :           call ModSphBessel( il(0:,i), kl(0:,i), input%preconditioning_param * atoms%rmsh(i,n), lmax )
     168              :         end do
     169              :       end if
     170        16477 :       do lh = 0, sphhar%nlh(nd)
     171        15523 :         l = sphhar%llh(lh,nd)
     172        15523 :         if ( potdenType == POTDEN_TYPE_POTYUK ) then
     173       156906 :           green_1(1:imax) = il(l,1:imax)
     174       156906 :           green_2(1:imax) = kl(l,1:imax)
     175          207 :           green_factor    = fpi_const * input%preconditioning_param
     176              :         else
     177     11053930 :           green_1(1:imax) = atoms%rmsh(1:imax,n) ** l
     178     11053930 :           green_2(1:imax) = 1.0 / ( green_1(1:imax) * atoms%rmsh(1:imax,n) )
     179        15316 :           green_factor    = fpi_const / ( 2 * l + 1 )
     180              :         end if
     181              : 
     182     11210836 :         integrand_1(1:imax) = green_1(1:imax) * rho(1:imax,lh,n)
     183     11210836 :         integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
     184        15523 :         if (.not.dosf) THEN
     185        15361 :          call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
     186        15361 :          call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
     187              :         else
     188          162 :            call sfint(integrand_1(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_1(1:imax))
     189          162 :            call sfint(integrand_2(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_2(1:imax))
     190              :         end if
     191        15523 :         termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
     192              :         vr(1:imax,lh,n) = green_factor * (   green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
     193     11210836 :                                            + green_2(1:imax) *            integral_1(1:imax)   )
     194        16136 :         IF (l_dfptvgen) THEN
     195              :            ! Integrate the imaginary part of the density perturbation as well.
     196            0 :            integrand_1(1:imax) = green_1(1:imax) * rhoIm(1:imax,lh,n)
     197            0 :            integrand_2(1:imax) = green_2(1:imax) * rhoIm(1:imax,lh,n)
     198            0 :            call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
     199            0 :            call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
     200            0 :            termsR = integral_2(imax) + ( AIMAG(vtl(lh,n)) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
     201              :            vrIm(1:imax,lh,n) = green_factor * (   green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
     202            0 :                                                 + green_2(1:imax) *            integral_1(1:imax)   )
     203              :         END IF
     204              :       end do
     205              :     end do
     206          341 :     if ( potdenType == POTDEN_TYPE_POTYUK ) then
     207            3 :       deallocate( il, kl )
     208              :     end if
     209              : 
     210          341 :     if ( potdenType /= POTDEN_TYPE_POTYUK .AND. potdenType /= POTDEN_TYPE_CRYSTALFIELD) then
     211          337 :       IF (.NOT.l_dfptvgen) THEN
     212          930 :          do n = 1, atoms%ntype
     213       412483 :          vr(1:atoms%jri(n),0,n) = vr(1:atoms%jri(n),0,n) - sfp_const * ( 1.0 / atoms%rmsh(1:atoms%jri(n),n) - 1.0 / atoms%rmt(n) ) * atoms%zatom(n)
     214              :          end do
     215            0 :       ELSE IF (juphon%l_phonon) THEN
     216            0 :          IF (.NOT.PRESENT(iDir2)) THEN
     217              :             ! DFPT(-phonon) case:
     218              :             ! l=1 contributions from the Coulomb singularity instead of l=0 (1/r -> 1/r^2)
     219            0 :             DO n = MERGE(1,iDtype,iDtype==0), MERGE(atoms%ntype,iDtype,iDtype==0)
     220            0 :                ptsym = sym%ntypsy(atoms%firstAtom(n))
     221            0 :                pref = MERGE(atoms%zatom(n),-atoms%zatom(n),iDtype==0)
     222            0 :                DO lh = 1, 3
     223            0 :                   l = sphhar%llh(lh, ptsym)
     224            0 :                   DO iMem = 1, sphhar%nmem(lh, ptsym)
     225            0 :                      m = sphhar%mlh(iMem, lh, ptsym)
     226            0 :                      lm = l*(l+1) + m + 1
     227              :                      vr(1:atoms%jri(n),lh,n) = vr(1:atoms%jri(n),lh,n) + &
     228              :                                                 conjg(sphhar%clnu(iMem, lh, ptsym)) * c_im(iDir, lm - 1) * pref * &
     229            0 :                                                 ( 1 - (atoms%rmsh(1:atoms%jri(n), n) / atoms%rmt(n))**3) / atoms%rmsh(1:atoms%jri(n),n)**2
     230              :                   END DO
     231              :                END DO
     232              :             END DO
     233              :          ELSE
     234              :             ! DFPT 2nd order case:
     235              :             ! l=2 contributions from the Coulomb singularity instead of l=0 (1/r -> 1/r^3)
     236            0 :             DO n = 1, atoms%ntype!MERGE(1,iDtype,iDtype==0), MERGE(atoms%ntype,iDtype,iDtype==0)
     237            0 :                ptsym = sym%ntypsy(atoms%firstAtom(n))
     238            0 :                pref = -atoms%zatom(n)!MERGE(atoms%zatom(n),-atoms%zatom(n),iDtype==0)
     239            0 :                DO lh = 4, 8
     240            0 :                   l = sphhar%llh(lh, ptsym)
     241            0 :                   DO iMem = 1, sphhar%nmem(lh, ptsym)
     242            0 :                      m = sphhar%mlh(iMem, lh, ptsym)
     243            0 :                      lm = l*(l+1) + m + 1
     244            0 :                      IF ((n.EQ.iDtype).OR.(0.EQ.iDtype)) vr(1:atoms%jri(n),lh,n) = vr(1:atoms%jri(n),lh,n) + &
     245              :                                                          conjg(sphhar%clnu(iMem, lh, ptsym)) * mat2ord(lm-4,iDir2,iDir) * pref * &
     246            0 :                                                          ( 1 - (atoms%rmsh(1:atoms%jri(n), n) / atoms%rmt(n))**5) / atoms%rmsh(1:atoms%jri(n),n)**3
     247              :                   END DO
     248              :                END DO
     249              :             END DO
     250              :          END IF
     251              :       END IF
     252              :     end if
     253              :     end if
     254              : 
     255         1364 :   end subroutine vmts
     256              : 
     257              : end module m_vmts
        

Generated by: LCOV version 2.0-1