Open orbeckst opened 11 months ago
From discord #users:
I'm trying to create a group using select_atoms and the cylayer keyword. However, reading the doc I've noticed that the selection argument has to be a group. Is there a way to use the cylayer keyword centered around a specific xyz coordinate in my box?
For high performance you might actually not use selections but do explicit calculations with the positions. Importantly, we have to take periodic boundaries into account with using MDAnalysis.lib.distances.minimize_vectors()
:
import numpy as np
import MDAnalysis as mda
from MDAnalysis.lib import distances
from MDAnalysisTests import datafiles as data
# make a universe to work with, this should really be YOUR data
u = mda.Universe(data.TPR, data.XTC)
# cyzone geometric parameters
# center should be located INSIDE the simulation box
center = np.array([-10.2, 0, 3.6])
zmax, zmin = 25, 5
rInner, rOuter = 5, 10
ag = u.select_atoms("resname SOL and name OW") # just oxygen for simplicity
# calculate distance from center WITH PBC
x = distances.minimize_vectors(ag.positions - center, ag.dimensions)
# make geometric selections
within_z_limits = (zmin < x[:, 2]) & (x[:, 2] < zmax)
radii = np.linalg.norm(x[:, :2], axis=1)
within_r = (rInner < radii) & (radii < rOuter)
within_cyzone = within_r & within_z_limits
# get atoms that fulfill geometric criteria
layer = ag[within_cyzone]
EDIT: **Unfortunately, this does not work as expected because the selections are still made inside the same universe.*** I tried out different positions and didn't manage to shift the selection.
~Yes, you're right, the selection argument in cylayer innerRadius externalRadius zMax zMin selection
has to be something that selects atoms. But it is more general: it can be any AtomGroup
so we can use the group
keyword to insert an arbitrary AtomGroup
as the reference point via~
u.select_atoms("cylayer {innerRadius} {externalRadius} {zMax} {zMin} group myGroup", myGroup=some_atomgroup)
~We can craft an AtomGroup that contains a fake atom with a fixed coordinate (using the very flexible Universe.empty()
:~
import MDAnalysis as mda
from MDAnalysisTests import datafiles as data
# make a universe to work with, this should really be YOUR data
u = mda.Universe(data.TPR, data.XTC)
# create a fake universe with a single atom as a fixed reference
fixed = mda.Universe.empty(n_atoms=1, trajectory=True)
# set the coordinates to what you need (this sets all atoms to the same value, but we only have one)
fixed.atoms.positions = [-10.2, 0, 3.6]
# select with the fixed reference group to anchor cyzone
# (you probably want updating=True if the contents of the selection should change
# along a trajectory)
layer = u.select_atoms("resname SOL and cylayer 5 10 25 5 group REFERENCE", REFERENCE=fixed.atoms, updating=True)
~Now layer
is your selected AtomGroup. Note that you could just provide reference.atoms as part of the selection and that these groups do not have to come from the same Universe! This gives you a lot of flexibility.~
Do double check (eg by writing out frames and visually checking) that this works as expected
I would want to update the FAQ with some of the issues that I repeatedly see new users trip up. I am listing a few points here and use the issue as a notepad.
Basics
The Introduction in the Quickstart is already laying out most of the key things to be aware of. It should be copied into the API docs and as a separate intro chapter for the rest of the user guide — it's really important to tell people up-front what the basic philosophy is.
Only one frame is loaded into memory. (in principle addresses in Why do the atom positions change over trajectories?)
MDAnalysis keeps track of the trajectory frame (keeps state). Think of it as a "cursor" that points to a frame on the trajectory on disk. Only this frame is "active".
Generally, figuring out when you need to make a copy is confusing (e.g. when using
Timestep
in the canonicalfor ts in u.trajectory
).Groups (Containers)
"fragments" are groups of particles that are all connected by bonds, and no bonded particles are excluded (like chemical molecules); fragments are not part of the Segment/Residue/Atom hierarchy (a fragment can span multiple residues or even segments)
"segments" are arbitrary but mutually exclusive collections of "Residue" instances.
AtomGroups are arbitrary collections of Atom instances.
attributes are "tags" for Atoms and containers (AtomGroups, Residues, ResidueGroups, Segments, SegmentGroups) that describe the Atom or container. attributes can be used for selections. Arbitrary attributes can be added.
Selecting
updating=True
) selections and really the need to be aware of it!Trajectory manipulation
"unwrap" means to "make molecules whole across periodic unit cell boundaries". "nojump" means to make trajectories continuous across periodic boundaries.
In order to do RMSD-superposition, you typically need to make a copy of the coordinates.