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

Run Details

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

# import MDAnalysis as mda
import numpy as np

from CodeEntropy import GeometricFunctions as GF


def select_levels(data_container, verbose):
    """
    Function to read input system and identify the number of molecules and the levels
    (i.e. united atom, residue and/or polymer) that should be used.
    The level refers to the size of the bead (atom or collection of atoms) that will be
    used in the entropy calculations.

    Input
    -----
       arg_DataContainer : MDAnalysis universe object containing the system of interest

    Returns
    -------
       number_molecules : integer
       levels : array of strings for each molecule
    """

    # fragments is MDAnalysis terminology for what chemists would call molecules
    number_molecules = len(data_container.atoms.fragments)
    print("The number of molecules is {}.".format(number_molecules))
    fragments = data_container.atoms.fragments
    levels = [[] for _ in range(number_molecules)]

    for molecule in range(number_molecules):
        levels[molecule].append("united_atom")  # every molecule has at least one atom

        atoms_in_fragment = fragments[molecule].select_atoms("not name H*")
        number_residues = len(atoms_in_fragment.residues)

        # if a fragment has more than one heavy atom assign residue level
        if len(atoms_in_fragment) > 1:
            levels[molecule].append("residue")

            # if assigned residue level and there is more than one residue assign
            # polymer level
            if number_residues > 1:
                levels[molecule].append("polymer")

    if verbose:
        print(levels)

    return number_molecules, levels


# END get_levels


def get_matrices(
    data_container, level, verbose, start, end, step, number_frames, highest_level
):
    """
    Function to create the force matrix needed for the transvibrational entropy
    calculation and the torque matrix for the rovibrational entropy calculation.

    Input
    -----
        data_container : MDAnalysis universe type with the information on the
        molecule of interest.
        level : string, which of the polymer, residue, or united atom levels
        are the matrices for.
        verbose : true/false, controlls how much is printed
        start : int, starting frame, default 0 (first frame)
        end : int, ending frame, default -1 (last frame)
        step : int, step for going through trajectories, default 1

    Returns
    -------
        force_matrix : force covariance matrix for transvibrational entropy
        torque_matrix : torque convariance matrix for rovibrational entropy
    """

    # Make beads
    list_of_beads = GF.get_beads(data_container, level)

    # number of beads and frames in trajectory
    number_beads = len(list_of_beads)

    # initialize force and torque arrays
    weighted_forces = [[0 for x in range(number_frames)] for y in range(number_beads)]
    weighted_torques = [[0 for x in range(number_frames)] for y in range(number_beads)]

    # Calculate forces/torques for each bead
    for bead_index in range(number_beads):
        for timestep in data_container.trajectory[start:end:step]:
            # Set up axes
            # translation and rotation use different axes
            # how the axes are defined depends on the level
            trans_axes, rot_axes = GF.get_axes(data_container, level, bead_index)

            # Sort out coordinates, forces, and torques for each atom in the bead
            weighted_forces[bead_index][timestep.frame] = GF.get_weighted_forces(
                data_container, list_of_beads[bead_index], trans_axes, highest_level
            )
            weighted_torques[bead_index][timestep.frame] = GF.get_weighted_torques(
                data_container, list_of_beads[bead_index], rot_axes
            )

    # Make covariance matrices - looping over pairs of beads
    # list of pairs of indices
    pair_list = [(i, j) for i in range(number_beads) for j in range(number_beads)]

    force_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]
    torque_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]

    for i, j in pair_list:
        # for each pair of beads
        # reducing effort because the matrix for [i][j] is the transpose of the one for
        # [j][i]
        if i <= j:
            # calculate the force covariance segment of the matrix
            force_submatrix[i][j] = GF.create_submatrix(
                weighted_forces[i], weighted_forces[j], number_frames
            )
            force_submatrix[j][i] = np.transpose(force_submatrix[i][j])

            # calculate the torque covariance segment of the matrix
            torque_submatrix[i][j] = GF.create_submatrix(
                weighted_torques[i], weighted_torques[j], number_frames
            )
            torque_submatrix[j][i] = np.transpose(torque_submatrix[i][j])

    # Tidy up
    # use np.block to make submatrices into one matrix
    force_matrix = np.block(
        [
            [force_submatrix[i][j] for j in range(number_beads)]
            for i in range(number_beads)
        ]
    )

    torque_matrix = np.block(
        [
            [torque_submatrix[i][j] for j in range(number_beads)]
            for i in range(number_beads)
        ]
    )

    # fliter zeros to remove any rows/columns that are all zero
    force_matrix = GF.filter_zero_rows_columns(force_matrix, verbose)
    torque_matrix = GF.filter_zero_rows_columns(torque_matrix, verbose)

    if verbose:
        with open("matrix.out", "a") as f:
            print("force_matrix \n", file=f)
            print(force_matrix, file=f)
            print("torque_matrix \n", file=f)
            print(torque_matrix, file=f)

    return force_matrix, torque_matrix


