Line data    Source code

       1             :       MODULE m_wann_uHu_radintsra
       2             :       CONTAINS
       3           0 :       SUBROUTINE wann_uHu_radintsra(jmtd,jri,rmsh,dx,
       4           0 :      >                          epar,vr,f,g,l,expect)
       5             : 
       6             :       USE m_intgr, ONLY : intgr3
       7             :       USE m_constants
       8             :       USE m_difcub
       9             : 
      10             :       IMPLICIT NONE
      11             :       
      12             :       INTEGER, INTENT(IN) :: jmtd
      13             :       REAL, INTENT(IN) :: f(jmtd,2),g(jmtd,2)
      14             :       REAL, INTENT(IN) :: vr(jmtd)
      15             :       REAL, INTENT(IN) :: epar
      16             :       REAL, INTENT(IN) :: rmsh(jmtd)
      17             :       REAL, INTENT(IN) :: dx
      18             :       REAL, INTENT(IN) :: expect
      19             :       INTEGER, INTENT(IN) :: jri
      20             :       INTEGER, INTENT(IN) :: l
      21             : 
      22           0 :       REAL, ALLOCATABLE :: x(:),dg(:,:),t(:,:),vv(:)
      23             :       REAL :: t11,t22,t12,t21,total,norm
      24             :       REAL :: mm,c,c2,cin2,cin
      25             :       REAL :: ll,xi,sfp
      26             :       INTEGER :: i,j
      27             : 
      28           0 :       c = c_light(1.)
      29           0 :       c2 = c*c
      30           0 :       cin = 1./c
      31           0 :       cin2 = cin*cin
      32           0 :       ll = l*(l+1)
      33             :       sfp = sqrt(4.0*pimach())
      34             : 
      35           0 :       allocate( x(jri), dg(jri,2), t(jri,2), vv(jri) )
      36             : 
      37             :       ! derivatives d/dr g for large and small component
      38           0 :       DO i=1,jri
      39           0 :        t(i,:) = g(i,:)/rmsh(i)
      40             :       ENDDO
      41             : 
      42           0 :       DO j = 1, 2
      43             :          ! derivative at 1st point
      44           0 :          dg(1,j) = difcub( rmsh(1),t(1,j),rmsh(1) )
      45             : 
      46             :          ! derivative at 2nd...(jri-2)th point
      47           0 :          DO i = 2, jri-2
      48           0 :             dg(i,j) = difcub( rmsh(i-1),t(i-1,j),rmsh(i) )
      49             :          ENDDO
      50             : 
      51             :          ! derivative at last two points
      52           0 :          dg(jri-1,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri-1) )
      53           0 :          dg(jri,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri) )
      54             :       ENDDO
      55             : 
      56           0 :       DO i=1,jri
      57           0 :        xi = rmsh(i)
      58           0 :        dg(i,:) = dg(i,:) * xi
      59           0 :        vv(i) = vr(i) / xi !* sfp
      60             :       ENDDO
      61             : 
      62             :       ! check normalization
      63           0 :       DO i = 1, jri
      64           0 :          x(i) = f(i,1)*g(i,1)+f(i,2)*g(i,2)
      65             :       ENDDO
      66           0 :       call intgr3(x,rmsh,dx,jri,norm)
      67           0 :       write(*,*)'norm:',norm
      68             : 
      69             : 
      70             :       ! compute matrix elements of semi-relativistic
      71             :       ! Hamiltonian [Eq.(3.54) in PhD thesis of P.Kurz]
      72           0 :       DO i = 1, jri
      73           0 :          mm = 1. + 0.5 * cin2 * ( epar - vv(i) )
      74             :          x(i) = f(i,1) * g(i,1)
      75           0 :      >        * ( 0.5 / mm * ll / rmsh(i) / rmsh(i) + vv(i) )
      76             :       ENDDO
      77           0 :       call intgr3(x,rmsh,dx,jri,t11) ! large-H-large
      78             : 
      79           0 :       DO i = 1, jri
      80           0 :          x(i) = f(i,2) * g(i,2) * ( -2. * c2 + vv(i) )
      81             :       ENDDO
      82           0 :       call intgr3(x,rmsh,dx,jri,t22) ! small-H-small
      83             : 
      84           0 :       DO i = 1, jri
      85           0 :          x(i) = f(i,1) * ( 2. * g(i,2) / rmsh(i) + dg(i,2) )
      86             :       ENDDO
      87           0 :       x = -c * x
      88           0 :       call intgr3(x,rmsh,dx,jri,t12) ! large-H-small
      89             : 
      90           0 :       DO i = 1, jri
      91           0 :          x(i) = f(i,2) * dg(i,1)
      92             :       ENDDO
      93           0 :       x = c * x
      94           0 :       call intgr3(x,rmsh,dx,jri,t21) ! small-H-large
      95             : 
      96           0 :       total = t11 + t22 + t12 + t21
      97           0 :       write(*,'(a,f16.12)')'expect:',expect
      98           0 :       write(*,'(a,f16.12)')'integr:',total
      99           0 :       if(abs(expect).gt.1e-12) then
     100           0 :        write(*,'(a,f7.3,a)')'differ:',abs(expect-total)/expect*100.,' %'
     101             :       endif
     102             : 
     103           0 :       deallocate( x, dg, t, vv )
     104             : 
     105           0 :       END SUBROUTINE wann_uHu_radintsra
     106             :       END MODULE m_wann_uHu_radintsra