How to compute the order parameter for the last carbon in lipid tail in membrane GROMACS?

User discussions
1

GROMACS version:5.1.2
GROMACS modification: No
In gromacs gmx order computes the order parameter per atom for carbon tails. For atom i the vector i-1, i+1 is used together with an axis. then we use this S mol = 1/2(3cos^2 (θ ) − 1),.

So to compute the molecular order parameter for atom for example C3. We need to have indexes for atom C2 and C4 (because we write a vector and check the angle).

But what we should do, when we have the last atom in the tail. For example, I have oleic acid and I want to compute the order parameter for C18? How to do this? I can make an index for C17, but I haven’t C19 (because this atom doesn’t exist in oleic acid), so how can I write a vector? Maybe I should make an index for H atom at the end of the tail???

I ask you because I want to compare results from the simulation of my membrane POPC to these from 1 H– 13 C solid-state nuclear magnetic resonance (NMR) spectroscopy and there they check C-H bond (so they don’t use vector Ci-1 Ci+1 for Ci), so they are able to check this parameter for the last carbon of the tail.

I think its impossible to compute last carbon, because “Note that the order program, by default, numbers the carbon atoms in each chain from 1. It is not possible to compute -SCD for terminal carbon atoms, as they lack neighboring atoms from which the local molecular axis is computed.” http://www.mdtutorials.com/gmx/membrane_protein/09_analysis.html
But I see that people use GROMACS and in the article they have the molecular order parameter for terminal carbon like here https://pubs.rsc.org/en/content/articlehtml/2013/cp/c2cp42738a?casa_token=my-RH0Q6BGwAAAAA:OeeBrgfrlOdwLVUHd7a6YLSuO0xvy2erDtADtmoN4x_nqXmrUwnBE9pcyholDOsPNwTEGXpheDYiTLU (Fig 3 red dots)

And I also have a question.

why in NMR they use CH vector and then in GROMACS they use vector Ci-1 - Ci+1? And they have similar results I don’t get it.

Thank you in advance for your help.

2

Hi,

There are a few different ways to compute the carbon-deuterium order parameters. The one implemented in gmx order is designed for united-atom force fields, so when you don’t have explicit hydrogen atoms in the tails and have to somehow predict where they would be to do the calculation. Because of the way it does this, you can’t use gmx order for the terminal methyl groups.

If you’re using an all-atom force field, you shouldn’t really use gmx order anyway, rather you should use a method where the actual position of the explicit hydrogen atoms are taken into account in the calculation. That said, the differences are fairly small for saturated carbon atoms but are larger for unsaturated carbons. You also get automatic averaging of the two different hydrogen atom (pro-R and pro-S) order parameters in a methylene group with gmx order that might not be desirable.

If you’re using a united-atom force field and want the order parameters for this terminal group, you should use a method that builds the hydrogen atoms back into the structure (e.g. protonate, babel, pdb2gmx, etc.) and then use an all-atom analysis method on that.

My paper from a few years ago has lots of details about all this ,including different tools available:

https://pubs.acs.org/doi/abs/10.1021/acs.jctc.7b00643

Also worth stressing is the gmx order doesn’t do the calculation correctly for unstaturated carbons (even with the -unsat option).

Cheers

Tom

3

Thank you so much.