[BUG] Gromacs `FreeEnergyMinimisation` and `FreeEnergyEquilibration` failing when using `restraint`

[mb2055]

Describe the bug
Minimising or equilibrating a perturbable system using a protocol defined by either BSS.Protocol.FreeEnergyEquilibration or BSS.Protocol.FreeEnergyMinimisation and run using BSS.Process.Gromacs will consistently fail when the restraint keyword is passed.

To Reproduce
Create a perturbable system (in this case the example from the BioSimSpace tutorials) and attempt to run with restraints:

import BioSimSpace as BSS

ethane = BSS.Parameters.gaff("CC").getMolecule()

methanol = BSS.Parameters.gaff("CO").getMolecule()

mapping = BSS.Align.matchAtoms(ethane, methanol)

ethane = BSS.Align.rmsdAlign(ethane, methanol, mapping)

merged = BSS.Align.merge(ethane, methanol, mapping)

vac = merged.toSystem()

minimise = BSS.Protocol.FreeEnergyMinimisation(restraint="heavy")

process_m1 = BSS.Process.Gromacs(vac, minimise, work_dir="Minimise_GROMACS")

process_m1.start()
process_m1.wait()
print(process_m1.getSystem())

this should raise the following exception:

RuntimeError: Unable to generate GROMACS binary run input file.
'gmx grompp' reported the following warnings:
WARNING 1 [file posre_0001.itp, line 3]:
  Some parameters for bonded interaction involving perturbed atoms are
  specified explicitly in state A, but not B - copying A to B

Use 'ignore_warnings' to ignore warnings.

Removing restraint="heavy" completely removes the issue.

(please complete the following information):

• OS: Ubuntu 22.04.5 LTS
• Version of Python: 3.12.4
• Version of BioSimSpace: 2024.3.0.dev+18.gdc6602e4 (latest version of devel)
• I confirm that I have checked this bug still exists in the latest released version of BioSimSpace: yes

Additional context
This issue is also present when using restraints to minimise/equilibrate perturbable systems using OpenMM (just using BSS.Protocol.Minimisation/BSS.Protocol.Equilibration and BSS.Process.OpenMM), so the issue certainly appears to be with the restraint logic, rather than any specific protocol.

In the case of openMM the issue is specifically with the following line: parmed.load_file('openmm.prm7', 'openmm_ref.rst7'), raising the error cannot reshape array of size 0 into shape (0,3)

[lohedges]

I think the two issues are different.

For the GROMACS one: can you see if it persists if you use the regular versions of the protocols, e.g. Minimisation rather than FreeEnergyMinimisation. In this case it should explicitly be using the lambda = 0 end state so there should only be state A in the topology file.

For the OpenMM one: can you try manually loading the prm7 and rst7 files written by BioSimSpace (the lambda = 0 state) with Parmed to see if it gives a more specific error. These will be in the process work_dir. The restraints are completely separate and are just applied directly using the OpenMM API.

The only difference to the last run that I can think of is the changed LJ sigma parameters, but I can't think why they would cause this issue.

[lohedges]

It might be that the restraint file needs parameters for both states. Presumably this is done in the sandpit, so maybe take a look there. If this is the case, this is another reason why the system doesn't work since the change should have also been raised and PR'ed to the core.

[mb2055]

Confirming that the process runs successfully when using Protocol.Minimisation with GROMACS, both with and without restraint.

Re-trying everything using functionality from the sandpit doesn't solve the problem, it raises the same error as above. Also, just using the regular Minimisation within the sandpit causes an error, but a completely different one (UnboundLocalError: cannot access local variable 'total_mass' where it is not associated with a value). This is assuming that sandpit functionality is the same as core and doesn't require any extra steps.

[mb2055]

Moved the openMM issue to a separate thread: #340

[lohedges]

Thanks for confirming. Could you check the GROMACS docs to see how restraints are handled for multi state topologies. It might be as easy as adding an additional column with the force constant for state B.

[mb2055]

Looking at the docstrings it seems that the positional restraints need to be explicitly defined for states A and B.

So the .itp file that defines the positional restraints needs the following columns
i funct fcx_A fcy_A fcz_A fcx_B fcy_B fcz_B,
currently it just has
i funct fcx fcy fcz.

This should be an easy fix since, as far as I can tell, for our use case the A and B columns can just be identical.

I'll take a look at implementing a fix and a unit test tomorrow.

[lohedges]

Great, that's nice and easy then. Yes, same for A and B (restraints aren't perturbed) and only for protocols that are instances of FreeEnergyMixin.