14059657073

Committed 25 Mar 2025 12:31PM UTC coverage: 38.633% (+8.9%) from 29.725%

Build # 14059657073

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #73 from CCPBioSim/34-implement-python-logging

Implement Python Logging

Run Details

137 of 211 new or added lines in 9 files covered. (64.93%)

5 existing lines in 1 file now uncovered.

243 of 629 relevant lines covered (38.63%)

1.16 hits per line

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

7.64

/CodeEntropy/calculations/GeometricFunctions.py

import logging

import numpy as np

logger = logging.getLogger(__name__)


def get_beads(data_container, level):
    """
    Function to define beads depending on the level in the hierarchy.

    Input
    -----
    data_container : the MDAnalysis universe
    level : the heirarchy level (polymer, residue, or united atom)

    Output
    ------
    list_of_beads : the relevent beads
    """

    if level == "polymer":
        list_of_beads = []
        atom_group = "all"
        list_of_beads.append(data_container.select_atoms(atom_group))

    if level == "residue":
        list_of_beads = []
        num_residues = len(data_container.residues)
        for residue in range(num_residues):
            atom_group = "resindex " + str(residue)
            list_of_beads.append(data_container.select_atoms(atom_group))

    if level == "united_atom":
        list_of_beads = []
        heavy_atoms = data_container.select_atoms("not name H*")
        for atom in heavy_atoms:
            atom_group = (
                "index "
                + str(atom.index)
                + " or (name H* and bonded index "
                + str(atom.index)
                + ")"
            )
            list_of_beads.append(data_container.select_atoms(atom_group))

    logger.debug(f"List of beads: {list_of_beads}")

    return list_of_beads


def get_axes(data_container, level, index=0):
    """
    Function to set the translational and rotational axes.
    The translational axes are based on the principal axes of the unit one level larger
    than the level we are interested in (except for the polymer level where there is no
    larger unit). The rotational axes use the covalent links between residues or atoms
    where possible to define the axes, or if the unit is not bonded to others of the
    same level the prinicpal axes of the unit are used.

    Input
    -----
    data_container : the information about the molecule and trajectory
    level : the level (united atom, residue, or polymer) of interest
    index : residue index (integer)

    Output
    ------
    trans_axes : translational axes
    rot_axes : rotational axes
    """
    index = int(index)

    if level == "polymer":
        # for polymer use principle axis for both translation and rotation
        trans_axes = data_container.atoms.principal_axes()
        rot_axes = data_container.atoms.principal_axes()

    if level == "residue":
        # Translation
        # for residues use principal axes of whole molecule for translation
        trans_axes = data_container.atoms.principal_axes()

        # Rotation
        # find bonds between atoms in residue of interest and other residues
        # we are assuming bonds only exist between adjacent residues
        # (linear chains of residues)
        # TODO refine selection so that it will work for branched polymers
        index_prev = index - 1
        index_next = index + 1
        atom_set = data_container.select_atoms(
            f"(resindex {index_prev} or resindex {index_next}) and bonded resid {index}"
        )
        residue = data_container.select_atoms(f"resindex {index}")

        if len(atom_set) == 0:
            # if no bonds to other residues use pricipal axes of residue
            rot_axes = residue.atoms.principal_axes()

        else:
            # set center of rotation to center of mass of the residue
            center = residue.atoms.center_of_mass()

            # get vector for average position of bonded atoms
            vector = get_avg_pos(atom_set, center)

            # use spherical coordinates function to get rotational axes
            rot_axes = get_sphCoord_axes(vector)

    if level == "united_atom":
        # Translation
        # for united atoms use principal axes of residue for translation
        trans_axes = data_container.residues.principal_axes()

        # Rotation
        # for united atoms use heavy atoms bonded to the heavy atom
        atom_set = data_container.select_atoms(f"not name H* and bonded index {index}")

        # center at position of heavy atom
        atom_group = data_container.select_atoms(f"index {index}")
        center = atom_group.positions[0]

        # get vector for average position of hydrogens
        vector = get_avg_pos(atom_set, center)

        # use spherical coordinates function to get rotational axes
        rot_axes = get_sphCoord_axes(vector)

    logger.debug(f"Translational Axes: {trans_axes}")
    logger.debug(f"Rotational Axes: {rot_axes}")

    return trans_axes, rot_axes


