Project

General

Profile

Bug #953

fix for g_chi omega angle calculation

Added by Alexandar Hansen almost 5 years ago. Updated over 3 years ago.

Status:
Closed
Priority:
Normal
Assignee:
Category:
analysis tools
Target version:
Affected version - extra info:
Affected version:
Difficulty:
uncategorized
Close

Description

I've been struggling with this for a while and have been able to calculate the omega angle properly (Ca(i-1) - CO - N - Ca(i) ) with an index, but doing it all with g_chi is preferable. Searching the the source, I found a fix. I don't know who or where this should go to, but hopefully someone will include this in a future version.

/* Modified by Alex Hansen, 1 June 2012 */
/* first the old version */
/*for(i=0; (i<nl); i++) 
  {*/
      /* Omega */
/*      if (has_dihedral(edOmega,&(dl[i])))
      {
          dl[i].j0[edOmega] = n/4;
          id[n++]=dl[i].atm.minO;
          id[n++]=dl[i].atm.minC;
          id[n++]=dl[i].atm.N;
          id[n++]=dl[i].atm.H;
      }
  }*/
/* and here's the fix. Note: we start at 1, not 0 */
  for(i=1; (i<nl); i++) 
  {
      /* Omega */
      if (has_dihedral(edOmega,&(dl[i])))
      {
          dl[i].j0[edOmega] = n/4;
          id[n++]=dl[i-1].atm.Cn[1];
          id[n++]=dl[i-1].atm.C;
          id[n++]=dl[i].atm.N;
          id[n++]=dl[i].atm.Cn[1];
      }
  }

g_chi_compare-TYR.agr (17.9 KB) Alexandar Hansen, 06/04/2012 06:05 PM

g_chi_compare-PRO.agr (16.8 KB) Alexandar Hansen, 06/04/2012 06:05 PM

g_chi_compare-ALA.agr (17.9 KB) Alexandar Hansen, 06/04/2012 06:05 PM

Cropped0.xtc (5.02 MB) Alexandar Hansen, 04/29/2013 10:31 PM

Im7Ex2x.tpr (872 KB) Alexandar Hansen, 04/29/2013 10:31 PM


Related issues

Related to GROMACS - Bug #1235: peptide dihedral angle definitions violate IUPAC New 04/30/2013
Related to GROMACS - Feature #1247: fix hardcoded references to atom names in analysis tools New 05/10/2013

Associated revisions

Revision b0dd1667 (diff)
Added by Mark Abraham almost 4 years ago

Fix g_chi -omega

Per IUPAC, we should calculate peptide omega angles from
Calpha-C-N-Calpha, and there are only n-1 such angles defined for a
polypeptide of n amino acid residues.

There are some other IUPAC violations in make_chi_ind(), but solving
those is deferred to 5.0, not least because there will be a
behavioural change when we conform to IUPAC.

The problem that led to the report had several causes. The hard-coding
of the string names for hydrogen bound to peptide nitrogen (which in
CHARMM27 is "HN") meant that an H atom was not found. has_dihedral()
then checked whether a Calpha was found, which was true. Then
make_chi_ind() computed the angle based on H without doing any further
check. So some garbage coordinate was looked up when computing the
angle.

A proper fix for this kind of stuff for 5.0 requires finding all the
places we try to guess what atoms are from their names and using
properly implemented selections.

Also removed orphaned header file.

Fixes #953, details and links there

Change-Id: I63257a1d6d8ed46d251b8bc0cb308dfbd2864a97

History

#1 Updated by Alexandar Hansen almost 5 years ago

Wanted to add some 'before' and 'after' figures for the trajectory I was having trouble with.

Strangely, when I use g_angle and define the dihedral as O(i-1)-CO-N(i)-HN, I get good values ~180 degrees. I really don't understand what's going wrong with g_chi as it was before.

Also, in the fix above, the superscript ones are 1's in brackets. The post's formatting changed them into supers.

#2 Updated by Mark Abraham about 4 years ago

  • Target version changed from 4.5.6 to 4.5.7

Your patch probably works if the input conforms to your assumptions and the previous residue has a C-alpha and C atom, but there has to be checking done in @has_dihedral() that is consistent with it or code might simpley break. I'm also not sure why the old code doesn't work, but if you'd like to upload some (small, compressed) input files then we have a chance of seeing it go wrong in action and determine the proper fix.

#3 Updated by Mark Abraham almost 4 years ago

  • Status changed from New to Accepted
  • Assignee changed from David van der Spoel to Mark Abraham
  • Target version changed from 4.5.7 to future
  • Affected version set to 4.5.5

#4 Updated by Alexandar Hansen almost 4 years ago

I apologize for the long delay. My project diverted from Gromacs for some time and I find I need to refamiliarize myself with the issue. The attached .xtc was generated from a Camshift MD simulation (ie. with chemical shift restraints augmenting the force field). As such, there may be some peculiarities with the structures/proton positions.

Regardless, the issue was with how the omega angle was defined. Where previously it was being calculated with the O, C, N, and H atoms of the plane, it should be defined as described above (and with my fix) according to IUPAC-IUB recommendations. http://www.chem.qmul.ac.uk/iupac/misc/ppep3.html#320. This also turns out to be more robust, based on my previous uploads.

#5 Updated by Mark Abraham almost 4 years ago

  • Description updated (diff)

#6 Updated by Mark Abraham almost 4 years ago

  • Description updated (diff)

#7 Updated by Mark Abraham almost 4 years ago

  • Status changed from Accepted to In Progress

Alexandar Hansen wrote:

I apologize for the long delay. My project diverted from Gromacs for some time and I find I need to refamiliarize myself with the issue. The attached .xtc was generated from a Camshift MD simulation (ie. with chemical shift restraints augmenting the force field). As such, there may be some peculiarities with the structures/proton positions.

Regardless, the issue was with how the omega angle was defined. Where previously it was being calculated with the O, C, N, and H atoms of the plane, it should be defined as described above (and with my fix) according to IUPAC-IUB recommendations. http://www.chem.qmul.ac.uk/iupac/misc/ppep3.html#320. This also turns out to be more robust, based on my previous uploads.

You are correct about how these dihedrals should be defined. The definitions used in g_chi were wrong, and also buggy because of the way CHARMM names the H on peptide N as "HN" and the way much of the GROMACS code base was not written in anticipation of this. This ought to have been immaterial if the correct definitions of the dihedrals were used. Sigh :-) Mentioning CHARMM at any point would have led to a 5-minute fix :-)

https://gerrit.gromacs.org/2348 is based on your suggestion and fixes some of these issues.

There should be further house cleaning later, noted in the related issue.

#8 Updated by Alexandar Hansen almost 4 years ago

Ah, thanks! That was doing something even more devious than I imagined! I'm glad I was able to find a real issue, and I wasn't merely fixing something that wasn't broken, as I usually do :-)

#9 Updated by Mark Abraham almost 4 years ago

Ja, I have "fixed" my share of "uh, but that's a feature" problems :-) Thanks for the interest

#10 Updated by Mark Abraham almost 4 years ago

  • Status changed from In Progress to Fix uploaded

#11 Updated by Mark Abraham almost 4 years ago

  • Target version changed from future to 4.6.3

#12 Updated by Mark Abraham almost 4 years ago

  • Status changed from Fix uploaded to Resolved
  • % Done changed from 0 to 100

#13 Updated by Rossen Apostolov over 3 years ago

  • Status changed from Resolved to Closed

Also available in: Atom PDF