What is gorder?

gorder is a command-line tool that aims to provide a comprehensive, reliable, simple, and fast way to calculate atomistic and coarse-grained lipid order parameters from Gromacs simulations.

The primary (and possibly naïve) goal of gorder is to completely replace the use of gmx order for atomistic and coarse-grained systems. gmx order is a tool still widely used in the Gromacs community despite its tedious usage and the well-known fact that it produces incorrect results for many systems [1].

gorder also offers many additional features, such as leaflet classification or construction of ordermaps.

gorder is written in Rust and is also available as a Rust crate.

Some theory

Lipid order parameters quantify the degree of ordering in lipid membranes. They can be obtained from both experimental measurements and simulations, making them quite useful.

The order parameter $S$ in atomistic simulations can be calculated as:

$S = \frac{1}{2} ⟨ 3 cos^{2} θ - 1 ⟩$

where $θ$ is the angle between the C-H bond vector and the bilayer normal, and the angular brackets represent molecular and temporal ensemble averages. The order parameter for atomistic systems is usually reported as the negative of this calculated value, $- S$ (often denoted as $- S_{C H}$ or $- S_{C D}$ ). Higher positive values of $- S$ indicate more ordered membrane structures. $- S$ of 0 corresponds to a disordered membrane.

For coarse-grained systems, the order parameter can be calculated using the same equation. However, here $θ$ corresponds to the angle between the vector connecting two consecutive beads of the lipid and the bilayer normal. Using this definition, the order parameter $S$ increases as the membrane structure becomes more ordered. $S$ of 0 corresponds to a disordered membrane. For coarse-grained systems, the order parameter is usually reported directly as $S$ to facilitate comparison with the atomistic $- S$ (which follows the same trend).

Installation

Installing the gorder tool consists of two simple steps:

Install Rust.
- Follow this guide for your operating system.
- If you already have Rust installed, you can skip this step. To check, run rustc --version.
Install gorder.
- Run cargo install gorder in your terminal.

You can verify that the installation was successful by running gorder --version. The current version of the gorder tool should be displayed in the terminal.

Atomistic order parameters

Preparing the input

Suppose we have a membrane composed of POPC, POPE, and POPG lipids.

To calculate atomistic order parameters, we need two Gromacs files:

A TPR file containing the system structure and topology (system.tpr).
An XTC trajectory file (md.xtc) whose frames will be analyzed.

Next, we create an input YAML file that specifies the options for the analysis:

structure: system.tpr
trajectory: md.xtc
analysis_type: !AAOrder
  heavy_atoms: "@membrane and name r'C3([2-9]|1[0-6])|C2([2-9]|1[0-8])'"
  hydrogens: "@membrane and element name hydrogen"
output: order.yaml

In the input YAML file, the analysis type AAOrder requires you to specify both heavy_atoms and hydrogens. gorder will then identify all bonds connecting the selected heavy atoms with the selected hydrogen atoms. The order parameters are calculated for all these identified bonds.

The atoms are selected using a query language called GSL. If you are familiar with the query language used in VMD, you'll find the basic syntax of GSL intuitive.

Here:

heavy_atoms are selected using the query @membrane and name r'C3([2-9]|1[0-6])|C2([2-9]|1[0-8])', which selects all palmitoyl and oleoyl carbons of the membrane lipids.
hydrogens are selected using the query @membrane and element name hydrogen.

@membrane is a GSL autodetection macro that selects all atoms of common membrane lipids. The block starting with r' is a regular expression which are also supported by GSL.

The results of the analysis will be saved in the order.yaml file.

Running the analysis

We save the input YAML file, for example, as analyze.yaml. Then, we run gorder as follows:

$ gorder analyze.yaml

During the analysis, we will see something like this (except colored):

>>> GORDER v0.2.0 <<<