# END get_matrices


def get_dihedrals(data_container, level):
    """
    Define the set of dihedrals for use in the conformational entropy function.
    If residue level, the dihedrals are defined from the atoms
    (4 bonded atoms for 1 dihedral).
    If polymer level, use the bonds between residues to cast dihedrals.
    Note: not using improper dihedrals only ones with 4 atoms/residues
    in a linear arrangement.

    Input
    -----
    data_container : system information
    level : level of the hierarchy (should be residue or polymer)

    Output
    ------
    dihedrals : set of dihedrals
    """
    # Start with empty array
    dihedrals = []

    # if united atom level, read dihedrals from MDAnalysis universe
    if level == "united_atom":
        # only use dihedrals made of heavy atoms
        heavy_atom_group = data_container.select_atoms("not name H*")
        dihedrals = heavy_atom_group.dihedrals

    # if residue level, looking for dihedrals involving residues
    if level == "residue":
        num_residues = len(data_container.residues)
        if num_residues < 4:
            print("no residue level dihedrals")

        else:
            # find bonds between residues N-3:N-2 and N-1:N
            for residue in range(4, num_residues + 1):
                # Using MDAnalysis selection,
                # assuming only one covalent bond between neighbouring residues
                # TODO not written for branched polymers
                atom_string = (
                    "resindex "
                    + str(residue - 4)
                    + " and bonded resindex "
                    + str(residue - 3)
                )
                atom1 = data_container.select_atoms(atom_string)

                atom_string = (
                    "resindex "
                    + str(residue - 3)
                    + " and bonded resindex "
                    + str(residue - 4)
                )
                atom2 = data_container.select_atoms(atom_string)

                atom_string = (
                    "resindex "
                    + str(residue - 2)
                    + " and bonded resindex "
                    + str(residue - 1)
                )
                atom3 = data_container.select_atoms(atom_string)

                atom_string = (
                    "resindex "
                    + str(residue - 1)
                    + " and bonded resindex "
                    + str(residue - 2)
                )
                atom4 = data_container.select_atoms(atom_string)

                atom_group = atom1 + atom2 + atom3 + atom4
                dihedrals.append(atom_group.dihedral)

    return dihedrals


# END get_dihedrals

