I've had a quick shot at running ABFEs with SOMD2. The overall experience was a massive step up from SOMD! The recent sire documentation was also super clear and comprehensive. However, I've had a couple of issues:

• Reproducing the results for the SOMD1 free leg
• Crashes at lambda = 1 during the free leg

Overall Set-Up

• The system is T4L lysozyme/ benzene, chosen since this is a standard ABFE test system, and I can get good results with 200 ps simulations with SOMD1. The input is the same as used for my SOMD1 calculations, and has been equilibrated for 5 ns.
• I've created the merged systems by mocking RBFE - I created a copy of the ligand, deleted the charge and LJ properties, and set the atom types to du. I then followed the typical RBFE procedure in BSS. I did this because a) SOMD2 checks for perturbable, and not decoupled molecules, and b) I ran into errors relating to missing properties if I just added the "is_perturbable" property to the decoupled ligand (prepared as standard for ABFE through BSS). See the relevant script in the input files below.
• I've used exactly the same Boresch restraints for both the SOMD1 and SOMD2 calculations, which were generated using the BSS algorithm.

To Reproduce

Input files and scripts: reproduce_somd2_abfe.tar.gz

1. Install my feature_abfe branch of SOMD2 (tiny modifications to allow a Boresch restraint to be passed in)
2. Uncompress input files above and cd reproduce_somd2_abfe
3. Run the free leg with python run_somd2.py free input/pert_system_free.bss, or bound with python run_somd2.py bound input/pert_system_bound.bss

Optionally, recreate the perturbed systems from the SOMD1 input with:

cd input
python prep_pert_systems.py free free.prm7 free.rst7 
python prep_pert_systems.py bound bound.prm7 bound.rst7

Because I ran my tests before the issue with installing the 2024 conda packages was fixed, I'm using:

biosimspace               2023.5.0                py311_0    openbiosim
sire                      2023.5.1        py311h0666ab4_0    openbiosim
python                    3.11.7          hab00c5b_1_cpython    conda-forge

Running on Ubuntu 22.04.

SOMD1 Results

