Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .travis.yml
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ install:

script:
- pytest --cov=./

- flake8

after_success:
- codecov
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file modified orientation/__pycache__/setup_system.cpython-37.pyc
Binary file not shown.
Binary file modified orientation/__pycache__/visuals.cpython-37.pyc
Binary file not shown.
40 changes: 22 additions & 18 deletions orientation/calc_angles.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import numpy as np


def get_com(universe, selection):

''' Return the centre of mass of a selection string '''
Expand All @@ -13,60 +14,64 @@ def get_principal_axes(universe, selection):

# Methodology taken from Beckstein
# https://stackoverflow.com/questions/49239475/
# how-to-use-mdanalysis-to-principal-axes-and-moment-of-inertia-with-a-group-of-at/49268247#49268247
# how-to-use-mdanalysis-to-principal-axes-and-
# moment-of-inertia-with-a-group-of-at/49268247#49268247

# select alpha carbons only
CA = universe.select_atoms(selection)

# calculate moment of inertia, and sort eigenvectors
I = CA.moment_of_inertia()
# PCB renamed I to Inertia (ambiguous variable name)
Inertia = CA.moment_of_inertia()

# UT = CA.principal_axes()
# U = UT.T

values, evecs = np.linalg.eigh(I)
values, evecs = np.linalg.eigh(Inertia)
indices = np.argsort(values)
U = evecs[:, indices]

## Below is for testing ##
# Lambda = U.T.dot(I.dot(U))
# Below is for testing ##
# Lambda = U.T.dot(Inertia.dot(U))
#
# print(Lambda)
# print(np.allclose(Lambda - np.diag(np.diagonal(Lambda)), 0))

#print("")
#print(U)
# print("")
# print(U)

return U


def dir_cosine(v1, v2):

'''
'''
Given two vectors (v1 and v2) work out the direction cosine between
them. For more information see: https://en.wikipedia.org/wiki/Direction_cosine)
them. For more information see:
https://en.wikipedia.org/wiki/Direction_cosine)
'''

direction_cosine = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2))
direction_cosine = np.dot(v1, v2) / (np.linalg.norm(v1) *
np.linalg.norm(v2))

return direction_cosine


def make_direction_cosine_matrix(ref, axes_set):

'''
Given a set of two axes (ref_option and axes_set) work out the direction cosine matrix
between them.
Given a set of two axes (ref_option and axes_set)
work out the direction cosine matrix between them.
'''
# Gather the individual vectors of our reference basis
ex = ref[0, :]
ey = ref[1, :]
ez = ref[2, :]

# principal axes calculated at this instance
u = axes_set[0, :] # 1st princpal axis (PA)
v = axes_set[1, :] # 2nd princpal axis (PA)
w = axes_set[2, :] # 3rd princpal axis (PA)
u = axes_set[0, :] # 1st princpal axis (PA)
v = axes_set[1, :] # 2nd princpal axis (PA)
w = axes_set[2, :] # 3rd princpal axis (PA)

# Work out the dot prod b/w the 1st PA and the unit vectors
u_alpha = dir_cosine(u, ex)
Expand All @@ -83,9 +88,8 @@ def make_direction_cosine_matrix(ref, axes_set):
w_beta = dir_cosine(w, ey)
w_gamma = dir_cosine(w, ez)

# make all of the above into a 3 X 3 matrix and return it
# make all of the above into a 3 X 3 matrix and return it
c_matrix = np.asarray([[u_alpha, u_beta, u_gamma],
[v_alpha, v_beta, v_gamma],
[w_alpha, w_beta, w_gamma]])

[w_alpha, w_beta, w_gamma]])
return c_matrix
Loading