!========================================================
program         anharmonic_elevels
!========================================================
 use    , intrinsic    :: iso_fortran_env
 implicit none
 integer, parameter    :: dp = real64
 !-------------------------------------------------------
 integer               :: DIM
 real(dp), allocatable :: H (:,:), X4(:,:)
 real(dp), allocatable :: X2(:,:), P2(:,:)
 real(dp), allocatable :: E(:)
 real(dp)              :: lambda
 !-------------------------------------------------------
 print *,'# ---------------------------------------------'
 print *,'# Enter Hilbert Space dimension:'
 read  *,              DIM
 print *,'# DIM    = ',DIM
 print *,'# Enter      lambda:'
 read  *,              lambda
 print *,'# lambda = ',lambda
 print *,'# ---------------------------------------------'

 ALLOCATE(H (DIM,DIM),X4(DIM,DIM),E(DIM))
 ALLOCATE(X2(DIM,DIM),P2(DIM,DIM))

 call calculate_X4
 call calculate_H
 call calculate_evs
 call calculate_obs

 print '(A,I8,20000G28.16)','EV ',DIM,lambda,E

!--------------------------------------------------------
contains
!--------------------------------------------------------

 !=======================================================
 subroutine     calculate_evs()
  !------------------------------------------------------
  real(dp), allocatable    :: WORK(:)
  integer                  :: LWORK
  real(dp)                 :: LWORKQUERY(1)
  character(1)             :: JOBZ,UPLO
  integer                  :: INFO
  !------------------------------------------------------
  JOBZ='V';UPLO='U'
  !------------------------------------------------------
  !Query for optimal LWORK:
  LWORK=-1
  call dsyev(JOBZ,UPLO,DIM,H,DIM,E,LWORKQUERY,LWORK,INFO)
  LWORK=LWORKQUERY(1) + 1
  !------------------------------------------------------
  ALLOCATE(WORK(LWORK))
  !------------------------------------------------------
  !Calculate eigenvalues and eigenvectors:
  call dsyev(JOBZ,UPLO,DIM,H,DIM,E,WORK      ,LWORK,INFO)
  if(INFO /= 0 ) stop 'ERROR: dsyev failed'
  !------------------------------------------------------
 end subroutine calculate_evs

 !=======================================================
 subroutine     calculate_H
  !------------------------------------------------------
  integer                  :: j
  !------------------------------------------------------
  H(:,:)  = lambda * X4(:,:)

  do j    = 1, DIM
   H(j,j) = H(j,j) + j - 0.5_dp
  end do
  !------------------------------------------------------

 end subroutine calculate_H

 !=======================================================
 subroutine     calculate_X4
  !------------------------------------------------------
  integer                  :: i, j, n, m
  real(dp), parameter      :: isqrt2=1.0_dp/sqrt(2.0_dp)
  real(dp)                 :: X(DIM,DIM), iP(DIM,DIM)
  !------------------------------------------------------
  !Position operator:
  X (:,:) = 0.0_dp
  iP(:,:) = 0.0_dp
  !Compute nonzero elements:
  do i = 1, DIM
   n=i-1 !indices 0,...,DIM-1
   ! The delta_{n,m+1} term, i.e. m=n-1
   m=n-1 !the energy level n -> i=n+1, m-> j=m+1
   j=m+1
   if(j >= 1  ) X (i,j) =  isqrt2*sqrt(DBLE(m+1))
   if(j >= 1  ) iP(i,j) = -isqrt2*sqrt(DBLE(m+1))
   ! The delta_{n,m-1} term, i.e. m=n+1
   m=n+1
   j=m+1
   if(j <= DIM) X (i,j) =  isqrt2*sqrt(DBLE(m))
   if(j <= DIM) iP(i,j) =  isqrt2*sqrt(DBLE(m))
  end do
 !------------------------------------------------------
  X2 =  MATMUL(X ,X )
  X4 =  MATMUL(X2,X2)
  P2 = -MATMUL(iP,iP)
 end subroutine calculate_X4

 !=======================================================
 subroutine     calculate_obs
  !------------------------------------------------------
  real(dp)                 :: OBS(DIM,DIM)
  real(dp),dimension(DIM)  :: avX2,avP2,avX4,DxDp!,norm
  integer                  :: i,j,k
  character(100)           :: fmt='(A5,I6,10000G25.15)'

 !norm  = 0.0_dp
  avX2  = 0.0_dp; avP2 = 0.0_dp; avX4 = 0.0_dp; DxDp = 0.0_dp
  do   i=1,DIM
   do  j=1,DIM
!   norm(i) = norm(i) + H(j,i)*H(j,i)
    do k=1,DIM
     avX2(i) = avX2(i) + H(j,i)*H(k,i)*X2(j,k)
     avX4(i) = avX4(i) + H(j,i)*H(k,i)*X4(j,k)
     avP2(i) = avP2(i) + H(j,i)*H(k,i)*P2(j,k)
    enddo
   enddo
  enddo
! norm = sqrt(norm)
  where( avX2 > 0.0_dp .and. avP2 > 0.0_dp) DxDp = sqrt(avX2*avP2)
! print fmt,'norm ',DIM,lambda,norm 
  print fmt,'avX2 ',DIM,lambda,avX2
  print fmt,'avX4 ',DIM,lambda,avX4
  print fmt,'avP2 ',DIM,lambda,avP2
  print fmt,'DxDp ',DIM,lambda,DxDp

 end subroutine calculate_obs
!========================================================
end  program    anharmonic_elevels
!========================================================
!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2021)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------