[gmx-users] Bug in free energy calculations in heme-containingsystems

[Berk Hess]

>From: Chris Oostenbrink <c.oostenbrink at few.vu.nl>
>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>To: gmx-users at gromacs.org
>CC: Chris Oostenbrink <c.oostenbrink at few.vu.nl>
>Subject: [gmx-users] Bug in free energy calculations in 
>heme-containingsystems
>Date: Tue, 25 Oct 2005 18:00:47 +0200
>
>Hello,
>An error occurs in free energy calculations (TI) of heme-containing 
>systems, even if the heme moiety itself is not involved in the 
>perturbation.
>In the gromos 43a1 force field, the heme iron atom is coordinated by four 
>nitrogen atoms (NA, NB, NC and ND), all in one plane and possibly by one 
>histidine Nitrogen, NE2, above this plane. Default force-field parameters 
>for the connection between His and Heme are defined in file ffG43a1bon.itp, 
>including the N[ABCD] - FE - NE2 angle of type ga_16.
>In a simulation of a heme-containing system, the parameters of the 
>intraheme N[ABCD] - FE - N[ABCD] angles are of type ga_1 (90 degrees) as 
>specified in the topology for state A. However, in the tpr-file generated 
>by grompp, these angles get the parameters of ga_16 (115 degrees) for state 
>B!!
>
>Effect
>Unperturbed simulations of such systems are probably unaffected, but as 
>soon as a thermodynamic integration calculation is performed two things 
>will happen (also if the Heme is not coordinated to a His): 1) The Heme 
>group will deform due to increasing angles around the iron, 2) a 
>corresponding free energy contribution of about 100 kJ/mol will be added to 
>the total free energy difference that is being calculated.
>
>Solution
>The reason for all this probably lies in grompp, where the standard 
>force-field parameters are incorporated into the system topology, resulting 
>in a faulty tpr-file. In the case where the Heme group is not coordinated 
>to a His, commenting out the default force-field parameters in file 
>ffG43a1bon.itp seems to help. Otherwise, explicitly stating all Heme 
>parameters in state B as well may also prevent the error from happening.
>
>Chris de Graaf and Chris Oostenbrink

This problem arises from the way grompp sets the B-state parameters.
The way it has been done up till now can lead to problems like this.
It also means that one had to set all B-state bonded parameters when
doing solvation free-energy studies.

Two weeks ago Erik has changed the way the B-state bonded parameters
are assigned. Default force-field parameters are now only used when
the A-state parameters are not set explicitly. When the A-state parameters
are specified in the molecule definition they will also be used for the 
B-state.
This will solve the HEME problem.
The change will be in 3.3.1 which will be released soon.

Berk.