Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

update FAQ / basic philosophy #341

Open
orbeckst opened this issue Nov 7, 2023 · 1 comment
Open

update FAQ / basic philosophy #341

orbeckst opened this issue Nov 7, 2023 · 1 comment

Comments

@orbeckst
Copy link
Member

orbeckst commented Nov 7, 2023

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?)

    A key idea is that only one frame of the trajectory is read from disk into memory. All objects from the universe then work with the data of the current frame. u.atoms.positions are the positions at the current frame. When you do u.trajectory[-1] you will move to the end of the trajectory. Now u.atoms.positions will contain a new set of coordinates.

  • 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 canonical for 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

  • Difference between static and dynamic (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.

@orbeckst
Copy link
Member Author

orbeckst commented Nov 9, 2023

From discord #users:

Q:How to make spatial selections relative to a FIXED point?

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?

A: Explicit selection with PBC

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]

NOT WORKING: Use a fake AtomGroup with a fixed dummy atom

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant