Question on the usage of multiconformational seeding

[mb2055]

Hi Lester,

I'm looking to use the perturbed_system option to seed some calculations with a pair of end states, but it's not clear to me from reading the code how exactly the option should be used.
I can see here that there's a number of atoms consistency check with the base system, and here that there's a link_to_perturbed and link_to_reference. Does this mean that the system given in the perturbed_system flag should be the exact reverse of the original system?

My specific scenario is that I have system1 and system2 (assume that they have different numbers of atoms), both of which have been minimised and equilibrated, and I now want to run production using the coordinates of system1 for lambda=0->0.5 and the coordinates of system2 for lambda=0. 5->1.0.
Should the seeding workflow be:

1. Align system2 to system1 then merge, taking note of the mapping. Call this merged_system_A.
2. Align system1 to system2 then merge using the exact reverse of the mapping from stage1. Call this merged_system_B.
3. Run Somd2 with merged_system_A as the primary system and merged_system_B as the perturbed_system

Or am I over-complicating things?

[lohedges]

Hi there,

Apologies this isn't clearer. The current process is just to work with the Cresset pipeline but could be made simpler. To be honest, we only need an array of coordinates, so I might make an option for this.

The assumption is that you are passing input that might have been pre-equilibrated with BioSimSpace or with Sire using the swap end states functionality. As such, the merged end states will be the normal (forward) way round, but the coordinates1 property of the ligand will specify the current conformer that you want to use, not coordinates0. Cresset equilibrate externally then set the required properties.

[lohedges]

The link_to_perturbed and link_to_reference is definitely confusing. I'll double check tomorrow, but I think this is because of the different setup pathway for repex and non repex simulations. One needs to use the regular system and copy across the positions, the other can use either system. (And it not really being as useful for non-repex anyway.)

[lohedges]

Hi again. I've checked the Cresset setup script again and this is the logic:

• Multiconformational seeding is only used for replica exchange.
• The setup system has the correct ligand coordinates as the coordinates1 property.
• In the runner setup, we get the coordinates of the perturbed system, copy them into lambda=0 state with setCoordinates, then update the box to match.
• The properties are then re-linked to the lambda=0 (reference) state, which now has the lambda=1 coordinates for the ligand.
• The perturbed system is then used to create a dynamics object for all lambda windows where lambda > 0.5

[mb2055]

Okay, I think I understand. Essentially the coordinates1 property of perturbed_system should contain the coordinates for seeding lambda=0.5->1.0? Assuming that all of the required conditions (repex, equal numbers of atoms etc.) are met.

[lohedges]

That's right. This is what you'd end up with if you setup the calculation using BioSimSpace in the forward direction, then ran standard minimisation and equilibration (not FreeEnergy ones) using the is_lambda=1 option. The issue with setting things up in reverse is that the atom order in the merged ligand is different, so you'd need to re-order to match the forward system. This is actually what Cresset end up doing, since they do generate both directions independently, then combine for the production stage. It's unfortunate that this option is a bit specific to their use case, since it would be easier to just pass in a numpy array for the positions and a matrix for the box vectors.

[mb2055]

Great, thanks for the help!