[*] Read config file 'analyze.yaml'.
[*] Will calculate all-atom order parameters.
[*] Membrane normal expected to be oriented along the z axis.
[*] Read molecular topology from 'system.tpr'.
[*] Detected 8768 heavy atoms using a query '@membrane and name r'C3([2-9]|1[0-6])|C2([2-9]|1[0-8])''.
[*] Detected 21592 hydrogen atoms using a query '@membrane and element name hydrogen'.
[*] Detecting molecule types...
[*] Detected 3 relevant molecule type(s).
[*] Molecule type POPE: 64 order bonds, 131 molecules.
[*] Molecule type POPC: 64 order bonds, 128 molecules.
[*] Molecule type POPG: 64 order bonds, 15 molecules.
[*] Will read trajectory file 'md.xtc' (start: 0 ps, end: inf ps, step: 1).
[*] Performing the analysis using 1 thread(s)...
[COMPLETED]   Step    225000000 | Time       450000 ps
[*] Writing the order parameters into a yaml file 'order.yaml'...
[✔] ANALYSIS COMPLETED

Note that the structure from the TPR file is not analyzed. The TPR file is only used to construct the system and obtain its topology.

The results of the analysis are saved in the order.yaml file. Here is an excerpt from the file:

# Order parameters calculated with 'gorder v0.2.0' using structure file 'system.tpr' and trajectory file 'md.xtc'.
- molecule: POPE
  order:
    POPE C22 (23):
      total: 0.1098
      bonds:
        POPE H2R (24):
          total: 0.0946
        POPE H2S (25):
          total: 0.125
    POPE C32 (32):
      total: 0.2341
      bonds:
        POPE H2X (33):
          total: 0.2438
        POPE H2Y (34):
          total: 0.2244
  # (...)
- molecule: POPC
  order:
    POPC C22 (32):
      total: 0.1094
      bonds:
        POPC H2R (33):
          total: 0.0953
        POPC H2S (34):
          total: 0.1235
    POPC C32 (41):
      total: 0.2325
      bonds:
        POPC H2X (42):
          total: 0.2405
        POPC H2Y (43):
          total: 0.2245
  # (...)
- molecule: POPG
  order:
    POPG C22 (25):
      total: 0.1081
      bonds:
        POPG H2R (26):
          total: 0.1024
        POPG H2S (27):
          total: 0.1138
    POPG C32 (34):
      total: 0.2028
      bonds:
        POPG H2X (35):
          total: 0.2131
        POPG H2Y (36):
          total: 0.1926

gorder automatically identified three molecule types and all relevant bonds. Order parameters are reported separately for each molecule type: for each bond type of each molecule type and for each heavy atom type of each molecule type. Order parameters for heavy atom types are obtained by averaging the order parameters of their bonds with hydrogens.

Let's take a closer look at a part of the YAML file:

- molecule: POPE         # name of the molecule
  order:                 # order parameters
    POPC C22 (23):       # heavy atom type: residue atom_name (relative_index)
      total: 0.1098      # order parameter of the heavy atom
      bonds:             # bonds of this heavy atom with hydrogens
        POPC H2R (24):   # hydrogen type: residue atom_name (relative_index)
          total: 0.0946  # order parameter of bond with this hydrogen
        POPC H2S (25):   # hydrogen type: residue atom_name (relative_index)
          total: 0.125   # order parameter of bond with this hydrogen

YAML files are easy to read programmatically and not completely human-unreadable. However, gorder also provides other output formats (XVG, CSV, human-readable table). See Output Formats for more information.

Using groups from an NDX file

gorder also supports using groups from NDX files. Let's suppose our NDX file already contains a group called TailCarbons, which specifies all the carbons to use. We no longer need to specify them using the complex regular expression from the previous example and can simply select the group. Our input YAML file will look like this:

structure: system.tpr
trajectory: md.xtc
index: index.ndx     # name of the NDX file
analysis_type: !AAOrder
  heavy_atoms: "TailCarbons"   # name of a group from NDX file
  hydrogens: "@membrane and element name hydrogen"
output: order.yaml

Then, we run the analysis in the same way as before.

Coarse-grained order parameters

Preparing the input

Suppose we have a Martini 3 membrane composed of POPC, POPE, and POPG lipids.

To calculate coarse-grained order parameters, we need two Gromacs files:

A TPR file containing the system structure and topology (system.tpr).
An XTC trajectory file (md.xtc) whose frames will be analyzed.

Next, we create an input YAML file that specifies the options for the analysis:

structure: system.tpr
trajectory: md.xtc
analysis_type: !CGOrder
  beads: "@membrane"
output: order.yaml

