13726590978

Committed 07 Mar 2025 06:07PM UTC coverage: 1.726%. First build

Build # 13726590978

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #45 from CCPBioSim/39-implement-ci-pipeline-pre-commit-hooks

Implement CI pipeline pre-commit hooks

Coverage Stats

27 of 1135 new or added lines in 23 files covered. (2.38%)

53 of 3071 relevant lines covered (1.73%)

0.05 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

0.0

/CodeEntropy/poseidon/extractData/generalFunctions.py

#!/usr/bin/env python

import logging
import math

import numpy as np
from numpy import linalg as LA

# import sys


def distance(x0, x1, dimensions):
    """
    calculates distance and accounts for PBCs.
    """

    x0 = np.array(x0)
    x1 = np.array(x1)
    delta = np.abs(x1 - x0)
    delta = np.where(delta > 0.5 * dimensions, delta - dimensions, delta)
    dist = np.sqrt((delta**2).sum(axis=-1))
    return dist


def vector(x0, x1, dimensions):
    """
    get vector of two coordinates over PBCs.
    """

    x0 = np.array(x0)
    x1 = np.array(x1)
    dimensions = np.array(dimensions)
    delta = x1 - x0
    delta2 = []
    for delt, box in zip(delta, dimensions):
        if delt < 0 and delt < 0.5 * box:
            delt = delt + box
        if delt > 0 and delt > 0.5 * box:
            delt = delt - box

        delta2.append(delt)

    delta2 = np.array(delta2)

    return delta2


def getMcoords(x0, x1, dimensions):
    """
    Get the middle (M) coords of two coords.
    """

    x0 = np.array(x0)
    x1 = np.array(x1)
    dimensions = np.array(dimensions)

    delta = x1 - x0

    delta2 = []
    for delt, box in zip(delta, dimensions):
        if delt < 0 and delt < 0.5 * box:
            delt = delt + box
        if delt > 0 and delt > 0.5 * box:
            delt = delt - box

        delta2.append(delt)

    delta2 = np.array(delta2)

    M_coords2 = 0.5 * delta2 + x0

    return M_coords2


def projectVectorToPlane(vector, plane_normal):
    """
    Project a vector onto a plane, need normal of plane as input.
    """

    plane_normal_magnitude = np.linalg.norm(plane_normal)

    project_vector = np.cross(vector, plane_normal)
    project_vector = np.divide(project_vector, plane_normal_magnitude)
    project_vector = np.cross(plane_normal, project_vector)
    projected_vector = np.divide(project_vector, plane_normal_magnitude)

    return projected_vector


def normaliseVector(x0, x1, dimensions):
    """
    Get vector of two coordinates over PBCs and then normalise by the
    magnitude of the vector
    """

    x0 = np.array(x0)
    x1 = np.array(x1)
    dimensions = np.array(dimensions)
    delta = x1 - x0

    delta2 = []
    for delt, box in zip(delta, dimensions):
        if delt < 0 and delt < 0.5 * box:
            delt = delt + box
        if delt > 0 and delt > 0.5 * box:
            delt = delt - box

        delta2.append(delt)

    delta = np.array(delta2)

    delta_magnitude = (delta[0] ** 2 + delta[1] ** 2 + delta[2] ** 2) ** 0.5
    delta_norm = np.divide(delta, delta_magnitude)

    return delta_norm


def angleBetweenVectors(x0, x1):
    """
    Get the angle between two vectors to determine how orthogonal they are.
    """

    top = np.dot(x0, x1)

    axis1_magnitude = np.linalg.norm(x0)  # magnitude/ distance
    axis2_magnitude = np.linalg.norm(x1)  # magnitude/ distance
    bottom = np.multiply(axis1_magnitude, axis2_magnitude)
    top_bottom = np.divide(top, bottom)
    if top_bottom > 1:
        top_bottom = 1
    if top_bottom < -1:
        top_bottom = -1
    theta = np.arccos(top_bottom)
    theta = np.degrees(theta)

    return theta


def angle2(a, b, c, dimensions):
    """
    Just used as a check. Doesn't work as periodic boundary conditions
    are not considered properly.
    """

    a = np.array(a)  # j
    b = np.array(b)  # i
    c = np.array(c)  # k
    dimensions = np.array(dimensions)
    ba = a - b
    bc = c - b
    ac = c - a
    ba = np.where(ba > 0.5 * dimensions, ba - dimensions, ba)
    bc = np.where(bc > 0.5 * dimensions, bc - dimensions, bc)
    ac = np.where(ac > 0.5 * dimensions, ac - dimensions, ac)

    cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))

    angle = np.arccos(cosine_angle)
    angle = np.degrees(angle)

    return angle


def angle(a, b, c, dimensions):

    a = np.array(a)  # j
    b = np.array(b)  # i
    c = np.array(c)  # k
    dimensions = np.array(dimensions)
    ba = np.abs(a - b)
    bc = np.abs(c - b)
    ac = np.abs(c - a)
    ba = np.where(ba > 0.5 * dimensions, ba - dimensions, ba)
    bc = np.where(bc > 0.5 * dimensions, bc - dimensions, bc)
    ac = np.where(ac > 0.5 * dimensions, ac - dimensions, ac)
    dist_ba = np.sqrt((ba**2).sum(axis=-1))
    dist_bc = np.sqrt((bc**2).sum(axis=-1))
    dist_ac = np.sqrt((ac**2).sum(axis=-1))

    cosine_angle = (dist_ac**2 - dist_bc**2 - dist_ba**2) / (-2 * dist_bc * dist_ba)

    return cosine_angle


