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

[Berk Hess]

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.