At first sight your script looks ok for the FEP method.
In the last point (run 400000) \lambda + \delta\lambda will be negative, which may not make sense for some parameters. You can compute this point separately, by simulating at 0.02 and perturbing to 0.00.
If you are using 'volume yes’ (NPT) then what you’re saving to bar*.lmp may not be good because you’ll have the Boltzmann factor multiplied by V.
I’ll try to run a similar calculation.
Agilio Padua