Ok, many thanks.
So you only have dV/dl contributions from restraints and the only user restraints you have defined are the ones from your DISANG file. That's obviously what you want. Assuming that the code is ok your results should be fine. You could try to check if applying the restraints to the other copy or modifying the force constants of the restraints make any difference but I guess that won't give you much more information.
The timasks are redundant in your case (probably you could just specify there whatever pmemd would accept to let you run the simulation) but needed to allow TI simulations.
I suspect I can't add much more to this discussion.
Cheers,
Hannes.
________________________________________
From: vladimir.palivec.marge.uochb.cas.cz [vladimir.palivec.marge.uochb.cas.cz]
Sent: 27 March 2015 15:46
To: AMBER Mailing List
Subject: Re: [AMBER] Free energy calculation: Thermodynamic integration with use of restraints
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
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Fri Mar 27 2015 - 10:00:03 PDT