def get_avg_pos(atom_set, center):
    """
    Function to get the average position of a set of atoms.

    Input
    -----
    atoms : MDAnalysis atom group
    center : position for center of rotation

    Output
    ------
    avg_position : three dimensional vector
    """
    # start with an empty vector
    avg_position = np.zeros((3))

    # get number of atoms
    number_atoms = len(atom_set.names)

    if number_atoms != 0:
        # sum positions for all atoms in the given set
        for atom_index in range(number_atoms):
            atom_position = atom_set.atoms[atom_index].position

            avg_position += atom_position

        avg_position /= number_atoms  # divide by number of atoms to get average

    else:
        # if no atoms in set the unit has no bonds to restrict its rotational motion,
        # so we can use a random vector to get the spherical coordinates axes
        avg_position = np.random.random(3)

    # transform the average position to a coordinate system with the origin at center
    avg_position = avg_position - center

    logger.debug(f"Average Position: {avg_position}")

    return avg_position


def get_sphCoord_axes(arg_r):
    """
    For a given vector in space, treat it is a radial vector rooted at 0,0,0 and
    derive a curvilinear coordinate system according to the rules of polar spherical
    coordinates
    """

    x2y2 = arg_r[0] ** 2 + arg_r[1] ** 2
    r2 = x2y2 + arg_r[2] ** 2

    # Check for division by zero
    if r2 == 0.0:
        raise ValueError("r2 is zero, cannot compute spherical coordinates.")

    if x2y2 == 0.0:
        raise ValueError("x2y2 is zero, cannot compute sin_phi and cos_phi.")

    # Check for non-negative values inside the square root
    if x2y2 / r2 < 0:
        raise ValueError(
            f"Negative value encountered for sin_theta calculation: {x2y2 / r2}. "
            f"Cannot take square root."
        )

    if x2y2 < 0:
        raise ValueError(
            f"Negative value encountered for sin_phi and cos_phi calculation: {x2y2}. "
            f"Cannot take square root."
        )

    if x2y2 != 0.0:
        sin_theta = np.sqrt(x2y2 / r2)
        cos_theta = arg_r[2] / np.sqrt(r2)

        sin_phi = arg_r[1] / np.sqrt(x2y2)
        cos_phi = arg_r[0] / np.sqrt(x2y2)

    else:
        sin_theta = 0.0
        cos_theta = 1

        sin_phi = 0.0
        cos_phi = 1

    # if abs(sin_theta) > 1 or abs(sin_phi) > 1:
    #     print('Bad sine : T {} , P {}'.format(sin_theta, sin_phi))

    # cos_theta = np.sqrt(1 - sin_theta*sin_theta)
    # cos_phi = np.sqrt(1 - sin_phi*sin_phi)

    # print('{} {} {}'.format(*arg_r))
    # print('Sin T : {}, cos T : {}'.format(sin_theta, cos_theta))
    # print('Sin P : {}, cos P : {}'.format(sin_phi, cos_phi))

    spherical_basis = np.zeros((3, 3))

    # r^
    spherical_basis[0, :] = np.asarray(
        [sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]
    )

    # Theta^
    spherical_basis[1, :] = np.asarray(
        [cos_theta * cos_phi, cos_theta * sin_phi, -sin_theta]
    )

    # Phi^
    spherical_basis[2, :] = np.asarray([-sin_phi, cos_phi, 0.0])

    logger.debug(f"Spherical Basis: {spherical_basis}")

    return spherical_basis


def get_weighted_forces(
    data_container, bead, trans_axes, highest_level, force_partitioning=0.5
):
    """
    Function to calculate the mass weighted forces for a given bead.

    Input
    -----
    bead : the part of the system to be considered
    trans_axes : the axes relative to which the forces are located

    Output
    ------
    weighted_force : the mass weighted sum of the forces in the bead
    """

    forces_trans = np.zeros((3,))

    # Sum forces from all atoms in the bead
    for atom in bead.atoms:
        # update local forces in translational axes
        forces_local = np.matmul(trans_axes, data_container.atoms[atom.index].force)
        forces_trans += forces_local

    if highest_level:
        # multiply by the force_partitioning parameter to avoid double counting
        # of the forces on weakly correlated atoms
        # the default value of force_partitioning is 0.5 (dividing by two)
        forces_trans = force_partitioning * forces_trans

    # divide the sum of forces by the mass of the bead to get the weighted forces
    mass = bead.total_mass()

    # Check that mass is positive to avoid division by 0 or negative values inside sqrt
    if mass <= 0:
        raise ValueError(
            f"Invalid mass value: {mass}. Mass must be positive to compute the square "
            f"root."
        )

    weighted_force = forces_trans / np.sqrt(mass)

    logger.debug(f"Weighted Force: {weighted_force}")

    return weighted_force