In the input YAML file, the analysis type CGOrder requires you to specify beads that will be considered for the analysis. gorder will then identify all bonds connecting the selected beads. The order parameters are calculated for all these identified bonds.

The atoms are selected using a query language called GSL. If you are familiar with the query language used in VMD, you'll find the basic syntax of GSL intuitive.

Here beads are selected using the query @membrane. That is a GSL autodetection macro that select all beads or atoms of common membrane lipids.

The results of the analysis will be saved in the order.yaml file.

Running the analysis

We save the input YAML file, for example, as analyze.yaml. Then, we run gorder as follows:

$ gorder analyze.yaml

During the analysis, we will see something like this (except colored):

>>> GORDER v0.2.0 <<<

[*] Read config file 'analyze.yaml'.
[*] Will calculate coarse-grained order parameters.
[*] Membrane normal expected to be oriented along the z axis.
[*] Read molecular topology from 'system.tpr'.
[*] Detected 6096 beads for order calculation using a query '@membrane'.
[*] Detecting molecule types...
[*] Detected 3 relevant molecule type(s).
[*] Molecule type POPC: 11 order bonds, 242 molecules.
[*] Molecule type POPE: 11 order bonds, 242 molecules.
[*] Molecule type POPG: 11 order bonds, 24 molecules.
[*] Will read trajectory file 'md.xtc' (start: 0 ps, end: inf ps, step: 1).
[*] Performing the analysis using 1 thread(s)...
[COMPLETED]   Step     50000000 | Time      1000000 ps
[*] Writing the order parameters into a yaml file 'order.yaml'...
[✔] ANALYSIS COMPLETED

Note that the structure from the TPR file is not analyzed. The TPR file is only used to construct the system and obtain its topology.

The results of the analysis are saved in the order.yaml file. Here is an excerpt from the file:

# Order parameters calculated with 'gorder v0.2.0' using structure file 'system.tpr' and trajectory file 'md.xtc'.
- molecule: POPC
  order:
    POPC NC3 (0) - POPC PO4 (1):
      total: -0.1362
    POPC PO4 (1) - POPC GL1 (2):
      total: 0.585
    POPC GL1 (2) - POPC GL2 (3):
      total: -0.1808
    POPC GL1 (2) - POPC C1A (4):
      total: 0.3835
  # (...)
- molecule: POPE
  order:
    POPE NH3 (0) - POPE PO4 (1):
      total: -0.1293
    POPE PO4 (1) - POPE GL1 (2):
      total: 0.6072
    POPE GL1 (2) - POPE GL2 (3):
      total: -0.1833
    POPE GL1 (2) - POPE C1A (4):
  # (...)
- molecule: POPG
  order:
    POPG GL0 (0) - POPG PO4 (1):
      total: 0.0004
    POPG PO4 (1) - POPG GL1 (2):
      total: 0.5917
    POPG GL1 (2) - POPG GL2 (3):
      total: -0.1807
    POPG GL1 (2) - POPG C1A (4):
      total: 0.395
  # (...)

gorder automatically identified three molecule types and all relevant bonds. Order parameters are reported for each bond type of each molecule type.

Let's take a closer look at a part of the YAML file:

- molecule: POPC                  # name of the molecule
  order:                          # order parameters
    POPC NC3 (0) - POPC PO4 (1):  # bond type (atom type 1 - atom type 2)
      total: -0.1362              # order parameter of this bond
    POPC PO4 (1) - POPC GL1 (2):  # atom types are specified as residue atom_name (relative_index)
      total: 0.585
    POPC GL1 (2) - POPC GL2 (3):
      total: -0.1808
    POPC GL1 (2) - POPC C1A (4):
      total: 0.3835

Using groups from an NDX file

gorder also supports using groups from NDX files. Let's suppose our NDX file already contains a group called LipidBeads, which specifies all the beads to use. We can then simply select the group. Our input YAML file will look like this:

structure: system.tpr
trajectory: md.xtc
index: index.ndx        # name of the NDX file
analysis_type: !CGOrder
  beads: "LipidBeads"   # name of a group from NDX file
output: order.yaml

Then, we run the analysis in the same way as before.

Output formats

The YAML output is always generated by gorder, as it contains all the calculated information. See Atomistic Order Parameters or Coarse-grained Order Parameters for more details about this format. However, gorder can also write output in other formats.

