LCOV - code coverage report
Current view: top level - wannier/uhu - wann_uHu_radintsra5.F (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 24 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : c*************************************************c
       2             : c   compute radial integrals of sph. Hamiltonian  c
       3             : c   in scalar relativistic approximation between  c
       4             : c   products of radial function u and spherical   c
       5             : c   Bessel functions                              c
       6             : c                                                 c
       7             : c   integral = < uj1 | H_sra(lzb) | uj2 >         c
       8             : c*************************************************c
       9             : c                       J.-P. Hanke, Dec. 2015    c
      10             : c*************************************************c
      11             :       MODULE m_wann_uHu_radintsra5
      12             :       CONTAINS
      13             : 
      14           0 :       SUBROUTINE wann_uHu_radintsra5(jmtd,jri,rmsh,dx,e,vr,
      15           0 :      >                           uj1,uj2,duj1,duj2,lmaxd,lzb,
      16             :      >                           integral,irank)
      17             : 
      18             :       USE m_intgr, ONLY : intgr3
      19             :       USE m_constants
      20             :       IMPLICIT NONE
      21             : 
      22             :       REAL difcub
      23             :       EXTERNAL difcub
      24             :       
      25             :       INTEGER, INTENT(IN) :: jri,jmtd,lmaxd,irank
      26             :       REAL, INTENT(OUT) :: integral
      27             :       REAL, INTENT(IN) :: vr(jmtd),rmsh(jmtd),dx,e
      28             :       REAL, INTENT(IN) :: uj1(jmtd,2)      ! u(b1)*j
      29             :       REAL, INTENT(IN) :: uj2(jmtd,2)      ! u(b2)*j
      30             :       REAL, INTENT(IN) :: duj1(jmtd,2)     ! d/dr (u(b1)j)
      31             :       REAL, INTENT(IN) :: duj2(jmtd,2)     ! d/dr (u(b2)j)
      32             :       INTEGER, INTENT(IN) :: lzb         ! l of zentrifugal barrier
      33             : 
      34           0 :       REAL, ALLOCATABLE :: x(:)
      35             :       REAL :: c,c2,cin2,cin
      36             :       REAL :: ll,xi,vv,mm
      37             :       REAL :: sfp,facr,facl
      38           0 :       REAL, ALLOCATABLE :: xx(:,:),intt(:)
      39             :       INTEGER :: i,j,symopt
      40             : 
      41           0 :       c = c_light(1.)
      42           0 :       c2 = c*c
      43           0 :       cin = 1./c
      44           0 :       cin2 = cin*cin
      45           0 :       ll = lzb*(lzb+1)
      46           0 :       sfp = sqrt(4.0*pimach())
      47           0 :       symopt=0
      48             :       
      49             :       if(symopt.eq.0) then ! symmetrize
      50           0 :        facr = 0.5
      51           0 :        facl = 0.5
      52             :       elseif(symopt.lt.0)then ! apply H to left
      53             :        facr = 0.0
      54             :        facl = 1.0
      55             :       else ! apply H to right
      56             :        facr = 1.0
      57             :        facl = 0.0
      58             :       endif
      59             : 
      60           0 :       allocate( x(jri) )
      61           0 :       allocate( xx(jri,8),intt(8) )
      62             : 
      63             :       ! compute matrix elements of semi-relativistic
      64             :       ! Hamiltonian [Eq.(3.54) in PhD thesis of P.Kurz]
      65           0 :       DO i = 1, jri
      66           0 :          xi = rmsh(i)
      67           0 :          vv = vr(i) / xi !* sfp   ! no need to correct for Y_00
      68           0 :          mm = 1. + 0.5 * cin2 * ( e - vv )
      69             :          x(i) = 
      70             :      >        ! large-H-large
      71             :      >          uj1(i,1) * uj2(i,1) * ( 0.5 / mm * ll / xi / xi + vv )
      72             :      >        ! small-H-small
      73             :      >        + uj1(i,2) * uj2(i,2) * ( -2. * c2 + vv )
      74             : c      ! NEW VERSION BELOW
      75             : c     >        - c / xi * ( uj1(i,1)*uj2(i,2) + uj1(i,2)*uj2(i,1) )
      76             : c     >        + c / 2. * ( uj1(i,2)*duj2(i,1) + duj1(i,2)*uj2(i,1) )
      77             : c     >        - c / 2. * ( uj1(i,1)*duj2(i,2) + duj1(i,1)*uj2(i,2) )
      78             :       ! OLD VERSION BELOW
      79             :      >        ! large-H-small (symmetrized)
      80             :      >        - facr * c * uj1(i,1) * ( 2. * uj2(i,2) / xi + duj2(i,2) )
      81             :      >        - facl * c * uj2(i,1) * ( 2. * uj1(i,2) / xi + duj1(i,2) )
      82             :      >        ! small-H-large (symmetrized)
      83             :      >        + facr * c * uj1(i,2) * duj2(i,1)
      84           0 :      >        + facl * c * uj2(i,2) * duj1(i,1)
      85             : 
      86             : c         xx(i,:) = 0.0
      87             : c
      88             : c         xx(i,1)=uj1(i,1) * uj2(i,1) * ( 0.5 / mm * ll / xi / xi )
      89             : c         xx(i,2)=uj1(i,1) * uj2(i,1) * ( vv )
      90             : c
      91             : c         xx(i,3)=uj1(i,2) * uj2(i,2) * ( -2. * c2 )
      92             : c         xx(i,4)=uj1(i,2) * uj2(i,2) * ( vv )
      93             : c
      94             : c         xx(i,5)=-uj1(i,1) * c * ( 2.*uj2(i,2)/xi )
      95             : c         xx(i,6)=-uj1(i,1) * c * ( duj2(i,2) )
      96             : c
      97             : c         xx(i,7)= c * uj1(i,2) * duj2(i,1)
      98             : c
      99             : c         x(i) = xx(i,1)+xx(i,2)+xx(i,3)+xx(i,4)
     100             : c     >         +xx(i,5)+xx(i,6)+xx(i,7)+xx(i,8)
     101             : 
     102             :       ENDDO
     103           0 :       call intgr3(x,rmsh,dx,jri,integral)
     104             : c      do i=1,8
     105             : c       call intgr3(xx(:,i),rmsh,dx,jri,intt(i))
     106             : c      enddo
     107             :       !integral = integral / sfp
     108             : 
     109             : c      if(irank.eq.0 .and. abs(integral).gt.1e-10) then
     110             : c       do i=1,8
     111             : c        write(*,'(a,i1,a,f8.2)')'int',i,':',intt(i)/integral*100.
     112             : c       enddo
     113             : c       write(*,*)'integral:',integral
     114             : c      endif
     115             : 
     116           0 :       deallocate( x )
     117           0 :       deallocate( xx,intt )
     118             : 
     119           0 :       END SUBROUTINE wann_uHu_radintsra5
     120             :       END MODULE m_wann_uHu_radintsra5

Generated by: LCOV version 1.13