determining angle

[Roja]

Hi,

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

Best regards,

Rose

[Tamas Hegedus]

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

--
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.

[jona...@barnoud.net]

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

[Roja]

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

--------

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

[Oliver Beckstein]

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.

--

Oliver Beckstein (he/his/him)

orbe...@gmail.com

[Roja]

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

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

oxygen_coordinate = u.select_atoms("name OW 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))

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

94.15836711676553

94.15836711676553

94.15836711676553

94.15836711676553

.

.

...

Best,

Rose

[Oliver Beckstein]

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

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)

email: orbe...@mdanalysis.org
twitter: @orbeckst

GitHub: @orbeckst

MDAnalysis – a NumFOCUS fiscally sponsored project

https://www.mdanalysis.org/

[Roja]

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

[Oliver Beckstein]

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.

[Roja]

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

water = u.select_atoms("resname SOL and sphzone 6.0 (resname ZnO)") #Is it also possible to write: water = u.select_atoms("resname SOL and sphzone 6.0 (npa_com)")   ?
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

NBINS = 1000
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

        dipole = 0.5*(H1.positions - H2.positions) - O.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)

        tmp_h_cos_theta, _ = np.histogram(cos_theta, bins=NBINS, range=ANGLE_RANGE, density=True)
        h_cos_theta += tmp_h_theta

# normalization of the density (average per frame)

h_cos_theta /= u.trajectory.n_frames

[Oliver Beckstein]

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.

[Oliver Beckstein]

One quick comment:

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

> 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().

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

[Roja]

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, edges = np.histogram([0], bins=NBINS, range=ANGLE_RANGE, density=True)
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)

        tmp_h_cos_alpha_a, _ = np.histogram(cos_alpha, bins=NBINS, range=ANGLE_RANGE, density=True)
        h_cos_alpha_a += tmp_h_cos_alpha_a

# normalization of the density (average per frame)

h_cos_alpha_a /= u.trajectory.n_frames
------------

[Oliver Beckstein]

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.

[Roja]

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

[Oliver Beckstein]

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.