13973993966

Committed 20 Mar 2025 04:08PM UTC coverage: 29.982%. Remained the same

Build # 13973993966

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #70 from CCPBioSim/49-remove-end-method-comments

Removed all instances of # End Method comments

Run Details

164 of 547 relevant lines covered (29.98%)

0.9 hits per line

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

44.44

/CodeEntropy/main_mcc.py

import argparse
import math
import os

import MDAnalysis as mda

# import numpy as np
import pandas as pd
import yaml

from CodeEntropy import EntropyFunctions as EF
from CodeEntropy import LevelFunctions as LF
from CodeEntropy import MDAUniverseHelper as MDAHelper

# from datetime import datetime

arg_map = {
    "top_traj_file": {
        "type": str,
        "nargs": "+",
        "help": "Path to Structure/topology file followed by Trajectory file(s)",
        "default": [],
    },
    "selection_string": {
        "type": str,
        "help": "Selection string for CodeEntropy",
        "default": "all",
    },
    "start": {
        "type": int,
        "help": "Start analysing the trajectory from this frame index",
        "default": 0,
    },
    "end": {
        "type": int,
        "help": "Stop analysing the trajectory at this frame index",
        "default": -1,
    },
    "step": {
        "type": int,
        "help": "Interval between two consecutive frames to be read index",
        "default": 1,
    },
    "bin_width": {
        "type": int,
        "help": "Bin width in degrees for making the histogram",
        "default": 30,
    },
    "temperature": {
        "type": float,
        "help": "Temperature for entropy calculation (K)",
        "default": 298.0,
    },
    "verbose": {
        "type": bool,
        "help": "True/False flag for noisy or quiet output",
        "default": False,
    },
    "thread": {"type": int, "help": "How many multiprocess to use", "default": 1},
    "outfile": {
        "type": str,
        "help": "Name of the file where the output will be written",
        "default": "outfile.out",
    },
    "resfile": {
        "type": str,
        "help": "Name of the file where the residue entropy output will be written",
        "default": "res_outfile.out",
    },
    "mout": {
        "type": str,
        "help": "Name of the file where certain matrices will be written",
        "default": None,
    },
    "force_partitioning": {"type": float, "help": "Force partitioning", "default": 0.5},
    "waterEntropy": {"type": bool, "help": "Calculate water entropy", "default": False},
}


def load_config(file_path):
    """Load YAML configuration file."""
    if not os.path.exists(file_path):
        raise FileNotFoundError(f"Configuration file '{file_path}' not found.")

    with open(file_path, "r") as file:
        config = yaml.safe_load(file)

        # If YAML content is empty, return an empty dictionary
        if config is None:
            config = {}

    return config


def setup_argparse():
    """Setup argument parsing dynamically based on arg_map."""
    parser = argparse.ArgumentParser(
        description="CodeEntropy: Entropy calculation with MCC method."
    )

    for arg, properties in arg_map.items():
        kwargs = {key: properties[key] for key in properties if key != "help"}
        parser.add_argument(f"--{arg}", **kwargs, help=properties.get("help"))

    return parser


def merge_configs(args, run_config):
    """Merge CLI arguments with YAML configuration."""
    if run_config is None:
        run_config = {}

    if not isinstance(run_config, dict):
        raise TypeError("run_config must be a dictionary or None.")

    # Step 1: Merge YAML configuration into args
    for key, value in run_config.items():
        if getattr(args, key, None) is None:
            setattr(args, key, value)

    # Step 2: Set default values for any missing arguments from `arg_map`
    for key, params in arg_map.items():
        if getattr(args, key, None) is None:
            setattr(args, key, params.get("default"))

    # Step 3: Override with CLI values if provided
    for key in arg_map.keys():
        cli_value = getattr(args, key, None)
        if cli_value is not None:
            run_config[key] = cli_value

    return args


