Subroutine Vhamil_2AP_CIV(E_hamil,r)
USE FileUnits_Module
IMPLICIT NONE
INTEGER,parameter:: nbasis=37,nenergy=16,nbasis2=22
INTEGER icase,iupper,jupper,info
INTEGER,ALLOCATABLE::iwork(:),ifail(:)
real*8 U(nenergy,nenergy),Hamil(nbasis,nbasis)
real*8 Hamil2(nbasis,nbasis)
real*8 Hamiltemp(nbasis2,nbasis2),temp,E_hamil(nbasis)
real*8 r(3),angle,E_p,anglec,angles
real*8 v(nenergy),xab,xbc,xac,a(8),b(8)
REAL*8,ALLOCATABLE::sd(:),subd(:),vtemp(:,:)
real*8 vli2,vli2_3sig,vli2_3siu2,vli2_3piu,vli2_3pig
INTEGER i,j,k,indexi,indexj,ngri,iflag,icount
data a(1)/ 0.14535747845035D-01/
data b(1)/0.24857019396506D+00/
data a(2)/0.14946553159223D-02/
data b(2)/0.74292801158851D-01 /
data a(3)/0.73180454277552D-03/
data b(3)/0.75185999079706D-01 /
data a(4)/0.0d0 /
data b(4)/1.0d0 /
data a(5)/0.0d0 /
data b(5)/1.0d0 /
data a(6)/0.0d0 /
data b(6)/1.0d0 /
data a(7)/0.0d0 /
data b(7)/1.0d0 /
data a(8)/0.25323569592969D+00 /
data b(8)/0.74167570269811D+00 /


!---------------------------------------------------------------
!this routine is written for the diabatic potential
!It is written by Xuan LI
!Contact : XUANLI@OU.EDU
!--------------------------------------------------------
!iflag=1  !for my potential
!iflag=2  !for Daniel Brue's program
!IFLAG=3 ! use mc input files: from molpro
!IFLAG=4 ! use Boss' spk data for 16 states
IFLAG=5 !use BOSS's SPK spdf data to get A' states
xab=r(1)
xbc=r(2)
xac=r(3)
Hamil=0.0d0
Hamiltemp=0.0d0






!--------------------------------------------------------------
!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^+
!-------------------------------------------------------------
!warning! there are 10 basis, however in each diatom we use 7 out of them,
! so be careful with the basis, the Hamiltonian is constructed as 10x10
! then we reduce it to 7x7 without the A'' states
!-----------------------------------------------
!basis notation
! |1> =sss
! |2> =ssx
! |3> =ssy
! |4> =ssz
! |5> =sxs
! |6> =sys
! |7> =szs
! |8> =xss
! |9> =yss
! |10>=zss
!----------------------------------------------


