[BUG] Residues in a topology and coordinates saved in a different order to input structure

[jasmin-guven]

Describe the bug
Hi! I'm trying to setup an MD system with custom residues and parameters. When I use BioSimSpace to setup the MD directory for a minimisation, some of the residues in the prm7 and pdb files produced are in a different order to the input.

To Reproduce
Steps to reproduce the behavior:

1. Load the files vim2_ligand_11_solv.prmtop and vim2_ligand_11_solv.inpcrd
2. Run the code:

import BioSimSpace as bss
import MDAnalysis as mda

system = bss.IO.readMolecules(
    [
        "vim2_ligand_11_solv.prmtop",
        "vim2_ligand_11.solv.inpcrd"
    ]
)

bss.IO.saveMolecules("test_bss.pdb", system, fileformat="PDB")

I tested with MDAnalysis like this:

universe = mda.Universe(
    "vim2_ligand_11_solv.prmtop",
    "vim2_ligand_11_solv.inpcrd",
    topology_format="PARM7"
)
universe.atoms.write("test_mda.pdb")

4. This is what goes wrong:
   The MH1 residue (reparameterised hydroxide) has a residue number 236, but is placed in between zinc 234 and the ligand ML1 with residue number 235. This causes issues down the line where I use Amber restraints (based on atom indexing) on the prm7 and rst7 generated by BioSimSpace in bss.Process.Amber

Expected behavior
Residues are not reordered from inputs.

Input files
bug.zip

(please complete the following information):

• OS: MacOS (Tahoe 26.3.1)
• Version of Python: 3.12.13
• Version of BioSimSpace: 2025.4.0
• I confirm that I have checked this bug still exists in the latest released version of BioSimSpace: [no]

I'm hoping I'm not missing anything obvious, thanks in advance!

[lohedges]

Thanks for reporting. I've transferred this across to Sire since it's most likely to be an issue with one of the parsers. I'll take a look and report back.

[jasmin-guven]

Thanks Lester, that makes sense!

[lohedges]

Okay, this is an issue with the AmberPrm parser. What's going on is that the bond graph merges ZN1 (protein molecule) with MH1 (separate molecule) because of the ZN1-MH1 covalent bond. I'll try to figure out how we can fix this in a clean way. I'm not sure if there are cases where you'd want things to be re-ordered this way, in which case, you might want to re-index based on the loaded system, not the original AMBER topology.

[lohedges]

The problem is that Sire has no concept of intermolecular bonds, so needs to pull all atoms/residues that are bonded into the same molecule. Each molecule in a system has a set of properties that define the force field parameters, e.g. bond, angle, dihedral, etc, and all of these are intra molecule.

In light of this, I've added a warning to the parser which triggers when residues are re-ordered. For example:

In [1]: import BioSimSpace as BSS

In [2]: system = BSS.IO.readMolecules(["bug/vim2_ligand_11_solv.inpcrd", "bug/vim2_ligand_11_solv.prmtop"])

WARNING: One or more residues have been reordered relative to the original topology when loading this file. This happens when covalent bonds cross AMBER's molecule boundaries (e.g. a metal ion bonded to a small molecule ligand), which forces Sire to group the bonded atoms into the same molecule. The following residues appear out of residue-number order in the loaded system: ML1(235). To avoid ordering issues, access residues by name rather than by position (e.g. mols.residues()["resname ML1"] instead of mols.residues()[1]).

In [3]: system.getIndex(system.search("resname ML1").residues()[0])
Out[3]: 235

How are you defining your restraints downstream?

[jasmin-guven]

Thanks Lester, the warning is definitely useful.

The distance restraints are defined and written out as atom indices in a restraints RST file, which is required by Amber:

&rst iat=3451,3483, r1=1.08, r2=1.58, r3=2.58, r4=3.08, rk2=48.5, rk3=48.5, /
&rst iat=3451,1318, r1=1.35, r2=1.85, r3=2.85, r4=3.35, rk2=20.0, rk3=20.0, /

I generate these restraints from MDAnalysis atom groups and my universe always matches the BioSimSpace/Sire input system being supplied into the BioSimSpace process. I could write up a work-around where I use parmed to create the restraints, but I think even then I'm not sure I would be able to avoid the reorganisation by Sire.

[lohedges]

Ah, I see. What happens if you instead create the universe from the toplogy/coordinate file written by Sire/BioSimSpace? Do the indices then agree? If they do, then this might be simplest workaround to make sure that things are consistent.

[jasmin-guven]

