14328256002

Committed 08 Apr 2025 08:00AM UTC coverage: 39.089% (-0.1%) from 39.213%

Build # 14328256002

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #82 from CCPBioSim/80-controlling-trajectory-size

Resolve issue #80 fixing indexing issue for start parameter affecting timesteps

Run Details

0 of 5 new or added lines in 2 files covered. (0.0%)

1 existing line in 1 file now uncovered.

249 of 637 relevant lines covered (39.09%)

1.17 hits per line

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

20.9

/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
            timestep_index = timestep.frame - start
            weighted_forces[bead_index][timestep_index] = GF.get_weighted_forces(
                data_container, list_of_beads[bead_index], trans_axes, highest_level
            )
            weighted_torques[bead_index][timestep_index] = 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
NEW 97	timestep_index = timestep.frame - start	×
NEW 98	weighted_forces[bead_index][timestep_index] = GF.get_weighted_forces(	×
99	data_container, list_of_beads[bead_index], trans_axes, highest_level
100	)
NEW 101	weighted_torques[bead_index][timestep_index] = GF.get_weighted_torques(	×
102	data_container, list_of_beads[bead_index], rot_axes
103	)
104
105	# Make covariance matrices - looping over pairs of beads
106	# list of pairs of indices
107	pair_list = [(i, j) for i in range(number_beads) for j in range(number_beads)]	×
108
109	force_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]	×
110	torque_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]	×
111
112	for i, j in pair_list:	×
113	# for each pair of beads
114	# reducing effort because the matrix for [i][j] is the transpose of the one for
115	# [j][i]
116	if i <= j:	×
117	# calculate the force covariance segment of the matrix
118	force_submatrix[i][j] = GF.create_submatrix(	×
119	weighted_forces[i], weighted_forces[j], number_frames
120	)
121	force_submatrix[j][i] = np.transpose(force_submatrix[i][j])	×
122
123	# calculate the torque covariance segment of the matrix
124	torque_submatrix[i][j] = GF.create_submatrix(	×
125	weighted_torques[i], weighted_torques[j], number_frames
126	)
127	torque_submatrix[j][i] = np.transpose(torque_submatrix[i][j])	×
128
129	# Tidy up
130	# use np.block to make submatrices into one matrix
131	force_matrix = np.block(	×
132	[
133	[force_submatrix[i][j] for j in range(number_beads)]
134	for i in range(number_beads)
135	]
136	)
137
138	torque_matrix = np.block(	×
139	[
140	[torque_submatrix[i][j] for j in range(number_beads)]
141	for i in range(number_beads)
142	]
143	)
144
145	# fliter zeros to remove any rows/columns that are all zero
146	force_matrix = GF.filter_zero_rows_columns(force_matrix, verbose)	×
147	torque_matrix = GF.filter_zero_rows_columns(torque_matrix, verbose)	×
148
149	logger.debug(f"Force Matrix: {force_matrix}")	×
150	logger.debug(f"Torque Matrix: {torque_matrix}")	×
151
152	return force_matrix, torque_matrix	×
153
154
155	def get_dihedrals(data_container, level):	3✔
156	"""
157	Define the set of dihedrals for use in the conformational entropy function.
158	If residue level, the dihedrals are defined from the atoms
159	(4 bonded atoms for 1 dihedral).
160	If polymer level, use the bonds between residues to cast dihedrals.
161	Note: not using improper dihedrals only ones with 4 atoms/residues
162	in a linear arrangement.
163
164	Input
165	-----
166	data_container : system information
167	level : level of the hierarchy (should be residue or polymer)
168
169	Output
170	------
171	dihedrals : set of dihedrals
172	"""
173	# Start with empty array
174	dihedrals = []	×
175
176	# if united atom level, read dihedrals from MDAnalysis universe
177	if level == "united_atom":	×
178	dihedrals = data_container.dihedrals	×
179
180	# if residue level, looking for dihedrals involving residues
181	if level == "residue":	×
182	num_residues = len(data_container.residues)	×
183	if num_residues < 4:	×
184	logger.debug("no residue level dihedrals")	×
185
186	else:
187	# find bonds between residues N-3:N-2 and N-1:N
188	for residue in range(4, num_residues + 1):	×
189	# Using MDAnalysis selection,
190	# assuming only one covalent bond between neighbouring residues
191	# TODO not written for branched polymers
192	atom_string = (	×
193	"resindex "
194	+ str(residue - 4)
195	+ " and bonded resindex "
196	+ str(residue - 3)
197	)
198	atom1 = data_container.select_atoms(atom_string)	×
199
200	atom_string = (	×
201	"resindex "
202	+ str(residue - 3)
203	+ " and bonded resindex "
204	+ str(residue - 4)
205	)
206	atom2 = data_container.select_atoms(atom_string)	×
207
208	atom_string = (	×
209	"resindex "
210	+ str(residue - 2)
211	+ " and bonded resindex "
212	+ str(residue - 1)
213	)
214	atom3 = data_container.select_atoms(atom_string)	×
215
216	atom_string = (	×
217	"resindex "
218	+ str(residue - 1)
219	+ " and bonded resindex "
220	+ str(residue - 2)
221	)
222	atom4 = data_container.select_atoms(atom_string)	×
223
224	atom_group = atom1 + atom2 + atom3 + atom4	×
225	dihedrals.append(atom_group.dihedral)	×
226
227	logger.debug(f"Dihedrals: {dihedrals}")	×
228
229	return dihedrals	×

CCPBioSim / CodeEntropy / 14328256002

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