Espaloma parameters are not loaded correctly.

[Kamecs]

Dear OpenFE Developer.

Thank you for your useful OpenFE code.

Recently, I am trying espaloma parameters to apply self-fine-tuned MLFF.
To begin with, I tried to use latest espaloma parameter (https://github.com/choderalab/espaloma/releases/tag/0.4.0).
As a result, following output was shown.

On the other hand, through the simple usage (e.g. https://docs.espaloma.org/en/latest/deploy.html), espaloma parameter could be applyed without such error log.

Could you please check if there are any issues with incorporating Espaloma into the OpenFF system build?

Best regards.

openfe quickrun network_setup/transformations/easy_rbfe_lig_ejm_31_complex_lig_ejm_42_complex.json -o results/easy_rbfe_lig_ejm_31_complex_lig_ejm_42_complex.json -d results/easy_rbfe_lig_ejm_31_complex_lig_ejm_42_complex
Loading file...
Planning simulations for this edge...
Starting the simulations for this edge...
'utf-8' codec can't decode byte 0x80 in position 64: invalid start byte
'utf-8' codec can't decode byte 0x80 in position 64: invalid start byte
'utf-8' codec can't decode byte 0x80 in position 64: invalid start byte
'utf-8' codec can't decode byte 0x80 in position 64: invalid start byte
'utf-8' codec can't decode byte 0x80 in position 64: invalid start byte
'utf-8' codec can't decode byte 0x80 in position 64: invalid start byte
Done with all simulations! Analyzing the results....
Here is the result:
        dG = None ± None

[Kamecs]

Incidentally, I'm using the Espaloma 0.4 libraries that I installed with conda.

[IAlibay]

Hi @Kamecs,

Could you provide more information on what you code you tried? It's not immediately clear to us how to reproduce this issue.

[Kamecs]

Dear @IAlibay

Thank you for your response.
I apologize for the poor quality of the information I provided.

In addition, the error that prevented execution was my own mistake. On the other hand, logs like
=== utf-8' codec can't decode byte 0x80 in position 64: invalid start byte ===
are still present. I would like you to confirm this point.

In an environment where openfe 1.11.1 and espaloma 0.4.0 are installed, creating an RBFE edge object and running FEP as described below will produce the above log:

Best regards.

import sys, os, re
import openfe
from openfe import SmallMoleculeComponent

from openfe.setup import LomapAtomMapper
from kartograf import KartografAtomMapper

from openfe.setup.ligand_network_planning import generate_radial_network
from openfe.setup.ligand_network_planning import generate_minimal_spanning_network
from openfe.setup.ligand_network_planning import generate_lomap_network
from openfe.utils.atommapping_network_plotting import plot_atommapping_network

from rdkit import Chem
supp = Chem.SDMolSupplier("tyk2_ligands.sdf", removeHs = False)
ligands = [SmallMoleculeComponent(sdf) for sdf in supp]

# Create a radial, choosing the first ligand as central ligand
ligand_network = generate_radial_network(
    ligands = ligands[1:],
    central_ligand = ligands[0],
    mappers = KartografAtomMapper()
)

from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol
from openff.units import unit

def create_setting():
    rbfe_settings = RelativeHybridTopologyProtocol.default_settings()

    # change the value
    rbfe_settings.thermo_settings.temperature = 300.0 * unit.kelvin
    rbfe_settings.protocol_repeats = 3
    rbfe_settings.integrator_settings.remove_com = True
    rbfe_settings.simulation_settings.minimization_steps  = 10000
    rbfe_settings.simulation_settings.equilibration_length = 1.0 * unit.nanosecond
    rbfe_settings.simulation_settings.production_length = 2.0 * unit.nanosecond
    rbfe_settings.simulation_settings.sampler_method = "repex"
    rbfe_settings.simulation_settings.real_time_analysis_interval = 0.01 * unit.nanosecond
    rbfe_settings.simulation_settings.real_time_analysis_minimum_time = 0.0 * unit.nanosecond 
    rbfe_settings.solvation_settings.box_shape = "octahedron"
    rbfe_settings.forcefield_settings.small_molecule_forcefield = "espaloma-latest.pt"
    rbfe_settings.output_settings.output_indices = "not water"
    rbfe_settings.output_settings.checkpoint_interval = 1.0 * unit.nanosecond
    return rbfe_settings


solvent = openfe.SolventComponent(
    positive_ion = "Na+",
    negative_ion = "Cl-",
    ion_concentration = 0.15 * unit.molar
)
protein = openfe.ProteinComponent.from_pdb_file("tyk2_protein.pdb")


transformations = []
for mapping in ligand_network.edges:
    for leg in ['solvent', 'complex']:
        # use the solvent and protein created above
        sysA_dict = {
            'ligand': mapping.componentA,
            'solvent': solvent}
        sysB_dict = {
            'ligand': mapping.componentB,
            'solvent': solvent
        }
        rbfe_settings = create_setting()

        rbfe_settings.solvation_settings.solvent_padding = 2.0 * unit.nanometer
        if leg == 'complex':
            sysA_dict['protein'] = protein
            sysB_dict['protein'] = protein
            rbfe_settings.solvation_settings.solvent_padding = 1.2 * unit.nanometer

        rbfe_protocol = RelativeHybridTopologyProtocol(
            settings = rbfe_settings
        )

        # we don't have to name objects, but it can make things (like filenames) more convenient
        sysA = openfe.ChemicalSystem(sysA_dict, name=f"{mapping.componentA.name}_{leg}")
        sysB = openfe.ChemicalSystem(sysB_dict, name=f"{mapping.componentB.name}_{leg}")

        prefix = "easy_rbfe_"  # prefix is only to exactly reproduce CLI

        transformation = openfe.Transformation(
            stateA = sysA,
            stateB = sysB,
            mapping = {'ligand': mapping},
            protocol = rbfe_protocol,  # use protocol created above
            name = f"{prefix}{sysA.name}_{sysB.name}"
        )
        transformations.append(transformation)

network = openfe.AlchemicalNetwork(transformations)


import pathlib
# first we create the directory
transformation_dir = pathlib.Path("network_setup/transformations")
transformation_dir.mkdir(exist_ok = True)

# then we write out each transformation
for transformation in network.edges:
    transformation.dump(transformation_dir / f"{transformation.name}.json")