13999062252

Committed 21 Mar 2025 07:03PM UTC coverage: 29.725%. Remained the same

Build # 13999062252

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #72 from CCPBioSim/63-conformational-entropy

Issue 63- Conformational Entropy Calculation. This merge closes #63 and resolves differences between the two versions of the code for aliphatic residues.

Run Details

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

2 existing lines in 1 file now uncovered.

162 of 545 relevant lines covered (29.72%)

0.89 hits per line

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

17.65

/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


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


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

1	# import MDAnalysis as mda
2	import numpy as np	3✔
3
4	from CodeEntropy import GeometricFunctions as GF	3✔
5
6
7	def select_levels(data_container, verbose):	3✔
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)	3✔
26	print("The number of molecules is {}.".format(number_molecules))	3✔
27	fragments = data_container.atoms.fragments	3✔
28	levels = [[] for _ in range(number_molecules)]	3✔
29
30	for molecule in range(number_molecules):	3✔
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:	3✔
46	print(levels)	×
47
48	return number_molecules, levels	3✔
49
50
51	def get_matrices(	3✔
52	data_container, level, verbose, start, end, step, number_frames, highest_level
53	):
54	"""
55	Function to create the force matrix needed for the transvibrational entropy
56	calculation and the torque matrix for the rovibrational entropy calculation.
57
58	Input
59	-----
60	data_container : MDAnalysis universe type with the information on the
61	molecule of interest.
62	level : string, which of the polymer, residue, or united atom levels
63	are the matrices for.
64	verbose : true/false, controlls how much is printed
65	start : int, starting frame, default 0 (first frame)
66	end : int, ending frame, default -1 (last frame)
67	step : int, step for going through trajectories, default 1
68
69	Returns
70	-------
71	force_matrix : force covariance matrix for transvibrational entropy
72	torque_matrix : torque convariance matrix for rovibrational entropy
73	"""
74
75	# Make beads
76	list_of_beads = GF.get_beads(data_container, level)	×
77
78	# number of beads and frames in trajectory
79	number_beads = len(list_of_beads)	×
80
81	# initialize force and torque arrays
82	weighted_forces = [[0 for x in range(number_frames)] for y in range(number_beads)]	×
83	weighted_torques = [[0 for x in range(number_frames)] for y in range(number_beads)]	×
84
85	# Calculate forces/torques for each bead
86	for bead_index in range(number_beads):	×
87	for timestep in data_container.trajectory[start:end:step]:	×
88	# Set up axes
89	# translation and rotation use different axes
90	# how the axes are defined depends on the level
91	trans_axes, rot_axes = GF.get_axes(data_container, level, bead_index)	×
92
93	# Sort out coordinates, forces, and torques for each atom in the bead
94	weighted_forces[bead_index][timestep.frame] = GF.get_weighted_forces(	×
95	data_container, list_of_beads[bead_index], trans_axes, highest_level
96	)
97	weighted_torques[bead_index][timestep.frame] = GF.get_weighted_torques(	×
98	data_container, list_of_beads[bead_index], rot_axes
99	)
100
101	# Make covariance matrices - looping over pairs of beads
102	# list of pairs of indices
103	pair_list = [(i, j) for i in range(number_beads) for j in range(number_beads)]	×
104
105	force_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]	×
106	torque_submatrix = [[0 for x in range(number_beads)] for y in range(number_beads)]	×
107
108	for i, j in pair_list:	×
109	# for each pair of beads
110	# reducing effort because the matrix for [i][j] is the transpose of the one for
111	# [j][i]
112	if i <= j:	×
113	# calculate the force covariance segment of the matrix
114	force_submatrix[i][j] = GF.create_submatrix(	×
115	weighted_forces[i], weighted_forces[j], number_frames
116	)
117	force_submatrix[j][i] = np.transpose(force_submatrix[i][j])	×
118
119	# calculate the torque covariance segment of the matrix
120	torque_submatrix[i][j] = GF.create_submatrix(	×
121	weighted_torques[i], weighted_torques[j], number_frames
122	)
123	torque_submatrix[j][i] = np.transpose(torque_submatrix[i][j])	×
124
125	# Tidy up
126	# use np.block to make submatrices into one matrix
127	force_matrix = np.block(	×
128	[
129	[force_submatrix[i][j] for j in range(number_beads)]
130	for i in range(number_beads)
131	]
132	)
133
134	torque_matrix = np.block(	×
135	[
136	[torque_submatrix[i][j] for j in range(number_beads)]
137	for i in range(number_beads)
138	]
139	)
140
141	# fliter zeros to remove any rows/columns that are all zero
142	force_matrix = GF.filter_zero_rows_columns(force_matrix, verbose)	×
143	torque_matrix = GF.filter_zero_rows_columns(torque_matrix, verbose)	×
144
145	if verbose:	×
146	with open("matrix.out", "a") as f:	×
147	print("force_matrix \n", file=f)	×
148	print(force_matrix, file=f)	×
149	print("torque_matrix \n", file=f)	×
150	print(torque_matrix, file=f)	×
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":	×
NEW 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	print("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	return dihedrals	×

CCPBioSim / CodeEntropy / 13999062252

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