[gmx-users] free energy with TIP4P bug fixed

[chris.neale at utoronto.ca]

I have done some tests that, for a single system, empirically indicate  
that the coulomb(SR) discrepancy occurs only with use of the free  
energy code and only when tip4p is in the same energygrps as any part  
of the molecule which is perturbed, even if the common atoms are not  
perturbed themselves.

This system is composed of:
SOL = 7180 tip4p
DPC = 49 detergents
DPN = perturbed detergent in which the charges are turned off.

general free energy mdp options in these runs:
init_lambda         = 0.00
sc_alpha            = 0.0
nsteps              = 0

Different mdp options and resulting Coulomb(SR) values are shown below:

## Reproduce the discrepancy:
free_energy=no; energygrps=System
-3.22653e+05
free_energy=yes; energygrps=System
-3.23064e+05

## Discrepancy is resolved with separate energygrps for tip4p,  
perturbed molecule, and rest of system:
free_energy=no; energygrps=SOL     DPC       DPN
-3.22653e+05
free_energy=yes; energygrps=SOL    DPC       DPN
-3.22653e+05

## Discrepancy is resolved when energygrps separate tip4p from the  
perturbed molecule, with the remainder of the system being flexible to  
join with either tip4p or the perturbed molecule:
free_energy=yes; energygrps=SOL         DPN_DPC
-3.22657e+05
free_energy=yes; energygrps=SOL_DPC     DPN
-3.22653e+05

## Reproduce the discrepancy when energygrps treat tip4p together with  
the perturbed molecule, even when the remainder of the system is  
separated:
free_energy=yes; energygrps=SOL_DPN       DPC
-3.23061e+05

## Reproduce the discrepancy when only 3 atoms of DPN are perturbed  
while those 3 perturbed atoms are in a separate energygrps section  
from tip4p, while tip4p is in the same energygrps section as the  
non-perturbed atoms from DPN (note that the perturbed .itp is  
different in this case in order to (i) perturb only 3 atoms, and (ii)  
to ensure that perturbed atoms are not in the same charge group as any  
non-perturbed atoms, which is required to pass grompp due to the  
energygrps selection):
free_energy=yes; energygrps=energygrps = SOL_DPC_DPN_&_!C1_C2_C7   
DPN_&_C1_C2_C7
-3.22886e+05

Again, this is not exhaustive in that it is only a single system, so I  
can not be sure that there are no other issues, but at least the  
discrepancies noted here probably hold for all systems.

Chris.

-- original message --

> Berk, David, and other developers,
>
> I believe that all the relevant information has already been posted,  
> but since I use TIP4P ubiquitously, I would like to get some  
> clarification regarding simulations that may be affected. I  
> currently understand that:
>
> 1. Utilization of TIP4P in the absence of the free_energy code is  
> currently believed to be without error in gromacs 3.3 through 4.0.3
> 2. Combining TIP4P and the free-energy code may lead to the errors noted.
> 3. Combining TIP4P and the free-energy code with at least 3 energy  
> groups (TIP4P, perturbed atoms, rest of system) will not lead to the  
> errors noted. This would only require 2 energy groups if the system  
> was composed of only TIP4P and the perturbed atoms.
>
> It would be very useful for me if somebody could validate the  
> statements above or state any incorrections in those statements.
>
> Thank you,
> Chris.