[gmx-developers] issue on benchmarking FEP calculations in Gromacs and NAMD

[Wang, Yuhang]

Dear Berk,

The equation I showed was for the self interaction contribution.The toy simulation I just did is simply a test of the PME calculation in GROMACS. Since the expected result is zero, any non-zero number reported by GROMACS will be an artifact of PME. Apparently GROMACS does not apply such self interaction correction.  Otherwise I should get zero instead of 11 kcal/mol.

By the way, I just found that NAMD doesn't need correction only because I set Decouple=on, which means the lambda parameter for the ion I put into the box is not changed for interaction with itself. Therefore the self interaction term is always on.

Thanks for mentioning Philippe Hünenberger. I will certainly read about his work. But right now, I am only concerned about why the two simulation packages produce different RAW results.


Steven W.
________________________________________
From: Berk Hess [berkhess at gmail.com] on behalf of Berk Hess [hess at kth.se]
Sent: Monday, June 11, 2012 3:35 PM
To: Wang, Yuhang; Discussion list for GROMACS development
Subject: Re: [gmx-developers] issue on benchmarking FEP calculations in Gromacs and NAMD

Hi,

Please keep the discussion on the list.
Do you really mean the single charge self interaction or the uniform
background charge contribution?
Gromacs takes all of this into account, so there is nothing to correct.
But you probably want the solvation free energy, not the cost of adding
or removing a charge.
Then the results depends on the reference state.
If you are not aware of the work of Philippe Hünenberger, please read
it. There are about five
corrections terms required to obtain the proper solvation free energy.

Cheers,

Berk


On 06/11/2012 10:31 PM, Wang, Yuhang wrote:
> Dear Berk,
>
> Thank you for the reply. The correction  term I applied is the PME self interaction correction based on the following equation:
> 0.5*(2.84/L)*331 (kcal/mol)
>
> By doing a simple simulation of a single ion in vacuum, I just figured out this:
> NAMD doesn't need this correction, but GROMACS does.
>
> Strangely, the correction term needs to be a negative constant (the equation says it should be positive). Anyway, with this correction, both GROMACS and NAMD give ~82 kcal/mol for the electrostatic component of SOD desolvation. Puzzle solved!
>
>
> Steven W.
>
>
>
> ________________________________________
> From: gmx-developers-bounces at gromacs.org [gmx-developers-bounces at gromacs.org] on behalf of Berk Hess [hess at kth.se]
> Sent: Monday, June 11, 2012 12:21 PM
> To: Discussion list for GROMACS development
> Subject: Re: [gmx-developers] issue on benchmarking FEP calculations in Gromacs and NAMD
>
> Hi,
>
> What is the correction you applied?
> I think Gromacs has all charge corrections in, except for the set of
> corrections
> required for absolute solvation free energies as described by Philippe
> Hünenberger.
>
> Cheers,
>
> Berk
>
> On 6/11/12 18:23 , Wang, Yuhang wrote:
>> Dear Floris,
>>
>> I actually did the FEP in steps: step 1, turn off electrostatics; step 2, turn off vdW. Parameters alchElecLambdaStart and alchvdWLambdaEnd were modified correspondingly. So that is not the problem. The real problem is I got two different results (~ 10kcal/mol difference) in step 1 (turning off electrostatics), although step 2 results matched quite well. Simulations in both packages used the same force field parameters for ion and water (TIP3P). PME was used in both cases. Then the question is why two packages differ from each other.
>>
>>
>> Steven W.
>> ________________________________________
>> From: gmx-developers-bounces at gromacs.org [gmx-developers-bounces at gromacs.org] on behalf of Floris Buelens [floris_buelens at yahoo.com]
>> Sent: Monday, June 11, 2012 3:33 AM
>> To: Discussion list for GROMACS development
>> Subject: Re: [gmx-developers] issue on benchmarking FEP calculations in Gromacs and NAMD
>>
>> I think you may be inadvertently coupling the van der Waals parameter of the ion to lambda in your NAMD configuration. Check in your trajectories whether you're seeing water overlap with the ion as you get to the end point, and read up on the values of alchVdwLambdaEnd and alchElecLambdaStart - the concept is a little different than for Gromacs, but I think what you need are values of alchVdwLambdaEnd=0 and alchElecLambdaStart = 0 to get the effect you're looking for (only decouple the electrostatic component without changing Lennard-Jones).
>>
>>
>>
>> ________________________________
>> From: "Wang, Yuhang"<ywang148 at illinois.edu>
>> To: "gmx-developers at gromacs.org"<gmx-developers at gromacs.org>
>> Sent: Friday, 8 June 2012, 18:21
>> Subject: [gmx-developers] issue on benchmarking FEP calculations in Gromacs and NAMD
>>
>>
>>
>> Dear Gromacs developers,
>>
>>
>> I have benchmarked the desolvation free energy calculation of Na+ ion using Gromacs and NAMD (for comparison). There is a large difference in the electrostatic desolvation free energy (see below and the attached figure):
>>
>> Gromacs: 82.10(raw)+11.65(Ewald correction) = 93.75 kcal/mol
>> NAMD:     93.62(raw)+11.23(Ewald correction) = 104.85 kcal/mol
>>
>> Both of them are FEP calculations with perturbation of electrostatic interactions and used the same Lennard-Jones parameters:
>> sigma=0.243 nm, epsilon=0.196 kJ/mol (0.0469 kca/mol)
>>
>> Ewald correction was calculated by: 0.5*(2.837297/L)*331 (unit: kcal/mol), "L" is the cubic box length.
>>
>> Question: how can I explain the difference? Does Gromacs have a different PME implementation than NAMD?
>>
>> P.S. my input scripts are in the attachments.
>>
>>
>> Steven W.
>>
>> --
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-developers-request at gromacs.org.
>> --
>> gmx-developers mailing list
>> gmx-developers at gromacs.org
>> http://lists.gromacs.org/mailman/listinfo/gmx-developers
>> Please don't post (un)subscribe requests to the list. Use the
>> www interface or send it to gmx-developers-request at gromacs.org.
> --
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-developers-request at gromacs.org.