Skip to content

Commit e4a41f4

Browse files
jkalayanjkalayan
committed
uncomment Sconf calcs
1 parent 612b5d0 commit e4a41f4

2 files changed

Lines changed: 137 additions & 141 deletions

File tree

‎CodeEntropy/dihedral_tools.py‎

Lines changed: 122 additions & 124 deletions
Original file line numberDiff line numberDiff line change
@@ -2,14 +2,13 @@
22

33
import numpy as np
44
from MDAnalysis.analysis.dihedrals import Dihedral
5-
6-
# from rich.progress import (
7-
# BarColumn,
8-
# Progress,
9-
# SpinnerColumn,
10-
# TextColumn,
11-
# TimeElapsedColumn,
12-
# )
5+
from rich.progress import (
6+
BarColumn,
7+
Progress,
8+
SpinnerColumn,
9+
TextColumn,
10+
TimeElapsedColumn,
11+
)
1312

1413
logger = logging.getLogger(__name__)
1514

@@ -47,122 +46,121 @@ def build_conformational_states(
4746
states_ua = {}
4847
states_res = [None] * number_groups
4948

50-
# SWITCH OFF SCONF
51-
# total_items = sum(
52-
# len(levels[mol_id]) for mols in groups.values() for mol_id in mols
53-
# )
54-
# with Progress(
55-
# SpinnerColumn(),
56-
# TextColumn("[bold blue]{task.fields[title]}", justify="right"),
57-
# BarColumn(),
58-
# TextColumn("[progress.percentage]{task.percentage:>3.1f}%"),
59-
# TimeElapsedColumn(),
60-
# ) as progress:
61-
62-
# task = progress.add_task(
63-
# "[green]Building Conformational States...",
64-
# total=total_items,
65-
# title="Starting...",
66-
# )
67-
68-
# for group_id in groups.keys():
69-
# molecules = groups[group_id]
70-
# mol = self._universe_operations.get_molecule_container(
71-
# data_container, molecules[0]
72-
# )
73-
# num_residues = len(mol.residues)
74-
# dihedrals_ua = [[] for _ in range(num_residues)]
75-
# peaks_ua = [{} for _ in range(num_residues)]
76-
# dihedrals_res = []
77-
# peaks_res = {}
78-
79-
# # Identify dihedral AtomGroups
80-
# for level in levels[molecules[0]]:
81-
# if level == "united_atom":
82-
# for res_id in range(num_residues):
83-
# selection1 = mol.residues[res_id].atoms.indices[0]
84-
# selection2 = mol.residues[res_id].atoms.indices[-1]
85-
# res_container = self._universe_operations.new_U_select_atom(
86-
# mol,
87-
# f"index {selection1}:" f"{selection2}",
88-
# )
89-
# heavy_res = self._universe_operations.new_U_select_atom(
90-
# res_container, "prop mass > 1.1"
91-
# )
92-
93-
# dihedrals_ua[res_id] = self._get_dihedrals(heavy_res, level)
94-
95-
# elif level == "residue":
96-
# dihedrals_res = self._get_dihedrals(mol, level)
97-
98-
# # Identify peaks
99-
# for level in levels[molecules[0]]:
100-
# if level == "united_atom":
101-
# for res_id in range(num_residues):
102-
# if len(dihedrals_ua[res_id]) == 0:
103-
# # No dihedrals means no histogram or peaks
104-
# peaks_ua[res_id] = []
105-
# else:
106-
# peaks_ua[res_id] = self._identify_peaks(
107-
# data_container,
108-
# molecules,
109-
# dihedrals_ua[res_id],
110-
# bin_width,
111-
# start,
112-
# end,
113-
# step,
114-
# )
115-
116-
# elif level == "residue":
117-
# if len(dihedrals_res) == 0:
118-
# # No dihedrals means no histogram or peaks
119-
# peaks_res = []
120-
# else:
121-
# peaks_res = self._identify_peaks(
122-
# data_container,
123-
# molecules,
124-
# dihedrals_res,
125-
# bin_width,
126-
# start,
127-
# end,
128-
# step,
129-
# )
130-
131-
# # Assign states for each group
132-
# for level in levels[molecules[0]]:
133-
# if level == "united_atom":
134-
# for res_id in range(num_residues):
135-
# key = (group_id, res_id)
136-
# if len(dihedrals_ua[res_id]) == 0:
137-
# # No conformational states
138-
# states_ua[key] = []
139-
# else:
140-
# states_ua[key] = self._assign_states(
141-
# data_container,
142-
# molecules,
143-
# dihedrals_ua[res_id],
144-
# peaks_ua[res_id],
145-
# start,
146-
# end,
147-
# step,
148-
# )
149-
150-
# elif level == "residue":
151-
# if len(dihedrals_res) == 0:
152-
# # No conformational states
153-
# states_res[group_id] = []
154-
# else:
155-
# states_res[group_id] = self._assign_states(
156-
# data_container,
157-
# molecules,
158-
# dihedrals_res,
159-
# peaks_res,
160-
# start,
161-
# end,
162-
# step,
163-
# )
164-
165-
# progress.advance(task)
49+
total_items = sum(
50+
len(levels[mol_id]) for mols in groups.values() for mol_id in mols
51+
)
52+
with Progress(
53+
SpinnerColumn(),
54+
TextColumn("[bold blue]{task.fields[title]}", justify="right"),
55+
BarColumn(),
56+
TextColumn("[progress.percentage]{task.percentage:>3.1f}%"),
57+
TimeElapsedColumn(),
58+
) as progress:
59+
60+
task = progress.add_task(
61+
"[green]Building Conformational States...",
62+
total=total_items,
63+
title="Starting...",
64+
)
65+
66+
for group_id in groups.keys():
67+
molecules = groups[group_id]
68+
mol = self._universe_operations.get_molecule_container(
69+
data_container, molecules[0]
70+
)
71+
num_residues = len(mol.residues)
72+
dihedrals_ua = [[] for _ in range(num_residues)]
73+
peaks_ua = [{} for _ in range(num_residues)]
74+
dihedrals_res = []
75+
peaks_res = {}
76+
77+
# Identify dihedral AtomGroups
78+
for level in levels[molecules[0]]:
79+
if level == "united_atom":
80+
for res_id in range(num_residues):
81+
selection1 = mol.residues[res_id].atoms.indices[0]
82+
selection2 = mol.residues[res_id].atoms.indices[-1]
83+
res_container = self._universe_operations.new_U_select_atom(
84+
mol,
85+
f"index {selection1}:" f"{selection2}",
86+
)
87+
heavy_res = self._universe_operations.new_U_select_atom(
88+
res_container, "prop mass > 1.1"
89+
)
90+
91+
dihedrals_ua[res_id] = self._get_dihedrals(heavy_res, level)
92+
93+
elif level == "residue":
94+
dihedrals_res = self._get_dihedrals(mol, level)
95+
96+
# Identify peaks
97+
for level in levels[molecules[0]]:
98+
if level == "united_atom":
99+
for res_id in range(num_residues):
100+
if len(dihedrals_ua[res_id]) == 0:
101+
# No dihedrals means no histogram or peaks
102+
peaks_ua[res_id] = []
103+
else:
104+
peaks_ua[res_id] = self._identify_peaks(
105+
data_container,
106+
molecules,
107+
dihedrals_ua[res_id],
108+
bin_width,
109+
start,
110+
end,
111+
step,
112+
)
113+
114+
elif level == "residue":
115+
if len(dihedrals_res) == 0:
116+
# No dihedrals means no histogram or peaks
117+
peaks_res = []
118+
else:
119+
peaks_res = self._identify_peaks(
120+
data_container,
121+
molecules,
122+
dihedrals_res,
123+
bin_width,
124+
start,
125+
end,
126+
step,
127+
)
128+
129+
# Assign states for each group
130+
for level in levels[molecules[0]]:
131+
if level == "united_atom":
132+
for res_id in range(num_residues):
133+
key = (group_id, res_id)
134+
if len(dihedrals_ua[res_id]) == 0:
135+
# No conformational states
136+
states_ua[key] = []
137+
else:
138+
states_ua[key] = self._assign_states(
139+
data_container,
140+
molecules,
141+
dihedrals_ua[res_id],
142+
peaks_ua[res_id],
143+
start,
144+
end,
145+
step,
146+
)
147+
148+
elif level == "residue":
149+
if len(dihedrals_res) == 0:
150+
# No conformational states
151+
states_res[group_id] = []
152+
else:
153+
states_res[group_id] = self._assign_states(
154+
data_container,
155+
molecules,
156+
dihedrals_res,
157+
peaks_res,
158+
start,
159+
end,
160+
step,
161+
)
162+
163+
progress.advance(task)
166164

167165
return states_ua, states_res
168166

‎CodeEntropy/entropy.py‎

Lines changed: 15 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -474,23 +474,21 @@ def _process_united_atom_entropy(
474474
t_matrix, "torque", self._args.temperature, highest
475475
)
476476

477-
# SWITCH OFF SCONF
478-
# # Get the relevant conformational states
479-
# values = states[key]
480-
# # Check if there is information in the states array
481-
# contains_non_empty_states = (
482-
# np.any(values) if isinstance(values, np.ndarray) else any(values)
483-
# )
484-
485-
# # Calculate the conformational entropy
486-
# # If there are no conformational states (i.e. no dihedrals)
487-
# # then the conformational entropy is zero
488-
# S_conf_res = (
489-
# ce.conformational_entropy_calculation(values)
490-
# if contains_non_empty_states
491-
# else 0
492-
# )
493-
S_conf_res = 0
477+
# Get the relevant conformational states
478+
values = states[key]
479+
# Check if there is information in the states array
480+
contains_non_empty_states = (
481+
np.any(values) if isinstance(values, np.ndarray) else any(values)
482+
)
483+
484+
# Calculate the conformational entropy
485+
# If there are no conformational states (i.e. no dihedrals)
486+
# then the conformational entropy is zero
487+
S_conf_res = (
488+
ce.conformational_entropy_calculation(values)
489+
if contains_non_empty_states
490+
else 0
491+
)
494492

495493
# Add the data to the united atom level entropy
496494
S_trans += S_trans_res

0 commit comments

Comments
 (0)