14191288198

Committed 01 Apr 2025 08:47AM UTC coverage: 39.213%. Remained the same

Build # 14191288198

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #79 from CCPBioSim/78-test-cases-for-input-arguments

Additional Test Cases for Input Configuration

Run Details

249 of 635 relevant lines covered (39.21%)

1.18 hits per line

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

42.07

/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.outfile = os.path.join(folder, args.outfile)

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

    # Initialize the logging system with the determined log level
    logging_config = LoggingConfig(folder, default_level=log_level)
    logger = logging_config.setup_logging()

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

    try:
        # 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)

                # 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
    if start == 0 and end == -1 and step == 1:
        number_frames = len(u.trajectory)
    elif end == -1:
        end = len(u.trajectory)
        number_frames = math.floor((end - start) / step) + 1
    else:
        number_frames = math.floor((end - start) / step) + 1
    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.outfile, "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.outfile
        )

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

CCPBioSim / CodeEntropy / 14191288198

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