14928438681

Committed 09 May 2025 12:01PM UTC coverage: 40.52% (+0.03%) from 40.491%

Build # 14928438681

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #87 from CCPBioSim/86-indexing-error-in-creating-dihedral-histograms

Fixing an off by one indexing error when slicing trajectories.

Coverage Stats

1 of 4 new or added lines in 1 file covered. (25.0%)

265 of 654 relevant lines covered (40.52%)

1.22 hits per line

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

42.51

/CodeEntropy/main_mcc.py

import logging
import math
import os
import re
import sys

import MDAnalysis as mda
import pandas as pd

from CodeEntropy.calculations import EntropyFunctions as EF
from CodeEntropy.calculations import LevelFunctions as LF
from CodeEntropy.calculations import MDAUniverseHelper as MDAHelper
from CodeEntropy.config.arg_config_manager import ConfigManager
from CodeEntropy.config.data_logger import DataLogger
from CodeEntropy.config.logging_config import LoggingConfig


def create_job_folder():
    """
    Create a new job folder with an incremented job number based on existing folders.
    """
    # Get the current working directory
    base_dir = os.getcwd()

    # List all folders in the base directory
    existing_folders = [
        f for f in os.listdir(base_dir) if os.path.isdir(os.path.join(base_dir, f))
    ]

    # Filter folders that match the pattern 'jobXXX'
    job_folders = [f for f in existing_folders if re.match(r"job\d{3}", f)]

    # Determine the next job number
    if job_folders:
        max_job_number = max([int(re.search(r"\d{3}", f).group()) for f in job_folders])
        next_job_number = max_job_number + 1
    else:
        next_job_number = 1

    # Format the new job folder name
    new_job_folder = f"job{next_job_number:03d}"
    new_job_folder_path = os.path.join(base_dir, new_job_folder)

    # Create the new job folder
    os.makedirs(new_job_folder_path, exist_ok=True)

    return new_job_folder_path


