Using OPC water models for RBFE calculations

[palominohernandez]

trafficstars

Dear all,

I have been trying to use other water models such as OPC (a 4 point model).

I thought that intuitively it should be as simple as:

rbfe_settings = RelativeHybridTopologyProtocol.default_settings()

rbfe_settings.forcefield_settings.forcefields=["amber/RNA.OL3.xml","amber/opc_standard.xml"]
rbfe_settings.solvation_settings.solvent_model="opc"

However, I get an error that says:

File /data/sw/mamba/envs/openfe_env/lib/python3.12/site-packages/pydantic/v1/main.py:384, in BaseModel.__setattr__(self, name, value)
    [382](https://vscode-remote+ssh-002dremote-002blinux.vscode-resource.vscode-cdn.net/data/sw/mamba/envs/openfe_env/lib/python3.12/site-packages/pydantic/v1/main.py:382) value, error_ = known_field.validate(value, dict_without_original_value, loc=name, cls=self.__class__)
    [383](https://vscode-remote+ssh-002dremote-002blinux.vscode-resource.vscode-cdn.net/data/sw/mamba/envs/openfe_env/lib/python3.12/site-packages/pydantic/v1/main.py:383) if error_:
--> [384](https://vscode-remote+ssh-002dremote-002blinux.vscode-resource.vscode-cdn.net/data/sw/mamba/envs/openfe_env/lib/python3.12/site-packages/pydantic/v1/main.py:384)     raise ValidationError([error_], self.__class__)
    [385](https://vscode-remote+ssh-002dremote-002blinux.vscode-resource.vscode-cdn.net/data/sw/mamba/envs/openfe_env/lib/python3.12/site-packages/pydantic/v1/main.py:385) else:
    [386](https://vscode-remote+ssh-002dremote-002blinux.vscode-resource.vscode-cdn.net/data/sw/mamba/envs/openfe_env/lib/python3.12/site-packages/pydantic/v1/main.py:386)     new_values[name] = value

ValidationError: 1 validation error for OpenMMSolvationSettings
solvent_model
  unexpected value; permitted: 'tip3p', 'spce', 'tip4pew', 'tip5p' (type=value_error.const; given=opc; permitted=('tip3p', 'spce', 'tip4pew', 'tip5p'))

I tried to go deep in the gufe rabbit hole yet it was a bit unclear.

I thought I could just set up another 4 point model such as tip4pew and run the system with it, as the opc parameters should override the tip4pew setup. Doing this, I could write the transformations to disk. However, when running the json files with openfe quickrun, I got the following output:

Done with all simulations! Analyzing the results....
Here is the result:
        dG = None ± None

Error: The protocol unit 'LIG1 to LIG2 repeat 0 generation 0' failed with the error message:
ValueError: Residue template U with the same override level 0 already exists.

Details provided in output.

For reference, the same simulation runs properly when using the default force field parameters.

Thanks, O.

[IAlibay]

Thanks @palominohernandez - unfortunately I don't think we've ever tested out simulations with the opc models.

My initial thought here is that openmm.Modeller.addSolvent just doesn't support it properly yet (hence why you are getting a validation error when trying to set the solvent_model settings to 'opc'.

There is likely a way to fix this, but we'll need to triage exactly what is going wrong & how much work fixing it will take.

[mikemhenry]

I believe that the issue is we only allow a subset of the solvation models that https://github.com/openmm/openmmforcefields supports. Looking at http://docs.openmm.org/latest/api-python/generated/openmm.app.modeller.Modeller.html#openmm.app.modeller.Modeller.addSolvent @IAlibay is correct that addSolvent only supports the models we allow in the settings:

model (str='tip3p') – the water model to use. Supported values are ‘tip3p’, ‘spce’, ‘tip4pew’, ‘tip5p’, and ‘swm4ndp’ (polarizable).

https://github.com/OpenFreeEnergy/openfe/blob/6aaf7e30eea4fd65255fc6660923b33c2156c6e3/openfe/protocols/openmm_utils/omm_settings.py#L46

There are other ways to solvate a system, but our current OpenMM protocols use Modeller.addSolvent AFAIK

[palominohernandez]

Thanks! Do I get it from the comments that there is no easy fix at the moment, correct?

[jthorton]

Just ran into this and I think the method described by @palominohernandez to use a similar supported water model like tip4pew to solvate the system and then the opc parameters to simulate after minimisation should be fine and follows the recommendations in the OpenMM docs:

There exist many different water models, many of which are very similar to each other. This method creates preequilibrated water boxes for a limited set of water models. In most cases, they work equally well for other models that involve the same number of particles. For example, to simulate a box of TIP3P-FB water, use this method to create a box of TIP3P water, construct a System using TIP3P-FB parameters, and perform a local energy minimization to correct for the small differences between the models. Likewise, a box of TIP4P-Ew water can be used for most four site water models.

This works with tip4p-fb in the following script, I think the issue with OPC is in the XML file in openmmtools. I noticed that this error only occurs if I load another force field with a residue named U. The amber/ff14SB.xml contains a residue named U which then clashes with some parameters in the opc ions file.

import string
import click
import tempfile
import pathlib
import gufe
from openff.units import unit
import openfe
from openfe.protocols.openmm_md.plain_md_methods import PlainMDProtocol
from rdkit import Chem


def get_settings():
    """
    Utility method for getting MDProtocol settings.

    These settings mostly follow defaults but use very short
    simulation times to avoid being too much of a burden on users' machines.
    """
    settings = openfe.protocols.openmm_md.plain_md_methods.PlainMDProtocol.default_settings()
    settings.simulation_settings.equilibration_length_nvt = 1 * unit.picosecond
    settings.simulation_settings.equilibration_length = 1 * unit.picosecond
    settings.simulation_settings.production_length = 1 * unit.picosecond
    settings.forcefield_settings.forcefields=[
        "amber/ff14SB.xml",  # ff14SB protein force field
        "amber/tip4pfb_standard.xml",  # TIP4P and recommended monovalent ion parameters
        "amber/phosaa10.xml",  # Handles THE TPO
        ]
    settings.solvation_settings.solvent_model="tip4pew"
    return settings


def run_md(dag):
    """
    Run a DAG and check it was ok.

    Parameters
    ----------
    dag : openfe.ProtocolDAG
      A ProtocolDAG to execute.

    Raises
    ------
    AssertionError
      If any of the simulation Units failed.
    """
    with tempfile.TemporaryDirectory() as tmpdir:
        workdir = pathlib.Path(tmpdir)
        dagres = gufe.protocols.execute_DAG(
            dag,
            shared_basedir=workdir,
            scratch_basedir=workdir,
            keep_shared=False,
            raise_error=True,
            n_retries=0,
        )

    # If everything is good then tell the user
    assert dagres.ok()
    print("SIMULATION COMPLETE")


@click.command
@click.option(
    '--pdb',
    type=click.Path(dir_okay=False, file_okay=True, path_type=pathlib.Path),
    required=True,
    help="Path to the prepared PDB file to validate",
)
@click.option(
    '--cofactors',
    type=click.Path(dir_okay=False, file_okay=True, path_type=pathlib.Path),
    default=None,
    help="Path to the prepared cofactors SDF file (optional)",
)
def run_inputs(pdb, cofactors):
    """
    Validate input files by running a short MD simulation

    Parameters
    ----------
    pdb : pathlib.Path
      A Path to a protein PDB file.
    cofactors : Optional[pathlib.Path]
      A Path to an SDF file containing the system's cofactors.
    """
    # Create the solvent and protein components
    solv = openfe.SolventComponent()
    prot = openfe.ProteinComponent.from_pdb_file(str(pdb))

    # Store there in a components dictionary
    components_dict = {
        'protein': prot,
        'solvent': solv,
    }

    # If we have cofactors, populate them and store them based on
    # an single letter index (we assume no more than len(alphabet) cofactors)
    if cofactors is not None:
        cofactors = [
            openfe.SmallMoleculeComponent(m)
            for m in Chem.SDMolSupplier(str(cofactors), removeHs=False)
        ]

        for cofactor, entry in zip(cofactors, string.ascii_lowercase):
            components_dict[entry] = cofactor

    # Create the ChemicalSystem
    system = openfe.ChemicalSystem(components_dict)

    # Get the settings and create the protocol
    settings = get_settings()
    protocol = PlainMDProtocol(settings=settings)

    # Now create the DAG and run it
    dag = protocol.create(stateA=system, stateB=system, mapping=None)
    run_md(dag)


if __name__ == "__main__":
    run_inputs()

[IAlibay]

@jthorton can I take from this that a locally modified version of that xml works?

[jthorton]

@jthorton can I take from this that a locally modified version of that xml works?

Yeah that works locally!