CSV format

CSV format may be useful if you want to work with the results using spreadsheet software such as LibreOffice Calc or Microsoft Excel, or when working with data analysis libraries like Pandas or Polars.

To write the results of the analysis in CSV format, add output_csv: PATH_TO_OUTPUT to your input YAML file.

Here is an excerpt from a CSV file:

molecule,residue,atom,relative index,total,hydrogen #1,hydrogen #2,hydrogen #3
POPE,POPE,C22,23,0.1098,0.0946,0.1250,
POPE,POPE,C32,32,0.2341,0.2438,0.2244,
POPE,POPE,C23,35,0.2172,0.2168,0.2175,
POPE,POPE,C24,38,0.2214,0.2193,0.2235,
POPE,POPE,C25,41,0.2333,0.2341,0.2325,
POPE,POPE,C26,44,0.2137,0.2144,0.2129,
POPE,POPE,C27,47,0.2040,0.2036,0.2044,
POPE,POPE,C28,50,0.1206,0.1193,0.1220,
POPE,POPE,C29,53,0.0748,0.0748,,
POPE,POPE,C210,55,0.0600,0.0600,,
(...)
POPE,POPE,C218,78,0.0324,0.0329,0.0321,0.0322
(...)
POPC,POPC,C22,32,0.1094,0.0953,0.1235,
POPC,POPC,C32,41,0.2325,0.2405,0.2245,
POPC,POPC,C23,44,0.2228,0.2232,0.2224,
POPC,POPC,C24,47,0.2228,0.2228,0.2228,
POPC,POPC,C25,50,0.2356,0.2366,0.2347,
POPC,POPC,C26,53,0.2150,0.2165,0.2134,
POPC,POPC,C27,56,0.2058,0.2055,0.2062,
POPC,POPC,C28,59,0.1220,0.1239,0.1201,
POPC,POPC,C29,62,0.0712,0.0712,,
POPC,POPC,C210,64,0.0648,0.0648,,
(...)
POPC,POPC,C218,87,0.0344,0.0351,0.0340,0.0341
(...)
POPG,POPG,C22,25,0.1081,0.1024,0.1138,
POPG,POPG,C32,34,0.2028,0.2131,0.1926,
POPG,POPG,C23,37,0.2107,0.2172,0.2041,
POPG,POPG,C24,40,0.2170,0.2189,0.2151,
POPG,POPG,C25,43,0.2367,0.2341,0.2392,
POPG,POPG,C26,46,0.2193,0.2191,0.2195,
POPG,POPG,C27,49,0.2101,0.2094,0.2108,
POPG,POPG,C28,52,0.1233,0.1181,0.1286,
POPG,POPG,C29,55,0.0806,0.0806,,
POPG,POPG,C210,57,0.0597,0.0597,,
(...)
POPG,POPG,C218,80,0.0351,0.0361,0.0345,0.0346

Each line corresponds to one heavy atom type (or one bond type in the case of coarse-grained systems). The first column corresponds to the molecule name, the second to the residue name, the third to the atom name, the fourth to the relative index of the atom type, the fifth to the order parameter calculated for this atom type, and the following columns correspond to the order parameters of the bonds the heavy atom type is involved in.

Table format

The Table format presents the results of the analysis in a clear, human-readable format.

To write the results of the analysis in Table format, add output_tab: PATH_TO_OUTPUT to your input YAML file.

Here is an excerpt from a Table file:

# Order parameters calculated with 'gorder v0.2.0' using structure file 'system.tpr' and trajectory file 'md.xtc'.

Molecule type POPE
           TOTAL   |   H #1   |   H #2   |   H #3   |
C22        0.1098  |  0.0946  |  0.1250  |          |
C32        0.2341  |  0.2438  |  0.2244  |          |
C23        0.2172  |  0.2168  |  0.2175  |          |
C24        0.2214  |  0.2193  |  0.2235  |          |
(...)

Molecule type POPC
           TOTAL   |   H #1   |   H #2   |   H #3   |
C22        0.1094  |  0.0953  |  0.1235  |          |
C32        0.2325  |  0.2405  |  0.2245  |          |
C23        0.2228  |  0.2232  |  0.2224  |          |
C24        0.2228  |  0.2228  |  0.2228  |          |
(...)