def com_vectors(coord_list, mass_list, dimensions):
    """
    doesnt work
    """
    cm = np.zeros(3)
    tot_mass = 0
    for coord, mass in zip(coord_list, mass_list):
        vec = vector(coord, [0, 0, 0], dimensions)
        cm += np.array(vec) * mass
        tot_mass += mass

    cm_weighted = np.array(cm) / float(tot_mass)

    return cm_weighted


def com(coord_list, mass_list):
    """ """
    x_cm = 0
    y_cm = 0
    z_cm = 0
    tot_mass = 0
    for coord, mass in zip(coord_list, mass_list):
        x_cm += mass * coord[0]
        y_cm += mass * coord[1]
        z_cm += mass * coord[2]
        tot_mass += mass
    x_cm = float(x_cm) / float(tot_mass)
    y_cm = float(y_cm) / float(tot_mass)
    z_cm = float(z_cm) / float(tot_mass)

    cm = (x_cm, y_cm, z_cm)

    return cm


def MOI(cm, coord_list, mass_list):
    """ """
    x_cm, y_cm, z_cm = cm[0], cm[1], cm[2]
    I_xx = 0
    I_xy = 0
    I_yx = 0

    I_yy = 0
    I_xz = 0
    I_zx = 0

    I_zz = 0
    I_yz = 0
    I_zy = 0

    for coord, mass in zip(coord_list, mass_list):
        I_xx += (abs(coord[1] - y_cm) ** 2 + abs(coord[2] - z_cm) ** 2) * mass
        I_xy += (coord[0] - x_cm) * (coord[1] - y_cm) * mass
        I_yx += (coord[0] - x_cm) * (coord[1] - y_cm) * mass

        I_yy += (abs(coord[0] - x_cm) ** 2 + abs(coord[2] - z_cm) ** 2) * mass
        I_xz += (coord[0] - x_cm) * (coord[2] - z_cm) * mass
        I_zx += (coord[0] - x_cm) * (coord[2] - z_cm) * mass

        I_zz += (abs(coord[0] - x_cm) ** 2 + abs(coord[1] - y_cm) ** 2) * mass
        I_yz += (coord[1] - y_cm) * (coord[2] - z_cm) * mass
        I_zy += (coord[1] - y_cm) * (coord[2] - z_cm) * mass

    inertia_tensor = np.array(
        [[I_xx, -I_xy, -I_xz], [-I_yx, I_yy, -I_yz], [-I_zx, -I_zy, I_zz]]
    )
    return inertia_tensor


def UAforceTorque(bonded_atoms_list, F_axes, Fmoi_axes, MI_axes):
    """
    Get UA force and torque from previously calculated Faxes and MIaxes
    """

    """
    Faxis_list = (Faxis1, Faxis2, Faxis3)
    f_list = (FO, FH1, FH2)
    MI_xyz = (axis1MI, axis2MI, axis3MI)
    water_torques = torque(c_list, m_list, Faxis_list, f_list, MI_xyz)
    """

    c_list = []
    f_list = []
    forces_summed = np.zeros(3)
    m_list = []

    for atoms in bonded_atoms_list:
        c_list.append(atoms.coords)
        f_list.append(atoms.forces)
        forces_summed += atoms.forces
        m_list.append(atoms.mass)

    cm = com(c_list, m_list)
    UA_torque = torque_old(cm, c_list, Fmoi_axes, f_list, MI_axes)

    F1 = np.dot(forces_summed, F_axes[0])
    F2 = np.dot(forces_summed, F_axes[1])
    F3 = np.dot(forces_summed, F_axes[2])
    mass_sqrt = sum(m_list) ** 0.5
    UA_force = (
        float(F1) / float(mass_sqrt),
        float(F2) / float(mass_sqrt),
        float(F3) / float(mass_sqrt),
    )

    return UA_force, UA_torque


def torque_old(cm, coord_list, axis_list, force_list, MI_coords):
    """ """
    x_cm, y_cm, z_cm = cm[0], cm[1], cm[2]
    MI_coords = np.sqrt(MI_coords)  # sqrt moi to weight torques

    O_torque = np.array([0, 0, 0])
    H1_torque = np.array([0, 0, 0])
    H2_torque = np.array([0, 0, 0])
    count_atom = 0
    for coord, force in zip(coord_list, force_list):
        count_atom += 1
        new_coords = []
        new_forces = []
        for axis in axis_list:
            new_c = (
                axis[0] * (coord[0] - x_cm)
                + axis[1] * (coord[1] - y_cm)
                + axis[2] * (coord[2] - z_cm)
            )
            new_f = axis[0] * force[0] + axis[1] * force[1] + axis[2] * force[2]
            new_coords.append(new_c)
            new_forces.append(new_f)

        torque_list = []
        torquex = float(
            new_coords[1] * new_forces[2] - new_coords[2] * new_forces[1]
        )  # / float(1e10)
        torquey = float(
            new_coords[2] * new_forces[0] - new_coords[0] * new_forces[2]
        )  # / float(1e10)
        torquez = float(
            new_coords[0] * new_forces[1] - new_coords[1] * new_forces[0]
        )  # / float(1e10)

        if count_atom == 1:
            O_torque = (torquex, torquey, torquez)
            O_torque = np.divide(O_torque, MI_coords)
        elif count_atom == 2:
            H1_torque = (torquex, torquey, torquez)
            H1_torque = np.divide(H1_torque, MI_coords)
        elif count_atom == 3:
            H2_torque = (torquex, torquey, torquez)
            H2_torque = np.divide(H2_torque, MI_coords)
        else:
            ("torques outs of range")
            break

    W_torque = O_torque + H1_torque + H2_torque

    return W_torque


