77 views

Skip to first unread message

Jul 6, 2022, 6:03:04 AMJul 6

to MDnalysis discussion

Hi,

How should we calculate the angle having three points by MDnalysis?

Best regards,

Rose

Jul 6, 2022, 6:09:03 AMJul 6

to mdnalysis-...@googlegroups.com

Dear Rose,

pls google, such as: mdanalysis angle

Second hit:

https://becksteinlab.physics.asu.edu/pages/courses/2013/SimBioNano/15/MDAnalysis/trajectories.html

See Exercises 4

Bests, tamas

pls google, such as: mdanalysis angle

Second hit:

https://becksteinlab.physics.asu.edu/pages/courses/2013/SimBioNano/15/MDAnalysis/trajectories.html

See Exercises 4

Bests, tamas

--

You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.

To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/7a6cf4b0-0bf3-4048-bc1d-cc7c0300fbd0n%40googlegroups.com.

Jul 6, 2022, 6:29:09 AMJul 6

to MDnalysis discussion

Hi Rose,

If you already have your 3 points, you can use calc_angles: https://docs.mdanalysis.org/stable/documentation_pages/lib/distances.html#MDAnalysis.lib.distances.calc_angles

With this function, you provide the 3 points as numpy arrays. These arrays can be either of shape (3,) if you only have one angle to compute, or of shape (3, N) if you have N angles.

If you work with periodic systems, you should provide the box, otherwise, you may have wrong results if one of the 3 atoms goes through the box.

Cheers,

Jonathan

Jul 26, 2022, 4:54:48 AMJul 26

to MDnalysis discussion

Hi,

Thank you for your answer. I have tried to do the exercise 4 for my system:

import numpy as np

import MDAnalysis as mda

from numpy.linalg import norm

u = mda.Universe("plumed-eq.tpr", "plumed-eq.gro")

npa = u.select_atoms("resname ZnO")

npa_com = npa.center_of_mass()

oxygen_coordinate = u.select_atoms("name OW and resid 3").positions

hydrogen_coordinate = u.select_atoms("name HW1 and resid 3").positions

BA = npa_com - oxygen_coordinate

BC = hydrogen_coordinate - oxygen_coordinate

theta = np.arccos(np.dot(BA, BC)/(norm(BA)*norm(BC)))

tetha

import MDAnalysis as mda

from numpy.linalg import norm

u = mda.Universe("plumed-eq.tpr", "plumed-eq.gro")

npa = u.select_atoms("resname ZnO")

npa_com = npa.center_of_mass()

oxygen_coordinate = u.select_atoms("name OW and resid 3").positions

hydrogen_coordinate = u.select_atoms("name HW1 and resid 3").positions

BA = npa_com - oxygen_coordinate

BC = hydrogen_coordinate - oxygen_coordinate

theta = np.arccos(np.dot(BA, BC)/(norm(BA)*norm(BC)))

tetha

--------

But I got this error;

ValueError Traceback (most recent call last)
/tmp/ipykernel_178608/45643382.py in <module>
12 BA = npa_com - oxygen_coordinate
13 BC = hydrogen_coordinate - oxygen_coordinate
---> 14 theta = np.arccos(np.dot(BA, BC)/(norm(BA)*norm(BC)))
15 tetha
16
<__array_function__ internals> in dot(*args, **kwargs)
ValueError: shapes (1,3) and (1,3) not aligned: 3 (dim 1) != 1 (dim 0)

How can I define the center of mass to be aligned with two other coordination?

Best,

Rose

Jul 27, 2022, 7:55:56 PMJul 27

to mdnalysis-discussion

Hi Rose,

The shape of the position difference is generally (N, 3) where N is the size of the AtomGroup. In your case the shape is (1, 3) because you only have one atom each. If you want to use np.dot(a, b) as a vector inner product, the two arguments a and b must have shape (3,). The easiest way to do this in your case is

BA = npa_com - oxygen_coordinate

BA = BA[0]

BC = hydrogen_coordinate - oxygen_coordinate

BC = BC[0]

or in one line each:

BA = (npa_com - oxygen_coordinate)[0]

BC = (hydrogen_coordinate - oxygen_coordinate)[0]