Molecule type POPG
           TOTAL   |   H #1   |   H #2   |   H #3   |
C22        0.1081  |  0.1024  |  0.1138  |          |
C32        0.2028  |  0.2131  |  0.1926  |          |
C23        0.2107  |  0.2172  |  0.2041  |          |
C24        0.2170  |  0.2189  |  0.2151  |          |
(...)

Each line corresponds to one heavy atom type (or one bond type in the case of coarse-grained systems). The first column corresponds to the name of the atom type, the second to the order parameter calculated for this atom type, and the following columns correspond to the order parameters of the bond the heavy atom type is involved in.

XVG format

XVG format may be useful if you want to quickly visualize the results as plots using xmgrace.

To write the results of the analysis in XVG format, add output_xvg: FILE_PATTERN to your input YAML file. One XVG file will be created for each molecule. The name of the molecule will be appended to the end of the file stem. For example, if your file pattern is order.xvg and the membrane contains POPC, POPE, and POPG lipid molecules, three XVG files will be generated: order_POPC.xvg, order_POPE.xvg, and order_POPG.xvg.

Here is an excerpt from an XVG file:

# Order parameters calculated with 'gorder v0.2.0' using structure file 'system.tpr' and trajectory file 'md.xtc'.
@    title "Atomistic order parameters for molecule type POPC"
@    xaxis label "Atom"
@    yaxis label "-Sch"
@    s0 legend "Full membrane"
@TYPE xy
# Atom C22:
1      0.1094 
# Atom C32:
2      0.2325 
# Atom C23:
3      0.2228 
# Atom C24:
4      0.2228 
# Atom C25:
5      0.2356 
# Atom C26:
6      0.2150 
# Atom C27:
7      0.2058 
# Atom C28:
8      0.1220 
(...)

XVG files provide the least complete information. The first column is the index of the analyzed heavy atom type (starting from 1), and the second column is the order parameter of the atom type.

Example of an input YAML file with all supported outputs requested

structure: system.tpr
trajectory: md.xtc
analysis_type: !AAOrder
  heavy_atoms: "@membrane and name r'C3([2-9]|1[0-6])|C2([2-9]|1[0-8])'"
  hydrogens: "@membrane and element name hydrogen"
output_yaml: order.yaml    # for yaml format, you can use either `output` or `output_yaml`
output_tab: order.tab
output_xvg: order.xvg
output_csv: order.csv

Order parameters for individual leaflets

gorder can calculate order parameters for the entire membrane as well as for the individual leaflets. To do this, you need to specify a method for classifying lipids into membrane leaflets. gorder assigns lipids to membrane leaflets independently for each frame, making it suitable even for the analysis of membranes where lipids flip-flop between leaflets.

There are three leaflet classification methods available in gorder: global, local, and individual.

Global method for leaflet classification

Fast and reliable. Recommended, especially good for disrupted membranes.

In this method, lipid molecules are assigned to membrane leaflets based on the position of their 'head identifier' relative to the global membrane center of geometry. The 'head identifier' is a single atom representing the head of the lipid. If the 'head identifier' is located "above" the membrane center, the lipid is assigned to the upper leaflet; if it is located "below", it is assigned to the lower leaflet.

To use this method, you must specify the 'head identifier' atoms and all atoms that form the membrane. GSL is used to define these selections.

leaflets: !Global
  membrane: "@membrane" 
  heads: "name P"

Here, we use autodetected membrane atoms to calculate the membrane center and select atoms named 'P' (phosphorus atoms of lipids) as head identifiers. Each analyzed lipid must have exactly one head identifier atom; otherwise, an error will occur.

Local method for leaflet classification

Very slow but reliable. Useful for some curved systems.

In this method, lipid molecules are assigned to membrane leaflets based on the position of their 'head identifier' relative to the local membrane center of geometry. The local membrane center is calculated using atoms in a cylinder around the 'head identifier'. If the 'head identifier' is located "above" the local center, the lipid is assigned to the upper leaflet; if "below", it is assigned to the lower leaflet.

For this method, you need to specify a selection of head identifiers, all atoms forming the membrane, and the radius of the cylinder used to define the local membrane.

leaflets: !Local
  membrane: "@membrane"
  heads: "name P"
  radius: 2.5

