simonboothroyd / absolv

Absolute solvation free energy calculations with OpenFF and OpenMM

ABsolute SOLVantion Free Energy Calculations

Absolute solvation free energy calculations using OpenMM

---

The absolv framework aims to offer a simple API for computing the change in free energy when transferring a solute from one solvent to another, or to vacuum in the case of solvation free energy calculations.

It offers two routes to this end: standard equilibrium calculations and non-equilibrium switching type calculations, where the latter will be the main focus of this framework.

Warning: This code is currently experimental and under active development. If you are using this it, please be aware that it is not guaranteed to provide correct results, the documentation and testing is incomplete, and the API can change without notice.

Installation

This package can be installed using conda (or mamba, a faster version of conda):

mamba install -c conda-forge absolv

If you are running with MPI on an HPC cluster, you may need to instruct conda to use your local installation depending on your setup

mamba install -c conda-forge absolv "openmpi=4.1.5=*external*"

Getting Started

To get started, see the usage guide.

absolv's Issues

ani2x compatibility?

Is there an easy way to integrate ani2x for ml/mm potentials in this workflow? I'm digging around and seeing some implementations in "plain" openmm, but not sure how thatd translate here.

runner.setup() does not accept target density for system

Hi Devs,

I noticed that the solvent boxes in absolv are always prepared with the same density, independent of the chosen solvent.

I found that absolv.setup.setup_system() would generally accept a target density for the solvent box. However, absolv.runner.setup() (which is the function that calls .setup.setup_system()) does not accept this argument.
(For clarity: .setup() calls a function _setup_solvent, which then calls .setup.setup_system()` here: line67)

As a results, users may not influence the density of the box that is assembled. I think it would be great to provide the possibility to influence this target density somehow, because depending on the solvent of interest the default value may actually be quite far off.

I have two suggestions, how to do that:
Either the target densities will be additional arguments for .runner.setup(). I guess that would be somewhat consistent.
Alternatively, (and I think that would be the better option), absolv.config.System() may accepttwo additiona keywords target_dens_a and target_dens_b.

In case, I may help with testing, I would be very much open to assist.

Pymbar 4 compatibility

After the recent API update in pymbar the CI has started failing. Happy to take a go at updating absolv to work with the new API if you like.

openff.interchange should be a dependency

Hi Devs,

I know the project is in an early stage, but I would already like to try it out. I installed it from conda-forge with micromamba, which seemed to have worked pretty well and also quickly.
However, when I tried to run the example from the user-guide

I received an error when initializing the runner:
prepared_system_a, prepared_system_b = absolv.runner.setup(system, config, force_field)
failed, because openff.interchange could not be imported.

I believe openff.interchange should be included in the dependencies, right?

Issues w/ tutorial NEQ run

I've tried running the NEQ example as outline in the documentation:

from openmm import unit

from absolv.models import (
    NonEquilibriumPfrom openff.toolkit.typing.engines.smirnoff import ForceField
force_field = ForceField("openff-2.0.0.offxml")

from absolv.runners.equilibrium import NonEquilibriumRunner
NonEquilibriumRunner.setup(schema, force_field)rotocol,
    SwitchingProtocol,
    SimulationProtocol,
    State,
    System,
    TransferFreeEnergySchema
)
from absolv.runners.nonequilibrium import NonEquilibriumRunner
from openff.toolkit.typing.engines.smirnoff import ForceField

schema = TransferFreeEnergySchema(
    system=System(solutes={"CCO": 1}, solvent_a=None, solvent_b={"O": 895}),
    # Define the state that the calculation will be performed at.
    state=State(temperature=298.15 * unit.kelvin, pressure=1.0 * unit.atmosphere),
    # Define the alchemical pathway to transform the solute along in vacuum ('solvent_a')
    # and water ('solvent_b')
    alchemical_protocol_a=NonEquilibriumProtocol(
        production_protocol=SimulationProtocol(
            n_steps_per_iteration=6250,
            n_iterations=160,
            timestep=2.0 * unit.femtoseconds
        ),
        switching_protocol=SwitchingProtocol(
            n_electrostatic_steps=60,
            n_steps_per_electrostatic_step=100,
            n_steric_steps=0,
            n_steps_per_steric_step=0,
            timestep=2.0*unit.femtoseconds
        )
    ),
    alchemical_protocol_b=NonEquilibriumProtocol(
        production_protocol=SimulationProtocol(
            n_steps_per_iteration=6250,
            n_iterations=160,
            timestep=2.0 * unit.femtoseconds
        ),
        switching_protocol=SwitchingProtocol(
            # Annihilate the electrostatic interactions over the first 12 ps
            n_electrostatic_steps=60,
            n_steps_per_electrostatic_step=100,
            # followed by decoupling the vdW interactions over the next 38 ps
            n_steric_steps=190,
            n_steps_per_steric_step=100,
        )
    )
)

force_field = ForceField("openff-2.0.0.offxml")
NonEquilibriumRunner.setup(schema, force_field)

However, the above gives the error:

Traceback (most recent call last):
  File "testRun.py", line 52, in <module>
    NonEquilibriumRunner.setup(schema, force_field)
  File "/home/spadavec/itx/research_projects/absolv/absolv/runners/_runners.py", line 171, in setup
    cls._setup_solvent(
  File "/home/spadavec/itx/research_projects/absolv/absolv/runners/_runners.py", line 96, in _setup_solvent
    topology, coordinates = PACKMOLCoordinateFactory.generate(components)
  File "/home/spadavec/itx/research_projects/absolv/absolv/factories/coordinate.py", line 147, in generate
    box_size = cls._approximate_box_size_by_density(components, target_density)
  File "/home/spadavec/itx/research_projects/absolv/absolv/factories/coordinate.py", line 68, in _approximate_box_size_by_density
    volume += molecule_volume * count
  File "/home/spadavec/miniconda3/envs/absolv2/lib/python3.8/site-packages/openmm/unit/quantity.py", line 235, in __add__
    raise TypeError('Cannot add two quantities with incompatible units "%s" and "%s".' % (self.unit, other.unit))
TypeError: Cannot add two quantities with incompatible units "angstrom**3" and "mole*milliliter/gram".

Is there something I need to change? I should also note in the documentation, the NEQ example needs to change the following line(s):

from

from absolv.runners.equilibrium import NonEquilibriumRunner

to

from absolv.runners.nonequilibrium import NonEquilibriumRunner

and SwitchingProtocol is never imported from absolv.models

Cheers