def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5):
    """
    Function to calculate the moment of inertia weighted torques for a given bead.

    This function computes torques in a rotated frame and then weights them using
    the moment of inertia tensor. To prevent numerical instability, it treats
    extremely small diagonal elements of the moment of inertia tensor as zero
    (since values below machine precision are effectively zero). This avoids
    unnecessary use of extended precision (e.g., float128).

    Additionally, if the computed torque is already zero, the function skips
    the division step, reducing unnecessary computations and potential errors.

    Parameters
    ----------
    data_container : object
        Contains atomic positions and forces.
    bead : object
        The part of the molecule to be considered.
    rot_axes : np.ndarray
        The axes relative to which the forces and coordinates are located.
    force_partitioning : float, optional
        Factor to adjust force contributions, default is 0.5.

    Returns
    -------
    np.ndarray
        The mass-weighted sum of the torques in the bead.
    """

    torques = np.zeros((3,))
    weighted_torque = np.zeros((3,))

    for atom in bead.atoms:

        # update local coordinates in rotational axes
        coords_rot = data_container.atoms[atom.index].position - bead.center_of_mass()
        coords_rot = np.matmul(rot_axes, coords_rot)
        # update local forces in rotational frame
        forces_rot = np.matmul(rot_axes, data_container.atoms[atom.index].force)

        # multiply by the force_partitioning parameter to avoid double counting
        # of the forces on weakly correlated atoms
        # the default value of force_partitioning is 0.5 (dividing by two)
        forces_rot = force_partitioning * forces_rot

        # define torques (cross product of coordinates and forces) in rotational axes
        torques_local = np.cross(coords_rot, forces_rot)
        torques += torques_local

    # divide by moment of inertia to get weighted torques
    # moment of inertia is a 3x3 tensor
    # the weighting is done in each dimension (x,y,z) using the diagonal elements of
    # the moment of inertia tensor
    moment_of_inertia = bead.moment_of_inertia()

    for dimension in range(3):
        # Skip calculation if torque is already zero
        if np.isclose(torques[dimension], 0):
            weighted_torque[dimension] = 0
            continue

        # Check for zero moment of inertia
        if np.isclose(moment_of_inertia[dimension, dimension], 0):
            raise ZeroDivisionError(
                f"Attempted to divide by zero moment of inertia in dimension "
                f"{dimension}."
            )

        # Check for negative moment of inertia
        if moment_of_inertia[dimension, dimension] < 0:
            raise ValueError(
                f"Negative value encountered for moment of inertia: "
                f"{moment_of_inertia[dimension, dimension]} "
                f"Cannot compute weighted torque."
            )

        # Compute weighted torque
        weighted_torque[dimension] = torques[dimension] / np.sqrt(
            moment_of_inertia[dimension, dimension]
        )

    logger.debug(f"Weighted Torque: {weighted_torque}")

    return weighted_torque


def create_submatrix(data_i, data_j, number_frames):
    """
    Function for making covariance matrices.

    Input
    -----
    data_i : values for bead i
    data_j : valuees for bead j

    Output
    ------
    submatrix : 3x3 matrix for the covariance between i and j
    """

    # Start with 3 by 3 matrix of zeros
    submatrix = np.zeros((3, 3))

    # For each frame calculate the outer product (cross product) of the data from the
    # two beads and add the result to the submatrix
    for frame in range(number_frames):
        outer_product_matrix = np.outer(data_i[frame], data_j[frame])
        submatrix = np.add(submatrix, outer_product_matrix)

    # Divide by the number of frames to get the average
    submatrix /= number_frames

    logger.debug(f"Submatrix: {submatrix}")

    return submatrix