def main():
    """
    Main function for calculating the entropy of a system using the multiscale cell
    correlation method.
    """
    folder = create_job_folder()
    data_logger = DataLogger()
    arg_config = ConfigManager()

    # Load configuration
    config = arg_config.load_config("config.yaml")
    if config is None:
        raise ValueError(
            "No configuration file found, and no CLI arguments were provided."
        )

    parser = arg_config.setup_argparse()
    args, unknown = parser.parse_known_args()
    args.output_file = os.path.join(folder, args.output_file)

    try:
        # Initialize the logging system once
        logging_config = LoggingConfig(folder)
        logger = logging_config.setup_logging()

        # Process each run in the YAML configuration
        for run_name, run_config in config.items():
            if isinstance(run_config, dict):
                # Merging CLI arguments with YAML configuration
                args = arg_config.merge_configs(args, run_config)

                # Determine logging level
                log_level = logging.DEBUG if args.verbose else logging.INFO

                # Update the logging level
                logging_config.update_logging_level(log_level)

                # Capture and log the command-line invocation
                command = " ".join(sys.argv)
                logging.getLogger("commands").info(command)

                # Ensure necessary arguments are provided
                if not getattr(args, "top_traj_file"):
                    raise ValueError(
                        "The 'top_traj_file' argument is required but not provided."
                    )
                if not getattr(args, "selection_string"):
                    raise ValueError(
                        "The 'selection_string' argument is required but not provided."
                    )

                # Log all inputs for the current run
                logger.info(f"All input for {run_name}")
                for arg in vars(args):
                    logger.info(f" {arg}: {getattr(args, arg) or ''}")
            else:
                logger.warning(f"Run configuration for {run_name} is not a dictionary.")
    except ValueError as e:
        logger.error(e)
        raise

    # Get topology and trajectory file names and make universe
    tprfile = args.top_traj_file[0]
    trrfile = args.top_traj_file[1:]
    u = mda.Universe(tprfile, trrfile)

    # Define bin_width for histogram from inputs
    bin_width = args.bin_width

    # Define trajectory slicing from inputs
    start = args.start
    if start is None:
        start = 0
    end = args.end
    if end is None:
        end = -1
    step = args.step
    if step is None:
        step = 1
    # Count number of frames, easy if not slicing
    # MDAnalysis trajectory slicing only includes up to end-1
    # This works the way we want it to if the whole trajectory is being included
    if start == 0 and end == -1 and step == 1:
        end = len(u.trajectory)
        number_frames = len(u.trajectory)
    elif end == -1:
        end = len(u.trajectory)
        number_frames = math.floor((end - start) / step)
    else:
        end = end + 1
        number_frames = math.floor((end - start) / step)
    logger.debug(f"Number of Frames: {number_frames}")

    # Create pandas data frame for results
    results_df = pd.DataFrame(columns=["Molecule ID", "Level", "Type", "Result"])
    residue_results_df = pd.DataFrame(
        columns=["Molecule ID", "Residue", "Type", "Result"]
    )

    # Reduce number of atoms in MDA universe to selection_string arg
    # (default all atoms included)
    if args.selection_string == "all":
        reduced_atom = u
    else:
        reduced_atom = MDAHelper.new_U_select_atom(u, args.selection_string)
        reduced_atom_name = f"{len(reduced_atom.trajectory)}_frame_dump_atom_selection"
        MDAHelper.write_universe(reduced_atom, reduced_atom_name)

    # Scan system for molecules and select levels (united atom, residue, polymer)
    # for each
    number_molecules, levels = LF.select_levels(reduced_atom, args.verbose)

    # Loop over molecules
    for molecule in range(number_molecules):
        # molecule data container of MDAnalysis Universe type for internal degrees
        # of freedom getting indices of first and last atoms in the molecule
        # assuming atoms are numbered consecutively and all atoms in a given
        # molecule are together
        index1 = reduced_atom.atoms.fragments[molecule].indices[0]
        index2 = reduced_atom.atoms.fragments[molecule].indices[-1]
        selection_string = f"index {index1}:{index2}"
        molecule_container = MDAHelper.new_U_select_atom(reduced_atom, selection_string)

        # Calculate entropy for each relevent level
        for level in levels[molecule]:
            if level == levels[molecule][-1]:
                highest_level = True
            else:
                highest_level = False

            if level == "united_atom":
                # loop over residues, report results per residue + total united atom
                # level. This is done per residue to reduce the size of the matrices -
                # amino acid resiudes have tens of united atoms but a whole protein
                # could have thousands. Doing the calculation per residue allows for
                # comparisons of contributions from different residues
                num_residues = len(molecule_container.residues)
                S_trans = 0
                S_rot = 0
                S_conf = 0

                for residue in range(num_residues):
                    # molecule data container of MDAnalysis Universe type for internal
                    # degrees of freedom getting indices of first and last atoms in the
                    # molecule assuming atoms are numbered consecutively and all atoms
                    # in a given molecule are together
                    index1 = molecule_container.residues[residue].atoms.indices[0]
                    index2 = molecule_container.residues[residue].atoms.indices[-1]
                    selection_string = f"index {index1}:{index2}"
                    residue_container = MDAHelper.new_U_select_atom(
                        molecule_container, selection_string
                    )
                    residue_heavy_atoms_container = MDAHelper.new_U_select_atom(
                        residue_container, "not name H*"
                    )  # only heavy atom dihedrals are relevant

                    # Vibrational entropy at every level
                    # Get the force and torque matrices for the beads at the relevant
                    # level

                    force_matrix, torque_matrix = LF.get_matrices(
                        residue_container,
                        level,
                        args.verbose,
                        start,
                        end,
                        step,
                        number_frames,
                        highest_level,
                    )

                    # Calculate the entropy from the diagonalisation of the matrices
                    S_trans_residue = EF.vibrational_entropy(
                        force_matrix, "force", args.temperature, highest_level
                    )
                    S_trans += S_trans_residue
                    logger.debug(f"S_trans_{level}_{residue} = {S_trans_residue}")
                    new_row = pd.DataFrame(
                        {
                            "Molecule ID": [molecule],
                            "Residue": [residue],
                            "Type": ["Transvibrational (J/mol/K)"],
                            "Result": [S_trans_residue],
                        }
                    )
                    residue_results_df = pd.concat(
                        [residue_results_df, new_row], ignore_index=True
                    )
                    data_logger.add_residue_data(
                        molecule, residue, "Transvibrational", S_trans_residue
                    )

                    S_rot_residue = EF.vibrational_entropy(
                        torque_matrix, "torque", args.temperature, highest_level
                    )
                    S_rot += S_rot_residue
                    logger.debug(f"S_rot_{level}_{residue} = {S_rot_residue}")
                    new_row = pd.DataFrame(
                        {
                            "Molecule ID": [molecule],
                            "Residue": [residue],
                            "Type": ["Rovibrational (J/mol/K)"],
                            "Result": [S_rot_residue],
                        }
                    )
                    residue_results_df = pd.concat(
                        [residue_results_df, new_row], ignore_index=True
                    )
                    data_logger.add_residue_data(
                        molecule, residue, "Rovibrational", S_rot_residue
                    )

                    # Conformational entropy based on atom dihedral angle distributions
                    # Gives entropy of conformations within each residue

                    # Get dihedral angle distribution
                    dihedrals = LF.get_dihedrals(residue_heavy_atoms_container, level)

                    # Calculate conformational entropy
                    S_conf_residue = EF.conformational_entropy(
                        residue_heavy_atoms_container,
                        dihedrals,
                        bin_width,
                        start,
                        end,
                        step,
                        number_frames,
                    )
                    S_conf += S_conf_residue
                    logger.debug(f"S_conf_{level}_{residue} = {S_conf_residue}")
                    new_row = pd.DataFrame(
                        {
                            "Molecule ID": [molecule],
                            "Residue": [residue],
                            "Type": ["Conformational (J/mol/K)"],
                            "Result": [S_conf_residue],
                        }
                    )
                    residue_results_df = pd.concat(
                        [residue_results_df, new_row], ignore_index=True
                    )
                    data_logger.add_residue_data(
                        molecule, residue, "Conformational", S_conf_residue
                    )

                # Print united atom level results summed over all residues
                logger.debug(f"S_trans_{level} = {S_trans}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Transvibrational (J/mol/K)"],
                        "Result": [S_trans],
                    }
                )

                results_df = pd.concat([results_df, new_row], ignore_index=True)

                data_logger.add_results_data(
                    molecule, level, "Transvibrational", S_trans
                )

                logger.debug(f"S_rot_{level} = {S_rot}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Rovibrational (J/mol/K)"],
                        "Result": [S_rot],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)

                data_logger.add_results_data(molecule, level, "Rovibrational", S_rot)
                logger.debug(f"S_conf_{level} = {S_conf}")

                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Conformational (J/mol/K)"],
                        "Result": [S_conf],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)

                data_logger.add_results_data(molecule, level, "Conformational", S_conf)

            if level in ("polymer", "residue"):
                # Vibrational entropy at every level
                # Get the force and torque matrices for the beads at the relevant level
                force_matrix, torque_matrix = LF.get_matrices(
                    molecule_container,
                    level,
                    args.verbose,
                    start,
                    end,
                    step,
                    number_frames,
                    highest_level,
                )

                # Calculate the entropy from the diagonalisation of the matrices
                S_trans = EF.vibrational_entropy(
                    force_matrix, "force", args.temperature, highest_level
                )
                logger.debug(f"S_trans_{level} = {S_trans}")

                # Create new row as a DataFrame for Transvibrational
                new_row_trans = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Transvibrational (J/mol/K)"],
                        "Result": [S_trans],
                    }
                )

                # Concatenate the new row to the DataFrame
                results_df = pd.concat([results_df, new_row_trans], ignore_index=True)

                # Calculate the entropy for Rovibrational
                S_rot = EF.vibrational_entropy(
                    torque_matrix, "torque", args.temperature, highest_level
                )
                logger.debug(f"S_rot_{level} = {S_rot}")

                # Create new row as a DataFrame for Rovibrational
                new_row_rot = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Rovibrational (J/mol/K)"],
                        "Result": [S_rot],
                    }
                )

                # Concatenate the new row to the DataFrame
                results_df = pd.concat([results_df, new_row_rot], ignore_index=True)

                data_logger.add_results_data(
                    molecule, level, "Transvibrational", S_trans
                )
                data_logger.add_results_data(molecule, level, "Rovibrational", S_rot)

                # Note: conformational entropy is not calculated at the polymer level,
                # because there is at most one polymer bead per molecule so no dihedral
                # angles.

            if level == "residue":
                # Conformational entropy based on distributions of dihedral angles
                # of residues. Gives conformational entropy of secondary structure

                # Get dihedral angle distribution
                dihedrals = LF.get_dihedrals(molecule_container, level)
                # Calculate conformational entropy
                S_conf = EF.conformational_entropy(
                    molecule_container,
                    dihedrals,
                    bin_width,
                    start,
                    end,
                    step,
                    number_frames,
                )
                logger.debug(f"S_conf_{level} = {S_conf}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Conformational (J/mol/K)"],
                        "Result": [S_conf],
                    }
                )
                results_df = pd.concat([results_df, new_row], ignore_index=True)
                data_logger.add_results_data(molecule, level, "Conformational", S_conf)

            # Orientational entropy based on network of neighbouring molecules,
            #  only calculated at the highest level (whole molecule)
            #    if highest_level:
            #        neigbours = LF.get_neighbours(reduced_atom, molecule)
            #        S_orient = EF.orientational_entropy(neighbours)
            #        print(f"S_orient_{level} = {S_orient}")
            #        new_row = pd.DataFrame({
            #            'Molecule ID': [molecule],
            #            'Level': [level],
            #            'Type':['Orientational (J/mol/K)'],
            #            'Result': [S_orient],})
            #        results_df = pd.concat([results_df, new_row], ignore_index=True)
            #        with open(args.output_file, "a") as out:
            #    print(molecule,
            #          "\t",
            #          level,
            #          "\tOrientational\t",
            #          S_orient,
            #          file=out)

        # Report total entropy for the molecule
        S_molecule = results_df[results_df["Molecule ID"] == molecule]["Result"].sum()
        logger.debug(f"S_molecule = {S_molecule}")
        new_row = pd.DataFrame(
            {
                "Molecule ID": [molecule],
                "Level": ["Molecule Total"],
                "Type": ["Molecule Total Entropy "],
                "Result": [S_molecule],
            }
        )
        results_df = pd.concat([results_df, new_row], ignore_index=True)

        data_logger.add_results_data(
            molecule, level, "Molecule Total Entropy", S_molecule
        )
        data_logger.save_dataframes_as_json(
            results_df, residue_results_df, args.output_file
        )

    logger.info("Molecules:")
    data_logger.log_tables()