def principalAxesMOI(bonded_atoms_list):
    """
    Calculate the principal axes of the MOI matrix, then calc forces
    """

    coord_list = []
    mass_list = []
    forces_summed = np.zeros(3)

    for atom in bonded_atoms_list:
        coord_list.append(atom.coords)
        mass_list.append(atom.mass)
        forces_summed += atom.forces

    cm = com(coord_list, mass_list)
    moI = MOI(cm, coord_list, mass_list)

    eigenvalues, eigenvectors = LA.eig(moI)
    # different values generated to Jon's code

    transposed = np.transpose(eigenvectors)  # turn columns to rows

    min_eigenvalue = abs(eigenvalues[0])
    if eigenvalues[1] < min_eigenvalue:
        min_eigenvalue = eigenvalues[1]
    if eigenvalues[2] < min_eigenvalue:
        min_eigenvalue = eigenvalues[2]

    max_eigenvalue = abs(eigenvalues[0])
    if eigenvalues[1] > max_eigenvalue:
        max_eigenvalue = eigenvalues[1]
    if eigenvalues[2] > max_eigenvalue:
        max_eigenvalue = eigenvalues[2]

    # print (min_eigenvalue * const, max_eigenvalue * const)
    # same as Jons

    FMIaxis1 = None
    FMIaxis2 = None
    FMIaxis3 = None

    axis1MI = None
    axis2MI = None
    axis3MI = None

    for i in range(0, 3):
        if eigenvalues[i] == max_eigenvalue:
            FMIaxis1 = transposed[i]
            axis1MI = eigenvalues[i]
        elif eigenvalues[i] == min_eigenvalue:
            FMIaxis3 = transposed[i]
            axis3MI = eigenvalues[i]
        else:
            FMIaxis2 = transposed[i]
            axis2MI = eigenvalues[i]

    Fm_axes = [FMIaxis1, FMIaxis2, FMIaxis3]
    MI_axes = [axis1MI, axis2MI, axis3MI]

    return Fm_axes, MI_axes


def calcAngleWithNearestNonlike(neighbour, nearest, dimensions):
    """
    Calc angle in plane and orthog to plane between any molecule
    and any other molecule.
    """

    PAxes = np.array(neighbour.WMprincipalAxis[0])
    PAxis = np.array(neighbour.WMprincipalAxis[1])
    COM = np.array(neighbour.WMprincipalAxis[2])

    # plane_normal = vector(COM, PAxes[0], dimensions)
    # plane_normal2 = vector(COM, PAxes[1], dimensions)
    # OM_vector = vector(COM, PAxes[2], dimensions)

    plane_normal = PAxes[0]
    plane_normal2 = PAxes[1]
    OM_vector = PAxes[2]

    OSolute_vector = vector(nearest.coords, COM, dimensions)

    projected_OSolute = projectVectorToPlane(OSolute_vector, plane_normal)
    projected_orthog_OSolute = projectVectorToPlane(OSolute_vector, plane_normal2)

    # print ('OM', OM_vector, oxygen.atom_num)

    soluteOMangle = angleBetweenVectors(OM_vector, projected_OSolute)
    soluteOorthogAngle = angleBetweenVectors(OM_vector, projected_orthog_OSolute)

    if (
        np.dot(OM_vector, projected_OSolute) > 0
        and np.dot(plane_normal2, projected_OSolute) > 0
    ):
        soluteOMangle = (90 - soluteOMangle) + 270

    if (
        np.dot(-OM_vector, projected_OSolute) > 0
        and np.dot(plane_normal2, projected_OSolute) > 0
    ):
        soluteOMangle = (180 - soluteOMangle) + 180

    if (
        np.dot(-plane_normal, projected_orthog_OSolute) > 0
        and np.dot(-OM_vector, projected_orthog_OSolute) > 0
    ):
        soluteOorthogAngle = (180 - soluteOorthogAngle) + 180

    if (
        np.dot(-plane_normal, projected_orthog_OSolute) > 0
        and np.dot(OM_vector, projected_orthog_OSolute) > 0
    ):
        soluteOorthogAngle = (90 - soluteOorthogAngle) + 270

    # cosAngle = angle(closest_H.coords, oxygen.coords, solute.coords, dimensions)
    # arccosAngle = np.arccos(cosAngle)
    # angleDegrees = np.degrees(arccosAngle)
    # if math.isnan(angleDegrees) is False:
    # angleDegrees = int(round(angleDegrees, 0))

    if math.isnan(soluteOMangle) is False:
        soluteOMangle = int(round(soluteOMangle, 0))
    # else:
    # print ('Error1: NaN for solute resid: %s, '\
    # 'water resid %s and angle %s.'\
    # % (nearest.resid, neighbour.resid, soluteOMangle))
    # continue

    if math.isnan(soluteOorthogAngle) is False:
        soluteOorthogAngle = int(round(soluteOorthogAngle, 0))
    # else:
    # print ('Error1: NaN for solute resid: %s, '\
    # 'water resid %s and angle %s.'\
    # % (nearest.resid, neighbour.resid, soluteOorthogAngle))
    # continue

    return soluteOMangle, soluteOorthogAngle


