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

21.21

/CodeEntropy/calculations/LevelFunctions.py

# import MDAnalysis as mda
import logging

import numpy as np

from CodeEntropy.calculations import GeometricFunctions as GF

logger = logging.getLogger(__name__)


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)
    logger.debug("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")

    logger.debug(f"levels {levels}")

    return number_molecules, 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)

    logger.debug(f"Force Matrix: {force_matrix}")
    logger.debug(f"Torque Matrix: {torque_matrix}")

    return force_matrix, torque_matrix


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":
        dihedrals = data_container.dihedrals

    # if residue level, looking for dihedrals involving residues
    if level == "residue":
        num_residues = len(data_container.residues)
        if num_residues < 4:
            logger.debug("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)

    logger.debug(f"Dihedrals: {dihedrals}")

    return dihedrals

1	# import MDAnalysis as mda
2	import logging	3✔
3
4	import numpy as np	3✔
5
6	from CodeEntropy.calculations import GeometricFunctions as GF	3✔
7
8	logger = logging.getLogger(__name__)	3✔
9
10
11	def select_levels(data_container, verbose):	3✔
12	"""
13	Function to read input system and identify the number of molecules and the levels
14	(i.e. united atom, residue and/or polymer) that should be used.
15	The level refers to the size of the bead (atom or collection of atoms) that will be
16	used in the entropy calculations.
17
18	Input
19	-----
20	arg_DataContainer : MDAnalysis universe object containing the system of interest
21
22	Returns
23	-------
24	number_molecules : integer
25	levels : array of strings for each molecule
26	"""
27
28	# fragments is MDAnalysis terminology for what chemists would call molecules
29	number_molecules = len(data_container.atoms.fragments)	3✔
30	logger.debug("The number of molecules is {}.".format(number_molecules))	3✔
31	fragments = data_container.atoms.fragments	3✔
32	levels = [[] for _ in range(number_molecules)]	3✔
33
34	for molecule in range(number_molecules):	3✔
35	levels[molecule].append("united_atom") # every molecule has at least one atom	×
36
37	atoms_in_fragment = fragments[molecule].select_atoms("not name H*")	×
38	number_residues = len(atoms_in_fragment.residues)	×
39
40	# if a fragment has more than one heavy atom assign residue level
41	if len(atoms_in_fragment) > 1:	×
42	levels[molecule].append("residue")	×
43
44	# if assigned residue level and there is more than one residue assign
45	# polymer level
46	if number_residues > 1:	×
47	levels[molecule].append("polymer")	×
48
49	logger.debug(f"levels {levels}")	3✔
50
51	return number_molecules, levels	3✔
52
53
54	def get_matrices(	3✔
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
97	weighted_forces[bead_index][timestep.frame] = GF.get_weighted_forces(	×
98	data_container, list_of_beads[bead_index], trans_axes, highest_level
99	)
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
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
117	force_submatrix[i][j] = GF.create_submatrix(	×
118	weighted_forces[i], weighted_forces[j], number_frames
119	)
120	force_submatrix[j][i] = np.transpose(force_submatrix[i][j])	×
121
122	# calculate the torque covariance segment of the matrix
123	torque_submatrix[i][j] = GF.create_submatrix(	×
124	weighted_torques[i], weighted_torques[j], number_frames
125	)
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
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
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
NEW 148	logger.debug(f"Force Matrix: {force_matrix}")	×
NEW 149	logger.debug(f"Torque Matrix: {torque_matrix}")	×
150
151	return force_matrix, torque_matrix	×
152
153
154	def get_dihedrals(data_container, level):	3✔
155	"""
156	Define the set of dihedrals for use in the conformational entropy function.
157	If residue level, the dihedrals are defined from the atoms
158	(4 bonded atoms for 1 dihedral).
159	If polymer level, use the bonds between residues to cast dihedrals.
160	Note: not using improper dihedrals only ones with 4 atoms/residues
161	in a linear arrangement.
162
163	Input
164	-----
165	data_container : system information
166	level : level of the hierarchy (should be residue or polymer)
167
168	Output
169	------
170	dihedrals : set of dihedrals
171	"""
172	# Start with empty array
173	dihedrals = []	×
174
175	# if united atom level, read dihedrals from MDAnalysis universe
176	if level == "united_atom":	×
177	dihedrals = data_container.dihedrals	×
178
179	# if residue level, looking for dihedrals involving residues
180	if level == "residue":	×
181	num_residues = len(data_container.residues)	×
182	if num_residues < 4:	×
NEW 183	logger.debug("no residue level dihedrals")	×
184
185	else:
186	# find bonds between residues N-3:N-2 and N-1:N
187	for residue in range(4, num_residues + 1):	×
188	# Using MDAnalysis selection,
189	# assuming only one covalent bond between neighbouring residues
190	# TODO not written for branched polymers
191	atom_string = (	×
192	"resindex "
193	+ str(residue - 4)
194	+ " and bonded resindex "
195	+ str(residue - 3)
196	)
197	atom1 = data_container.select_atoms(atom_string)	×
198
199	atom_string = (	×
200	"resindex "
201	+ str(residue - 3)
202	+ " and bonded resindex "
203	+ str(residue - 4)
204	)
205	atom2 = data_container.select_atoms(atom_string)	×
206
207	atom_string = (	×
208	"resindex "
209	+ str(residue - 2)
210	+ " and bonded resindex "
211	+ str(residue - 1)
212	)
213	atom3 = data_container.select_atoms(atom_string)	×
214
215	atom_string = (	×
216	"resindex "
217	+ str(residue - 1)
218	+ " and bonded resindex "
219	+ str(residue - 2)
220	)
221	atom4 = data_container.select_atoms(atom_string)	×
222
223	atom_group = atom1 + atom2 + atom3 + atom4	×
224	dihedrals.append(atom_group.dihedral)	×
225
NEW 226	logger.debug(f"Dihedrals: {dihedrals}")	×
227
228	return dihedrals	×

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