Re: LJcorrection option in NAMD

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

On 06/21/2016 01:03 PM, Grace Brannigan wrote:
> Hi Brian,
> As a decade-long user of the alchemy code, I certainly welcome
> improvements to consistency and accuracy, and appreciate you working
> to implement fixes.
> It is alarming, though, to hear that changes are being committed
> casually and without documentation. A few months ago, for instance,
> my postdoc Reza Salari observed discrepancies of ~30,000 kcal/mol in
> FEP output from the development version. That was hard to miss, but
> naturally makes us wonder about less obvious bugs.
Fortunately, although I don't think I outlined this well enough, the
error was actually quite small in its extent, since the core simulation
dynamics (i.e. the forces and energies) were still correct (that's what
my error testing did verify!). The only bug (now fixed) was in the FEP
/analysis/ code. Post processing of the trajectories would have lead to
the expected result (I should also note that multistate techniques such
as WHAM are now quite standard in the field).

Again, not to minimize the mistake, but the behavior before the bug was
actually quite a bit more pernicious in my opinion, because alchemical
scaling of bonded terms was entirely excluded and could (at least
heuristically) lead to unsatisfying behavior in certain circumstances.
There was, however, no indication from the code output that this kind of
situation was possible.

> Even undocumented fixes can introduce unexpected discrepancies
> compared to previous results, which can take a user some time to track
> down if they don't know where to look. You might know it's an
> improvement - but the user doesn't.
This is true and of course to be avoided. Unfortunately undocumented
features aren't always even apparent to all of the developers! Such are
the perils of academic software.

> So, I'd like to request that you be considerate to the numerous
> research groups using the FEP features when choosing to
> commit changes. And let me know if you need some unusual test cases...
Duly noted. Candidates for regression testing would be quite welcome, as
this is something that we need to move towards in order to prevent
problems in the future. A community based effort could go a long way in
making this happen and further boost the other ongoing developments to
the code base.


> Thanks,
> --
> Grace Brannigan, Ph.D.
> Assistant Professor
> Center for Computational and Integrative Biology (CCIB) &
> Department of Physics
> Rutgers University, Camden, NJ
> (856)225-6780 <tel:%28856%29225-6780>
> <>
> On Tue, Jun 21, 2016 at 10: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.
> 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.
> 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 <tel:773%2F834-2812>
> email:
> --
> <>
> Grace Brannigan, Ph.D.
> Assistant Professor
> Center for Computational and Integrative Biology (CCIB) &
> Department of Physics
> Rutgers University, Camden, NJ
> (856)225-6780
> <>

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