How to input yank to run hydration free energy by using openmm serialized system (.xml) files

[nividic]

@jchodera I was trying to run hydration free energy calculation of phenol. I would like to input Yank by using openmm files of the serialized systems (vacuum and solvated phases). I made a .yaml script file to run but I'm having troubles. I would appreciate if you could write an example of this specific case please, I have not found a lot in the documentation. I have attached my files as well. Thanks a lot

hydration.zip

[jchodera]

@andrrizzi prepared input files for SMIRNOFF hydration free energies this way. Perhaps he can add one of the example input files to the yank-examples repo when he returns from vacation? @davidlmobley is also working to set up a repo to put all his example calculations for the SMIRNOFF paper into---he should already have the input scripts for that.

[nividic]

@jchodera. @davidlmobley provided me with a script example I was able to start running my testing. Thanks a lot

[jchodera]

For reference, that repo is https://github.com/mobleylab/SMIRNOFF_paper_code

[leelasd]

@andrrizzi Can you please provide a sample yaml script for running hydration free energy and binding free energy with openmm serialized systems ? I tried following the manual. But seems like I am missing something.

[andrrizzi]

Hi @leelasd. Here is an example:

---
options:
  minimize: yes
  verbose: yes
  temperature: 298.15*kelvin
  pressure: 1*atmosphere
  number_of_iterations: 5000
  anisotropic_dispersion_cutoff: auto

systems:
  hydration-system:
    phase1_path: [input/solvent1.xml, input/solvent1.pdb]
    phase2_path: [input/solvent2.xml, input/solvent2.pdb]

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]
        lambda_sterics:        [1.00, 1.00, 1.00, 1.00, 1.00]

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

Hope this helps!

[leelasd]

Thanks @andrrizzi !! But how to define which atoms to be annihilated? Does it figure it out automatically?

[andrrizzi]

Yes, with that script it modifies anything that is not water or ions. If your 2 phases have uncommon solvents (or if your water molecules have an uncommon residue name), you can specify what is the solvent like this:

systems:
  hydration-system:
    phase1_path: [input/solvent1.xml, input/solvent1.pdb]
    phase2_path: [input/solvent2.xml, input/solvent2.pdb]
    solvent_dsl: resname WAT and resname Ca+

[leelasd]

Thanks for the quick response, much appreciated. I had to add number_of_iterations: 1 to your script to get it to run.

Can you also please post a script with GROMACS input as an example ?

[andrrizzi]

I had to add number_of_iterations: 1 to your script to get it to run.

Thanks, I've updated the example above. With Gromacs files it's pretty much the same except you need to specify the long-range treatment (which is already included in the OpenMM xml file).

solvents:
  PME:
    nonbonded_method: PME
    nonbonded_cutoff: 11*angstroms
    switch_distance: 10*angstroms
    ewald_error_tolerance: 1.0e-4

systems:
  hydration-system:
    phase1_path: [input/solvent1.top, input/solvent1.gro]
    phase2_path: [input/solvent2.top, input/solvent2.gro]
    solvent: PME
    gromacs_include_dir: OPTIONAL/path/to/extra/parameters/

[leelasd]

Thanks, @andrrizzi.

Just one more question. Can Yank support OPLS-AA force field style LJ combination rules (style: 3 for in gromacs input file ) ?

[jchodera]

YANK supports whatever OpenMM's [GromacsTopFile]() supports, but it looks like this currently only supports combination rule types 1 and 2. You can create an issue on the OpenMM project page requesting support for type 3---I don't think it would be hard to add!

[andrrizzi]

Also, interactions involving alchemical atoms currently support only a single type of combination rule (see this alchemy.py snippet). We plan to add more, but unfortunately it's not on the high priority list right now. All the interactions involving only solvent atoms (and receptor atoms for binding free energies) support everything that OpenMM supports.

[jchodera]

Oh, good point. This will have to wait for general support for metaprogramming CustomNonbondedForce alchemical modifications: https://github.com/choderalab/openmmtools/issues/172