where you just select the first element at index 0 from the [[x,y,z]] array.

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/4087b58c-56ca-49f9-a9c8-f93d3e09000an%40googlegroups.com.

Message has been deleted

Jul 28, 2022, 5:16:47 AMJul 28

to MDnalysis discussion

Hi Oliver,

Thank you so much for your help.

What
if I want to calculate over the trajectory? and this time don't limit
to only one water but the whole waters in each trajectory. For each
trajectory, I want to calculate the oxygen_coordinat and
hydrogen_coordinate per each water molecule. To be more specific, this
is what I finally want to calculate:

1- I want to calculate the histogram of (over the trajectories)angle consisting 3 points : first one is
the center of mass of nanoparticle(np) second one is the oxygen atom of
water molecule and third is H atom of water.

So in each
trajectory, the center of mass of nanoparticle (np_COM) is one point but
there are 68545(residue 3 to residue 68547) water molecules which I
need to have the hydrogen_ coordinate and oxygen_coordinate per each
water molecule. So I will have 68547 theta per each trajectory and at
the end I want to have the distribution of angles for whole simulation
time.

So I did this but I get the same angle:

------------------------------

import numpy as np

import MDAnalysis as mda

from numpy.linalg import norm

import MDAnalysis as mda

from numpy.linalg import norm

traj = mda.Universe("plumed-eq.xtc")

u = mda.Universe("plumed-eq.tpr", "plumed-eq.gro")

npa = u.select_atoms("resname ZnO")

npa_com = npa.center_of_mass()

npa_com = npa.center_of_mass()

oxygen_coordinate = u.select_atoms("name OW and resid 3-68547").positions

hydrogen_coordinate = u.select_atoms("name HW1 and resid 3-68547").positions

hydrogen_coordinate = u.select_atoms("name HW1 and resid 3-68547").positions

BA = (npa_com - oxygen_coordinate)[0]

BC = (hydrogen_coordinate - oxygen_coordinate)[0]

theta = np.arccos(np.dot(BA, BC)/(norm(BA)*norm(BC)))

theta

np.rad2deg(theta)

for ts in traj.trajectory:

print (np.rad2deg(theta))

theta

np.rad2deg(theta)

for ts in traj.trajectory:

print (np.rad2deg(theta))

----------------------

94.15836711676553

94.15836711676553

94.15836711676553

94.15836711676553

.

.

...

Best,

Rose

Jul 28, 2022, 7:56:46 PMJul 28

to mdnalysis-discussion

Hi Rose,

You need to **do all the calculations *inside* the loop over frames**:

import numpy as np

import MDAnalysis as mda

from numpy.linalg import norm

import MDAnalysis as mda

from numpy.linalg import norm

