Calculation angle between two plane

[MOHD HOMAIDUR RAHMAN]

Dear MDAnalysis User's

I am trying to calculate the angle between two normal vectors of Imidazolium ring plane. I did for cross product and then did dot product to get all angle between each residue.

For this, I wrote a script as

residue_ids = np.arange(0, 1024, 2)
residues = traj.atoms.residues[residue_ids]
pi_angle = []
for res in residues:
    C1_C = traj.select_atoms("name C1").center_of_mass()
    C2_B = traj.select_atoms("name C2").center_of_mass()
    C3_A = traj.select_atoms("name C4").center_of_mass()    
    BA = C3_A - C2_B
    BC = C1_C - C2_B
    pi_angle.append(np.cross(BA , BC))
pi_angle = np.array(pi_angle)
theta = []
k=0
for i in range(5):
    for j in range(i+1, 512):
        k +=1
        theta.append(np.arccos(np.dot(pi_angle[i],pi_angle[j])/(norm(pi_angle[i])*norm(pi_angle[j]))))

In my simulation box I have total 1024 molecules (512 imidazolium cation and 512 anions). I selected three atoms from imidazolium ring and did cross product of selected vector to get a vector normal to the ring. I got the value for all 512 imidazolium cation but all the value are same. To cross check my value I used ambertool to calculate corrplane of my imidazolium cation for one molecule, I got a different value than this.

So I think I am doing some mistake in selections of atoms. I have total 5 ring atoms name are C1 , N6 , C2 , C4 , N7 .

Please guide me to resolve this issue or if any module is there in MDAnalysis pakage to direct print corrplane  of selected molecules or groups of atoms.

Thank you  in advance

Regards

Rahman

[Oliver Beckstein]

Hi Rahman,

Currently you are selecting *all* atoms of the name C1, all atoms C2 etc, i.e. you’re averaging over all you residues.

Instead, I assume you want to do this individually?

> On 13 Jun, 2016, at 17:14, MOHD HOMAIDUR RAHMAN <rahm...@gmail.com> wrote:
>
> for res in residues:
> C1_C = traj.select_atoms("name C1").center_of_mass()
> C2_B = traj.select_atoms("name C2").center_of_mass()
> C3_A = traj.select_atoms(“name C4”).center_of_mass()
> BA = C3_A - C2_B
> BC = C1_C - C2_B
> pi_angle.append(np.cross(BA , BC))
>

I’d write a function

def ring_plane_vector(res):
C1_C = res.atoms.select_atoms(“name C1”).positions[0]
C2_B = res.atoms.select_atoms(“name C2”).positions[0]
C3_A = res.atoms.select_atoms(“name C4”).positions[0]

BA = C3_A - C2_B
BC = C1_C - C2_B

return np.cross(BA , BC)

and then do something like

pi_angle = [ring_plane_vector(res) for res in residues]


There’s also different way to access atoms using the so-called “quick selectors”:

def ring_plane_vector(res):
C1_C = res.C1.position
C2_B = res.C2.position
C3_A = res.C4.position

BA = C3_A - C2_B
BC = C1_C - C2_B

return np.cross(BA , BC)

Should also work but is a bit less flexible and can give rather cryptic error messages if your residue does not contain the atom name.

Oliver

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst * orbe...@gmail.com

[MOHD HOMAIDUR RAHMAN]

Thank you Oliver for your quik reply

For my next doubt more specifically, could I run a for loop over all cation residue and groups of atom.
Like in my simulation box I have 512 cation and I want to select only ring atoms of cation with center of geometry or center of mass

cog=[]
For res in resid:
    ring = res.atoms.select_atoms(“name C1 N6 C2 C4 N7 ”).center_of_geometry
    cog.append(ring)

It is possible to use box information in this calculation. My aim is to calculate distance between two ring.

[Richard Gowers]

Hi Rahman,

If you're trying to calculate the distances between the CoG of each ring....

from MDAnalysis.lib.distances import distance_array

# your above code to calculate the CoG of each ring
cog = np.array(cog, dtype=np.float32)

# pass in Universe.dimensions as the box kwarg to use box information
# alternatively, you can just pass in box=np.array([10, 10, 10, 90, 90, 90]) for a 10x10x10 cubic box
dist_matrix = distance_array(cog, cog, box=u.dimensions)

Then dist_matrix should be a 512x512 array of the distance between each ring

Richard

[MOHD HOMAIDUR RAHMAN]

Thanks Richard

I read about this module from tutorials  at [link] and my confusion is about atom selection as I mention earlier.
I am loading amber topology and coordinate file. The trajectory box information are present in the coordinate file for each time step. Can I extract the box information over time step and the coordinate of selected atoms ( only rings atoms ) by center of geometry . So I can give this as input in distance array module.

I want to select res id and atoms  C1 N6 C2 C4 N7 by center of geometry.

Thanks

Rahman

[Oliver Beckstein]

> On 15 Jun, 2016, at 05:59, MOHD HOMAIDUR RAHMAN <rahm...@gmail.com> wrote:
>
> I am loading amber topology and coordinate file. The trajectory box information are present in the coordinate file for each time step. Can I extract the box information over time step and the coordinate of selected atoms ( only rings atoms ) by center of geometry . So I can give this as input in distance array module.

You can get the box information from the Timestep’s dimension attribute (provided the trajectory has unit cell information):

for ts in u.trajectory:
box = ts.dimensions
dist = distance_array(…., box=box)

[Oliver Beckstein]

P.S.: Richard already gave you another way to use the dimension information: Just use Universe.dimension:

> On 15 Jun, 2016, at 04:25, Richard Gowers <richard...@gmail.com> wrote:
>
> dist_matrix = distance_array(cog, cog, box=u.dimensions)
>

The dimensions attribute automatically updates whenever you iterate through the trajectory.