def main():
    """
    Main function for calculating the entropy of a system using the multiscale cell
    correlation method.
    """
    try:
        config = load_config("config.yaml")

        if config is None:
            raise ValueError(
                "No configuration file found, and no CLI arguments were provided."
            )

        parser = setup_argparse()
        args, unknown = parser.parse_known_args()

        # 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 = 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."
                    )

                # REPLACE INPUTS
                print(f"Printing all input for {run_name}")
                for arg in vars(args):
                    print(f" {arg}: {getattr(args, arg) or ''}")
            else:
                print(f"Run configuration for {run_name} is not a dictionary.")
    except ValueError as e:
        print(e)
        raise

    # startTime = datetime.now()

    # 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
    print(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"]
    )

    # printing headings for output files
    with open(args.outfile, "a") as out:
        print("Molecule\tLevel\tType\tResult (J/mol/K)\n", file=out)

    with open(args.resfile, "a") as res:
        print("Molecule\tResidue\tType\tResult (J/mol/K)\n", file=res)

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

                    # 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
                    print(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
                    )
                    with open(args.resfile, "a") as res:
                        print(
                            molecule,
                            "\t",
                            residue,
                            "\tTransvibration\t",
                            S_trans_residue,
                            file=res,
                        )

                    S_rot_residue = EF.vibrational_entropy(
                        torque_matrix, "torque", args.temperature, highest_level
                    )
                    S_rot += S_rot_residue
                    print(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
                    )
                    with open(args.resfile, "a") as res:
                        #  print(new_row, file=res)
                        print(
                            molecule,
                            "\t",
                            residue,
                            "\tRovibrational \t",
                            S_rot_residue,
                            file=res,
                        )

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

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

                    # Calculate conformational entropy
                    S_conf_residue = EF.conformational_entropy(
                        residue_container,
                        dihedrals,
                        bin_width,
                        start,
                        end,
                        step,
                        number_frames,
                    )
                    S_conf += S_conf_residue
                    print(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
                    )
                    with open(args.resfile, "a") as res:
                        print(
                            molecule,
                            "\t",
                            residue,
                            "\tConformational\t",
                            S_conf_residue,
                            file=res,
                        )

                # Print united atom level results summed over all residues
                print(f"S_trans_{level} = {S_trans}")
                new_row = pd.DataFrame(
                    {
                        "Molecule ID": [molecule],
                        "Level": [level],
                        "Type": ["Transvibrational (J/mol/K)"],
                        "Result": [S_trans],
                    }
                )
                with open(args.outfile, "a") as out:
                    print(
                        molecule, "\t", level, "\tTransvibration\t", S_trans, file=out
                    )

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

                print(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)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)

                print(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)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)

            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
                )
                print(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)
                with open(args.outfile, "a") as out:
                    print(
                        molecule, "\t", level, "\tTransvibrational\t", S_trans, file=out
                    )

                S_rot = EF.vibrational_entropy(
                    torque_matrix, "torque", args.temperature, highest_level
                )
                print(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)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)

                # 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,
                )
                print(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)
                with open(args.outfile, "a") as out:
                    print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)

            # 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()
        print(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)
        with open(args.outfile, "a") as out:
            print(molecule, "\t Molecule\tTotal Entropy\t", S_molecule, file=out)


if __name__ == "__main__":

    main()