!---------------------------------------------
!cosntruc h(m,m')_A+h(m,m')_B+h(m,m')_C
 IF(IFLAG.eq.5)E_p=(0.19603613d0-0.12834886d0)
 Hamil(1,1)= 0.0d0
 Hamil(2,2)= -E_p
 Hamil(3,3)= -E_p
 Hamil(4,4)= -E_p
 Hamil(5,5)= -E_p
 Hamil(6,6)= -E_p
 Hamil(7,7)= -E_p
 Hamil(8,8)= -E_p
 Hamil(9,9)= -E_p
 Hamil(10,10)= -E_p
 DO i=11,nbasis
 HAMIL(i,i)=-2*E_p
 ENDDO
!-----------------------------------------------
!now add the h(m,m')_AB
!-----------------------------------------------
IF(IFLAG.eq.4.or.IFLAG.eq.5)THEN
!call pot16(xab,1,v(1))
!USE Flavio's codes
v(1)=vli2(xab)
call pot16(xab,2,v(2))
call pot16(xab,3,v(3))
call pot16(xab,4,v(4))
v(5)=v(4)
call pot16(xab,6,v(6))
v(7)=v(6)
call pot16(xab,8,v(8))
call pot16(xab,9,v(9))
call pot16(xab,10,v(10))
v(11)=v(10)
call pot16(xab,12,v(12))
call pot16(xab,13,v(13))
v(14)=v(13)
call pot16(xab,15,v(15))
v(16)=v(15)
ELSE
stop 'wrong FLAG'
ENDIF



!rotational angles
! c for cos
! s for sin
    anglec=(xac*xac-xab*xab-xbc*xbc)/(2.0d0*xab*xbc)
    IF(abs(anglec-1.0d0).le.1D-10)THEN
    angles=0.0d0
    anglec=1.0d0
    ELSE
    angles=sqrt(abs(1.0000000d0-anglec*anglec))
    ENDIF

     call U_matrix_2(nbasis,nenergy,U,1,anglec,angles,xab,a,b)

hamil2=0.0d0
indexi=0
DO i=1,nbasis
    IF(i.eq.1.or.(i.ge.5.and.i.le.10).or.(i.ge.29))THEN
      indexi=indexi+1
      indexj=0
      DO j=1,nbasis
         IF(j.eq.1.or.(j.ge.5.and.j.le.10).or.(j.ge.29))THEN
         indexj=indexj+1
         temp=0.0d0
            DO k=1,nenergy
              temp=temp+v(k)*U(indexi,k)*U(indexj,k)
            ENDDO
         Hamil2(i,j)=temp
          ENDIF
      ENDDO
    ENDIF
ENDDO
   DO i=2,4
    Hamil2(i,i)=HAMIL2(1,1)
   ENDDO
   DO i=11,13
    Hamil2(i,i)=HAMIL2(5,5)
   ENDDO
   DO i=14,16
    Hamil2(i,i)=HAMIL2(6,6)
   ENDDO
   DO i=17,19
    Hamil2(i,i)=HAMIL2(7,7)
   ENDDO
   DO i=20,22
    Hamil2(i,i)=HAMIL2(8,8)
   ENDDO
   DO i=23,25
    Hamil2(i,i)=HAMIL2(9,9)
   ENDDO
   DO i=26,28
    Hamil2(i,i)=HAMIL2(10,10)
   ENDDO

   DO i=1,nbasis
     DO j=1,nbasis
       HAMIL(i,j)=HAMIL(i,j)+HAMIL2(i,j)
     ENDDO
   ENDDO
!-----------------------------------------------
!now add the h(m,m')_BC
!-----------------------------------------------

IF(iflag.eq.4.or.iflag.eq.5)THEN
!call pot16(xbc,1,v(1))
!USE Flavio's codes
v(1)=vli2(xbc)
call pot16(xbc,2,v(2))
call pot16(xbc,3,v(3))
call pot16(xbc,4,v(4))
v(5)=v(4)
call pot16(xbc,6,v(6))
v(7)=v(6)
call pot16(xbc,8,v(8))
call pot16(xbc,9,v(9))
call pot16(xbc,10,v(10))
v(11)=v(10)
call pot16(xbc,12,v(12))
call pot16(xbc,13,v(13))
v(14)=v(13)
call pot16(xbc,15,v(15))
v(16)=v(15)
ELSE
stop 'wrong FLAG'
ENDIF


    call U_matrix_2(nbasis,nenergy,U,2,0.0d0,0.0d0,xbc,a,b)
hamil2=0.0d0
indexi=0
DO i=1,nbasis
    IF((i.ge.1.and.i.le.7).or.(i.ge.11.and.i.le.19))THEN
      indexi=indexi+1
      indexj=0
      DO j=1,nbasis
         IF((j.ge.1.and.j.le.7).or.(j.ge.11.and.j.le.19))THEN
         indexj=indexj+1
         temp=0.0d0
            DO k=1,nenergy
              temp=temp+v(k)*U(indexi,k)*U(indexj,k)
            ENDDO
         Hamil2(i,j)=temp
          ENDIF
      ENDDO
    ENDIF
ENDDO

   DO i=8,10
    Hamil2(i,i)=HAMIL2(1,1)
   ENDDO
   DO i=20,26,3
    Hamil2(i,i)=HAMIL2(2,2)
   ENDDO
   DO i=21,27,3
    Hamil2(i,i)=HAMIL2(3,3)
   ENDDO
   DO i=22,28,3
    Hamil2(i,i)=HAMIL2(4,4)
   ENDDO
   DO i=29,35,3
    Hamil2(i,i)=HAMIL2(5,5)
   ENDDO
   DO i=30,36,3
    Hamil2(i,i)=HAMIL2(6,6)
   ENDDO
   DO i=31,37,3
    Hamil2(i,i)=HAMIL2(7,7)
   ENDDO

   DO i=1,nbasis
     DO j=1,nbasis
       HAMIL(i,j)=HAMIL(i,j)+HAMIL2(i,j)
     ENDDO
   ENDDO
!-----------------------------------------------
!now add the h(m,m')_AC
!-----------------------------------------------
IF(iflag.eq.4.or.iflag.eq.5)THEN
!call pot16(xac,1,v(1))
!USE Flavio's codes
v(1)=vli2(xac)
call pot16(xac,2,v(2))
call pot16(xac,3,v(3))
call pot16(xac,4,v(4))
v(5)=v(4)
call pot16(xac,6,v(6))
v(7)=v(6)
call pot16(xac,8,v(8))
call pot16(xac,9,v(9))
call pot16(xac,10,v(10))
v(11)=v(10)
call pot16(xac,12,v(12))
call pot16(xac,13,v(13))
v(14)=v(13)
call pot16(xac,15,v(15))
v(16)=v(15)
ELSE
stop 'wrong FLAG'
ENDIF

    anglec=(-xab*xab+xbc*xbc+xac*xac)/(2.0d0*xac*xbc)
    IF(abs(anglec-1.0d0).le.1D-10)THEN
    angles=0.0d0
    ELSE
    angles=sqrt(abs(1.0000000d0-anglec*anglec))
    ENDIF
    call U_matrix_2(nbasis,nenergy,U,3,anglec,angles,xac,a,b)

hamil2=0.0d0
indexi=0
DO i=1,nbasis
    IF((i.ge.1.and.i.le.4).or.(i.ge.8.and.i.le.10).or.(i.ge.20.and.i.le.28))THEN
      indexi=indexi+1
      indexj=0
      DO j=1,nbasis
         IF((j.ge.1.and.j.le.4).or.(j.ge.8.and.j.le.10).or.(j.ge.20.and.j.le.28))THEN
         indexj=indexj+1
         temp=0.0d0
            DO k=1,nenergy
              temp=temp+v(k)*U(indexi,k)*U(indexj,k)
            ENDDO
         Hamil2(i,j)=temp
          ENDIF
      ENDDO
    ENDIF
ENDDO

   DO i=5,7
    Hamil2(i,i)=HAMIL2(1,1)
   ENDDO
   DO i=11,17,3
    Hamil2(i,i)=HAMIL2(2,2)
   ENDDO
   DO i=12,18,3
    Hamil2(i,i)=HAMIL2(3,3)
   ENDDO
   DO i=13,19,3
    Hamil2(i,i)=HAMIL2(4,4)
   ENDDO
   DO i=29,31
    Hamil2(i,i)=HAMIL2(8,8)
   ENDDO
   DO i=32,34
    Hamil2(i,i)=HAMIL2(9,9)
   ENDDO
   DO i=35,37
    Hamil2(i,i)=HAMIL2(10,10)
   ENDDO

     DO i=1,nbasis
     DO j=1,nbasis
       HAMIL(i,j)=HAMIL(i,j)+HAMIL2(i,j)
     ENDDO
   ENDDO

299 FORMAT(E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x)
145 CONTINUE
!-------------------------------------------------------------
!now diagonalize the Hamil to get the potential
!-------------------------------------------------------------
!first, pick out the 22x22 submatrix


indexi=0
DO i=1,nbasis

   indexj=0
 IF(i.ne.3.and.i.ne.6.and.i.ne.9.and.i.ne.12.and.i.ne.18.and.i.ne.21.and.i.ne.23.and.i.ne.16)THEN
 IF(i.ne.25.and.i.ne.27.and.i.ne.30.and.i.ne.32.and.i.ne.34.and.i.ne.36.and.i.ne.14)THEN
   indexi=indexi+1
   DO j=1,nbasis
   IF(j.ne.3.and.j.ne.6.and.j.ne.9.and.j.ne.12.and.j.ne.18.and.j.ne.21.and.j.ne.23.and.j.ne.16)THEN
    IF(j.ne.25.and.j.ne.27.and.j.ne.30.and.j.ne.32.and.j.ne.34.and.j.ne.36.and.j.ne.14)THEN
     indexj=indexj+1
     HAMILtemp(indexi,indexj)=hamil(i,j)
    ENDIF
   ENDIF
   ENDDO
  ENDIF
  ENDIF
ENDDO
     

298 FORMAT(E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x,E8.2,1x)
297 FORMAT(7(E11.5,1x))

ngri=nbasis2
k=0


ALLOCATE(sd(ngri*(ngri+1)/2))
DO  j=1,ngri
   DO  i=j,ngri
      k=k+1
      sd(k)=Hamiltemp(i,j)
   ENDDO
ENDDO


Hamiltemp=0.0d0


iupper=min(4000,ngri)        !iupper is the upper limit of the output eigenvalue numbers
jupper=iupper
ALLOCATE(vtemp(ngri,ngri))
ALLOCATE(subd(8*ngri),iwork(5*ngri),ifail(ngri))

CALL dspevx('V','I','L',ngri,sd,0.d0,0.d0,1,iupper,0.0d-12,jupper,E_hamil,vtemp,ngri,subd,iwork,ifail,info)

DEALLOCATE(sd)
DEALLOCATE(subd,iwork,ifail,vtemp)
return

END Subroutine
