SUBROUTINE MLR_Potential(rau, Potential)
USE Numeric_Kinds_Module
USE Li2_a3Sigma_MLR
IMPLICIT NONE
INTEGER i
REAL(KIND=WP_Kind) r, rau, alpha, Potential, ULR_Pot, ULR_Pot_re, ULR_Beta_term, ULR_re_term, ULR_y_term
REAL(KIND=WP_Kind) D6, D8, D10, Beta_Infinity, Deau, ypterm, yqterm

!Convert the radial distance from au to Angstroms
r=rau*0.529176
!Calculate the long range potential at r
D6=(1.d0-exp(-(bds*rho*r+cds*(rho*r)**2)/6))**5
D8=(1.d0-exp(-(bds*rho*r+cds*(rho*r)**2)/8))**7
D10=(1.d0-exp(-(bds*rho*r+cds*(rho*r)**2)/10))**9
ULR_Pot=(d6*c6/r**6+d8*c8/r**8+d10*c10/r**10)
ULR_Pot=ULR_Pot

!Calculate the long range potential at re
D6=(1.d0-exp(-(bds*rho*re+cds*(rho*re)**2)/6))**5
D8=(1.d0-exp(-(bds*rho*re+cds*(rho*re)**2)/8))**7
D10=(1.d0-exp(-(bds*rho*re+cds*(rho*re)**2)/10))**9
ULR_Pot_re=(d6*C6/re**6+d8*C8/re**8+d10*C10/re**10)
ULR_Pot_re=ULR_Pot_re

!Calculate the y_term
ULR_y_Term=(r**p-re**p)/(r**p+re**p)
ypterm=(r**p-rref**p)/(r**p+rref**p)
yqterm=(r**q-rref**q)/(r**q+rref**q)

!Calculate the Beta_Term
Beta_Infinity=log(2.d0*de/ULR_Pot_re)
ULR_Beta_term=Beta_Infinity*ypterm
DO i=0,3
    ULR_Beta_term=ULR_Beta_term+(1.d0-ypterm)*Beta(i)*yqterm**i
ENDDO

Potential=De*(1.d0-(ULR_Pot/ULR_Pot_re)*EXP(-ULR_Beta_term*ULR_y_Term))**2

!Convert the potential energy from inverse centimeters to bohr
Potential=Potential/2.1947e+05

RETURN
ENDSUBROUTINE MLR_Potential