Subroutine U_matrix_2(nbasis,nenergy,U,icase,anglec,angles,r,a,b)
IMPLICIT NONE
INTEGER nbasis,nenergy,icase
real*8 U(nenergy,nenergy),Rmatrix(nenergy,nenergy),npower
real*8 Pmatrix(nenergy,nenergy),Gmatrix(nenergy,nenergy)
real*8 Cmatrix(nenergy,nenergy),Mtemp(nenergy,nenergy),Tmatrix(nenergy,nenergy)
INTEGER i,j,k
real*8 angle,anglec,angles,temp,r,a(8),b(8),phi,phic,phis,phi2
real*8 anglec2,angles2,angleM,temp2,temp3
REAL*8 CTEMP1(4,4),Ctemp2(4,4),ctemp(4,4)
LOGICAL DEBUG
data DEBUG /.FALSE./
data npower /2d0/
!--------------------------------------------------------------
!define the basis
! |1> =ss
! |2> =sx
! |3> =sy
! |4> =sz
! |5> =xs
! |6> =ys
! |7> =zs
! |8> =xx
! |9> =xy
! |10> =xz
! |11> =yx
! |12> =yy
! |13> =yz
! |14> =zx
! |15> =zy
! |16> =zz
!--------------------------------------------------------------
!--------------------------------------------------------------
!define the energy for XY diatom
!E1=3sigma_u^+
!E2=3sigma_g^+(sz case)
!E3=2 3sigma_u^+(sz case)
!E4=E5=3pi_g
!E6=E7=3pi_u
!E8=3 3sigma_u^+
!E9=3sigma_g^-
!E10=E11=2 3pi_g
!E12=4 3sigma_u^+
!E13=E14=2 3pi_u
!E15=E16=3delta_u
!--------------------------------------------------------------
!U=R*P*T*G*C
!in the triplet Li case
!--------------------------------------------------------------
IF(icase.eq.1)THEN
!--------------------------------------------------------------
!we do the AB case
!--------------------------------------------------------------
anglec2=anglec*anglec
angles2=angles*angles
angleM=angles*anglec
Rmatrix=0.0d0
Rmatrix(1,1)=1.0d0
Rmatrix(2,2)=anglec
Rmatrix(2,4)=angles
Rmatrix(3,3)=1.0d0
Rmatrix(4,2)=-angles
Rmatrix(4,4)=anglec
Rmatrix(5,5)=anglec
Rmatrix(5,7)=angles
Rmatrix(6,6)=1.0d0
Rmatrix(7,5)=-angles
Rmatrix(7,7)=anglec
Rmatrix(8,8)=anglec2
Rmatrix(8,10)=angleM
Rmatrix(8,14)=angleM
Rmatrix(8,16)=angles2
Rmatrix(9,9)=anglec
Rmatrix(9,15)=angles
Rmatrix(10,8)=-angleM
Rmatrix(10,10)=anglec2
Rmatrix(10,14)=-angles2
Rmatrix(10,16)=angleM
Rmatrix(11,11)=anglec
Rmatrix(11,13)=angles
Rmatrix(12,12)=1.0d0
Rmatrix(13,11)=-angles
Rmatrix(13,13)=anglec
Rmatrix(14,8)=-angleM
Rmatrix(14,10)=-angles2
Rmatrix(14,14)=anglec2
Rmatrix(14,16)=angleM
Rmatrix(15,9)=-angles
Rmatrix(15,15)=anglec
Rmatrix(16,8)=angles2
Rmatrix(16,10)=-angleM
Rmatrix(16,14)=-angleM
Rmatrix(16,16)=anglec2



u=Rmatrix
!CHECK
!write(154,*)'Rmatrix'
298 FORMAT(E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x)
!----------------------------------
! test for Umatrix
!---------------------------------
IF(DEBUG)THEN
DO i=1,nenergy
   DO j=1,nenergy
      temp=0.0d0
      DO k=1,nenergy
      temp=temp+U(k,i)*U(k,j)
      ENDDO
      IF(i.eq.j)THEN
           IF(abs(temp-1.0d0).gt.1.0d-8)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 1'
           ENDIF
      ELSE
           IF(abs(temp).gt.1.0d-8)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 1'
           ENDIF
      ENDIF
   ENDDO
