[For Science!] Reconstructing rotational data with MEL.

#1

So I actually posted at the start of the week at the Autodesk Forum, but got no replies there. So I hope this community can help me out:

As a life sciences student I’m exploring the capabilities of using Maya to make clear visualizations of protein dynamics simulations. (PELE in this case - protein energy landscape exploration).

Specifically I’m working on a hormone molecule/ligand that can move through a protein over the simulated trajectories.

Maya is very useful in this regard as it can interpolate between simulation time-steps (frames), allowing for smooth animations.

The ligands were imported by script, having their translation and rotation all at 0.

I have been successful in generating an animation of the simulation output by MEL script using the position data of the pivots of the ligands.

Now the animation would be perfect if I could also assign a rotation to each molecule, as they are all exactly the same. Thus it might be possible to write a script for this, that rotates until 2 molecules are perfectly aligned, writing the rotational value to a text file?

I’m really curious if you guys can think of any way this might be done, if possible at all. See the attached picture for the scene we are talking about.

It has to be scripted, in the picture I show only 5 objects (all with Unknown rotation aka 0),
however my full data set comprises of no less then 7170 instances… (which I handle in small chunks).

Translation data is calculated by float $pos[] = xform -q -ws -piv $ligandName;

Polygons only scene (no mMaya needed):
https://www.dropbox.com/s/ybqn0r7qgkiz6gz/1A_1-10.mb?dl=1

Used mMaya plugin: https://clarafi.com/tools/mmaya/

I’m really curious to what you guys think. Madness or might it be feasible?
If you know of any other programs at all which might be able to do these alignments of polygons and query rotational values it would already be a huge help.

Thanks a lot for your input!

Cheers,
Wilglide

#2

Well, this should be doable. Fist you need a rotation for every object. This could be done by:

[ul]
[li]calculate the center position of each molecule group create a locator there[/li][li]grab two arbitrary atoms in this molecule and create locators there[/li][li]now you have three points. Calculate the vector between center and atom 1 and then between center and atom 2[/li][li]with the two vectors you can calculate the cross vector and a normal[/li][li]with these data you can create a matrix.[/li][/ul]Or…

[ul]
[li]use the three locators and create an aim contraint from the center to atom 1 and use atom b as up object[/li][li]Now you have a rotated center[/li][li]parent the molecule below the locator[/li][/ul]Once you have rotated objects you can animate them.

#3

Genius! :slight_smile:
Thanks a lot, this sounds feasible indeed.
The mathematics I should be able to do. However, the implementation in Maya might be quite a ‘rabbit hole’.
I think I will also need a try function that minimizes the rotation angle depending on where your normals are assigned then…

Some questions I still have:
[ul]
[li]Do you know if I can use this plug-in below to build the rotation node?
[/li]http://download.autodesk.com/us/maya/2010help/api/build_rotation_node_8cpp-example.html#_a11
[li]How can I assign a rotation, without actually transforming the node?
[/li][li]Do know a best practice to build these normals in Maya on the 2 vector plane?
[/li][/ul]

I would be very grateful if you would be willing to help me break down the C++ code in the link…
Do I understand correctly that it assumes upX, upY, upZ, forwardX, forwardY and forwardZ to be known for a given obj?
These are required to be orthogonal already, so I guess this could be 1 vector (of the 2 vector plane) and a determined orthogonal normal?

Thanks again. :bowdown:

#4

I’m quite busy at the moment, give me some days.

#5

Sure I understand, no rush, just for long run fun that is also very useful.
also still have 3 reports to finish by the end of next week.
I’m lucky this week if I have 2 hours free-time in the evening…

So I will make a start in messing around with this only in 2 weeks and will try to post an update then.
Part will be outside Maya because the atom coordinates are in these .pdb files such as.

https://www.dropbox.com/s/typ0h3x1mobj2s8/traj_COP1A_.1.pdb_Sep_1.txt?dl=1
https://www.dropbox.com/s/bc1xannaoz5qazr/traj_COP1A_.1.pdb_Sep_2.txt?dl=1

(says .txt but is .pdb format)

I worked with Unix in Ubuntu a bit before to manage pdb files, so probable awk and grep method.
Or might play around in Python.

Cheers!

#6

i don’t want to go back to MEL, to there is what I have with python and pymel (just for practice):

from pymel.core.datatypes import Vector, Matrix, Point
import pymel.core as pm

def matrixFromNormal(front, side):
    # normalize vectors
    front.normalize()
    side.normalize()

    #get the up axis as the cross vector   
    up = Vector.cross(front, side)
    up.normalize()
    #cross again in case front and side were not orthogonal
    side = Vector.cross(up, front)
    side.normalize()
    
    #matrix
    return Matrix (
        side.x, side.y, side.z, 0, 
        up.x, up.y, up.z, 0, 
        front.x, front.y, front.z, 0, 
        0, 0, 0, 1)
          

