Line data    Source code

       1             :       MODULE m_wann_uHu_radintsra2
       2             :       CONTAINS
       3             : 
       4             :       ! < p | H(l) | q >
       5           0 :       SUBROUTINE wann_uHu_radintsra2(jmtd,jri,rmsh,dx,
       6           0 :      >                           e,vr,
       7           0 :      >                           p,q,l,
       8             :      >                           integral)
       9             : 
      10             :       USE m_intgr, ONLY : intgr3
      11             :       USE m_constants
      12             :       USE m_difcub
      13             : 
      14             :       IMPLICIT NONE
      15             :       
      16             :       INTEGER, INTENT(IN) :: jmtd
      17             :       REAL, INTENT(IN) :: p(jmtd,2)
      18             :       REAL, INTENT(IN) :: q(jmtd,2)
      19             :       REAL, INTENT(IN) :: vr(jmtd)
      20             :       REAL, INTENT(IN) :: e
      21             :       REAL, INTENT(IN) :: rmsh(jmtd)
      22             :       REAL, INTENT(IN) :: dx
      23             :       INTEGER, INTENT(IN) :: jri
      24             :       INTEGER, INTENT(IN) :: l
      25             :       REAL, INTENT(OUT) :: integral
      26             : 
      27           0 :       REAL, ALLOCATABLE :: x(:),dq(:,:),t(:,:)
      28             :       REAL :: c,c2,cin2,cin
      29             :       REAL :: ll,xi,vv,mm,sfp
      30             :       INTEGER :: i,j
      31             : 
      32           0 :       c = c_light(1.)
      33           0 :       c2 = c*c
      34           0 :       cin = 1./c
      35           0 :       cin2 = cin*cin
      36           0 :       ll = l*(l+1)
      37             :       sfp = sqrt(4.0*pimach())
      38             : 
      39           0 :       allocate( x(jri), dq(jri,2), t(jri,2) )
      40             : 
      41           0 :       DO i=1,jri
      42           0 :        t(i,:) = q(i,:) / rmsh(i)
      43             :       ENDDO
      44             : 
      45             :       ! derivatives d/dr for large and small component
      46           0 :       DO j = 1, 2
      47             :          ! derivative at 1st point
      48           0 :          dq(1,j) = difcub( rmsh(1),t(1,j),rmsh(1) )
      49             : 
      50             :          ! derivative at 2nd...(jri-2)th point
      51           0 :          DO i = 2, jri-2
      52           0 :             dq(i,j) = difcub( rmsh(i-1),t(i-1,j),rmsh(i) )
      53             :          ENDDO
      54             : 
      55             :          ! derivative at last two points
      56           0 :          dq(jri-1,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri-1) )
      57           0 :          dq(jri,j) = difcub( rmsh(jri-3),t(jri-3,j),rmsh(jri) )
      58             :       ENDDO
      59             : 
      60           0 :       DO i=1,jri
      61           0 :        dq(i,:) = dq(i,:)*rmsh(i)
      62             :       ENDDO
      63             : 
      64             :       ! compute matrix elements of semi-relativistic
      65             :       ! Hamiltonian [Eq.(3.54) in PhD thesis of P.Kurz]
      66           0 :       DO i = 1, jri
      67           0 :          xi = rmsh(i)
      68           0 :          vv = vr(i) / xi !* sfp
      69           0 :          mm = 1. + 0.5 * cin2 * ( e - vv )
      70             :          x(i) = 
      71             :      >        ! large-H-large
      72             :      >          p(i,1) * q(i,1) * ( 0.5 / mm * ll / xi / xi + vv )
      73             :      >        ! small-H-small
      74             :      >        + p(i,2) * q(i,2) * ( -2. * c2 + vv )
      75             :      >        ! large-H-small (not symmetrized)
      76             :      >        - c * p(i,1) * (2. * q(i,2) / xi + dq(i,2) )
      77             :      >        ! small-H-large (not symmetrized)
      78           0 :      >        + c * p(i,2) * dq(i,1)
      79             :       ENDDO
      80           0 :       call intgr3(x,rmsh,dx,jri,integral)
      81             : 
      82           0 :       deallocate( x, dq, t )
      83             : 
      84           0 :       END SUBROUTINE wann_uHu_radintsra2
      85             :       END MODULE m_wann_uHu_radintsra2