13769110208

Committed 10 Mar 2025 03:57PM UTC coverage: 1.774% (+0.05%) from 1.726%

Build # 13769110208

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #53 from CCPBioSim/7-nan-error-handling

Add error handling for NaN and invalid values in calculations:

Run Details

2 of 33 new or added lines in 3 files covered. (6.06%)

3 existing lines in 2 files now uncovered.

55 of 3101 relevant lines covered (1.77%)

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)

    if np.any(delta < 0):
        negative_indices = np.where(delta < 0)
        negative_values = delta[negative_indices]
        raise ValueError(
            f"Negative values encountered in 'delta' at indices {negative_indices} "
            f"with values {negative_values}. "
            f"Cannot take square root of negative values."
        )

    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)

    if np.any(ba < 0):
        negative_indices = np.where(ba < 0)
        negative_values = ba[negative_indices]
        raise ValueError(
            f"Negative values encountered in 'ba' at indices {negative_indices} "
            f"with values {negative_values}. "
            f"Cannot take square root of negative values."
        )

    dist_ba = np.sqrt((ba**2).sum(axis=-1))

    if np.any(bc < 0):
        negative_indices = np.where(bc < 0)
        negative_values = bc[negative_indices]
        raise ValueError(
            f"Negative values encountered in 'bc' at indices {negative_indices} "
            f"with values {negative_values}. "
            f"Cannot take square root of negative values."
        )

    dist_bc = np.sqrt((bc**2).sum(axis=-1))

    if np.any(ac < 0):
        negative_indices = np.where(ac < 0)
        negative_values = ac[negative_indices]
        raise ValueError(
            f"Negative values encountered in 'ac' at indices {negative_indices} "
            f"with values {negative_values}. "
            f"Cannot take square root of negative values."
        )

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

CCPBioSim / CodeEntropy / 13769110208

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