Hi Thomas,
thank you very much for your prompt answer !
Although your reply confirmed some of my hypothesis it also
generated couple of additional questions but maybe I would like first
understand fully the "processing" of the force-constants which are
supplied in RST file (in &rst records) for the given TI MD run using the
given value of lambda.
During my first "explorational" trial with "T4 lysozyme + benzene" I used
just one RST file common for all
lambdas (0.1,0.2,...,0.9 ) with all force constants = 20 kcal/(mol A2) or
kcal/(mol rad2) and have found
in output files informations that they were rescaled:
#lambda 0.1
Scaling up weight of 20.000 by 1.111
New weight 22.222
#lambda 0.2
Scaling up weight of 20.000 by 1.250
New weight 25.000
#lambda 0.3
Scaling up weight of 20.000 by 1.429
New weight 28.571
#lambda 0.4
Scaling up weight of 20.000 by 1.667
New weight 33.333
#lambda 0.5
Scaling up weight of 20.000 by 2.000
New weight 40.000
#lambda 0.6
Scaling up weight of 20.000 by 2.500
New weight 50.000
#lambda 0.7
Scaling up weight of 20.000 by 3.333
New weight 66.667
#lambda 0.8
Scaling up weight of 20.000 by 5.000
New weight 100.000
#lambda 0.9
Scaling up weight of 20.000 by 10.000
New weight 200.000
meanwhile I "decoded" how these scaling constants are connected with
lambdas values. They are simply 1/(1-lambda) in other words each reported
scaled value of the given force constant in output file noted as the "New
weight ..." is simply the relevant value from the given RST file (&rst
record)
divided by the (1-lambda) for example 22.222 = 20/(1-0.1), 25.000 =
20/(1-0.2) etc. Here an interesting thing is that in Amber11 manual is
written (page 101, comment of dvdl_norest):
-------------------------------------------------------------------------------------------------------------
Please note that the force constants of restraints are divided by lambda
if soft core potentials are switched
on. This results in the restraint being applied in full to the
disappearing atoms at any lambda.
-------------------------------------------------------------------------------------------------------------
So here are two information:
a)
"restraints are divided by lambda" which is not in agreement with my just
commented experience from which is clear that the restraints are divided
by (1-lambda).
b)
The scaling of the restraints (that reported in output files) is present
probably just in case of usage of softcore potentials not in the case when
the softcore potentials are not used
and just simply linear mixing is used instead. Am I right ?
Anyway if I understood well the purpose of the scaling of the restraint
force constant "K" by factor "X" is just because during the consequent MD
calculation it will be
multiplied by factor 1/X so the scaling ensures, that the final restraint
applied to the ligand for the given lambda will be exactly that which is
defined in RST file (&rst record).
So for example this notice in output file:
-------------------------------------------------------
Scaling up weight of 20.000 by 2.000
New weight 40.000
--------------------------------------------------------
is really just informative and the final restraint applied to the ligand
will be still original value of 20 in this case. Am I right ?
If this is true, the first step of the decoupling of the restrained ligand:
dA*
(RL)s -----> (R)s***(L)v
(i.e. gradually turn on restraints while el./vdw interactions are
simultaneously gradually turned off using softcore potentials)
of the whole process
dA* dAr
(RL)s -----> (R)s***(L)v----->(R)s + (L)v (1)
might be in principle realised in one Amber TI step by providing for each
lambda different RST file (let say RST_lambda) where I use
different values of force constants in each RST_lambda (i.e. smoothly
growing with increasing lambda in my particular case).
Of course immediately arises question how exactly smoothly increase
restraints values in RST_lambda files to obtain converged results (when
simultaneously decoupling both native
interactions using SC potentials.). In ideal case this should somehow
"copy" the speed of vanishing of the native (vdw/el) interactions.
Just in this context I understand your comment:
"with restraints, I think there is no true *one* step solution yet. You
would activate your restraints
first (no softcore required), then disappear the whole ligand (eel&vdw)
while restraints remain in place"
which perhaps means to divide/rewrite step:
dA*
(RL)s -----> (R)s***(L)v
like this:
dAa* dAb*
(RL)s ----->(R***L)s -----> (R)s***(L)v
where (R***L)s is the new artificial itermediate state: fully natively
(vdw/el) interacting R with L (complex) in solvent but with already
applied target restraints
which are gradually introduced during the "dAa*-step". This has the
advantage (comparing to one step approach) that the way of smooth
introduction of restraits
during "dAa*-step" should not affect much the convergence as the position
of the ligand is during this WHOLE step sufficienly "fixed" with native
interactions.
But as dA* = dAa* + dAb* and the dA* is the change in total free energy
while going from "(RL)s" to "(R)s***(L)v" I still think that also the
dv/dl contributions
from the restraint energy response should be included into the final dv/dl
for the given lambda which implies dvdl_norest=0 for both subprocesses
(dAa*, dAb*). It seems to me (if I understood everything well) that you
agree here with me just for the substep "dAa*" and in case of
"dAb*-step" (gradual vanishing of vdw/el interactions of the fully
restrained ligand) you suggest dvdl_norest=1.
If this is true, please explain why (what is your argumentation for
excluding restraint dv/dl contribution) ?
I am really sorry for such long contributions but if I should be able to
use Amber TI to obtain reliable results I need to understand well at least
the necessary things. BTW I hope that this discussion will be useful
(now/in future or for the eventual future TI tutorial) also for the others.
Thank you very much in advance for your (or anyones) comments !
Best wishes,
Marek
Dne Mon, 07 May 2012 09:49:31 +0200 <steinbrt.rci.rutgers.edu> napsal/-a:
> Hi,
>
>> I suppose that the restraints (distance/angle) are zero in state 1 and
>> are gradually turned on during the first step (1--->2) with
>> maximum/final
>> values in state 2
>> (these values (force constants) are then used in analyt. calc. of dAr)
>> as
>> the whole process (1--->2--->3) should connect two real states. I hope
>> this assumption is OK.
>
> that sounds correct.
>
>> So my question is how to do this first step with Amber assuming "One
>> step"
>> approach (using simultaneously softcore potentials for vdw and el.
>> interactions).
>
> with restraints, I think there is no true *one* step solution yet. You
> would activate your restraints first (no softcore required), then
> disappear the whole ligand (eel&vdw) while restraints remain in place (no
> vacuum step is then required). That's still two steps less than the old
> setup :-)
>
>> Restraints for the ligand are defined (in proper RST file) just in case
>> of
>> the "V0 system" as the ligand is not present in
>> the "V1 system". Am I right ?
>
> yes. That setup comes from Ambers dual topology approach in which V1
> doesnt even know that the ligand exists.
>
>> So dvdl_norest=0 or dvdl_norest=1 in this particular case ?
>
> yes, dvdl_norest should only be on when the restraints shall not be
> counted, e.g. while disappearing the ligand.
>
>> #lambda 0.1
>> Scaling up weight of 20.000 by 1.111
>
> This is done to prevent your restraint force from being scaled down when
> you want it to stay in place. The alternative would be to supply a scaled
> restraint in V0 at each lambda. Since the restraint only exists in V0, it
> will be multiplied by lambda later during the MD, this scaling just
> effectively prevents that from happening.
>
> The goal of this is that if you set a 5.0 kcal restraint, sander will
> apply that force at any lambda with no adjustment needed by the user. It
> is done in this odd way because a restraint is really like an *internal*
> potential of the ligand, but it cant be treated like one in the code.
>
> Thomas
>
> Dr. Thomas Steinbrecher
> formerly at the
> BioMaps Institute
> Rutgers University
> 610 Taylor Rd.
> Piscataway, NJ 08854
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
> __________ Informace od ESET NOD32 Antivirus, verze databaze 7116
> (20120507) __________
>
> Tuto zpravu proveril ESET NOD32 Antivirus.
>
> http://www.eset.cz
>
>
>
--
Tato zpráva byla vytvořena převratným poštovním klientem Opery:
http://www.opera.com/mail/
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Mon May 07 2012 - 08:30:03 PDT
