module diabmodel
  use dim_parameter,only:qn,ndiab,pst
  use accuracy_constants, only:dp,idp 
  implicit none 
  logical :: debug=.false.
  contains

  subroutine diab(ex,ey,n,x1,x2,p)
      use ctrans_mod, only:ctrans
      integer,intent(in),optional :: n ! number of parameter & nmbr of points \
      integer id
      integer key,i,j
      double precision, intent(in)::x1(qn),x2(qn)
      double precision, contiguous,intent(in):: p(:)! array containing parameters
      double precision, intent(out)::ex(ndiab,ndiab),ey(ndiab,ndiab)
      key =87
     call diab_x(ex,x1,x2,key,p)
           !ey=0.0d0
      call diab_y(ey,x1,x2,key,p)
  end subroutine

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  subroutine diab_x(e,q,t,key,p)

    
    real(dp),intent(in)::q(qn),t(qn)
    real(dp),intent(out)::e(:,:)
    integer(idp),intent(in)::key
    real(dp),intent(in),contiguous::p(:)
    integer(idp) id,i,j 
    real(dp) tmp_v,xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
    xs=q(2)
    ys=q(3)
    xb=q(4)
    yb=q(5)
    a=q(1)
    b=q(6)

    ss=xs**2+ys**2 ! totaly symmetric term
    sb=xb**2+yb**2

    v3_vec( 1) = xs*(xs**2-3*ys**2)
    v3_vec( 2) = xb*(xb**2-3*yb**2)
    v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys
    v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb
    v3_vec( 5) = ys*(3*xs**2-ys**2)
    v3_vec( 6) = yb*(3*xb**2-yb**2)
    v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
    v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb 



    e=0.0d0


   
    id=key   !1
    ! V-term 
    ! order 1
    e(1,1)=e(1,1)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb 
    id=id+1   !2
    e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb 
    e(3,3)=e(3,3)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb 
    id=id+1   !3
    e(4,4)=e(4,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
    ! order 2
    id=id+1   !4
    e(1,1)=e(1,1)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)&
                 +p(pst(1,id)+2)*(xs*xb-ys*yb) 
    id=id+1   !5
    e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) +&
                  p(pst(1,id)+2)*(xs*xb-ys*yb) 
    e(3,3)=e(3,3)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) +&
                  p(pst(1,id)+2)*(xs*xb-ys*yb) 
    id=id+1   !6
    e(4,4)=e(4,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) + &
                   p(pst(1,id)+2)*(xs*xb-ys*yb) 
     ! order 3
    id=id+1   !7
    e(1,1)=e(1,1)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb + b**2* &
                  (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)
    id=id+1   !8
    e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb+ b**2* &
                  (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)
    e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb+ b**2* &
                  (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)
    id=id+1   !9
    e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb + b**2* &
                  (p(pst(1,id)+2)*xs +p(pst(1,id)+3)*xb)

    ! JAHN TELLER COUPLING W AND Z
    ! order 0
    id=id+1   !10
    e(2,2)=e(2,2)+p(pst(1,id))
    e(3,3)=e(3,3)-p(pst(1,id))
    ! order 1
    id=id+1   !11
    e(2,2)=e(2,2)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb 
    e(3,3)=e(3,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb 
    e(2,3)=e(2,3)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb
    ! order 2
    id=id+1   !12
    e(2,2)=e(2,2)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                  +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb + &
                  b**2*(p(pst(1,id)+5) +b**4*(p(pst(1,id)+6)) + b**6*(p(pst(1,id)+7)) + & 
                  b**8*(p(pst(1,id)+8)))
    e(3,3)=e(3,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                  +p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb) -&
                  b**8*(p(pst(1,id)+5))
    e(2,3)=e(2,3)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+ &
                  p(pst(1,id)+2)*(xs*yb+xb*ys) 
    ! order 3

    id=id+1   !13
    do i=1,4
      j=i-1
      e(2,2)=e(2,2)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
      e(3,3)=e(3,3)-(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
      e(2,3)=e(2,3)+(-p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i+4)
    enddo 

    e(2,2)=e(2,2)+p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb 
    e(3,3)=e(3,3)-(p(pst(1,id)+8)*xs*ss+p(pst(1,id)+9)*xb*sb)
    e(2,3)=e(2,3)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb 
    ! PSEUDO JAHN TELLER 

    ! A2 ground state coupled with E 
    ! ###################################################
    ! ###################################################

    ! order 0
    id=id+1   !14
    e(1,2)=e(1,2)+b*p(pst(1,id))
    
    ! order 1 
    id=id+1   !15 
    e(1,2)=e(1,2)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
    e(1,3)=e(1,3)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb)
    ! order 2
    id=id+1   !16
    e(1,2)=e(1,2)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)&
                 +p(pst(1,id)+2)*(xs*xb-ys*yb))
    e(1,3)=e(1,3)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)&
                 +p(pst(1,id)+2)*(xs*yb+xb*ys))

    !! THE COUPLING OF A2 WITH A1
    !####################################################
    !####################################################
    ! order 1
    id=id+1   !17
    e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
    id=id+1   !18
    e(1,4)=e(1,4)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)&
                 +p(pst(1,id)+2)*(xs*xb-ys*yb)) 

      
    !!! THE COUPLING OF A1 WITH E 
    !!####################################################
    !####################################################
    ! order 0
    id=id+1   !19
    e(2,4)=e(2,4)+p(pst(1,id))

    ! order 1
    id=id+1   !20
    e(2,4)=e(2,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
    e(3,4)=e(3,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb

    ! order 2
    id=id+1   !21
    e(2,4)=e(2,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                 +p(pst(1,id)+2)*(xs*xb-ys*yb)
    e(3,4)=e(3,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) &
                 -p(pst(1,id)+2)*(xs*yb+xb*ys)

    !! End of the model 

    e(2,1)=e(1,2)
    e(3,1)=e(1,3) 
    e(3,2)=e(2,3)
    e(4,1)=e(1,4)
    e(4,2)=e(2,4)
    e(4,3)=e(3,4)
  end subroutine diab_x

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! THE Y COMPONENT OF DIPOLE 
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine diab_y(e,q,t,key,p)
    !integer(idp), intent(in)::npar 
    real(dp),intent(in)::q(qn),t(qn)
    real(dp),intent(out)::e(:,:)
    integer(idp),intent(in):: key
    real(dp),intent(in),contiguous::p(:)
    integer(idp) id,i,j 
    real(dp) tmp_v,ys,xb,a,b,xs,yb,ss,sb,v3_vec(8) 
    xs=q(2)
    ys=q(3)
    xb=q(4)
    yb=q(5)
    a=q(1)
    b=q(6)

    ss=xs**2+ys**2 ! totaly symmetric term
    sb=xb**2+yb**2

    v3_vec( 1) = xs*(xs**2-3*ys**2)
    v3_vec( 2) = xb*(xb**2-3*yb**2)
    v3_vec( 3) = xb*(xs**2-ys**2) - 2*yb*xs*ys
    v3_vec( 4) = xs*(xb**2-yb**2) - 2*ys*xb*yb
    v3_vec( 5) = ys*(3*xs**2-ys**2)
    v3_vec( 6) = yb*(3*xb**2-yb**2)
    v3_vec( 7) = yb*(xs**2-ys**2)+2*xb*xs*ys
    v3_vec( 8) = ys*(xb**2-yb**2)+2*xs*xb*yb 





    e=0.0d0
    ! V-term 
    id=key   !1
    e(1,1)=e(1,1)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb

    id=id+1   !2
    e(2,2)=e(2,2)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb 
    e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb 
    id=id+1   !3
    e(4,4)=e(4,4)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb 
    ! order 2
    id=id+1   !4
    e(1,1)=e(1,1)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) &
                 -p(pst(1,id)+2)*(xs*yb+xb*ys)
    id=id+1   !5
    e(2,2)=e(2,2)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) &
                 -p(pst(1,id)+2)*(xs*yb+xb*ys)

   
    e(3,3)=e(3,3)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) &
                 -p(pst(1,id)+2)*(xs*yb+xb*ys)
    id=id+1   !6
    e(4,4)=e(4,4)-p(pst(1,id))*(2*xs*ys)-p(pst(1,id)+1)*(2*xb*yb) &
                 -p(pst(1,id)+2)*(xs*yb+xb*ys)
    ! order 3
    id=id+1   !7
    e(1,1)=e(1,1)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* &
                  (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)
    id=id+1   !8
    e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb+b**2* &
                  (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)
    e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* &
                  (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)
    id=id+1   !9
    e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb +b**2* &
                  (p(pst(1,id)+2)*ys +p(pst(1,id)+3)*yb)

    ! V- term + totally symmetric coord a

    ! JAHN TELLER COUPLING TERM 
     ! order 0
    id=id+1   !10
    e(2,3)=e(2,3)+p(pst(1,id))
    ! order 1

    id=id+1   !11
    e(2,2)=e(2,2)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb
    e(3,3)=e(3,3)+p(pst(1,id))*ys+p(pst(1,id)+1)*yb
    e(2,3)=e(2,3)-p(pst(1,id))*xs-p(pst(1,id)+1)*xb
    !id=id+1   !12
    ! order 2 
    id=id+1   !12
    e(2,2)=e(2,2)+p(pst(1,id))*2*xs*ys+p(pst(1,id)+1)*2*xb*yb+p(pst(1,id)+2)*(xs*yb+xb*ys)
    e(3,3)=e(3,3)-p(pst(1,id))*2*xs*ys-p(pst(1,id)+1)*2*xb*yb-p(pst(1,id)+2)*(xs*yb+xb*ys)
    e(2,3)=e(2,3)-(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)) &
                  -p(pst(1,id)+2)*(xs*xb-ys*yb)+p(pst(1,id)+3)*ss+p(pst(1,id)+4)*sb +& 
                  b**8*(p(pst(1,id)+5))
    ! order 3 
    id=id+1   !13 
    do i=1,4
        j=i-1
        e(2,2)=e(2,2)+(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4)
        e(3,3)=e(3,3)-(p(pst(1,id)+j)-p(pst(1,id)+j+4))*v3_vec(i+4)
        e(2,3)=e(2,3)+(p(pst(1,id)+j)+p(pst(1,id)+j+4))*v3_vec(i)
    enddo 
    e(2,2)=e(2,2)-p(pst(1,id)+8)*ys*ss-p(pst(1,id)+9)*yb*sb
    e(3,3)=e(3,3)+p(pst(1,id)+8)*ys*ss+p(pst(1,id)+9)*yb*sb
    e(2,3)=e(2,3)-p(pst(1,id)+8)*xs*ss-p(pst(1,id)+1)*xb*sb 
    ! PSEUDO JAHN TELLER 
    ! ORDER 0 
    ! THE COUPLING OF A2 GROUND STATE WITH E
    ! ###################################################
    ! ################################################### 
    ! order 0
    id=id+1   !14 
    e(1,3)=e(1,3)-b*(p(pst(1,id)))
    ! order 1
    id=id+1   !15
    e(1,2)=e(1,2)-b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb)
    e(1,3)=e(1,3)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)

    ! order 2
    id=id+1   !16
    e(1,2)=e(1,2)+b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)&
                 +p(pst(1,id)+2)*(xs*yb+xb*ys))
    e(1,3)=e(1,3)+b*(p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2)&
                 +p(pst(1,id)+2)*(xs*xb-ys*yb))
    ! THE COUPLING OF A2 WITH A1
    !####################################################
    !####################################################
    ! order 1
    id=id+1   !17
    e(1,4)=e(1,4)+b*(p(pst(1,id))*ys+p(pst(1,id)+1)*yb)
    ! order 2
    id=id+1   !18
    e(1,4)=e(1,4)-b*(p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb)&
                 +p(pst(1,id)+2)*(xs*yb+xb*ys))

    
    ! THE COUPLING OF A1 WITH E
    !####################################################
    !####################################################
    ! order 0
    id=id+1   !19
    e(3,4)=e(3,4)-p(pst(1,id))  
    ! order 1 
    id=id+1   !20
    e(2,4)=e(2,4)-p(pst(1,id))*ys-p(pst(1,id)+1)*yb
    e(3,4)=e(3,4)+p(pst(1,id))*xs+p(pst(1,id)+1)*xb
    ! order 2 
    id=id+1   !21
    e(2,4)=e(2,4)+p(pst(1,id))*(2*xs*ys)+p(pst(1,id)+1)*(2*xb*yb) &
                 +p(pst(1,id)+2)*(xs*yb+xb*ys)
    e(3,4)=e(3,4)+p(pst(1,id))*(xs**2-ys**2)+p(pst(1,id)+1)*(xb**2-yb**2) &
                 +p(pst(1,id)+2)*(xs*xb-ys*yb)  
    ! end of the model 
    e(2,1)=e(1,2)
    e(3,1)=e(1,3) 
    e(3,2)=e(2,3)
    e(4,1)=e(1,4)
    e(4,2)=e(2,4)
    e(4,3)=e(3,4)
end subroutine diab_y
 subroutine copy_2_lower_triangle(mat)
    real(dp), intent(inout) :: mat(:, :)
    integer :: m, n
!     write lower triangle of matrix symmetrical
    do n = 2, size(mat, 1)
      do m = 1, n - 1
        mat(n, m) = mat(m, n)
      end do
    end do
  end subroutine copy_2_lower_triangle

end module diabmodel