def redoTransform(node = None):
    ss = pm.ls(sl=1)
    node = pm.ls(node)[0] if node else ss[0]
    if node:
        pm.makeIdentity(node, apply=1, t=1, r=1, s=1);
        f = 0
        pp = []

        for k in range(3): # find three first shells of mesh, and their centers
            pm.select(node + (".f[%i]" % f), r=1)
            pm.mel.eval("polyConvertToShell")
            ff = pm.filterExpand(sm=34)
            f += len(ff)
        
            bb = pm.xform(q=1, bb=1, ws=1) # world space
            center = (Point(bb[0:3]) - Point(bb[3:0])) / 2.0
            pp.append(center)
        
        tm = matrixFromNormal(pp[1]-pp[0], pp[2]-pp[0]) # make matrix using three points (plane)
        node.setMatrix(tm.inverse())
        pm.makeIdentity(node, apply=1, r=1, s=1);
        node.setMatrix(tm)
        pm.select(ss)
        return tm

redoTransform()

i couldn’t open example file, so I expect these molecules are single-mesh polys. hopefully three central pieces are made first.

to try ->>> select a molecule and run the script. also you can pass a name to redoTransform function.

redoTransform('ligandName')
#7

Hi denisT,

I am not used to pyMel yet, but your code gives an impression of how to build the matrix using side, up and front.

I have tried to test your script, however I’m unable to split the atoms in separate polygons.
It is an nParticle to polygon conversion, maybe someone else knows if these nParticles can be converted individually.

My own idea is to just work outside of Maya on the text PDB files outputting matrix data to a text file, because the atom coordinates there match those in Maya and the individual atoms of the ligand have their center position in there (the molecules are just point (atom) clouds in PDB format, in contrast to Maya).

Both approaches are worthwhile exploring though.

Cheers,
Wilglide

#8

i can’t open test file, because it’s binary and probably saved in a different to my version of maya. try to post ZIP of maya ASCII file.

#9

Ok I didn’t know, here is the ASCII file for the polygon ligands:
https://www.dropbox.com/s/gtv8gzjgivwgxii/thisPDBPolygonTestFile_ASCII.zip?dl=1

However as I said, the atom coordinate information is lost here, the coordinates of the first 2 ligands is still in the PDB files I also linked however.

The nParticle format requires mMaya to be installed, but I don’t presume the atom coordinates is preserved in there either…

(I use Maya 2015 to be sure)

#10

So I am back a bit later then anticipated, this worked out for me.

Here is a proper method to calculate rotation(x,y,z) that I used:


#CC: EulerAngles.c by Ken Shoemake, and deals with the case where cos(y) is close to zero
#https://afni.nimh.nih.gov/pub/dist/src/pkundu/meica.libs/nibabel/eulerangles.py

import math

cy_thresh = None
M = np.asarray(A1_Rot)
if cy_thresh is None:
    try:
         cy_thresh = np.finfo(M.dtype).eps * 4
    except ValueError:
        cy_thresh = _FLOAT_EPS_4
r11, r12, r13, r21, r22, r23, r31, r32, r33 = M.flat
# cy: sqrt((cos(y)*cos(z))**2 + (cos(x)*cos(y))**2)
cy = math.sqrt(r33*r33 + r23*r23)
if cy > cy_thresh: # cos(y) not close to zero, standard form
    z = math.atan2(-r12,  r11) # atan2(cos(y)*sin(z), cos(y)*cos(z))
    y = math.atan2(r13,  cy) # atan2(sin(y), cy)
    x = math.atan2(-r23, r33) # atan2(cos(y)*sin(x), cos(x)*cos(y))
else: # cos(y) (close to) zero, so x -> 0.0 (see above)
    # so r21 -> sin(z), r22 -> cos(z) and
    z = math.atan2(r21,  r22)
    y = math.atan2(r13,  cy) # atan2(sin(y), cy)
    x = 0.0

#x,y,z are messed up, so be careful
#I found -x,y,z to correspond actually to z,y,x
#I am too lazy to adjust the formula above, order of input cross vector or order of A1_rot (see below)
A1_Rot_vec=np.array([-x,y,z])

Rotations are calculated respective to the X,Y,Z system (0,0,1) (0,1,0) (1,0,0).
Because my first molecule is always at the same position and rotation of each simulation trial,
I substract the rotation of that first molecule from the 2nd … Nth.

In this way I can leave the rotation of the first at rotX=0 rotY=0 rotZ=0 and adjust the follow up molecules respective to the first molecule (rotated) coordinate system.

I use 3 atoms to calculate the 3 vector coordinate system:

 
import numpy as np

#molecule 1 [[L1][L2][L3]], atom1 atom2 atom3
M1 = np.array([[40.672,23.105,93.316],[39.990,23.781,92.284],[41.900,23.610,93.791]])

L1 = M1[0]
L2 = M1[1]
L3 = M1[2]

# L2 as vector
# L2-L1

# L3 as vector
# L3-L1

#L4 (cross right-hand rule)
L4 = np.cross(  L3-L1, L2-L1 )
#print L4

#new L3' (cross right-hand rule)
L33 = np.cross( L2-L1, L4 )
#print L33

#A1_rotational_centre
A1_Rot = np.array([L2-L1,L33,L4])

Now I am going to write a general method that can run over all my coordinate PDB files and extract and write the [xRot,yRot,zRot] to a text file, that is relative to the rotation of the first molecule of each trial.

Maybe someone can use this.
Thanks for the help guys.

If there are any questions I am happy to try to answer them.