SUBROUTINE pot(r, hstep, v0, v1, v2, jvalue, mu)
USE Numeric_Kinds_Module
USE FileUnits_Common_Module
USE FileUnits_OneDim_Module
USE Numbers_Module
USE CriticalDist_Module
USE QState_Module
USE OneDim_Parms_Module
IMPLICIT NONE
CHARACTER(LEN=1) xa
INTEGER jvalue
REAL(KIND=WP_Kind), PARAMETER :: SodiumMass=22.989767_WP_Kind , CesiumMass=132.905429_WP_Kind
REAL(KIND=WP_Kind) r, hstep, h, v0, v1, v2, mu, cent0, cent1, cent2, VNaCs

IF(SodiumMass/=MassA.or.CesiumMass/=MassB)THEN
   WRITE(Msg_Unit,*)"SodiumMass/=MassA.or.CesiumMass/=MassB"
   WRITE(Msg_Unit,*)SodiumMass, MassA, CesiumMass, MassB
   WRITE(Out_Unit,*)"SodiumMass/=MassA.or.CesiumMass/=MassB"
   WRITE(Out_Unit,*)SodiumMass, MassA, CesiumMass, MassB
ENDIF

IF(Elec==0)THEN
   xa="X"
   IF(r==0.d0)THEN
      v0=5.0d+5
      v1=-6.0d+6
      v2=2.0d+8
      RETURN
   ENDIF
ELSEIF(Elec==1)THEN
   xa="A"
   IF(r==0.d0)THEN
      v0=1.9d+2
      v1=-7.0d+2
      v2=2.5d+2
      RETURN
   ENDIF
ELSE
   STOP "Elec is not correct"
ENDIF
mu=1822.888506d0*SodiumMass*CesiumMass/(SodiumMass+CesiumMass)
v0=VNaCs(r,xa)
h=1.0d-3
v1=(VNaCs(r+h,xa)-VNaCs(r-h,xa))/h
v2=(VNaCs(r+h,xa)-Two*VNaCs(r,xa)+VNaCs(r-h,xa))/(Two*h**2)
IF(r.eq.0)THEN
   v1=-68.0d0
   v2=1.86d0
ENDIF

CALL CenTerm(r, jvalue, mu, cent0, cent1, cent2)
v0=v0+cent0
v1=v1+cent1
v2=v2+cent2

RETURN
END SUBROUTINE pot