ENDDO
ENDIF
!---------------------------------------------------------------
Pmatrix=0.0d0
temp=1.0d0/sqrt(2.0d0)
Pmatrix(1,1)=1.d0
Pmatrix(2,2)=-temp
Pmatrix(2,4)=temp
Pmatrix(3,2)=-temp
Pmatrix(3,4)=-temp
Pmatrix(4,3)=1.0d0
Pmatrix(5,5)=-temp
Pmatrix(5,7)=temp
Pmatrix(6,5)=-temp
Pmatrix(6,7)=-temp
Pmatrix(7,6)=1.0d0

Pmatrix(8,8)=0.5d0
Pmatrix(8,10)=-0.5d0
Pmatrix(8,14)=-0.5d0
Pmatrix(8,16)=0.5d0
Pmatrix(9,8)=0.5d0
Pmatrix(9,10)=0.5d0
Pmatrix(9,14)=-0.5d0
Pmatrix(9,16)=-0.5d0
Pmatrix(10,9)=-temp
Pmatrix(10,15)=temp
Pmatrix(11,8)=0.5d0
Pmatrix(11,10)=-0.5d0
Pmatrix(11,14)=0.5d0
Pmatrix(11,16)=-0.5d0
Pmatrix(12,8)=0.5d0
Pmatrix(12,10)=0.5d0
Pmatrix(12,14)=0.5d0
Pmatrix(12,16)=0.5d0
Pmatrix(13,9)=-temp
Pmatrix(13,15)=-temp
Pmatrix(14,11)=-temp
Pmatrix(14,13)=temp
Pmatrix(15,11)=-temp
Pmatrix(15,13)=-temp
Pmatrix(16,12)=1.0d0