def calcWaterSoluteAngle(solute, oxygen, H1, H2, dimensions):
    """
    Calc angle in plane and orthog to plane between any water
    and any solute atom.
    """

    dist_solute_H1, dist_solute_H2 = None, None

    for atomDist in solute.nearest_all_atom_array:
        if atomDist[0] == H1.atom_num:
            dist_solute_H1 = atomDist[1]
        elif atomDist[0] == H2.atom_num:
            dist_solute_H2 = atomDist[1]
        else:
            continue

    if dist_solute_H1 is None:
        dist_solute_H1 = distance(solute.coords, H1.coords, dimensions)
    if dist_solute_H2 is None:
        dist_solute_H2 = distance(solute.coords, H2.coords, dimensions)
    closest_H = None
    if dist_solute_H1 < dist_solute_H2:
        closest_H = H1
    elif dist_solute_H1 > dist_solute_H2:
        closest_H = H2
    else:
        closest_H = H1

    H1O_vector = vector(oxygen.coords, H1.coords, dimensions)
    H2O_vector = vector(oxygen.coords, H2.coords, dimensions)
    plane_normal = np.cross(H1O_vector, H2O_vector)
    M_coords = getMcoords(H1.coords, H2.coords, dimensions)
    OM_vector = vector(M_coords, oxygen.coords, dimensions)
    plane_normal2 = np.cross(OM_vector, plane_normal)
    OSolute_vector = vector(solute.coords, oxygen.coords, dimensions)
    projected_OSolute = projectVectorToPlane(OSolute_vector, plane_normal)
    projected_orthog_OSolute = projectVectorToPlane(OSolute_vector, plane_normal2)

    # print ('OM', OM_vector, oxygen.atom_num)

    soluteOMangle = angleBetweenVectors(OM_vector, projected_OSolute)
    soluteOorthogAngle = angleBetweenVectors(OM_vector, projected_orthog_OSolute)

    if (
        np.dot(OM_vector, projected_OSolute) > 0
        and np.dot(plane_normal2, projected_OSolute) > 0
    ):
        soluteOMangle = (90 - soluteOMangle) + 270
        # only go up to 180 degrees rather than 360 as above
        # if soluteOMangle > 180:
        # soluteOMangle =

    if (
        np.dot(-OM_vector, projected_OSolute) > 0
        and np.dot(plane_normal2, projected_OSolute) > 0
    ):
        soluteOMangle = (180 - soluteOMangle) + 180

    if (
        np.dot(-plane_normal, projected_orthog_OSolute) > 0
        and np.dot(-OM_vector, projected_orthog_OSolute) > 0
    ):
        soluteOorthogAngle = (180 - soluteOorthogAngle) + 180

    if (
        np.dot(-plane_normal, projected_orthog_OSolute) > 0
        and np.dot(OM_vector, projected_orthog_OSolute) > 0
    ):
        soluteOorthogAngle = (90 - soluteOorthogAngle) + 270

    cosAngle = angle(closest_H.coords, oxygen.coords, solute.coords, dimensions)
    arccosAngle = np.arccos(cosAngle)
    angleDegrees = np.degrees(arccosAngle)
    if math.isnan(angleDegrees) is False:
        angleDegrees = int(round(angleDegrees, 0))
        soluteOMangle = int(round(soluteOMangle, 0))
        soluteOorthogAngle = int(round(soluteOorthogAngle, 0))
    else:
        logging.error(
            "NaN for solute resid: %s, water resid %s and angle %s."
            % (solute.resid, oxygen.resid, angleDegrees)
        )

    return soluteOMangle, soluteOorthogAngle

