Re: LJcorrection option in NAMD

From: Brian Radak (
Date: Wed Jun 22 2016 - 17:04:02 CDT

On 06/21/2016 09:48 AM, Jérôme Hénin wrote:
> Brian, thanks for outlining those issues, in particular missing docs.
> I can't recommend enough that the docs be updated at the same time as
> the features themselves. It's vital for users of development versions,
> and it's the best time to do it as a developer as well - not least, it
> removes the risk that some things never get documented at all, like
> the WCA features.
This is on the TODO list. It is, however, a bit of double duty when one
writes new code alongside undocumented (even in the source code) features.

> Interactions with the pressure control, both through the kinetic
> pressure term for vanishing atoms (an issue you mentioned earlier) or
> long-range corrections, are among the subtleties that make it so
> complicated to maintain a solid alchemy code. Phrases like "should
> work more or less ok now" are pretty worrying for a user to read. They
> are also a reminder that we lack test cases for this code.
The virial debate is one that I raised in a previous thread and seems to
be more a confusion in the construction of "dual-topology" thermodynamic
cycles than anything else. I don't think anyone is really advocating
wide use of NpT alchemy anyway.

I should apologize for the loose language and I can see how that was a
bit alarming. This is, however, to be contrasted with "doesn't work as
expected", which is what the LJcorrection code was doing before the
updates. The "more or less" comes from the fact that compounded rounding
errors mean that the alchemical code doesn't exactly reproduce the
expected alchemical endpoint (agreement within <0.1 kcal/mol though).
This is suboptimal, but well within the reasonable statistical error of
free energy calculations. I also would expect that any serious
researcher doing free energy calculations requiring high precision would
be able enough to check this (that's how I originally came across the
unexpected behavior).

> I know that you need to get some science done, but I hope you can find
> some time to clean up the existing code, and build tests to ensure its
> robustness.
> Jerome
> On 21 June 2016 at 14:36, Brian Radak <
> <>> wrote:
> On 06/20/2016 05:09 PM, Yue Qian wrote:
> Dear NAMD users,
> I am currently using NAMD to carry out free energy of binding
> calculations of a protein-ligand system. I have encountered
> several problems regrading the long-range LJ corrections. I
> have googled for certain issues but couldn't find relevant,
> informative discussions.
> This is what I found on the manual:
> LJcorrection: Description: Apply an analytical correction to
> the reported vdW energy and virial that is equal to the amount
> lost due to switching and cutoff of the LJ potential. The
> correction will use the average of vdW parameters for all
> particles in the system and assume a constant, homogeneous
> distribution of particles beyond the switching distance. See
> [60] for details (the equations used in the NAMD
> implementation are slightly different due to the use of a
> different switching function). Periodic boundary conditions
> are required to make use of tail corrections. LJcorrection as
> implemented is inconsistent with vdwForceSwitching.
> Here are my questions:
> 1. To turn on the long-range LJ correction, is "LJcorrection
> yes" is the only thing I need to do? (From the manual, I
> noticed that "LJcorrection as implemented is inconsistent with
> vdwForceSwitching." vdwForceSwitching was kept off.) I have
> tested this protocol in unbound state calculation with only
> the solute molecule and TIP3 water molecules with LJcorrection
> on and off. However, I did not observe significant difference,
> which was expected to be around 0.5 ~ 0.6 kcal/mol. But with
> "Ewald off", the calculation result was ~0.8kcal/mol more
> positive. This patten applies to different cutoffs, from 8 Ang
> to 12 Ang.
> Correct, "LJCorrection yes" or "LJCorrection on" will apply the
> correction. As of 2.11 there actually is now a vdwForceSwitching
> compatible correction, I just haven't updated the documentation
> yet (sorry!). If a warning message still states this
> incompatibility, then it is not present, because I also removed
> that message in the update.
> The new code also improves consistency in conjunction with
> alchemical endpoints, but might not work as expected with the
> (also mostly undocumented) WCA FEP code. Standard FEP and TI
> should work more or less ok now.
> In general, such a correction will not significantly change free
> energies unless the solute rather very large. Even then I might
> only expect sub kcal/mol changes.
> The main changes I expected when I did all of this was to get
> different densities during an NpT equilibration. In general the
> correction increases density by ~5%.
> 2. According to the manual, the LJcorrction option assumes a
> homogeneous, isotropic environment which only applies to
> unbound state calculation. To be able to apply to bound state
> calculation, heterogenous, nonisotropic environment needs to
> be considered. The cited paper (J. Phys. Chem. B, Vol. 111,
> No. 45, 2007) covered both scenarios. But GROMACS was used in
> their work. I was wondering if there's anyway to incorporate
> such methods in NAMD.
> The main options here are to apply a Lennard-Jones Ewald approach
> (not available in NAMD) or else to use the perturbation-like
> method by re-evaluating the energy with a larger cutoff. The
> latter is doable but has never been automated in NAMD to my
> knowledge (I suppose we could give this a try if there is interest).
> Any clue is appreciated. Thank you very much.
> Best,
> Yue
> --
> Yue Qian
> --
> Brian Radak
> Postdoctoral Scholar
> Gordon Center for Integrative Science, W323A
> Department of Biochemistry & Molecular Biology
> University of Chicago
> 929 E. 57th St.
> Chicago, IL 60637-1454
> Tel: 773/834-2812
> email: <>

Brian Radak
Postdoctoral Appointee
Leadership Computing Facility
Argonne National Laboratory
9700 South Cass Avenue, Bldg. 240
Argonne, IL 60439-4854
(630) 252-8643

This archive was generated by hypermail 2.1.6 : Tue Dec 27 2016 - 23:22:16 CST