if __name__ == "__main__":

    main()

1	import logging	3✔
2	import math	3✔
3	import os	3✔
4	import re	3✔
5	import sys	3✔
6
7	import MDAnalysis as mda	3✔
8	import pandas as pd	3✔
9
10	from CodeEntropy.calculations import EntropyFunctions as EF	3✔
11	from CodeEntropy.calculations import LevelFunctions as LF	3✔
12	from CodeEntropy.calculations import MDAUniverseHelper as MDAHelper	3✔
13	from CodeEntropy.config.arg_config_manager import ConfigManager	3✔
14	from CodeEntropy.config.data_logger import DataLogger	3✔
15	from CodeEntropy.config.logging_config import LoggingConfig	3✔
16
17
18	def create_job_folder():	3✔
19	"""
20	Create a new job folder with an incremented job number based on existing folders.
21	"""
22	# Get the current working directory
23	base_dir = os.getcwd()	3✔
24
25	# List all folders in the base directory
26	existing_folders = [	3✔
27	f for f in os.listdir(base_dir) if os.path.isdir(os.path.join(base_dir, f))
28	]
29
30	# Filter folders that match the pattern 'jobXXX'
31	job_folders = [f for f in existing_folders if re.match(r"job\d{3}", f)]	3✔
32
33	# Determine the next job number
34	if job_folders:	3✔
35	max_job_number = max([int(re.search(r"\d{3}", f).group()) for f in job_folders])	×
36	next_job_number = max_job_number + 1	×
37	else:
38	next_job_number = 1	3✔
39
40	# Format the new job folder name
41	new_job_folder = f"job{next_job_number:03d}"	3✔
42	new_job_folder_path = os.path.join(base_dir, new_job_folder)	3✔
43
44	# Create the new job folder
45	os.makedirs(new_job_folder_path, exist_ok=True)	3✔
46
47	return new_job_folder_path	3✔
48
49
50	def main():	3✔
51	"""
52	Main function for calculating the entropy of a system using the multiscale cell
53	correlation method.
54	"""
55	folder = create_job_folder()	3✔
56	data_logger = DataLogger()	3✔
57	arg_config = ConfigManager()	3✔
58
59	# Load configuration
60	config = arg_config.load_config("config.yaml")	3✔
61	if config is None:	3✔
62	raise ValueError(	3✔
63	"No configuration file found, and no CLI arguments were provided."
64	)
65
66	parser = arg_config.setup_argparse()	3✔
67	args, unknown = parser.parse_known_args()	3✔
68	args.output_file = os.path.join(folder, args.output_file)	3✔
69
70	try:	3✔
71	# Initialize the logging system once
72	logging_config = LoggingConfig(folder)	3✔
73	logger = logging_config.setup_logging()	3✔
74
75	# Process each run in the YAML configuration
76	for run_name, run_config in config.items():	3✔
77	if isinstance(run_config, dict):	3✔
78	# Merging CLI arguments with YAML configuration
79	args = arg_config.merge_configs(args, run_config)	3✔
80
81	# Determine logging level
82	log_level = logging.DEBUG if args.verbose else logging.INFO	3✔
83
84	# Update the logging level
85	logging_config.update_logging_level(log_level)	3✔
86
87	# Capture and log the command-line invocation
88	command = " ".join(sys.argv)	3✔
89	logging.getLogger("commands").info(command)	3✔
90
91	# Ensure necessary arguments are provided
92	if not getattr(args, "top_traj_file"):	3✔
93	raise ValueError(	×
94	"The 'top_traj_file' argument is required but not provided."
95	)
96	if not getattr(args, "selection_string"):	3✔
97	raise ValueError(	×
98	"The 'selection_string' argument is required but not provided."
99	)
100
101	# Log all inputs for the current run
102	logger.info(f"All input for {run_name}")	3✔
103	for arg in vars(args):	3✔
104	logger.info(f" {arg}: {getattr(args, arg) or ''}")	3✔
105	else:
106	logger.warning(f"Run configuration for {run_name} is not a dictionary.")	×
107	except ValueError as e:	×
108	logger.error(e)	×
109	raise	×
110
111	# Get topology and trajectory file names and make universe
112	tprfile = args.top_traj_file[0]	3✔
113	trrfile = args.top_traj_file[1:]	3✔
114	u = mda.Universe(tprfile, trrfile)	3✔
115
116	# Define bin_width for histogram from inputs
117	bin_width = args.bin_width	3✔
118
119	# Define trajectory slicing from inputs
120	start = args.start	3✔
121	if start is None:	3✔
122	start = 0	×
123	end = args.end	3✔
124	if end is None:	3✔
125	end = -1	×
126	step = args.step	3✔
127	if step is None:	3✔
128	step = 1	×
129	# Count number of frames, easy if not slicing
130	# MDAnalysis trajectory slicing only includes up to end-1
131	# This works the way we want it to if the whole trajectory is being included
132	if start == 0 and end == -1 and step == 1:	3✔
133	end = len(u.trajectory)	3✔
134	number_frames = len(u.trajectory)	3✔
135	elif end == -1:	×
136	end = len(u.trajectory)	×
NEW 137	number_frames = math.floor((end - start) / step)	×
138	else:
NEW 139	end = end + 1	×
NEW 140	number_frames = math.floor((end - start) / step)	×
141	logger.debug(f"Number of Frames: {number_frames}")	3✔
142
143	# Create pandas data frame for results
144	results_df = pd.DataFrame(columns=["Molecule ID", "Level", "Type", "Result"])	3✔
145	residue_results_df = pd.DataFrame(	3✔
146	columns=["Molecule ID", "Residue", "Type", "Result"]
147	)
148
149	# Reduce number of atoms in MDA universe to selection_string arg
150	# (default all atoms included)
151	if args.selection_string == "all":	3✔
152	reduced_atom = u	3✔
153	else:
154	reduced_atom = MDAHelper.new_U_select_atom(u, args.selection_string)	×
155	reduced_atom_name = f"{len(reduced_atom.trajectory)}_frame_dump_atom_selection"	×
156	MDAHelper.write_universe(reduced_atom, reduced_atom_name)	×
157
158	# Scan system for molecules and select levels (united atom, residue, polymer)
159	# for each
160	number_molecules, levels = LF.select_levels(reduced_atom, args.verbose)	3✔
161
162	# Loop over molecules
163	for molecule in range(number_molecules):	3✔
164	# molecule data container of MDAnalysis Universe type for internal degrees
165	# of freedom getting indices of first and last atoms in the molecule
166	# assuming atoms are numbered consecutively and all atoms in a given
167	# molecule are together
168	index1 = reduced_atom.atoms.fragments[molecule].indices[0]	×
169	index2 = reduced_atom.atoms.fragments[molecule].indices[-1]	×
170	selection_string = f"index {index1}:{index2}"	×
171	molecule_container = MDAHelper.new_U_select_atom(reduced_atom, selection_string)	×
172
173	# Calculate entropy for each relevent level
174	for level in levels[molecule]:	×
175	if level == levels[molecule][-1]:	×
176	highest_level = True	×
177	else:
178	highest_level = False	×
179
180	if level == "united_atom":	×
181	# loop over residues, report results per residue + total united atom
182	# level. This is done per residue to reduce the size of the matrices -
183	# amino acid resiudes have tens of united atoms but a whole protein
184	# could have thousands. Doing the calculation per residue allows for
185	# comparisons of contributions from different residues
186	num_residues = len(molecule_container.residues)	×
187	S_trans = 0	×
188	S_rot = 0	×
189	S_conf = 0	×
190
191	for residue in range(num_residues):	×
192	# molecule data container of MDAnalysis Universe type for internal
193	# degrees of freedom getting indices of first and last atoms in the
194	# molecule assuming atoms are numbered consecutively and all atoms
195	# in a given molecule are together
196	index1 = molecule_container.residues[residue].atoms.indices[0]	×
197	index2 = molecule_container.residues[residue].atoms.indices[-1]	×
198	selection_string = f"index {index1}:{index2}"	×
199	residue_container = MDAHelper.new_U_select_atom(	×
200	molecule_container, selection_string
201	)
202	residue_heavy_atoms_container = MDAHelper.new_U_select_atom(	×
203	residue_container, "not name H*"
204	) # only heavy atom dihedrals are relevant
205
206	# Vibrational entropy at every level
207	# Get the force and torque matrices for the beads at the relevant
208	# level
209
210	force_matrix, torque_matrix = LF.get_matrices(	×
211	residue_container,
212	level,
213	args.verbose,
214	start,
215	end,
216	step,
217	number_frames,
218	highest_level,
219	)
220
221	# Calculate the entropy from the diagonalisation of the matrices
222	S_trans_residue = EF.vibrational_entropy(	×
223	force_matrix, "force", args.temperature, highest_level
224	)
225	S_trans += S_trans_residue	×
226	logger.debug(f"S_trans_{level}_{residue} = {S_trans_residue}")	×
227	new_row = pd.DataFrame(	×
228	{
229	"Molecule ID": [molecule],
230	"Residue": [residue],
231	"Type": ["Transvibrational (J/mol/K)"],
232	"Result": [S_trans_residue],
233	}
234	)
235	residue_results_df = pd.concat(	×
236	[residue_results_df, new_row], ignore_index=True
237	)
238	data_logger.add_residue_data(	×
239	molecule, residue, "Transvibrational", S_trans_residue
240	)
241
242	S_rot_residue = EF.vibrational_entropy(	×
243	torque_matrix, "torque", args.temperature, highest_level
244	)
245	S_rot += S_rot_residue	×
246	logger.debug(f"S_rot_{level}_{residue} = {S_rot_residue}")	×
247	new_row = pd.DataFrame(	×
248	{
249	"Molecule ID": [molecule],
250	"Residue": [residue],
251	"Type": ["Rovibrational (J/mol/K)"],
252	"Result": [S_rot_residue],
253	}
254	)
255	residue_results_df = pd.concat(	×
256	[residue_results_df, new_row], ignore_index=True
257	)
258	data_logger.add_residue_data(	×
259	molecule, residue, "Rovibrational", S_rot_residue
260	)
261
262	# Conformational entropy based on atom dihedral angle distributions
263	# Gives entropy of conformations within each residue
264
265	# Get dihedral angle distribution
266	dihedrals = LF.get_dihedrals(residue_heavy_atoms_container, level)	×
267
268	# Calculate conformational entropy
269	S_conf_residue = EF.conformational_entropy(	×
270	residue_heavy_atoms_container,
271	dihedrals,
272	bin_width,
273	start,
274	end,
275	step,
276	number_frames,
277	)
278	S_conf += S_conf_residue	×
279	logger.debug(f"S_conf_{level}_{residue} = {S_conf_residue}")	×
280	new_row = pd.DataFrame(	×
281	{
282	"Molecule ID": [molecule],
283	"Residue": [residue],
284	"Type": ["Conformational (J/mol/K)"],
285	"Result": [S_conf_residue],
286	}
287	)
288	residue_results_df = pd.concat(	×
289	[residue_results_df, new_row], ignore_index=True
290	)
291	data_logger.add_residue_data(	×
292	molecule, residue, "Conformational", S_conf_residue
293	)
294
295	# Print united atom level results summed over all residues
296	logger.debug(f"S_trans_{level} = {S_trans}")	×
297	new_row = pd.DataFrame(	×
298	{
299	"Molecule ID": [molecule],
300	"Level": [level],
301	"Type": ["Transvibrational (J/mol/K)"],
302	"Result": [S_trans],
303	}
304	)
305
306	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
307
308	data_logger.add_results_data(	×
309	molecule, level, "Transvibrational", S_trans
310	)
311
312	logger.debug(f"S_rot_{level} = {S_rot}")	×
313	new_row = pd.DataFrame(	×
314	{
315	"Molecule ID": [molecule],
316	"Level": [level],
317	"Type": ["Rovibrational (J/mol/K)"],
318	"Result": [S_rot],
319	}
320	)
321	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
322
323	data_logger.add_results_data(molecule, level, "Rovibrational", S_rot)	×
324	logger.debug(f"S_conf_{level} = {S_conf}")	×
325
326	new_row = pd.DataFrame(	×
327	{
328	"Molecule ID": [molecule],
329	"Level": [level],
330	"Type": ["Conformational (J/mol/K)"],
331	"Result": [S_conf],
332	}
333	)
334	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
335
336	data_logger.add_results_data(molecule, level, "Conformational", S_conf)	×
337
338	if level in ("polymer", "residue"):	×
339	# Vibrational entropy at every level
340	# Get the force and torque matrices for the beads at the relevant level
341	force_matrix, torque_matrix = LF.get_matrices(	×
342	molecule_container,
343	level,
344	args.verbose,
345	start,
346	end,
347	step,
348	number_frames,
349	highest_level,
350	)
351
352	# Calculate the entropy from the diagonalisation of the matrices
353	S_trans = EF.vibrational_entropy(	×
354	force_matrix, "force", args.temperature, highest_level
355	)
356	logger.debug(f"S_trans_{level} = {S_trans}")	×
357
358	# Create new row as a DataFrame for Transvibrational
359	new_row_trans = pd.DataFrame(	×
360	{
361	"Molecule ID": [molecule],
362	"Level": [level],
363	"Type": ["Transvibrational (J/mol/K)"],
364	"Result": [S_trans],
365	}
366	)
367
368	# Concatenate the new row to the DataFrame
369	results_df = pd.concat([results_df, new_row_trans], ignore_index=True)	×
370
371	# Calculate the entropy for Rovibrational
372	S_rot = EF.vibrational_entropy(	×
373	torque_matrix, "torque", args.temperature, highest_level
374	)
375	logger.debug(f"S_rot_{level} = {S_rot}")	×
376
377	# Create new row as a DataFrame for Rovibrational
378	new_row_rot = pd.DataFrame(	×
379	{
380	"Molecule ID": [molecule],
381	"Level": [level],
382	"Type": ["Rovibrational (J/mol/K)"],
383	"Result": [S_rot],
384	}
385	)
386
387	# Concatenate the new row to the DataFrame
388	results_df = pd.concat([results_df, new_row_rot], ignore_index=True)	×
389
390	data_logger.add_results_data(	×
391	molecule, level, "Transvibrational", S_trans
392	)
393	data_logger.add_results_data(molecule, level, "Rovibrational", S_rot)	×
394
395	# Note: conformational entropy is not calculated at the polymer level,
396	# because there is at most one polymer bead per molecule so no dihedral
397	# angles.
398
399	if level == "residue":	×
400	# Conformational entropy based on distributions of dihedral angles
401	# of residues. Gives conformational entropy of secondary structure
402
403	# Get dihedral angle distribution
404	dihedrals = LF.get_dihedrals(molecule_container, level)	×
405	# Calculate conformational entropy
406	S_conf = EF.conformational_entropy(	×
407	molecule_container,
408	dihedrals,
409	bin_width,
410	start,
411	end,
412	step,
413	number_frames,
414	)
415	logger.debug(f"S_conf_{level} = {S_conf}")	×
416	new_row = pd.DataFrame(	×
417	{
418	"Molecule ID": [molecule],
419	"Level": [level],
420	"Type": ["Conformational (J/mol/K)"],
421	"Result": [S_conf],
422	}
423	)
424	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
425	data_logger.add_results_data(molecule, level, "Conformational", S_conf)	×
426
427	# Orientational entropy based on network of neighbouring molecules,
428	# only calculated at the highest level (whole molecule)
429	# if highest_level:
430	# neigbours = LF.get_neighbours(reduced_atom, molecule)
431	# S_orient = EF.orientational_entropy(neighbours)
432	# print(f"S_orient_{level} = {S_orient}")
433	# new_row = pd.DataFrame({
434	# 'Molecule ID': [molecule],
435	# 'Level': [level],
436	# 'Type':['Orientational (J/mol/K)'],
437	# 'Result': [S_orient],})
438	# results_df = pd.concat([results_df, new_row], ignore_index=True)
439	# with open(args.output_file, "a") as out:
440	# print(molecule,
441	# "\t",
442	# level,
443	# "\tOrientational\t",
444	# S_orient,
445	# file=out)
446
447	# Report total entropy for the molecule
448	S_molecule = results_df[results_df["Molecule ID"] == molecule]["Result"].sum()	×
449	logger.debug(f"S_molecule = {S_molecule}")	×
450	new_row = pd.DataFrame(	×
451	{
452	"Molecule ID": [molecule],
453	"Level": ["Molecule Total"],
454	"Type": ["Molecule Total Entropy "],
455	"Result": [S_molecule],
456	}
457	)
458	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
459
460	data_logger.add_results_data(	×
461	molecule, level, "Molecule Total Entropy", S_molecule
462	)
463	data_logger.save_dataframes_as_json(	×
464	results_df, residue_results_df, args.output_file
465	)
466
467	logger.info("Molecules:")	3✔
468	data_logger.log_tables()	3✔
469
470
471	if __name__ == "__main__":	3✔
472
473	main()	×

CCPBioSim / CodeEntropy / 14928438681

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