_add_position_restraints change topology

[kexul]

Dear BioSimSpace Developers: Recently I found that the .top file would change after _add_position_restraints during Gromacs setup, which would make the simulation fail to proceed. https://github.com/michellab/BioSimSpace/blob/3b8f38c4feec2280e1ea89c7f16b28631b1de4d8/python/BioSimSpace/Process/_gromacs.py#L2005
Here is the code to reproduce:

import BioSimSpace as BSS
from pathlib import Path

work_dir = '.'

p0 = str(list(Path(work_dir).glob(f'*0.prm7'))[0])
p1 = str(list(Path(work_dir).glob(f'*1.prm7'))[0])
r0 = str(list(Path(work_dir).glob(f'*0.rst7'))[0])
r1 = str(list(Path(work_dir).glob(f'*1.rst7'))[0])

morph = BSS.IO.readPerturbableSystem(p0, r0, p1, r1)

sys = BSS.Solvent.tip3p(morph, [4*BSS.Units.Length.nanometer]*3, ion_conc=0.154, work_dir=f'solv')

protocol = BSS.Protocol.Minimisation()
process = BSS.Process.Gromacs(sys, protocol, work_dir='mini', ignore_warnings=True)
process.addArgs({'-ntmpi':1, '-ntomp':'8', '-gpu_id':0})
process.start()
process.wait()
sys = process.getSystem()

protocol = BSS.Protocol.Equilibration(restraint='all')  # without restraint, the simulation could proceed with no error. 
process = BSS.Process.Gromacs(sys, protocol, work_dir='equ', ignore_warnings=True)
process.addArgs({'-ntmpi':1, '-ntomp':'8', '-gpu_id':0})
process.start()
process.wait()
sys = process.getSystem()

Input files: temp.zip

Differences in [pairs]:

Additional information

This is a perturbation that has ring breaking, [16, 37] is the bond that would break during the perturbation. Though this type of perturbation is not fully supported currently, could we have an easy fix? Thanks!

[lohedges]

Hi there. I'm a bit confused as to where you are seeing differences in the topology file? I've written the topology of the minimised system and compared it to the one written by the second process when applying the equilibration protocol. The only differences I see are the time stamp and the inclusion of a position restraint file. What files are you comparing when looking at the [pairs] records?

In my case I do still see Constraint error in algorithm Lincs at step 0, so the constraint is indeed invalid given the starting topology.

[lohedges]

Hmmm, I've checked the topology with and without constraints in the equilibration protocol and it looks like it's extracting the lambda=0 state in one case, and not in the other. This should be consistent, so I'll try to work out why that's different.

[kexul]

Hi @lohedges , thanks for checking ! I just find that adding system = self._checkPerturbable(system) here would help. https://github.com/michellab/BioSimSpace/blob/38457685dae4861ee360bf1661bdddc37df77e4e/python/BioSimSpace/Process/_gromacs.py#L2033

[lohedges]

Yes, that's exactly the fix that I've just implemented :-)

[lohedges]

I've added it on my feature branch which is still to be merged, so I'd suggest adding to your local copy for now.