Re: LJcorrection option in NAMD

From: JC Gumbart (gumbart_at_physics.gatech.edu)
Date: Mon Jun 27 2016 - 00:23:12 CDT

Hi Brian,

> 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.

There has a been a community effort surrounding this for a few years now, although incomplete in some respects: http://www.alchemistry.org/wiki/Test_System_Repository <http://www.alchemistry.org/wiki/Test_System_Repository>

For example, no NAMD files (because none of us have found the time to contribute…).

Thanks for your great work! Development is HARD, for myriad reasons.

Best,
JC

> On Jun 22, 2016, at 5:22 PM, Brian Radak <bradak_at_anl.gov> wrote:
>
>
>
> 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.
>
> Regards,
> Brian
>
>>
>> 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>
>> http://branniganlab.wordpress.com <http://branniganlab.wordpress.com/>
>>
>>
>>
>>
>>
>>
>>
>> On Tue, Jun 21, 2016 at 10:48 AM, Jérôme Hénin <jerome.henin_at_ibpc.fr <mailto:jerome.henin_at_ibpc.fr>> 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 < <mailto:brian.radak.accts_at_gmail.com>brian.radak.accts_at_gmail.com <mailto:brian.radak.accts_at_gmail.com>> 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: radak_at_uchicago.edu
>>
>> --
>> <mailto:radak_at_uchicago.edu>Grace Brannigan, Ph.D.
>> Assistant Professor
>> Center for Computational and Integrative Biology (CCIB) &
>> Department of Physics
>> Rutgers University, Camden, NJ
>> (856)225-6780
>> <mailto:radak_at_uchicago.edu>http://branniganlab.wordpress.com <http://branniganlab.wordpress.com/>
> --
> 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>

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