From: Brian Radak (bradak_at_anl.gov)
Date: Wed Jan 11 2017 - 10:20:21 CST
I think you need to be careful how you compute your derivative in this
case - just changing the radii means that you have a highly non-linear
coupling. What about epsilon?
On 01/11/2017 12:39 AM, Abhishek Kumar Singh wrote:
> I'm working with a protocol that is different from NAMD dual topology
> procedure. I'm transforming K+ to Na+ by different parameter (λ= 1.0,
> 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1 and 0.0), decreasing vdw
> radii. Then reading the dcd file with 2 parameters ( λ1=K+, λ0=Na+)
> for each window and have calculated energy difference. The area
> under the plot (derivate vs λ) gives me = -23.40 but in principle it
> should be = -18.60 kcal/mol.
>
>
> On Wed, Jan 11, 2017 at 1:58 AM, Radak, Brian K <bradak_at_anl.gov
> <mailto:bradak_at_anl.gov>> wrote:
>
> I ran this calculation very quickly and got something close to the
> reported value. The transformation is from SOD -> POT in a TIP3Pm
> box (which are the default parameters shipped with
> par_all36_prot.prm, for example). I added a 100 kcal/mol-A^2 bond
> with 0 length between the ions using extraBonds, vdwForceSwitching
> on, and NO softcore potential (i.e. alchVdwShiftCoeff 0.0).
>
> lambda delta f (in kT)
> ------------------------------
> 0.0 0.0000
> 0.1 7.3138
> 0.2 12.7031
> 0.4 19.8229
> 0.6 24.7088
> 0.8 28.3123
> 0.9 29.7617
> 1.0 31.0306 +/- 0.306
>
> 31.0306 --> 18.4 +/- 0.2 kcal/mol
>
> This was using WHAM. Using a "quick and dirty" TI protocol I got
> 31.3467 -->18.6 kcal/mol
>
> Brian Radak
>
> Postdoctoral Appointee
> Leadership Computing Facility
> Argonne National Laboratory
>
> 9700 South Cass Avenue, Bldg. 240
> Argonne, IL 60439-4854
> (630) 252-8643
> brian.radak_at_anl.gov <mailto:brian.radak_at_anl.gov>
>
> ________________________________________
> From: owner-namd-l_at_ks.uiuc.edu <mailto:owner-namd-l_at_ks.uiuc.edu>
> [owner-namd-l_at_ks.uiuc.edu <mailto:owner-namd-l_at_ks.uiuc.edu>] on
> behalf of Radak, Brian K [bradak_at_anl.gov <mailto:bradak_at_anl.gov>]
> Sent: Tuesday, January 10, 2017 12:27 PM
> To: namd-l_at_ks.uiuc.edu <mailto:namd-l_at_ks.uiuc.edu>; Hannes Loeffler
> Subject: RE: namd-l: Solvation free energy difference of
> alchemical transformation of K+ into Na+
>
> IMO, the "correct" way to do this kind of calculation is to have a
> harmonic bond between the two ions - this is weakly suggested in
> the manual by referring to an old and probably not widely read
> paper by Axelsen and Li. This does not affect the free energy that
> is observed, whereas Cartesian restrains on the atoms separately does.
>
> There might also be some discrepancy due to lack of a
> LJcorrection, which is off by default. The use of
> vdwForceSwitching might also cause a difference. Both of these are
> much more common options in CHARMM, where I suspect the reference
> calculations were performed.
>
> The WCA code is indeed present and largely undocumented. There
> have been some efforts to remedy this, but they are going slowly.
> Because of some slightly unfortunate code design, many of the
> recent improvements to the alchemy module were not extended to the
> WCA code - I don't foresee this changing because I'm the one that
> did the alchemy changes and I don't currently have time or
> resources to do the WCA changes.
>
> I will try to reproduce the calculation myself and see what I come
> up with.
>
> Cheers,
> Brian
>
> Brian Radak
> Postdoctoral Appointee
> Leadership Computing Facility
> Argonne National Laboratory
>
> 9700 South Cass Avenue, Bldg. 240
> Argonne, IL 60439-4854
> (630) 252-8643
> brian.radak_at_anl.gov <mailto:brian.radak_at_anl.gov>
>
> ________________________________________
> From: owner-namd-l_at_ks.uiuc.edu <mailto:owner-namd-l_at_ks.uiuc.edu>
> [owner-namd-l_at_ks.uiuc.edu <mailto:owner-namd-l_at_ks.uiuc.edu>] on
> behalf of Hannes Loeffler [Hannes.Loeffler_at_stfc.ac.uk
> <mailto:Hannes.Loeffler_at_stfc.ac.uk>]
> Sent: Tuesday, January 10, 2017 9:59 AM
> To: namd-l_at_ks.uiuc.edu <mailto:namd-l_at_ks.uiuc.edu>
> Subject: Re: namd-l: Solvation free energy difference of
> alchemical transformation of K+ into Na+
>
> On Tue, 10 Jan 2017 15:09:52 +0000
> "B.W.J. Irwin" <bwji2_at_cam.ac.uk <mailto:bwji2_at_cam.ac.uk>> wrote:
>
> > This .pdf is useful (sorry if you already know this)
> >
> http://www.ks.uiuc.edu/Training/Tutorials/namd/FEP/tutorial-FEP.pdf
> <http://www.ks.uiuc.edu/Training/Tutorials/namd/FEP/tutorial-FEP.pdf>
> >
> > Particularly the heading:
> > Why a soft-core potential should always be included
>
> Except that this here may very well be a case where that should not be
> necessary. The dual-topology approach in NAMD means that there is
> always a part vanishing and one appearing. Softcores however are
> really only effective/needed at short distances(*) but Abhishek seems
> to keep both ions at the same coordinates. This means there may
> always
> be enough excluded volume available to ensure water molecules not to
> come too close to the ion(s). Actually, I would think that a
> dual-topology approach in this paticular case is not really a clever
> idea because a simple linear transformation of the three non-bonded
> parameters in a single-topology approach would not require the
> restraint (there is already the constraint of a common coordinate
> set).
> But that's how FEP/TI is currently implemented in NAMD(**).
>
> (*) The code appears to have some WCA functionality included but that
> does not seem to be documented. It is functional?
> (**) I heard there were plans to work on this?
>
>
> > If you haven't included one.
> >
> > The fepout fileshould show the contribution from each λ-window
> > "Delta Delta A" and a measure of uncertainty.
>
> Abishek said TI not FEP.
>
>
> > If a particular window has a large contribution you could try
> > splitting it into two smaller windows.
> >
> >
> >
> > On 2017-01-10 12:15, Hannes Loeffler wrote:
> > > It would be useful it you kept the discussion in one thread.
> > >
> > > Your attachment suggests that you plot λ vs <U> but your
> really need
> > > <∂U/∂λ>. Is that what it is? How have you extracted this?
> > >
> > > Maybe you send all relevant input files (control file, parmfile,
> > > coordinates/alchfile) and output files (alchOutFile) so that
> > > somebody on
> > > the list can take a look.
> > >
> > >
> > > On Tue, 10 Jan 2017 16:39:01 +0530
> > > Abhishek Kumar Singh <bo13m1002_at_iith.ac.in
> <mailto:bo13m1002_at_iith.ac.in>> wrote:
> > >
> > >> Hello,
> > >>
> > >> I have alchemically transformed K+ to Na+ at the center of water
> > >> boxes (40, 60 and 80 A edge length respectively) using NAMD
> > >> software and CHARMM36 forcefield. A harmonic spring of force
> > >> constant = 10 kcal/mol.A^-2 has been applied throughout the
> > >> transformation. Temp and Pressure has been kept at 310K and 1 atm
> > >> respectively. Free energy of alchemical transformation has been
> > >> calculated using thermodynamic integaration. The derivative has
> > >> been computed as (∂G/∂λ) = ⟨ U(λ=1)- U)(λ=0) ⟩λ, where λ is the
> > >> coupling parameter such that λ=1 and λ=0 are K+ and Na+ vdw
> > >> parameters respectively. Free energy derivatives has been
> computed
> > >> for 11 λvdw values ( 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2,
> > >> 0.1 and 0.0). The area under the plot (derivate vs λ) gives me =
> > >> -23.40, -23.46 and -23.4 8
> > >> kcal/mol for the box sizes of 40A, 60A and 80A respectively. I
> > >> am attaching alongwith the derivate vsλ values for 40A water box.
> > >>
> > >> The standard CHARMM36 prm file actually states that the value
> > >> should be -18.6 kcal/mol.
> > >> So what can be reason for this difference. Is there any
> problem in
> > >> this procedure ?
> > >> Your suggestions are most welcome.
> > >>
> > >> Sincerely
> > >> Abhishek
> >
>
>
>
>
>
>
>
> --
> Thanks & Regards,
> Abhishek
-- Brian Radak Postdoctoral Appointee Leadership Computing Facility Argonne National Laboratory 9700 South Cass Avenue, Bldg. 240 Argonne, IL 60439-4854 (630) 252-8643 brian.radak_at_anl.gov
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:01 CST