!========================================================
program         anharmonic_elevels
!========================================================
 use    , intrinsic    :: iso_fortran_env
 implicit none
 integer, parameter    :: dp = real64
 !-------------------------------------------------------
 integer               :: DIM
 real(dp), allocatable :: H(:,:), X(:,:), X4(:,:)
 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),X(DIM,DIM),X4(DIM,DIM),E(DIM))

 call calculate_X4
 call calculate_H
 call calculate_evs

 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                  :: LDA ,INFO,i,j
  !------------------------------------------------------
  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'
  !------------------------------------------------------
  !Print eigenvectors:
  print *,'# ********************** EVEC ***************'
  do j=1,DIM
   print '(A,20000G20.6)',' # EVEC ',lambda,H(:,j)
  enddo
  print *,'# ********************** EVEC ***************'
 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
  !------------------------------------------------------
  print *,'# ********************** H ******************'
  do j    = 1, DIM
   print '(A,20000G20.6)',' # HH ',H(:,j)
  end do
  print *,'# ********************** H ******************'

 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)                 :: X2(DIM,DIM)
  !------------------------------------------------------
  !Position operator:
  X(:,:) = 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))
   ! 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))
  end do
 !------------------------------------------------------
 !X2 = X . X  operator:
  X2 = MATMUL(X ,X )
 !X4 = X2. X2 operator:
  X4 = MATMUL(X2,X2)
 end subroutine calculate_X4

!========================================================
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/>.
!  -----------------------------------------------------------------------