def filter_zero_rows_columns(arg_matrix, verbose):
    """
    function for removing rows and columns that contain only zeros from a matrix

    Input
    -----
    arg_matrix : matrix

    Output
    ------
    arg_matrix : the reduced size matrix
    """

    # record the initial size
    init_shape = np.shape(arg_matrix)

    zero_indices = list(
        filter(
            lambda row: np.all(np.isclose(arg_matrix[row, :], 0.0)),
            np.arange(np.shape(arg_matrix)[0]),
        )
    )
    all_indices = np.ones((np.shape(arg_matrix)[0]), dtype=bool)
    all_indices[zero_indices] = False
    arg_matrix = arg_matrix[all_indices, :]

    all_indices = np.ones((np.shape(arg_matrix)[1]), dtype=bool)
    zero_indices = list(
        filter(
            lambda col: np.all(np.isclose(arg_matrix[:, col], 0.0)),
            np.arange(np.shape(arg_matrix)[1]),
        )
    )
    all_indices[zero_indices] = False
    arg_matrix = arg_matrix[:, all_indices]

    # get the final shape
    final_shape = np.shape(arg_matrix)

    if init_shape != final_shape:
        logger.debug(
            "A shape change has occurred ({},{}) -> ({}, {})".format(
                *init_shape, *final_shape
            )
        )

    logger.debug(f"arg_matrix: {arg_matrix}")

    return arg_matrix

