[gmx-developers] problems with parameter assignement to CSH and COH groups when "-vsite hydrogen" flag is used

[Karmen Condic-Jurkic]

Hi everyone,


I'm a new Gromacs user and I bumped into a problem with the parameter assignment to the constraint bonds in COH and CSH groups in serines and cysteines when flag -vsite hydrogens is being used in the topology preparation. Namely, after running an energy minimisation on the system built with virtual sites, I noticed that CSH angles in cysteines ended up with a strange pointy form corresponding to 75 degrees, instead of the equilibrium value of 96. After a bit of a digging, it turned out that the problem is created during the grompp run, where only a single function is created to represent the constraint bond that prevents angle bending for both COH and CSH groups. There parameter used in that particular function corresponds to the COH group (0.198 nm), which is too short for CSH (should be 0.237 nm) and results with the distorted CSH angle. The reason why this happens seems to be the vague definition of the constraint bond given through the atom types (CH2 and H), which are unfortunately the same for both CYS and SER. In addition, the CH2-H distance parameter for CSH is actually completely omitted from the current ffG54a7bon.file! There is no problem with other COH groups, which are defined through different atom types, as shown below in the excerpt from ffG54a7bon.itp file with parts of interest highlighted:

 *************************************************************************
; get the constraint distances for dummy atom constructions

#include "ff_dum.itp"

[ constrainttypes ]
; now the constraints for the rigid NH3 groups
 MNH3    C    2   DC_MNC1
 MNH3  CH1    2   DC_MNC2
 MNH3  CH2    2   DC_MNC2
 MNH3 MNH3    2   DC_MNMN
; and the angle-constraints for OH and SH groups in proteins:
  CH2    H    2   DC_CO
  CH1    H    2   DC_CO
    C    H    2   DC_CO
    P    H    2   DC_PO
  ????????????????????? missing DC_CS value
**************************************************************************

I'm not really sure if this belongs to the developers or to the users mailing list, but I gues that depends on the solution of the problem. It can be probably done on a more general level by redefining the constraints using more unique atom descriptors or just by calculating the atom distance for each constraint of interest by looking at the atoms directly and not reading in a fixed value. As I have no experience in coding in C, this goes way beyond my programming skills so I came up with an easier solution - introducing additional atom type for hydrogen bound to sulphur named HS. With this, the grompp was able to override this issue and this seems to be working:

**************************************************************************
; get the constraint distances for dummy atom constructions

#include "ff_dum.itp"

[ constrainttypes ]
; now the constraints for the rigid NH3 groups
 MNH3    C    2   DC_MNC1
 MNH3  CH1    2   DC_MNC2
 MNH3  CH2    2   DC_MNC2
 MNH3 MNH3    2   DC_MNMN
; and the angle-constraints for OH and SH groups in proteins:
  CH2    H    2   DC_CO
  CH1    H    2   DC_CO
    C    H    2   DC_CO
    P    H    2   DC_PO
  CH2   HS    2   DC_CS

************************************************************************


I duplicated all the parameters given for atom type H and reassigned them to HS type. The modified force field files (version G54a7) are in the attachment and have suffix _HS.

I was trying to find if this problem was already known, but I couldn't find anything related to this issue. In any case, it would be good to fix this as many proteins contain cysteines, but in case it has already been fixed, I apologize for spamming.

All the best,

Karmen



-----------------------------
Dr. Karmen Condic-Jurkic
Postdoctoral Researcher
School of Chemistry and Molecular Biosciences
The University of Queensland
Brisbane St. Lucia, QLD 4072
Australia
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20140115/fbea5b5b/attachment.html>