What is the proper way to quantify PMF uncertainty calculated from gmx wham?

GROMACS version: 2021.4
GROMACS modification: No

I used gmx wham to perform PMF calculation of my umbrella sampling simulations. I supplied the option -nBootstrap 200:

gmx wham -it tpr.dat -if pullf.dat -o hist -unit kCal -nBootstrap 200

(I prepared the dat files of list of mdrun outputs and assume they are correct)

I did bootstrap to calculate the uncertainty of my PMF. Is the third column of the output file “bsResult.xvg” the standard deviation that I needed:

@ title “Average and stddev from bootstrapping”
@ xaxis label “\xx\f{} (nm)”
@ yaxis label “E (kcal mol\S-1\N)”
@TYPE xy
@TYPE xydy
3.814324e+00 0.000000e+00 0.000000e+00
3.827767e+00 9.938912e-01 7.105413e-03
3.841210e+00 -3.882630e-02 6.887042e-01
3.854653e+00 -1.222086e-01 6.931874e-01

When I parse the file and plot it using matplotlib it looks like this:

image

Should I report uncertainties this way or as SEM (i.e. by dividing by sqrt(200))?

If this is the average profile and the standard deviation associated by the calculations of gmx wham, then they should already be the correct errors and they do not need further treatment!