1	import logging	3✔
2
3	import numpy as np	3✔
4
5	logger = logging.getLogger(__name__)	3✔
6
7
8	def get_beads(data_container, level):	3✔
9	"""
10	Function to define beads depending on the level in the hierarchy.
11
12	Input
13	-----
14	data_container : the MDAnalysis universe
15	level : the heirarchy level (polymer, residue, or united atom)
16
17	Output
18	------
19	list_of_beads : the relevent beads
20	"""
21
22	if level == "polymer":	×
23	list_of_beads = []	×
24	atom_group = "all"	×
25	list_of_beads.append(data_container.select_atoms(atom_group))	×
26
27	if level == "residue":	×
28	list_of_beads = []	×
29	num_residues = len(data_container.residues)	×
30	for residue in range(num_residues):	×
31	atom_group = "resindex " + str(residue)	×
32	list_of_beads.append(data_container.select_atoms(atom_group))	×
33
34	if level == "united_atom":	×
35	list_of_beads = []	×
36	heavy_atoms = data_container.select_atoms("not name H*")	×
37	for atom in heavy_atoms:	×
38	atom_group = (	×
39	"index "
40	+ str(atom.index)
41	+ " or (name H* and bonded index "
42	+ str(atom.index)
43	+ ")"
44	)
45	list_of_beads.append(data_container.select_atoms(atom_group))	×
46
NEW 47	logger.debug(f"List of beads: {list_of_beads}")	×
48
49	return list_of_beads	×
50
51
52	def get_axes(data_container, level, index=0):	3✔
53	"""
54	Function to set the translational and rotational axes.
55	The translational axes are based on the principal axes of the unit one level larger
56	than the level we are interested in (except for the polymer level where there is no
57	larger unit). The rotational axes use the covalent links between residues or atoms
58	where possible to define the axes, or if the unit is not bonded to others of the
59	same level the prinicpal axes of the unit are used.
60
61	Input
62	-----
63	data_container : the information about the molecule and trajectory
64	level : the level (united atom, residue, or polymer) of interest
65	index : residue index (integer)
66
67	Output
68	------
69	trans_axes : translational axes
70	rot_axes : rotational axes
71	"""
72	index = int(index)	×
73
74	if level == "polymer":	×
75	# for polymer use principle axis for both translation and rotation
76	trans_axes = data_container.atoms.principal_axes()	×
77	rot_axes = data_container.atoms.principal_axes()	×
78
79	if level == "residue":	×
80	# Translation
81	# for residues use principal axes of whole molecule for translation
82	trans_axes = data_container.atoms.principal_axes()	×
83
84	# Rotation
85	# find bonds between atoms in residue of interest and other residues
86	# we are assuming bonds only exist between adjacent residues
87	# (linear chains of residues)
88	# TODO refine selection so that it will work for branched polymers
89	index_prev = index - 1	×
90	index_next = index + 1	×
91	atom_set = data_container.select_atoms(	×
92	f"(resindex {index_prev} or resindex {index_next}) and bonded resid {index}"
93	)
94	residue = data_container.select_atoms(f"resindex {index}")	×
95
96	if len(atom_set) == 0:	×
97	# if no bonds to other residues use pricipal axes of residue
98	rot_axes = residue.atoms.principal_axes()	×
99
100	else:
101	# set center of rotation to center of mass of the residue
102	center = residue.atoms.center_of_mass()	×
103
104	# get vector for average position of bonded atoms
105	vector = get_avg_pos(atom_set, center)	×
106
107	# use spherical coordinates function to get rotational axes
108	rot_axes = get_sphCoord_axes(vector)	×
109
110	if level == "united_atom":	×
111	# Translation
112	# for united atoms use principal axes of residue for translation
113	trans_axes = data_container.residues.principal_axes()	×
114
115	# Rotation
116	# for united atoms use heavy atoms bonded to the heavy atom
117	atom_set = data_container.select_atoms(f"not name H* and bonded index {index}")	×
118
119	# center at position of heavy atom
120	atom_group = data_container.select_atoms(f"index {index}")	×
121	center = atom_group.positions[0]	×
122
123	# get vector for average position of hydrogens
124	vector = get_avg_pos(atom_set, center)	×
125
126	# use spherical coordinates function to get rotational axes
127	rot_axes = get_sphCoord_axes(vector)	×
128
NEW 129	logger.debug(f"Translational Axes: {trans_axes}")	×
NEW 130	logger.debug(f"Rotational Axes: {rot_axes}")	×
131
132	return trans_axes, rot_axes	×
133
134
135	def get_avg_pos(atom_set, center):	3✔
136	"""
137	Function to get the average position of a set of atoms.
138
139	Input
140	-----
141	atoms : MDAnalysis atom group
142	center : position for center of rotation
143
144	Output
145	------
146	avg_position : three dimensional vector
147	"""
148	# start with an empty vector
149	avg_position = np.zeros((3))	×
150
151	# get number of atoms
152	number_atoms = len(atom_set.names)	×
153
154	if number_atoms != 0:	×
155	# sum positions for all atoms in the given set
156	for atom_index in range(number_atoms):	×
157	atom_position = atom_set.atoms[atom_index].position	×
158
159	avg_position += atom_position	×
160
161	avg_position /= number_atoms # divide by number of atoms to get average	×
162
163	else:
164	# if no atoms in set the unit has no bonds to restrict its rotational motion,
165	# so we can use a random vector to get the spherical coordinates axes
166	avg_position = np.random.random(3)	×
167
168	# transform the average position to a coordinate system with the origin at center
169	avg_position = avg_position - center	×
170
NEW 171	logger.debug(f"Average Position: {avg_position}")	×
172
173	return avg_position	×
174
175
176	def get_sphCoord_axes(arg_r):	3✔
177	"""
178	For a given vector in space, treat it is a radial vector rooted at 0,0,0 and
179	derive a curvilinear coordinate system according to the rules of polar spherical
180	coordinates
181	"""
182
183	x2y2 = arg_r[0] 2 + arg_r[1] 2	×
184	r2 = x2y2 + arg_r[2] ** 2	×
185
186	# Check for division by zero
187	if r2 == 0.0:	×
188	raise ValueError("r2 is zero, cannot compute spherical coordinates.")	×
189
190	if x2y2 == 0.0:	×
191	raise ValueError("x2y2 is zero, cannot compute sin_phi and cos_phi.")	×
192
193	# Check for non-negative values inside the square root
194	if x2y2 / r2 < 0:	×
195	raise ValueError(	×
196	f"Negative value encountered for sin_theta calculation: {x2y2 / r2}. "
197	f"Cannot take square root."
198	)
199
200	if x2y2 < 0:	×
201	raise ValueError(	×
202	f"Negative value encountered for sin_phi and cos_phi calculation: {x2y2}. "
203	f"Cannot take square root."
204	)
205
206	if x2y2 != 0.0:	×
207	sin_theta = np.sqrt(x2y2 / r2)	×
208	cos_theta = arg_r[2] / np.sqrt(r2)	×
209
210	sin_phi = arg_r[1] / np.sqrt(x2y2)	×
211	cos_phi = arg_r[0] / np.sqrt(x2y2)	×
212
213	else:
214	sin_theta = 0.0	×
215	cos_theta = 1	×
216
217	sin_phi = 0.0	×
218	cos_phi = 1	×
219
220	# if abs(sin_theta) > 1 or abs(sin_phi) > 1:
221	# print('Bad sine : T {} , P {}'.format(sin_theta, sin_phi))
222
223	# cos_theta = np.sqrt(1 - sin_theta*sin_theta)
224	# cos_phi = np.sqrt(1 - sin_phi*sin_phi)
225
226	# print('{} {} {}'.format(*arg_r))
227	# print('Sin T : {}, cos T : {}'.format(sin_theta, cos_theta))
228	# print('Sin P : {}, cos P : {}'.format(sin_phi, cos_phi))
229
230	spherical_basis = np.zeros((3, 3))	×
231
232	# r^
233	spherical_basis[0, :] = np.asarray(	×
234	[sin_theta * cos_phi, sin_theta * sin_phi, cos_theta]
235	)
236
237	# Theta^
238	spherical_basis[1, :] = np.asarray(	×
239	[cos_theta * cos_phi, cos_theta * sin_phi, -sin_theta]
240	)
241
242	# Phi^
243	spherical_basis[2, :] = np.asarray([-sin_phi, cos_phi, 0.0])	×
244
NEW 245	logger.debug(f"Spherical Basis: {spherical_basis}")	×
246
247	return spherical_basis	×
248
249
250	def get_weighted_forces(	3✔
251	data_container, bead, trans_axes, highest_level, force_partitioning=0.5
252	):
253	"""
254	Function to calculate the mass weighted forces for a given bead.
255
256	Input
257	-----
258	bead : the part of the system to be considered
259	trans_axes : the axes relative to which the forces are located
260
261	Output
262	------
263	weighted_force : the mass weighted sum of the forces in the bead
264	"""
265
266	forces_trans = np.zeros((3,))	×
267
268	# Sum forces from all atoms in the bead
269	for atom in bead.atoms:	×
270	# update local forces in translational axes
271	forces_local = np.matmul(trans_axes, data_container.atoms[atom.index].force)	×
272	forces_trans += forces_local	×
273
274	if highest_level:	×
275	# multiply by the force_partitioning parameter to avoid double counting
276	# of the forces on weakly correlated atoms
277	# the default value of force_partitioning is 0.5 (dividing by two)
278	forces_trans = force_partitioning * forces_trans	×
279
280	# divide the sum of forces by the mass of the bead to get the weighted forces
281	mass = bead.total_mass()	×
282
283	# Check that mass is positive to avoid division by 0 or negative values inside sqrt
284	if mass <= 0:	×
285	raise ValueError(	×
286	f"Invalid mass value: {mass}. Mass must be positive to compute the square "
287	f"root."
288	)
289
290	weighted_force = forces_trans / np.sqrt(mass)	×
291
NEW 292	logger.debug(f"Weighted Force: {weighted_force}")	×
293
294	return weighted_force	×
295
296
297	def get_weighted_torques(data_container, bead, rot_axes, force_partitioning=0.5):	3✔
298	"""
299	Function to calculate the moment of inertia weighted torques for a given bead.
300
301	This function computes torques in a rotated frame and then weights them using
302	the moment of inertia tensor. To prevent numerical instability, it treats
303	extremely small diagonal elements of the moment of inertia tensor as zero
304	(since values below machine precision are effectively zero). This avoids
305	unnecessary use of extended precision (e.g., float128).
306
307	Additionally, if the computed torque is already zero, the function skips
308	the division step, reducing unnecessary computations and potential errors.
309
310	Parameters
311	----------
312	data_container : object
313	Contains atomic positions and forces.
314	bead : object
315	The part of the molecule to be considered.
316	rot_axes : np.ndarray
317	The axes relative to which the forces and coordinates are located.
318	force_partitioning : float, optional
319	Factor to adjust force contributions, default is 0.5.
320
321	Returns
322	-------
323	np.ndarray
324	The mass-weighted sum of the torques in the bead.
325	"""
326
327	torques = np.zeros((3,))	×
328	weighted_torque = np.zeros((3,))	×
329
330	for atom in bead.atoms:	×
331
332	# update local coordinates in rotational axes
333	coords_rot = data_container.atoms[atom.index].position - bead.center_of_mass()	×
334	coords_rot = np.matmul(rot_axes, coords_rot)	×
335	# update local forces in rotational frame
336	forces_rot = np.matmul(rot_axes, data_container.atoms[atom.index].force)	×
337
338	# multiply by the force_partitioning parameter to avoid double counting
339	# of the forces on weakly correlated atoms
340	# the default value of force_partitioning is 0.5 (dividing by two)
341	forces_rot = force_partitioning * forces_rot	×
342
343	# define torques (cross product of coordinates and forces) in rotational axes
344	torques_local = np.cross(coords_rot, forces_rot)	×
345	torques += torques_local	×
346
347	# divide by moment of inertia to get weighted torques
348	# moment of inertia is a 3x3 tensor
349	# the weighting is done in each dimension (x,y,z) using the diagonal elements of
350	# the moment of inertia tensor
351	moment_of_inertia = bead.moment_of_inertia()	×
352
353	for dimension in range(3):	×
354	# Skip calculation if torque is already zero
355	if np.isclose(torques[dimension], 0):	×
356	weighted_torque[dimension] = 0	×
357	continue	×
358
359	# Check for zero moment of inertia
360	if np.isclose(moment_of_inertia[dimension, dimension], 0):	×
361	raise ZeroDivisionError(	×
362	f"Attempted to divide by zero moment of inertia in dimension "
363	f"{dimension}."
364	)
365
366	# Check for negative moment of inertia
367	if moment_of_inertia[dimension, dimension] < 0:	×
368	raise ValueError(	×
369	f"Negative value encountered for moment of inertia: "
370	f"{moment_of_inertia[dimension, dimension]} "
371	f"Cannot compute weighted torque."
372	)
373
374	# Compute weighted torque
375	weighted_torque[dimension] = torques[dimension] / np.sqrt(	×
376	moment_of_inertia[dimension, dimension]
377	)
378
NEW 379	logger.debug(f"Weighted Torque: {weighted_torque}")	×
380
381	return weighted_torque	×
382
383
384	def create_submatrix(data_i, data_j, number_frames):	3✔
385	"""
386	Function for making covariance matrices.
387
388	Input
389	-----
390	data_i : values for bead i
391	data_j : valuees for bead j
392
393	Output
394	------
395	submatrix : 3x3 matrix for the covariance between i and j
396	"""
397
398	# Start with 3 by 3 matrix of zeros
399	submatrix = np.zeros((3, 3))	×
400
401	# For each frame calculate the outer product (cross product) of the data from the
402	# two beads and add the result to the submatrix
403	for frame in range(number_frames):	×
404	outer_product_matrix = np.outer(data_i[frame], data_j[frame])	×
405	submatrix = np.add(submatrix, outer_product_matrix)	×
406
407	# Divide by the number of frames to get the average
408	submatrix /= number_frames	×
409
NEW 410	logger.debug(f"Submatrix: {submatrix}")	×
411
412	return submatrix	×
413
414
415	def filter_zero_rows_columns(arg_matrix, verbose):	3✔
416	"""
417	function for removing rows and columns that contain only zeros from a matrix
418
419	Input
420	-----
421	arg_matrix : matrix
422
423	Output
424	------
425	arg_matrix : the reduced size matrix
426	"""
427
428	# record the initial size
429	init_shape = np.shape(arg_matrix)	×
430
431	zero_indices = list(	×
432	filter(
433	lambda row: np.all(np.isclose(arg_matrix[row, :], 0.0)),
434	np.arange(np.shape(arg_matrix)[0]),
435	)
436	)
437	all_indices = np.ones((np.shape(arg_matrix)[0]), dtype=bool)	×
438	all_indices[zero_indices] = False	×
439	arg_matrix = arg_matrix[all_indices, :]	×
440
441	all_indices = np.ones((np.shape(arg_matrix)[1]), dtype=bool)	×
442	zero_indices = list(	×
443	filter(
444	lambda col: np.all(np.isclose(arg_matrix[:, col], 0.0)),
445	np.arange(np.shape(arg_matrix)[1]),
446	)
447	)
448	all_indices[zero_indices] = False	×
449	arg_matrix = arg_matrix[:, all_indices]	×
450
451	# get the final shape
452	final_shape = np.shape(arg_matrix)	×
453
NEW 454	if init_shape != final_shape:	×
NEW 455	logger.debug(	×
456	"A shape change has occurred ({},{}) -> ({}, {})".format(
457	init_shape, final_shape
458	)
459	)
460
NEW 461	logger.debug(f"arg_matrix: {arg_matrix}")	×
462
463	return arg_matrix	×

CCPBioSim / CodeEntropy / 14059657073

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