Autodetected membrane atoms will be used to calculate the membrane center. Only atoms within a cylinder of radius 2.5 nm (with infinite height) centered on the 'head identifier' and oriented along the membrane normal will be used for the local center calculation. The atoms named 'P' (phosphorus atoms of lipids) are used as 'head identifiers'.

Individual method for leaflet classification

Fast but less reliable. Suitable for large membranes.

In this method, lipid molecules are assigned to membrane leaflets based on the position of their 'head identifier' relative to their 'tail ends'. 'Tail ends' refer to the last heavy atoms or beads of the lipid tails. Each lipid molecule may have multiple 'tail ends', but only one 'head identifier'. If the 'head identifier' is located "above" the 'tail ends', the lipid is assigned to the upper leaflet; if it is located "below", it is assigned to the lower leaflet.

To use this method, you must specify selections for the 'head identifiers' and the 'tail ends'.

leaflets: !Individual
  heads: "name P"
  methyls: "name C218 C316"

In this example, atoms named 'P' (phosphorus atoms of lipids) are used as head identifiers, and 'C218' or 'C316' atoms (the last carbons of oleoyl and palmitoyl chains) are used as tail ends.

Example output YAML file

When the leaflet classification method is specified, gorder will calculate order parameters for both the entire membrane and for the individual leaflets. Here is an excerpt from an output YAML file containing results for the individual membrane leaflets:

# Order parameters calculated with 'gorder v0.2.0' using structure file 'system.tpr' and trajectory file 'md.xtc'.
- molecule: POPE
  order:
    POPE C22 (23):
      total: 0.1098
      upper: 0.1101
      lower: 0.1095
      bonds:
        POPE H2R (24):
          total: 0.0946
          upper: 0.0938
          lower: 0.0955
        POPE H2S (25):
          total: 0.125
          upper: 0.1265
          lower: 0.1235
    POPE C32 (32):
      total: 0.2341
      upper: 0.2309
      lower: 0.2372
      bonds:
        POPE H2X (33):
          total: 0.2438
          upper: 0.2409
          lower: 0.2467
        POPE H2Y (34):
          total: 0.2244
          upper: 0.2209
          lower: 0.2278
  #(...)
- molecule: POPC
  order:
    POPC C22 (32):
      total: 0.1094
      upper: 0.1065
      lower: 0.1124
      bonds:
        POPC H2R (33):
          total: 0.0953
          upper: 0.0928
          lower: 0.0978
        POPC H2S (34):
          total: 0.1235
          upper: 0.1201
          lower: 0.1269
    POPC C32 (41):
      total: 0.2325
      upper: 0.2363
      lower: 0.2287
      bonds:
        POPC H2X (42):
          total: 0.2405
          upper: 0.2455
          lower: 0.2354
        POPC H2Y (43):
          total: 0.2245
          upper: 0.2272
          lower: 0.2219
  #(...)
- molecule: POPG
  order:
    POPG C22 (25):
      total: 0.1081
      upper: 0.0975
      lower: 0.1203
      bonds:
        POPG H2R (26):
          total: 0.1024
          upper: 0.0934
          lower: 0.1128
        POPG H2S (27):
          total: 0.1138
          upper: 0.1016
          lower: 0.1278
    POPG C32 (34):
      total: 0.2028
      upper: 0.2174
      lower: 0.1862
      bonds:
        POPG H2X (35):
          total: 0.2131
          upper: 0.2315
          lower: 0.192
        POPG H2Y (36):
          total: 0.1926
          upper: 0.2033
          lower: 0.1803
  #(...)

Maps of order parameters

Basic analysis

gorder can create maps of order parameters ("ordermaps") for the individual bond types and atom types. To construct an ordermap, the membrane plane is mapped onto a grid with bins of a specified size. Order parameters calculated for individual bonds are projected into the respective bins based on the position of the bond center. The average order parameter is then computed for each bin in the grid.

To enable the calculation of ordermaps, add the following section to your input YAML file:

ordermaps:
  output_directory: "ordermaps"

gorder will generate one ordermap per bond type and save it in the ordermaps directory. If you calculate atomistic order parameters, ordermaps for individual heavy atom types will also be created as averages of the order parameters of heavy atom - hydrogen bonds. If the calculation of order parameters for the individual membrane leaflets is requested (see Order parameters for individual leaflets), additional ordermaps for the upper and lower leaflets will be constructed.