1	#!/usr/bin/env python
2
NEW 3	import logging	×
4	import math	×
5
6	import numpy as np	×
7	from numpy import linalg as LA	×
8
9	# import sys
10
11
12	def distance(x0, x1, dimensions):	×
13	"""
14	calculates distance and accounts for PBCs.
15	"""
16
17	x0 = np.array(x0)	×
18	x1 = np.array(x1)	×
19	delta = np.abs(x1 - x0)	×
20	delta = np.where(delta > 0.5 * dimensions, delta - dimensions, delta)	×
NEW 21	dist = np.sqrt((delta**2).sum(axis=-1))	×
22	return dist	×
23
24
25	def vector(x0, x1, dimensions):	×
26	"""
27	get vector of two coordinates over PBCs.
28	"""
29
30	x0 = np.array(x0)	×
31	x1 = np.array(x1)	×
32	dimensions = np.array(dimensions)	×
33	delta = x1 - x0	×
34	delta2 = []	×
35	for delt, box in zip(delta, dimensions):	×
36	if delt < 0 and delt < 0.5 * box:	×
37	delt = delt + box	×
38	if delt > 0 and delt > 0.5 * box:	×
39	delt = delt - box	×
40
41	delta2.append(delt)	×
42
43	delta2 = np.array(delta2)	×
44
45	return delta2	×
46
47
48	def getMcoords(x0, x1, dimensions):	×
49	"""
50	Get the middle (M) coords of two coords.
51	"""
52
53	x0 = np.array(x0)	×
54	x1 = np.array(x1)	×
55	dimensions = np.array(dimensions)	×
56
57	delta = x1 - x0	×
58
59	delta2 = []	×
60	for delt, box in zip(delta, dimensions):	×
61	if delt < 0 and delt < 0.5 * box:	×
62	delt = delt + box	×
63	if delt > 0 and delt > 0.5 * box:	×
64	delt = delt - box	×
65
66	delta2.append(delt)	×
67
68	delta2 = np.array(delta2)	×
69
70	M_coords2 = 0.5 * delta2 + x0	×
71
72	return M_coords2	×
73
74
75	def projectVectorToPlane(vector, plane_normal):	×
76	"""
77	Project a vector onto a plane, need normal of plane as input.
78	"""
79
80	plane_normal_magnitude = np.linalg.norm(plane_normal)	×
81
82	project_vector = np.cross(vector, plane_normal)	×
83	project_vector = np.divide(project_vector, plane_normal_magnitude)	×
84	project_vector = np.cross(plane_normal, project_vector)	×
85	projected_vector = np.divide(project_vector, plane_normal_magnitude)	×
86
87	return projected_vector	×
88
89
90	def normaliseVector(x0, x1, dimensions):	×
91	"""
92	Get vector of two coordinates over PBCs and then normalise by the
93	magnitude of the vector
94	"""
95
96	x0 = np.array(x0)	×
97	x1 = np.array(x1)	×
98	dimensions = np.array(dimensions)	×
99	delta = x1 - x0	×
100
101	delta2 = []	×
102	for delt, box in zip(delta, dimensions):	×
103	if delt < 0 and delt < 0.5 * box:	×
104	delt = delt + box	×
105	if delt > 0 and delt > 0.5 * box:	×
106	delt = delt - box	×
107
108	delta2.append(delt)	×
109
110	delta = np.array(delta2)	×
111
112	delta_magnitude = (delta[0] 2 + delta[1] 2 + delta[2] 2) 0.5	×
113	delta_norm = np.divide(delta, delta_magnitude)	×
114
115	return delta_norm	×
116
117
118	def angleBetweenVectors(x0, x1):	×
119	"""
120	Get the angle between two vectors to determine how orthogonal they are.
121	"""
122
123	top = np.dot(x0, x1)	×
124
NEW 125	axis1_magnitude = np.linalg.norm(x0) # magnitude/ distance	×
NEW 126	axis2_magnitude = np.linalg.norm(x1) # magnitude/ distance	×
127	bottom = np.multiply(axis1_magnitude, axis2_magnitude)	×
128	top_bottom = np.divide(top, bottom)	×
129	if top_bottom > 1:	×
130	top_bottom = 1	×
131	if top_bottom < -1:	×
132	top_bottom = -1	×
133	theta = np.arccos(top_bottom)	×
134	theta = np.degrees(theta)	×
135
136	return theta	×
137
138
139	def angle2(a, b, c, dimensions):	×
140	"""
141	Just used as a check. Doesn't work as periodic boundary conditions
142	are not considered properly.
143	"""
144
NEW 145	a = np.array(a) # j	×
NEW 146	b = np.array(b) # i	×
NEW 147	c = np.array(c) # k	×
148	dimensions = np.array(dimensions)	×
149	ba = a - b	×
150	bc = c - b	×
151	ac = c - a	×
152	ba = np.where(ba > 0.5 * dimensions, ba - dimensions, ba)	×
153	bc = np.where(bc > 0.5 * dimensions, bc - dimensions, bc)	×
154	ac = np.where(ac > 0.5 * dimensions, ac - dimensions, ac)	×
155
156	cosine_angle = np.dot(ba, bc) / (np.linalg.norm(ba) * np.linalg.norm(bc))	×
157
158	angle = np.arccos(cosine_angle)	×
159	angle = np.degrees(angle)	×
160
161	return angle	×
162
163
164	def angle(a, b, c, dimensions):	×
165
NEW 166	a = np.array(a) # j	×
NEW 167	b = np.array(b) # i	×
NEW 168	c = np.array(c) # k	×
169	dimensions = np.array(dimensions)	×
170	ba = np.abs(a - b)	×
171	bc = np.abs(c - b)	×
172	ac = np.abs(c - a)	×
173	ba = np.where(ba > 0.5 * dimensions, ba - dimensions, ba)	×
174	bc = np.where(bc > 0.5 * dimensions, bc - dimensions, bc)	×
175	ac = np.where(ac > 0.5 * dimensions, ac - dimensions, ac)	×
NEW 176	dist_ba = np.sqrt((ba**2).sum(axis=-1))	×
NEW 177	dist_bc = np.sqrt((bc**2).sum(axis=-1))	×
NEW 178	dist_ac = np.sqrt((ac**2).sum(axis=-1))	×
179
NEW 180	cosine_angle = (dist_ac2 - dist_bc2 - dist_ba*2) / (-2 dist_bc * dist_ba)	×
181
182	return cosine_angle	×
183
184
185	def com_vectors(coord_list, mass_list, dimensions):	×
186	"""
187	doesnt work
188	"""
189	cm = np.zeros(3)	×
190	tot_mass = 0	×
191	for coord, mass in zip(coord_list, mass_list):	×
192	vec = vector(coord, [0, 0, 0], dimensions)	×
193	cm += np.array(vec) * mass	×
194	tot_mass += mass	×
195
196	cm_weighted = np.array(cm) / float(tot_mass)	×
197
198	return cm_weighted	×
199
200
201	def com(coord_list, mass_list):	×
202	""" """
203	x_cm = 0	×
204	y_cm = 0	×
205	z_cm = 0	×
206	tot_mass = 0	×
207	for coord, mass in zip(coord_list, mass_list):	×
208	x_cm += mass * coord[0]	×
209	y_cm += mass * coord[1]	×
210	z_cm += mass * coord[2]	×
211	tot_mass += mass	×
212	x_cm = float(x_cm) / float(tot_mass)	×
213	y_cm = float(y_cm) / float(tot_mass)	×
214	z_cm = float(z_cm) / float(tot_mass)	×
215
216	cm = (x_cm, y_cm, z_cm)	×
217
218	return cm	×
219
220
221	def MOI(cm, coord_list, mass_list):	×
222	""" """
223	x_cm, y_cm, z_cm = cm[0], cm[1], cm[2]	×
224	I_xx = 0	×
225	I_xy = 0	×
226	I_yx = 0	×
227
228	I_yy = 0	×
229	I_xz = 0	×
230	I_zx = 0	×
231
232	I_zz = 0	×
233	I_yz = 0	×
234	I_zy = 0	×
235
236	for coord, mass in zip(coord_list, mass_list):	×
NEW 237	I_xx += (abs(coord[1] - y_cm) 2 + abs(coord[2] - z_cm) 2) * mass	×
238	I_xy += (coord[0] - x_cm) * (coord[1] - y_cm) * mass	×
239	I_yx += (coord[0] - x_cm) * (coord[1] - y_cm) * mass	×
240
NEW 241	I_yy += (abs(coord[0] - x_cm) 2 + abs(coord[2] - z_cm) 2) * mass	×
242	I_xz += (coord[0] - x_cm) * (coord[2] - z_cm) * mass	×
243	I_zx += (coord[0] - x_cm) * (coord[2] - z_cm) * mass	×
244
NEW 245	I_zz += (abs(coord[0] - x_cm) 2 + abs(coord[1] - y_cm) 2) * mass	×
246	I_yz += (coord[1] - y_cm) * (coord[2] - z_cm) * mass	×
247	I_zy += (coord[1] - y_cm) * (coord[2] - z_cm) * mass	×
248
NEW 249	inertia_tensor = np.array(	×
250	[[I_xx, -I_xy, -I_xz], [-I_yx, I_yy, -I_yz], [-I_zx, -I_zy, I_zz]]
251	)
NEW 252	return inertia_tensor	×
253
254
255	def UAforceTorque(bonded_atoms_list, F_axes, Fmoi_axes, MI_axes):	×
256	"""
257	Get UA force and torque from previously calculated Faxes and MIaxes
258	"""
259
NEW 260	"""	×
261	Faxis_list = (Faxis1, Faxis2, Faxis3)
262	f_list = (FO, FH1, FH2)
263	MI_xyz = (axis1MI, axis2MI, axis3MI)
264	water_torques = torque(c_list, m_list, Faxis_list, f_list, MI_xyz)
265	"""
266
267	c_list = []	×
268	f_list = []	×
269	forces_summed = np.zeros(3)	×
270	m_list = []	×
271
272	for atoms in bonded_atoms_list:	×
273	c_list.append(atoms.coords)	×
274	f_list.append(atoms.forces)	×
275	forces_summed += atoms.forces	×
276	m_list.append(atoms.mass)	×
277
278	cm = com(c_list, m_list)	×
279	UA_torque = torque_old(cm, c_list, Fmoi_axes, f_list, MI_axes)	×
280
281	F1 = np.dot(forces_summed, F_axes[0])	×
282	F2 = np.dot(forces_summed, F_axes[1])	×
283	F3 = np.dot(forces_summed, F_axes[2])	×
284	mass_sqrt = sum(m_list) ** 0.5	×
NEW 285	UA_force = (	×
286	float(F1) / float(mass_sqrt),
287	float(F2) / float(mass_sqrt),
288	float(F3) / float(mass_sqrt),
289	)
290
291	return UA_force, UA_torque	×
292
293
294	def torque_old(cm, coord_list, axis_list, force_list, MI_coords):	×
295	""" """
296	x_cm, y_cm, z_cm = cm[0], cm[1], cm[2]	×
NEW 297	MI_coords = np.sqrt(MI_coords) # sqrt moi to weight torques	×
298
299	O_torque = np.array([0, 0, 0])	×
300	H1_torque = np.array([0, 0, 0])	×
301	H2_torque = np.array([0, 0, 0])	×
302	count_atom = 0	×
303	for coord, force in zip(coord_list, force_list):	×
NEW 304	count_atom += 1	×
305	new_coords = []	×
306	new_forces = []	×
307	for axis in axis_list:	×
NEW 308	new_c = (	×
309	axis[0] * (coord[0] - x_cm)
310	+ axis[1] * (coord[1] - y_cm)
311	+ axis[2] * (coord[2] - z_cm)
312	)
NEW 313	new_f = axis[0] * force[0] + axis[1] * force[1] + axis[2] * force[2]	×
314	new_coords.append(new_c)	×
315	new_forces.append(new_f)	×
316
317	torque_list = []	×
NEW 318	torquex = float(	×
319	new_coords[1] * new_forces[2] - new_coords[2] * new_forces[1]
320	) # / float(1e10)
NEW 321	torquey = float(	×
322	new_coords[2] * new_forces[0] - new_coords[0] * new_forces[2]
323	) # / float(1e10)
NEW 324	torquez = float(	×
325	new_coords[0] * new_forces[1] - new_coords[1] * new_forces[0]
326	) # / float(1e10)
327
328	if count_atom == 1:	×
329	O_torque = (torquex, torquey, torquez)	×
330	O_torque = np.divide(O_torque, MI_coords)	×
331	elif count_atom == 2:	×
332	H1_torque = (torquex, torquey, torquez)	×
333	H1_torque = np.divide(H1_torque, MI_coords)	×
334	elif count_atom == 3:	×
335	H2_torque = (torquex, torquey, torquez)	×
336	H2_torque = np.divide(H2_torque, MI_coords)	×
337	else:
338	("torques outs of range")	×
339	break	×
340
341	W_torque = O_torque + H1_torque + H2_torque	×
342
NEW 343	return W_torque	×
344
345
346	def principalAxesMOI(bonded_atoms_list):	×
347	"""
348	Calculate the principal axes of the MOI matrix, then calc forces
349	"""
350
351	coord_list = []	×
352	mass_list = []	×
353	forces_summed = np.zeros(3)	×
354
355	for atom in bonded_atoms_list:	×
356	coord_list.append(atom.coords)	×
357	mass_list.append(atom.mass)	×
358	forces_summed += atom.forces	×
359
360	cm = com(coord_list, mass_list)	×
361	moI = MOI(cm, coord_list, mass_list)	×
362
NEW 363	eigenvalues, eigenvectors = LA.eig(moI)	×
364	# different values generated to Jon's code
365
NEW 366	transposed = np.transpose(eigenvectors) # turn columns to rows	×
367
368	min_eigenvalue = abs(eigenvalues[0])	×
369	if eigenvalues[1] < min_eigenvalue:	×
370	min_eigenvalue = eigenvalues[1]	×
371	if eigenvalues[2] < min_eigenvalue:	×
372	min_eigenvalue = eigenvalues[2]	×
373
374	max_eigenvalue = abs(eigenvalues[0])	×
375	if eigenvalues[1] > max_eigenvalue:	×
376	max_eigenvalue = eigenvalues[1]	×
377	if eigenvalues[2] > max_eigenvalue:	×
378	max_eigenvalue = eigenvalues[2]	×
379
380	# print (min_eigenvalue * const, max_eigenvalue * const)
381	# same as Jons
382
383	FMIaxis1 = None	×
384	FMIaxis2 = None	×
385	FMIaxis3 = None	×
386
387	axis1MI = None	×
388	axis2MI = None	×
389	axis3MI = None	×
390
NEW 391	for i in range(0, 3):	×
392	if eigenvalues[i] == max_eigenvalue:	×
393	FMIaxis1 = transposed[i]	×
394	axis1MI = eigenvalues[i]	×
395	elif eigenvalues[i] == min_eigenvalue:	×
396	FMIaxis3 = transposed[i]	×
397	axis3MI = eigenvalues[i]	×
398	else:
399	FMIaxis2 = transposed[i]	×
400	axis2MI = eigenvalues[i]	×
401
402	Fm_axes = [FMIaxis1, FMIaxis2, FMIaxis3]	×
403	MI_axes = [axis1MI, axis2MI, axis3MI]	×
404
405	return Fm_axes, MI_axes	×
406
407
408	def calcAngleWithNearestNonlike(neighbour, nearest, dimensions):	×
409	"""
410	Calc angle in plane and orthog to plane between any molecule
411	and any other molecule.
412	"""
413
414	PAxes = np.array(neighbour.WMprincipalAxis[0])	×
415	PAxis = np.array(neighbour.WMprincipalAxis[1])	×
416	COM = np.array(neighbour.WMprincipalAxis[2])	×
417
418	# plane_normal = vector(COM, PAxes[0], dimensions)
419	# plane_normal2 = vector(COM, PAxes[1], dimensions)
420	# OM_vector = vector(COM, PAxes[2], dimensions)
421
422	plane_normal = PAxes[0]	×
423	plane_normal2 = PAxes[1]	×
424	OM_vector = PAxes[2]	×
425
426	OSolute_vector = vector(nearest.coords, COM, dimensions)	×
427
428	projected_OSolute = projectVectorToPlane(OSolute_vector, plane_normal)	×
NEW 429	projected_orthog_OSolute = projectVectorToPlane(OSolute_vector, plane_normal2)	×
430
431	# print ('OM', OM_vector, oxygen.atom_num)
432
433	soluteOMangle = angleBetweenVectors(OM_vector, projected_OSolute)	×
NEW 434	soluteOorthogAngle = angleBetweenVectors(OM_vector, projected_orthog_OSolute)	×
435
NEW 436	if (	×
437	np.dot(OM_vector, projected_OSolute) > 0
438	and np.dot(plane_normal2, projected_OSolute) > 0
439	):
440	soluteOMangle = (90 - soluteOMangle) + 270	×
441
NEW 442	if (	×
443	np.dot(-OM_vector, projected_OSolute) > 0
444	and np.dot(plane_normal2, projected_OSolute) > 0
445	):
446	soluteOMangle = (180 - soluteOMangle) + 180	×
447
NEW 448	if (	×
449	np.dot(-plane_normal, projected_orthog_OSolute) > 0
450	and np.dot(-OM_vector, projected_orthog_OSolute) > 0
451	):
452	soluteOorthogAngle = (180 - soluteOorthogAngle) + 180	×
453
NEW 454	if (	×
455	np.dot(-plane_normal, projected_orthog_OSolute) > 0
456	and np.dot(OM_vector, projected_orthog_OSolute) > 0
457	):
458	soluteOorthogAngle = (90 - soluteOorthogAngle) + 270	×
459
460	# cosAngle = angle(closest_H.coords, oxygen.coords, solute.coords, dimensions)
461	# arccosAngle = np.arccos(cosAngle)
462	# angleDegrees = np.degrees(arccosAngle)
463	# if math.isnan(angleDegrees) is False:
464	# angleDegrees = int(round(angleDegrees, 0))
465
466	if math.isnan(soluteOMangle) is False:	×
467	soluteOMangle = int(round(soluteOMangle, 0))	×
468	# else:
469	# print ('Error1: NaN for solute resid: %s, '\
470	# 'water resid %s and angle %s.'\
471	# % (nearest.resid, neighbour.resid, soluteOMangle))
472	# continue
473
474	if math.isnan(soluteOorthogAngle) is False:	×
475	soluteOorthogAngle = int(round(soluteOorthogAngle, 0))	×
476	# else:
477	# print ('Error1: NaN for solute resid: %s, '\
478	# 'water resid %s and angle %s.'\
479	# % (nearest.resid, neighbour.resid, soluteOorthogAngle))
480	# continue
481
482	return soluteOMangle, soluteOorthogAngle	×
483
484
NEW 485	def calcWaterSoluteAngle(solute, oxygen, H1, H2, dimensions):	×
486	"""
487	Calc angle in plane and orthog to plane between any water
488	and any solute atom.
489	"""
490
491	dist_solute_H1, dist_solute_H2 = None, None	×
492
493	for atomDist in solute.nearest_all_atom_array:	×
494	if atomDist[0] == H1.atom_num:	×
495	dist_solute_H1 = atomDist[1]	×
496	elif atomDist[0] == H2.atom_num:	×
497	dist_solute_H2 = atomDist[1]	×
498	else:
499	continue	×
500
NEW 501	if dist_solute_H1 is None:	×
502	dist_solute_H1 = distance(solute.coords, H1.coords, dimensions)	×
NEW 503	if dist_solute_H2 is None:	×
504	dist_solute_H2 = distance(solute.coords, H2.coords, dimensions)	×
505	closest_H = None	×
506	if dist_solute_H1 < dist_solute_H2:	×
507	closest_H = H1	×
508	elif dist_solute_H1 > dist_solute_H2:	×
509	closest_H = H2	×
510	else:
511	closest_H = H1	×
512
NEW 513	H1O_vector = vector(oxygen.coords, H1.coords, dimensions)	×
NEW 514	H2O_vector = vector(oxygen.coords, H2.coords, dimensions)	×
515	plane_normal = np.cross(H1O_vector, H2O_vector)	×
516	M_coords = getMcoords(H1.coords, H2.coords, dimensions)	×
NEW 517	OM_vector = vector(M_coords, oxygen.coords, dimensions)	×
518	plane_normal2 = np.cross(OM_vector, plane_normal)	×
NEW 519	OSolute_vector = vector(solute.coords, oxygen.coords, dimensions)	×
520	projected_OSolute = projectVectorToPlane(OSolute_vector, plane_normal)	×
NEW 521	projected_orthog_OSolute = projectVectorToPlane(OSolute_vector, plane_normal2)	×
522
523	# print ('OM', OM_vector, oxygen.atom_num)
524
525	soluteOMangle = angleBetweenVectors(OM_vector, projected_OSolute)	×
NEW 526	soluteOorthogAngle = angleBetweenVectors(OM_vector, projected_orthog_OSolute)	×
527
NEW 528	if (	×
529	np.dot(OM_vector, projected_OSolute) > 0
530	and np.dot(plane_normal2, projected_OSolute) > 0
531	):
532	soluteOMangle = (90 - soluteOMangle) + 270	×
533	# only go up to 180 degrees rather than 360 as above
534	# if soluteOMangle > 180:
535	# soluteOMangle =
536
NEW 537	if (	×
538	np.dot(-OM_vector, projected_OSolute) > 0
539	and np.dot(plane_normal2, projected_OSolute) > 0
540	):
541	soluteOMangle = (180 - soluteOMangle) + 180	×
542
NEW 543	if (	×
544	np.dot(-plane_normal, projected_orthog_OSolute) > 0
545	and np.dot(-OM_vector, projected_orthog_OSolute) > 0
546	):
547	soluteOorthogAngle = (180 - soluteOorthogAngle) + 180	×
548
NEW 549	if (	×
550	np.dot(-plane_normal, projected_orthog_OSolute) > 0
551	and np.dot(OM_vector, projected_orthog_OSolute) > 0
552	):
553	soluteOorthogAngle = (90 - soluteOorthogAngle) + 270	×
554
NEW 555	cosAngle = angle(closest_H.coords, oxygen.coords, solute.coords, dimensions)	×
556	arccosAngle = np.arccos(cosAngle)	×
557	angleDegrees = np.degrees(arccosAngle)	×
558	if math.isnan(angleDegrees) is False:	×
559	angleDegrees = int(round(angleDegrees, 0))	×
560	soluteOMangle = int(round(soluteOMangle, 0))	×
561	soluteOorthogAngle = int(round(soluteOorthogAngle, 0))	×
562	else:
NEW 563	logging.error(	×
564	"NaN for solute resid: %s, water resid %s and angle %s."
565	% (solute.resid, oxygen.resid, angleDegrees)
566	)
567
568	return soluteOMangle, soluteOorthogAngle	×

CCPBioSim / CodeEntropy / 13726590978

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous