I think I'm in a deep mis-understanding now.

As a result of NVT calculations, we can easily draw RDF(radial distribution function),

PMF(potential of mean-force) = -kB(boltzman-constant)*T*ln(RDF). Is this formula true?

If it is, easily find PMF with a simple operation,

and then, can I think of the difference between the two energy valleys as Gibbs free energy difference(ΔG)?

When Gibbs free energy is approximated by Helmholtz energy, ΔA, can I easily get the pKa through pKa = ΔG/2.303RT ?

If my assumptions are correct, why are you doing Umbrella sampling and special(difficulty) processing?