phononevents_tetrahedron.f90 Source File


Source Code

!> Tetrahedron integration weights for isotope scattering, from one q-point
subroutine iso_tetrahedron_fft_oneqp(qp,dr,scq,gi1,scsigma,blochlcorrections)
    !> qpoint mesh
    !type(lo_fft_mesh), intent(in) :: qp
    class(lo_qpoint_mesh), intent(in) :: qp
    !> phonon dispersoins
    type(lo_phonon_dispersions), intent(in) :: dr
    !> current q-point
    type(lo_iso2), intent(inout) :: scq
    !> current grid-index
    integer, intent(in) :: gi1
    !> baseline smearing parameter
    real(flyt), intent(in) :: scsigma
    !> blochl corrections?
    logical, intent(in) :: blochlcorrections
    !
    real(flyt), parameter :: thres_weight=lo_tiny
    integer, dimension(4) :: tetqpoints
    integer :: gi2
    integer :: i,j,ii,b1,b2,i_gamma
    real(flyt), dimension(4) :: cval,wts1,wts2
    real(flyt) :: om1,omthres,minc,maxc
    real(flyt), dimension(:,:,:), allocatable :: qpw
    real(flyt), dimension(:,:), allocatable :: omr2
    real(flyt), dimension(:), allocatable :: omr1
    
    lo_allocate(omr1(dr%nb))
    lo_allocate(omr2(4,dr%nb))
    lo_allocate(qpw(dr%nb,dr%nb,qp%nq_tot))
    omthres=dr%omega_min*0.2_flyt


    ! Calculate all integration weights for the relevant tetrahedra
    scq%gi1=gi1
    qpw=0.0_flyt
    omr1=dr%aq(gi1)%omega
    tetloop: do i=1,qp%ntet
        ! Grab frequencies
        do j=1,4
            gi2=qp%tet(i)%gridind(j)
            omr2(j,:)=dr%aq(gi2)%omega
            tetqpoints(j)=gi2
        enddo
        
        ! Add integration weights
        do b1=1,dr%nb
        do b2=1,dr%nb
            minc=lo_huge
            maxc=-lo_huge
            do j=1,4
                cval(j)=omr2(j,b2)
                minc=min(minc,cval(j))
                maxc=max(maxc,cval(j))
            enddo
            ! add a small delta for safety
            minc=minc-3*scsigma/100.0_flyt
            maxc=maxc+3*scsigma/100.0_flyt
            ! check if it's relevant
            om1=omr1(b1)
            if ( om1 .gt. minc .and. om1 .lt. maxc .and. om1 .gt. omthres ) then
                ! Take care if it's completely degenerate
                if ( lo_stddev(cval)/om1 .lt. 1E-3_flyt ) then
                    call lo_integration_weights_for_one_tetrahedron(&
                    qp%tet(i),cval,om1-scsigma/200.0_flyt,wts1,scsigma,blochlcorrections)
                    call lo_integration_weights_for_one_tetrahedron(&
                    qp%tet(i),cval,om1+scsigma/200.0_flyt,wts2,scsigma,blochlcorrections)
                else
                    call lo_integration_weights_for_one_tetrahedron(&
                    qp%tet(i),cval,om1-scsigma/200.0_flyt,wts1,scsigma,blochlcorrections)
                    call lo_integration_weights_for_one_tetrahedron(&
                    qp%tet(i),cval,om1+scsigma/200.0_flyt,wts2,scsigma,blochlcorrections)
                endif
                ! sum up integration weights
                qpw(b1,b2,tetqpoints)=qpw(b1,b2,tetqpoints)+&
                (wts1+wts2)*0.5_flyt
            endif
        enddo
        enddo
    enddo tetloop
   
    ! Set the integration weights for the acoustic branches to zero
    i_gamma=-1
    do i=1,qp%nq_tot
        if ( lo_sqnorm(qp%ap(i)%w) .lt. lo_sqtol ) then
            i_gamma=i
            exit
        endif
    enddo
    if ( i_gamma .lt. 0 ) call lo_stop_gracefully(['FFT mesh does not contain gamma'],lo_exitcode_symmetry,__FILE__,__LINE__)

    do b2=1,dr%nb
    do b1=1,dr%nb
        if ( dr%aq(i_gamma)%omega(b1) .lt. lo_freqtol .or. dr%aq(i_gamma)%omega(b2) .lt. lo_freqtol ) then
            qpw(b1,b2,i_gamma)=0.0_flyt
        endif
    enddo
    enddo
 
    ! Store weights per q-point, and only the relevant q-points
    do b1=1,dr%nb
    do b2=1,dr%nb
        ! count q-points
        ii=0
        do i=1,qp%nq_tot
            if ( abs(qpw(b1,b2,i)) .gt. thres_weight ) then
                ii=ii+1
            endif
        enddo
        ! make some space
        scq%band(b1,b2)%n=ii
        if ( scq%band(b1,b2)%n .gt. 0 ) then
           lo_allocate(scq%band(b1,b2)%e( ii ))
        endif
        ! store weights and indices
        ii=0
        do i=1,qp%nq_tot
            if ( abs(qpw(b1,b2,i)) .gt. thres_weight ) then
                ii=ii+1
                scq%band(b1,b2)%e(ii)%gi2=i
                scq%band(b1,b2)%e(ii)%deltafunction=qpw(b1,b2,i)
            endif
        enddo
    enddo
    enddo
    
    lo_deallocate(omr1)
    lo_deallocate(omr2)
    lo_deallocate(qpw)