1	# import MDAnalysis as mda
NEW 2	import numpy as np	×
3
4	from CodeEntropy import GeometricFunctions as GF	×
5
6
7	def select_levels(data_container, verbose):	×
8	"""
9	Function to read input system and identify the number of molecules and the levels
10	(i.e. united atom, residue and/or polymer) that should be used.
11	The level refers to the size of the bead (atom or collection of atoms) that will be
12	used in the entropy calculations.
13
14	Input
15	-----
16	arg_DataContainer : MDAnalysis universe object containing the system of interest
17
18	Returns
19	-------
20	number_molecules : integer
21	levels : array of strings for each molecule
22	"""
23
24	# fragments is MDAnalysis terminology for what chemists would call molecules
25	number_molecules = len(data_container.atoms.fragments)	×
26	print("The number of molecules is {}.".format(number_molecules))	×
27	fragments = data_container.atoms.fragments	×
28	levels = [[] for _ in range(number_molecules)]	×
29
30	for molecule in range(number_molecules):	×
NEW 31	levels[molecule].append("united_atom") # every molecule has at least one atom	×
32
33	atoms_in_fragment = fragments[molecule].select_atoms("not name H*")	×
34	number_residues = len(atoms_in_fragment.residues)	×
35
36	# if a fragment has more than one heavy atom assign residue level
37	if len(atoms_in_fragment) > 1:	×
38	levels[molecule].append("residue")	×
39
40	# if assigned residue level and there is more than one residue assign
41	# polymer level
42	if number_residues > 1:	×
43	levels[molecule].append("polymer")	×
44
45	if verbose:	×
46	print(levels)	×
47
48	return number_molecules, levels	×
49
50
51	# END get_levels
52
53
NEW 54	def get_matrices(	×
55	data_container, level, verbose, start, end, step, number_frames, highest_level
56	):
57	"""
58	Function to create the force matrix needed for the transvibrational entropy
59	calculation and the torque matrix for the rovibrational entropy calculation.
60
61	Input
62	-----
63	data_container : MDAnalysis universe type with the information on the
64	molecule of interest.
65	level : string, which of the polymer, residue, or united atom levels
66	are the matrices for.
67	verbose : true/false, controlls how much is printed
68	start : int, starting frame, default 0 (first frame)
69	end : int, ending frame, default -1 (last frame)
70	step : int, step for going through trajectories, default 1
71
72	Returns
73	-------
74	force_matrix : force covariance matrix for transvibrational entropy
75	torque_matrix : torque convariance matrix for rovibrational entropy
76	"""
77
78	# Make beads
79	list_of_beads = GF.get_beads(data_container, level)	×
80
81	# number of beads and frames in trajectory
82	number_beads = len(list_of_beads)	×
83
84	# initialize force and torque arrays
85	weighted_forces = [[0 for x in range(number_frames)] for y in range(number_beads)]	×
86	weighted_torques = [[0 for x in range(number_frames)] for y in range(number_beads)]	×
87
88	# Calculate forces/torques for each bead
89	for bead_index in range(number_beads):	×
90	for timestep in data_container.trajectory[start:end:step]:	×
91	# Set up axes
92	# translation and rotation use different axes
93	# how the axes are defined depends on the level
94	trans_axes, rot_axes = GF.get_axes(data_container, level, bead_index)	×
95
96	# Sort out coordinates, forces, and torques for each atom in the bead
NEW 97	weighted_forces[bead_index][timestep.frame] = GF.get_weighted_forces(	×
98	data_container, list_of_beads[bead_index], trans_axes, highest_level
99	)
NEW 100	weighted_torques[bead_index][timestep.frame] = GF.get_weighted_torques(	×
101	data_container, list_of_beads[bead_index], rot_axes
102	)
103
104	# Make covariance matrices - looping over pairs of beads
105	# list of pairs of indices
NEW 106	pair_list = [(i, j) for i in range(number_beads) for j in range(number_beads)]	×
107
108	force_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]	×
109	torque_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]	×
110
111	for i, j in pair_list:	×
112	# for each pair of beads
113	# reducing effort because the matrix for [i][j] is the transpose of the one for
114	# [j][i]
115	if i <= j:	×
116	# calculate the force covariance segment of the matrix
NEW 117	force_submatrix[i][j] = GF.create_submatrix(	×
118	weighted_forces[i], weighted_forces[j], number_frames
119	)
NEW 120	force_submatrix[j][i] = np.transpose(force_submatrix[i][j])	×
121
122	# calculate the torque covariance segment of the matrix
NEW 123	torque_submatrix[i][j] = GF.create_submatrix(	×
124	weighted_torques[i], weighted_torques[j], number_frames
125	)
NEW 126	torque_submatrix[j][i] = np.transpose(torque_submatrix[i][j])	×
127
128	# Tidy up
129	# use np.block to make submatrices into one matrix
NEW 130	force_matrix = np.block(	×
131	[
132	[force_submatrix[i][j] for j in range(number_beads)]
133	for i in range(number_beads)
134	]
135	)
136
NEW 137	torque_matrix = np.block(	×
138	[
139	[torque_submatrix[i][j] for j in range(number_beads)]
140	for i in range(number_beads)
141	]
142	)
143
144	# fliter zeros to remove any rows/columns that are all zero
145	force_matrix = GF.filter_zero_rows_columns(force_matrix, verbose)	×
146	torque_matrix = GF.filter_zero_rows_columns(torque_matrix, verbose)	×
147
148	if verbose:	×
149	with open("matrix.out", "a") as f:	×
150	print("force_matrix \n", file=f)	×
151	print(force_matrix, file=f)	×
152	print("torque_matrix \n", file=f)	×
153	print(torque_matrix, file=f)	×
154
155	return force_matrix, torque_matrix	×
156
157
158	# END get_matrices
159
160
161	def get_dihedrals(data_container, level):	×
162	"""
163	Define the set of dihedrals for use in the conformational entropy function.
164	If residue level, the dihedrals are defined from the atoms
165	(4 bonded atoms for 1 dihedral).
166	If polymer level, use the bonds between residues to cast dihedrals.
167	Note: not using improper dihedrals only ones with 4 atoms/residues
168	in a linear arrangement.
169
170	Input
171	-----
172	data_container : system information
173	level : level of the hierarchy (should be residue or polymer)
174
175	Output
176	------
177	dihedrals : set of dihedrals
178	"""
179	# Start with empty array
180	dihedrals = []	×
181
182	# if united atom level, read dihedrals from MDAnalysis universe
183	if level == "united_atom":	×
184	# only use dihedrals made of heavy atoms
NEW 185	heavy_atom_group = data_container.select_atoms("not name H*")	×
186	dihedrals = heavy_atom_group.dihedrals	×
187
188	# if residue level, looking for dihedrals involving residues
189	if level == "residue":	×
190	num_residues = len(data_container.residues)	×
191	if num_residues < 4:	×
192	print("no residue level dihedrals")	×
193
194	else:
195	# find bonds between residues N-3:N-2 and N-1:N
NEW 196	for residue in range(4, num_residues + 1):	×
197	# Using MDAnalysis selection,
198	# assuming only one covalent bond between neighbouring residues
199	# TODO not written for branched polymers
NEW 200	atom_string = (	×
201	"resindex "
202	+ str(residue - 4)
203	+ " and bonded resindex "
204	+ str(residue - 3)
205	)
206	atom1 = data_container.select_atoms(atom_string)	×
207
NEW 208	atom_string = (	×
209	"resindex "
210	+ str(residue - 3)
211	+ " and bonded resindex "
212	+ str(residue - 4)
213	)
214	atom2 = data_container.select_atoms(atom_string)	×
215
NEW 216	atom_string = (	×
217	"resindex "
218	+ str(residue - 2)
219	+ " and bonded resindex "
220	+ str(residue - 1)
221	)
222	atom3 = data_container.select_atoms(atom_string)	×
223
NEW 224	atom_string = (	×
225	"resindex "
226	+ str(residue - 1)
227	+ " and bonded resindex "
228	+ str(residue - 2)
229	)
230	atom4 = data_container.select_atoms(atom_string)	×
231
232	atom_group = atom1 + atom2 + atom3 + atom4	×
233	dihedrals.append(atom_group.dihedral)	×
234
235	return dihedrals	×
236
237
238	# END get_dihedrals

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