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/GeometricFunctions.py

import numpy as np


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))

    return list_of_beads


# END


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)

    return trans_axes, rot_axes


# END


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

    return avg_position


# END


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

    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])

    return spherical_basis


# END


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()

    weighted_force = forces_trans / np.sqrt(mass)

    return weighted_force


# END


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.

    Input
    -----
    bead : the part of the molecule to be considered
    rot_axes : the axes relative to which the forces and coordinates are located
    frame : the frame number from the trajectory

    Output
    ------
    weighted_torque : 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):
        # cannot divide by zero
        if np.isclose(moment_of_inertia[dimension, dimension], 0):
            weighted_torque[dimension] = torques[dimension]
        else:
            weighted_torque[dimension] = torques[dimension] / np.sqrt(
                moment_of_inertia[dimension, dimension]
            )

    return weighted_torque


# END


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

    return submatrix


# END


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 verbose and init_shape != final_shape:
        print(
            "A shape change has occured ({},{}) -> ({}, {})".format(
                *init_shape, *final_shape
            )
        )

    return arg_matrix


# END

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

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