Experimental free energy: -5.19 +/- 0.16 kcal mol-1 (https://pubs.acs.org/doi/epdf/10.1021/ct060037v)

All uncertainties are 95 % t-based CIs from 5 replicate runs.
Used default lambda values shown here
Used these SOMD settings. 12 A cutoff with RF.
Analysed with SOMD "native" MBAR script, no subsampling.

Free energy of releasing restraint to standard state volume: -7.08 kcal mol-1
Symmetry correction for restraint: - kTln2 = -0.41 kcal mol-1 (as rotation around axis of 6-fold symmetry occurs as the restraints are turned on, but flipping of the ring is not sampled).

200 ps per window, no equilibration

Free: -5.07 +/- 0.08 kcal mol-1
Bound: 6.86 +/-0.30 kcal mol-1
Overall = free - bound - restraint_corr + symmetry_corr = -5.07 - 6.86 + 7.08 - 0.41 = -5.26 +/- 0.35 kcal mol-1

In good agreement with experiment.

30 ns per window, 10 ns of which discarded to equilibration

Free: -5.18 +/- 0.02 kcal mol-1
Bound: 5.52 +/- 0.80 kcal mol-1
Overall: -4.03 +/- 0.80 kcal mol-1

Larger error and drift to more positive free energies due to water starting to move into the binding site for some of the runs towards lambda = 1 in the bound vanish leg.

SOMD2 Results

Bound Leg

SOMD2 defaults other than RF with 12 A cutoff and 100 ps runs

Repeat 1: 6.2 +/- 0.3 kcal mol-1
Repeat 2: 6.9 +/- 0.4 kcal mol-1
Repeat 3: 7.1 +/- 0.4 kcal mol-1

In good agreement with the SOMD1 results for 200 ps, at least within the large errors.

The overlap is more than enough:

and the ligand clearly stays in the binding site when fully decoupled.

As above but without restraints

-1.8 +/- 0.2 kcal mol-1

The with_restraints - without_restraints difference is in the expected order of magnitude.

SOMD2 defaults other than RF with 12 A cutoff, 5 ns runs, and 10 ps energy frequency

Repeat 1: 5.5 +/- 0.2 kcal mol-1
Repeat 2: 5.3 +/- 0.2 kcal mol-1
Repeat 3: 5.3 +/- 0.2 kcal mol-1

A slight decrease compared to the shorter runs, as expected, and in good agreement with the SOMD1 30 ns runs.

Free Leg

SOMD2 defaults other than RF with 12 A cutoff and 100 ps runs

Repeat 1: -9.1 +/- 0.6 kcal mol-1
Repeat 2: -7.8 +/- 0.6 kcal mol-1
Repeat 3: -8.3 +/- 0.6 kcal mol-1

Substantially more negative than for SOMD1, producing a substantially more negative free energy of binding compared to experiment. Also, substantially greater variability in results than expected. As expected from the bound leg, the overlap was good:

SOMD2 defaults other than RF with 12 A cutoff, 5 ns runs, and 10 ps energy frequency

Repeat 1: -8.8 kcal mol-1

Attempted 5 more repeats, but all failed at lambda = 1.0 with:

or

Tried:

• Dropping the timestep to 0.5 fs (still crashed)
• Running with perturbable_constraints="none" and 0.5 fs timestep (still crashed)

Investigating lambda = 1 crashes with bound leg

• Completed 5 further 5 ns runs with the bound leg at lambda = 1 - no crashes observed
• Turned off the restraints - this introduced the lambda = 1 crashes (over three repeats, one crash and two warnings as above).

So we get crashes when the ligand is non-interacting which are prevented by coupling the non-interacting ligand to interacting particles - seems like energy is building up in the decoupled ligand. To try and dissipate this energy, I tried running tried free leg 5 ns runs with a Brownian integrator - this seemed to reduce the frequency of crashes - as 4 / 6 runs completed.

Questions and Comments on Software

A few questions and things I ran into while setting up and running the simulations, not all of which are related to SOMD2:

• In sire, I got an error when trying to convert atoms in the perturbed molecule to vectors. This makes sense, but it would be helpful if the error message stated that this was why.
• Might be helpful to have a more flexible BoreschRestraints object in sire which allows you to input the equilibrium values, rather than measuring them from the atoms supplied
• How best to check the type of the restraints object in the SOMD2 config? BoreschRestraints didn't seem to be an instance of Restraints. Is there a better way than checking against all possible types of restraints?
• Properties set by BSS for ABFEs. The standard Exscientia sandpit ABFE workflow sets the ligand property is_decoupled rather than is_perturbed and sets the charges, LJ terms, and atom types as end-state 1 properties in the ligand. However, as mentioned above, SOMD2 only checks for is_perturbed and I got a sire error (not all required properties present) when I just set the is_perturbed property in the decoupled ligand returned by the standard BSS workflow for ABFE. Given that the ABFE and RBFE workflows are very similar in BSS (handled by the same class) maybe it would be be better to have the properties set in a more similar way (e.g. have is_perturbed for the ABFE ligand).
• Where should the restraint corrections be calculated? Given that sire can now hold Boresch restraints, maybe the correction calculation code should be moved from BSS to sire.
• How will complex lambda schedules and Boresch restraints be input into the SOMD2 config?

Thanks!

Describe the bug

Sodium-in-a-box (a single sodium atom being perturbed in a water box) system triggers a HMR crash, because the code currently tries to find perturbable hydrogen atoms in the system, which here don't exist.

Exact error reported is:

╭─────────────────────────────── Traceback (most recent call last) ────────────────────────────────╮
│ /home/akalpokas/mambaforge/envs/somd2_fast/bin/somd2:8 in <module>                               │
│                                                                                                  │
│   5 from somd2.app import cli                                                                    │
│   6 if __name__ == '__main__':                                                                   │
│   7 │   sys.argv[0] = re.sub(r'(-script\.pyw|\.exe)?$', '', sys.argv[0])                         │
│ ❱ 8 │   sys.exit(cli())                                                                          │
│   9                                                                                              │
│                                                                                                  │
│ /home/akalpokas/code/git/somd2/src/somd2/app/run.py:94 in cli                                    │
│                                                                                                  │
│   91 │   _logger.info(f"sire version: {sire_version}+{sire_revisionid}")                         │
│   92 │                                                                                           │
│   93 │   # Instantiate a Runner object to run the simulation.                                    │
│ ❱ 94 │   runner = Runner(system, config)                                                         │
│   95 │                                                                                           │
│   96 │   # Run the simulation.                                                                   │
│   97 │   runner.run()                                                                            │
│                                                                                                  │
│ /home/akalpokas/code/git/somd2/src/somd2/runner/_runner.py:179 in __init__                       │
│                                                                                                  │
│   176 │   │   ]                                                                                  │
│   177 │   │                                                                                      │
│   178 │   │   # Work out the current hydrogen mass factor.                                       │
│ ❱ 179 │   │   h_mass_factor = self._get_h_mass_factor(self._system)                              │
│   180 │   │                                                                                      │
│   181 │   │   # HMR has already been applied.                                                    │
│   182 │   │   from math import isclose                                                           │
│                                                                                                  │
│ /home/akalpokas/code/git/somd2/src/somd2/runner/_runner.py:609 in _get_h_mass_factor             │
│                                                                                                  │
│   606 │   │   expected_h_mass = Element("H").mass().value()                                      │
│   607 │   │                                                                                      │
│   608 │   │   # Get the hydrogen mass.                                                           │
│ ❱ 609 │   │   h_mass = system.molecules("property is_perturbable")["element H"][0].mass()        │
│   610 │   │                                                                                      │
│   611 │   │   # Work out the current hydrogen mass factor. We round to 3dp due to                │
│   612 │   │   # the precision of atomic masses loaded from text files.                           │
│                                                                                                  │
│ /home/akalpokas/mambaforge/envs/somd2_fast/lib/python3.11/site-packages/sire/mol/__init__.py:505 │
│ in __fixed__getitem__                                                                            │
│                                                                                                  │
│    502 │   │   │   # try to search for the object - this will raise                              │
│    503 │   │   │   # a SyntaxError if this is not a search term                                  │
│    504 │   │   │   # (and is instead a name)                                                     │
│ ❱  505 │   │   │   return __from_select_result(obj.search(key, map=map))                         │
│    506 │   │   except SyntaxError:                                                               │
│    507 │   │   │   # ignore SyntaxErrors as this is a name                                       │
│    508 │   │   │   pass                                                                          │
│                                                                                                  │
│ /home/akalpokas/mambaforge/envs/somd2_fast/lib/python3.11/site-packages/sire/mol/__init__.py:364 │
│ in __from_select_result                                                                          │
│                                                                                                  │
│    361 │   │   )                                                                                 │
│    362 │                                                                                         │
│    363 │   if obj.list_count() == 0:                                                             │
│ ❱  364 │   │   raise KeyError("Nothing matched the search.")                                     │
│    365 │                                                                                         │
│    366 │   typ = obj.get_common_type()                                                           │
│    367                                                                                           │
╰──────────────────────────────────────────────────────────────────────────────────────────────────╯
KeyError: 'Nothing matched the search.'