u=matmul(Rmatrix,Pmatrix)
MTEMP=u
!write(154,*)'Pmatrix'
!DO i=1,7
!write(154,298)(Pmatrix(i,j),j=1,7)
!ENDDO
!----------------------------------
! test for Umatrix
!---------------------------------
IF(DEBUG)THEn
DO i=1,nenergy
   DO j=1,nenergy
      temp=0.0d0
      DO k=1,nenergy
      temp=temp+U(k,i)*U(k,j)
      ENDDO
      IF(i.eq.j)THEN
           IF(abs(temp-1.0d0).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 2'
           ENDIF
      ELSE
           IF(abs(temp).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 2'
           ENDIF
      ENDIF
   ENDDO
ENDDO
ENDIF

!--------------------------------------------
!T matrix
Tmatrix=0.0d0
DO i=1,7
Tmatrix(i,i)=1.0d0
ENDDO
temp=1.0d0/sqrt(2.0d0)
temp2=1.0d0/sqrt(3.0d0)
temp3=1.0d0/sqrt(6.0d0)

Tmatrix(8,15)=1.0d0
Tmatrix(9,10)=temp
Tmatrix(9,13)=temp
Tmatrix(10,8)=temp2
Tmatrix(10,9)=temp
Tmatrix(10,12)=temp3
Tmatrix(11,10)=-temp
Tmatrix(11,13)=temp
Tmatrix(12,8)=-temp2
Tmatrix(12,12)=sqrt(2d0/3d0)
Tmatrix(13,11)=temp
Tmatrix(13,14)=temp
Tmatrix(14,8)=temp2
Tmatrix(14,9)=-temp
Tmatrix(14,12)=temp3
Tmatrix(15,11)=-temp
Tmatrix(15,14)=temp
Tmatrix(16,16)=1.0d0

u=matmul(MTEMP,Tmatrix)
Mtemp=u
!---------------------------------------------------------------
Gmatrix=0.0d0
Gmatrix(1,1)=1.d0
Gmatrix(2,4)=temp
Gmatrix(2,6)=-temp
Gmatrix(3,2)=temp
Gmatrix(3,3)=-temp
Gmatrix(4,5)=temp
Gmatrix(4,7)=-temp
Gmatrix(5,4)=temp
Gmatrix(5,6)=temp
Gmatrix(6,2)=temp
Gmatrix(6,3)=temp
Gmatrix(7,5)=temp
Gmatrix(7,7)=temp

DO i=8,nenergy
Gmatrix(i,i)=1.0d0
ENDDO


u=matmul(Mtemp,Gmatrix)
Mtemp=u
!write(154,*)'Gmatrix'
!DO i=1,7
!write(154,298)(Gmatrix(i,j),j=1,7)
!ENDDO
!----------------------------------
! test for Umatrix
!---------------------------------
IF(DEBUG)THEN
DO i=1,nenergy
   DO j=1,nenergy
      temp=0.0d0
      DO k=1,nenergy
      temp=temp+U(k,i)*U(k,j)
      ENDDO
      IF(i.eq.j)THEN
           IF(abs(temp-1.0d0).gt.1.0d-8)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 3'
           ENDIF
      ELSE
           IF(abs(temp).gt.1.0d-8)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 3'
           ENDIF
      ENDIF
   ENDDO
ENDDO
ENDIF

!--------------------------------------------------------------


!1,3 states are both 3sigma_u^+ so they mix
IF(abs(a(1)).lt.1d-10.or.abs(exp(b(1)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(1)*r**(npower))*exp(b(1)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP1=0.0d0
CTEMP1(1,1)=phic
CTEMP1(1,2)=phis
CTEMP1(2,1)=-phis
CTEMP1(2,2)=phic
CTEMP1(3,3)=1.0d0
CTEMP1(4,4)=1.0d0

IF(abs(a(2)).lt.1d-10.or.abs(exp(b(2)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(2)*r**(npower))*exp(b(2)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(1,1)=phic
CTEMP2(1,3)=phis
CTEMP2(3,1)=-phis
CTEMP2(3,3)=phic
CTEMP2(2,2)=1.0d0
CTEMP2(4,4)=1.0d0


CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp
IF(abs(a(3)).lt.1d-10.or.abs(exp(b(3)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(3)*r**(npower))*exp(b(3)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(1,1)=phic
CTEMP2(1,4)=phis
CTEMP2(4,1)=-phis
CTEMP2(4,4)=phic
CTEMP2(2,2)=1.0d0
CTEMP2(3,3)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(4)).lt.1d-10.or.abs(exp(b(4)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(4)*r**(npower))*exp(b(4)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(2,2)=phic
CTEMP2(2,3)=phis
CTEMP2(3,2)=-phis
CTEMP2(3,3)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(4,4)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(5)).lt.1d-10.or.abs(exp(b(5)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(5)*r**(npower))*exp(b(5)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(2,2)=phic
CTEMP2(2,4)=phis
CTEMP2(4,2)=-phis
CTEMP2(4,4)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(3,3)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(6)).lt.1d-10.or.abs(exp(b(6)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(6)*r**(npower))*exp(b(6)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(3,3)=phic
CTEMP2(3,4)=phis
CTEMP2(4,3)=-phis
CTEMP2(4,4)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(2,2)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)



CMATRIX=0.0d0
CMATRIX(1,1)=CTEMP(1,1)
CMATRIX(1,3)=CTEMP(1,2)
CMATRIX(1,8)=CTEMP(1,3)
CMATRIX(1,12)=CTEMP(1,4)
CMATRIX(3,1)=CTEMP(2,1)
CMATRIX(3,3)=CTEMP(2,2)
CMATRIX(3,8)=CTEMP(2,3)
CMATRIX(3,12)=CTEMP(2,4)
CMATRIX(8,1)=CTEMP(3,1)
CMATRIX(8,3)=CTEMP(3,2)
CMATRIX(8,8)=CTEMP(3,3)
CMATRIX(8,12)=CTEMP(3,4)
CMATRIX(12,1)=CTEMP(4,1)
CMATRIX(12,3)=CTEMP(4,2)
CMATRIX(12,8)=CTEMP(4,3)
CMATRIX(12,12)=CTEMP(4,4)

!Pi_g coupling
IF(abs(a(7)).lt.1d-10.or.abs(exp(b(7)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(7)*r**(npower))*exp(b(7)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CMATRIX(4,4)=phic
CMATRIX(4,10)=phis
CMATRIX(5,5)=phic
CMATRIX(5,11)=-phis
CMATRIX(10,4)=-phis
CMATRIX(10,10)=phic
CMATRIX(11,5)=phis
CMATRIX(11,11)=phic

!Pi_u coupling
IF(abs(a(8)).lt.1d-10.or.abs(exp(b(8)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(8)*r**(npower))*exp(b(8)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF


CMATRIX(6,6)=phic
CMATRIX(6,13)=phis
CMATRIX(13,6)=-phis
CMATRIX(13,13)=phic
CMATRIX(7,7)=phic
CMATRIX(7,14)=-phis
CMATRIX(14,7)=phis
CMATRIX(14,14)=phic

CMATRIX(2,2)=1.0d0
CMATRIX(9,9)=1.0d0
CMATRIX(15,15)=1.0d0
CMATRIX(16,16)=1.0d0

IF(DEBUG)THEN
DO i=1,nenergy
DO j=1,nenergy
!IF(i.eq.j.and.abs(CMATRIX(i,j)-1.0d0).gt.1d-8)THEN
!print*,'I,I=',i,j,CMATRIX(i,j)
!ELSEIF(i.ne.j.and.abs(CMATRIX(i,j).gt.1d-8))THEN
!print*,'I,J=',i,j,CMATRIX(i,j)
!ENDIF
ENDDO
ENDDO
ENDIF

u=matmul(Mtemp,Cmatrix)
!write(154,*)'Umatrix'
!DO i=1,7
!write(154,298)(U(i,j),j=1,7)
!ENDDO
!stop
!print*,'U=',u
!print*,'before test1'
!----------------------------------
! test for Umatrix
!---------------------------------
IF(DEBUG)THEN
DO i=1,nenergy
   DO j=1,nenergy
      temp=0.0d0
      DO k=1,nenergy
      temp=temp+U(k,i)*U(k,j)
      ENDDO
!      print*,'i,j',i,j,temp
      IF(i.eq.j)THEN
           IF(abs(temp-1.0d0).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 4'
           ENDIF
      ELSE
           IF(abs(temp).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 4'
           ENDIF
      ENDIF
   ENDDO
ENDDO
print*,'Umatrix check is done'
!stop
ENDIF

ELSEIF(icase.eq.2)THEN
!--------------------------------------------------------------
!we do the BC case
!--------------------------------------------------------------
Pmatrix=0.0d0
temp=1.0d0/sqrt(2.0d0)
Pmatrix(1,1)=1.d0
Pmatrix(2,2)=-temp
Pmatrix(2,4)=temp
Pmatrix(3,2)=-temp
Pmatrix(3,4)=-temp
Pmatrix(4,3)=1.0d0
Pmatrix(5,5)=-temp
Pmatrix(5,7)=temp
Pmatrix(6,5)=-temp
Pmatrix(6,7)=-temp
Pmatrix(7,6)=1.0d0
Pmatrix(8,8)=0.5d0
Pmatrix(8,10)=-0.5d0
Pmatrix(8,14)=-0.5d0
Pmatrix(8,16)=0.5d0
Pmatrix(9,8)=0.5d0
Pmatrix(9,10)=0.5d0
Pmatrix(9,14)=-0.5d0
Pmatrix(9,16)=-0.5d0
Pmatrix(10,9)=-temp
Pmatrix(10,15)=temp
Pmatrix(11,8)=0.5d0
Pmatrix(11,10)=-0.5d0
Pmatrix(11,14)=0.5d0
Pmatrix(11,16)=-0.5d0
Pmatrix(12,8)=0.5d0
Pmatrix(12,10)=0.5d0
Pmatrix(12,14)=0.5d0
Pmatrix(12,16)=0.5d0
Pmatrix(13,9)=-temp
Pmatrix(13,15)=-temp
Pmatrix(14,11)=-temp
Pmatrix(14,13)=temp
Pmatrix(15,11)=-temp
Pmatrix(15,13)=-temp
Pmatrix(16,12)=1.0d0
!--------------------------------------------
!T matrix
Tmatrix=0.0d0
DO i=1,7
Tmatrix(i,i)=1.0d0
ENDDO
temp=1.0d0/sqrt(2.0d0)
temp2=1.0d0/sqrt(3.0d0)
temp3=1.0d0/sqrt(6.0d0)

Tmatrix(8,15)=1.0d0
Tmatrix(9,10)=temp
Tmatrix(9,13)=temp
Tmatrix(10,8)=temp2
Tmatrix(10,9)=temp
Tmatrix(10,12)=temp3
Tmatrix(11,10)=-temp
Tmatrix(11,13)=temp
Tmatrix(12,8)=-temp2
Tmatrix(12,12)=sqrt(2d0/3d0)
Tmatrix(13,11)=temp
Tmatrix(13,14)=temp
Tmatrix(14,8)=temp2
Tmatrix(14,9)=-temp
Tmatrix(14,12)=temp3
Tmatrix(15,11)=-temp
Tmatrix(15,14)=temp
Tmatrix(16,16)=1.0d0

u=matmul(PMATRIX,Tmatrix)
Mtemp=u





!---------------------------------------------------------------
Gmatrix=0.0d0
Gmatrix(1,1)=1.d0
Gmatrix(2,4)=temp
Gmatrix(2,6)=-temp
Gmatrix(3,2)=temp
Gmatrix(3,3)=-temp
Gmatrix(4,5)=temp
Gmatrix(4,7)=-temp
Gmatrix(5,4)=temp
Gmatrix(5,6)=temp
Gmatrix(6,2)=temp
Gmatrix(6,3)=temp
Gmatrix(7,5)=temp
Gmatrix(7,7)=temp
DO i=8,nenergy
Gmatrix(i,i)=1.0d0
ENDDO

u=matmul(Mtemp,Gmatrix)
Mtemp=u
!--------------------------------------------------------------
Cmatrix=0.0d0
!1,3 states are both 3sigma_u^+ so they mix
IF(abs(a(1)).lt.1d-10.or.abs(exp(b(1)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(1)*r**(npower))*exp(b(1)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP1=0.0d0
CTEMP1(1,1)=phic
CTEMP1(1,2)=phis
CTEMP1(2,1)=-phis
CTEMP1(2,2)=phic
CTEMP1(3,3)=1.0d0
CTEMP1(4,4)=1.0d0

IF(abs(a(2)).lt.1d-10.or.abs(exp(b(2)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(2)*r**(npower))*exp(b(2)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(1,1)=phic
CTEMP2(1,3)=phis
CTEMP2(3,1)=-phis
CTEMP2(3,3)=phic
CTEMP2(2,2)=1.0d0
CTEMP2(4,4)=1.0d0


CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp
IF(abs(a(3)).lt.1d-10.or.abs(exp(b(3)*r)).gt.1d20)THEN

phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(3)*r**(npower))*exp(b(3)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(1,1)=phic
CTEMP2(1,4)=phis
CTEMP2(4,1)=-phis
CTEMP2(4,4)=phic
CTEMP2(2,2)=1.0d0
CTEMP2(3,3)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(4)).lt.1d-10.or.abs(exp(b(4)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(4)*r**(npower))*exp(b(4)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(2,2)=phic
CTEMP2(2,3)=phis
CTEMP2(3,2)=-phis
CTEMP2(3,3)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(4,4)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(5)).lt.1d-10.or.abs(exp(b(5)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(5)*r**(npower))*exp(b(5)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(2,2)=phic
CTEMP2(2,4)=phis
CTEMP2(4,2)=-phis
CTEMP2(4,4)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(3,3)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(6)).lt.1d-10.or.abs(exp(b(6)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(6)*r**(npower))*exp(b(6)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(3,3)=phic
CTEMP2(3,4)=phis
CTEMP2(4,3)=-phis
CTEMP2(4,4)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(2,2)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)

CMATRIX=0.0d0
CMATRIX(1,1)=CTEMP(1,1)
CMATRIX(1,3)=CTEMP(1,2)
CMATRIX(1,8)=CTEMP(1,3)
CMATRIX(1,12)=CTEMP(1,4)
CMATRIX(3,1)=CTEMP(2,1)
CMATRIX(3,3)=CTEMP(2,2)
CMATRIX(3,8)=CTEMP(2,3)
CMATRIX(3,12)=CTEMP(2,4)
CMATRIX(8,1)=CTEMP(3,1)
CMATRIX(8,3)=CTEMP(3,2)
CMATRIX(8,8)=CTEMP(3,3)
CMATRIX(8,12)=CTEMP(3,4)
CMATRIX(12,1)=CTEMP(4,1)
CMATRIX(12,3)=CTEMP(4,2)
CMATRIX(12,8)=CTEMP(4,3)
CMATRIX(12,12)=CTEMP(4,4)

!Pi_g coupling
IF(abs(a(7)).lt.1d-10.or.abs(exp(b(7)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(7)*r**(npower))*exp(b(7)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF

CMATRIX(4,4)=phic
CMATRIX(4,10)=phis
CMATRIX(5,5)=phic
CMATRIX(5,11)=-phis
CMATRIX(10,4)=-phis
CMATRIX(10,10)=phic
CMATRIX(11,5)=phis
CMATRIX(11,11)=phic

!Pi_u coupling
IF(abs(a(8)).lt.1d-10.or.abs(exp(b(8)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(8)*r**(npower))*exp(b(8)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF

CMATRIX(6,6)=phic
CMATRIX(6,13)=phis
CMATRIX(13,6)=-phis
CMATRIX(13,13)=phic
CMATRIX(7,7)=phic
CMATRIX(7,14)=-phis
CMATRIX(14,7)=phis
CMATRIX(14,14)=phic

CMATRIX(2,2)=1.0d0
CMATRIX(9,9)=1.0d0
CMATRIX(15,15)=1.0d0
CMATRIX(16,16)=1.0d0

u=matmul(Mtemp,Cmatrix)

!----------------------------------
! test for Umatrix
!---------------------------------
IF(DEBUG)THEN
DO i=1,nenergy
   DO j=1,nenergy
      temp=0.0d0
      DO k=1,nenergy
      temp=temp+U(k,i)*U(k,j)
      ENDDO
!      print*,'i,j',i,j,temp
      IF(i.eq.j)THEN
           IF(abs(temp-1.0d0).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 2_4'
           ENDIF
      ELSE
           IF(abs(temp).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 2_4'
           ENDIF
      ENDIF
   ENDDO
ENDDO
print*,'Umatrix check is done'
!stop
ENDIF

ELSEIF(icase.eq.3)THEN
!--------------------------------------------------------------
!we do the AC case
!--------------------------------------------------------------
anglec2=anglec*anglec
angles2=angles*angles
angleM=angles*anglec

Rmatrix=0.0d0
Rmatrix(1,1)=1.0d0
Rmatrix(2,2)=anglec
Rmatrix(2,4)=angles
Rmatrix(3,3)=1.0d0
Rmatrix(4,2)=-angles
Rmatrix(4,4)=anglec
Rmatrix(5,5)=anglec
Rmatrix(5,7)=angles
Rmatrix(6,6)=1.0d0
Rmatrix(7,5)=-angles
Rmatrix(7,7)=anglec
Rmatrix(8,8)=anglec2
Rmatrix(8,10)=angleM
Rmatrix(8,14)=angleM
Rmatrix(8,16)=angles2
Rmatrix(9,9)=anglec
Rmatrix(9,15)=angles
Rmatrix(10,8)=-angleM
Rmatrix(10,10)=anglec2
Rmatrix(10,14)=-angles2
Rmatrix(10,16)=angleM
Rmatrix(11,11)=anglec
Rmatrix(11,13)=angles
Rmatrix(12,12)=1.0d0
Rmatrix(13,11)=-angles
Rmatrix(13,13)=anglec
Rmatrix(14,8)=-angleM
Rmatrix(14,10)=-angles2
Rmatrix(14,14)=anglec2
Rmatrix(14,16)=angleM
Rmatrix(15,9)=-angles
Rmatrix(15,15)=anglec
Rmatrix(16,8)=angles2
Rmatrix(16,10)=-angleM
Rmatrix(16,14)=-angleM
Rmatrix(16,16)=anglec2

!---------------------------------------------------------------
Pmatrix=0.0d0
temp=1.0d0/sqrt(2.0d0)
Pmatrix(1,1)=1.d0
Pmatrix(2,2)=-temp
Pmatrix(2,4)=temp
Pmatrix(3,2)=-temp
Pmatrix(3,4)=-temp
Pmatrix(4,3)=1.0d0
Pmatrix(5,5)=-temp
Pmatrix(5,7)=temp
Pmatrix(6,5)=-temp
Pmatrix(6,7)=-temp
Pmatrix(7,6)=1.0d0
Pmatrix(8,8)=0.5d0
Pmatrix(8,10)=-0.5d0
Pmatrix(8,14)=-0.5d0
Pmatrix(8,16)=0.5d0
Pmatrix(9,8)=0.5d0
Pmatrix(9,10)=0.5d0
Pmatrix(9,14)=-0.5d0
Pmatrix(9,16)=-0.5d0
Pmatrix(10,9)=-temp
Pmatrix(10,15)=temp
Pmatrix(11,8)=0.5d0
Pmatrix(11,10)=-0.5d0
Pmatrix(11,14)=0.5d0
Pmatrix(11,16)=-0.5d0
Pmatrix(12,8)=0.5d0
Pmatrix(12,10)=0.5d0
Pmatrix(12,14)=0.5d0
Pmatrix(12,16)=0.5d0
Pmatrix(13,9)=-temp
Pmatrix(13,15)=-temp
Pmatrix(14,11)=-temp
Pmatrix(14,13)=temp
Pmatrix(15,11)=-temp
Pmatrix(15,13)=-temp
Pmatrix(16,12)=1.0d0
u=matmul(Rmatrix,Pmatrix)
Mtemp=u
!---------------------------------------------------------------
!--------------------------------------------
!T matrix
Tmatrix=0.0d0
DO i=1,7
Tmatrix(i,i)=1.0d0
ENDDO
temp=1.0d0/sqrt(2.0d0)
temp2=1.0d0/sqrt(3.0d0)
temp3=1.0d0/sqrt(6.0d0)

Tmatrix(8,15)=1.0d0
Tmatrix(9,10)=temp
Tmatrix(9,13)=temp
Tmatrix(10,8)=temp2
Tmatrix(10,9)=temp
Tmatrix(10,12)=temp3
Tmatrix(11,10)=-temp
Tmatrix(11,13)=temp
Tmatrix(12,8)=-temp2
Tmatrix(12,12)=sqrt(2d0/3d0)
Tmatrix(13,11)=temp
Tmatrix(13,14)=temp
Tmatrix(14,8)=temp2
Tmatrix(14,9)=-temp
Tmatrix(14,12)=temp3
Tmatrix(15,11)=-temp
Tmatrix(15,14)=temp
Tmatrix(16,16)=1.0d0

u=matmul(MTEMP,Tmatrix)
Mtemp=u




!---------------------------------------------------------------
Gmatrix=0.0d0
Gmatrix(1,1)=1.d0
Gmatrix(2,4)=temp
Gmatrix(2,6)=-temp
Gmatrix(3,2)=temp
Gmatrix(3,3)=-temp
Gmatrix(4,5)=temp
Gmatrix(4,7)=-temp
Gmatrix(5,4)=temp
Gmatrix(5,6)=temp
Gmatrix(6,2)=temp
Gmatrix(6,3)=temp
Gmatrix(7,5)=temp
Gmatrix(7,7)=temp
DO i=8,nenergy
Gmatrix(i,i)=1.0d0
ENDDO

u=matmul(Mtemp,Gmatrix)
Mtemp=u
!--------------------------------------------------------------
Cmatrix=0.0d0
!1,3 states are both 3sigma_u^+ so they mix
IF(abs(a(1)).lt.1d-10.or.abs(exp(b(1)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(1)*r**(npower))*exp(b(1)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP1=0.0d0
CTEMP1(1,1)=phic
CTEMP1(1,2)=phis
CTEMP1(2,1)=-phis
CTEMP1(2,2)=phic
CTEMP1(3,3)=1.0d0
CTEMP1(4,4)=1.0d0

IF(abs(a(2)).lt.1d-10.or.abs(exp(b(2)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(2)*r**(npower))*exp(b(2)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(1,1)=phic
CTEMP2(1,3)=phis
CTEMP2(3,1)=-phis
CTEMP2(3,3)=phic
CTEMP2(2,2)=1.0d0
CTEMP2(4,4)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(3)).lt.1d-10.or.abs(exp(b(3)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(3)*r**(npower))*exp(b(3)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(1,1)=phic
CTEMP2(1,4)=phis
CTEMP2(4,1)=-phis
CTEMP2(4,4)=phic
CTEMP2(2,2)=1.0d0
CTEMP2(3,3)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(4)).lt.1d-10.or.abs(exp(b(4)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(4)*r**(npower))*exp(b(4)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(2,2)=phic
CTEMP2(2,3)=phis
CTEMP2(3,2)=-phis
CTEMP2(3,3)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(4,4)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(5)).lt.1d-10.or.abs(exp(b(5)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(5)*r**(npower))*exp(b(5)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(2,2)=phic
CTEMP2(2,4)=phis
CTEMP2(4,2)=-phis
CTEMP2(4,4)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(3,3)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)
CTEMP1=Ctemp

IF(abs(a(6)).lt.1d-10.or.abs(exp(b(6)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(6)*r**(npower))*exp(b(6)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CTEMP2=0.0d0
CTEMP2(3,3)=phic
CTEMP2(3,4)=phis
CTEMP2(4,3)=-phis
CTEMP2(4,4)=phic
CTEMP2(1,1)=1.0d0
CTEMP2(2,2)=1.0d0
CTEMP=matmul(CTEMP1,CTEMP2)

CMATRIX=0.0d0
CMATRIX(1,1)=CTEMP(1,1)
CMATRIX(1,3)=CTEMP(1,2)
CMATRIX(1,8)=CTEMP(1,3)
CMATRIX(1,12)=CTEMP(1,4)
CMATRIX(3,1)=CTEMP(2,1)
CMATRIX(3,3)=CTEMP(2,2)
CMATRIX(3,8)=CTEMP(2,3)
CMATRIX(3,12)=CTEMP(2,4)
CMATRIX(8,1)=CTEMP(3,1)
CMATRIX(8,3)=CTEMP(3,2)
CMATRIX(8,8)=CTEMP(3,3)
CMATRIX(8,12)=CTEMP(3,4)
CMATRIX(12,1)=CTEMP(4,1)
CMATRIX(12,3)=CTEMP(4,2)
CMATRIX(12,8)=CTEMP(4,3)
CMATRIX(12,12)=CTEMP(4,4)

!Pi_g coupling
IF(abs(a(7)).lt.1d-10.or.abs(exp(b(7)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(7)*r**(npower))*exp(b(7)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CMATRIX(4,4)=phic
CMATRIX(4,10)=phis
CMATRIX(5,5)=phic
CMATRIX(5,11)=-phis
CMATRIX(10,4)=-phis
CMATRIX(10,10)=phic
CMATRIX(11,5)=phis
CMATRIX(11,11)=phic

!Pi_u coupling
IF(abs(a(8)).lt.1d-10.or.abs(exp(b(8)*r)).gt.1d20)THEN
phic=1.0d0
phis=0.0d0
ELSE
phi=1.0d0/(a(8)*r**(npower))*exp(b(8)*r)
phi2=phi*phi
phic=sqrt(phi2/(1.0d0+phi2))
phis=sqrt(1.0d0/(1.0d0+phi2))
ENDIF
CMATRIX(6,6)=phic
CMATRIX(6,13)=phis
CMATRIX(13,6)=-phis
CMATRIX(13,13)=phic
CMATRIX(7,7)=phic
CMATRIX(7,14)=-phis
CMATRIX(14,7)=phis
CMATRIX(14,14)=phic

CMATRIX(2,2)=1.0d0
CMATRIX(9,9)=1.0d0
CMATRIX(15,15)=1.0d0
CMATRIX(16,16)=1.0d0
u=matmul(Mtemp,Cmatrix)

!----------------------------------
! test for Umatrix
!---------------------------------
IF(DEBUG)THEN
DO i=1,nenergy
   DO j=1,nenergy
      temp=0.0d0
      DO k=1,nenergy
      temp=temp+U(k,i)*U(k,j)
      ENDDO
!      print*,'i,j',i,j,temp
      IF(i.eq.j)THEN
           IF(abs(temp-1.0d0).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 3_4'
           ENDIF
      ELSE
           IF(abs(temp).gt.1.0d-6)THEN
           print*,'wrong 2dbasis=',i,j,temp
           stop 'error 3_4'
           ENDIF
      ENDIF
   ENDDO
ENDDO
print*,'Umatrix check is done'
stop
ENDIF
ELSE
print*,'unable to do more than 3 atoms now'
stop 'Umatrix in lib'
ENDIF
return
END subroutine
