module pesmod
    !use dim_parameter,only:qn,ndiab,pst
    use iso_fortran_env, only: idp => int32,  dp => real64
    use dip_param

    use nn_params
    use nncommons, only: len_in
    implicit none
    
    integer(idp),parameter:: ndiab=4
    private 
    public :: potential_nh3_low
    contains 


    subroutine potential_nh3_low(e,q,nn_out)
    implicit none 
    real(dp),intent(in)::q(maxnin)
    real(dp),intent(out)::e(:,:)
    real(dp),intent(inout):: nn_out(maxnout)
    !logical,intent(in):: sym
    !integer(idp),intent(in)::key
    integer(idp) id,i, ii
    real(dp) tmp,xs,xb,ys,yb,a,b,ss,sb
    logical, parameter:: switch = .false. ! parameter to switch when you want to know Ex and Ey
    logical, parameter:: sym = .false. ! to include r**2 in the Genetic fitting 
    real(dp),parameter:: tol = 1.0d-9
    real(dp),dimension(16):: scal, shift 
    !real(dp),dimension(17):: shi,scalee
    !include "shift-scale-70N.incl"
    !include "sh-scal-50.incl"
    xs=q(len_in+2)
    ys=q(len_in+3)
    xb=q(len_in+4)
    yb=q(len_in+5)
    a=q(len_in+1)
    b=q(len_in+6)

    ss= xs**2+ys**2
    sb= xb**2+yb**2
    call init_dip_planar_data()
    shift =1.0d0 
    scal = 1.0d-3
    !scal = 1.0d-3
    !scal(1:3) = 1.0d-3
    ii = 1
    do i =1 ,27  ! pes parameter key index 
      if( abs(p(i)) .gt. tol ) then 
        p(i) = p(i)*(shift(ii)+ scal(ii)*nn_out(ii))
        ii = ii + 1

      endif 
    enddo 
    !write(6,*)" PES non zero", ii-1
      

    e = 0.0_dp
    id = 1 ! 1 
    ! v -term 
    ! oth order only  
     ! A2"
    e(1,1) = e(1,1) + p(pst(1,id)) !- p(pst(1,1)) 
    if (sym) then 
      e(1,1)=e(1,1) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*sb
    endif 
    id = id +1 !2
    ! E block 
    e(2,2) = e(2,2)+ p(pst(1,id))! - p(pst(1,1)) 
    e(3,3) = e(3,3)+ p(pst(1,id))! - p(pst(1,1)) 
    if (sym) then 
      e(2,2)=e(2,2)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2
      e(3,3)=e(3,3)+p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2
    endif 

    id = id +1 ! 3 

    !A1'
    e(4,4) = e(4,4) + p(pst(1,id))! - p(pst(1,1)) 
    if (sym) then 
      e(4,4)=e(4,4) + p(pst(1,id)+1)*ss +p(pst(1,id)+2)*ss**2 &
                    + p(pst(1,id)+3)*ss**3
    endif 

    ! W and z of JT block 

    ! first order 
    id = id +1 !4 

    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 ! 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)
    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)

    ! the coupling of A2' and E'
    ! order  1
    id = id +1 !6

    e(1,2) = e(1,2) + b*(p(pst(1,id))*xs + p(pst(1,id)+1)*yb)
    e(1,3) = e(1,3) - b*(p(pst(1,id))*ys + p(pst(1,id)+1)*yb)
    ! order 2 
    id =id +1  ! 7
    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 A1' and E' 
    id = id +1 ! 8
    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
    if (sym) then 
      e(2,4) = e(2,4)+ p(pst(1,id)+2)*xs*ss
      e(3,4) = e(3,4)- p(pst(1,id)+2)*ys*ss
    endif

    ! order 2 
    id = id +1 ! 9
    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) 
    ! the couplign between  A
    id = id +1 ! 10 
    e(1,4) = e(1,4) + b*p(pst(1,id))
    !write(*,*) 'Id =  ',id

    call copy_2_lower_triangle(e)
    if (switch) then
      write( 6,*)" EX AND EY ARE SWITCHED" 
      tmp = e(2,2)
      e(2,2) = e(3,3)
      e(3,3) = tmp

      ! the coupling 

      tmp = e(2,4)
      e(2,4) = e(3,4)
      e(3,4) = tmp
      e(4,2) = e(2,4)
      e(4,3) = e(3,4) 
    endif 
    end subroutine potential_nh3_low 


    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