Re: Solvation free energy difference of alchemical transformation of K+ into Na+

From: Abhishek Kumar Singh (bo13m1002_at_iith.ac.in)
Date: Wed Jan 11 2017 - 00:39:11 CST

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> 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
>
> ________________________________________
> From: owner-namd-l_at_ks.uiuc.edu [owner-namd-l_at_ks.uiuc.edu] on behalf of
> Radak, Brian K [bradak_at_anl.gov]
> Sent: Tuesday, January 10, 2017 12:27 PM
> To: 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
>
> ________________________________________
> From: owner-namd-l_at_ks.uiuc.edu [owner-namd-l_at_ks.uiuc.edu] on behalf of
> Hannes Loeffler [Hannes.Loeffler_at_stfc.ac.uk]
> Sent: Tuesday, January 10, 2017 9:59 AM
> To: 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> wrote:
>
> > This .pdf is useful (sorry if you already know this)
> > 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> 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

This archive was generated by hypermail 2.1.6 : Sun Dec 31 2017 - 23:20:59 CST