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 : sfp = sqrt(4.0*pimach())
47 0 : symopt=0
48 :
49 : if(symopt.eq.0) then ! symmetrize
50 : facr = 0.5
51 : 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
|