end subroutine

!> Tetrahedron integration weights for threephonon scattering, from one q-point
subroutine threephonon_tetrahedron_fft_oneqp(qp,dr,scq,gi1,scsigma,blochlcorrections)
    !> qpoint mesh
    !type(lo_fft_mesh), intent(in) :: qp
    class(lo_qpoint_mesh), intent(in) :: qp
    !> phonon dispersions
    type(lo_phonon_dispersions), intent(in) :: dr
    !> q-point in question
    type(lo_3phqp2), intent(inout) :: scq
    !> grid-index in question
    integer, intent(in) :: gi1
    !> baseline smearing
    real(flyt), intent(in) :: scsigma
    !> blochl corrections?
    logical, intent(in) :: blochlcorrections

    integer :: i,j,ii,jj,b1,b2,b3
    integer :: gi2,gi3,i_gamma
    integer, dimension(3) :: dims
    integer, dimension(4) :: tetqpoints
    real(flyt), dimension(:,:,:,:), allocatable :: qpwp,qpwm
    real(flyt), dimension(:,:), allocatable :: omr2,omr3
    real(flyt), dimension(:), allocatable :: omr1
    real(flyt), dimension(4) :: cval_p,cval_m,wts1,wts2
    real(flyt) :: om1,omthres,minc_m,minc_p,maxc_m,maxc_p,wthres,sigmatol,deltaeps

    ! Some temporary space
    lo_allocate(omr1(dr%nb))
    lo_allocate(omr2(4,dr%nb))
    lo_allocate(omr3(4,dr%nb))
    lo_allocate(qpwp(dr%nb,dr%nb,dr%nb,qp%nq_tot))
    lo_allocate(qpwm(dr%nb,dr%nb,dr%nb,qp%nq_tot))
    ! smallest frequency allowed, should allow everything except acoustic modes at Gamma
    omthres=dr%omega_min*0.5_flyt
    ! smallest weight to care about
    wthres=1E-30_flyt !lo_tiny ! 1E-40_flyt
    ! Dimensions of q-point grid
    select type(qp)
    class is(lo_fft_mesh)
        dims=qp%griddensity
    class default
        write(*,*) 'Really need an fft grid for this'
        stop
    end select
    ! The reference q-point, q1 in q1+q2+q3=G
    scq%gi1=gi1
    ! Some tolerances
    sigmatol=scsigma/30.0_flyt
    deltaeps=scsigma/1E2_flyt
    qpwp=0.0_flyt
    qpwm=0.0_flyt
    omr1=dr%aq(gi1)%omega
    tetloop: do i=1,qp%ntet
        ! Grab frequencies        
        do j=1,4
            gi2=qp%tet(i)%gridind(j)
            gi3=fft_third_grid_index(gi1,gi2,dims)
            omr2(j,:)=dr%aq(gi2)%omega
            omr3(j,:)=dr%aq(gi3)%omega
            tetqpoints(j)=gi2
        enddo
        ! Add integration weights        
        do b1=1,dr%nb
        do b2=1,dr%nb
        do b3=1,dr%nb
            minc_p=lo_huge
            maxc_p=-lo_huge
            minc_m=lo_huge
            maxc_m=-lo_huge
            do j=1,4
                cval_p(j)=omr3(j,b3)-omr2(j,b2)   ! plus events
                cval_m(j)=omr3(j,b3)+omr2(j,b2)   ! minus events
                minc_p=min(minc_p,cval_p(j))
                minc_m=min(minc_m,cval_m(j))
                maxc_p=max(maxc_p,cval_p(j))
                maxc_m=max(maxc_m,cval_m(j))
            enddo
            ! add a small delta for safety, to avid pathological cases
            minc_p=minc_p-sigmatol
            minc_m=minc_m-sigmatol
            maxc_p=maxc_p+sigmatol
            maxc_m=maxc_m+sigmatol
            ! check if it's relevant
            om1=omr1(b1)
            if ( om1 .gt. minc_p .and. om1 .lt. maxc_p ) then !.and. om1 .gt. omthres ) then
                if ( blochlcorrections ) then
                    wts1=lo_LV_tetrahedron_weights(cval_p,om1,lo_freqtol,scsigma)
                    qpwp(b1,b2,b3,tetqpoints)=qpwp(b1,b2,b3,tetqpoints)+wts1*qp%tet(i)%weight
                else
                    wts1=0.0_flyt
                    wts2=0.0_flyt
                    ! this plus event is relevant, store weights
                    call lo_integration_weights_for_one_tetrahedron(qp%tet(i),cval_p,om1-deltaeps,wts1,scsigma,blochlcorrections)
                    call lo_integration_weights_for_one_tetrahedron(qp%tet(i),cval_p,om1+deltaeps,wts2,scsigma,blochlcorrections)
                    qpwp(b1,b2,b3,tetqpoints)=qpwp(b1,b2,b3,tetqpoints)+(wts1+wts2)*0.5_flyt
                endif
            endif
            if ( om1 .gt. minc_m .and. om1 .lt. maxc_m ) then !.and. om1 .gt. omthres ) then
                if ( blochlcorrections ) then
                    wts1=lo_LV_tetrahedron_weights(cval_m,om1,lo_freqtol,scsigma)
                    qpwm(b1,b2,b3,tetqpoints)=qpwm(b1,b2,b3,tetqpoints)+wts1*qp%tet(i)%weight
                else
                ! this minus event is relevant, store weights
                wts1=0.0_flyt
                wts2=0.0_flyt
                call lo_integration_weights_for_one_tetrahedron(qp%tet(i),cval_m,om1-deltaeps,wts1,scsigma,blochlcorrections)
                call lo_integration_weights_for_one_tetrahedron(qp%tet(i),cval_m,om1+deltaeps,wts2,scsigma,blochlcorrections)
                qpwm(b1,b2,b3,tetqpoints)=qpwm(b1,b2,b3,tetqpoints)+(wts1+wts2)*0.5_flyt
                endif
            endif
        enddo
        enddo
        enddo        
    enddo tetloop

    i_gamma=-1
    do i=1,qp%nq_tot
        if ( lo_sqnorm(qp%ap(i)%w) .lt. lo_sqtol ) then
            i_gamma=i
            exit
        endif
    enddo
    if ( i_gamma .lt. 0 ) call lo_stop_gracefully(['FFT mesh does not contain gamma'],lo_exitcode_symmetry,__FILE__,__LINE__)
    do b3=1,dr%nb
    do b2=1,dr%nb
    do b1=1,dr%nb
        if ( dr%aq(i_gamma)%omega(b1) .lt. lo_freqtol .or. &
             dr%aq(i_gamma)%omega(b2) .lt. lo_freqtol .or. &
             dr%aq(i_gamma)%omega(b3) .lt. lo_freqtol ) then
            qpwp(b1,b2,b3,i_gamma)=0.0_flyt
            qpwm(b1,b2,b3,i_gamma)=0.0_flyt
        endif
    enddo
    enddo
    enddo
    qpwp(1:3,1:3,1:3,i_gamma)=0.0_flyt
    qpwm(1:3,1:3,1:3,i_gamma)=0.0_flyt
    
    ! Store weights per q-point    
    do b1=1,dr%nb
    do b2=1,dr%nb
    do b3=1,dr%nb
        ! count q-points
        ii=0
        jj=0
        do i=1,qp%nq_tot
            if ( abs(qpwp(b1,b2,b3,i)) .gt. wthres ) then
                ii=ii+1
            endif
            if ( abs(qpwm(b1,b2,b3,i)) .gt. wthres ) then
                jj=jj+1
            endif
        enddo
        ! make some space
        scq%plus(b1,b2,b3)%n=ii
        scq%minus(b1,b2,b3)%n=jj
        if ( scq%plus(b1,b2,b3)%n .gt. 0 ) then
           lo_allocate(scq%plus(b1,b2,b3)%e( ii ))
        endif
        if ( scq%minus(b1,b2,b3)%n .gt. 0 ) then
           lo_allocate(scq%minus(b1,b2,b3)%e( jj ))
        endif
        ! store weights and indices
        ii=0
        jj=0
        do i=1,qp%nq_tot
            if ( abs(qpwp(b1,b2,b3,i)) .gt. wthres ) then
                ii=ii+1
                gi2=i
                scq%plus(b1,b2,b3)%e(ii)%gi2=gi2
                scq%plus(b1,b2,b3)%e(ii)%gi3=fft_third_grid_index(gi1,gi2,dims)
                scq%plus(b1,b2,b3)%e(ii)%deltafunction=qpwp(b1,b2,b3,i)
            endif
            if ( abs(qpwm(b1,b2,b3,i)) .gt. wthres ) then
                jj=jj+1
                gi2=i
                scq%minus(b1,b2,b3)%e(jj)%gi2=gi2
                scq%minus(b1,b2,b3)%e(jj)%gi3=fft_third_grid_index(gi1,gi2,dims)
                scq%minus(b1,b2,b3)%e(jj)%deltafunction=qpwm(b1,b2,b3,i)
            endif
        enddo
    enddo
    enddo
    enddo

    ! Cleanup    
    lo_deallocate(qpwp)
    lo_deallocate(qpwm)
    lo_deallocate(omr1)
    lo_deallocate(omr2)
    lo_deallocate(omr3)    
end subroutine