How is the cylinder defined in cyzone ?

[peabreu]

Hello there,
I have been using MDAnalysis to calculate the number of molecules above and below a molecule (graphene oxide flake) and I have been using cyzone to select these molecules. But I can not understand how the cylinder axis is defined. According to what I found in the manual ( https://docs.mdanalysis.org/2.9.0/documentation_pages/selections.html ):
cyzone externalRadius zMax zMin selection
"selects all atoms within a cylindric zone centered in the center of geometry (COG) of a given selection,
[...]
Changed in version 0.10.0: keywords cyzone and cylayer now take zMax and zMin to be relative to the COG of selection, instead of absolute z-values in the box."
The manual never discusses the orientation of the created cylinder. Maybe I am missing something or misunderstanding the selection method. If the COG has coordinates (x_g, y_g, z_g) does this means that the selection is just for entities that have
z < z_g + zMax and z > z_g + Zmin and are within a distance less than externalRadius from the COG ?
I have tried to read the relevant code but I do not quite understand how it works.
As my molecule is rotating in my simulation I just want o be sure if I am using the correct selection to select molecules above and below my molecule ( the flake shape can be approximated very well by a circle ) but I am not sure if cyzone is the answer to my problem.
Thanks.
Paulo

[thisisshivamtiwari]

Your reading is almost right, with one important difference: the cylinder axis is always the box z-axis (lab frame). It is not tied to the orientation of your molecule.

From CylindricalSelection._apply:

vecs = group.positions - sel.center_of_geometry()            # L477 (minimum-imaged if periodic)
mask = (vecs[:, 2] > self.zmin) & (vecs[:, 2] < self.zmax)   # L516
radii = vecs[:, 0] ** 2 + vecs[:, 1] ** 2                    # L522
mask = radii < self.exRadius**2                              # L523

So with COG = (x_g, y_g, z_g), cyzone R zMax zMin sel selects atoms with

• z_g + zMin < z < z_g + zMax (strict), and
• (x - x_g)² + (y - y_g)² < R². This is the distance in the xy-plane from the axis, not the 3D distance from the COG.

A quick check on a toy universe (MDAnalysis 2.11) confirms it: an atom 5.9 Å to the side and 3 Å above the COG (3D distance 6.6 Å) is selected by cyzone 6 8 -8, an atom exactly at dz = zMax is not, and tilting the reference molecule by 45° does not change the selection.

That means cyzone only does what you want while the flake stays flat in the xy-plane. For a rotating flake you can build the cylinder along the flake normal yourself, frame by frame:

import numpy as np
from MDAnalysis.lib.distances import minimize_vectors

def flake_cylinder(flake, others, radius, zmax, zmin):
    n = flake.principal_axes()[0]        # largest moment of inertia -> plate normal
    vecs = others.positions - flake.center_of_geometry()
    if others.dimensions is not None:
        vecs = minimize_vectors(vecs, others.dimensions)
    h = vecs @ n                                         # signed height along the normal
    r = np.linalg.norm(vecs - np.outer(h, n), axis=1)   # distance from the axis
    keep = (h > zmin) & (h < zmax) & (r < radius)
    return others[keep], h[keep]

flake = u.select_atoms("resname GO")
water = u.select_atoms("resname SOL and name OW")
for ts in u.trajectory:
    inside, h = flake_cylinder(flake, water, 15, 8, -8)
    n_above, n_below = (h > 0).sum(), (h < 0).sum()

Two notes:

• principal_axes() sorts by eigenvalue (largest first; see core/topologyattrs.py). For a flat disc the largest moment is about the normal, so [0] is the normal. Its sign is arbitrary. If "above" and "below" need to stay consistent over time, flip n to agree with the previous frame (if n @ n_prev < 0: n = -n).
• principal_axes() uses masses, so the topology needs them. If the flake is broken across the periodic boundary, use principal_axes(wrap=True) or make the flake whole first.

[peabreu]

Thanks a lot for the detailed explanation. The code was also very helpful.