Re: Question on restraint behavior when part of system decouples

From: Brian Radak (brian.radak_at_gmail.com)
Date: Mon Jan 29 2018 - 11:19:00 CST

Jerome has already given a great response, but I think I can add a little
bit regarding alchemical bonds and the use of extraBonds. These features
were added recently (by me) and are probably not well known, nonetheless,
they were not really intended for ligand binding, so you probably won't
need to change any of them.

As of NAMD 2.11 bonds in the PSF (or added via extraBonds) can be affected
by alchemical coupling if desired, but not be default. You can set
alchBondLambdaEnd in much the same way as alchVdwLambdaEnd to control the
region on which this scaling occurs and how quickly. There is also the
alchBondDecouple keyword, which is somewhat similar to alchDecouple.

There are two modes:

alchBondDecouple = False (default) - Only bonds involving atoms *outside
the alchemical region* are scaled, bonds between atoms in the alchemical
region are not scaled.

*If you are just decoupling a ligand, then this mode does nothing
different.*
alchBondDecouple = True - All bonds involving alchemical atoms are
decoupled. *This means that the alchemical region should blow apart into
non-interacting ideal gas atoms*. You probably don't want this.

There is also generally a lot of confusion surrounding the alchDecouple
keyword, which decouples non-bonded interactions between alchemical atoms
when off (default), but keeps them when true. With the default bonded
settings, alchDecouple = True means that the non-interacting ligand
endpoint behaves just like the molecule in the gas phase, but with a
periodic boundary if using PME. The use of this keyword can totally change
what thermodynamic cycle is appropriate and what assumptions are being made
regarding the unbound ligand.

HTH,
Brian

P.S. WARNING - I just discovered a subtle bug when using LJcorrection and
alchVdwLambdaEnd = 0.0 - DO NOT USE THIS COMBINATION WITH CONSTANT
PRESSURE. I'll hopefully have a patch up soon.

On Mon, Jan 29, 2018 at 10:53 AM, Jérôme Hénin <jerome.henin_at_ibpc.fr> wrote:

> Hi Randy,
>
> On 29 January 2018 at 16:19, Randy J. Zauhar <r.zauhar_at_usciences.edu>
> wrote:
>
>> Hi, I have always been troubled by the following aspect of the way
>> constraints are applied in the typical ‘alchemical’ transformation when
>> evaluating free energy of binding of ligand to protein.
>>
>> 1) If ligand is in complex with protein with restraints, and some
>> restraints involve both ligand and protein atoms (say dihedral angle
>> involving protein and ligand) and others only ligand atoms (say ligand
>> RMSD), and the ligand is decoupled from the complex - do all restraint
>> energy contributions ‘disappear’ during decoupling, or only the
>> protein-ligand part, or do NONE of the contributions disappear? In the end,
>> is the ligand effectively unrestrained in vacuum, or are the restraints
>> still in place?
>>
>
> If the restraints are implemented within the Colvars module, then they are
> unaffected by the alchemical lambda and will persist in the decoupled
> system. If they are implemented as extraBonds, I am not entirely sure,
> depending on the options that allow for perturbing bonded terms.
>
>
>> 2) Similarly, when the ligand is decoupled from bulk solvent, why does
>> the recommended protocol only restraint the conformation, why not position
>> and and orientation as well, to mimic the loss of positional and rotational
>> freedom that occurs upon complex formation?
>
>
> I'm not sure which recommended protocol you have in mind, but
> translational and rotational restraints when decoupling from isotropic bulk
> solution would have zero effect: they do not change the excess free energy
> that is being calculated. Restraint in the bulk calculation will not mimic
> the entropic effects of binding. More typically, restraints in the
> "protein" leg of the alchemical transformation are there to improve
> numerical convergence by increasing the overlap between the end-states. And
> then the free energy of those restraints needs to be calculated to account
> for those rotational and translational degrees of freedom.
>
> And again, if ligand in bulk decouples from solvent, is the end result a
>> small molecule in vacuum, or a small molecule in vacuum with restraints?
>> Those are very different situations!
>>
>
> They are, but the effect of those restraints in vacuum is exactly the same
> as in bulk solution, so their impact on the deltaG of decoupling is zero.
>
> I had assumed that the decoupling ONLY involved the nonbonded interactions
>> between parts of the system,
>
>
> Typically, yes.
>
>
>> and that you needed to handle everything involving restraints via
>> integration;
>
>
> Or analytically if possible.
>
>
>> but if that is so, my mental picture of how the thermodynamic cycle works
>> is really off.
>>
>
> There is more than one way to design a thermodynamic cycle for double
> decoupling, but I think the statements I've made above are general.
>
> Jerome
>
>
>

This archive was generated by hypermail 2.1.6 : Wed Sep 18 2019 - 23:19:36 CDT