module diabmodel
  use iso_fortran_env, only: idp => int32, dp => real64
  use dip_param
  implicit none 
  include "nnparams.incl"
  integer(idp),parameter:: ndiab=4
  logical :: debug=.false.
  contains

  subroutine diab_x(e,q,nn_out)
    real(dp),intent(in)::q(maxnin)
    real(dp),intent(inout):: nn_out(maxnout)
    real(dp),intent(out)::e(ndiab,ndiab)
    integer(idp) id,i,j ,ii
    real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
    real(dp),dimension(12):: shift,scal
    xs=q(5)
    ys=q(6)
    xb=q(7)
    yb=q(8)
    a=q(4)
    b=q(9)

    call init_dip_planar_data()

    ! modify the parametr
    !shift=1.0_dp
    !scal=1.0d-2
    scal = [0.69753,0.44797,51.14259,2.76924,1.45300,91.58246, & 
            14.07390,1.02550,1.68623,4.80804,7.69958,0.97871]
    shift = [0.23479,0.26845,35.05940,2.27175,-0.33017,117.48895, & 
            -1.68211,0.79418,-1.60443,-10.41309,8.47695,1.25334]
    ! V term of A2''
    ii=1
    do i =1,np 
      if (p(i) .ne. 0.0d0) then 
        p(i) =p(i)*(shift(ii) + scal(ii)*nn_out(ii) )
        ii=ii+1
      else
        p(i) = p(i)
      endif 
    enddo

    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=1   !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
    id=id+1   !8
    e(2,2)=e(2,2)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
    e(3,3)=e(3,3)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb
    id=id+1   !9
    e(4,4)=e(4,4)+p(pst(1,id))*xs*ss+p(pst(1,id)+1)*xb*sb 

    ! 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
    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)
    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) + p(pst(1,id)+3)*(xs**2+ys**2))
    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))
    ! order 3 
    id =id+1 ! 17

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

                 

    !! THE COUPLING OF A2 WITH A1
    !####################################################
    !####################################################
    ! order 1
    id=id+1   !18
    e(1,4)=e(1,4)+b*(p(pst(1,id))*xs+p(pst(1,id)+1)*xb)
    id=id+1   !19
    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   !20
    e(2,4)=e(2,4)+p(pst(1,id))

    ! order 1
    id=id+1   !21
    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   !22
    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) +p(pst(1,id)+3)*(xs**2+ys**2)
    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) 
    ! order 3
    id=id+1 !23
    do i=1,4
      e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
      e(3,4)=e(3,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
    enddo 

    !! 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,nn_out)
    real(dp),intent(in)::q(maxnin)
    real(dp),intent(inout):: nn_out(maxnout)
    real(dp),intent(out)::e(ndiab,ndiab)
    integer(idp) id,i,j, ii 
    real(dp) xs,xb,ys,yb,a,b,ss,sb,v3_vec(8)
     real(dp),dimension(12):: shift,scal
    xs=q(5)
    ys=q(6)
    xb=q(7)
    yb=q(8)
    a=q(4)
    b=q(9)

    call init_dip_planar_data()
    ! modify the parametr
    !shift=1.0_dp
    !scal=1.0d-3
    scal = [0.69753,0.44797,51.14259,2.76924,1.45300,91.58246, & 
            14.07390,1.02550,1.68623,4.80804,7.69958,0.97871]
    shift = [0.23479,0.26845,35.05940,2.27175,-0.33017,117.48895, & 
            -1.68211,0.79418,-1.60443,-10.41309,8.47695,1.25334]

    ! V term of A2''
    ii=1
    do i =1,np 
      if (p(i) .ne. 0) then 
        p(i) =p(i)*(shift(ii) + scal(ii)*nn_out(ii) )
        ii=ii+1
      else 
        p(i) = p(i)
      endif 
    enddo



    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=1   !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
    id=id+1   !8
    e(2,2)=e(2,2)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb
    e(3,3)=e(3,3)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb 
    id=id+1   !9
    e(4,4)=e(4,4)+p(pst(1,id))*ys*ss+p(pst(1,id)+1)*yb*sb

    ! 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
    ! 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) - p(pst(1,id)+3)*(xs**2+ys**2))

     ! order 3
    id = id+1 ! 17
    do i=1,4
      e(1,2)=e(1,2)+b*(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
      e(1,3)=e(1,3)-b*(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
    enddo 
          
                 
    ! 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) - p(pst(1,id)+3)*(xs**2+ys**2)

    id =id+1 ! 23
    ! order 3
    do i=1,4
      e(2,4)=e(2,4)+(p(pst(1,id)+(i-1))-p(pst(1,id)+(i+3)))*v3_vec(i+4)
      e(3,4)=e(3,4)-(p(pst(1,id)+(i-1))+p(pst(1,id)+(i+3)))*v3_vec(i)
    enddo        
    ! 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