Include custom restraints in the Relative Free Energy Protocol

Lately, I've been working with some proteins that are missing several large loops, which I believe are causing instabilities during the simulation [1].

Of course, one possible solution would be to model the missing loops. But I'm wondering how difficult it would be to include custom restraints in the terminal residues for cases where loops are very large or they are far from the binding site.

I haven't seen any information on this in previous issues or the documentation, but if I missed it and you could point me in the right direction, that would be very helpful.

[1] For this protein I mentioned, I was finding the infamous Potential energy is NaN error for some edges:

...
ERROR:openmmtools.mcmc:Potential energy is NaN after 19 attempts of integration with move LangevinDynamicsMove Trying to reinitialize Context as a last-resort restart attempt...
ERROR:openmmtools.mcmc:Potential energy is NaN after 20 attempts of integration with move LangevinDynamicsMove
...

I checked the ligands alignment and mapping, the forcefield_settings.hydrogen_mass setting, run plain MD at the end points, increased equilibration length, changed the small molecule FF and finally checked the forces per particle in the nan-error-logs (I didn't try reducing the timestep, as the system is big and already runs slowly). I couldn't get rid of the error but found that the atoms experiencing the largest forces (300 times larger than the average force per particle) were consistently located in the terminal residues of the unmodeled missing loops.

Start by reading the Relative Free Energy Protocol implementation and its documentation to determine how terminal residues and restraints are currently handled. Reproduce the reported Potential energy is NaN behavior from the issue context, then define the supported custom-restraint interface and verify it on systems with unmodeled loops.