To reproduce

Run the provided input file with current default SOMD2 parameters.

Input files

inputs.tar.gz

Environment information:

SOMD2 version: somd2 version: 0.1.dev281+g32ab743.d20240319
Sire version: sire version: 2024.1.0.dev+d971cfd - compiled with OpenMM performance enhancements enabled

The option restart=True is currently non-functional, throwing the error SireBase::missing_property: None of the contained forcefields have a property called config. Available properties are [ energy_trajectory, time, space, ensemble ]. This results from SOMD2 attempting to compare the streamed checkpoint file with its current configuration (needed to ensure that the requested restart from checkpoint is sensible) and being unable to find either the old configuration or the old lambda value in the porperties of the streamed system.

Should be a simple fix - both lambda and config need to be added to the shared property of the system before it is saved to checkpoint.

Describe the bug

Occasional simulation crashes seem to be caused by what the engine is reporting as a periodic box shrinking. This crash tends to happen between λ values of 0.6-1.0 for the system provided (GLU --> GLY mutation on a tripeptide with co-alc water), after around 1ns of equilibration time. I have also seen this crash happen on two other systems (same system without co-alc water and a large 770 residue system), with variety of box sizes. I don't believe that the issue is necessarily with barostat itself, but rather the system blowing up and triggering this error somehow (suggested previously here). Currently, I am unable to reproduce this crash in sire, which is interesting since I tried to keep the simulation parameters as similar as possible.