If the output directory (here ordermaps) does not exist, it will be created. If it already exists, it will be backed up unless the overwrite option is enabled.

By default:

The ordermap is placed in the plane perpendicular to the membrane normal (usually the xy plane since the membrane normal is typically along the z-axis).
The dimensions of the ordermaps are derived from the simulation box dimensions in the input structure file.
The bin size is set to 0.1×0.1 nm.
At least one sample per bin is required to calculate the order parameter for that bin.

You can manually override these default settings, see below.

Options for ordermaps

Plane of the ordermaps

To specify the plane for constructing the ordermaps, use the plane option:

ordermaps:
  output_directory: "ordermaps"
  plane: xz

The ordermaps will be created in the xz plane. For the xz plane, the first dimension corresponds to the x-axis and the second to the z-axis. For the yz plane, the first dimension corresponds to the z-axis and the second to the y-axis.

Dimensions of the ordermaps

To define the size (span) of the ordermaps, use the dim option:

ordermaps:
  output_directory: "ordermaps"
  dim:
    - !Manual { start: 5.0, end: 10.0 }
    - !Auto

In this example, the ordermaps are placed in the xy plane. They will cover an area spanning 5–10 nm along the x-axis, while the size along the y-axis will be automatically set based on the simulation box size.

Size of the bins

To set the bin size (granularity of the ordermaps), use the bin_size option:

ordermaps:
  output_directory: "ordermaps"
  bin_size: [0.05, 0.2]

Here, the bins are 0.05 nm wide along the x-axis and 0.2 nm wide along the y-axis.

Required number of samples

To specify the required number of samples per bin, use the min_samples option:

ordermaps:
  output_directory: "ordermaps"
  min_samples: 50

For this example, at least 50 samples must be collected in a bin to calculate the order parameter for it. If the collected sample count is less than 50, the order parameter for that bin will be recorded as NaN.

Examples of ordermaps

Ordermaps are saved in a custom two-dimensional gorder format:

# Map of average order parameters calculated for the atom type POPC-C22-32.
# Calculated with 'gorder v0.2.0'.
@ xlabel x-dimension [nm]
@ ylabel y-dimension [nm]
@ zlabel order parameter ($-S_{CH}$)
@ zrange -1 1 0.2
$ type colorbar
$ colormap seismic_r
0.0000 0.0000 0.1025
0.0000 0.1000 0.0610
0.0000 0.2000 0.0922
0.0000 0.3000 0.0798
0.0000 0.4000 0.1215
0.0000 0.5000 0.1275
0.0000 0.6000 0.1120
0.0000 0.7000 0.0732
0.0000 0.8000 0.0723
0.0000 0.9000 0.1040
0.0000 1.0000 0.0711
(...)
0.1000 0.0000 0.0825
0.1000 0.1000 0.0918
0.1000 0.2000 0.1077
0.1000 0.3000 0.1022
0.1000 0.4000 0.1100
0.1000 0.5000 0.1124
0.1000 0.6000 0.0917
0.1000 0.7000 0.1038
0.1000 0.8000 0.1128
0.1000 0.9000 0.1064
0.1000 1.0000 0.1046
(...)
0.2000 0.0000 0.0603
0.2000 0.1000 0.0873
0.2000 0.2000 0.1058
0.2000 0.3000 0.1097
0.2000 0.4000 0.1014
0.2000 0.5000 0.1147
0.2000 0.6000 0.1199
0.2000 0.7000 0.1315
0.2000 0.8000 0.1123
0.2000 0.9000 0.1012
0.2000 1.0000 0.1207
(...)

This example ordermap corresponds to the atom type POPC-C22-32 (C22 of POPC). Lines beginning with # are comments, while lines starting with @ and $ define plotting parameters. Other lines contain the calculated order parameters. The first column represents the first coordinate of the bin (here x-dimension), the second column represents the second coordinate (here y-dimension), and the third column contains the average order parameter for the bin. The data is stored in row-major order. Note that the order parameter may be NaN; ensure your plotting software can handle such cases.

When visualized, the ordermap for this atom type might look like this:

Example of an ordermap

