Recreate P(Wf) & P(-Wr) from the raw data

GROMACS version:2023.3
GROMACS modification: No
“Can I recreate P(Wf) & P(-Wr) by myself from the raw data?”
When I originally run the alchemical simulation (solute, solvent ,and ions), I get dgdl.xvg file for nonequilibium free energy calculations from a starting state A (I just choose the first one) to state B.
gmx mdrun … -dhdl dgdl.xvg
These individual xvg file contains time axis and dU/dl(?) in the unit of kJ/mol.

The way I understand is that we compute W by
The P(wf) in the forward process range around -27kJ/mol.

I somehow get the close number by simply taking the average of the 2nd column in the dgdl.xvg.
If I believe the formula is the way to compute W in the program, then I can use the system energy(U) stored in edr file at each step (80ps 1fs step) and compute W by using above formula. Unfortunately, the number looked completely off(these dgdl.xvg & energy.xvg files are taken from the same simulation). This suggest that I am misunderstanding the calculation of the W term. What is the second column in the dgdl.xvg? Is there a way to manually get the W from edr file? Thank you very much for the clarification.