1	import argparse	3✔
2	import math	3✔
3	import os	3✔
4
5	import MDAnalysis as mda	3✔
6
7	# import numpy as np
8	import pandas as pd	3✔
9	import yaml	3✔
10
11	from CodeEntropy import EntropyFunctions as EF	3✔
12	from CodeEntropy import LevelFunctions as LF	3✔
13	from CodeEntropy import MDAUniverseHelper as MDAHelper	3✔
14
15	# from datetime import datetime
16
17	arg_map = {	3✔
18	"top_traj_file": {
19	"type": str,
20	"nargs": "+",
21	"help": "Path to Structure/topology file followed by Trajectory file(s)",
22	"default": [],
23	},
24	"selection_string": {
25	"type": str,
26	"help": "Selection string for CodeEntropy",
27	"default": "all",
28	},
29	"start": {
30	"type": int,
31	"help": "Start analysing the trajectory from this frame index",
32	"default": 0,
33	},
34	"end": {
35	"type": int,
36	"help": "Stop analysing the trajectory at this frame index",
37	"default": -1,
38	},
39	"step": {
40	"type": int,
41	"help": "Interval between two consecutive frames to be read index",
42	"default": 1,
43	},
44	"bin_width": {
45	"type": int,
46	"help": "Bin width in degrees for making the histogram",
47	"default": 30,
48	},
49	"temperature": {
50	"type": float,
51	"help": "Temperature for entropy calculation (K)",
52	"default": 298.0,
53	},
54	"verbose": {
55	"type": bool,
56	"help": "True/False flag for noisy or quiet output",
57	"default": False,
58	},
59	"thread": {"type": int, "help": "How many multiprocess to use", "default": 1},
60	"outfile": {
61	"type": str,
62	"help": "Name of the file where the output will be written",
63	"default": "outfile.out",
64	},
65	"resfile": {
66	"type": str,
67	"help": "Name of the file where the residue entropy output will be written",
68	"default": "res_outfile.out",
69	},
70	"mout": {
71	"type": str,
72	"help": "Name of the file where certain matrices will be written",
73	"default": None,
74	},
75	"force_partitioning": {"type": float, "help": "Force partitioning", "default": 0.5},
76	"waterEntropy": {"type": bool, "help": "Calculate water entropy", "default": False},
77	}
78
79
80	def load_config(file_path):	3✔
81	"""Load YAML configuration file."""
82	if not os.path.exists(file_path):	3✔
83	raise FileNotFoundError(f"Configuration file '{file_path}' not found.")	×
84
85	with open(file_path, "r") as file:	3✔
86	config = yaml.safe_load(file)	3✔
87
88	# If YAML content is empty, return an empty dictionary
89	if config is None:	3✔
90	config = {}	3✔
91
92	return config	3✔
93
94
95	def setup_argparse():	3✔
96	"""Setup argument parsing dynamically based on arg_map."""
97	parser = argparse.ArgumentParser(	3✔
98	description="CodeEntropy: Entropy calculation with MCC method."
99	)
100
101	for arg, properties in arg_map.items():	3✔
102	kwargs = {key: properties[key] for key in properties if key != "help"}	3✔
103	parser.add_argument(f"--{arg}", **kwargs, help=properties.get("help"))	3✔
104
105	return parser	3✔
106
107
108	def merge_configs(args, run_config):	3✔
109	"""Merge CLI arguments with YAML configuration."""
110	if run_config is None:	3✔
111	run_config = {}	×
112
113	if not isinstance(run_config, dict):	3✔
114	raise TypeError("run_config must be a dictionary or None.")	3✔
115
116	# Step 1: Merge YAML configuration into args
117	for key, value in run_config.items():	3✔
118	if getattr(args, key, None) is None:	3✔
119	setattr(args, key, value)	3✔
120
121	# Step 2: Set default values for any missing arguments from `arg_map`
122	for key, params in arg_map.items():	3✔
123	if getattr(args, key, None) is None:	3✔
124	setattr(args, key, params.get("default"))	3✔
125
126	# Step 3: Override with CLI values if provided
127	for key in arg_map.keys():	3✔
128	cli_value = getattr(args, key, None)	3✔
129	if cli_value is not None:	3✔
130	run_config[key] = cli_value	3✔
131
132	return args	3✔
133
134
135	def main():	3✔
136	"""
137	Main function for calculating the entropy of a system using the multiscale cell
138	correlation method.
139	"""
140	try:	3✔
141	config = load_config("config.yaml")	3✔
142
143	if config is None:	3✔
144	raise ValueError(	3✔
145	"No configuration file found, and no CLI arguments were provided."
146	)
147
148	parser = setup_argparse()	3✔
149	args, unknown = parser.parse_known_args()	3✔
150
151	# Process each run in the YAML configuration
152	for run_name, run_config in config.items():	3✔
153	if isinstance(run_config, dict):	3✔
154	# Merging CLI arguments with YAML configuration
155	args = merge_configs(args, run_config)	3✔
156
157	# Ensure necessary arguments are provided
158	if not getattr(args, "top_traj_file"):	3✔
159	raise ValueError(	×
160	"The 'top_traj_file' argument is required but not provided."
161	)
162	if not getattr(args, "selection_string"):	3✔
163	raise ValueError(	×
164	"The 'selection_string' argument is required but not provided."
165	)
166
167	# REPLACE INPUTS
168	print(f"Printing all input for {run_name}")	3✔
169	for arg in vars(args):	3✔
170	print(f" {arg}: {getattr(args, arg) or ''}")	3✔
171	else:
172	print(f"Run configuration for {run_name} is not a dictionary.")	×
173	except ValueError as e:	3✔
174	print(e)	3✔
175	raise	3✔
176
177	# startTime = datetime.now()
178
179	# Get topology and trajectory file names and make universe
180	tprfile = args.top_traj_file[0]	3✔
181	trrfile = args.top_traj_file[1:]	3✔
182	u = mda.Universe(tprfile, trrfile)	3✔
183
184	# Define bin_width for histogram from inputs
185	bin_width = args.bin_width	3✔
186
187	# Define trajectory slicing from inputs
188	start = args.start	3✔
189	if start is None:	3✔
190	start = 0	×
191	end = args.end	3✔
192	if end is None:	3✔
193	end = -1	×
194	step = args.step	3✔
195	if step is None:	3✔
196	step = 1	×
197	# Count number of frames, easy if not slicing
198	if start == 0 and end == -1 and step == 1:	3✔
199	number_frames = len(u.trajectory)	3✔
200	elif end == -1:	×
201	end = len(u.trajectory)	×
202	number_frames = math.floor((end - start) / step) + 1	×
203	else:
204	number_frames = math.floor((end - start) / step) + 1	×
205	print(number_frames)	3✔
206
207	# Create pandas data frame for results
208	results_df = pd.DataFrame(columns=["Molecule ID", "Level", "Type", "Result"])	3✔
209	residue_results_df = pd.DataFrame(	3✔
210	columns=["Molecule ID", "Residue", "Type", "Result"]
211	)
212
213	# printing headings for output files
214	with open(args.outfile, "a") as out:	3✔
215	print("Molecule\tLevel\tType\tResult (J/mol/K)\n", file=out)	3✔
216
217	with open(args.resfile, "a") as res:	3✔
218	print("Molecule\tResidue\tType\tResult (J/mol/K)\n", file=res)	3✔
219
220	# Reduce number of atoms in MDA universe to selection_string arg
221	# (default all atoms included)
222	if args.selection_string == "all":	3✔
223	reduced_atom = u	3✔
224	else:
225	reduced_atom = MDAHelper.new_U_select_atom(u, args.selection_string)	×
226	reduced_atom_name = f"{len(reduced_atom.trajectory)}_frame_dump_atom_selection"	×
227	MDAHelper.write_universe(reduced_atom, reduced_atom_name)	×
228
229	# Scan system for molecules and select levels (united atom, residue, polymer)
230	# for each
231	number_molecules, levels = LF.select_levels(reduced_atom, args.verbose)	3✔
232
233	# Loop over molecules
234	for molecule in range(number_molecules):	3✔
235	# molecule data container of MDAnalysis Universe type for internal degrees
236	# of freedom getting indices of first and last atoms in the molecule
237	# assuming atoms are numbered consecutively and all atoms in a given
238	# molecule are together
239	index1 = reduced_atom.atoms.fragments[molecule].indices[0]	×
240	index2 = reduced_atom.atoms.fragments[molecule].indices[-1]	×
241	selection_string = f"index {index1}:{index2}"	×
242	molecule_container = MDAHelper.new_U_select_atom(reduced_atom, selection_string)	×
243
244	# Calculate entropy for each relevent level
245	for level in levels[molecule]:	×
246	if level == levels[molecule][-1]:	×
247	highest_level = True	×
248	else:
249	highest_level = False	×
250
251	if level == "united_atom":	×
252	# loop over residues, report results per residue + total united atom
253	# level. This is done per residue to reduce the size of the matrices -
254	# amino acid resiudes have tens of united atoms but a whole protein
255	# could have thousands. Doing the calculation per residue allows for
256	# comparisons of contributions from different residues
257	num_residues = len(molecule_container.residues)	×
258	S_trans = 0	×
259	S_rot = 0	×
260	S_conf = 0	×
261	for residue in range(num_residues):	×
262	# molecule data container of MDAnalysis Universe type for internal
263	# degrees of freedom getting indices of first and last atoms in the
264	# molecule assuming atoms are numbered consecutively and all atoms
265	# in a given molecule are together
266	index1 = molecule_container.residues[residue].atoms.indices[0]	×
267	index2 = molecule_container.residues[residue].atoms.indices[-1]	×
268	selection_string = f"index {index1}:{index2}"	×
269	residue_container = MDAHelper.new_U_select_atom(	×
270	molecule_container, selection_string
271	)
272
273	# Vibrational entropy at every level
274	# Get the force and torque matrices for the beads at the relevant
275	# level
276	force_matrix, torque_matrix = LF.get_matrices(	×
277	residue_container,
278	level,
279	args.verbose,
280	start,
281	end,
282	step,
283	number_frames,
284	highest_level,
285	)
286
287	# Calculate the entropy from the diagonalisation of the matrices
288	S_trans_residue = EF.vibrational_entropy(	×
289	force_matrix, "force", args.temperature, highest_level
290	)
291	S_trans += S_trans_residue	×
292	print(f"S_trans_{level}_{residue} = {S_trans_residue}")	×
293	new_row = pd.DataFrame(	×
294	{
295	"Molecule ID": [molecule],
296	"Residue": [residue],
297	"Type": ["Transvibrational (J/mol/K)"],
298	"Result": [S_trans_residue],
299	}
300	)
301	residue_results_df = pd.concat(	×
302	[residue_results_df, new_row], ignore_index=True
303	)
304	with open(args.resfile, "a") as res:	×
305	print(	×
306	molecule,
307	"\t",
308	residue,
309	"\tTransvibration\t",
310	S_trans_residue,
311	file=res,
312	)
313
314	S_rot_residue = EF.vibrational_entropy(	×
315	torque_matrix, "torque", args.temperature, highest_level
316	)
317	S_rot += S_rot_residue	×
318	print(f"S_rot_{level}_{residue} = {S_rot_residue}")	×
319	new_row = pd.DataFrame(	×
320	{
321	"Molecule ID": [molecule],
322	"Residue": [residue],
323	"Type": ["Rovibrational (J/mol/K)"],
324	"Result": [S_rot_residue],
325	}
326	)
327	residue_results_df = pd.concat(	×
328	[residue_results_df, new_row], ignore_index=True
329	)
330	with open(args.resfile, "a") as res:	×
331	# print(new_row, file=res)
332	print(	×
333	molecule,
334	"\t",
335	residue,
336	"\tRovibrational \t",
337	S_rot_residue,
338	file=res,
339	)
340
341	# Conformational entropy based on atom dihedral angle distributions
342	# Gives entropy of conformations within each residue
343
344	# Get dihedral angle distribution
345	dihedrals = LF.get_dihedrals(residue_container, level)	×
346
347	# Calculate conformational entropy
348	S_conf_residue = EF.conformational_entropy(	×
349	residue_container,
350	dihedrals,
351	bin_width,
352	start,
353	end,
354	step,
355	number_frames,
356	)
357	S_conf += S_conf_residue	×
358	print(f"S_conf_{level}_{residue} = {S_conf_residue}")	×
359	new_row = pd.DataFrame(	×
360	{
361	"Molecule ID": [molecule],
362	"Residue": [residue],
363	"Type": ["Conformational (J/mol/K)"],
364	"Result": [S_conf_residue],
365	}
366	)
367	residue_results_df = pd.concat(	×
368	[residue_results_df, new_row], ignore_index=True
369	)
370	with open(args.resfile, "a") as res:	×
371	print(	×
372	molecule,
373	"\t",
374	residue,
375	"\tConformational\t",
376	S_conf_residue,
377	file=res,
378	)
379
380	# Print united atom level results summed over all residues
381	print(f"S_trans_{level} = {S_trans}")	×
382	new_row = pd.DataFrame(	×
383	{
384	"Molecule ID": [molecule],
385	"Level": [level],
386	"Type": ["Transvibrational (J/mol/K)"],
387	"Result": [S_trans],
388	}
389	)
390	with open(args.outfile, "a") as out:	×
391	print(	×
392	molecule, "\t", level, "\tTransvibration\t", S_trans, file=out
393	)
394
395	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
396
397	print(f"S_rot_{level} = {S_rot}")	×
398	new_row = pd.DataFrame(	×
399	{
400	"Molecule ID": [molecule],
401	"Level": [level],
402	"Type": ["Rovibrational (J/mol/K)"],
403	"Result": [S_rot],
404	}
405	)
406	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
407	with open(args.outfile, "a") as out:	×
408	print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)	×
409
410	print(f"S_conf_{level} = {S_conf}")	×
411	new_row = pd.DataFrame(	×
412	{
413	"Molecule ID": [molecule],
414	"Level": [level],
415	"Type": ["Conformational (J/mol/K)"],
416	"Result": [S_conf],
417	}
418	)
419	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
420	with open(args.outfile, "a") as out:	×
421	print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)	×
422
423	if level in ("polymer", "residue"):	×
424	# Vibrational entropy at every level
425	# Get the force and torque matrices for the beads at the relevant level
426	force_matrix, torque_matrix = LF.get_matrices(	×
427	molecule_container,
428	level,
429	args.verbose,
430	start,
431	end,
432	step,
433	number_frames,
434	highest_level,
435	)
436
437	# Calculate the entropy from the diagonalisation of the matrices
438	S_trans = EF.vibrational_entropy(	×
439	force_matrix, "force", args.temperature, highest_level
440	)
441	print(f"S_trans_{level} = {S_trans}")	×
442	new_row = pd.DataFrame(	×
443	{
444	"Molecule ID": [molecule],
445	"Level": [level],
446	"Type": ["Transvibrational (J/mol/K)"],
447	"Result": [S_trans],
448	}
449	)
450	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
451	with open(args.outfile, "a") as out:	×
452	print(	×
453	molecule, "\t", level, "\tTransvibrational\t", S_trans, file=out
454	)
455
456	S_rot = EF.vibrational_entropy(	×
457	torque_matrix, "torque", args.temperature, highest_level
458	)
459	print(f"S_rot_{level} = {S_rot}")	×
460	new_row = pd.DataFrame(	×
461	{
462	"Molecule ID": [molecule],
463	"Level": [level],
464	"Type": ["Rovibrational (J/mol/K)"],
465	"Result": [S_rot],
466	}
467	)
468	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
469	with open(args.outfile, "a") as out:	×
470	print(molecule, "\t", level, "\tRovibrational \t", S_rot, file=out)	×
471
472	# Note: conformational entropy is not calculated at the polymer level,
473	# because there is at most one polymer bead per molecule so no dihedral
474	# angles.
475
476	if level == "residue":	×
477	# Conformational entropy based on distributions of dihedral angles
478	# of residues. Gives conformational entropy of secondary structure
479
480	# Get dihedral angle distribution
481	dihedrals = LF.get_dihedrals(molecule_container, level)	×
482	# Calculate conformational entropy
483	S_conf = EF.conformational_entropy(	×
484	molecule_container,
485	dihedrals,
486	bin_width,
487	start,
488	end,
489	step,
490	number_frames,
491	)
492	print(f"S_conf_{level} = {S_conf}")	×
493	new_row = pd.DataFrame(	×
494	{
495	"Molecule ID": [molecule],
496	"Level": [level],
497	"Type": ["Conformational (J/mol/K)"],
498	"Result": [S_conf],
499	}
500	)
501	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
502	with open(args.outfile, "a") as out:	×
503	print(molecule, "\t", level, "\tConformational\t", S_conf, file=out)	×
504
505	# Orientational entropy based on network of neighbouring molecules,
506	# only calculated at the highest level (whole molecule)
507	# if highest_level:
508	# neigbours = LF.get_neighbours(reduced_atom, molecule)
509	# S_orient = EF.orientational_entropy(neighbours)
510	# print(f"S_orient_{level} = {S_orient}")
511	# new_row = pd.DataFrame({
512	# 'Molecule ID': [molecule],
513	# 'Level': [level],
514	# 'Type':['Orientational (J/mol/K)'],
515	# 'Result': [S_orient],})
516	# results_df = pd.concat([results_df, new_row], ignore_index=True)
517	# with open(args.outfile, "a") as out:
518	# print(molecule,
519	# "\t",
520	# level,
521	# "\tOrientational\t",
522	# S_orient,
523	# file=out)
524
525	# Report total entropy for the molecule
526	S_molecule = results_df[results_df["Molecule ID"] == molecule]["Result"].sum()	×
527	print(f"S_molecule = {S_molecule}")	×
528	new_row = pd.DataFrame(	×
529	{
530	"Molecule ID": [molecule],
531	"Level": ["Molecule Total"],
532	"Type": ["Molecule Total Entropy "],
533	"Result": [S_molecule],
534	}
535	)
536	results_df = pd.concat([results_df, new_row], ignore_index=True)	×
537	with open(args.outfile, "a") as out:	×
538	print(molecule, "\t Molecule\tTotal Entropy\t", S_molecule, file=out)	×
539
540
541	if __name__ == "__main__":	3✔
542
543	main()	×

CCPBioSim / CodeEntropy / 13973993966

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