Peptide ligand is not alchemically modified

[juliebehr]

I am trying to run yank on a system where the ligand is a peptide, and it does not seem to recognize the peptide for alchemical modification, such that all of the "alchemical" states are exactly the same. The calculated energies for each replica in each of the different possible states does not vary.

Input files here: peptide-ligand-yank.zip

[jchodera]

Do you have an output dump you could send along too?

[andrrizzi]

Ah! I think I understand what happens. The current code relies on a DSL string to identify the ligand atoms, but YamlBuilder does not construct the string correctly. Does it works if you use the command line with the correct DSL string for the parameter --ligand?

[juliebehr]

(Output files are still zipping) I think the same thing happened with the command line interface but I might have been using it incorrectly?

[juliebehr]

Update: the files are too big to be logical to zip / upload anyway, but can be found on src in /home/behrj/experiments/

I'm retrying the command line now

[juliebehr]

My replicas are currently minimizing and presumably will be for a while (hours) before I'll be able to confirm whether the energies are still not changing, but for the record I'm running this:

mkdir output
yank prepare binding amber --setupdir=setup/systems/1sb0_peptide_RF/ --ligand='chainid 1' --store=output --iterations=2000 --nbmethod=CutoffPeriodic --cutoff=1*nanometer --temperature=310*kelvin --pressure=1*atmosphere --minimize -v
yank run --store=output --verbose

and have confirmed that 'chainid 1' selects the peptide:

import mdtraj as md
trj = md.load_pdb("setup/systems/1sb0_peptide_RF/complex.pdb")
ligand = trj.top.subset(trj.top.select('chainid 1'))
print([residue for residue in ligand.residues])
>> [LYS88, ARG89, ILE90, LYS91, GLU92, LEU93, GLU94, LEU95, LEU96, LEU97, MET98, SER99, THR100, GLU101, ASN102, GLU103, LEU104, LYS105]

[juliebehr]

Andrea has pointed out that inside yank it's actually loading from the prmtop, which condenses the system into 1 chain and therefore the ligand cannot be identified via 'chainid 1' -- not sure how to differentiate it from the receptor.

[juliebehr]

I think this is solved -- I'm able to run using --ligand='protein and resi >= 87' and I should be able to manually find what this needs to be modified to for my different setups. Maybe in the future it would be helpful if there was a way to tell tleap to keep the chains separate, so there could be a way to differentiate the ligand besides resname (since that doesn't work for peptides)?

[andrrizzi]

Thanks for figure things out!

Maybe in the future it would be helpful if there was a way to tell tleap to keep the chains separate

Yeah, that sounds like the best solution. I'll check out if that's possible on the ambertools manual.

[jchodera]

We'll need to figure out a way to have the ligand: section of the YAML file automagically identified when we construct these systems. Maybe we can just count atom indices and identify the ligand that way?

[andrrizzi]

For small molecules ligands the pipeline already automagically identifies it through the residue name (because tleap doesn't change it when it create the complex system). The main problem is the protein.

Maybe we can just count atom indices

This would work for the protein! We could even count residues instead.

[Lnaden]

I don't think there will be time to generalize the "ligand" selection process for 1.0, but this should be something we plan to fix for a future cut 1.X, X > 0

1.0 touch