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

Selecting atom indices within a particular distance #1600

Closed
mayankb3192 opened this issue Nov 9, 2020 · 6 comments
Closed

Selecting atom indices within a particular distance #1600

mayankb3192 opened this issue Nov 9, 2020 · 6 comments
Labels

Comments

@mayankb3192
Copy link

Hello,
Is there a way to select atom indices within a particular distance of another atom?
For example to select water molecules within 0.3 nm of protein for a trajectory t, could we use

t.topology.select('water and within 0.3 of protein')
However this brings up issue of 'Cannot use literals as truth'
I have had issues with a distance based atom selection. Any help appreciated.

@sroet
Copy link
Contributor

sroet commented Nov 10, 2020

So there is an inherent issue with selection based on a distance from a trajectory; it is not well-defined.
Say; A trajectory with 2 frames and in frame 1 atoms [1,2,3] are within the cutoff, and in frame 2 atoms [3,4].

Would you then expect this to return "any atom that is within the cutoff at any point in the trajectory" so in this case [1,2,3,4] or "any atom that is within the cutoff for the whole trajectory" so in this case [3].?

What would you like to achieve with this selection?

@mayankb3192
Copy link
Author

Ideally I would like to update it every frame, but even if not can I at least select it for a particular frame? Like VMD gives two options for something like this - update every frame and second keep it constant from the first frame.

@sroet
Copy link
Contributor

sroet commented Nov 12, 2020

mdtraj.compute_neighbors Should work for you then.

Something like this should work (DISCLAIMER: I did not test or run this code)

import mdtraj as md
t = md.load('trajectory.extension')
query = t.topology.select('protein')
haystack = t.topology.select('water')

#This should be a list of 1D np arrays that contain all the water atom indices that are within the cutoff per frame
waters_in_contact_per_frame = md.compute_neighbors(t, cutoff=0.3, query_indices=query, haystack_indices=haystack)

If you want to do more in-depth analysis of these contact, I would like you to point to another package called:
ContactMapExplorer Which is still in development, but reaching maturity ;)

@stale
Copy link

stale bot commented Jun 2, 2021

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

@stale stale bot added the wontfix label Jun 2, 2021
@stale stale bot closed this as completed Jul 8, 2021
@mmfarrugia
Copy link

So there is an inherent issue with selection based on a distance from a trajectory; it is not well-defined. Say; A trajectory with 2 frames and in frame 1 atoms [1,2,3] are within the cutoff, and in frame 2 atoms [3,4].

Would you then expect this to return "any atom that is within the cutoff at any point in the trajectory" so in this case [1,2,3,4] or "any atom that is within the cutoff for the whole trajectory" so in this case [3].?

What would you like to achieve with this selection?

This is from a long time ago, but would anyone know a way to do this but with the first behavior where any atoms (or residues, which would presumably just require byres) which come within a cutoff of a residue at any point in the trajectory are included?

I can see that one way to do this is to use the solution proposed by @sroet and then consolidate the output programmatically, but I have a total of 16 microseconds of MD for which I would do this, so I would really rather not do this so inefficiently if there is another way. I have checked into mdtraj.compute_neighbors but haven't seen anything with quite this behavior, it looks like contact_map is the only one with this functionality.

Ultimately, I am hoping to feed this as an mdtraj selection string for more intensive analysis (in pyemma) so doing it with mdtraj itself or mdtraj selection syntax would be great if possible, but not strictly necessary so I thought that I'd post here as a suggestion for future functionality if there is not currently another method.

Thank you!

@sroet
Copy link
Contributor

sroet commented Jun 21, 2023

I can see that one way to do this is to use the solution proposed by @sroet and then consolidate the output programmatically, but I have a total of 16 microseconds of MD for which I would do this, so I would really rather not do this so inefficiently if there is another way. I have checked into mdtraj.compute_neighbors but haven't seen anything with quite this behavior, it looks like contact_map is the only one with this functionality.

The way we did (something similar) in ContactMapExplorer (which is a mostly stalled project by now) was by running these lines on every frame. If you just want a list of atoms out of it you can probably rewrite the loop to:
(DISCLAIMER the following code has not been tested/run and assumes you can load the whole trajectory in memory)

import mdtraj as md
query = set([1,2,3]) # atom numbers of the residues you want to query
cutoff = 0.45 #cutoff in nm
traj = md.load('traj_file.xtc')
n_frames = len(traj)
atoms = set([])
for frame_number in range(n_frames):
    neighborlist = md.compute_neighborlist(traj, cutoff, frame_number)
    for atom_idx in query:
        neighbor_atoms = set(neighborlist[atom_idx])
        neighbor_atoms -= query #ignore the query atoms, comment out if you want to include them in the output
        atoms |= neighbor_atoms
print(atoms) 

The reason it uses md.compute_neighborlist is because it was significantly faster than the base md.compute_neighbors (If I remeber correctly it was 1 or 2 orders of magnitude faster)

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

No branches or pull requests

3 participants