Doesn't work means that it says that I have to specify at least one atom
in the timask. I have printed just averaged values. Here is output for
lambda = 0.5
| TI region 1
A V E R A G E S O V E R 200 S T E P S
NSTEP = 200000 TIME(PS) = 400.000 TEMP(K) = 300.09 PRESS =
-24.6
Etot = -88964.0251 EKtot = 19723.7890 EPtot =
-108687.8141
BOND = 897.2671 ANGLE = 2388.0587 DIHED =
3215.9752
1-4 NB = 1172.9601 1-4 EEL = 14646.7577 VDWAALS =
15882.7941
EELEC = -146892.0553 EHBOND = 0.0000 RESTRAINT =
0.4282
EAMBER (non-restraint) = -108688.2423
EKCMT = 8074.5637 VIRIAL = 8271.1655 VOLUME =
344621.5907
Density =
0.9622
DV/DL = 0.8564
Ewald error estimate: 0.5405E-04
------------------------------------------------------------------------------
NMR restraints: Bond = 2.066 Angle = 0.002 Torsion = 0.005
===============================================================================
| TI region 2
A V E R A G E S O V E R 200 S T E P S
NSTEP = 200000 TIME(PS) = 400.000 TEMP(K) = 300.09 PRESS =
-24.6
Etot = -88964.0251 EKtot = 19723.7890 EPtot =
-108687.8141
BOND = 897.2671 ANGLE = 2388.0587 DIHED =
3215.9752
1-4 NB = 1172.9601 1-4 EEL = 14646.7577 VDWAALS =
15882.7941
EELEC = -146892.0553 EHBOND = 0.0000 RESTRAINT =
0.4282
EAMBER (non-restraint) = -108688.2423
EKCMT = 8074.5637 VIRIAL = 8271.1655 VOLUME =
344621.5907
Density =
0.9622
DV/DL = 0.8564
Ewald error estimate: 0.5405E-04
------------------------------------------------------------------------------
NMR restraints: Bond = 2.066 Angle = 0.002 Torsion = 0.005
===============================================================================
DV/DL, AVERAGES OVER 200 STEPS
NSTEP = 200000 TIME(PS) = 400.000 TEMP(K) = 1.61 PRESS =
0.0
Etot = 0.0000 EKtot = 0.0000 EPtot =
0.8564
BOND = 0.0000 ANGLE = 0.0000 DIHED =
0.0000
1-4 NB = 0.0000 1-4 EEL = 0.0000 VDWAALS =
-0.0000
EELEC = 0.0000 EHBOND = 0.0000 RESTRAINT =
0.8564
EAMBER (non-restraint) = 0.0000
EKCMT = 0.0000 VIRIAL = 0.0000 VOLUME =
0.0000
Density =
0.0000
DV/DL = 0.8564
Ewald error estimate: 0.0000E+00
------------------------------------------------------------------------------
NMR restraints on final step:
------------------------------------------------------------------------------
Final Restraint Analysis for coords: 3.restrt
Restraints, deviations, and energy contributions: pencut = 0.10
------------------------------------------------------------------------------
First atom Last atom curr. value target deviation penalty
------------------------------------------------------------------------------
CA GLU 34 -- C3 DPN 9320: 15.984 16.630 0.646 4.170 d
0: 0
Total distance penalty: 4.170
Total angle penalty: 0.004
Total torsion penalty: 0.010
------------------------------------------------------------------------------
----------------------------------------------------------------------------
################################
In between are rms fluctuations which I didn't copy. Input file is like this:
icfe=1, clambda = 0.5,
ifsc=0,
timask1=':9321', scmask1='',
timask2=':9320', scmask2='',
nmropt = 1,
where 9320, 9321 is exactly the same molecule. DISANG file contain 3
dihedral, 2 angle, and 1 dist restraint. I used harmonic potential with a
constant of value 10 for both angles and distances.
> "doesn't work" means exactly what?
>
> Your output file should give you a summary under the heading "DV/DL,
> AVERAGES OVER"
> I would expect all component terms to be zero except for RESTRAINT. Can
> you paste this section?
>
> I would have set the masks something like timask1=':1' and timask2=':2'
> assuming that the two ligand copies are residue 1 and 2 but I suppose
> that's what you mean, also. implying that the restraints are only applied
> to one copy.
>
> Can you provide the complete input file?
>
> ________________________________________
> From: vladimir.palivec.marge.uochb.cas.cz
> [vladimir.palivec.marge.uochb.cas.cz]
> Sent: 26 March 2015 15:22
> To: AMBER Mailing List
> Subject: Re: [AMBER] Free energy calculation: Thermodynamic integration
> with use of restraints for
>
> Hi,
>
> suggestion 1 doesn't work, however, suggestion 2. looks promising.
>
> Here is the setup:
> timask1='not_restrained', scmask1='',
> timask2='restrained', scmask2=' ,
>
> for lamda = 1
> I got :
> RESTRAINT = 0.5716
> DV/DL = 0.5716
>
> for lambda = 0.5
> RESTRAINT = 0.5010
> DV/DL = 1.0019
> so it it twice the restraint energy
>
> or lambda = 0.0
> RESTRAINT = 0.0000
> DV/DL = 5.0038
>
> with these three points, this results in dG_restrains_on = +1.9 kcal/mol,
> which seems reasonable. However, I am still not sure, whether this is the
> right way to do it.
>
> Vlada
>
>> Hi again,
>>
>> I have not done this by myself so I can only make suggestions how I
>> _think_ this may work. So you will have to experiment unless you get an
>> answer form someone more knowledgeable.
>>
>> Set nmropt=1 and define the necessary &wt namelist to define your
>> restraints. I assume you can define multiple restraints simultaneously.
>> If not you would have to do it one-by-one... At the end you are
>> interested in scaling only those restraints with lambda and nothing
>> else.
>>
>> Suggestion 1:
>> Create a prmtop file with only one ligand. I assume that's exactly the
>> prmtop you already have for the annihilation step (mutate ligand to
>> nothing). Set icfe=1, ifsc=0 and the timasks and scmasks to empty
>> strings. The hope here is that pmemd would only lambda scale the NMR
>> restraints. I do not know if using empty strings would switch off the
>> TI
>> code or this doesn't work for some other reason. If so
>>
>> Suggestion 2:
>> Create a prmtop file with two identical copies of the ligand (that's how
>> you would setup relative free energy simulations). Set icfe=1, ifsc=0
>> and
>> timask1 to the first ligand and timask2 to the second ligand. The
>> scmasks
>> stay empty (which also implies ifsc=0). In this way ligand 1 would be
>> transformed into ligand 2 but since all force field parameters are
>> identical the resulting free energy should be zero. Not very sufficient
>> obviously so I hope suggestion 1 just works out for you.
>>
>> Please report back if any of this works.
>>
>>
>> Thanks.,
>> Hannes.
>>
>>
>> ________________________________________
>> From: vladimir.palivec.marge.uochb.cas.cz
>> [vladimir.palivec.marge.uochb.cas.cz]
>> Sent: 26 March 2015 08:34
>> To: AMBER Mailing List
>> Subject: Re: [AMBER] Free energy calculation: Thermodynamic integration
>> with use of restraints
>>
>> Hello,
>>
>> thank you for an advice. This might be relly the way, however, I do not
>> know how to perform this step. I dont know how to couple a system
>> between
>> free and restrained state (how exatly to do it, as in the manual there
>> is
>> not a word about this).
>>
>> Please, can you give me any suggestions?
>>
>>
>>> I just realise that in both sander and pmemd, NMR restraints (nmropt=1)
>>> are scaled too as per the link in your original email. This may be the
>>> better option to pursue.
>>>
>>> ________________________________________
>>> From: Loeffler, Hannes (STFC,DL,SC)
>>> Sent: 25 March 2015 19:12
>>> To: AMBER Mailing List
>>> Subject: RE: [AMBER] Free energy calculation: Thermodynamic integration
>>> with use of restraints
>>>
>>> So, what you want to do is an "absolute" free energy calculation with
>>> the
>>> TI code. You also want to only couple a set of restraints to lambda.
>>> If
>>> I understand the Boresch' method correctly you would need to define six
>>> such restraints: one bond, two angles, three dihedrals. Right?
>>>
>>> Sorry, I'm not too well at the moment so I may get this wrong but
>>> wouldn't
>>> it be possible to define those restraints in terms of existing atoms
>>> and
>>> introduce additional bonds, angles and dihedrals between them. In the
>>> input file you would then define the TI region by including those atoms
>>> only while the softcore mask is empty. Maybe that's a plan but I would
>>> have to rethink it. The only other idea I can come up with is to
>>> create
>>> multiple topology files with those additional terms pre-scaled by the
>>> current lambda but then you would need to do a EXP/BAR estimation for
>>> the
>>> associated free energy as I can't see at the moment how this could be
>>> fit
>>> into the TI formalism.
>>>
>>> Cheers,
>>> Hannes.
>>>
>>> ________________________________________
>>> From: vladimir.palivec.marge.uochb.cas.cz
>>> [vladimir.palivec.marge.uochb.cas.cz]
>>> Sent: 25 March 2015 17:55
>>> To: amber.ambermd.org
>>> Subject: [AMBER] Free energy calculation: Thermodynamic integration
>>> with
>>> use of restraints
>>>
>>> Hello,
>>>
>>> this question is similar to this one:
>>>
>>> http://archive.ambermd.org/201205/0093.html
>>>
>>> however, I think I quite did not get it.
>>>
>>> I would like to do a thermodynamic integration simulations to calculate
>>> dG
>>> of binding of a ligand to a protein. After search through a literature
>>> I
>>> decided to do it in multiple steps:
>>>
>>> 1. constrain a ligand in a certain position using 6 restraints
>>> - here i need to calculate dG of proposing restraints
>>>
>>> 2. turn off el. i.
>>>
>>> 3. using soft core, turn off vww interactions
>>>
>>> 4. analytically remove restraints
>>>
>>>
>>> I am not sure how to perform the step 1. I thought it is possible to
>>> use
>>> pmemd with single topology and gradually propose restraints and thus
>>> get
>>> the dG (with use of nmropt=1 and DIANG file). However, I have no found
>>> any
>>> way to do it - either in manual or forum. Please, can somebody give me
>>> a
>>> hint how this calculation should be performed? How to set up the
>>> calculation and how to obtain energy of proposing restraints?
>>>
>>> Is it even possible to do it with one topology file?
>>>
>>> Thank you very much for any help!
>>>
>>> Vlada
>>>
>>>
>>>
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> AMBER mailing list
>>> AMBER.ambermd.org
>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>
>>> _______________________________________________
>>> AMBER mailing list
>>> AMBER.ambermd.org
>>> http://lists.ambermd.org/mailman/listinfo/amber
>>>
>>
>>
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>> _______________________________________________
>> AMBER mailing list
>> AMBER.ambermd.org
>> http://lists.ambermd.org/mailman/listinfo/amber
>>
>
>
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
> _______________________________________________
> AMBER mailing list
> AMBER.ambermd.org
> http://lists.ambermd.org/mailman/listinfo/amber
>
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Mar 27 2015 - 09:00:05 PDT
