choderalab / yank

An open, extensible Python framework for GPU-accelerated alchemical free energy calculations.
http://getyank.org
MIT License
181 stars 72 forks source link

Implementing RMSD or other restraints in solvent phase #1208

Open jeffcomer opened 4 years ago

jeffcomer commented 4 years ago

Currently, the yank documentation says: "The RMSD restraints must also be applied in the solvent phase of the calculation if they are on in the decoupled state, but this is not implemented yet."

An alternative method would be provided by restraining the ligand torsions as described in: https://github.com/choderalab/yank/pull/845

Either way, the restraint contribution to the free energy needs to be calculated during the solvent phase. I'm thinking about trying to implement this myself. Has this not been implemented because it would be extremely difficult or simply because the developers have other priorities?

I'm thinking to implement it by optionally allowing "lambda_restraints" during the solvent phase. Do the developers have any suggestions or advice before I attempt this?

andrrizzi commented 4 years ago

Hi @jeffcomer . We didn't implement this yet because we had other priorities. We were actually hoping to integrate this into a large overhaul of the restraints that would also allow us to run multiple restraints. The changes were almost finished in a branch, but it got stuck.

This being said, I don't expect this to be extremely difficult to implement, but it won't be trivial either. This feature would require changing the YAML syntax, the system/alchemical phase preparation code, and a few bits of the analysis.

jchodera commented 4 years ago

@jeffcomer : To clarify @andrrizzi's point, once we introduce the ability to have multiple restraints, the RMSD restraint can be introduced only at intermediate lambda values, with the harmonic, flat-bottom, or Boresch restraints used at the non-interacting endpoint. This would avoid the need to introduce RMSD restraints in the solvent phase.

We haven't ended up using the RMSD restraints much because we really wanted to instead restraint the RMSD of the ligand once we aligned the system on the receptor, which is a feature that hasn't been implemented in OpenMM yet. While computing the RMSD this way is straightforward, the gradient computation is much less straightforward. We think this would be much more useful than overall RMSD, however.