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

Berk Hess gmx3 at hotmail.com
Wed Oct 26 11:47:14 CEST 2005

>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 
>Date: Tue, 25 Oct 2005 18:00:47 +0200
>An error occurs in free energy calculations (TI) of heme-containing 
>systems, even if the heme moiety itself is not involved in the 
>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 
>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.
>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 
This will solve the HEME problem.
The change will be in 3.3.1 which will be released soon.


More information about the gromacs.org_gmx-users mailing list