`PeriodicNonequilibriumIntegrator` docstring example storing work

[ijpulidos]

PeriodicNonequilibriumIntegrator seems to not be using the nsteps_eq parameter to perform any equilibration. Instead, users need to manually call integrator.step(nsteps_eq) to run an equilibration before the non-equilibrium part.

The signature of the class suggests that it should do this automatically.

[jchodera]

I'm not sure this is unexpected behavior. The PeriodicNonequilibriumIntegrator is meant to set up an integrator that cycles through the following sequence:

    eq 0 for nsteps_eq | neq 0->1 over nsteps_neq | eq 1 for nsteps_eq | neq 1->0 for nsteps_neq

However, we should update the docstring example, since it doesn't show how to correctly measure forward and reverse work.

Here's the example usage for how to do that:

# Run one periodic cycle: (eq 0) > (neq 0->1) > (eq 1) > (neq 1->0)
# Store dimensionless forward and reverse protocol work
# Run equilibration at lambda=0
integrator.step(nsteps_eq)
# Run forward process and measure dimensionless protocol work
initial_work = integrator.get_protocol_work(dimensionless=True) # in kT
integrator.step(nsteps_neq)
forward_work = integrator.get_protocol_work(dimensionless=True) - initial_work # in kT
# Run equilibration at lambda=1
integrator.step(nsteps_eq)
# Run reverse process and measure protocol work
initial_work = integrator.get_protocol_work(dimensionless=True) # in kT
integrator.step(nsteps_neq)
reverse_work = integrator.get_protocol_work(dimensionless=True) - initial_work # in kT

The signature of the class suggests that it should do this automatically.

It does run the equilibrations automatically if you use

integrator.step(2*nsteps_eq + 2*nsteps_neq)

or

integrator.step(nsteps_per_period)

However, it doesn't save the intermedate work values, which we need for free energy calculations.