From: Radak, Brian K (bradak_at_anl.gov)
Date: Tue Jan 10 2017 - 14:28:57 CST
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
>
This archive was generated by hypermail 2.1.6 : Mon Dec 31 2018 - 23:20:00 CST