Exact error reported is:
Minimisation/dynamics at λ = 1.0 failed with the following exception The periodic box size has decreased to less than twice the nonbonded cutoff.

To reproduce

SOMD2

Extract the provided tar.gz file and run the SOMD2 input file with:

somd2 e2g_box_30.bss --log-level debug --timestep 4fs --log-file output.log --cutoff-type rf --runtime 500ps --num-lambda 33 --cutoff 10A --equilibration-time 10000ps --equilibration-timestep 2fs --energy-frequency 1ps --frame-frequency 200ps --checkpoint-frequency 5ns

output.log file is also provided which contains the logged errors from a previous run.

Sire

import sire as sr
from sire import morph
import numpy as np

system = sr.stream.load("e2g_box_30.bss")
system = morph.link_to_reference(system)
system = morph.repartition_hydrogen_masses(system, 1.5)

# run params
timestep = "2fs"
equil_time = "10000ps"
run_time = "500ps"
cutoff_type = "rf"
energy_frequency = "1ps"
lambda_value = 0.65625 # problematic lambda value
cutoff = "10A"
constraint = "h-bonds"
perturbable_constraint = "h-bonds-not-perturbed"

# miminmisation
min_system = (
    system.minimisation(
        lambda_value=lambda_value,
        constraint=constraint,
        perturbable_constraint="none",
    )
    .run()
    .commit()
)

lambda_values = list(np.linspace(0.0, 1.0, 33))

d = min_system.dynamics(
    timestep="2fs",
    temperature="25oC",
    cutoff="10A",
    cutoff_type=cutoff_type,
    lambda_value=0.65625,
    shift_delta="2A",
    pressure="1atm",
    constraint=constraint,
    perturbable_constraint=perturbable_constraint,
)

# generate random velocities
d.randomise_velocities()

# equilibrate, not saving anything
d.run(equil_time, save_frequency=0)

lambda_index = lambda_values.index(0.65625)
lambda_windows = lambda_values[max(lambda_index - 1, 0) : min(len(lambda_values),lambda_index  + 2)]

# run the dynamics, saving the energy every 1ps
d.run(
    run_time,
    energy_frequency=energy_frequency,
    frame_frequency=0,
    lambda_windows=lambda_windows,
)

Input files

inputs.tar.gz

Environment information:

SOMD2 version: 0.1.dev261+g6e0da22.d20240307
Sire version: 2024.1.0.dev+d971cfd - compiled with OpenMM performance enhancements enabled

It would make sense to update README.md to provide instructions to set CUDA_VISIBLE_DEVICES in order to run jobs on a CUDA platform.

The logic for applying a somd1 compatible perturbation by modifying bonded terms assumes that the containers preserve order, which isn't the case. I'll just fully re-implement the pert file logic from BioSimSpace.