I can't think of a way to do that as the topology and coordinates for the process do not get generated until after the process is started (as far as I'm aware). This is what I tried in my meze code :

    def _run(
            self,
            system: Optional[bssSystem], 
            recipe: MezeRecipe,
            protocol: bss.Protocol,
            process_name: Optional[str] = "meze-run",
            config_options: Optional[dict] = None,
            namelist_options: Optional[list] = None,
            is_gpu: bool = True,
            distance_restraints: Optional[list] = None,
            distance_write_frequency: Optional[int] = 100,
    ):
        input_system = system or self.system        
        run_directory = os.path.join(recipe.workdir, process_name)
        os.makedirs(run_directory, exist_ok=True)

        namelist_options = namelist_options or []

        process = bss.Process.Amber(
            system = input_system,
            protocol = protocol,
            work_dir=run_directory,
            name=process_name,
            extra_options=config_options,
            extra_lines=namelist_options,
            is_gpu=is_gpu,
            exe=recipe.path_to_engine
        )
        self.topology = os.path.join(
            run_directory,
            process_name,
            f"{process_name}.prm7"
        )
        self.coordinates = os.path.join(
            run_directory,
            process_name,
            f"{process_name}.rst7"
        )
        self.universe = mda.Universe(
                    self.topology,
                    self.coordinates,
                    topology_format="PARM7"
        )

        if self.recipe.model == 0:
            coordination_restraints = self._prepare_distance_restraints()
            angle_restraints = self._prepare_angle_restraints()
            distance_restraints = coordination_restraints + angle_restraints
        
        if self.restraint_file and os.path.isfile(self.restraint_file):
            with open(self.restraint_file, "r") as file:
                added_distance_restraints = file.readlines()
            distance_restraints = (distance_restraints or []) + added_distance_restraints

        if distance_restraints:
            self.write_restrained_atoms_pdb(
                output_path=f"{run_directory}/{process_name}_restrained_atoms.pdb",
                restraints=distance_restraints
            )
            config_file = process._config_file
            restraint_file = os.path.join(recipe.workdir, "restraints.RST")

            with open(restraint_file, "w") as file:
                file.writelines(distance_restraints)
            
            step_restraint_file = os.path.join(run_directory, "restraints.RST")
            shutil.copyfile(restraint_file, step_restraint_file)
            
            with open(config_file, "a") as file:
                file.write(f"&wt TYPE='DUMPFREQ', istep1={distance_write_frequency} /\n")

            if not namelist_options:
                with open(config_file, "a") as file:
                    file.write(f"&wt TYPE=\"END\", /\n")
        
            with open(config_file, "a") as file:
                file.write("\n")
                file.write(f"DISANG=restraints.RST\n")
                file.write(f"DUMPAVE=distances.out\n")


        process.start()
        process.wait()

        new_system = process.getSystem()
        topology, new_coordinates = bss.IO.saveMolecules(
            f"{run_directory}/next", system=new_system, fileformat=["prm7", "rst7"]
        )

        return dataclasses.replace(
            self,
            topology=topology,
            coordinates=new_coordinates,
            recipe=recipe
        )

I get a FileNotFoundError when I try to create the MDAnalysis universe with self.universe = mda.Universe(...)

I also tried to do it in a test notebook but now for some reason I'm not able to read in the inpcrd file:

system = bss.IO.readMolecules(
    [
        "vim2_ligand_11_solv.prmtop",
        "vim2_ligand_11.solv.inpcrd"
    ]
)
protocol = bss.Protocol.Minimisation(
    steps=1000
)

process = bss.Process.Amber(
    system = system,
    protocol = protocol,
    work_dir="bug",
    name="test"
)

And I get the error:

╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ in <module>:11                                                                                   │
│                                                                                                  │
│    8 │   steps=1000                                                                              │
│    9 )                                                                                           │
│   10                                                                                             │
│ ❱ 11 process = bss.Process.Amber(                                                                │
│   12 │   system = system,                                                                        │
│   13 │   protocol = protocol,                                                                    │
│   14 │   work_dir="bug",                                                                         │
│                                                                                                  │
│ /Users/af25016/miniforge3/envs/obss-dev/lib/python3.12/site-packages/BioSimSpace/Process/_amber. │
│ py:119 in __init__                                                                               │
│                                                                                                  │
│    116 │   │   from ..Protocol._free_energy_mixin import _FreeEnergyMixin                        │
│    117 │   │                                                                                     │
│    118 │   │   # Call the base class constructor.                                                │
│ ❱  119 │   │   super().__init__(                                                                 │
│    120 │   │   │   system,                                                                       │
│    121 │   │   │   protocol,                                                                     │
│    122 │   │   │   reference_system=reference_system,                                            │
│                                                                                                  │
│ /Users/af25016/miniforge3/envs/obss-dev/lib/python3.12/site-packages/BioSimSpace/Process/_proces │
│ s.py:207 in __init__                                                                             │
│                                                                                                  │
│    204 │   │   for mol in system.getMolecules():                                                 │
│    205 │   │   │   if not mol.isPerturbable():                                                   │
│    206 │   │   │   │   if not mol._sire_object.has_property(prop):                               │
│ ❱  207 │   │   │   │   │   raise _IncompatibleError(                                             │
│    208 │   │   │   │   │   │   "System object contains molecules without coordinates!"           │
│    209 │   │   │   │   │   )                                                                     │
│    210                                                                                           │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
IncompatibleError: System object contains molecules without coordinates!

[lohedges]

It looks like the issue is that your run_directory already includes the process name, e.g.:

        input_system = system or self.system        
        run_directory = os.path.join(recipe.workdir, process_name)
        os.makedirs(run_directory, exist_ok=True)

Then, when creating the topology and coordinate file names later, the process_name is added again:

        self.topology = os.path.join(
            run_directory,
            process_name,
            f"{process_name}.prm7"
        )
        self.coordinates = os.path.join(
            run_directory,
            process_name,
            f"{process_name}.rst7"
        )

I think you just want:

        self.topology = os.path.join(
            run_directory,
            f"{process_name}.prm7"
        )
        self.coordinates = os.path.join(
            run_directory,
            f"{process_name}.rst7"
        )

[jasmin-guven]

Oh yeah, oops 😬

I fixed that and reorganised my code around a little bit and it looks like now the correct atoms are selected. Thanks Lester!

[lohedges]

Great, glad that worked. I'll close this for now since I think the warning and workaround should be sufficient. When I get time, I'll also think about doing something similar in the other topology parsers, i.e. I imagine the same thing happens for the GroTop parser too.