SUBROUTINE POT16(r,Istate,E)
! an equivalent of subroutine POT7 ..and we have 16 states instead of 7 states
IMPLICIT NONE
real*8 r,E,tab(3),EP
INTEGER Istate,ithcall,iop(2),i
INTEGER,parameter:: nr1=51
Real*8 V(nr1,16),RV(nr1),temp
real*8 aa(nr1),bb(nr1),cc(nr1)
REAL*8 SECOND_1(nr1),SECOND_2(nr1),SECOND_3(nr1),SECOND_4(nr1),SECOND_6(nr1),shiftss,shiftsp,shiftpp
REAL*8 SECOND_8(nr1),SECOND_9(nr1),SECOND_10(nr1),SECOND_12(nr1),SECOND_13(nr1),SECOND_15(nr1)
data ithcall /0/
!NEW FCI data with a new p basis
data shiftss /-0.39207225d0/
data shiftsp /-0.32438499d0/
data shiftpp /-0.25669772d0/
data EP/0.06768727d0/
!--------------------------------------------------------------
!define the energy for XY diatom
!E1=3sigma_u^+
!E2=3sigma_g^+(sz case)
!E3=3sigma_u^+(sz case)
!E4=E5=3pi_g^+
!E6=E7=3pi_u^+
!---------------------------------------------------------------


!there are 33 lines of mcI files
IF(ithcall.eq.0)THEN
ithcall=1
!state 1,3,8,12
OPEN(unit=910,file='Potentials/WEB_LI3/3siu5.fci.table')
DO i=1,nr1
READ(910,*)RV(i),V(i,1),V(i,3),V(i,8),V(i,12)
ENDDO
CLOSE(910)

! state 2
OPEN(unit=910,file='Potentials/WEB_LI3/3sig1.fci.table')
DO i=1,nr1
READ(910,*)RV(i),V(i,2)
ENDDO
CLOSE(910)
! state (4,5),(10,11)
OPEN(unit=910,file='Potentials/WEB_LI3/3pig6.fci.table')
DO i=1,nr1
READ(910,*)RV(i),V(i,4),V(i,10)
ENDDO
CLOSE(910)

!state (6,7) and (13,14)
OPEN(unit=910,file='Potentials/WEB_LI3/3piu2.fci.table')
DO i=1,nr1
READ(910,*)RV(i),V(i,6),V(i,13)
ENDDO
CLOSE(910)

!state (15,16)
OPEN(unit=910,file='Potentials/WEB_LI3/3del8.fci.table')
DO i=1,nr1
READ(910,*)RV(i),V(i,15)
ENDDO
CLOSE(910)

!state 9
OPEN(unit=910,file='Potentials/WEB_LI3/3del4.fci.table')
DO i=1,nr1
READ(910,*)RV(i),V(i,9)
ENDDO
CLOSE(910)

CALL spl1d1 (nr1, RV(1), V(1,1), SECOND_1(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,2), SECOND_2(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,3), SECOND_3(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,4), SECOND_4(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,6), SECOND_6(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,8), SECOND_8(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,9), SECOND_9(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,10), SECOND_10(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,12), SECOND_12(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,13), SECOND_13(1), iop, 1, aa (1), bb (1), cc (1) )
CALL spl1d1 (nr1, RV(1), V(1,15), SECOND_15(1), iop, 1, aa (1), bb (1), cc (1) )
print*,'Initialization DONE'
ENDIF

IF(Istate.eq.1)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,1), SECOND_1(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,1)
ELSE
E=V(nr1,1)
ENDIF
E=E-shiftss

return
ELSEIF(Istate.eq.2)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,2), SECOND_2(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,2)
ELSE
E=V(nr1,2)
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.3)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,3), SECOND_3(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,3)
ELSE
E=V(nr1,3)
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.4)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,4), SECOND_4(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,4)
ELSE
!E=V(nr1,4)
E=EP+(V(nr1,4)-shiftss-EP)*exp(-2d0*(r-RV(nr1)))+shiftss
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.6)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,6), SECOND_6(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,6)
ELSE
!E=V(nr1,6)
E=EP+(V(nr1,6)-shiftss-EP)*exp(-2d0*(r-RV(nr1)))+shiftss
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.8)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,8), SECOND_8(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,8)
ELSE
E=V(nr1,8)
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.9)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,9), SECOND_9(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,9)
ELSE
E=V(nr1,9)
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.10)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,10), SECOND_10(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,10)
ELSE
E=V(nr1,10)
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.12)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,12), SECOND_12(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,12)
ELSE
E=V(nr1,12)
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.13)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,13), SECOND_13(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,13)
ELSE
E=V(nr1,13)
ENDIF
E=E-shiftss
return

ELSEIF(Istate.eq.15)THEN
IF(r.ge.RV(1).and.r.le.RV(nr1))THEN
CALL spl1d2 (nr1, RV(1), V(1,15), SECOND_15(1), 1, r, tab)
E=tab(1)
ELSEIF(r.lt.RV(1))THEN
E=V(1,15)
ELSE
E=V(nr1,15)
ENDIF
E=E-shiftss
return

ELSE 
stop 'WRONG istate'
ENDIF

END
