Re: [AMBER] Fw: [ Calculation of ligand restraints in TI]

From: Josiah Bones <>
Date: Wed, 30 Jun 2021 05:36:44 +0000

Hi Taisung,

Thanks again for this information.

I have since run some TI calculations attempting absolute binding free energies (with separate runs for water+ligand as is more common).

I have utilised default smoothstep settings, and have implemented a set of simulations where the restraints are turned on (no softcore, V0 is ligand in complex with no restraints, V1 is ligand in complex +restraints).

However, I have realised that it is possible my simulation is not correctly dealing with the main restrained leg of the run, as you mentioned "TypeRestBA", and this is default in my simulations. I believe the correct behaviour in my case would be to have the restraints maintained through the entire main leg of decoupling the restrained ligand from the protein - Even when decoupled, these restraints should be fully active I think? Am I right in thinking I have made a mistake by keeping the default behaviour where the restraints are gradually turned off during decoupling? However, I expect the default behaviour would be desirable for the phase where restraints are turned on?

Any comments about the impact of this variable on how to correctly implement Boresch restraints would be much appreciated. Thanks!


From: <>
Sent: Tuesday, 22 June 2021 3:46 AM
To: 'AMBER Mailing List' <>; Josiah Bones <>
Subject: RE: [AMBER] Fw: [ Calculation of ligand restraints in TI]

Hi Josiah,

        Here is what I know/believe:

        The 6DOF restraints on the bound ligand at its dummy state (i.e.,
lambda=1; not interacting w/ protein except the restraints), if correctly
set up, are similar to restraining the COM of the ligand in bulk solution.
Your ligand has 6DOF to move without changing anything in its dynamics or
free energy state. However, at the real state (i.e., the ligand is
interacting w/ protein in addition to the restraints), the restraints will
have effects on its dynamics and free energy state and hence you need to do
some corrections.

        When turning on gti_output=1, you will get the restraint energy for
each window. However, it is still not the correction term. Remember, what
you miss is the free energy of restraining ligand in the protein
environment. One quick thing you could do is: check the printed restraint
energy of the real state (lambda=0). If the standard deviation of the
restraint energy distribution is small, e.g., less than 2kcal/mol, you can
directly calculate the restraint free energy through Zwanzig equation. In
most cases, you should be able to get a good estimate.

        A classic reference to start with is " Absolute Binding Free
Energies: A Quantitative Approach for Their Calculation" 2003 J Phys Chem B
107(35). If I find good recent references, I will let you know.



-----Original Message-----
From: Josiah Bones <>
Sent: Thursday, June 17, 2021 6:09 PM
Subject: [AMBER] Fw: [ Calculation of ligand
restraints in TI]

Hi Taisung, thanks for your reply.

I'll give a very quick overview of what I'm trying to do, to try and get an
idea if I am on the right track.

I'm attempting a single simulation where V0=ligand in bulk solution, and
V1=ligand bound to protein. This is due to charges on the ligands. I intend
to use the 6DOF restraints on the bound ligand, while only restraining the
COM of the ligand in bulk solution to ensure it doesn't interact with the
nearby protein, but not limit its internal degrees of freedom. I believe the
COM restraints don't need to be accounted for in any special way in the
thermodynamic cycle?

In this particular case, I believe you are saying I would still need an
additional simulation to calculate the dG of turning on the Boresch
restraints for the bound state, and then use the analytical equation for
these Boresch restraints. I refer to the recent paper, regarding these

"Additionally, Boresch restraints have been implemented, which can be used
in an automated fashion for ABFE simulations with many diverse ligands."


"The implementation in AMBER20 also permits these restraints to be included
in the overall alchemical transformation such that the component of the free
energy arising from the restraints in the bound state can be computed in the
same way as other force field terms."

I saw the option gti_output in the Amber20 manual, and wondered if this
provided the contribution of the restraints within the single calculation.
If this is not the case, would you have any recommendations for what
settings I need to include in the different steps in order to capture the
thermodynamic cycle correctly, and to automate the use of the Boresch
restraints as detailed in the recent paper? I am comfortable with the
overall running of the primary simulation (capturing dG of binding in one
step) but not with the additional method required for capturing the Boresch
restraints. Any help would be greatly appreciated.

Thanks again,


From: <>
Sent: Thursday, 17 June 2021 10:22 PM
To: <>; Josiah Bones
Cc: 'David A Case' <>;
<>; 'Darrin York' <>
Subject: RE: [ [AMBER] Calculation of ligand
restraints in TI]

Hi Josiah,

        In the current AMBER20 TI implementation, restraints involving atoms
in the scmasks will be scaled with respect to lambda, i.e., if an atom
disappears at lambda=1 state, all restraints involving this atom will become
zero at lambda=1 state. This implementation was designed for relative
binding free energy--so that Ligand A and Ligand B can keep independent sets
of restraints when doing A-->B.

        This behavior, however, does not change what you need to do for
absolute binding free energy. I believe that a proper 6-DOF-restraint will
not have any effect on the free energy of a non-interacting ligand. You
still need to calculate the free energy of turning on the restraints of the
ligand to the protein

        You also can recover the pre-AMBER20 behavior through changing the
restraint lambda-dependence by enabling lambda-scheduling of "TypeRestBA"
(AMBER manual Section 23.1.7). Let me know if you want to do so.



Date: Thu, 17 Jun 2021 01:13:46 +0000
From: Josiah Bones <>
To: "" <>
Subject: [AMBER] Calculation of ligand restraints in TI

Hi there!

I am just reaching out for some help as to the behaviour of AMBER20
pmemd.cuda in how it deals with the 6DoF restraints commonly employed in
absolute binding free energy calculations, where the ligand's movement is
restrained in a binding pocket. These calculations typically utilise two
main simulations (one to decouple ligand from protein, and another for
decoupling the ligand in water). From my understanding, there is then the
requirement of running an additional simulation to capture the free energy
of turning on the restraints of the ligand to the protein, and finally an
analytical calculation to remove them.

>From reading recent works examining the changes in AMBER20 to the TI
capabilities, I am suspecting, but not certain, that the restraints are
already calculated in the ligand+complex simulation, removing the need for
an additional simulation. Is this assumption correct, and is anything needed
to be done to ensure the correct behaviour? And does this contract the
requirements of such an ABFE calculation in a single smoothstep process
using only two simulations: complex+ligand(restrained),
water+ligand(unrestrained), and then the analytical calculation to account
for the restraints?

Thank you kindly for any assistance,

AMBER mailing list

----- End forwarded message -----

AMBER mailing list

AMBER mailing list
Received on Tue Jun 29 2021 - 23:00:02 PDT
Custom Search