Numerical calculation for problem set 6, problem 3 (nuclear fermion gas)

Clear["Global`*"]

hbarSI = 1.05 10^-34 (* kg m^2 /s *) ;
mSI = 1.67 10^-27 (* kg *) ;
r0 = 1.4 10^(-13) 10^-2 ; (* m *) ;
rho[r0_] = 3/(4 Pi r0^3) ;
eF[rho_] = (hbarSI)^2/(2 mSI) ((3 Pi^2/2) rho)^(2/3) ;

jInMeV = 6.241509×(10^12)  ;

ef = eF[rho[r0]] (* J *) 

ef jInMeV 
3 ef jInMeV /5

Junk

(*mEvInGrams = QuantityForm[(QuantityUnit[1.782662`*^-27]),"LongForm"]

(* MeV *)

mProtonMeV = ParticleData["Proton", "Mass"] 
(*WolframAlpha["mass of proton",{{"UnitConversion",2},"ComputableData"}]*)
mNeutronMeV = ParticleData["Neutron", "Mass"] *)

(*WolframAlpha["mass of neutron",{{"UnitConversion",2},"ComputableData"}]*)

(*
mProton = mProtonMeV mEvInGrams  ;
mNeutron = mNeutronMeV mEvInGrams  ;

(*hbar = Quantity[ 1, "hbar"]*)
cm = Quantity[1, "cm"] ;
hbar = QuantityForm[(QuantityUnit[1]),"LongForm"]/(2 Pi) ;

eF[m_, s_] = ((hbar^2)/(2 m)) (6 Pi^2 (rho/2)/(2 s + 1))^(2/3)
eF[mProton, 1/2] // UnitSimplify
eF[mNeutron, 1/2]// UnitSimplify*)