Ordermaps are especially useful for analyzing simulations involving transmembrane proteins. For such analyses, ensure the protein is centered and, ideally, RMSD-fitted in the trajectory.

Example of an ordermap with transmembrane protein

In this map of coarse-grained order for a particular lipid bond, the bright white region in the center represents the area occupied by the transmembrane protein. In this zone, no lipids (or an insufficient number to meet the min_samples threshold of 50) were detected. The membrane is disrupted in the vicinity of the protein.

Analyzing a part of the trajectory

By default, gorder analyzes the entire provided trajectory, meaning it processes all frames from the start to the end of the trajectory.

If you want to analyze only a specific part of the trajectory, you can use the begin, end, and step options to control the time range and frequency of analysis.

begin: Specifies the time (in ps) from which the trajectory will start being read.
end: Specifies the time (in ps) at which the trajectory analysis will stop.
step: Specifies the step size for reading trajectory frames (i.e., how often a frame will be analyzed).

Example: Calculating coarse-grained order parameters from a specific part of a trajectory

structure: system.tpr
trajectory: md.xtc
analysis_type: !CGOrder
  beads: "@membrane"
begin: 100000.0
end: 300000.0
step: 5
output: order.yaml

In this example, the analysis starts at time 100,000 ps (100 ns) and ends at time 300,000 ps (300 ns). Every 5th frame of the trajectory will be analyzed.

Multithreaded analysis

The trajectory analysis can be performed using multiple threads. Utilizing more threads is generally much faster than using just one, although the speed gain also depends on the read performance of your storage device.

By default, gorder uses only 1 thread to perform the analysis. If you want to speed up the process by using multiple threads, you can specify the number of threads in your input YAML file:

n_threads: 4

In the example above, the analysis will be performed using 4 threads.

gorder guarantees binary identical output, regardless of the number of threads used.

Other options

Membrane normal

By default, gorder uses the z-axis as the membrane normal, assuming that the membrane is built in the xy-plane. If your membrane is built in a different plane, you can manually specify the membrane normal in the input YAML file:

membrane_normal: x

In this example, the membrane normal is oriented along the x-axis, meaning that the membrane is built in the yz-plane.

Min samples

By default, gorder calculates order parameters for all bonds for which it collects at least one sample. If you want to increase the number of required samples, you can specify the minimum number of samples manually in the input YAML file:

min_samples: 10000

In this example, at least 10,000 samples are required for the bond to have a valid order parameter. If the number of samples is lower than this value, the order parameter will be NaN.

Overwrite

By default, gorder does not overwrite output files. Before writing an output file (or creating an output directory for the ordermaps), any existing file (or directory) with the same name is backed up. You can change this behavior by specifying the overwrite option in the input YAML file:

overwrite: true

In this example, no backups will be created, and files will be overwritten directly.

You can also specify this option on the command line when calling gorder:

$ gorder INPUT_YAML --overwrite

Silent

By default, gorder writes a lot of information about the progress of the analysis to the terminal (standard output). If you prefer not to see any of this information, you can specify it in the input YAML file:

silent: true

In this example, gorder will not write anything to standard output. Output files will still be generated, and the analysis will proceed as normal. If an error occurs, information about the error will be written to standard error output.

You can also specify this option on the command line when calling gorder:

$ gorder INPUT_YAML --silent

Using `gorder` as a Rust crate

gorder can also be used as a Rust crate, allowing you to call it directly from your Rust code.

To use gorder in your project, first add it as a dependency:

$ cargo add gorder

Next, import the crate into your Rust code:

#![allow(unused)]
fn main() {
use gorder::prelude::*;
}

Once imported, you can access all the options and functionality described in this manual.

For instance, the following input yaml file:

structure: system.tpr
trajectory: md.xtc
analysis_type: !CGOrder
  beads: "@membrane"
output: order.yaml

can also by written as the following Rust code:

#![allow(unused)]
fn main() {
Analysis::new()
    .structure("system.tpr")
    .trajectory("md.xtc")
    .analysis_type(AnalysisType::cgorder("@membrane"))
    .output("order.yaml")
    .build()
    .unwrap()
}

For some more examples on how to specify the analysis options and how to perform the analysis in Rust, refer to the corresponding docs.rs.

gorder manual