u = mda.Universe(TPR, "plumed-eq.xtc”) # load a TPR file for the topology!!

# make selections (only needs to be done once!)

npa = u.select_atoms("resname ZnO”)

water = u.select_atoms("resid 3-68547”) # or “resname SOL” if water is standard solvent in GROMACS

ow = water.select_atoms("name OW")

hw1 = water.select_atoms("name HW1")

# create an empty histogram

NBINS = 91

ANGLE_RANGE = (0, 180)

h, edges = np.histogram([0], bins=NBINS, range=ANGLE_RANGE, density=True)

h *= 0 # reset counts

for ts in u.trajectory:

npa_com = npa.center_of_mass()

BA = npa_com - ow.positions

BC = hw1.positions - ow.positions

theta =np.rad2deg( np.arccos(np.einsum('...i,...i', BA, BC)/(norm(BA, axis=1)*norm(BC, axis=1))) )

print(f"frame={ts.frame} theta={theta}º")

# add to the histogram (or rather, add to the probability density)

tmp_h, _ = np.histogram(theta, bins=NBINS, range=ANGLE_RANGE, density=True)

h += tmp_h

# normalization of the density (average per frame)

h /= u.trajectory.n_frames

# now plot your histogram

import matplotlib.pyplot as plt

midpoints = 0.5*(edges[:-1] + edges[1:])

plt.step(midpoints, h, where="mid")

plt.xlabel("angle (degree)")

plt.ylabel("probability density")

I added code to build up the histogram. We first initialize and empty histogram and then histogram the data for each frame and add it to the histogram. (I’m calculating the histogram as a density, which I can still add like counts but I need to remember to divide everything by the total number of frames at the end.)

Also note that I wrote the code so that you're now doing ALL theta calculations via numpy. In order to do the “dot product for all pairs BA[i] @ BC[i]” I used the einsum() function. Similarly, using norm(…, axis=1) calculates the length of each BA[i] and BA[j].

theta is an ARRAY of size N where N == len(ow). If the above is not quite clear to you, run parts interactively and look at the shape of the arrays, e.g.,

ow.positions.shape

BA.shape

norm(BA, axis=1).shape

BAdotBC = np.einsum('...i,...i', BA, BC)

BAdotBC.shape

BAdotBC[0]

BA[0] @ BC[0] # should be the same as BAdotBC[0]

c = np.einsum('...i,...i', BA, BC)/(norm(BA, axis=1)*norm(BC, axis=1))

c.shape

I hope this helps.

Best,

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/76a0b5be-31a7-4b52-978c-525913f688ebn%40googlegroups.com.

--

Oliver Beckstein (he/his/him)

GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

Jul 31, 2022, 6:41:36 PMJul 31

to MDnalysis discussion

Thank you so much for your help. Is this the same for calculating cos theta? or the water dipole vector instead of -OH vector?

Best regards

Jul 31, 2022, 7:37:52 PMJul 31

to mdnalysis-discussion

On Jul 31, 2022, at 3:41 PM, Roja <rose....@gmail.com> wrote:

Thank you so much for your help. Is this the same for calculating cos theta? or the water dipole vector instead of -OH vector?

In principle, yes. The same ideas apply. You need to do the calculations INSIDE the loop over frames. AtomGroups should be constructed OUTSIDE the loop for better performance (but if they would change such as AROUND selections, then make them updating=True).

The actual calculations for cos theta would simply remove the arccos() function call. For the water dipole vector you have calculate the appropriate vector first. Most of these calculations can be written elegantly with numpy (which takes some time figuring out) but you can always start with a slow version that include loops and then speed it up using numpy with broadcasting (and tricks like einsum()).

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/c769c5c9-fce0-4c49-ba82-547ff6e2b3f9n%40googlegroups.com.

Aug 1, 2022, 10:02:30 AMAug 1

to MDnalysis discussion

So I have done some modefications to calculate the cos.angle between dipole vector of water and the other vector(BA) but I get this error:

ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (0,3)->(0,3) (66396,3)->(66396,3)

u = mda.Universe("nvt.tpr", "nvt.xtc")

# make selections

npa = u.select_atoms("resname ZnO")

ow = water.select_atoms("name OW and sphzone 6.0 (resname ZnO)")

hw1 = water.select_atoms("name HW1 and sphzone 6.0 (resname ZnO)")

# create an empty histogram

ANGLE_RANGE = (-1,1)

h_cos_theta, edges = np.histogram([0], bins=NBINS, range=ANGLE_RANGE, density=True)

h_cos_theta *= 0 # reset counts

for ts in u.trajectory:

npa_com = npa.center_of_mass()

BA = npa_com - ow.positions

cos_theta =np.einsum('...i,...i', BA, dipole)/(norm(BA, axis=1)*norm(dipole, axis=1))

print(f"frame={ts.frame} cos_theta={cos_theta}º")

# add to the histogram (or rather, add to the probability density)

h_cos_theta += tmp_h_theta

# normalization of the density (average per frame)

Aug 1, 2022, 9:02:15 PMAug 1

to mdnalysis-...@googlegroups.com

Hi Rose,

1) I don’t know where your error originates. If in doubt, include full error outputs. You might have to dig deeper into how to do the calculations with numpy. Or someone else knows immediately where the problems occur.

2) You can use the sphzone selection but you need to make it an updating AtomGroup with updating=True. Only an UpdatingAtomGroup will properly recalculate the waters inside the sphere at each step. (This will also make your code slower but there’s not much you can do about it if you want to use AtomGroups.)

Oliver

Am 8/1/22 um 07:02 schrieb Roja <rose....@gmail.com>:

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/e650b765-7f2a-4b4f-a105-ce26e6e15be6n%40googlegroups.com.

Aug 2, 2022, 12:27:46 PMAug 2

to mdnalysis-...@googlegroups.com

One quick comment:

Or if the error comes from einsum(), see if you can just change the call to the function — it takes a while to learn einsum() but it’s worth it.

Oliver

> So I have done some modefications to calculate the cos.angle between dipole vector of water and the other vector(BA) but I get this error:

> ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (0,3)->(0,3) (66396,3)->(66396,3)

The shape of the first operand should probably be (1,3) — try not.reshape().
> ValueError: operands could not be broadcast together with remapped shapes [original->remapped]: (0,3)->(0,3) (66396,3)->(66396,3)

Or if the error comes from einsum(), see if you can just change the call to the function — it takes a while to learn einsum() but it’s worth it.

Oliver

Sep 20, 2022, 10:02:19 AM (6 days ago) Sep 20

to MDnalysis discussion

Hi,

I am trying to calculate the distribution of angles between surface atoms of ZnS nanoparticle (type ZnB or type ZnC or type SB or type SC) and OH vector of water. But when I change the hw1 to hw2 in BA(hw1.positions - ow.positions or hw2.positions - ow.positions) I get different answers. Is there any way to specify that the code takes the hw vector closer to the surface (type ZnB or type ZnC or type SB or type SC)?

Best,

Rose

------------------

npa = u.select_atoms("resname ZnS")

water = u.select_atoms("same residue as name HW1 and around 2.38 (type ZnB or type ZnC or type SB or type SC)", periodic = True, updating=True) #Is it also possible to write: water = u.select_atoms("resname SOL and sphzone 6.0 (npa_com)") ?ow = water.select_atoms("name OW", periodic = True, updating=True)

hw1 = water.select_atoms("name HW1", periodic = True, updating=True)

hw2 = water.select_atoms("name HW2", periodic = True, updating=True)

# create an empty histogram

NBINS = 1000

ANGLE_RANGE = (-1,1)

h_cos_alpha_a *= 0 # reset counts

for ts in u.trajectory:

npa_com = npa.center_of_mass(pbc = True)

BA = ow.positions - npa_com

BC = hw1.positions - ow.positions

cos_alpha =np.einsum('...i,...i', BA, BC)/(norm(BA, axis=1)*norm(BC, axis=1))# add to the histogram (or rather, add to the probability density)

h_cos_alpha_a += tmp_h_cos_alpha_a

# normalization of the density (average per frame)

------------

Sep 20, 2022, 8:43:32 PM (6 days ago) Sep 20

to mdnalysis-discussion

Hi Rose,

On Sep 20, 2022, at 7:02 AM, Roja <rose....@gmail.com> wrote:Hi,I am trying to calculate the distribution of angles between surface atoms of ZnS nanoparticle (type ZnB or type ZnC or type SB or type SC) and OH vector of water. But when I change the hw1 to hw2 in BA(hw1.positions - ow.positions or hw2.positions - ow.positions) I get different answers. Is there any way to specify that the code takes the hw vector closer to the surface (type ZnB or type ZnC or type SB or type SC)?

Do you want to find the HW1 or HW2 that is closer to the nano particle and then use that HW* for the vector? If so, you’ll need to calculate the distance of your HW* to the npa atoms and then choose the closer one. You can then create a vector of the correct HW1 and HW2 coordinates

BC = hw.positions - ow.positions

There are different ways to do this. I played around with arrays and came up with the following: NOTE: not really tested, just as an idea:

npa_com = npa.center_of_mass(unwrap=True) # I think you want unwrap=True, see docs

water = u.select_atoms(f"resname SOL and sphzone {npa_com[0]} {npa_com[1]} {npa_com[2]}”)

hw1 = water.select_atoms("name HW1")

hw2 = water.select_atoms("name HW2”)

ow =

# distance array is compute intensive; hopefully, there aren’t too many waters in the group

d1 = mda.lib.distances.distance_array(hw1.positions, protein.positions, box=u.dimensions)

d2 = mda.lib.distances.distance_array(hw2.positions, protein.positions, box=u.dimensions)

# create a boolean array that contains True when HW1 is closer than HW2

hw1_close = d1.min(axis=1) - d2.min(axis=1) < 0

# and the opposite for HW2

hw2_close = ~hw1_close

# now make a group of the closest hydrogens

hw_closest = hw1[hw1_close] + hw2[hw2_close]

# the order has changed so we need to match the OWs

ow_reordered = ow [hw1_close] + ow[hw2_close]

# calculate the vectors

BC = hw_closest,positions - ow_reordered.positions

BA = ow_reordered.positions - npa_com

...

The idea in the code above is to work with arrays and AtomGroups as much as possible and avoid a loop over individual water molecules.

The distance_array() is slow and one could use capped_distance() instead, but then you have to do more book keeping. I assume that the around selection is already expensive and reduced the number of waters sufficiently.

Best,Rose------------------npa = u.select_atoms("resname ZnS")water = u.select_atoms("same residue as name HW1 and around 2.38 (type ZnB or type ZnC or type SB or type SC)", periodic = True, updating=True) #Is it also possible to write: water = u.select_atoms("resname SOL and sphzone 6.0 (npa_com)") ?

If you want to use the distance to the COM then you can use a static selection where you use a f-string to introduce the npa COM but then you make the selection INSIDE the loops over frames. You’ll also need to calculate the npa.center_of_mass() for every time step (as you already do).

I don’t know why 2.38 vs 6 so you need to try it and compare what you get, e.g., by initial printing the number of waters at each ts.

Best,

Oliver

--

You received this message because you are subscribed to the Google Groups "MDnalysis discussion" group.

To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-discus...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/07ee47bc-39d1-44f6-aa32-539b745e04d9n%40googlegroups.com.

Sep 22, 2022, 8:49:01 AM (4 days ago) Sep 22

to MDnalysis discussion

Hi,

Thank you for your answer. I have tested this code and got an error for this part (resname SOL and sphzone {npa_com[0]} {npa_com[1]} {npa_com[2]}). So I changed it to (same residue as name SOL and point {npa_com[0]} {npa_com[1]} {npa_com[2]} 32.7). I guess the result is reasonable. But now the problem is that when I want to calculate dipole vector of water (and these waters are the ones whose OH orientations have already been calculated by the code you sent me in the previous message). So the question is how should I select water molecules so that it takes the SOL molecules which the orientaion of their OH vectors has already been calculated? As you know the waters on the border (32.7 Å) are the problem right now. Hope I have explained it well:)

And a we expect the difference between peaks in these two plots(OH and dipole vector) should have almost 52 degrees difference (104.7/2).

Thank you so much for your help,

Rose

This was an old comment to myself, sorry :)

Sep 22, 2022, 1:01:11 PM (4 days ago) Sep 22

to mdnalysis-discussion

Hi Rose,

Good that you’ve been able to fix my code ;-).

On Sep 22, 2022, at 5:49 AM, Roja <rose....@gmail.com> wrote:

Hi,

Thank you for your answer. I have tested this code and got an error for this part (resname SOL and sphzone {npa_com[0]} {npa_com[1]} {npa_com[2]}). So I changed it to (same residue as name SOL and point {npa_com[0]} {npa_com[1]} {npa_com[2]} 32.7). I guess the result is reasonable.

Better not guess. You’re ultimately responsible for your research. Be sure. Test it. Break it. Don’t trust me.

But now the problem is that when I want to calculate dipole vector of water (and these waters are the ones whose OH orientations have already been calculated by the code you sent me in the previous message). So the question is how should I select water molecules so that it takes the SOL molecules which the orientaion of their OH vectors has already been calculated? As you know the waters on the border (32.7 Å) are the problem right now. Hope I have explained it well:)

You have the atomgroup hw_closest.You can now ask for the residue that these hydrogens belong to:

water_closest = hw_closest.residues.atoms

This is the atomgroup of the closest waters. Now do your dipole calculation with them.

Oliver

To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-discussion/a1775bf3-19cd-45ba-b9f1-26427d64ec25n%40googlegroups.com.

Reply all

Reply to author

Forward

0 new messages