g_order is incorrect for unsaturated carbons
The problem (which exists in 4.5.4), and a patch, are located here
I haven't checked 4.6, but there is no mention of a fix to g_order in the 4.6 release notes:
Please note that this patch does not fix g_order for all options, but should be useful for anybody who wants to fully implement the fix.
I realize that it is most desirable that I upload a fully fixed version of g_order, but I don't have time at the moment. I have therefore opened this issue so that at least it doesn't get forgotten about. A quick fix would be to have the current code throw an error when encountering unsaturated carbons, since these values are wrong (though widely reported in the literature based on g_order).
#1 Updated by Reid Van Lehn almost 6 years ago
I encountered the bug that Chris refers to here and wanted to chime in on its origin in the existing g_order code. The source code as stands works by following an approximation for Scd nicely derived in this paper: http://scitation.aip.org/content/aip/journal/jcp/109/6/10.1063/1.476823. Namely, for saturated carbons, the code defines a new set of molecular axes on the carbon atom C(N) for which the order parameter is calculated, defining a z-axis from C(N-1) to C(N+1), a y-axis that bisects the H-C-H bond angle, and a x-axis orthogonal to the C(N-1)-C(N)-C(N+1) plane. It then makes use of the relation of Scd = 2/3 Sxx + 1/3 Syy, where Scd is the desired order parameter that captures the orientation of the C-D vector for this carbon with respect to a reference axis and Sxx and Syy are order parameters reflecting the orientation of the molecular X and Y axes defined above relative to the reference axes. However, the prefactors in this expression (2/3 and 1/3) emerge from the assumed tetrahedral geometry of the central carbon and are related to the approximate cos^2 and sin^2 of 109.5 degrees, the presumed H-C-H bond angle. This expression is thus NOT valid for unsaturated carbons due to a different bond angle geometry. Similarly, the expression Scd = 2/3 Sxx + 1/3 Syy also implicitly cancels out a cross term related to Sxy due to the assumed equivalence of two hydrogen/deuterium atoms bonded to the central carbon; again this cross term cannot be ignored for an unsaturated bond with only a single bonded hydrogen.
The correct case for handling the an unsaturated bond is derived in the reference above, leading to the relation: Scd = 1/4 Szz + 3/4 Syy - sqrt(3)/2 Syz, where the signs in front of the different terms depends on the geometry of the double bond and how the axis Sz is defined. These prefactors assume an idealized 120 degree bond angle appropriate for the sp2 carbons included in a double bond. In the version of the g_order code I looked at (v. 4.6.1), the molecular axis Sz is correctly redefined for unsaturated bonds for the C9 carbon ONLY, but the incorrect Scd expression is then used. The reason C9 results look so wrong compared to C10 (in DOPC/POPC) is likely because the errors of incorrectly defining the Sz axis for C10 and using the incorrect Scd definition lead to a fortuitous cancellation of errors, yielding somewhat correct results for C10 but not C9. I've patched my own local version of g_order to use the correct expression for Scd and am able to reproduce results similar to those of Poger and Mark (http://pubs.acs.org/doi/abs/10.1021/ct900487a) and similar to the POPC results reported by Piggot et al (http://pubs.acs.org/doi/abs/10.1021/ct3003157). I haven't looked at Chris' patch which I'm sure works fine, but the solution I'm suggesting is in the spirit of the original implementation of g_order and thus might be easy to implement.
The solution to this fix is thus to fix the definition of Scd for unsaturated carbons, and ensure that the molecular axes are always correct - specifically, Sz always needs to be parallel to the double bond (currently only correct for the C9 carbon as defined in the GROMOS force field, for example). I have a heavily modified version of g_order that makes these changes amongst others, but is probably too different from the current version to be useful as is. If anyone wants I can clean up the code to remove unnecessary changes and upload it to at least serve as a template for fixing this bug in the future.
#2 Updated by David van der Spoel almost 6 years ago
- Affected version set to 4.5.5
thanks for loking into this issue. your description seems to be rather discouraging as we can not have tools that implicitly assume we are using a certain force field, and it even sounds like the lipid types are hardcoded. we used to have more such tools but are trying to phase them out. doescthe paper you mention give any formulae that are not dependent on knowing the geometry?
#3 Updated by Erik Lindahl about 5 years ago
David: g_order already has an option -unsat that should be used for unsaturated carbons, and the fix here seems to only concern that case. However, the present patch is a bit too dirty to apply right away, so this will have to wait until somebody has time to clean it up (or rewrite g_order properly).
#7 Updated by Thomas Piggot over 2 years ago
Some of my recent work has been looking at comparing different tools for the calculation of these order parameters. The modified version of g_order here from Chris still doesn't perform this calculation correctly for the double bond. Reid sent me his code and with some slight modifications I have a version that works correctly for these united-atom unsaturated double bonds. There are also some further modifications that calculate the individual hydrogen atom order parameters in united-atom methylene groups.
As per the comments above by Reid, this version of g_order: is quite different to the original version (e.g. it takes in a .dat file that defines things for the calculation a la g_density/gmx density); is just designed to calculate deuterium order parameters (both for a complete membrane and radially outwards from an object); and it is based upon g_order 4.6.x and so is still written in C (and my coding skills are not good enough to easily make the conversion to C++). However, I can provide this version here, if that will help towards solving this bug?