The current behaviour of the restart flag doesn't properly capture the scenario in which --restart is specified, but one or more checkpoint files are missing.
A warning of Unable to load checkpoint file for λ={value}, starting from scratch. is correctly raised when this scenario is encountered, however the simulation never actually starts from scratch and an error of Exception raised for λ = {value}: cannot access local variable 'system' where it is not associated with a value is raised due to the absence of a system.

Describe the bug

GPU utilisation seems to be unstable during FEP runs for the bound leg (protein + ligand + water, ~63K atoms) for a sample simulation in TYK2 benchmark. This results in slow MD simulations for this leg as the GPU usage jumps between 0 to 100% rapidly which can be seen in this nvtop output screenshot:

For the free leg (ligand + water only, ~3.5K atoms) GPU is properly utilised (drops indicate points at which the engine starts a new simulation with a new λ value:

To reproduce

Run provided input files with:

export NUMEXPR_MAX_THREADS=8
export CUDA_VISIBLE_DEVICES=0

somd2 ejm31-ejm50.bss --log-level trace

Full reproducible example also available in somd2 examples repository.

Things I have tried so far

1. Increasing the number of CPU threads to 30 - doesn't seem to change much
2. Run a multi-GPU simulation - same behaviour on each GPU
3. Run on a different machine - on my laptop the GPU utilisation is spiking between 50-100%, instead of 0-100%

Input files

inputs.tar.gz and somd2 examples repository.

Environment information:

SOMD2 version: 0.1.dev169+ga642d47
BioSimSpace version: 2023.5.0.dev
Sire version: 2023.5.0.dev

Describe the bug

Currently streaming to serialized checkpoint (.s3) files causes crashes when doing FEP on large systems due to the machine running out of memory. For example, this particular system of around 770 residues causes memory crashes on a 32GB RAM machine after around 1ns of dynamics when the frame saving frequency is 100ps. The memory usage of the SOMD2 process can also be seen to steadily climb as the simulation progresses. While it is possible to circumvent this issue by either severely decreasing the frame saving frequency or disabling checkpointing entirely, neither serve as a good long term solution to the problem.

To reproduce

Extract the provided tar.gz file and run the SOMD2 input file with:

somd2 perturbablemin_system.bss --timestep 1fs --cutoff-type rf --equilibration-timestep 1fs --equilibration-time 500ps --checkpoint-frequency 500ps --frame-frequency 100ps

Given the timestep used and the size of the system, it might take awhile to run the FEP simulation to the point where it crashes. Also I believe running this on multiple GPUs at once is more problematic than just a single one, since with multiple GPU runs there will be multiple trajectories being stored in the memory at a given time.

Input files

perturbablemin_system.tar.gz

Environment information:

SOMD2 version: 0.1.dev266+ge2e7a97
Sire version: 2024.1.0.dev+d971cfd

A warning should be raised when HMR factor is set to be below 1.0 in the configuration, since this might lead to unexpected simulation behaviour.

Equilibration currently has a few issues:

• Lack of minimisation between equilibration and production runs, with the change in timestep between the two presenting a source of instability due to bond length constraints.
• HMR is performed prior to equilibration, it should be performed after equilibration (if a 4fs timestep is set).
• Equilibration can be applied upon restart from a checkpoint file, this is likely to cause more confusion than it solves. Restarting from checkpoint should perform only a brief minimisation before continuing with production.

Describe the bug

Applying HMR and increasing the timestep to 4fs seems to crash SOMD2 simulations with TYK2 dataset when using the most recent version of SOMD2 (0.1.dev220+geb89275.d20231205). The particular example runs fine at 2fs timestep with/without HMR applied, but the dynamics crash immediately once a higher timestep is used. I was able to run these simulations with HMR and 4fs timestep using a previous version of SOMD2 (0.1.dev169+ga642d47).

Interestingly, I can use the sire faster FEP script to run same the perturbation and the input file with HMR applied at 4fs timestep in the same conda envrionment that this version of SOMD2 is installed, which makes me believe that there might be an issue with the way SOMD2 is either applying HMR or passing the system to sire with HMR applied to do the dynamics?

I tried to debug by printing out the masses before and after HMR is applied for the perturbable molecule, and I get the same masses between the sire script and SOMD2, so my guess it's the latter, although there could be something else going on that I haven't considered entirely (like the recent updates to the constraints code).

To reproduce

To test SOMD2 crashing

Run provided input file with:

somd2 lig_ejm31-lig_ejm42.bss --timestep 4fs --h-mass-factor 1.5 --cutoff-type rf --num-lambda 21 --log-level debug --runtime 500ps --equilibration-time 2ps --energy-frequency 0.1ps --no-save-trajectories

Which should crash immediately. Now reduce the timestep to 2fs, and this runs fine.

To show sire can run the system at 4fs timestep with HMR:

Run the following code, which is adapted from sire faster FEP script:

import sire as sr

mols = sr.load(f"lig_ejm31-lig_ejm42.bss")

for mol in mols.molecules("molecule property is_perturbable"):
    mol = mol.perturbation().link_to_reference().commit()
    mol = sr.morph.repartition_hydrogen_masses(mol, mass_factor=1.5)
    mols.update(mol)

mols = mols.minimisation().run().commit()

timestep = "4fs"
energy_frequency = "0.1ps"

constraint = "h-bonds"
perturbable_constraint = None

equil_time = "2ps"
run_time = "500ps"

lambda_values = [x / 100.0 for x in range(0, 101, 5)]

print("\nWater leg")

for i, lambda_value in enumerate(lambda_values):
    print(f"Simulating lambda={lambda_value:.2f}")
    # minimise the system at this lambda value
    min_mols = mols.minimisation(lambda_value=lambda_value,
                                 constraint=constraint,
                                 perturbable_constraint="none").run().commit()

    # create a dynamics object for the system
    d = min_mols.dynamics(
        timestep=timestep,
        temperature="25oC",
        lambda_value=lambda_value,
        constraint=constraint,
        perturbable_constraint=perturbable_constraint,
        cutoff_type="RF"
    )

    # generate random velocities
    d.randomise_velocities()

    # equilibrate, not saving anything
    d.run(equil_time, save_frequency=0)
    print("Equilibration complete")
    print(d)

    # get the values of lambda for neighbouring windows
    lambda_windows = lambda_values[
        max(i - 1, 0) : min(len(lambda_values), i + 2)
    ]

    # run the dynamics, saving the energy every 0.1 ps
    d.run(
        run_time,
        energy_frequency=energy_frequency,
        frame_frequency=0,
        lambda_windows=lambda_windows,
    )
    print("Dynamics complete")
    print(d)

    # stream the EnergyTrajectory to a sire save stream object
    sr.stream.save(
        d.commit().energy_trajectory(), f"energy_water_{lambda_value:.2f}.s3"
    )

Input files

inputs.tar.gz

Environment information:

SOMD2 version: 0.1.dev220+geb89275.d20231205
BioSimSpace version: 2023.5.0.dev
Sire version: 2023.5.0.dev (106d783) as far as I can tell

We currently check that no HMR has been applied to the input system so that we don't re-apply if the h_mass_factor is not 1. What we should really do is check the current HMR factor and only apply the new one if it is different. If there is a difference, then we can either invert the existing HMR and re-apply, or apply the partial factor.

e.g. when integrator=verlet is passed to the config, querying the underlying dynamics object shows that the sire default of langevin_middle is still being used.

To make the code more flexible and compatible with adaptive lambda scheduling we should add the ability to set lambda values using a list.
The simplest way to do this is to add an additional option to the config (for example lambda_values) that overwrites the evenly-spaced values generated by the num_lambda logic.