openbiosim / somd2 Goto Github PK
View Code? Open in Web Editor NEWSOMD2 molecular dynamics for free energy calculations
License: GNU General Public License v3.0
SOMD2 molecular dynamics for free energy calculations
License: GNU General Public License v3.0
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:
Input files and scripts: reproduce_somd2_abfe.tar.gz
feature_abfe
branch of SOMD2 (tiny modifications to allow a Boresch restraint to be passed in)cd reproduce_somd2_abfe
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.
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 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.
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:
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.
A few questions and things I ran into while setting up and running the simulations, not all of which are related to SOMD2:
BoreschRestraints
object in sire which allows you to input the equilibrium values, rather than measuring them from the atoms suppliedis_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).Thanks!
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.'
Run the provided input file with current default SOMD2 parameters.
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.
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.
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.
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,
)
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.
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:
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.
inputs.tar.gz and somd2 examples repository.
SOMD2 version: 0.1.dev169+ga642d47
BioSimSpace version: 2023.5.0.dev
Sire version: 2023.5.0.dev
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.
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.
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:
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).
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.
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"
)
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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.