SUBROUTINE Morse_Potential(r, n, mu, v, vp, vpp)
USE Numeric_Kinds_Module
USE Numbers_Module
IMPLICIT NONE
SAVE
INTEGER IFirst, NuMax
REAL(KIND=WP_Kind) :: De   = 3.880453448111275E-002
REAL(KIND=WP_Kind) :: re   = 5.05240520316459 
REAL(KIND=WP_Kind) :: k    = 8.099321394117176E-003
REAL(KIND=WP_Kind) :: Beta = 0.323048760974807
REAL(KIND=WP_Kind) :: vasy =-1.393390000000000D-177
INTEGER, PARAMETER :: Morse_Unit=8
INTEGER :: n, i, nu
REAL(KIND=WP_Kind) :: mu, w0, wnu, Eigens, arg
REAL(KIND=WP_Kind) :: r, v, vp, vpp
DATA IFirst/-1/

IF(IFirst==-1)THEN
  w0=Sqrt(Two*k/mu)
  WRITE( Morse_Unit,*)
  WRITE( Morse_Unit,*)"Morse Parameters"
  WRITE( Morse_Unit,*)"  De=",De
  WRITE( Morse_Unit,*)"  re=",re
  WRITE( Morse_Unit,*)"   k=",k
  WRITE( Morse_Unit,*)"Beta=",Beta
  WRITE( Morse_Unit,*)"  w0=",w0
  OPEN(unit= Morse_Unit,file="E:\ParkerE\Potentials\Diatomics\MorseEigen.txt")
  numax=Two*De/w0-One/Two
  DO nu=0,numax
    wnu=w0*(One*nu+Half)
    Eigens=wnu*(One-wnu/(Four*De))+vasy-De
    WRITE( Morse_Unit,*) nu, Eigens
  ENDDO
  CLOSE(unit= Morse_Unit)
ENDIF

arg=Beta*(re-r)
v=De*(One-exp(arg))**2-De+vasy
vp=Two*Beta*De*exp(arg)*(One-exp(arg))
vpp=Two*Beta**2*De*exp(arg)*(Two*exp(arg)-One)

RETURN
END SUBROUTINE Morse_Potential