!------------------------------------------------------------------------------
SUBROUTINE calc_pot !(vpot)
USE Numeric_Kinds_Module
USE CommonInfo_Module
USE FileUnits_Module
USE APH_Module
USE Aph_Ham_Module
USE Masses_Module
USE Convrsns_Module
USE Cheby_Module
IMPLICIT NONE
!==============================================================================
LOGICAL   :: const_rho, const_chi, const_theta
INTEGER   :: irho, itheta, ichi, icase
REAL (dp) :: rho_min, rho_max, rho, chi,theta, upot
REAL (dp) :: theta_min, theta_max, chi_min, chi_max, vmin
REAL (dp) :: temp(ntheta,ntheta)!, vpot (nrho, ntheta, nchi+1)

NAMELIST / pot2 / rho_min,rho_max,const_rho,rho,const_chi,chi,const_theta,theta

!==============================================================================
! Read in atom masses and APH coordinate information

OPEN (unit = In_Unit, file =TRIM(inputfile2))
READ (In_Unit, mass)
CLOSE (In_Unit)

OPEN (unit = In_Unit, file =TRIM(inputfile2))
READ (In_Unit, pot2)
CLOSE (In_Unit)

!==============================================================================
! Allocate APH grid and potential arrays

ALLOCATE(rho_val(nrho), theta_val(ntheta), chi_val(nchi+1))
ALLOCATE(v_pot(nrho,ntheta,nchi+1))

!==============================================================================
! Calculate mass dependent parameters

CALL massfac (.true.)

!==============================================================================
! Output to log and print to screen

 WRITE(*,20) 'amass',amass/amutoau
 WRITE(*,20) 'bmass',bmass/amutoau
 WRITE(*,20) 'cmass',cmass/amutoau

 WRITE(Out_Unit,21) 'Mass Info:','amu','au'
 WRITE(Out_Unit,22) 'amass', amass/amutoau, amass
 WRITE(Out_Unit,22) 'bmass', bmass/amutoau, bmass
 WRITE(Out_Unit,22) 'cmass', cmass/amutoau, cmass
 WRITE(Out_Unit,23) 'reduced mass', usys/amutoau, usys

 20 FORMAT(1X,a,10('.'),F10.6,' au')
 21 FORMAT(/1x,a,T21,a,T34,a,/1X,50('-'))
 22 FORMAT(1x,a,10("."),F10.6,3x,F9.4)
 23 FORMAT(1x,a, 3("."),F10.6,3x,F9.4)

!==============================================================================
! Calculate APH grid values

! RHO GRID=========
IF (const_rho) THEN
   nrho = 1
   rho_val (1) = rho
ELSE
   deltarho = (rho_max - rho_min) / (nrho - 1)
   DO irho = 1, nrho
      rho_val (irho) = (irho - 1) * deltarho + rho_min
   ENDDO
ENDIF

! THETA GRID=======
theta_min = 0.d0
theta_max = 2.d0 * atan (1.d0)
IF (const_theta) THEN
  ntheta = 1
  theta_val (1) = theta
ELSE
  CALL DVR_kinetic(ntheta,temp,theta_val)
ENDIF

! CHI GRID=========
chi_min = 0.d0
chi_max = 8.d0 * atan (1.d0)
IF (const_chi) THEN
   nchi = 1
   chi_val (1) = chi
ELSE
   deltachi = (chi_max - chi_min) / (nchi)
   DO ichi = 1, nchi
      chi_val (ichi) = (ichi - 1) * deltachi
   ENDDO
ENDIF

! POTENTIAL ENERGY=
DO ichi = 1, nchi+1
   chi = chi_val (ichi)
   DO itheta = 1, ntheta
      theta = theta_val (itheta)
      DO irho = 1, nrho
         rho = rho_val (irho)
         CALL potaph (rho, theta, chi, upot)
         !IF(upot.gt.vcut) upot=vcut
         !upot=upot*((15.d0/(4.d0*usys2))*rhom2 (irho))
         v_pot (irho, itheta, ichi) = upot
      ENDDO
   ENDDO
ENDDO


!==============================================================================
! Output to log and print to screen

 WRITE(Out_Unit,31) 'APH Grid Info:','au'
 WRITE(Out_Unit,33) 'nrho  ', nrho,  'rhomin  ',rho_val(1)
 WRITE(Out_Unit,33) 'ntheta', ntheta,'rhomax  ',rho_val(nrho)
 WRITE(Out_Unit,33) 'nchi  ', nchi,  'deltarho',deltarho
 WRITE(Out_Unit,34) 'ndim  ', ndim,  'deltachi',deltachi

 31 FORMAT(/1X,a,T44,a,/1x,50('-'))
 33 FORMAT(1X,a,9('.'),I4,6X,a,7('.'),f6.3)
 34 FORMAT(1X,a,9('.'),I7,3X,a,7('.'),f6.3)

RETURN
END SUBROUTINE calc_pot