19675730718

Committed 25 Nov 2025 03:56PM UTC coverage: 90.133% (-9.9%) from 100.0%

Build # 19675730718

Build Type

Pull #200

github

Committed by

web-flow

Commit Message

Merge 444701d78 into 51f6db64b

Pull Request Pull Request #200: Update Averaging over Groups of Molecules for Conformational Entropy

Coverage Stats

99 of 210 new or added lines in 5 files covered. (47.14%)

1014 of 1125 relevant lines covered (90.13%)

0.9 hits per line

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

29.41

/CodeEntropy/dihedral_tools.py

import logging

import numpy as np
from MDAnalysis.analysis.dihedrals import Dihedral
from rich.progress import (
    BarColumn,
    Progress,
    SpinnerColumn,
    TextColumn,
    TimeElapsedColumn,
)

logger = logging.getLogger(__name__)


class DihedralAnalysis:
    """
    Functions for finding dihedral angles and analysing them to get the
    states needed for the conformational entropy functions.
    """

    def __init__(self, universe_operations=None):
        """
        Initialise with placeholders.
        """
        self._universe_operations = universe_operations
        self.data_container = None
        self.states_ua = None
        self.states_res = None

    def build_conformational_states(
        self,
        data_container,
        levels,
        groups,
        start,
        end,
        step,
        bin_width,
    ):
        """
        Build the conformational states descriptors based on dihedral angles
        needed for the calculation of the conformational entropy.
        """
        number_groups = len(groups)
        states_ua = {}
        states_res = [None] * number_groups

        total_items = sum(
            len(levels[mol_id]) for mols in groups.values() for mol_id in mols
        )

        with Progress(
            SpinnerColumn(),
            TextColumn("[bold blue]{task.fields[title]}", justify="right"),
            BarColumn(),
            TextColumn("[progress.percentage]{task.percentage:>3.1f}%"),
            TimeElapsedColumn(),
        ) as progress:

            task = progress.add_task(
                "[green]Building Conformational States...",
                total=total_items,
                title="Starting...",
            )

        for group_id in groups.keys():
            molecules = groups[group_id]
            mol = self._universe_operations.get_molecule_container(data_container, 0)
            num_residues = len(mol.residues)
            dihedrals_ua = [[] for _ in range(num_residues)]
            peaks_ua = [{} for _ in range(num_residues)]
            dihedrals_res = []
            peaks_res = {}

            # Identify dihedral AtomGroups
            for level in levels[molecules[0]]:
                if level == "united_atom":
                    for res_id in range(num_residues):
                        selection1 = mol.residues[res_id].atoms.indices[0]
                        selection2 = mol.residues[res_id].atoms.indices[-1]
                        res_container = self._universe_operations.new_U_select_atom(
                            mol,
                            f"index {selection1}:" f"{selection2}",
                        )
                        heavy_res = self._universe_operations.new_U_select_atom(
                            res_container, "prop mass > 1.1"
                        )

                        dihedrals_ua[res_id] = self._get_dihedrals(heavy_res, level)

                elif level == "residue":
                    dihedrals_res = self._get_dihedrals(mol, level)

            # Identify peaks
            for level in levels[molecules[0]]:
                if level == "united_atom":
                    for res_id in range(num_residues):
                        if len(dihedrals_ua[res_id]) == 0:
                            # No dihedrals means no histogram or peaks
                            peaks_ua[res_id] = []
                        else:
                            peaks_ua[res_id] = self._identify_peaks(
                                data_container,
                                molecules,
                                dihedrals_ua[res_id],
                                bin_width,
                                start,
                                end,
                                step,
                            )

                elif level == "residue":
                    if len(dihedrals_res) == 0:
                        # No dihedrals means no histogram or peaks
                        peaks_res = []
                    else:
                        peaks_res = self._identify_peaks(
                            data_container,
                            molecules,
                            dihedrals_res,
                            bin_width,
                            start,
                            end,
                            step,
                        )

            # Assign states for each group
            for level in levels[molecules[0]]:
                if level == "united_atom":
                    for res_id in range(num_residues):
                        key = (group_id, res_id)
                        if len(dihedrals_ua[res_id]) == 0:
                            # No conformational states
                            states_ua[key] = []
                        else:
                            states_ua[key] = self._assign_states(
                                data_container,
                                molecules,
                                dihedrals_ua[res_id],
                                peaks_ua[res_id],
                                start,
                                end,
                                step,
                            )

                elif level == "residue":
                    if len(dihedrals_res) == 0:
                        # No conformational states
                        states_res[group_id] = []
                    else:
                        states_res[group_id] = self._assign_states(
                            data_container,
                            molecules,
                            dihedrals_res,
                            peaks_res,
                            start,
                            end,
                            step,
                        )

            progress.advance(task)

        return states_ua, states_res

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

        Args:
          data_container (MDAnalysis.Universe): system information
          level (str): level of the hierarchy (should be residue or polymer)

        Returns:
           dihedrals (array): set of dihedrals
        """
        # Start with empty array
        dihedrals = []
        atom_groups = []

        # if united atom level, read dihedrals from MDAnalysis universe
        if level == "united_atom":
            dihedrals = data_container.dihedrals
            num_dihedrals = len(dihedrals)
            for index in range(num_dihedrals):
                atom_groups.append(dihedrals[index].atoms)

        # if residue level, looking for dihedrals involving residues
        if level == "residue":
            num_residues = len(data_container.residues)
            logger.debug(f"Number Residues: {num_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
                    atom_groups.append(atom_group)

        logger.debug(f"Level: {level}, Dihedrals: {atom_groups}")

        return atom_groups

    def _identify_peaks(
        self,
        data_container,
        molecules,
        dihedrals,
        bin_width,
        start,
        end,
        step,
    ):
        """
        Build a histogram of the dihedral data and identify the peaks.
        This is to give the information needed for the adaptive method
        of identifying dihedral states.
        """
        peak_values = [] * len(dihedrals)
        for dihedral_index in range(len(dihedrals)):
            phi = []
            # get the values of the angle for the dihedral
            # loop over all molecules in the averaging group
            # dihedral angle values have a range from -180 to 180
            for molecule in molecules:
                mol = self._universe_operations.get_molecule_container(
                    data_container, molecule
                )
                number_frames = len(mol.trajectory)
                R = Dihedral(dihedrals).run()
                for timestep in range(number_frames):
                    value = R.results.angles[timestep][dihedral_index]

                    # We want postive values in range 0 to 360 to make
                    # the peak assignment.
                    # works using the fact that dihedrals have circular symetry
                    # (i.e. -15 degrees = +345 degrees)
                    if value < 0:
                        value += 360
                    phi.append(value)

            # create a histogram using numpy
            number_bins = int(360 / bin_width)
            popul, bin_edges = np.histogram(a=phi, bins=number_bins, range=(0, 360))
            bin_value = [
                0.5 * (bin_edges[i] + bin_edges[i + 1]) for i in range(0, len(popul))
            ]

            # identify "convex turning-points" and populate a list of peaks
            # peak : a bin whose neighboring bins have smaller population
            # NOTE might have problems if the peak is wide with a flat or
            # sawtooth top in which case check you have a sensible bin width

            peaks = []
            for bin_index in range(number_bins):
                # if there is no dihedrals in a bin then it cannot be a peak
                if popul[bin_index] == 0:
                    pass
                # being careful of the last bin
                # (dihedrals have circular symmetry, the histogram does not)
                elif (
                    bin_index == number_bins - 1
                ):  # the -1 is because the index starts with 0 not 1
                    if (
                        popul[bin_index] >= popul[bin_index - 1]
                        and popul[bin_index] >= popul[0]
                    ):
                        peaks.append(bin_value[bin_index])
                else:
                    if (
                        popul[bin_index] >= popul[bin_index - 1]
                        and popul[bin_index] >= popul[bin_index + 1]
                    ):
                        peaks.append(bin_value[bin_index])

            peak_values.append(peaks)

            logger.debug(f"Dihedral: {dihedral_index}, Peak Values: {peak_values}")

        return peak_values

    def _assign_states(
        self,
        data_container,
        molecules,
        dihedrals,
        peaks,
        start,
        end,
        step,
    ):
        """
        Turn the dihedral values into conformations based on the peaks
        from the histogram.
        Then combine these to form states for each molecule.
        """
        conformations = []
        states = []

        # get the values of the angle for the dihedral
        # dihedral angle values have a range from -180 to 180
        for molecule in molecules:
            mol = self._universe_operations.get_molecule_container(
                data_container, molecule
            )
            number_frames = len(mol.trajectory)
            R = Dihedral(dihedrals).run()
            for dihedral_index in range(len(dihedrals)):
                conformation = []
                for timestep in range(number_frames):
                    value = R.results.angles[timestep][dihedral_index]

                    # We want postive values in range 0 to 360 to make
                    # the peak assignment.
                    # works using the fact that dihedrals have circular symetry
                    # (i.e. -15 degrees = +345 degrees)
                    if value < 0:
                        value += 360

                    # Find the turning point/peak that the snapshot is closest to.
                    distances = [abs(value - peak) for peak in peaks[dihedral_index]]
                    conformation.append(np.argmin(distances))

                    logger.debug(
                        f"Dihedral: {dihedral_index} Conformations: {conformation}"
                    )
                conformations.append(conformation)

            # for all the dihedrals available concatenate the label of each
            # dihedral into the state for that frame
            mol_states = [
                state
                for state in (
                    "".join(
                        str(int(conformations[d][f])) for d in range(len(dihedrals))
                    )
                    for f in range(number_frames)
                )
                if state
            ]

            if states is None:
                states = mol_states
            else:
                states.extend(mol_states)

        logger.debug(f"States: {states}")

        return states

1	import logging	1✔
2
3	import numpy as np	1✔
4	from MDAnalysis.analysis.dihedrals import Dihedral	1✔
5	from rich.progress import (	1✔
6	BarColumn,
7	Progress,
8	SpinnerColumn,
9	TextColumn,
10	TimeElapsedColumn,
11	)
12
13	logger = logging.getLogger(__name__)	1✔
14
15
16	class DihedralAnalysis:	1✔
17	"""
18	Functions for finding dihedral angles and analysing them to get the
19	states needed for the conformational entropy functions.
20	"""
21
22	def __init__(self, universe_operations=None):	1✔
23	"""
24	Initialise with placeholders.
25	"""
26	self._universe_operations = universe_operations	1✔
27	self.data_container = None	1✔
28	self.states_ua = None	1✔
29	self.states_res = None	1✔
30
31	def build_conformational_states(	1✔
32	self,
33	data_container,
34	levels,
35	groups,
36	start,
37	end,
38	step,
39	bin_width,
40	):
41	"""
42	Build the conformational states descriptors based on dihedral angles
43	needed for the calculation of the conformational entropy.
44	"""
NEW 45	number_groups = len(groups)	×
NEW 46	states_ua = {}	×
NEW 47	states_res = [None] * number_groups	×
48
NEW 49	total_items = sum(	×
50	len(levels[mol_id]) for mols in groups.values() for mol_id in mols
51	)
52
NEW 53	with Progress(	×
54	SpinnerColumn(),
55	TextColumn("[bold blue]{task.fields[title]}", justify="right"),
56	BarColumn(),
57	TextColumn("[progress.percentage]{task.percentage:>3.1f}%"),
58	TimeElapsedColumn(),
59	) as progress:
60
NEW 61	task = progress.add_task(	×
62	"[green]Building Conformational States...",
63	total=total_items,
64	title="Starting...",
65	)
66
NEW 67	for group_id in groups.keys():	×
NEW 68	molecules = groups[group_id]	×
NEW 69	mol = self._universe_operations.get_molecule_container(data_container, 0)	×
NEW 70	num_residues = len(mol.residues)	×
NEW 71	dihedrals_ua = [[] for _ in range(num_residues)]	×
NEW 72	peaks_ua = [{} for _ in range(num_residues)]	×
NEW 73	dihedrals_res = []	×
NEW 74	peaks_res = {}	×
75
76	# Identify dihedral AtomGroups
NEW 77	for level in levels[molecules[0]]:	×
NEW 78	if level == "united_atom":	×
NEW 79	for res_id in range(num_residues):	×
NEW 80	selection1 = mol.residues[res_id].atoms.indices[0]	×
NEW 81	selection2 = mol.residues[res_id].atoms.indices[-1]	×
NEW 82	res_container = self._universe_operations.new_U_select_atom(	×
83	mol,
84	f"index {selection1}:" f"{selection2}",
85	)
NEW 86	heavy_res = self._universe_operations.new_U_select_atom(	×
87	res_container, "prop mass > 1.1"
88	)
89
NEW 90	dihedrals_ua[res_id] = self._get_dihedrals(heavy_res, level)	×
91
NEW 92	elif level == "residue":	×
NEW 93	dihedrals_res = self._get_dihedrals(mol, level)	×
94
95	# Identify peaks
NEW 96	for level in levels[molecules[0]]:	×
NEW 97	if level == "united_atom":	×
NEW 98	for res_id in range(num_residues):	×
NEW 99	if len(dihedrals_ua[res_id]) == 0:	×
100	# No dihedrals means no histogram or peaks
NEW 101	peaks_ua[res_id] = []	×
102	else:
NEW 103	peaks_ua[res_id] = self._identify_peaks(	×
104	data_container,
105	molecules,
106	dihedrals_ua[res_id],
107	bin_width,
108	start,
109	end,
110	step,
111	)
112
NEW 113	elif level == "residue":	×
NEW 114	if len(dihedrals_res) == 0:	×
115	# No dihedrals means no histogram or peaks
NEW 116	peaks_res = []	×
117	else:
NEW 118	peaks_res = self._identify_peaks(	×
119	data_container,
120	molecules,
121	dihedrals_res,
122	bin_width,
123	start,
124	end,
125	step,
126	)
127
128	# Assign states for each group
NEW 129	for level in levels[molecules[0]]:	×
NEW 130	if level == "united_atom":	×
NEW 131	for res_id in range(num_residues):	×
NEW 132	key = (group_id, res_id)	×
NEW 133	if len(dihedrals_ua[res_id]) == 0:	×
134	# No conformational states
NEW 135	states_ua[key] = []	×
136	else:
NEW 137	states_ua[key] = self._assign_states(	×
138	data_container,
139	molecules,
140	dihedrals_ua[res_id],
141	peaks_ua[res_id],
142	start,
143	end,
144	step,
145	)
146
NEW 147	elif level == "residue":	×
NEW 148	if len(dihedrals_res) == 0:	×
149	# No conformational states
NEW 150	states_res[group_id] = []	×
151	else:
NEW 152	states_res[group_id] = self._assign_states(	×
153	data_container,
154	molecules,
155	dihedrals_res,
156	peaks_res,
157	start,
158	end,
159	step,
160	)
161
NEW 162	progress.advance(task)	×
163
NEW 164	return states_ua, states_res	×
165
166	def _get_dihedrals(self, data_container, level):	1✔
167	"""
168	Define the set of dihedrals for use in the conformational entropy function.
169	If united atom level, the dihedrals are defined from the heavy atoms
170	(4 bonded atoms for 1 dihedral).
171	If residue level, use the bonds between residues to cast dihedrals.
172	Note: not using improper dihedrals only ones with 4 atoms/residues
173	in a linear arrangement.
174
175	Args:
176	data_container (MDAnalysis.Universe): system information
177	level (str): level of the hierarchy (should be residue or polymer)
178
179	Returns:
180	dihedrals (array): set of dihedrals
181	"""
182	# Start with empty array
183	dihedrals = []	1✔
184	atom_groups = []	1✔
185
186	# if united atom level, read dihedrals from MDAnalysis universe
187	if level == "united_atom":	1✔
188	dihedrals = data_container.dihedrals	1✔
189	num_dihedrals = len(dihedrals)	1✔
190	for index in range(num_dihedrals):	1✔
191	atom_groups.append(dihedrals[index].atoms)	1✔
192
193	# if residue level, looking for dihedrals involving residues
194	if level == "residue":	1✔
195	num_residues = len(data_container.residues)	1✔
196	logger.debug(f"Number Residues: {num_residues}")	1✔
197	if num_residues < 4:	1✔
198	logger.debug("no residue level dihedrals")	1✔
199
200	else:
201	# find bonds between residues N-3:N-2 and N-1:N
202	for residue in range(4, num_residues + 1):	1✔
203	# Using MDAnalysis selection,
204	# assuming only one covalent bond between neighbouring residues
205	# TODO not written for branched polymers
206	atom_string = (	1✔
207	"resindex "
208	+ str(residue - 4)
209	+ " and bonded resindex "
210	+ str(residue - 3)
211	)
212	atom1 = data_container.select_atoms(atom_string)	1✔
213
214	atom_string = (	1✔
215	"resindex "
216	+ str(residue - 3)
217	+ " and bonded resindex "
218	+ str(residue - 4)
219	)
220	atom2 = data_container.select_atoms(atom_string)	1✔
221
222	atom_string = (	1✔
223	"resindex "
224	+ str(residue - 2)
225	+ " and bonded resindex "
226	+ str(residue - 1)
227	)
228	atom3 = data_container.select_atoms(atom_string)	1✔
229
230	atom_string = (	1✔
231	"resindex "
232	+ str(residue - 1)
233	+ " and bonded resindex "
234	+ str(residue - 2)
235	)
236	atom4 = data_container.select_atoms(atom_string)	1✔
237
238	atom_group = atom1 + atom2 + atom3 + atom4	1✔
239	atom_groups.append(atom_group)	1✔
240
241	logger.debug(f"Level: {level}, Dihedrals: {atom_groups}")	1✔
242
243	return atom_groups	1✔
244
245	def _identify_peaks(	1✔
246	self,
247	data_container,
248	molecules,
249	dihedrals,
250	bin_width,
251	start,
252	end,
253	step,
254	):
255	"""
256	Build a histogram of the dihedral data and identify the peaks.
257	This is to give the information needed for the adaptive method
258	of identifying dihedral states.
259	"""
NEW 260	peak_values = [] * len(dihedrals)	×
NEW 261	for dihedral_index in range(len(dihedrals)):	×
NEW 262	phi = []	×
263	# get the values of the angle for the dihedral
264	# loop over all molecules in the averaging group
265	# dihedral angle values have a range from -180 to 180
NEW 266	for molecule in molecules:	×
NEW 267	mol = self._universe_operations.get_molecule_container(	×
268	data_container, molecule
269	)
NEW 270	number_frames = len(mol.trajectory)	×
NEW 271	R = Dihedral(dihedrals).run()	×
NEW 272	for timestep in range(number_frames):	×
NEW 273	value = R.results.angles[timestep][dihedral_index]	×
274
275	# We want postive values in range 0 to 360 to make
276	# the peak assignment.
277	# works using the fact that dihedrals have circular symetry
278	# (i.e. -15 degrees = +345 degrees)
NEW 279	if value < 0:	×
NEW 280	value += 360	×
NEW 281	phi.append(value)	×
282
283	# create a histogram using numpy
NEW 284	number_bins = int(360 / bin_width)	×
NEW 285	popul, bin_edges = np.histogram(a=phi, bins=number_bins, range=(0, 360))	×
NEW 286	bin_value = [	×
287	0.5 * (bin_edges[i] + bin_edges[i + 1]) for i in range(0, len(popul))
288	]
289
290	# identify "convex turning-points" and populate a list of peaks
291	# peak : a bin whose neighboring bins have smaller population
292	# NOTE might have problems if the peak is wide with a flat or
293	# sawtooth top in which case check you have a sensible bin width
294
NEW 295	peaks = []	×
NEW 296	for bin_index in range(number_bins):	×
297	# if there is no dihedrals in a bin then it cannot be a peak
NEW 298	if popul[bin_index] == 0:	×
NEW 299	pass	×
300	# being careful of the last bin
301	# (dihedrals have circular symmetry, the histogram does not)
NEW 302	elif (	×
303	bin_index == number_bins - 1
304	): # the -1 is because the index starts with 0 not 1
NEW 305	if (	×
306	popul[bin_index] >= popul[bin_index - 1]
307	and popul[bin_index] >= popul[0]
308	):
NEW 309	peaks.append(bin_value[bin_index])	×
310	else:
NEW 311	if (	×
312	popul[bin_index] >= popul[bin_index - 1]
313	and popul[bin_index] >= popul[bin_index + 1]
314	):
NEW 315	peaks.append(bin_value[bin_index])	×
316
NEW 317	peak_values.append(peaks)	×
318
NEW 319	logger.debug(f"Dihedral: {dihedral_index}, Peak Values: {peak_values}")	×
320
NEW 321	return peak_values	×
322
323	def _assign_states(	1✔
324	self,
325	data_container,
326	molecules,
327	dihedrals,
328	peaks,
329	start,
330	end,
331	step,
332	):
333	"""
334	Turn the dihedral values into conformations based on the peaks
335	from the histogram.
336	Then combine these to form states for each molecule.
337	"""
NEW 338	conformations = []	×
NEW 339	states = []	×
340
341	# get the values of the angle for the dihedral
342	# dihedral angle values have a range from -180 to 180
NEW 343	for molecule in molecules:	×
NEW 344	mol = self._universe_operations.get_molecule_container(	×
345	data_container, molecule
346	)
NEW 347	number_frames = len(mol.trajectory)	×
NEW 348	R = Dihedral(dihedrals).run()	×
NEW 349	for dihedral_index in range(len(dihedrals)):	×
NEW 350	conformation = []	×
NEW 351	for timestep in range(number_frames):	×
NEW 352	value = R.results.angles[timestep][dihedral_index]	×
353
354	# We want postive values in range 0 to 360 to make
355	# the peak assignment.
356	# works using the fact that dihedrals have circular symetry
357	# (i.e. -15 degrees = +345 degrees)
NEW 358	if value < 0:	×
NEW 359	value += 360	×
360
361	# Find the turning point/peak that the snapshot is closest to.
NEW 362	distances = [abs(value - peak) for peak in peaks[dihedral_index]]	×
NEW 363	conformation.append(np.argmin(distances))	×
364
NEW 365	logger.debug(	×
366	f"Dihedral: {dihedral_index} Conformations: {conformation}"
367	)
NEW 368	conformations.append(conformation)	×
369
370	# for all the dihedrals available concatenate the label of each
371	# dihedral into the state for that frame
NEW 372	mol_states = [	×
373	state
374	for state in (
375	"".join(
376	str(int(conformations[d][f])) for d in range(len(dihedrals))
377	)
378	for f in range(number_frames)
379	)
380	if state
381	]
382
NEW 383	if states is None:	×
NEW 384	states = mol_states	×
385	else:
NEW 386	states.extend(mol_states)	×
387
NEW 388	logger.debug(f"States: {states}")	×
389
NEW 390	return states	×

CCPBioSim / CodeEntropy / 19675730718

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