Using yank for calculation of mixing free energy

[Eli-VDB]

hello,

I would like to use the YANK program to calculate mixing free energies for a polymer-drug system. I figured that the mixing free energy could be calculated in similar fashion as the binding free energy for a host-guest system or solvation free energies.

To be more precise, I have 15 polymer chains + n drug molecules, with n being an integer ranging from 1 to 476. Here I have 2 options either I remove all drug molecule at once to calculate a binding free energy or in my case a mixing free energy, or I just remove one drug molecule (which I think is in some way equivalent to the cases studied by your group). The solvent should be vaccuum in both cases.

The systems are constructed in openmm (I am hence using pdb and xml files) but I am struggling to figure out the correct format of the yaml file for these type of calculations. Hence my question if you could provide me with some help on how to set-up a correct yaml file for my calculations.

At the moment the following yaml file is used:

options:
checkpoint_interval: 50
minimize: yes
verbose: yes
temperature: 298.15kelvin
pressure: 1atmosphere
anisotropic_dispersion_cutoff: auto

systems:
hydration-system:
phase1_path: [polymer_10drugmolecules.xml, polymer_10drugmolecules.pdb]
phase2_path: [drug_molecule_system.xml, drug_molecule.pdb]
solvent_dsl: resname ETO or resname TER or resname INI or (resname UNL and resSeq > 1)

protocols:
hydration-protocol:
solvent1:
alchemical_path:
lambda_electrostatics: [1.00, 0.75, 0.50, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
lambda_sterics: [1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.80, 0.70, 0.60, 0.50, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00]
solvent2:
alchemical_path:
lambda_electrostatics: [1.00, 0.75, 0.50, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
lambda_sterics: [1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.50, 0.40, 0.30, 0.20, 0.10, 0.00]

experiments:
system: hydration-system
protocol: hydration-protocol

As you can see from this file I selected the polymer residues (ETO TER INI) and all drug-molecules except for the 1st one as my solvent, in the other case (solvent 1) I have 1 single molecule which should be in vaccuum, Assuming I am doing everything correct here.

Could you provide me with some feedback about what is missing, what can be done differently or ...

Thank you in advance for the help! Feel free to ask for more information, hopefully everything is clear. Kind regards, Elias

I have also posted this at choderalab/yank-examples, but I believe this is the correct location to ask my question, sorry for the inconvenience

[andrrizzi]

Hi @Eli-VDB . If I understand correctly, the polymer will play the role of the receptor, not the solvent (which is vacuum) so you can leave solvent_dsl out of the YAML. You'll have to specify ligand_dsl, however, to specify the molecules you'll want to alchemically decouple.

If you want to alchemically decouple a single molecule, this should be pretty standard. But if you want to decouple 10 molecules at the same time you'll have problems with restraints (see the tutorial) as you'll likely need one restraint per molecule for good performance, and the YAML interface does not support this kind of setup. You might be able to make it work if you are willing to dive into the Python code, however.

[Eli-VDB]

Hi @andrrizzi, thank you for the reply.

I will attempt to use a complex-solvent system instead for a single drug-molecule first. I have changed the yaml file to the following:

systems:
  mixing-system:
    phase1_path: [polymer_10drugmolecules.xml, polymer_10drugmolecules.pdb]
    phase2_path: [drug_molecule_system.xml, drug_molecule.pdb]
    ligand_dsl: resname UNL and resSeq 1

protocols:
  mixing-protocol:
    complex: 
      alchemical_path:
        lambda_electrostatics: [1.00, 0.75, 0.50, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.80, 0.70, 0.60, 0.50, 0.40, 0.35, 0.30, 0.25, 0.20, 0.15, 0.10, 0.05, 0.00]
    solvent:
      alchemical_path:
        lambda_electrostatics: [1.00, 0.75, 0.50, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00, 0.95, 0.90, 0.85, 0.80, 0.75, 0.70, 0.65, 0.60, 0.50, 0.40, 0.30, 0.20, 0.10, 0.00]

experiments:
  system: mixing-system
  protocol: mixing-protocol

2 questions here

1. Am I still correct to specify the second phase as a single ligand molecule? and the first phase as a complex
2. I get the following error when running the yank simulation: '''python !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! CRITICAL: Experiment NaN ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! The following experiment threw a NaN! It should NOT be considered! ''' As system xml files and pdb files should normally be correct because I am able to run openmm simulations with them, I expect that something is wrong with my system definition, or am I perhaps messing things up with periodic/non-periodic boundary conditions? (i.e. my complex system is periodic, the ligand should be aperiodic but perhaps I have to specify this explicitly?)

Furthermore is not clear to me how Yank will treat this solvent in case the solvent is vacuum (in the protocol section). Furthermore,

thank you in advance! kind regards, Elias

[andrrizzi]

1. Yes, that sounds correct to me.
2. YANK should save in the experiment folder the state of the simulation right before it NaNned. It might help debugging what is going on.

If your OpenMM XML file is using vacuum, then YANK will treat the molecules in vacuum as well.